c substract 2 file
	character*80 fa,fb,fc, ahead(72),bhead(72)
	real*8 a8(8),b8(8),bcoef(8)
	real a(2048,2048), b(2048,2048)

        real starx(400),stary(400),starv(400)
	real pz0(100),pz1(100)
        integer*2 px(100),py(100)

	logical logi
	real*8 cda,tda,sdb,cdb
	common /dec/cda,tda,sdb,cdb

c read a8, b8, bfile->a; xy--(a8)-->ad; ad--(b8)-->xy, a--xy-->b
c read         afile->a; find scale,
c do a - b

	if(iargc().lt.3)then
	  write(*,*)
	  write(*,*)'         ******** substarct 2048*2048 Pic. ***31217'
	  write(*,*)
	  write(*,*)'    hh1 file1.fit file2.fit out.fit'
	  write(*,*)
	  write(*,*)'    ex:  hh a b c'
	  write(*,*)
	  write(*,*)'    notice: program will call "sht coord"'
	  write(*,*)'            will call many files to work out SCALE'
	  write(*,*)
	  stop
	endif

	call getarg(1,fa)
        k=index(fa,'.')-1
	if(k.le.0)k=lnblnk(fa)
	fa=fa(1:k)//'.fit'
	inquire(file=fa,exist=logi)
        if(.not.logi)then
	  write(*,*)fa,' not found!'
	  stop
	endif
	call getarg(2,fb)
        k=index(fb,'.')-1
	if(k.le.0)k=lnblnk(fb)
	fb=fb(1:k)//'.fit'
	inquire(file=fb,exist=logi)
        if(.not.logi)then
	  write(*,*)fb,' not found!'
	  stop
	endif

	call gethead(fa,ahead)
	call geta8(ahead,a8)
        ipos=indexpos(ahead,'VOLT2   ')
        if(ahead(ipos)(42:46).ne."PIXEL")then
          write(*,*)fa(1:25)," do statsky first !"
          stop
        endif

        read(ahead(ipos)(22:),*)awhite

	call gethead(fb,bhead)
        ipos=indexpos(bhead,'VOLT2   ')
        if(bhead(ipos)(42:46).ne."PIXEL")then
          write(*,*)fb(1:25)," do statsky first !"
          stop
        endif

	call system("sht "//fa(1:lnblnk(fa)+1)//fb(1:lnblnk(fb))//" hh.tmp")
        fc="hh.tmp"	
	call gethead(fc,bhead)
	call geta8(bhead,b8)
        read(bhead(ipos)(22:),*)bwhite
	call xy2ad(b8,bcoef,a8)
	call readata(a,fc)
	call unlink("hh.tmp")
	write(*,*)"Wait"

        cda=dcos(a8(8))
        tda=dtan(a8(8))
        sdb=dsin(bcoef(8))
        cdb=dcos(bcoef(8))

	do 10 j=1,2048
	do 10 i=1,2048
	call xy_rad_xy(float(i),float(j),a8,bcoef,x,y)
	b(i,j)=-9999.
	if(x.lt.1.or.x.gt.2048.)goto 10
	if(y.lt.1.or.y.gt.2048.)goto 10
        iy0=y
        iy1=iy0+1
        wy=y-iy0
        zy=1.-wy
        ix0=x
        ix1=ix0+1
        wx=x-ix0
        zx=1.-wx
        x1=zx*zy
        x2=wx*zy
        x3=zx*wy
        x4=wx*wy
        b(i,j)=a(ix0,iy0)*x1+a(ix1,iy0)*x2+a(ix0,iy1)*x3+a(ix1,iy1)*x4
10	continue

111	call readata(a,fa)

c        nstar=400
c        call getstar(a,starx,stary,starv,nstar,awhite)
c        write(*,*)'nstar: ',nstar
c        k1=nstar/3.
c        k2=k1+99
c        if(k2.gt.nstar)k2=nstar
c        np=0
c        do 110 i=k1,k2
c        np=np+1
c        px(np)=starx(i)+0.5
c110     py(np)=stary(i)+0.5
c        call star100(a,awhite,px,py,pz0,np)
c        call star100(b,bwhite,px,py,pz1,np)
c        do 130 j=1,np
c130     pz1(j)=pz1(j)/pz0(j)
c        scale=xmedian(pz1,np)

	call system("j_find20 "//fa)
	call system("j_phot "//fa)
	call system("j_xy2048ad "//fa)
	call system("j_find20 "//fb)
	call system("j_phot "//fb)
	call system("j_xy2048ad "//fb)
	call system("all2048 "//fa//" "//fb)
	open(3,file='all2048.tmp',status='old')
	read(3,*)scale
	close(3)
        write(*,*)fb(1:22),'scale:',scale

	do 20 j=1,2048
	do 20 i=1,2048
	if(b(i,j).ne.-9999.)then
          a(i,j)=(a(i,j)-awhite)-(b(i,j)-bwhite)/scale
	else 
          a(i,j)=awhite
	endif
20	continue

	call getarg(3,fc)
        k=index(fc,'.')-1
	if(k.le.0)k=lnblnk(fc)
	fc=fc(1:k)//'.fit'
	write(*,*)' out file: '//fc(1:25)
	k=2048*2048*4+5760
	call swap4(a,2048*2048*4)
	open(50,file=fc,status='unknown',access='direct',recl=k)
	write(50,rec=1)ahead,a
	close(50)
	end

        subroutine readata(a,f1)
        integer c(720),a(1)
        character*80 f1,ch*2880
        equivalence (c,ch)
        nn=720         
        open(1,file=f1,status='old',access='direct',recl=nn*4)
        read(1,rec=1)ch
        n1=2048   
        n2=2048
        n=1   
55      k=index(ch,'END          ')
        if(k.eq.0)then
          n=n+1    
          read(1,rec=n)ch
          goto 55  
        endif   
        i=0
        k1=n+1
        k2=n+(n1*n2-1)/nn
        do 60 k=k1,k2
        if(mod(k,200).eq.0)call dispdot()
        read(1,rec=k)ch
        do 60 j=1,nn
        i=i+1
60      a(i)=c(j)
        close (1)
	call swap4(a,n1*n2*4)
        end
 
	subroutine gethead(f1,head)
	character*80 f1,head(72)
	open(52,file=f1,status='old',access='direct',recl=5760)
	read(52,rec=1)head
	close (52)
	end
	
	subroutine geta8(head,a8)
	real*8 a8(8)
	character*80 head(72)
	ipos=indexpos(head,'A81     ')-1
	do 10 i=1,8
10	read(head(i+ipos)(12:),*)a8(i)
	end

        function indexpos(head,f1)
        character*80 head(72),f1*8
        do 10 indexpos=1,72
10      if(head(indexpos)(1:8).eq.f1)return
        end

	subroutine xy_rad_xy(x,y,a8,b8,xx,yy)
        real x,y,xx,yy
        real*8 a8(8),alpha,xi,xn,tmp, b8(8),td,co
	real*8 cda,tda,sdb,cdb
	common /dec/cda,tda,sdb,cdb
        xi=a8(1)*x+a8(3)*y+a8(5)
        xn=a8(2)*x+a8(4)*y+a8(6)
	tmp=1d0-xn*tda
        alpha=datan(xi/cda/tmp)
c        delta=datan((xn+tda)*dcos(alpha)/tmp)
         td=(xn+tda)*dcos(alpha)/tmp
        alpha=alpha+a8(7)
c        td=dtan(delta)
        co=dcos(alpha-b8(7))
        tmp=sdb*td+cdb*co
        xi=dsin(alpha-b8(7))/tmp
        xn=(cdb*td-sdb*co)/tmp
        xx=b8(1)*xi+b8(3)*xn+b8(5)
        yy=b8(2)*xi+b8(4)*xn+b8(6)
        end

        subroutine xy2ad(x,a,y)
c 2 set 8 plate coef convert
        real*8 x(8),a(8),y(8),z
        z=x(1)*x(4)-x(3)*x(2)
        a(1)= x(4)/z
        a(2)=-x(2)/z
        a(3)=-x(3)/z
        a(4)= x(1)/z
        a(5)=(x(3)*x(6)-x(5)*x(4))/z
        a(6)=(x(5)*x(2)-x(1)*x(6))/z
	a(7)=2.*y(7)-x(7)
	a(8)=2.*y(8)-x(8)
        end

        function xmedian(a,n)
        real a(n)
        do 10 i=1,n-1
        do 10 j=i+1,n
        if(a(i).gt.a(j))goto 10
        x=a(i)
        a(i)=a(j)
        a(j)=x
10      continue
        m=n/2
        xmedian=a(m+1)
        if(m+m.ne.n)return
        xmedian=(a(m)+a(m+1))/2
        end


	subroutine star100(a,white,px,py,pz,n)
	integer*2 ix,iy,px(1),py(1)
	real a(2048,2048),pz(1)
	do 10 i=1,n
	z=0.
	do 20 ii=px(i)-7,px(i)+7
	do 20 jj=py(i)-7,py(i)+7
20	z=z+a(ii,jj)-white
10	pz(i)=z
	end

	subroutine getstar(a,starx,stary,starv,nstar,white)
	real      a(2048,2048)
	integer*2 b(2048,2048)
	real starx(1),stary(1),starv(1)
	do 5 j=1,2048
	do 5 i=1,2048
	k=a(i,j)+0.5
	if(k.lt.0)k=0
	if(k.gt.32767)k=32767
5	b(i,j)=k
	kt=2000
	kk=0
	ktw=kt+white
	do 10 i=30,2018
	do 10 j=30,2018
	if(b(i,j).gt.ktw)then
	  call center9(i,j,b,ktw,ix,iy,iv,n9)
	  if(iv.lt.ktw*2)goto 10
	  x=ix
	  y=iy
	  call star(a,2048,2048,x,y,4,ierr,white)
	  if(ierr.ne.0)x=-99.
	  z=(ix-x)**2+(iy-y)**2
	  if(z.gt.5.)goto 10
	  if(kk.eq.nstar)goto 11
	  kk=kk+1
	  starx(kk)=x
	  stary(kk)=y
	  starv(kk)=iv
	  if(mod(kk,10).eq.1)call dispdot()
	endif
10	continue
11	nstar=kk
	do 20 i=1,nstar-1
	do 20 j=i+1,nstar
	if(starv(i).ge.starv(j))goto 20
	x=starx(i)
	starx(i)=starx(j)
	starx(j)=x
	x=stary(i)
	stary(i)=stary(j)
	stary(j)=x
	x=starv(i)
	starv(i)=starv(j)
	starv(j)=x
20	continue
	if(nstar.gt.300)nstar=300
	end

        subroutine star(map,n1,n2,xx,yy,ir,ierr,sky)
c when ir=4, 15*15 for center
        real      map(n1,n2)
        integer*2 maps(-25:25,-25:25)
        real b(-25:25),xmap(-25:25,-25:25)

c sky first
        n25=25
        x=xx
        y=yy
        kx=nint(x)
        ky=nint(y)
        if(kx.le.n25)kx=n25+1
        if(kx.gt.n1-n25)kx=n1-n25
        if(ky.le.n25)ky=n25+1
        if(ky.gt.n2-n25)ky=n2-n25
        do 10 j=-n25,n25
        do 10 i=-n25,n25
10      maps(i,j)=map(kx+i,ky+j)
        n51=n25+n25+1
        ierr=1                                  ! error =1
        i5=ir*1.1+2
        do 1000 ll=1,7
         if(x.lt.1..or.x.gt.n1)return
         if(y.lt.1..or.y.gt.n2)return
        ix=nint(x)
        iy=nint(y)
cccccc x center    weight_center method
        do 30 i=ix-i5,ix+i5
        z=0.
        do 40 j=iy-i5,iy+i5
40      z=z+map(i,j)-sky          
30      b(i-ix)=z                       ! for weight_center cal.
        z=0.
        s=0.
        do 50 i=-i5,i5
        z=z+b(i)
50      s=s+b(i)*i
        s=s/z
        x=ix+s
cccccc y center
        do 130 j=iy-i5,iy+i5
        z=0.
        do 140 i=ix-i5,ix+i5
140     z=z+map(i,j)-sky
130     b(j-iy)=z
        z=0.
        s=0.
        do 150 i=-i5,i5
        z=z+b(i)
150     s=s+b(i)*i
        s=s/z
        y=iy+s
1000    continue
        ix=nint(x)
        iy=nint(y)
        do 101 j=-n25,n25
        do 101 i=-n25,n25
101     xmap(i,j)=map(ix+i,iy+j)
        x=ix-25+center_1(xmap,sky,51,1)
        y=iy-25+center_1(xmap,sky,51,0)
         if(x.lt.1..or.x.gt.n1)return
         if(y.lt.1..or.y.gt.n2)return
        ierr=0
          xx=x
          yy=y
        return
        end
 
        function center_1(a,sky,n,ip)
        real a(51,51)
        integer b(51)
        skyn=sky*n
        if(ip.eq.1)then
        do 10 i=1,n
        x=0.
        do 20 j=1,n
20      x=x+a(i,j)
10      b(i)=x-skyn
        else
        do 110 i=1,n
        x=0.
        do 120 j=1,n
120     x=x+a(j,i)
110     b(i)=x-skyn
        endif
        call histat(b,mode,maxh,xmean,xpeak,fsigma,gsigma,n)
        center_1=xpeak
        end