## Plot the density for the observed data

plotchi<-function(rmsea, n){
  filename<-paste('allres-',rmsea,'-',n,'-chi.txt', sep='')
  temp<-read.table(filename)
  temp<-temp[!is.na(temp[,3]),3]
  filename1<-paste(rmsea,'-',n,'.txt', sep='')
  write.table(temp,file=filename1,row.names=F, col.names=F)

  nn<-length(temp)

  ## now write the contral file for FPDE
  filename2<-paste(rmsea,'-',n,'.ctr', sep='')
  filename3<-paste(rmsea,'-',n,'.out', sep='')
  cat('filein=',filename1,';\n', file=filename2)
  cat('fileps = plots.eps;\n', file=filename2, append=T)
  cat('options=(Nsample=', nn, ', nboot=1, ntarget=5, Npolydeg=5, Valpar3=3);\n', file=filename2, append=T)
  cat('***\n', file=filename2, append=T)

  ## write the exectuable file
  filename4<-paste(rmsea,'-',n,'.run', sep='')
  cat(filename2, '\n', file=filename4)
  cat(filename3, '\n', file=filename4, append=T)

  ## 
  filename5<-paste('run.bat', sep='')
  cmd <- paste('FPDEAdjust2.exe <', filename4)
  cat(cmd, file=filename5)
  system("run.bat")
}

for (rmsea in c('r0','p05', 'p08', 'p12')){
  for (n in c('n100', 'n200', 'n300', 'n400', 'n500')) {
    plotchi(rmsea, n)
  } 
}

plotchi('p05','n100')

## function to plot the density according the the results from FPDE
plot.res<-function(rmsea,n){
  filename1<-paste(rmsea,'-',n,'.txt', sep='')
  
  filename2<-paste(rmsea,'-',n,'.R', sep='')
  cat('library(stats)\n', file=filename2)
  cat('data_in = scan("', filename1, '")\n', sep='', file=filename2, append=T)
  cat('seq_in = seq(min(data_in),max(data_in),length=1000)\n', file=filename2, append=T)
  cat('a = 0\n', file=filename2, append=T)
  cat('b = 1\n', file=filename2, append=T)
  cat('sdata_in = sort(data_in)\n', file=filename2, append=T)
  cat('nlength = length(data_in)\n', file=filename2, append=T)
  cat('pdata_in = seq(from=1, to=nlength, by=1) / nlength - 0.5 / nlength\n', file=filename2, append=T)
  cat('tdata_in=qchisq(pdata_in, df=    0.3000000000E+01)\n', file=filename2, append=T)

  filename3<-paste(rmsea,'-',n,'.out', sep='')

  temp<-readLines(filename3)

  nl<-grep("========================================= stage k=2 ==============================================", temp)

  for (i in (nl[1]+2):(nl[1]+6)){
    cat(temp[i],'\n',file=filename2, append=T)
  }

  ## Now write the codes for plot
  filename4<-paste(rmsea,'-',n,'.pdf', sep='')
  cat('pdf(file="', filename4, '")\n', sep='', file=filename2, append=T)

  ee<-as.numeric(substr(rmsea,2,3))/100
  nn<-as.numeric(substr(n,2,4))
  ncp<-ee^2*3*nn*2
  
  lmt<-max(dchisq(1:100, 3, ncp))+.01

  cat('plot(seq_in,den_2  ,type="l", lwd=2, main="",xlab="chi-square", ylab="density", ylim=c(0,', lmt,'))\n',file=filename2, append=T)
  cat('hist(data_in,freq=F, add=T)\n',file=filename2, append=T)

  cat('curve(dchisq(x, df=3,ncp=', ncp, '),add=T, lty=2, lwd=2)\n',file=filename2, append=T)
  cat('legend("topright", c("Approximation","Observed"), lty=c(2,1), lwd=c(2,2))\n',file=filename2, append=T)
  cat('dev.off()\n', sep='', file=filename2, append=T)
  
  source(filename2)
}

for (rmsea in c('r0','p05', 'p08', 'p12')){
  for (n in c('n100', 'n200', 'n300', 'n400', 'n500')) {
    plot.res(rmsea, n)
  } 
}

plot.res('p08','n200')