c  add auto select max. error of star in determine 6 coef.
c         use file "gsc.tmp" produced by "gsc" to determine 6 plate coef.
c         use left_button pick > 3 stars, then right_button
c         calaulate 6 plate ceof
c same principle of gsc, current unit in hour and degree
c do not conside epoch, epoch is taken from gsc.tmp result
c modify 94,4,14
c if need to add seconds to TIME, file "time.par" format examples:
c +43			! comment here
c 1994,5,27 for sun
c 1994,10,4 add *.fit filename
c 94,11,17 fitshead
c 94,11,27 auto match stars.
c 94,11,28 all step is freehand
c 96,1,19, add seeing
	character*30 f1,cc*1
	integer*2 b(512*512)
 	integer*2 a(512*512*2)		! -32 change to 16 when cal.
	character*80 head(72)
	real*8 a8(8)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        pi=4.*atan(1.0)
	cx=pi/12.
	cy=pi/180.
  	phi=(40.+23/60.+36/3600.)*pi/180.	! XL latitude
  	sphi=sin(phi)
	cphi=cos(phi)
	k=iargc()
        if(k.eq.0)stop 'usage:  coord filename [!] 2003.11'  ! p means print
c check file
        call getarg(1,f1)
        if(inqh(f1,n1,n2,head,ibit).eq.0)stop 'DATA not found'
	if(n1.ne.512.or.n2.ne.512)stop 'not 512*512 CCD'
	call getpar(head)				! get s0,iy,im,id,ut
	idone=0
	if(indexpos(head,'A81     ').ne.73)then
	  idone=1
	  ipos=indexpos(head,'A81     ')-1
          do 11 i=1,8
11        read(head(ipos+i)(11:),*)a8(i)
	endif
	iprint=0
	if(k.ge.2)then
	  call getarg(2,cc)
	  if(cc.eq.'p')iprint=1
	  if(cc.eq.'!')idone=0
	endif
	if(idone.eq.1.and.iprint.eq.0)stop 'coord already done!'
	if(kgetade(head).ne.0)stop 'bad coordination of fits_head'
c read data
10	read(head(3)(65:),'(i5)')ii
	ii=ii/2+1
	call readata(a,b,f1,ibit,ii)
	iturn=0
	if((iy-1900)*10+im.ge.945.or.iy.lt.1990)iturn=1
c turn image
        m1=512
        m2=512

	ipos=indexpos(head,'INSTRUME')
	if(ipos.ne.73)then
          do i=12,20
	  if(head(ipos)(i:i).eq.'3')iturn=0
	  if(head(ipos)(i:i).eq.'4')iturn=4
	  enddo
	endif
        if(iturn.eq.0)call turnx2(b,m1,m2,-180)
        if(iturn.eq.1)call turnx2(b,m1,m2,-360)
        if(iturn.eq.4)call turnx2(b,m1,m2,180)

c main_job
	call main_job(head,a,a(512*512+1),a8,b,f1,iturn,
     *                match,rms,idone,iprint,seeing)
c write something to fits head
	call put_a8(head,a8,match,rms,seeing)	! 45-52	 write down a8
	call put_gr(head)		! 53-57  write gal_co. redden
	call put_xmass(head)		! 15-17, 58 write center,h_angle,airmass
	call put_aa(head)		! 59,60  write Azimuth, Altitude
	call put_moon(head,a8)		! 61-65  phase, position, pos_angle
	close (1)
	if(iprint.eq.1)then
          call writetail(head)
          call system('print pgplot.hp')
	endif
	if(idone.eq.1)stop
	write(*,*) ' re_write fits head !'
	i=lnblnk(f1)+1
	f1(i:i)=char(0)
	call tvcoordc(f1,head)
	end
	
	function kgetade(head)
	character*80 head(72)
	character*1 c
	character*20 line
	common /platec/alpha_c,delta_c,epoch
	kgetade=1
	ipos=indexpos(head,'RA      ')
	if(ipos.eq.73)return
	line=head(ipos)(12:31)
	line(20:20)=char(0)
	do 15 i=1,19
	if(line(i:i).eq.':')line(i:i)=','
15	if(ichar(line(i:i)).eq.39)line(i:i)=' '
	read(line(1:),*,err=10)i,j,x
	alpha_c=i+j/60.+x/3600.
c	write(*,*)'alpha',alpha_c
	
c	i=19
c	if(head(ipos)(28:28).eq.'.')i=20
c	if(head(ipos)(29:29).eq.'.')i=21
c	read(head(ipos)(i:),'(i2,1x,i2,1x,f4.1)',err=10)i,j,x
c	alpha_c=i+j/60.+x/3600.
c	write(*,*)'alpha',alpha_c

	ipos=indexpos(head,'DEC     ')
	if(ipos.eq.73)return
	line=head(ipos)(12:31)
	line(20:20)=char(0)
	c=' '
	do 20 i=1,19
	if(line(i:i).eq.'-')c='-'
	if(line(i:i).eq.'-')line(i:i)=' '
	if(line(i:i).eq.':')line(i:i)=','
20	if(ichar(line(i:i)).eq.39)line(i:i)=' '
	read(line(1:),*,err=10)i,j,x
	delta_c=i+j/60.+x/3600.
	if(c.eq.'-')delta_c=-delta_c
c	write(*,*)'delta',delta_c

c	i=19
c	if(head(ipos)(29:29).eq.'.')i=20
c	if(head(ipos)(i+1:i+1).eq.'-')then
c	read(head(ipos)(i+1:),'(a,i1,1x,i2,1x,f4.1)',err=10)c,i,j,x
c	else
c	read(head(ipos)(i:),'(a,i2,1x,i2,1x,f4.1)',err=10)c,i,j,x
c	endif
c	delta_c=i+j/60.+x/3600.
c	if(c.eq.'-')delta_c=-delta_c
c	write(*,*)'delta',delta_c

	ipos=indexpos(head,'EPOCH   ')
	if(ipos.eq.73)write(*,*)' no keywords EPOCH  ='
	if(ipos.eq.73)then
	  epoch=1995.1
	else
	i=24
	  if(head(ipos)(29:29).eq.'.')i=25
	  read(head(ipos)(i:),'(f6.1)',err=10)epoch
	  if(epoch.lt.1950.or.epoch.gt.2050.)return
	endif
	kgetade=0
	return
10	write(*,*)' wrong R.A. or DEC data'
	end

	subroutine putdata2(a,b)
	integer*2 a(1),b(1)
	do 10 i=1,512*512
10	b(i)=a(i)
	end

        subroutine center9i(i1,j1,a,kt,ix,iy,iv,n1)
        integer t1(200)
        integer*2 c1(200),a(512,512)
        iv=0
        n1=0
        i=i1
        j=j1
10      k=0
        iu=j
        id=j
        if(a(i,j).lt.kt)goto 100        ! find end
        k=k+a(i,j)
        a(i,j)=-a(i,j)

20      if(a(i,iu-1).lt.kt)goto 30      ! up find end
        iu=iu-1
          if(iu.lt.1)goto 30
        k=k+a(i,iu)
        a(i,iu)=-a(i,iu)
        goto 20
30      if(a(i,id+1).lt.kt)goto 40      ! down find end
        id=id+1
          if(id.gt.512)goto 40
        k=k+a(i,id)
        a(i,id)=-a(i,id)
        goto 30
40      iv=iv+k
        if(n1.lt.200)n1=n1+1
        t1(n1)=k
        c1(n1)=iu+id
        i=i+1
        j=c1(n1)*.5
        if(a(i,j).lt.a(i,j+1))j=j+1
        goto 10

100     iu=t1(1)
        id=1
        do 110 k=2,n1
        if(t1(k).gt.iu)then
          iu=t1(k)
          id=k
        endif
110     continue
        ix=i1+id-1
        iu=0
        do 120 k=1,n1
120     iu=iu+c1(k)
        iy=iu/2/n1
        end


	subroutine getstar(a,b,iturn,starx,stary,starv,nstar,white,seeing)
	integer*2 a(512,512),b(512*512)
	real starx(1),stary(1),starv(1)
	real see400(400)
c	call pgbegin(1,'/su',1,1)
c	call pgenv(512.,1.,1.,512.,1,1)
	kt=3000
