common mpsf,nx,ny,np,psfm,wp,xp,yp,gp
common comp,r,x1,y1,t,cx,cy,ax,ay,fo,sa,clx,cly,sax,say,ashx,ashy,tx,ty
common done,done,crve,iscr
common numb,n
common filen,fname,logfil
common spher,csph
common camer,camera,lam
;
done=0 & crve=1 & iscr=1
;
print,'filename for data,fname' & fname=' '&   read,fname
logfil=3 & openw,logfil,fname+'.log'
print,' Input 0 for batch and 1 for interactive input' &inter=0 &read,inter
;
n=256&camera='WFFULL'&lam=0.368&csph=0.0
if(inter eq 1) then begin
   print,' Input n, lam csph crve (eg 256 0.547 0.0 1)'
   read,n,lam,csph,crve 
   infil=4 & openw,infil,fname+'.in'
   printf,infil,n,lam,csph,crve 
   print,'Input camera (WF/PC/FC + FULL/SEC)' & read,camera
   printf,infil,camera
   scale,m,pm               
   print,'Oversampling :',m,'Pupil Magnification :',pm,' Maximum psf size :',n/m
endif else begin
   infil=4 & openr,infil,fname+'.in'
   readf,infil,n,lam,csph,crve 
   readf,infil,camera
   scale,m,pm
   print,'Oversampling :',m,'Pupil Magnification :',pm,' Maximum psf size :',n/m
endelse
print,' n,camera,csph,crve,iscr,lam=',n,camera,csph,crve,iscr,lam
printf,logfil,' n,camera,csph,crve,iscr,lam=',n,camera,csph,crve,iscr,lam
;
getpsf,inter
;       Set starting guess
if(inter eq 1) then begin   
   seta,ac
   inloop: printa,ac,0
   print,'OK? (1-->Y)' & iok=1 &read,iok
   if(iok ne 1) then begin
      print,' input index and correction'
      index=0 & correct=0. & read,index,correct & ac(index)=correct
      goto , inloop
   endif
   printf,infil,'$(i3)',n_elements(ac)
   for i=0,n_elements(ac)-1 do begin
      temp=ac(i)*lam/0.547
      printf,infil,temp
   endfor
endif else begin
   nc=0& readf,infil,nc
   ac=fltarr(nc)
   for i=0,n_elements(ac)-1 do begin
      readf,infil,temp
      ac(i)=temp*0.547/lam
   endfor
endelse
;       Print it out
print,' Starting position:'
printf,logfil,' Starting position:'
printa,ac,1
nc=n_elements(ac)
;
;       Curvefit
;
;for icsph=0,5 do begin
;csph=-0.5176+icsph*0.1
if(crve) then begin
   if(nx*ny*np lt 0)then nx=long(nx) & y=fltarr(nx*ny*np)
   for ip=0,np-1 do for j=0,ny-1 do y(j*nx+ip*nx*ny)=psfm(*,j,ip)
   x=fltarr(nx*ny*np) & w=x
   for ip=0,np-1 do w(ip*nx*ny)=replicate(wp(ip),nx*ny)
   w=w/(y+0.03*max(y))
   sigmaa=replicate(1.,nc) & for ip=0,2*np-1 do sigmaa(ip)=0
   for ip=2*np,nc-1 do sigmaa(ip)=1 & for ip=0,2*np-1 do sigmaa(ip)=1
   print,'Fitting everything ...csph=',csph
   yfit=CRVEFIT(X,Y,W,Ac,SIGMAA)
;
;       Simplex
;
endif else begin
   openw,68,'fitbsimp.out'
   da=replicate(0.1,n_elements(ac))
   nloops=2000
   simplex,ac,da,nloops
endelse
;       Print result out
print,' Ending position:',csph
printf,logfil,' Ending position:',csph
printa,ac,1
w=166/nx & psprint,max([np*n,w*nx*np]),n+3*w*ny,fname
if(np ge 2) then begin
   psfc=fltarr(nx,ny,np)
   for ip=0,np-1 do for j=0,ny-1 do $
       psfc(0,j,ip)=yfit(j*ny+ip*nx*ny:(j+1)*ny-1+ip*nx*ny)
endif else begin
   psfc=fltarr(nx,ny)
   for j=0,ny-1 do psfc(0,j)=yfit(j*ny:(j+1)*ny-1)
endelse
print,'best fit is in psfc, measurement is in psfm'
;endfor
stop
end