;FIT_CORRELATIONS
;
;PURPOSE:
;	To try and fit the autocorrelation of the ShackHartman grid
;
;INPUT:
;	The function fit_function.pro is used


pro fit_simul_wind_correlations,Date=Date,stats_directory=stats_directory,$
        x_center_file=x_center_file,y_center_file=y_center_file,$
        depth=depth,x_off=x_off,y_off=y_off


If (not keyword_set(Date)) then Date='10'
If (not keyword_set(stats_directory)) then stats_directory='stats_1_31'
If (not keyword_set(postscript_dir)) then postscript_dir='distorted_postscripts'

;This is the factor that converts between an pixel spacing on the ShackHartmann and an
;angle on the primary mirror
Ang_factor=float(((24e-6/60e-3)*(3.5e-3))^2)
stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/'+stats_directory+'/'
postscript_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/'+postscript_dir+'/'
cor_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$
          '/ShackHartman/correlations/'

common func_simul_wind_xy, dis, xtrans, xlong, ytrans, ylong

center_sub=1
if center_sub eq 0 then begin
readcol,stats_directory+'Depth_0Correlations.txt',$
	xlong,xtrans,ylong,ytrans,format='(f,f,f,f)'
endif else begin
readcol,stats_directory+'Depth_0_Cen_Sub_Correlations.txt',$
        xlong,xtrans,ylong,ytrans,format='(f,f,f,f)'
endelse
dis=findgen(25)*.17
weight=[.03,.03,.03,.03,.03,.03,.03,.03,.03,.03,$
	.03,.03,.03,.03,.03,.03,.03,.03,.03,.03,$
	.03,.03,.03,.03]

xtrans=Ang_factor*xtrans
xtranskeep=xtrans
xlong=Ang_factor*xlong
xlongkeep=xlong
ytrans=Ang_factor*ytrans
ytranskeep=ytrans
ylong=Ang_factor*ylong
ylongkeep=ylong     
rs=findgen(6)
                                                                                          
;The function we are trying to fit is
;y=((d+(k0/kl)^(1/(b-4))*D)/po))^(b-2)
;
;p(0)=ko
;p(1)=wind
;p(2)=kl
;p(3)=kt
;p(4)=b
;p(5)=po
 
;Minimize the difference
if date eq '10' or date eq '12' then begin
rs=amoeba(1.0e-5,scale=[0.1,0.4,.2,.3,.3,.2],p0=[0.2,1,.2,.1,11./3.,.2],$
	function_value=fval,function_name='fit_simul_wind_function',ncalls=tries,$
	nmax=1000000)
endif else begin
rs=amoeba(1.0e-6,scale=[0.2,0.1,.1,.3,.3,.2],p0=[0.1,1,.1,.3,11./3.,.3],$
        function_value=fval,function_name='fit_simul_wind_function',ncalls=tries,$
        nmax=1000000)
endelse

a=float(1./(rs(4)-4))
b=float(rs(4)-2)
koxl=float(rs(0)/rs(2))
koyl=float(rs(0)/rs(2))
koxt=float(rs(0)/rs(3))
koyt=float(rs(0)/rs(3))

fitxt=fltarr(N_elements(dis))
fitxt=rs(3)*(((dis+(koxt^a)*.17)/rs(5))^b*(5e-7/(dis+(koxt^a)*.17))^2)
fitxl=fltarr(N_elements(dis))
fitxl=rs(1)*rs(2)*(((dis+(koxl^a)*.17)/rs(5))^b*(5e-7/(dis+(koxl^a)*.17))^2)
fityt=fltarr(N_elements(dis))
fityt=rs(1)*rs(3)*(((dis+(koyt^a)*.17)/rs(5))^b*(5e-7/(dis+(koyt^a)*.17))^2)
fityl=fltarr(N_elements(dis))
fityl=rs(2)*(((dis+(koyl^a)*.17)/rs(5))^b*(5e-7/(dis+(koyl^a)*.17))^2)


set_plot,'ps'
if Center_sub eq 0 then begin
  device,filename=Postscript_dir+'Correlations_Fitted_Wind_Simul.ps'
endif else begin
  device,filename=Postscript_dir+'Correlations_Cen_Sub_Fitted_Wind_Simul.ps'
endelse                                                            
loadct,39                       ;load color table 39
device,/color                   ;allow color on the postscript
device,ysize=8.5,/inches        ;Height of plot in y
device,xsize=10.0,/inches        ;Width of plot in x
device,yoffset=1.0,/inches      ;Y position of lower left corner
                                                                                
white='FFFFFF'x
black='000000'x
!P.CHARSIZE=1.2
!P.THICK=4.

sym_arr=[2,0]
items=['Data Points','Fit']