5	kk=0
	ktw=kt+white
	call putdata2(b,a)
	do 10 i=10,502
        if(i.lt.50.and.kk.gt.10.and.kt.eq.3000)goto 30
	do 10 j=10,502
	if(a(i,j).gt.ktw)then
	  call center9i(i,j,a,ktw,ix,iy,iv,n9)
	  if(iv.lt.ktw*2)goto 10
	  x=ix
	  y=iy
	  call star(b,512,512,x,y,4,ierr,0,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
	  see400(kk)=see(b,512,512,x,y,white)
	  if(iturn.eq.1)then
	    starx(kk)=513.-x
	    stary(kk)=513.-y
	  endif
	  if(iturn.eq.4)then
	    starx(kk)=513.-x
	    stary(kk)=513.-y
	  endif
	  starv(kk)=iv
	  if(mod(kk,10).eq.1)call dispdot()
c	  write(*,*)kk,ix,iy,iv,n9,x,y,z
c	  n9=alog(iv/float(kt))+20
c	  if(n9.lt.20)n9=20
c	  if(n9.gt.27)n9=27
c	  call pgpoint(1,float(ix),float(iy),n9)
c	  call pgpoint(1,x,y,5)
	endif
10	continue
11	write(*,*)
	if(kk.lt.20.and.kt.eq.3000)then
	  kt=1000
	  if(kk.lt.10)kt=300
	  goto 5
	endif
        if(kk.lt.10)stop 'Fail:  too less stars'
	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
	seeing=xmedian(see400,nstar)
	if(nstar.gt.300)nstar=300
        return
30      kt=8000
        write(*,*)
        goto 5
	end

	function see(map,n1,n2,x,y,white)
c (a81+a84)/2 *180*3600/pi = arcsec / pixel
c schmidt is 1.709"/pixel
c seeing= 2.355*sigma*1.709
        integer*2 map(n1,n2)
        integer b(113)        ! (15-1)*8+1
        ix=x+0.5
        iy=y+0.5
        j=1
        do 10 i=ix-7,ix+7
        b(j)=map(i,iy)-white
10      j=j+8
        do 20 j=1,112,8
20      b(j+4)=(b(j)+b(j+8))*0.5+0.5
        do 30 j=1,112,4
30      b(j+2)=(b(j)+b(j+4))*0.5+0.5
        do 40 j=1,112,2
40      b(j+1)=(b(j)+b(j+2))*0.5+0.5
        call histat(b,mode,maxh,xmean,xpeak,fsigma1,gsigma,113)
        j=1
        do 110 i=iy-7,iy+7
        b(j)=map(ix,i)-white
110     j=j+8
        do 120 j=1,112,8
120     b(j+4)=(b(j)+b(j+8))*0.5+0.5
        do 130 j=1,112,4
130     b(j+2)=(b(j)+b(j+4))*0.5+0.5
        do 140 j=1,112,2
140     b(j+1)=(b(j)+b(j+2))*0.5+0.5
        call histat(b,mode,maxh,xmean,xpeak,fsigma2,gsigma,113)
c       write(*,'(a,f9.2,a)')'seeing:',(fsigma1+fsigma2)/16.*2.355*1.709,'"'
        see=(fsigma1+fsigma2)/16.*2.355*1.709
        end

	subroutine getmxy(a8,x,y,z,phase,pl,pb,xx)
	real*8 a8(8),b8(12),jd
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	call xytoad(a8,b8)
	j=(iy+((im-1)*30.4+id)/365.)*10.+0.5
	epoch_g=j*0.1
	x=a8(7)/cx
	y=a8(8)/cy
	call astprs(x,y,2000.,alpha_g,delta_g,epoch_g)
	call ad_azl(alpha_g,delta_g,a1,d1)
	alpha_g=alpha_g*cx
	delta_g=delta_g*cy
	call azl_ad(a1,d1-0.1,x,y)
	call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	a3=b8(1)*xi+b8(3)*xn+b8(5)
	d3=b8(2)*xi+b8(4)*xn+b8(6)
	call azl_ad(a1,d1+0.1,x,y)
	call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	a2=b8(1)*xi+b8(3)*xn+b8(5)
	d2=b8(2)*xi+b8(4)*xn+b8(6)
	z=atan2(a2-a3,d2-d3)/cy
	call fjd(jd,iy,im,id)
	jd=jd+ut/24.
	call moon2(jd,a2,d2)				! get monn l,b
	call lb_ad(a2,d2,x,y)				! get moon alpha delta
	phase=amod(a2-sunlamda(jd)+360.,360.)/360.*29.53  ! get moon phase
	call ad_azl(x,y,a2,d2)				! get moon A, H
	pl=a2
	pb=d2
	x=sin(d1*cy)*sin(d2*cy)+cos(d1*cy)*cos(d2*cy)*cos(a1*cy-a2*cy)
        y=acos(x)                                    
	xx=y/cy
	cosA=(sin(d2*cy)-sin(d1*cy)*x)/cos(d1*cy)/sin(y)
	if(cosa.gt.0.999999)cosa=0.999999
	if(cosa.lt.-.999999)cosa=-.999999
	y=acos(cosA)/cy
	x=a2-a1
	if(abs(x).gt.180.)x=a1-a2
	if(x.lt.0.)y=-y
	z=270.+z-y
	z=amod(z,360.)
	if(z.lt.0.)z=z+360.
	x=cos(z*cy)
	y=sin(z*cy)
	end

	subroutine azl_ad(a,h,alpha,delta)	!dd->hd
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	aa=a*cy
	hh=h*cy
	sina=sin(aa)
	cosa=cos(aa)
	sinh=sin(hh)
	cosh=cos(hh)
	sinde=sphi*sinh-cphi*cosa*cosh
	cosde=sqrt(1.-sinde*sinde)
	sint=cosh*sina/cosde
	cost=(cosa*cosh+cphi*sinde)/(sphi*cosde)
        delta=atan(sinde/cosde)/cy
        t=atan(sint/cost)
        if(cost.lt.0.)t=t+pi
        alpha=s0-t/cx	
	if(alpha.lt.0.)alpha=alpha+24.
	if(alpha.gt.24.)alpha=alpha-24.
	end

	subroutine ad_azl(alpha,delta,a,h)     ! hd,dd
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        hour=(s0-alpha)*cx
        delt=delta*cy
        s_delta=sin(delt)
	c_delta=cos(delt)
        s_hour=sin(hour)
	c_hour=cos(hour)
        sinh=sphi*s_delta+cphi*c_delta*c_hour
        cosh=sqrt(1.-sinh*sinh)
        sinA=c_delta*s_hour/cosh
        cosA=(sphi*c_delta*c_hour-cphi*s_delta)/cosh
        h=atan(sinh/cosh)/cy
        a=atan(sinA/cosA)
        if(cosA.lt.0.)a=a+pi
	if(a.lt.0.)a=a+pi+pi
	if(a.gt.pi+pi)a=a-pi-pi
	a=a/cy
	end

	subroutine put_aa(head)     ! 59,60
	character*80 head(72)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        ipos=indexpos(head,'HA      ')
	read(head(ipos)(19:),'(i2,1x,i2,1x,f5.2)')i,j,x
	t=i+j/60.+x/3600.
	ipos=indexpos(head,'DEC     ')
	read(head(ipos)(20:),'(i2,1x,i2,1x,f4.1)')i,j,x
	d=i+j/60.+x/3600.
	if(head(ipos)(19:19).eq.'-')d=-d
	call ad_azl(s0-t,d,a,h) 		    ! hd,dd
	ipos=indexpos(head,'A81     ')+14           ! 59 row
        head(ipos)(1:32)='AZIMUTH =                      /'
        head(ipos+1)(1:32)='ALTITUDE=                      /'
 	head(ipos)(34:76)='Azimuth is measured from south through west'
 	head(ipos+1)(34:54)='both units in degrees'
	write(head(ipos)(21:30),'(f10.2)')a        ! in degree
	write(head(ipos+1)(21:30),'(f10.2)')h
	end

	subroutine put_moon(head,a8)		! 61-65
	character*80 head(72)
	real*8 a8(8)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ')+16           ! 61 row
        head(ipos)(1:32)='MPHASE  =                      /'
        head(ipos+1)(1:32)='MAZIMUTH=                      /'
        head(ipos+2)(1:32)='MALITIUD=                      /'
        head(ipos+3)(1:32)='MANGLE  =                      /'
        head(ipos+4)(1:32)='MDIRECT =                      /'
	head(ipos)(34:71)='Define the length of Moon_Month 29.53d'
 	head(ipos+1)(34:78)='Moon Azimuth is measured from south thr. west'
 	head(ipos+2)(34:54)='both units in degrees'
	head(ipos+3)(34:75)='Position angle of Moon to CCD_Field center'
	head(ipos+4)(34:77)='Moon direction in the X_Y plane of CCD frame'
	call getmxy(a8,x,y,z,phase,pl,pb,xx)
	write(head(ipos)(21:30),'(f10.1)')phase
	write(head(ipos+1)(21:30),'(f10.2)')pl
	write(head(ipos+2)(21:30),'(f10.2)')pb
	write(head(ipos+3)(21:30),'(f10.2)')xx		
        write(head(ipos+4)(21:30),'(f10.2)')z	
	end
	
	function sunlamda(jd)			! in degree
c sun lamda in degree, beta==0
	real*8 jd,n,l,g,x
	n=jd-2451545.d0
	l=280.460d0+0.9856474d0*n
	g=357.528d0+0.9856993d0*n
	x=dmod(g,360d0)*datan(1d0)/45d0
	x=l+1.915d0*dsin(x)+0.020d0*dsin(x+x)
	x=dmod(x,360d0)
	if(x.lt.0d0)x=x+360d0
	sunlamda=x
	end

	subroutine moon2(jd,pl,pb)		!in degree
	real*8 jd,t,g,r,x1,x2,x3,x4,x5,x6
	g=datan(1d0)/45d0
	t=(jd-2451545.d0)/36525.d0
	r=218.32d0+dmod(481267.833d0*t,360d0)
	x1=(134.9d0+dmod(477198.85d0*t,360d0))*g
	x2=(259.2d0-dmod(413335.38d0*t,360d0))*g
	x3=(235.7d0+dmod(890534.23d0*t,360d0))*g
	x4=(269.9d0+dmod(954397.70d0*t,360d0))*g
	x5=(357.5d0+dmod( 35999.05d0*t,360d0))*g
	x6=(186.6d0+dmod(966404.05d0*t,360d0))*g
	r=r+6.29d0*dsin(x1)-1.27d0*dsin(x2)+0.66d0*dsin(x3)
     *     +0.21d0*dsin(x4)-0.19d0*dsin(x5)-0.11d0*dsin(x6)
	r=dmod(r,360d0)
	if(r.lt.0d0)r=r+360d0
	pl=r
	x1=( 93.3d0+dmod(483202.03d0*t,360d0))*g
	x2=(228.2d0+dmod(960400.87d0*t,360d0))*g
	x3=(318.3d0+dmod(  6003.18d0*t,360d0))*g
	x4=(217.6d0-dmod(407332.20d0*t,360d0))*g
	r=5.13d0*dsin(x1)+0.28d0*dsin(x2)-0.28d0*sin(x3)-0.17d0*dsin(x4)
	r=dmod(r,360d0)
	if(r.gt.90d0)r=r-360d0
	pb=r
	end

	subroutine lb_ad(pl,pb,alpha,delta)       !dd-->hd
	cy=atan(1.)/45.
	xl=cos(pb*cy)*cos(pl*cy)
	xm=0.9175*cos(pb*cy)*sin(pl*cy)-0.3978*sin(pb*cy)
	xn=0.3978*cos(pb*cy)*sin(pl*cy)+0.9175*sin(pb*cy)
	alpha=atan2(xm,xl)/cy/15.
	if(alpha.lt.0.)alpha=alpha+24.
	delta=asin(xn)/cy
	end

	subroutine ad_lb(al,de,pl,pb)		!hd-->dd
	cy=atan(1.)/45.
	cx=cy*15.
	xl=cos(de*cy)*cos(al*cx)
	xm=0.9175*cos(de*cy)*sin(al*cx)+0.3978*sin(de*cy)
	xn=-.3978*cos(de*cy)*sin(al*cx)+0.9175*sin(de*cy)
	pl=atan2(xm,xl)/cy
	if(pl.lt.0.)pl=pl+360.
	pb=asin(xn)/cy
	end

	subroutine put_a8(head,a8,match,rms,seeing)      !45-52
	character*80 head(72),ac*12,dc*12
	real*8 a8(8)
        common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'VOLT3   ')
	if(ipos.ne.73)then
          write(head(ipos)(24:30),'(f7.2)')seeing
	  head(ipos)(34:49)='SEEING  (arcsec)'
	endif
	ipos=indexpos(head,'A81     ')
	if(ipos.eq.73)ipos=indexpos(head,'        ')
	if(ipos.gt.51)ipos=indexpos(head,'COMMEMT ')
	if(ipos.gt.51)ipos=51
	head(ipos)(1:46)='A81     =                      / 6 plate coef.'
	head(ipos+1)(1:46)='A82     =                      / matched star:'
	head(ipos+2)(1:46)='A83     =                      / RMS (arcsec):'
	head(ipos+3)(1:32)='A84     =                      /'
	head(ipos+4)(1:32)='A85     =                      /'
	head(ipos+5)(1:32)='A86     =                      /'
	head(ipos+6)(1:46)='A87     =                      / p_c (2000.0) '
	head(ipos+7)(1:46)='A88     =                      / in rad.      '
	do 11 i=1,6
