character*80 line
      character*1 all
      character*40 imgname
      character*2880 head1,head2,head3
      call getarg(1,imgname)
      call getarg(1,all)
      open(19,file=imgname)
      read(19,'(3a)') head1,head2,head3
      close(19)
      kh=index(head3,'FOCUSVAL')
      read(head3(kh+11:),'(3(1x,f5.3))') fa0,fb0,fc0
c     FOCUSVAL= '*0.581*0.864*0.844' / Focus 

      k=index(imgname,'.fit')
      call system('fkey '//imgname//' '//imgname(k-4:k-1))

      call focus(imgname(k-4:k-1)//'_1',all)
      call focus(imgname(k-4:k-1)//'_2',all)
      call focus(imgname(k-4:k-1)//'_3',all)
      call focus(imgname(k-4:k-1)//'_4',all)


c      call system('rm -f '//imgname(1:k-1)//'_?.fit')

c     convert equation focus plan to filed center to be 0,0
c    CCD1   CCD3
c        0,0
c    CCD2   CCD4

C horizon gap: 5.5mm = 5500um/15um = 367 pixels
c vertical gap: 2.0mm = 2000um/15um = 133 pixels
c  then
      dy=366.6/2.0
      dx=133.3/2.0

      k=index(imgname,'.fit')
      open(19,file=imgname(k-4:k-1)//'_1-0.bat',status='old')
      read(19,*)
      read(19,*)
      read(19,'(a)') line
      k=index(line,'*x')
      read(line(k-8:k-1),*) aa1
      k=index(line,'*y')
      read(line(k-8:k-1),*) bb1
      read(line(K+3:k+13),*) cc1
      write(*,*)
      write(*,*)
      write(*,*)

      write(*,'(2f9.5,f10.5)') aa1,bb1,cc1

c    z=aa1*x+bb1*y+cc1
c
c      N
c    E   W
c      S
c Then set S-W corner to be 0,0

c===>
c      xx=x-4096-dx
c      yy=y+dy
c => 
c   x=xx+4096+dx
c   y=yy-dy
c   z=aa1*(4096+xx+dx)+bb1*(yy-dy)+cc1
c   z=aa1*xx+aa1*(4096+dx)+bb1*yy-bb1*dy+cc1
c   z=aaa1*xx+bbb1*yy+ccc1
c     aaa1=aa1; bbb1=bb1; ccc1=cc1+aa1*(4096+dx)-bb1*dy

      aaa1=aa1
      bbb1=bb1
      ccc1=cc1+aa1*(4096+dx)-bb1*dy

      close(19)

      k=index(imgname,'.fit')
      open(19,file=imgname(k-4:k-1)//'_2-0.bat',status='old')
      read(19,*)
      read(19,*)
      read(19,'(a)') line
      k=index(line,'*x')
      read(line(k-8:k-1),*) aa2
      k=index(line,'*y')
      read(line(k-8:k-1),*) bb2
      read(line(K+3:k+13),*) cc2
      write(*,'(2f9.5,f10.5)') aa2,bb2,cc2
      close(19)

c    z=aa2*x+bb2*y+cc2
c
c      N
c    E   W
c      S
c
c Then set N-W corner to be 0,0
c
c===>
c      xx=x-4096-dx
c      yy=y-4032-dy
c => 
c   x=4096+dx+xx
c   y=4032+dy+yy
c   z=aa2(4096+dx+xx)+bb2(4032+dy+yy)+cc2
c   z=aa2*xx+bb2*yy+aa2*(4096+dx)+bb2*(4032+dy)+cc2
c   z=aaa2*xx+bbb2*yy+ccc2
c     aaa2=aa2; bbb2=bb2; ccc2=cc2+aa2*(4096+dx)+bb2*(4032+dy)

      aaa2=aa2
      bbb2=bb2
      ccc2=cc2+aa2*(4096+dx)+bb2*(4032+dy)



      k=index(imgname,'.fit')
      open(19,file=imgname(k-4:k-1)//'_3-0.bat',status='old')
      read(19,*)
      read(19,*)
      read(19,'(a)') line
      k=index(line,'*x')
      read(line(k-8:k-1),*) aa3
      k=index(line,'*y')
      read(line(k-8:k-1),*) bb3
      read(line(K+3:k+13),*) cc3
      write(*,'(2f9.5,f10.5)') aa3,bb3,cc3
      close(19)

c    z=aa3*x+bb3*y+cc3
c
c      N
c    E   W
c      S
c
c S-E corner is 0,0

c===>
c      xx=x+dx
c      yy=y+dy
c => 
c   x=yy-dx
c   y=xx-dy
c   z=aa3*(xx-dx)+bb3*(yy-dy)+cc3
c   z=aa3*xx+bb3*yy+cc3-aa3*dx-bb3*dy
c   z=aaa3*xx+bbb3*yy+ccc3
c   aaa3=bb3; bbb3=aa3; ccc3=cc3-aa3*dx-bb3*dy
      aaa3=aa3
      bbb3=bb3
      ccc3=cc3-aa3*dx-bb3*dy


      k=index(imgname,'.fit')
      open(19,file=imgname(k-4:k-1)//'_4-0.bat',status='old')
      read(19,*)
      read(19,*)
      read(19,'(a)') line
      k=index(line,'*x')
      read(line(k-8:k-1),*) aa4
      k=index(line,'*y')
      read(line(k-8:k-1),*) bb4
      read(line(K+3:k+13),*) cc4
      write(*,'(2f9.5,f10.5)') aa4,bb4,cc4
      close(19)
c    z=aa4*x+bb4*y+cc4
c
c      N
c    E   W
c      S
c
c Then set N-E corner to be 0,0
c
c===>
c      xx=x+dx
c      yy=y-4032-dy
c => 
c   x=-xx-dx
c   y=4032+dy+yy
c   z=aa4*(xx-dx)+bb4(4032+dy+yy)+cc4
c   z=aa4*xx+bb4*yy+bb4*4032+cc4-aa4*dx+bb4*dy
c   z=aaa4*xx+bbb4*yy+ccc4
c     aaa4=aa4; bbb4=bb4; ccc4=bb4*4032+cc4-aa4*dx+bb4*dy

      aaa4=aa4
      bbb4=bb4
      ccc4=bb4*4032+cc4-aa4*dx+bb4*dy

      write(*,*) '=============================================='
      write(*,*) "In common coordinate system with 0,0 in field center"
      write(*,'(2f9.5,f10.5)') aaa1,bbb1,ccc1
      write(*,'(2f9.5,f10.5)') aaa2,bbb2,ccc2
      write(*,'(2f9.5,f10.5)') aaa3,bbb3,ccc3
      write(*,'(2f9.5,f10.5)') aaa4,bbb4,ccc4
      write(*,*) '=============================================='

      aaa=(aaa1+aaa2+aaa3+aaa4)/4.0
      bbb=(bbb1+bbb2+bbb3+bbb4)/4.0
      ccc=(ccc1+ccc2+ccc3+ccc4)/4.0
      write(*,*) 'The mean plan is:'
      write(*,'(f9.5,a3,f9.5,a3,f10.5)') aaa,'*x+',bbb,'*y+',ccc-30.0

c The coordinates of the 3 actuators are:
c                           x     y
c A    333.3mm,  32.05mm  22220.0   2136.7
c B  -138.44mm,-305.26mm  -9229.3 -20350.7
c C  -194.45mm, 272.40mm -12963.3  18160.0

c 333300/15
c  32050/15
c -138440/15
c -305260/15
c -194450/15
c  272400/15

      fa= 22220.0*aaa+2136.7*bbb+ccc-30.0
      fb=-9229.3*aaa-20350.7*bbb+ccc-30.0
      fc=-12963.3*aaa+18160.0*bbb+ccc-30.0

      write(*,'(a,5f6.3)') 'present focus:',fa0,fb0,fc0
      write(*,*) 'Actuator A =',fa
      write(*,*) 'Actuator B =',fb
      write(*,*) 'Actuator C =',fc
      write(*,'(a10,3f10.3)') imgname(1:10),fa,fb,fc
c clearing
      fa1=fa0-fa*0.0005
      fb1=fb0-fb*0.0005
      fc1=fc0-fc*0.0005
      write(*,'(a,5f6.3)') 'Expect focus:',fa1,fb1,fc1
      open(20,file=imgname(1:10)//'.txt')
      write(20,'(4f6.0,9f6.3)') ccc-30,fa,fb,fc,fa0,fb0,fc0,fa1,fb1,fc1,
     +     fa1/fa1,fb1/fa1,fc1/fa1
      close(20)

      call system('rm -f *_?.*')
      call system('rm -f *-?.bat')
      call system('rm -f test.cat')

      stop
      end

      subroutine focus(logname,all)
      character*1 all
      parameter (nstar=2000)
      parameter (np=20)
      real x(nstar,np),y(nstar,np),see(nstar,np)
      real x0(nstar),y0(nstar),z0(nstar),zz(nstar),xx(np),yy(np)
      real x1(nstar),y1(nstar)
      integer nsp(nstar),ns
      character*80 line,logname,outname
      character*1 note

      note='n'

      outname=logname(1:6)//'.dat'
      open(19,file=logname(1:6),status='old')
      i=0
      
40    read(19,'(a)',end=50) line
      if(line(1:2).eq.'##') goto 40
      if(line(1:1).eq.'#') then
         if(nsp(i).lt.7) i=i-1
         j=0         
         i=i+1
         read(line(2:),*) x0(i),y0(i)
         read(19,*)
      endif

      if(line(1:1).eq.'#') then
         goto 40
      else
         j=j+1
         nsp(i)=j
      endif
      k=8-j
      read(line,*) x(i,k),y(i,k),adu,sky,see(i,k),ring
      goto 40
 50   close(19)
      ns=i
      write(*,*) 'number of stars =',ns
      nstep=1
      if(all.ne.'!'.and.ns.gt.20) nstep=ns/20
      write(*,*) 'number of stars select =',ns/nstep
      open(20,file=outname)
      ii=0
      do i=1,ns,nstep
         vmin=999999.0
         do j=1,nsp(i)
            yy(j)=see(i,j)
            vmin=min(vmin,yy(j))
            if(vmin.eq.yy(j)) nn=j
         enddo
c         write(*,*) 'vmin=',vmin

c         s(1)=1.97
c         s(2)=1.84
c         s(3)=1.60
c         s(4)=1.69
c         s(5)=1.90
c         s(6)=2.17
c         s(7)=2.74
         vmin=9999999.0
         do nc=40,250
            ccc=nc/50.0
            do nb=-40,80
               bbb=nb/10.0
               do na=1,30
                  aaa=na/50.0
                  sum=0
                  do j=1,nsp(i)
                     sum=sum+(see(i,j)-aaa*(j-bbb)**2.0-ccc)**2.0
                  enddo
                  vmin=min(vmin,sum)
                  if(vmin.eq.sum) then
                     aaaa=aaa
                     bbbb=bbb
                     cccc=ccc
                  endif
               enddo
            enddo
         enddo

         z0(i)=bbbb*30.0
c         write(*,'(2f9.3,3f9.5)') x0(i),y0(i),aaaa,bbbb,cccc
         if(z0(i).gt.-50.0.and.z0(i).lt.260.0) then
            ii=ii+1
            x0(ii)=x0(i)
            y0(ii)=y0(i)
            z0(ii)=z0(i)
            write(20,'(2f7.1,f9.3)') x0(ii),y0(ii),z0(ii)
            x1(ii)=(x0(ii)-2048)**2.0/2048
            y1(ii)=(y0(ii)-2048)**2.0/2048
         endif
      enddo
      close(20)

c fitting the telescope focus with flat plan
      ns=ii
       nn=ns
       a=0
       b=0
       do j=1,10
          do i=1,nn
             zz(i)=z0(i)-b*y0(i)
          enddo
          call medfit(x0,zz,nn,aa,bb,abdev)
          a=bb
          do i=1,nn
             zz(i)=z0(i)-a*x0(i)
          enddo
          call medfit(y0,zz,nn,aa,bb,abdev)
          b=bb
       enddo
       write(*,*) '-------------------------------'
       write(*,'(2(a5,f8.5),a5,f10.5)') 'z = ',a,' x + ',b,' y + ',aa
       write(*,'(a2,f9.5)') '+-',abdev
       write(*,*) '-------------------------------'

c if note ne 'o' then do not give any output for display 
       open(20,file=logname(1:6)//'-0.bat')
       write(20,*) 'set terminal gif'
       write(20,*) 'set output "'//logname(1:kk-1)//'-flat.gif"'
       write(20,'(2(a,f8.5),a,f10.5a)') 'splot [0:4096][0:4096]'//
     +      '[00:200] ',
     +      a,'*x+',b,'*y+',aa,' , "'//outname(1:kk+3)//'"'
       close(20)
       if(note.ne.'o') goto 80

       write(*,'(2(a,f8.5),a,f10.5,a)') 'splot [0:4096][0:4096]'//
     +      '[00:200] ',
     +      a,'*x+',b,'*y+',aa,' , "'//outname(1:kk+3)//'"'
       call system('gnuplot '//logname(1:kk-1)//'-0.bat')

       open(20,file=logname(1:kk-1)//'-1.bat')
       write(20,*) 'set terminal gif'
       do i=1,360,5
          ir=i+30
          if(ir.ge.360) ir=ir-360
          write(20,'(a,i3)') 'set view ,',ir
          ii=1000+i
          write(20,'(a,i4,a)') 'set output "focus',ii,'.gif"'
         
       write(20,'(2(a,f8.5),a,f9.5,a)') 'splot [0:4596][0:4096]'//
     +         '[00:200]',a,
     +      '*x+',b,'*y+',aa,' , "'//outname(1:kk+3)//'"'
       enddo
c 'splot [0:4096] [0:4096]  'a,'*x+',b,'*y+',aa,' , "'//outname(1:kk+3)//'"'
       close(20)
       call system('gnuplot '//logname(1:kk-1)//'-1.bat')
       CALL system('convert -adjoin -delay 10 -loop 0 
     +      focus1???.gif '//logname(1:kk-1)//'-anime.gif')
       call system('/bin/rm focus1???.gif')

 

c     fitting the telescope focus surface with z = a x**2 + b y**2 + c
       nn=ns
       do j=1,10
          do i=1,nn
             zz(i)=z0(i)-b*y1(i)
          enddo
          call medfit(x1,zz,nn,aa,bb,abdev)
          a=bb
          do i=1,nn
             zz(i)=z0(i)-a*x1(i)
          enddo
          call medfit(y1,zz,nn,aa,bb,abdev)
          b=bb
       enddo
       write(*,*) '-------------------------------'
       write(*,'(2(a5,f8.5),a5,f10.5)') 'z = ',a,' x*x  ',b,
     +      ' y*y + ',aa
       write(*,*) 'aa=',aa
       write(*,'(a2,f9.5)') '+-',abdev
       write(*,*) '-------------------------------'
 80    continue
       return
      stop
      end

 
C  (C) Copr. 1986-92 Numerical Recipes Software ,2:-5K#R..
      SUBROUTINE medfit(x,y,ndata,a,b,abdev)
      INTEGER ndata,NMAX,ndatat
      PARAMETER (NMAX=20100)
      REAL a,abdev,b,x(ndata),y(ndata),arr(NMAX),xt(NMAX),yt(NMAX),aa,
     *abdevt
      COMMON /arrays/ xt,yt,arr,aa,abdevt,ndatat
CU    USES rofunc
      INTEGER j
      REAL b1,b2,bb,chisq,del,f,f1,f2,sigb,sx,sxx,sxy,sy,rofunc
      sx=0.
      sy=0.
      sxy=0.
      sxx=0.
      do 11 j=1,ndata
        xt(j)=x(j)
        yt(j)=y(j)
        sx=sx+x(j)
        sy=sy+y(j)
        sxy=sxy+x(j)*y(j)
        sxx=sxx+x(j)**2
11    continue
      ndatat=ndata
      del=ndata*sxx-sx**2
      aa=(sxx*sy-sx*sxy)/del
      bb=(ndata*sxy-sx*sy)/del
      chisq=0.
      do 12 j=1,ndata
        chisq=chisq+(y(j)-(aa+bb*x(j)))**2
12    continue
      sigb=sqrt(chisq/del)
      b1=bb
      f1=rofunc(b1)
      b2=bb+sign(3.*sigb,f1)
      f2=rofunc(b2)
1     if(f1*f2.gt.0.)then
       bb=2.*b2-b1
        b1=b2
        f1=f2
        b2=bb
        f2=rofunc(b2)
        goto 1
      endif
      sigb=0.01*sigb
2     if(abs(b2-b1).gt.sigb)then
        bb=0.5*(b1+b2)
        if(bb.eq.b1.or.bb.eq.b2)goto 3
        f=rofunc(bb)
        if(f*f1.ge.0.)then
          f1=f
          b1=bb
        else
          f2=f
          b2=bb
        endif
        goto 2
      endif
3     a=aa
      b=bb
      abdev=abdevt/ndata
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software ,2:-5K#R..
      FUNCTION rofunc(b)
      INTEGER NMAX
      REAL rofunc,b,EPS
      PARAMETER (NMAX=20100,EPS=1.e-7)
CU    USES select
      INTEGER j,ndata
      REAL aa,abdev,d,sum,arr(NMAX),x(NMAX),y(NMAX),select
      COMMON /arrays/ x,y,arr,aa,abdev,ndata
      do 11 j=1,ndata
        arr(j)=y(j)-b*x(j)
11    continue
      if (mod(ndata,2).eq.0) then
        j=ndata/2
        aa=0.5*(select(j,ndata,arr)+select(j+1,ndata,arr))
      else
        aa=select((ndata+1)/2,ndata,arr)
      endif
      sum=0.
      abdev=0.
      do 12 j=1,ndata
        d=y(j)-(b*x(j)+aa)
        abdev=abdev+abs(d)
        if (y(j).ne.0.) d=d/abs(y(j))
        if (abs(d).gt.EPS) sum=sum+x(j)*sign(1.0,d)
12    continue
      rofunc=sum
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software ,2:-5K#R..
      FUNCTION select(k,n,arr)
      INTEGER k,n
      REAL select,arr(n)
      INTEGER i,ir,j,l,mid
      REAL a,temp
      l=1
      ir=n
1     if(ir-l.le.1)then
        if(ir-l.eq.1)then
          if(arr(ir).lt.arr(l))then
            temp=arr(l)
            arr(l)=arr(ir)
            arr(ir)=temp
          endif
        endif
        select=arr(k)
        return
      else
        mid=(l+ir)/2
        temp=arr(mid)
        arr(mid)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        if(j.ge.k)ir=j-1
        if(j.le.k)l=i
      endif
      goto 1
      
      END