pro fun,x,a,f
common comp,r,x1,y1,t,cx,cy,ax,ay,fo,sa,clx,cly,sax,say,ashx,ashy,tx,ty
common mpsf,nx,ny,np,psfm,wp,xp,yp,gp
common done,done,crve,iscr
common numb,n
common filen,fname,logfil
common prtt,wave,mfc
common iter,titer
common spher,csph
common camer,camera,lam
w=166/nx   
o=2*np-2         
scale,m,pm
if(not done) then begin
   print,'************** n=',n
   done=1 
   nscl=1. & if(n gt 256) then nscl=n/256
   if(iscr)then tvwindow,0,xsize=max([np*n/nscl,w*nx*np]),ysize=n/nscl+3*w*ny
   print,'First time caller'
   x1=indgen(n)#replicate(1.,n)
   y1=replicate(1.,n)#indgen(n)
   x1=(x1-n/2)*(4./n)/pm
   y1=(y1-n/2)*(4./n)/pm
   r=sqrt(x1^2+y1^2)
   t=atan(y1,x1)
   tx =1.8992573*r*cos(t)
   ty =1.8992573*r*sin(t)
   fo =3.8874443*(r^2-0.554450)
   ax =2.3137662*r^2*cos(2*t)
;   ax =74.782446*(r^6-1.663350*r^4+0.80313*r^2 -0.104406)
   ay =2.3137662*r^2*sin(2*t)     
   cx =8.3345629*(r^3-0.673796*r)*cos(t)
   cy =8.3345629*(r^3-0.673796*r)*sin(t)
   sa =16.895979*(r^4-1.108900*r^2+0.241243)                  
   if(n_elements(a) ge 10+o) then begin
      clx=2.6701691*r^3*cos(3*t)
      cly=2.6701691*r^3*sin(3*t)
      sax=12.033645*(r^4-0.750864*r^2)*cos(2*t)
      say=12.033645*(r^4-0.750864*r^2)*sin(2*t)
      ashx=29.851527*r^4*cos(4*t)
      ashy=29.851527*r^4*sin(4*t)
   endif
endif
i=complex(0.,1.)
scl=max(psfm)
os=max([n,w*nx])
if(crve)then begin
   if(titer eq -1) then begin
      f=fltarr(nx*ny*np)
      mfc=fltarr(np)
   endif
   mf=0
endif else begin
   titer=-1
   mfc=fltarr(np)
   f=0
   mf=0
endelse
if((titer eq -1) or (titer ge 2*np)) then begin
   wave=a(2+o)*fo+a(3+o)*cx+a(4+o)*cy
   if(csph ne 0.) then begin
      if(n_elements(a) ge 8+o)  then begin
         wave=wave+a(5+o)*ax+a(6+o)*ay+csph*sa
         a(7+o)=csph
      endif else begin
         wave=wave+csph*sa
      endelse
   endif else begin            
      if(n_elements(a) ge 8+o)  then wave=wave+a(5+o)*ax+a(6+o)*ay+a(7+o)*sa
   endelse
   if(n_elements(a) ge 10+o) then wave=wave+a(8+o)*clx+a(9+o)*cly
   if(n_elements(a) ge 14+o) then wave=wave+$
                             a(10+o)*sax+a(11+o)*say+a(12+o)*ashx+a(13+o)*ashy
endif
if(titer eq -1) then begin
   print,'tx1,ty1,fo,cx,cy',a(0),a(1),a(2+o),a(3+o),a(4+o)
   if(n_elements(a) ge 8+o) then print,'ax,ay,sa',a(5+o),a(6+o),a(7+o)
   if(n_elements(a) ge 10+o) then print,'clx,cly',a(8+o),a(9+o) 
   if(n_elements(a) ge 14+o) then print,'sax,say,ashx,ashy',$
                                    a(10+o),a(11+o),a(12+o),a(13+o)