;**************PLOT XTRANS AND FIT
plot,fitxt,psym=0,title="Fit for Transverse X Covariance (5-"+Date+')',$
	xtitle="n",$
	ytitle="Radians Squared"

oplot,xtranskeep,psym=1

legend,items,psym=sym_arr
xyouts,.7,.75,'W='+strtrim(string(rs(1)),2),/normal
xyouts,.7,.7,'k!l0!N='+strtrim(string(rs(0)),2),/normal
xyouts,.7,.65,'k!lT!N='+strtrim(string(rs(3)),2),/normal
xyouts,.7,.6,'r!lo!N='+strtrim(string(rs(5)),2),/normal
xyouts,.7,.55,'!4b!6=='+strtrim(string(rs(4)),2),/normal

;*************PLOT XLONG AND FIT
plot,fitxl,psym=0,title="Fit for Longitudinal X Covariance (5-"+Date+')',$
	 xtitle="n",$
         ytitle="Radians Squared"

oplot,xlongkeep,psym=1
legend,items,psym=sym_arr
xyouts,.7,.75,'W='+strtrim(string(rs(1)),2),/normal
xyouts,.7,.7,'k!l0!N='+strtrim(string(rs(0)),2),/normal
xyouts,.7,.65,'k!lL!N='+strtrim(string(rs(2)),2),/normal
xyouts,.7,.6,'r!lo!N='+strtrim(string(rs(5)),2),/normal
xyouts,.7,.55,'!4b!6='+strtrim(string(rs(4)),2),/normal

;*************PLOT YTRANS AND FIT
plot,fityt,psym=0,title="Fit for Transverse Y Covariance (5-"+Date+')',$
	 xtitle="n",$
         ytitle="Radians Squared"

oplot,ytranskeep,psym=1
                                                                                        
legend,items,psym=sym_arr
xyouts,.7,.75,'W='+strtrim(string(rs(1)),2),/normal
xyouts,.7,.7,'k!l0!N='+strtrim(string(rs(0)),2),/normal
xyouts,.7,.65,'k!lT!N='+strtrim(string(rs(3)),2),/normal
xyouts,.7,.6,'r!lo!N='+strtrim(string(rs(5)),2),/normal
xyouts,.7,.55,'!4b!6='+strtrim(string(rs(4)),2),/normal

;*************PLOT YTRANS AND FIT
plot,fityl,psym=0,title="Fit for Longitudinal Y Covariance (5-"+Date+')',$
	 xtitle="n",$
         ytitle="Radians Squared"

oplot,ylongkeep,psym=1
                                                                                        
legend,items,psym=sym_arr
xyouts,.7,.75,'W='+strtrim(string(rs(1)),2),/normal
xyouts,.7,.7,'k!l0!N='+strtrim(string(rs(0)),2),/normal
xyouts,.7,.65,'k!lL!N='+strtrim(string(rs(2)),2),/normal
xyouts,.7,.6,'r!lo!N='+strtrim(string(rs(5)),2),/normal
xyouts,.7,.55,'!4b!6='+strtrim(string(rs(4)),2),/normal


device,/close
set_plot,'X'

set_plot,'z'
erase
device, set_font='Courier'
device,set_resolution=[700,600]
!p.charsize=.8
!p.charthick=1.2
!x.thick=2
!y.thick=2
!p.thick=1
!p.noerase=0
red='FF0000'x
colors=[70,210,140,0]
sym_arr=[2,5,4,1]
items=['X Longitudinal','Y Longitudinal','Y Transverse','X Transverse']
plot,fitxt,psym=0,title="Fits for 5/"+Date+'/2005',$
        ytitle='Radians Squared',xtitle="n",$
        yrange=[0,2.5e-12],$
        position=[0.15,.1,.9,.9],/normal,$
        background=white,color=black
oplot,xtranskeep,psym=1,color=0         
oplot,fitxl,psym=0,color=black                          
oplot,xlongkeep,psym=2,color=70
oplot,fityt,psym=0,color=black         
oplot,ytranskeep,psym=4,color=140                  
oplot,fityl,psym=0,color=black 
oplot,ylongkeep,psym=5,color=210
legend,items,color=colors,psym=sym_arr,textcolors=[0,0,0,0]
jpgimg=tvrd()
tvlct,reds,greens,blues,/get
if Center_sub eq 0 then begin
   write_png,'/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_report_figs/Wind_No_Dainty_Correlations_5-'+Date+'-2005.png',jpgimg,reds,greens,blues
endif else begin
   write_png,'/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_report_figs/Wind_No_Dainty_Correlations_Cen_Sub_5-'+Date+'-2005.png',jpgimg,reds,greens,blues
endelse
stop

end