c combine 2048*2048 pic.
 	parameter (nfile=24)	
	character*80 f0,f2,head(72)
	character*12 ac,dC
	integer c(8)
	real b8(8)
	data b8/8.,64.,6.,64.,8.,62.,6.,62./
	
	if(iargc().lt.1)then
	  write(*,*)
	  write(*,*)'    ******** calculate shift pixel *******' 
	  write(*,*)
	  write(*,*)'    frame batch_file'
	  write(*,*)
	  stop
	endif

	call getarg(1,f0)
	open(51,file=f0,status='old')
	i=0
20	i=i+1
	read(51,'(a)',end=10)f2
	goto 20
10	close(51)
	m=i-1
	if(m.gt.nfile)m=nfile

	open(51,file=f0,status='old')
	do 30 i=1,m
	read(51,'(a)')f2
	call gethead(f2,head)
c	write(*,*)f2(1:9)
	call dis_4(head,1,2048,1,2048,b8)
	read(head(4)(65:),*)ic1,ic2
	read(head(5)(65:),*)ir1,ir2
c	write(*,*)ic1,ic2,ir1,ir2
	call dis_4(head,ic1,ic2,ir1,ir2,b8)
30	continue
	close(51)
c--------------------
	write(*,*)'-------------------'
	do 40 i=1,8,2
	call toms(b8(i),ac,1)
	call toms(b8(i+1),dc,0)
40	write(*,*)ac,dc

	open(51,file=f0,status='old')
	do 50 i=1,m
	read(51,'(a)')f2
	call gethead(f2,head)
	read(head(4)(65:),*)ic1,ic2
	read(head(5)(65:),*)ir1,ir2
c	write(*,*)f2(1:9),ic1,ic2,ir1,ir2
	call dis_4a(head,b8,c)
	jc1=max(c(1),c(5))
        jc2=min(c(3),c(7))
        jr1=max(c(2),c(4))
        jr2=min(c(6),c(8))
	write(*,*)f2(1:9),jc1,jc2,jr1,jr2,'   min'
	jc1=min(c(1),c(5))
        jc2=max(c(3),c(7))
        jr1=min(c(2),c(4))
        jr2=max(c(6),c(8))
	write(*,*)'         ',jc1,jc2,jr1,jr2,'   max'
	j=25
        jc1=jc1-j
	jc2=jc2+j
	jr1=jr1-j
	jr2=jr2+j
	if(jc1.lt.ic1)jc1=ic1
	if(jc2.gt.ic2)jc2=ic2
	if(jr1.lt.ir1)jr1=ir1
	if(jr2.gt.ir2)jr2=ir2
c	write(*,*)'   max+25',jc1,jc2,jr1,jr2
50	continue
	close(51)

	end

        subroutine dis_4a(head,b8,c)
        character*80 head(72),ac*12,dc*12
        real*8 a8(8),aa(2,6),adcoef(2,6)
	real b8(8)
	integer c(8)
	pi=3.14159265
	cx=pi/12.
	cy=pi/180.
        ipos=indexpos(head,'A81     ')-1
        do 10 i=1,8
10      read(head(ipos+i)(11:),*)a8(i)
	aa(1,1)=a8(1)
	aa(2,1)=a8(2)
	aa(1,2)=a8(3)
	aa(2,2)=a8(4)
	aa(1,3)=a8(5)
	aa(2,3)=a8(6)
	aa(1,4)=a8(7)
	aa(2,4)=a8(8)
	call xytoad(aa,adcoef)
	alpha_c=a8(7)
	delta_c=a8(8)
	do 20 i=1,8,2
        x_a=b8(i)*cx
        x_d=b8(i+1)*cy
        call standc(alpha_c,delta_c,x_a,x_d,xi,xn)
        x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
        y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
c20	write(*,*)nint(x),nint(y)
	c(i)=nint(x)
	c(i+1)=nint(y)
        call xy_rad(x,y,a8,alpha,delta)
        call toms(alpha,ac,1)
        call toms(delta,dc,0)
c20	write(*,*)ac,dc
20	continue
	end

        subroutine dis_4(head,ic1,ic2,ir1,ir2,b8)
        character*80 head(72),ac*12,dc*12
        real*8 a8(8)
	real b8(8)
	c1=ic1
	c2=ic2
	r1=ir1
	r2=ir2
        ipos=indexpos(head,'A81     ')-1
        do 10 i=1,8
10      read(head(ipos+i)(11:),*)a8(i)

        call xy_rad(c1,r1,a8,alpha,delta)
        call toms(alpha,ac,1)
        call toms(delta,dc,0)
	if(alpha.lt.b8(1))b8(1)=alpha
	if(delta.lt.b8(2))b8(2)=delta
c	write(*,*)'L_u:',ac,dc

        call xy_rad(c2,r1,a8,alpha,delta)
        call toms(alpha,ac,1)
        call toms(delta,dc,0)
	if(alpha.gt.b8(3))b8(3)=alpha
	if(delta.lt.b8(4))b8(4)=delta
c	write(*,*)'R_u:',ac,dc

        call xy_rad(c1,r2,a8,alpha,delta)
        call toms(alpha,ac,1)
        call toms(delta,dc,0)
	if(alpha.lt.b8(5))b8(5)=alpha
	if(delta.gt.b8(6))b8(6)=delta
c	write(*,*)'L_d:',ac,dc

        call xy_rad(c2,r2,a8,alpha,delta)
        call toms(alpha,ac,1)
        call toms(delta,dc,0)
	if(alpha.gt.b8(7))b8(7)=alpha
	if(delta.gt.b8(8))b8(8)=delta
c	write(*,*)'R_d:',ac,dc
	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(x,y,a8,alpha,delta)
c 2048*2048 ccd's X Y (in pixel) -------------> alpha delta (in h,d)
        real x,y,alpha,delta
        real*8 a8(8)                            ! from fits_head "A81"--"A88"
        pi=4.*atan(1.)
        xi=a8(1)*x+a8(3)*y+a8(5)
        xn=a8(2)*x+a8(4)*y+a8(6)
        rc=a8(7)                                ! center of ccd (in rad.)
        dc=a8(8)
        cd=cos(dc)
        td=tan(dc)
        alpha=atan(xi/cd/(1.-xn*td))
        delta=atan((xn+td)*cos(alpha)/(1.-xn*td))
        alpha=alpha+rc
        if(abs(dc-delta).gt.2.)then
          delta=-delta
          alpha=alpha+pi
        endif
        alpha=alpha*12./pi
        delta=delta*180./pi
        end