11	write(head(ipos-1+i)(11:30),'(E20.12)')a8(i)
	write(head(ipos+6)(11:30),'(f20.7)')a8(7)
	write(head(ipos+7)(11:30),'(f20.7)')a8(8)
	write(head(ipos+1)(47:52),'(i6)')match
	write(head(ipos+2)(47:52),'(f6.2)')rms
        x=a8(7)/cx                              ! put lamda beta
        y=a8(8)/cy
        call ad_lb(x,y,alpha,delta)
        call toms(alpha,ac,2)
        call toms(delta,dc,0)
        head(ipos+3)(57:63)='l_sun ='
        head(ipos+4)(57:63)='b_sun ='
        head(ipos+3)(65:73)=ac(1:9)
        head(ipos+4)(65:73)=dc(1:9)
	end

	subroutine put_gr(head)          !51:53-57
	character*80 head(72)
	character*12 ac,dc
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ')+6
	read(head(ipos)(21:),*)xa			! center plate in rad.
	read(head(ipos+1)(21:),*)xd
	xa=xa/cx
	xd=xd/cy				! to h,d
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	head(ipos+2)(1:46)="RA2000  = '                  ' / plate center "
	head(ipos+3)(1:32)="DEC2000 = '                  ' /"
	head(ipos+2)(18:29)=ac
	head(ipos+3)(19:29)=dc(1:11)
        call astprs(xa,xd,2000.,x,y,1950.)
        call galactic(x*cx,y*cy,xa,xd)
        head(ipos+4)(1:32)='GALLONG =                      /'
        head(ipos+4)(34:62)='galactic coordinates (1950.0)'
        head(ipos+5)(1:32)='GALLATI =                      /'
        head(ipos+5)(34:54)='both units in degrees'
	write(head(ipos+4)(21:30),'(f10.2)')xa        ! in degree
	write(head(ipos+5)(21:30),'(f10.2)')xd
	head(ipos+6)(1:32)='EXTINCT =                      /'
	head(ipos+6)(34:63)='BH Bmag, EXTINCTION = 4*E(B-V)'
	xa=redden(xa,xd)
	write(head(ipos+6)(21:30),'(f10.3)')xa
	end

	subroutine put_xmass(head)          !15,16,17,58
	character*80 head(72)
	character*12 ac,dc
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ')+6
	read(head(ipos)(21:),*)xa			! center plate in rad.
	read(head(ipos+1)(21:),*)xd
	xa=xa/cx
	xd=xd/cy					! to h,d
        j=(iy+((im-1)*30.4+id)/365.)*10.+0.5
	epoch=j*0.1
        call astprs(xa,xd,2000.,x,y,epoch)
	call toms(x,ac,1)
	call toms(y,dc,0)
	ipos=indexpos(head,'RA      ')
	head(ipos)(11:30)="'                  '"
	head(ipos)(18:29)=ac
	write(head(ipos)(50:),"('(',f6.1,')')")epoch
	ipos=indexpos(head,'EPOCH   ')
	if(ipos.ne.73)head(ipos)(11:30)="'                  '"
	if(ipos.ne.73)write(head(ipos)(24:29),'(f6.1)')epoch
	ipos=indexpos(head,'DEC     ')
	head(ipos)(11:30)="'                  '"
	head(ipos)(19:29)=dc(1:11)
	write(head(ipos)(50:),"('(',f6.1,')')")epoch
	h=amod(s0-x,24.)				! get HA, here I wast
        if(h.lt.0.)h=h+24.				! more time, x->xa
	call toms(h,ac,1)
	ipos=indexpos(head,'HA      ')
	head(ipos)(11:30)="'                  '"
	head(ipos)(18:29)=ac
	xmm=xmass(h*cx,xd*cy)				! middle time
	ipos=indexpos(head,'EXPOSURE')
	read(head(ipos)(20:),*)k				! exp time
	xmb=xmass((h-k/7200.)*cx,xd*cy)			! begin time
	xme=xmass((h+k/7200.)*cx,xd*cy)			! end time
	x=(xmb+4.*xmm+xme)/6.
	ipos=indexpos(head,'A81     ')+13
	head(ipos)(1:32)='AIRMASS =                      /'
	write(head(ipos)(21:30),'(f10.3)')x
	end

	subroutine getpar(head)
c get s0 iy im id ut
	character*80 head(72),c9*9
	real*8 d1,d
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'DATE-OBS')
        read(head(ipos)(22:),'(i2,1x,i2,1x,i2)')id,im,iy        ! read obs date
	iy=iy+1900
	if(iy.lt.1950)iy=iy+100
	ipos=indexpos(head,'TIME    ')
	read(head(ipos)(20:),'(i2,1x,i2,1x,f4.1)')i,j,x
	ipos=indexpos(head,'EXPOSURE')
	read(head(ipos)(20:),*)k				! exp time
	ut=i+j/60.+x/3600.+k/7200.			! center UT time, in h
        call fjd(d1,1992,12,31)			! 1992,12,31
        d2=(6+38/60.+40.1954/3600.)/24.         ! 6 38 40.1954
        x1=7+50/60.+18.344/3600.                ! 7:50:18.344 lamda of XL
        call fjd(d,iy,im,id)
        d=(d-d1)*1.0027379093d0+d2
        d=dmod(d,1d0)
        if(d.lt.0d0)d=d+1d0
        s0=d*24d0
        s0=s0+x1+ut*1.002738
	end

        subroutine fjd(d,iy,im,id)
        real*8 d
        i=(im+9)*0.09
        i=(i+iy-1900)*1.75+0.01
        j=30.56*im
        j=(iy-1950)*367+id+j-i
        d=2433338.5d0+j
        end

	function xmass(angle,delta)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
  	x=sin(delta)*sphi+cos(delta)*cos(angle)*cphi
c  	z=90.-atan2(x,sqrt(1.-x*x))*180./pi      ! Zenth
  	x=1./x-1.
  	xmass=1.+x-x*(0.0018167+x*(0.002875+x*0.0008083))
	end

	subroutine put_a2(head)
	character*80 head(72)
	character*35 dis
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ')
	if(ipos.eq.73)ipos=indexpos(head,'        ')
	if(ipos.gt.51)ipos=indexpos(head,'COMMEMT ')
	if(ipos.gt.51)ipos=51
	head(ipos)(1:32)='A81     =                      /'
        head(ipos)(34:53)='6 plate coefficients'
	head(ipos+1)(1:32)='A82     =                      /'
	head(ipos+2)(1:32)='A83     =                      /'
	head(ipos+3)(1:32)='A84     =                      /'
	head(ipos+4)(1:32)='A85     =                      /'
	head(ipos+5)(1:32)='A86     =                      /'
	head(ipos+6)(1:46)='A87     =                      / p_c (2000.0) '
	head(ipos+7)(1:46)='A88     =                      / in radians   '
10	write(*,*)'input center R.A, Decl, [1950.0]: ex. 3:46:02.1 23:0'
	read(*,'(a)')dis
	if(itohd(dis,x,y,epoch).eq.0)goto 10
	if(epoch.eq.0.)epoch=1950.
	call astprs(x,y,epoch,x,y,2000.)
	write(head(ipos+6)(11:30),'(f20.7)')x*cx
	write(head(ipos+7)(11:30),'(f20.7)')y*cy
	end

        function redden(al,b)