endif
for ip=0,np-1 do begin
   if((titer eq -1) or (titer ge 2*np) or (titer/2 eq ip)) then begin
      pupil=[0.]
      if(camera eq 'PCFULL')then begin
         sens=0.194*n*pm/(4.*sqrt(2)*400) & fx=xp(ip)-400 & fy=yp(ip)-400
         pupil=(r lt 1.) and (r gt 0.33) 
         pupil(where((abs(x1) lt 0.0105) or (abs(y1) lt 0.0105)))=0.
         pupil(where((shift(r,0.0,0.890*n*pm/4) lt 0.071) or $
                     (shift(r,0.760*n*pm/4,-0.463*n*pm/4) lt 0.071) or $
                     (shift(r,-0.760*n*pm/4,-0.463*n*pm/4) lt 0.071))) =0. 
         pupil=rot_int(pupil,90+90*gp(ip))
         pupil=pupil and (shift(r,fx*sens,fy*sens) gt 0.414)
         pupil(where((abs(x1-fx*sens*4/(n*pm)) lt 0.0205) or $
                    ((abs(y1-fy*sens*4/(n*pm)) lt 0.0205) and $
                     (x1-fx*sens*4/(n*pm) lt 0.))))=0
      endif
      if(camera eq 'PCSEC')then begin
         sens=0.194*n*pm/(4.*sqrt(2)*400) & fx=xp(ip)-400 & fy=yp(ip)-400
         pupil=(r lt 1.) and (r gt 0.33) and $
                (shift(r,fx*sens,fy*sens) gt 0.414)
      endif
      if(camera eq 'WFFULL')then begin
         sens=0.391*n*pm/(4.*sqrt(2)*400) & fxr=xp(ip)-400 & fyr=yp(ip)-400
         fx=fxr/sqrt(2)-fyr/sqrt(2) & fy=fxr/sqrt(2)+fyr/sqrt(2) 
         pupil=(r lt 1.) and (r gt 0.33) 
         pupil(where((abs(x1) lt 0.0105) or (abs(y1) lt 0.0105)))=0.
         pupil(where((shift(r,0.0,0.890*n*pm/4) lt 0.071) or $
                     (shift(r,0.760*n*pm/4,-0.463*n*pm/4) lt 0.071) or $
                     (shift(r,-0.760*n*pm/4,-0.463*n*pm/4) lt 0.071))) =0. 
         pupil=rot_int(pupil,90+gp(ip)*90)
         pupil=pupil and (shift(r,fx*sens,fy*sens) gt 0.421)
         pupil(where((abs(x1-fx*sens*4/(n*pm)) lt 0.0205) or $
                    ((abs(y1-fy*sens*4/(n*pm)) lt 0.0205) and $
                     (x1-fx*sens*4/(n*pm) lt 0.))))=0
         pupil=rot_int(pupil,45)
      endif
      if(camera eq 'WFSEC')then begin
         sens=0.391*n*pm/(4.*sqrt(2)*400) & fx=xp(ip)-400 & fy=yp(ip)-400
         pupil=(r lt 1.) and (r gt 0.33) and $
                (shift(r,fx*sens,fy*sens) gt 0.421)
         pupil=rot_int(pupil,135)
      endif
      if(camera eq 'FCFULL') then begin
         pupil=(r lt 1.) and (r gt 0.33) 
         pupil(where((abs(x1) lt 0.011) or (abs(y1) lt 0.011)))=0.
         pupil(where((shift(r,0.0,0.889*n*pm/4) lt 0.063) or $
                     (shift(r,-0.826*n*pm/4,-0.340*n*pm/4) lt 0.063) or $
                     (shift(r,0.762*n*pm/4,-0.467*n*pm/4) lt 0.063))) =0. 
         pupil=rot_int(pupil,33-180.)
      endif
      if(camera eq 'FCSEC')then begin
         pupil=(r lt 1.) and (r gt 0.33) 
      endif
      if(camera eq 'F288') then begin
         pupil=(r lt 0.938) and (r gt 0.414) 
         pupil(where((abs(x1) lt 0.041) or (abs(y1) lt 0.041)))=0.
         pupil(where((shift(r,0.0,0.891*n*pm/4) lt 0.134) or $
                     (shift(r,-0.761*n*pm/4,-0.464*n*pm/4) lt 0.134) or $
                     (shift(r,0.760*n*pm/4,-0.466*n*pm/4) lt 0.134))) =0. 
         pupil=rot_int(pupil,33-180.)
      endif
      if(n_elements(pupil) eq 1) then begin
        print,' invalid camera type',camera
        stop
      endif     
      wav=(wave+a(0+2*ip)*tx+a(1+2*ip)*ty)*pupil
      psfc=shift(abs(fft(pupil*exp(2*!pi*i*wav),1))^2,n/2,n/2)
      if(m ge 2)then psfc=smooth(psfc,m)
      psfc=congrid(psfc(n/2-m*nx/2:n/2+m*nx/2,n/2-m*ny/2:n/2+m*ny/2),nx,ny)
      if(camera eq 'F288')then psfc=smooth(psfc,2)
      psfc=psfc/total(psfc) 
      psfa=psfc &for j=0,ny-1 do psfa(0,j)=psfm(*,j,ip)
      if(iscr and (titer eq -1))then begin
         tvscl,congrid(wav(n/4:3*n/4,n/4:3*n/4),n/nscl,n/nscl),ip*os,3*w*ny
         tv,congrid(psfa*255/scl,w*nx,w*ny),ip*os,2*w*ny
         tv,congrid(psfc*255/scl,w*nx,w*ny),ip*os,1*w*ny
         tv,congrid(abs(psfa-psfc)*255/scl,w*nx,w*ny),ip*os,0
      endif 
      dmf=total((psfc-psfa)^2*wp(ip)/(psfa+0.03*scl))
      if(titer eq -1) then begin
;         print,        'star, wmin, wmax, error',ip+1,min(wav),max(wav),dmf
;         printf,logfil,'star, wmin, wmax, error',ip+1,min(wav),max(wav),dmf
         mfc(ip)=dmf
      endif
      if(crve)then begin
         for j=0,ny-1 do f(j*nx+ip*nx*ny)=psfc(*,j)
         mf=mf+dmf
      endif else begin
         f=f+dmf
      endelse
   endif else begin
      if(crve)then begin
         mf=mf+mfc(ip)
      endif else begin
         f=f+mfc(ip)
      endelse
   endelse
endfor
if(titer eq -1) then begin
   if(crve)then begin
      print,'f=sqrt(total((psfc-psfa)^2)*W)',sqrt(mf)
      printf,logfil,'f=sqrt(total((psfc-psfa)^2)*W)',sqrt(mf)
      print,'M-components',mfc
      printf,logfil,'Merit function components',mfc
   endif else begin
      print,'f=sqrt(total((psfc-psfa)^2)*w)',sqrt(f)
      printf,logfil,'f=sqrt(total((psfc-psfa)^2)*w)',sqrt(f)
      print,'M-components',mfc
      printf,logfil,'Merit function components',mfc
   endelse
endif
return
end