k=4
  tn=4

  library(package="orthopolynom")
  (leg4coef <- legendre.polynomials(n=k-1, normalized=TRUE))

  #for x,y axis in the plot ***************************
    #this for the 4 environment meausred for the phenotypes
    ev=c(0,15,30,75)
    ev2=ev
    ev2=ev2-ev[1]
    x=(ev2-ev2[tn]/2)/(ev2[tn]/2)

    #this is for infinite points between the first and last environment
    #st=-1 
    #fn=1 
    #x=seq(st,fn,by=(fn-st)/(tn-1))
  #****************************************************

  leg4 <- as.matrix(as.data.frame(polynomial.values(polynomials=leg4coef,x=scaleX(x, u=-1, v=1))))

  colnames(leg4) <- seq(0,k-1)
  leg4 <- leg4[, 1:ncol(leg4)]   #same as phi in page 31

  v1=read.table("legp_plot.in")  #same as K in page 31

  km=matrix(0,k,k)
  for (i in 1:k) {
    km[i,i]=v1$V2[i]
  }
  ki=0
  for (i in 2:k) {
    for (j in 1:(i-1)) {
      ki=ki+1
      km[i,j]=v1$V2[k+ki]
      km[j,i]=v1$V2[k+ki]
    }
  }

  vcv=leg4%*%km%*%t(leg4)