c this subroutine written by burstein.
      INTEGER*4 IRED(1200),IHI(201)
      XRAD = 57.2958
      ALIN = AL/0.3 +  0.5
      ILIN = ALIN + 0.51
      BB = ABS(B)
      redden= -0.99
      if(bb.lt.10.)  return
      OPEN(UNIT=11,FILE='/EOD/bur/redsouth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=4800)
      OPEN(UNIT=12,FILE='/EOD/bur/rednorth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=4800)
      OPEN(UNIT=13,FILE='/EOD/bur/hinorth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=804)
      OPEN(UNIT=14,FILE='/EOD/bur/hisouth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=804)
      IF (B) 11,11,18
   11 IF (B + 62.) 13,12,12
   12 AREC = -(B+10.)/0.6 + 1.
      KZ = AREC + 0.51
      READ (11,rec=KZ) IRED
	call swap4(ired,4800)
      KQ = KZ -1
      IEBV = IRED(ILIN)
      GO TO 19
   13 ALC = AL/XRAD
      AREC = 101. + SIN(ALC)*(90.+B)/0.3
      ALIN = 101. + COS(ALC)*(90.+B)/0.3
      KW = AREC + 0.51
      ILIN = ALIN + 0.51
      READ (14,rec=KW) IHI
	call swap4(ihi,804)
      AHI = IHI(ILIN)
      BMV = -0.0372 + 0.357*AHI/10000.
      IEBV = BMV*1000. + 0.5
      KQ = KW - 1
      GO TO 19
   18 IF (B - 62.) 17,17,16
   17 AREC = (B-10.)/0.6 + 1.
      KY = AREC + 0.51
      READ (12,rec=KY) IRED
	call swap4(ired,4800)
      KQ = KY - 1
      IEBV = IRED(ILIN)
      GO TO 19
   16 ALC = AL/XRAD
      AREC = 101. + SIN(ALC)*(90.-B)/0.3
      ALIN = 101. + COS(ALC)*(90.-B)/0.3
      KX = AREC + 0.51
      ILIN = ALIN + 0.51
      READ (13,rec=KX) IHI
	call swap4(ihi,804)
      AHI = IHI(ILIN)
      BMV = -0.0372 + 0.357*AHI/10000.
      IEBV = BMV*1000. + 0.5
      KQ = KX - 1
   19 ABV = IEBV
      redden= ABV*4./1000. + 0.005
	close(11)
	close(12)
	close(13)
	close(14)
      return
      END

        subroutine galactic(alpha,delta,xl,xb)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        a0=282.25*cy
        theta=62.6*cy
        Xl0=33.0*cy
        x1=cos(delta)*cos(alpha-a0)
        x2=cos(delta)*sin(alpha-a0)*cos(theta)+sin(delta)*sin(theta)
        xl=atan2(x2,x1)+xl0
        if(xl.lt.0.)xl=xl+2.0*pi
        xb=asin(sin(delta)*cos(theta)-
     *          cos(delta)*sin(alpha-a0)*sin(theta))
        xl=xl/cy
        xb=xb/cy
        end

	subroutine readata(a,b,f1,ibit,ii)
	integer c(720),a(1),b(1)
	character*30 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=512
	n2=512
	n=1
55	k=index(ch,'END          ')
	if(k.eq.0)then
	  n=n+1
	  read(1,rec=n)ch
	  goto 55
	endif
c	write(*,"(' reading data ')")
	i=0
	if(ibit.lt.0)ibit=-ibit
	k1=n+1
	k2=n+(n1*n2*ibit/32-1)/nn
	do 60 k=k1,k2
c	if(mod(k,120).eq.0)call dispdot()
	read(1,rec=k)ch
	do 60 j=1,nn
	i=i+1
60	a(i)=c(j)
	close (1)
	if(ibit.eq.32)then
	  call swap4(a,n1*n2*4)
	  call shrink4(a,n1,n2,a,ii)
	else
	  call swap2(a,n1*n2*2)
	endif
 	call shrink2(a,n1,n2,b)		! 2048*2048-->512*512
c	write(*,*)
	end

	subroutine shrink2(a,n1,n2,b)
	integer*2 a(n1,n2)
	integer*2 b(512*512)
	k=0
	m=1
	do 10 j=1,n2,m
	do 10 i=1,n1,m
	k=k+1
10	b(k)=a(i,j)
	end

	subroutine shrink4(a,n1,n2,b,ii)
	real a(n1*n2)
	integer*2 b(n1*n2)
	do 10 i=1,n1*n2
	x=a(i)/ii
        if(x.gt.32767.)x=32767.
        if(x.lt.0.)x=0.
10	b(i)=x
	end

	function inqh(f1,n1,n2,head,ibit)
	character*(*) f1
	character*80 head(72)
	logical inqh
	k=index(f1,'.')
	if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
	inquire(file=f1,exist=inqh)
	if(.not.inqh)return
	open(1,file=f1,status='old',
     *       access='direct',recl=5760)
	read(1,rec=1)head
	read(head(2)(22:),*)ibit
	read(head(4)(22:),*)n1	!NAXIS1
	read(head(5)(22:),*)n2	!NAXIS2
c	ipos=indexpos(head,'OBJECT  ')
c	write(*,'(1x,a)')head(ipos)(1:70)       	! object
	ipos=indexpos(head,'DATE-OBS')
	head(ipos)(63:72)='(dd/mm/yy)'
c	do 10 i=ipos,ipos+5				! da.ti.ex.ra.de.ha
c10	write(*,'(1x,a)')head(i)(1:78)
c	ipos=indexpos(head,'EPOCH   ')
c        if(ipos.ne.73)write(*,'(1x,a)')head(ipos)(1:78)
        ipos=indexpos(head,'TELESCOP')
        if(ipos.ne.73)head(ipos)(12:28)='BAO 60-90 SCHIMDT'
        ipos=indexpos(head,'BZERO   ')
        if(ipos.ne.73)head(ipos)(1:1)='P'
        ipos=indexpos(head,'BSCALE  ')
        if(ipos.ne.73)head(ipos)(1:1)='P'
	ipos=indexpos(head,'TIME    ')
	if(ipos.ne.73)head(ipos)(42:57)='OF START OF OBS.'
	ipos=indexpos(head,'HA      ')
	if(ipos.ne.73)head(ipos)(45:61)='OF MIDDLE OF OBS.'
	close (1)
	end

	subroutine main_job(head,map,map1,a8,ma,f1,iturn,
     *                      match,rms,idone,iprint,seeing)
	character*80 head(72)
	integer*2 map(512,512),ma(512,512)	! map buffer
	integer*2 map1(512,512)		! 2nd buffer
	real*8 a8(8)
	character*12 ac,dc,c3*3			! alpha delta display line
	character*50 ns,dis*60			! for display
	real x9(99),y9(99)			! al_del of starard posi star
	real sx(4),sy(4)			! scale on right_bottom of pic
	real*8 adcoef(2,6),xycoef(2,6),dummy8		! 2 set 6 plte coef
	real starx(400),stary(400),starv(400)   ! take ccd stars
	parameter (ngsc=3999)		        ! No. of gsc stars.
	real aagsc(ngsc),ddgsc(ngsc),mmgsc(ngsc)
	character f1*30
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	common /platec/alpha_c,delta_c,epoch
	call whitexblack(ma,512,512,white,x)
	black=white+20.*x
	n1=512
	n2=512
	ir=4
	if(idone.eq.1)goto 1001
******** auto match
c	write(*,*)'wait, auto match stars. Now, pick out CCD stars'
	nstar=400
	call getstar(map1,map,iturn,starx,stary,starv,nstar,white,seeing)
	ipp=0
        jpp=1
	if(nstar.ge.300)jpp=20
1012    if(jpp.ne.1)call ijpp(starx,stary,starv,nstar,jpp)
c******** take out gsc stars
	epoch1=epoch
	alpha1=alpha_c
	delta1=delta_c
1011	call fgsc(alpha_c,delta_c,epoch,15.,aagsc,ddgsc,mmgsc,ng)
        if(jpp.ne.1)call ijpp(aagsc,ddgsc,mmgsc,ng,jpp)
	call astprs(alpha_c,delta_c,epoch,xa,xd,2000.)
	epoch=2000.
        xycoef(1,1)=-0.83e-5
	xycoef(2,1)=-0.10e-6
	xycoef(1,2)= 0.10e-6
	xycoef(2,2)=-0.83e-5
	xycoef(1,3)= 0.21e-2
	xycoef(2,3)= 0.21e-2
         if(iturn.eq.4)then
           xycoef(1,1)=-0.830e-5
   	   xycoef(2,1)=-0.225e-6
	   xycoef(1,2)=-0.217e-6
	   xycoef(2,2)= 0.829e-5
	   xycoef(1,3)= 0.872e-2
	   xycoef(2,3)=-0.826e-2
         endif
        
	alpha_c=xa*cx
	delta_c=xd*cy
	ng=30
	rms=100.
	match=kgsc(aagsc,ddgsc,mmgsc,ng,
     *	       starx,stary,starv,ng,		
     *		xycoef,alpha_c,delta_c,rms)
	if(match.eq.0)then
	  call shift_c(starx,stary,aagsc,ddgsc,ng,xycoef,alpha_c,delta_c)
	  match=kgsc(aagsc,ddgsc,mmgsc,ng,
     *	       starx,stary,starv,ng,	
     *		xycoef,alpha_c,delta_c,rms)
	  if(match.eq.0)then
	    if(ipp.eq.2)stop 'match fail_1(position_error)!'
	    ipp=ipp+1
	    i=epoch1+0.1
	    if(i.eq.1950)epoch=2000.0
	    if(i.ge.1995)epoch=1950.0
	    alpha_c=alpha1
	    delta_c=delta1
	    if(ipp.eq.2)then
	      epoch=epoch1
	      if(jpp.ne.1)stop 'match fail_2(position_error)!'
	      if(nstar.gt.60)jpp=20
	      goto 1012
	    endif
	    write(*,*)'I assume epoch =',epoch
	    goto 1011
	  endif
	endif
	call fgsc(alpha_c/cx,delta_c/cy,epoch,20.,aagsc,ddgsc,mmgsc,ng)
        if(jpp.ne.1)call ijpp(aagsc,ddgsc,mmgsc,ng,jpp)
	rms=5.
	match=kgsc(aagsc,ddgsc,mmgsc,ng,
     *	       starx,stary,starv,nstar,
     *		xycoef,alpha_c,delta_c,rms)
	if(match.eq.0)stop 'match fail (star too faint) !'
********* save a8
	a8(1)=xycoef(1,1)
	a8(2)=xycoef(2,1)
	a8(3)=xycoef(1,2)
	a8(4)=xycoef(2,2)
	a8(5)=xycoef(1,3)
	a8(6)=xycoef(2,3)
******95,12,19 add
        if(alpha_c.lt.0.)alpha_c=alpha_c+pi+pi
	a8(7)=alpha_c
	a8(8)=delta_c			! 2000 epoch
********* save a8 end
1001	if(iprint.eq.0)return
	xycoef(1,1)=a8(1)
	xycoef(2,1)=a8(2)
	xycoef(1,2)=a8(3)
	xycoef(2,2)=a8(4)
	xycoef(1,3)=a8(5)
	xycoef(2,3)=a8(6)
	alpha_c=a8(7)
	delta_c=a8(8)
	epoch=2000.
	call toms(alpha_c/cx,ac,1)
	call toms(delta_c/cy,dc,0)
	write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
	call xytoad(xycoef,adcoef)	! 2 set of 6 coef convert
	xycoef(1,4)=alpha_c
	xycoef(2,4)=delta_c
        y=(abs(xycoef(1,1))+abs(xycoef(2,2)))*.5
        x=y*180.*60./pi*n1                      ! X size in minute
        i=nint(x/4.)
        k4=nint(i/x*n1)/4
        write(c3(1:),"(i2)")i
        c3(3:3)=char(39)
c hardcopy ****************
	call unlink('pgplot.hp')
        call pgbegin(0,'/hp',1,1)
        call pgslw(2)
        if(iturn.eq.4)then
          xycoef(1,2)=-xycoef(1,2)
          xycoef(2,2)=-xycoef(2,2)
          dummy8=xycoef(2,3)
          xycoef(2,3)=xycoef(1,3)
          xycoef(1,3)=-dummy8
          call xytoad(xycoef,adcoef)
        endif
 
c physical size of laser_printer, 300 dots/inch, 2336 dots in width
c                                 3200 dots in length
 	xp=17.06
        dy=300./xp
        x1=xp/2.54*300./2336.
        y1=xp/2.54*300./3200.
c create pg window with input scale
        x=(1.-x1)*0.5
        x1=(1.+x1)*0.5
        y=0.95
        y1=y-y1
        call pgadvance
        call pgvport(x,x1,y1,y)
        call pgwindow(1.,512.,1.,512.)
        call pgbox('bc',0,0,'bc',0,0)
  	call pgxgray(ma,512,512,1,512,1,512,white,black,2)
	call pgsch(0.8)
        call pglabel('R.A  South','Decl.  East',f1//ns)
	sx(1)=500
	sx(2)=sx(1)
	sx(3)=sx(1)-k4
	sx(4)=sx(3)
	sy(1)=15
	sy(2)=10
	sy(3)=sy(2)
	sy(4)=sy(1)
	call pgline(4,sx,sy)
	call pgtext(450.,15.,c3)
        y1=-4.*dy                       ! for print display line on paper
        x1=(1.-17./xp)*256.             ! in 2 coloumn, cal x position first
        call pgsch(.6)
c draw grid, if to near pole, not draw
        x=delta_c/cy
        x60=60.
        if(x.gt.60.)x60=30.
        if(x.gt.80.)x60=12.
        if(x.gt.86.)x60=6.
        if(x.gt.88.)x60=3.
        if(x.gt.89.)x60=2.
	call pgslw(1)
        call xy_rad(1.,1.,xycoef,a1,d1)
        call xy_rad(1.,512.,xycoef,a2,d2)
        call xy_rad(512.,1.,xycoef,a3,d3)
        call xy_rad(512.,512.,xycoef,a4,d4)
	x=amin1(a1,a2,a3,a4)
	y=amax1(a1,a2,a3,a4)
	i1=x*x60
	i2=y*x60+1
	x=amin1(d1,d2,d3,d4)
	y=amax1(d1,d2,d3,d4)
	j1=x*6.
	if(x.lt.0)j1=j1-1
	j2=y*6.
	if(y.gt.0)j2=j2+1
	do 153 i=i1,i2
	xa=i/x60*cx
	k=0
	do 152 j=j1,j2
	k=k+1
	xd=j/6.*cy
	call standc(alpha_c,delta_c,xa,xd,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)
	x9(k)=x*0.25
152	y9(k)=(513.-y)*0.25
	call pgline(k,x9,y9)
	call toms(i/x60,ac,1)
	ns=ac(2:3)//'\\uh\\d'//ac(5:6)//'\\um\\d'
153	if(x9(1).gt.0..and.x9(1).lt.512.)call pgtext(x9(1)-15.,-9.,ns)
	do 154 j=j1,j2
	xd=j/6.*cy
	k=0
	do 155 i=i1,i2
	k=k+1
	xa=i/x60*cx
	call standc(alpha_c,delta_c,xa,xd,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)
	x9(k)=x
155	y9(k)=(513.-y)
	call pgline(k,x9,y9)
	call toms(j/6.,dc,0)
	ns=dc(1:3)//'\\uo\\d'//dc(5:6)//char(39)
154	if(y9(k).gt.0..and.y9(k).lt.512.)call pgtext(-33.,y9(k),ns)
c draw AZIMUTH & ALTITUDE grid, if to near pole, not draw
159     j=(iy+((im-1)*30.4+id)/365.)*10.+0.5
        epoch_g=j*0.1
        call astprs(a1,d1,epoch,a1,d1,epoch_g)
        call astprs(a2,d2,epoch,a2,d2,epoch_g)
        call astprs(a3,d3,epoch,a3,d3,epoch_g)
        call astprs(a4,d4,epoch,a4,d4,epoch_g)
	call astprs(alpha_c/cx,delta_c/cy,epoch,alpha_g,delta_g,epoch_g)
        call ad_azl(a1,d1,az1,al1)
        call ad_azl(a2,d2,az2,al2)
        call ad_azl(a3,d3,az3,al3)
        call ad_azl(a4,d4,az4,al4)
        call ad_azl(alpha_g,delta_g,a1,d1)
        alpha_g=alpha_g*cx
        delta_g=delta_g*cy
        x60=6.
        if(d1.gt.30.)x60=4.
        if(d1.gt.50.)x60=3.
        if(d1.gt.70.)x60=2.
        if(d1.gt.80.)x60=1.
	if(d1.gt.84.)x60=0.5
        if(d1.gt.86.)x60=0.4
        if(d1.gt.88.)x60=0.2
        if(d1.gt.89.)x60=0.1
	call pgsci(3)
	call pgsls(3)
c judge direction of A,H
c center (1024.5,1024.5)--a1,l1
c another point, a2=a1, l2=l1-1
	call azl_ad(a1,d1-1.,x,y)
	call standc(alpha_g,delta_g,x*cx,y*cy,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)
	a3=x-1024.5
	d3=y-1024.5
	x=a3/sqrt(a3*a3+d3*d3)
	x=asin(x)/cy
	if(x.lt.0.)x=-x
	iah=0
	if(x.gt.48.)iah=1
c iah=0:     x direction is (A); y direction is (H)
c iah=1:     x direction is (H); y direction is (A)
	x=amin1(az1,az2,az3,az4)
	y=amax1(az1,az2,az3,az4)
        if(y-x.gt.200.)then
          xa=x
          x=y-360-5./x60
          y=xa+5./x60
        endif
	i1=x*x60
	i2=y*x60+1
	x=amin1(al1,al2,al3,al4)
	y=amax1(al1,al2,al3,al4)
	j1=x*6.
	j2=y*6.+1
	kk=0
	do 163 i=i1,i2
	xa=i/x60
	k=0
	do 162 j=j1,j2
	k=k+1
	xd=j/6.
	call azl_ad(xa,xd,x,y)
	call standc(alpha_g,delta_g,x*cx,y*cy,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)
	x9(k)=x
162	y9(k)=(513.-y)
	call pgline(k,x9,y9)
	xa=xa+0.0001
	j=xa
	k=(xa-j)*60.
	write(ac(1:),'(i3,1x,i2.2)')j,k
	ns=ac(1:3)//'\\uo\\d'//ac(5:6)//char(39)
	if(iah.eq.0)then
	  x=x9(1)-(x9(1)-x9(2))*(y9(1)-512.)/(y9(1)-y9(2))
	  if(x.gt.0..and.x.lt.512.)call pgtext(x-15.,516.,ns)
	else
	  y=y9(1)-(y9(1)-y9(2))*(x9(1)-512.)/(x9(1)-x9(2))
	  if(y.gt.0.and.y.lt.512)call pgtext(512.,y-3.,ns)
  	  if(y.gt.200.and.y.lt.350.and.kk.eq.0)then
	    call pgtext(520.,y+40.,'(A)')
	    kk=1
	  endif
	endif
163	continue
	if(iah.eq.1)call pgtext(490.,525.,'(H)')
	if(iah.eq.0)call pgtext(490.,525.,'(A)')
	do 164 j=j1,j2
	xd=j/6.
	k=0
	do 165 i=i1,i2
	k=k+1
	xa=i/x60
	call azl_ad(xa,xd,x,y)
	call standc(alpha_g,delta_g,x*cx,y*cy,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)
	x9(k)=x
165	y9(k)=(513.-y)

	call pgline(k,x9,y9)
	call toms(xd,dc,0)
	ns=dc(1:3)//'\\uo\\d'//dc(5:6)//char(39)
	if(iah.eq.0)then
	  if(x9(1).gt.x9(2))then
	    do 166 i=1,k
	    if(x9(i).lt.512.)goto 167
166	    continue
	  else
	    do 168 i=1,k
	    if(x9(i).gt.512.)goto 167
168	    continue
	  endif
167	  if(i.gt.k)goto 164
	  y=y9(i)-(y9(i)-y9(i-1))*(x9(i)-512.)/(x9(i)-x9(i-1))
	  if(y.gt.0..and.y.lt.512.)then
	    call pgtext(512.,y-3.,ns)
  	    if(y.gt.200.and.y.lt.350.and.kk.eq.0)then
	      call pgtext(520.,y+50.,'(H)')
	      kk=1
	    endif
	  endif
	else
	  if(y9(1).gt.y9(2))then
	    do 266 i=1,k
	    if(y9(i).lt.0.)goto 267
266	    continue
	  else
	    do 268 i=1,k
	    if(y9(i).gt.0.)goto 267
268	    continue
	  endif
267	  if(i.gt.k)goto 164
	  x=x9(i)-(x9(i)-x9(i-1))*(y9(i)-512.)/(y9(i)-y9(i-1))
	  if(x.gt.0..and.x.lt.512.)call pgtext(x-15.,516.,ns)
	endif
164	continue
c/* draw somthing about moon */
 	call getmxy(a8,x,y,z,phase,pl,pb,xx)
	xa=2100.
	xd=250.
	a1=xa+60.
	d1=xd+60.
	a2=xa+60.
	d2=xd-60.
	a3=xa-60.
	d3=xd-60.
	a4=xa-60.
	d4=xd+60.
	call grlin2(a1,d1,a2,d2)
	call grlin2(a2,d2,a3,d3)
	call grlin2(a3,d3,a4,d4)
	call grlin2(a4,d4,a1,d1)
	if(pb.gt.0)then
	  a1=x*150.+xa
	  d1=y*150.+xd
	  call grlin2(a1,d1,xa,xd)
	  dis='moon(??.?) to ccd '
	  write(dis(6:9),'(f4.1)')phase
	  i=xx+0.5
	  write(dis(19:),'(i3)')i
	  dis=dis(1:21)//'\\uo'
	  call pgtext(440.,-230.,dis)
	else
	  call pgtext(500.,-230.,"no moon")
	endif
        call pgiden
        call pgend
c submit "print map job" to system, in queeue, print map automatically.
	
	end
	
	subroutine star(map,n1,n2,xx,yy,ir,ierr,iturn,sky)
c when ir=4, 15*15 for center
	integer*2 map(n1,n2),maps(-25:25,-25:25)
	real b(-25:25),xmap(-25:25,-25:25)
c sky first
	n25=25
	x=xx
	y=yy
          if(iturn.eq.1)x=513.-x
          if(iturn.eq.1)y=513.-y
          if(iturn.eq.4)x=513.-x
          if(iturn.eq.4)y=513.-y
	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
c	call whitexblack(maps,n51,n51,sky,sigma)
	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		! must minus sky, other wise wrong
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
          if(iturn.eq.1)xx=513.-xx
          if(iturn.eq.1)yy=513.-yy
          if(iturn.eq.4)xx=513.-xx
          if(iturn.eq.4)yy=513.-yy
	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

        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

	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
		
c	parameter (ngsc=400,nstar=300)
c	real*8 a(8)
c	real aa(ngsc),dd(ngsc),gg(ngsc)
c	real xx(nstar),yy(nstar),zz(nstar)
c	A(1)    =  -0.828736234594E-05 
c        A(2)    =  -0.239285576687E-06
c        A(3)    =   0.231124845397E-06 
c        A(4)    =  -0.827521226147E-05 
c        A(5)    =   0.825358366881E-02 
c        A(6)    =   0.872310202561E-02
c        A(7)    =            2.0749257
c        A(8)    =            0.9818901
c	read(13)ng
c	read(13)(aa(i),i=1,ng)
c	read(13)(dd(i),i=1,ng)
c	read(13)(gg(i),i=1,ng)
c	read(12)nc
c	read(12)(xx(i),i=1,nc)
c	read(12)(yy(i),i=1,nc)
c	read(12)(zz(i),i=1,nc)
c	i=kgsc(aa,dd,gg,ng,xx,yy,zz,nc,a)
c	if(i.eq.0)write(*,*)' ok!'
c	end

        subroutine ijpp(a,b,c,n,j)
        real a(n),b(n),c(n)
        n=n-j
        do 10 i=1,n
        a(i)=a(i+j)
        b(i)=b(i+j)
10      c(i)=c(i+j)
        end
	
	function kgsc(aa,dd,gg,ng,xx,yy,zz,nc,xycoef,alpha_c,delta_c,rms)
	parameter (ngsc=400,nstar=300)
	real aa(1),dd(1),gg(1),xx(1),yy(1),zz(1)
c	real*8 cf8(8)
	real wa(ngsc),wd(ngsc),gr(nstar),gd(nstar)
	character*12 ac,dc
c	character*12 ac1,dc1,ns*70,ch*1
	byte xyw(nstar)
	real*8 xixn(2,nstar),xy(2,nstar),xycoef(2,6),adcoef(2,6)
	n12=12
	if(nc.le.30)n12=6
	pi=4.*atan(1.0)
        cx=pi/12.
        cy=pi/180.
	epoch=2000.
	n1=512
	n2=512
	call xytoad(xycoef,adcoef)
	do 100 i=1,ng
        xa=aa(i)*cx
        xd=dd(i)*cy
        call standc(alpha_c,delta_c,xa,xd,xi,xn)
        wa(i)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
100     wd(i)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)

c100	write(*,*)aa(i),dd(i),wa(i),wd(i)
c	write(*,*)xycoef,alpha_c,delta_c
c	call pgbegin(0,'/xw',1,1)
c	call pgenv(1.,512.,512.,1.,1,0)
c	do 101 i=1,ng
c          k=36.-gg(i)
c101     call pgpoint(1,wa(i),wd(i),k)
c	call pgsci(3)
c	do 102 i=1,nc
c          n9=alog(zz(i)/4000.)+20
c          if(n9.lt.20)n9=20
c          if(n9.gt.27)n9=27
c102       call pgpoint(1,xx(i),yy(i),n9)
c	
c	  read(*,*)

	do 105 i=1,nc
105	xyw(i)=0
c	m=0
	do 110 i=1,ng
	x=rms
	do 120 j=1,nc
	z1=wa(i)-xx(j)
	z2=wd(i)-yy(j)
	z=abs(z1)+abs(z2)
	if(z.lt.x)then
	k=j
	x=z
	endif
120	continue
	if(x.ne.rms)then
	xyw(k)=1
	gr(k)=aa(i)
	gd(k)=dd(i)
c	else
c	m=m+1
c	call toms(aa(i),ac1,1)
c	call toms(dd(i),dc1,0)
c		write(*,'(i5,a,2f7.1,2(2x,a))')m,
c     *  ' a gsc star lost',wa(i),wd(i),ac1,dc1
	endif
110	continue

	j=0
	do 152 i=1,nc
	if(xyw(i).ne.0)then
	j=j+1
	xx(j)=xx(i)
	yy(j)=yy(i)
	xy(1,j)=xx(i)	
	xy(2,j)=yy(i)
	gr(j)=gr(i)
	gd(j)=gd(i)
	xyw(j)=1
        xa=gr(j)*cx
        xd=gd(j)*cy
        call standc(alpha_c,delta_c,xa,xd,xi,xn)
        xixn(1,j)=xi
        xixn(2,j)=xn
	endif
152	continue
	n=j
	k=0
	kgsc=0
	if(n.lt.n12)return	! Fail, too less star match
151	k=k+1
	if(n-k.lt.n12)return	! Fail, too less star left
        call plate(xixn,xy,xyw,n,xycoef,3)
        call toms(alpha_c/cx,ac,1)
        call toms(delta_c/cy,dc,0)
c       write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
c	call pgsch(.7)
c	call pgtext(0.,-45.,ns)

c        write(*,2)epoch,ac,dc
c2	format(' **Frame Center:  (',f6.1,')',/5x,2a)
c        write(*,*)(xycoef(1,i),i=1,3)
c        write(*,*)(xycoef(2,i),i=1,3)
        sig=0.
        m=0.
	z=0.
c	call pgsci(6)
        do 150 i=1,n
        xi=xycoef(1,1)*xy(1,i)+xycoef(1,2)*xy(2,i)+xycoef(1,3)
        xn=xycoef(2,1)*xy(1,i)+xycoef(2,2)*xy(2,i)+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,xa,xd)
        xa=xa/cx
        xd=xd/cy
        sigma=((gr(i)-xa)*15.*cos(delta_c))**2+(gd(i)-xd)**2
        sigma=sqrt(sigma)*3600.
        if(xyw(i).eq.1)then
	if(sigma.gt.z)then
	  z=sigma
	  j=i
	endif	
        sig=sig+sigma
        m=m+1
        endif
c        call toms(gr(i),ac1,1)
c        call toms(gd(i),dc1,0)
c        call toms(xa,ac,1)
c        call toms(xd,dc,0)
c        if(i.eq.1)write(*,449)epoch,epoch
c449      format(9x,'GSC (',f6.1,')',12x,'Result (',f6.1,')'
c     *  /2x,24('-'),3x,24('-'))

c       ch=' '
c       if(xyw(i).eq.1)ch='*'
c	write(ns(1:),'(i2)')i
c	call pgsci(2)
c	call pgsch(.7)
c	call pgtext(xx(i)-15.,yy(i)-15.,ns(1:2))
c	call pgsch(1.)
c	if(xyw(i).ne.1)call pgsci(3)
c	call pgpoint(1,xx(i),yy(i),24)
c       write(*,'(i3,2a,2x,2a,g12.2,,4x,a)')i,ac1,dc1,ac,dc,sigma,ch
150	continue
	x=sig/m
c        write(*,'(62x,g12.2)')x
c 3 choice here, after rather good result, must adjest plate center, use 0
c60	write(*,*)'1,2....n  delete or recover this star to re_cal 6 coef.'
c        write(*,*)' 0        calculate center of frame'
c        write(*,*)'-1        quit'
c	read(*,*,err=60)k
	if(x.gt.6.0.and.k.gt.9)return	! Fail, too big error
	if(z.gt.2.0)then
	   xyw(j)=-xyw(j)
	   goto 151
	endif
	k=-k
        xa=(n1+1)*0.5
        xd=(n2+1)*0.5
        xi=xycoef(1,1)*xa+xycoef(1,2)*xd+xycoef(1,3)
        xn=xycoef(2,1)*xa+xycoef(2,2)*xd+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,xa,xd)
        alpha_c=xa
        delta_c=xd
        do 141 i=1,n
        xa=gr(i)*cx
        xd=gd(i)*cy
        call standc(alpha_c,delta_c,xa,xd,xi,xn)
        xixn(1,i)=xi
141     xixn(2,i)=xn
	if(k.lt.0)goto 151
c	cf8(1)=xycoef(1,1)
c	cf8(2)=xycoef(2,1)
c	cf8(3)=xycoef(1,2)
c	cf8(4)=xycoef(2,2)
c	cf8(5)=xycoef(1,3)
c	cf8(6)=xycoef(2,3)
c epoch is 2000
c	cf8(7)=alpha_c
c	cf8(8)=delta_c
c	write(*,*)'GSC:',ng,' CCD:',nc,' Matched:',n
	if(ng.gt.30)write(*,3)n-k,x
3	format(9x,'By using',i4,' stars to calculate 8 ceof, RMS:',
     *  f5.2,' (arcsec)')
	kgsc=n-k
	rms=x
	return
	end

c	parameter (ngsc=3999)		! No. of local area stars.
c	real aa(ngsc),dd(ngsc),gscm(ngsc)
c	character*30 ns
c
c	ep=2000.
c6	write(*,*)'ccd frame center: (ex: 3:12 -23:45 [2000.])'
c	read(*,'(a)')ns
c	if(itohd(ns,a_c,d_c,ep).eq.0)goto 6
c	if(ep.ne.0.)epoch=ep
c
c******** take out gsc stars
c	size=60.
c	call fgsc(a_c,d_c,epoch,size,aa,dd,gscm,ng)
c	call pgbegin(0,'/su',1,1)
c	d1=d_c-0.5
c	d2=d_c+0.5
c	a1=a_c+1.0/15.
c	a2=a_c-0.9/15.
c	call pgenv(a1,a2,d1,d2,0,0)
c	do 10 i=1,ng
c	k=36.-gscm(i)
c10	call pgpoint(1,aa(i),dd(i),k)
c	call pgend
c	write(13)ng
c	write(13)(aa(i),i=1,ng)
c	write(13)(dd(i),i=1,ng)
c	write(13)(gscm(i),i=1,ng)
c	end

	subroutine fgsc(a_c,d_c,e_p,s_z,aa,dd,xmag,n)
	parameter (ngsc=3999)		! No. of local area stars.
	character*6 dir,ns*30			! gsc directory name
	character*8 file1,file2,file3,file4	! gsc file name
	real aa(1),dd(1),xmag(1)	! alpha delta for gsc
	integer*2 ee(ngsc),ingsc(ngsc),e2
	byte w910(2)
	character col(ngsc),color,color1
	logical logi
	equivalence (w910,e2)
	nn=ngsc
	pi=4.*atan(1.)
	cx=pi/12.
	cy=pi/180.
	ns='/EOD/GSC/regions.bin'
	inquire(file=ns,exist=logi)
	if(.not.logi)stop 'GSC cataloge not found !'

	alph_a=a_c
	delt_a=d_c
	if(delta.gt. 89.5)delta= 89.5		! near 90 is too large cal.err
	if(delta.lt.-89.5)delta=-89.5
	call astprs(alph_a,delt_a,e_p,alpha,delta,2000.)

c get 4 corner position of pic. unit is hour & degree
	hw_de=s_z/120.
	up_d=delta+hw_de
	dn_d=delta-hw_de
	hs=hw_de*cy
        hs2=511./hs*0.5
	x=min(cos(up_d*cy),cos(dn_d*cy))
	hw_al=hw_de/x/15.
	if(hw_al.gt.12..or.x.lt.0.)hw_al=12.
	up_a=alpha+hw_al
	dn_a=alpha-hw_al	
c for R.A. to joint 24 hour and 0 hour
c iraf=0 normal, iraf=1 center R.A. is large, iraf=-1 center R.A. near 0 hour
	iaf=0
	if(up_a.ge.24.) iaf=1
	if(up_a.ge.24.) up_a=up_a-24.
	if(dn_a.lt.0.) iaf=-1
	if(dn_a.lt.0.) dn_a=dn_a+24.
	
**************************FIND STAR***********************
	n=0
c if delta>83.5 use all N8230 directory files
	if(abs(delta).gt.83.5)then
	if(delta.gt.0.)then
	dir='N8230/'
	do 20 i=4615,4662		! file name is 4615.BAO~~4662.BAO
	write(file1(1:),"(i4,'.BAO')")i
20	call find(dir,file1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf)
	else
	dir='S8230/'
	do 21 i=9490,9537		! as north
	write(file1(1:),"(i4,'.BAO')")i
21	call find(dir,file1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf)
	endif
	else	
c if delta <83.5 use 4 files or lease.
c according to 4 corner to index which file should be read.
	if(jindex(up_a,up_d,dir,file1).ne.0)goto 30
	call find(dir,file1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf)
c	write(*,'(a)')file1
30	if(jindex(up_a,dn_d,dir,file2).ne.0)goto 31
	if(file2.eq.file1)goto 31
	call find(dir,file2,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf)
c	write(*,'(a)')file2
31	if(jindex(dn_a,up_d,dir,file3).ne.0)goto 32
	if(file3.eq.file1.or.file3.eq.file2)goto 32
	call find(dir,file3,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf)
c	write(*,'(a)')file3
32	if(jindex(dn_a,dn_d,dir,file4).ne.0)goto 33
	if(file4.eq.file1.or.file4.eq.file2.or.file4.eq.file3)goto 33
	call find(dir,file4,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf)
c	write(*,'(a)')file4
33	continue
	endif
c find star end. but in some case a star appear several time in
c different files, find it, keep only one, ignore repeat one.
c GSC says, a star at more appears 7 times in different pltes.
c according to positions' error to detect the repeat one.
	nn=n
c 1. determine how many stars
	cosz=15.*cos(delta*cy)
	n=1
	do 34 i=1,nn-1
	ingsc(i)=n
	x=((aa(i)-aa(i+1))*cosz)**2
	y=(dd(i)-dd(i+1))**2
34	if(x+y.gt.5e-6)n=n+1
c accordind to test, 1e-6 or 1e-5 all OK, so choice middle value
	ingsc(nn)=n
c 2. acording color, change positions of alpha and delta

	color1='V'
	i=0
3410	i=i+1
	if(i.ge.nn)goto 3430
	if(ingsc(i).eq.ingsc(i+1))then
	if(col(i).eq.color1)then
	j=i+1
	x=aa(j)
	y=dd(j)
	k=ee(j)
	color=col(j)
	aa(j)=aa(i)
	dd(j)=dd(i)
	ee(j)=ee(i)
	col(j)=col(i)
	aa(i)=x
	dd(i)=y
	ee(i)=k
	col(i)=color
	endif
3420    if(i.ge.nn-1)goto 3430
        if(ingsc(i+1).eq.ingsc(i+2))then
        i=i+1
        goto 3420
        endif
        endif
        goto 3410
3430    continue

c 3. put each stars in front, mult stars in last 
	do 341 i=1,n
	if(ingsc(i).ne.i)then
	do 342 j=i+1,nn
	if(ingsc(j).eq.i)goto 343
342	continue
343	x=aa(j)
	y=dd(j)
	k=ee(j)
c	color=col(j)
	aa(j)=aa(i)
	dd(j)=dd(i)
	ee(j)=ee(i)
c	col(j)=col(i)
	ingsc(j)=ingsc(i)
	aa(i)=x
	dd(i)=y
	ee(i)=k
c	col(i)=color
	ingsc(i)=i
	endif
341	continue
c 4. delete out of frame star
	k=0
	alph_a=alpha*cx                        ! center of pic.
        delt_a=delta*cy  
	do 344 i=1,n
        xa=aa(i)*cx
        xd=dd(i)*cy
        call standc(alph_a,delt_a,xa,xd,xi,xn)
        ix=hs2*(hs-xi)+1.5
        iy=hs2*(xn-hs)+512.5
        if(inpic(ix,iy).eq.0)goto 344
c release last 2 byte message
        e2=ee(i)
        mag=w910(1)
        if(mag.lt.0)x=mag*0.1-5.
        if(mag.gt.0)x=mag*0.1+5.
        if(mag.lt.0)x=-x
c	if(mag.lt.0)col(i)=col(i)+32
	k=k+1
	aa(k)=aa(i)
	dd(k)=dd(i)
	xmag(k)=x
c	col(k)=col(i)
c	ingsc(k)=ingsc(i)
344	continue
c	write(*,*)nn,k
c 5. 2nd delete same position star
	n=1
	do 400 i=2,k
	do 410 j=1,n
	x=((aa(i)-aa(j))*cosz)**2
	y=(dd(i)-dd(j))**2
	if(x+y.lt.5e-6)then               ! keep bright star
	  if(xmag(j).le.xmag(i))goto 400
	  xmag(j)=xmag(i)
	  aa(j)=aa(i)
	  dd(j)=dd(i)
	  goto 400
	endif
410	continue
	n=n+1
	aa(n)=aa(i)
	dd(n)=dd(i)
	xmag(n)=xmag(i)
400	continue
c	write(*,*)nn,n

	do 500 i=1,n-1
	do 500 j=i+1,n
	if(xmag(i).lt.xmag(j))goto 500
	x=xmag(i)
	xmag(i)=xmag(j)	
	xmag(j)=x
	x=aa(i)
	aa(i)=aa(j)
	aa(j)=x
	x=dd(i)
	dd(i)=dd(j)
	dd(j)=x
500	continue
	if(n.gt.400)n=400	! control gsc number 400, star 300
	end

	subroutine find(dir,f1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf)
c find gsc
c dir file: is cataloge of star
c aa dd ee: 10 byte of gsc star
c nn n :    max array and now size of array
c up.. dn..: 4 corner of chart, in hour & degree
c iaf     : R.A. indicator,  for join 24~0 hour data
	real aa(nn),dd(nn)
	integer*2 ee(nn)
	character*1 col(nn),color,dir1
	character dir*6,f1*8
	integer*2 k512,mag
	byte u(512),v(10,51),w(10)
	equivalence (u,v),(u(511),k512),(alpha,w),(delta,w(5)),(mag,w(9))
c   1    2     3    4     5    6    7    8    9     10
c  ------------------   ------------------  - ---  -- --
c      alpha                 delta          c mag  pe me
c c~classify in first bit, 0~stellar 1~non-stellar
c mag  7 bit, (mag-5)*10
c pe~ position error in ", pe=pe*10, if pe > 15 then pe=0, up 4 bit
c me~ magnitude error in m, me=me*10, if me > 15 then me=0, low 4 bit
	open(1,file='/EOD/GSC/'//dir//f1,status='old',
     *       access='direct',recl=512)
	nrec=0
10	nrec=nrec+1
	read(1,rec=nrec,err=100)u
	do 20 i=1,51
	do 30 j=1,10
30	w(j)=v(j,i)
	call swap4(w,8)
	if(w(9).eq.0)goto 100
c add color identfy
	color='V'
	if(dir1.eq.'S')color='B'
	if(alpha.lt.0.)then
	alpha=-alpha
	color='B'
	if(dir1.eq.'S')color='V'
	endif
	if(delta.gt.90.)then
	delta=delta-90.
	color='R'
	endif
	if(delta.gt.up_d.or.delta.lt.dn_d)goto 20
	if(iaf) 40,50,60
c iaf=-1, ex.   center in 1 hour, dn_a in 23 hour.  24 23...........3 2 1
40	if(alpha.gt.up_a.and.alpha.lt.dn_a)goto 20
		if(alpha.gt.dn_a)alpha=alpha-24.     ! value around 0 pi
	goto 70	
c iaf=0 normal ex.  .........8 7 6 5 .................
50	if(alpha.gt.up_a.or.alpha.lt.dn_a)goto 20
	goto 70
c iaf=1   ex,   center in 23 hour, up_a in 1 hour. 23 22 21.........1 0 
60	if(alpha.gt.up_a.and.alpha.lt.dn_a)goto 20
		if(alpha.lt.up_a)alpha=alpha+24.	! around 2*pi
c bug 96,12,11, fixed
c70	n=n+1                           
c	if(n.gt.nn)return

70	if(n+1.gt.nn)return
	n=n+1
	aa(n)=alpha
	dd(n)=delta
	ee(n)=mag
	col(n)=color
20	continue
	goto 10
100	close (1)
	end

	function jindex(alpha,delta,dir,file)
c given a alpha delta, output a dir & file name
c in this file, the given star must be there.
	character dir*6,file*8
c 1991,1,5
	real a1(9537),a2(9537),d1(9537),d2(9537)
	character*5 cc(24)
	integer*2 ii(24)
	data ii/ 593,1177,1728,2258,2780,3245,3651,4013,
     *          4293,4491,4614,4662,5259,5837,6411,6988,
     *          7522,8021,8463,8839,9133,9345,9489,9537/
	data ip/0/
	data cc/'N0000','N0730','N1500','N2230','N3000','N3730',
     *          'N4500','N5230','N6000','N6730','N7500','N8230',
     *          'S0000','S0730','S1500','S2230','S3000','S3730',
     *          'S4500','S5230','S6000','S6730','S7500','S8230'/
	if(ip.ne.0)goto 99
c read index file of gsc only once
c because slow operate at optical disk, so search file first at /usr/local
	open(1,file='/EOD/GSC/regions.bin',status='old',form='unformatted')
	read(1)a1,a2,d1,d2
	close (1)
	call swap4(a1,9537*4)
	call swap4(a2,9537*4)
	call swap4(d1,9537*4)
	call swap4(d2,9537*4)
99	ip=1
	jindex=0
	do 20 i=1,9537
	if(alpha.lt.a1(i).or.alpha.gt.a2(i))goto 20
	dd1=d1(i)
	dd2=d2(i)
	if(delta.lt.0.)dd1=dd2
	if(delta.lt.0.)dd2=d1(i)
	if(delta.lt.dd1.or.delta.gt.dd2)goto 20
	goto 30
20	continue
	jindex=1
	return
30	write(file(1:),"(i4.4,'.BAO')")i
	do 40 j=1,24
	if(i.le.ii(j))goto 50
40	continue
50	dir=cc(j)//'/'
	end

        function inpic(ix,iy)
c judge ix,iy in frame or not
        inpic=0
        if(ix.lt.1.or.ix.gt.512)return
        if(iy.lt.1.or.iy.gt.512)return
        inpic=1
        end

	subroutine shift_c(gx,gy,aa,dd,n,xycoef,alpha_c,delta_c)
c n===30
	real gx(n),gy(n),aa(n),dd(n)
	real cx(30),cy(30)
	real czr1(30),czr2(30)
	integer*2 w(60)
	integer cz1(30),cz2(30)

	real*8 xycoef(2,6),adcoef(2,6)

	pi=4.*atan(1.0)
        x1=pi/12.
        x2=pi/180.

	call xytoad(xycoef,adcoef)
	do 100 i=1,n
        xa=aa(i)*x1
        xd=dd(i)*x2
        call standc(alpha_c,delta_c,xa,xd,xi,xn)
        cx(i)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
100     cy(i)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)

c	call pgbegin(0,'?',1,1)
c	call pgenv(0.,512.,512.,0.,1,0)
c	call pgpoint(n,cx,cy,5)
c	call pgsci(3)
c	call pgpoint(n,gx,gy,6)
	m=0.
	do 30 i=1,n
	r=500
	cz1(i)=0
	do 40 j=1,n
	rr=sqrt((cx(i)-gx(j))**2+(cy(i)-gy(j))**2)
	if(rr.lt.r)then
          r=rr
	  cz1(i)=j
	  czr1(i)=r
	endif
40	continue
	r=500
	cz2(i)=0
	do 50 j=1,n
	rr=sqrt((cx(i)-gx(j))**2+(cy(i)-gy(j))**2)
	if(rr.lt.r.and.rr.ne.czr1(i))then
          r=rr
	  cz2(i)=j
	  czr2(i)=r
	endif
50	continue
	if(cz1(i).ne.0)then
	  m=m+1
	  w(m)=czr1(i)*10
	endif
	if(cz2(i).ne.0)then
	  m=m+1
	  w(m)=czr2(i)*10
	endif
c	write(*,*)i,czr1(i),czr2(i)
30	continue
	call whitexblack(w,m,1,peak,s)
	peak=peak*0.1
	s=3
c	write(*,*)'peak:  ',peak
	m=0
	do 60 i=1,n
        if(cz1(i).eq.0.and.cz2(i).eq.0)goto 60
70	 if(cz1(i).eq.0)then
	  if(abs(czr2(i)-peak).gt.s)goto 60
	  j=cz2(i)
c	  write(*,*)i,cx(i)-gx(j),cy(i)-gy(j)
	  m=m+1
	  czr1(m)=cx(i)-gx(j)
	  czr2(m)=cy(i)-gy(j)
	  goto 60
	endif
	if(cz2(i).eq.0)then
	  if(abs(czr1(i)-peak).gt.s)goto 60
	  j=cz1(i)
c	  write(*,*)i,cx(i)-gx(j),cy(i)-gy(j)
	  m=m+1
	  czr1(m)=cx(i)-gx(j)
	  czr2(m)=cy(i)-gy(j)
	  goto 60
	endif
	x1=abs(czr1(i)-peak)
	x2=abs(czr2(i)-peak)
	if(x1.gt.x2)cz1(i)=0
	if(x2.gt.x1)cz2(i)=0
	goto 70
60	continue
c	do 80 i=1,m
c80	write(*,*)i,czr1(i),czr2(i)
	x1=xmedian(czr1,m)
	x2=xmedian(czr2,m)
	s=-sqrt(x1*x1+x2*x2)*xycoef(1,1)*10800/pi
	write(*,11)s
11	format(' Auto shift',f5.1,' arc_minutes')
c	call pgend
	alpha_c=alpha_c+x1*xycoef(1,1)
	delta_c=delta_c+x2*xycoef(1,1)
	end

	function xmedian(a,n)
	real a(n)
	do 10 i=1,n-1
	do 10 j=i+1,n
	if(a(i).lt.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
