c 2010.3 for bok
c deal with FLOAT data only.
c check key_word A81 1st,
c  yes: as normal, check !, (redo coord)
c   no: turn data, new fits head, write_down whole data, filename app "u"
c data a is real array, turn_data in "readata" using c_sub

	parameter (n4096=2048)     ! 111:36:01.6 W
	parameter (n512=512)       !  31:57:46.5 N
	                           !  2071       H
	character*60 f1,cc*1
	integer*2    b(n512*n512)
 	real         a(n4096*n4096)
	character*80 head(108)
        integer      nh            ! using head length
        integer      itype         ! 0: no;  1: yes;  2: yes+!
	real*8 a8(8)
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut

	pi=4.*atan(1.0)
	cx=pi/12.
	cy=pi/180.
  	phi=(31.+57./60.+46.5/3600.)*pi/180.	! 31:57:46.5 BOB latitude
  	sphi=sin(phi)
	cphi=cos(phi)
c ! means: force to do a8(8) again 
	k=iargc()
        if(k.eq.0)stop 'usage:  tv_coordlc filename [!] (2k*2k)'
c check file
        call getarg(1,f1)
        if(inqh(f1,n1,n2,head,ibit).eq.0)stop 'DATA not found'
        if(n1.ne.n2.or.n1.ne.n4096)stop 'Not delinha data'
        if(ibit.ne.-32)stop 'deal with float data only!'
c check itype
        nh=indexpos1(head,"END     ",108)
        if(nh.gt.108)nh=108
        i= indexpos1(head,"A81     ",nh)
        itype=1
	if(i.gt.nh)itype=0             ! not found A81
        if(itype.eq.1 .and. k.ge.2)then  
	  call getarg(2,cc)
	  if(cc.eq.'!')itype=2
	endif
        if(itype.eq.1)call dis_5(head,a8)      ! display & get a8[]
        if(itype.eq.0)call newhead(head,nh)    ! produce new fits_head
c cal astronomal parameter
	call getpar(head)			! get s0,iy,im,id,ut
c read data
        nh=72
	call readata(a,b,f1,itype,n1,nh)
c main_job
	call main_job(head,a,a8,itype,b,f1,match,rms,seeing,sky,n1)
	if(itype.eq.1)stop
c write something to fits head
	call put_a8(head,a8,match,rms,seeing,sky)  ! 49-58 write sky seeing,a8
	call put_gr(head)		! 59-63  write gal_co. redden
	call put_xmass(head)		! 15-17, 64 write center,h_angle,airmass
	call put_aa(head)		! 65,66  write Azimuth, Altitude
	call put_moon(head,a8)		! 67-71  phase, position, pos_angle
        call put_ds9(head,a8)           ! 40-48
	i=lnblnk(f1)+1
	f1(i:i)=char(0)
	if(itype.eq.2)then
	write(*,*) ' re_write fits head !'
	call tvcoordc(f1,head)
	endif
	if(itype.eq.0)then
        f1(1:1)='p'
	write(*,'(2a)')' produce new_fits_file: ',f1
	call swap4(a,n1*n1*4)
        call writefits(f1,head,80*72,1)
        call writefits(f1,a,n1*n1*4,0)
	endif	
	end

	subroutine newhead(head,nh)
	character*80 head(1),ohead(108)
	do 10 i=1,nh
10	ohead(i)=head(i)
	do 20 i=6,72
        do 20 j=1,80
20	head(i)(j:j)=char(32)
c        call puthead( 6,head,ohead,"CCDNAME ",nh)
        call puthead( 7,head,ohead,"OBJNAME ",nh)
c        call puthead( 8,head,ohead,"OBSERVER",nh)
c        call puthead( 9,head,ohead,"INSTRUME",nh)
c        call puthead(10,head,ohead,"TELESCOP",nh)
c        call puthead(11,head,ohead,"SITELONG",nh)
c        call puthead(12,head,ohead,"SITELAT ",nh)
c        call puthead(13,head,ohead,"SITEELEV",nh)
c        call puthead(14,head,ohead,"FILTER  ",nh)
c        call puthead(15,head,ohead,"GAIN8   ",nh)
c        call puthead(16,head,ohead,"RDNOISE8",nh)
c        call puthead(17,head,ohead,"DARKCUR ",nh)
c        call puthead(18,head,ohead,"CAMTEMP ",nh)
c        call puthead(19,head,ohead,"DEWTEMP ",nh)
        call puthead(20,head,ohead,"DATE-OBS",nh)
        call puthead(21,head,ohead,"TIME-OBS",nh)
        call puthead(22,head,ohead,"EXPTIME ",nh)
        call puthead(23,head,ohead,"RA      ",nh)
        call puthead(24,head,ohead,"DEC     ",nh)
c        call puthead(25,head,ohead,"EPOCH   ",nh)
c        call puthead(26,head,ohead,"HA      ",nh)
	head(11)(12:23)="111:36:01.6'"
	head(12)(12:23)=" 31:57:46.5'"
	head(13)(12:23)="       2071'"

        call puthead(72,head,ohead,"END     ",nh)
	end

	function indexpos(head,f1,nh)
	character*80 head(1),f1*8
	do 10 indexpos=1,nh       
10	if(head(indexpos)(1:8).eq.f1)return
        write(*,*)'wrong key_word: ',f1
        stop
	end

	function indexpos1(head,f1,nh)
	character*80 head(1),f1*8
	do 10 indexpos1=1,nh     
10	if(head(indexpos1)(1:8).eq.f1)return
	end
		
	subroutine puthead(k,head,ohead,c8,nh)
        character*80 head(1),ohead(1),c8*8
	ipos=indexpos(ohead,c8,nh)
	head(k)=ohead(ipos)
	end

	subroutine getstar(b,starx,stary,starv,nstar,white,seeing)
	real a(2048,2048),b(2048,2048)
	real starx(1),stary(1),starv(1)
	real see400(500)
	n3000=2000
	kt=n3000
5	kk=0
	ktw=kt+white
	do 6 i=1,2048
	do 6 j=1,2048
6	a(i,j)=b(i,j)
	do 10 i=10,4054
        if(i.lt.100.and.kk.gt.20.and.kt.eq.n3000)goto 30
	do 10 j=10,4054
	if(a(i,j).gt.ktw)then
	  call center_9(i,j,a,ktw,ix,iy,iv,n9)
c  	                 if(iv.lt.kt*5)goto 10
	  x=ix
	  y=iy
	  call star(b,2048,x,y,4,ierr,0)
	  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,2048,x,y,white)
	  starv(kk)=iv
            if(iv.lt.4*kt)kk=kk-1
	    call pgsci(8)
	    iv=(iv-2000.)*6./50000.+21
	    if(iv.lt.20)iv=20
	    if(iv.gt.27)iv=27
	    call pgpoint(1,x,y,iv)
	  if(mod(kk,10).eq.1)call dispdot()
	endif
10	continue
11	write(*,*)
	if(kk.lt.40.and.kt.eq.n3000)then
	  kt=kt/4
	  goto 5
	endif
	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=kt*4
	write(*,*)
	goto 5
	end

        subroutine center_9(i1,j1,a,kt,ix,iy,iv,n1)
        integer t1(200)
        integer*2 c1(200)
	real  a(2048,2048)
        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.2048)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

        function see(map,n1,x,y,white)
c (a81+a84)/2 *180*3600/pi = arcsec / pixel
c schmidt is 1.709"/pixel
c bok is 0.45"/pixel
c seeing= 2.355*sigma*1.709
        real map(n1,n1)
        integer b(113)        ! (29-1)*4+1
        ix=x+0.5 
        iy=y+0.5
        j=1
        do 10 i=ix-14,ix+14
        b(j)=map(i,iy)-white
10      j=j+4    
        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-14,iy+14
        b(j)=map(ix,i)-white
110     j=j+4
        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)/8.*2.355*0.45        ! for BOK
        end

	subroutine getmxy(a8,x,y,z,phase,pl,pb,xx)
	real*8 a8(8),b8(12),jd
	common /bok/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 dis_5(head,a8)
	character*80 head(72),ac*12,dc*12
	real*8 a8(8)
	ipos=indexpos(head,'A81     ',72)-1
	do 10 i=1,8
10	read(head(ipos+i)(11:),*)a8(i)
        call xy_rad(1.,1.,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,1)
1	format(5x,'Letf_up corner:  ',35x,'Right_up corner:')
	write(*,2)ac,dc
2	format(a,1x,a,'                          ',$)
        call xy_rad(2048.,1.,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,3)ac,dc
3	format(a,1x,a/28x,'plate center: (2000.0)')
        call xy_rad(2032.5,2032.5,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,4)ac,dc
4	format(26x,a,1x,a/
     * 4x,'Left_down corner:',35x,'Right_down corner:')
        call xy_rad(1.,2048.,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,2)ac,dc
        call xy_rad(2048.,2048.,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,'(a,1x,a)')ac,dc
	end

	subroutine azl_ad(a,h,alpha,delta)	!dd->hd
	common /bok/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 /bok/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)     ! 65,66
	character*80 head(72)
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        ipos=indexpos(head,'HA      ',72)
	read(head(ipos)(19:),'(i2,1x,i2,1x,f5.2)')i,j,x
	t=i+j/60.+x/3600.
	ipos=indexpos(head,'DEC     ',72)
	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     ',72)+14           ! 65 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)		! 66-71
	character*80 head(72)
	real*8 a8(8)
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',72)+16           
        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_ds9(head,a8)                        ! 40-48
	real*8 a8(8)
	character*80 head(72)
	ipos=40
	head(ipos)(1:47)="CTYPE1  = 'RA---ARC'           / DS9 parameters"
	head(ipos+1)(1:32)="CTYPE2  = 'DEC--ARC'           /"
	head(ipos+2)(1:32)='CROTA2  =                      /'
	head(ipos+3)(1:32)='CRVAL1  =                      /'
	head(ipos+4)(1:32)='CRVAL2  =                      /'
	head(ipos+5)(1:32)='CDELT1  =                      /'
	head(ipos+6)(1:32)='CDELT2  =                      /'
	head(ipos+7)(1:32)='CRPIX1  =    2.032500e+03      /'
	head(ipos+8)(1:32)='CRPIX2  =    2.032500e+03      /'
	rac=180./3.141592654
        jpos=indexpos(head,"RA      ",72)
	read(head(jpos)(19:),'(i2,1x,i2,1x,f5.2)')i,j,y
        x=(i+j/60.+y/3600.)*15
	write(head(ipos+3)(11:25),'(e15.6)')x
        write(head(ipos+5)(11:25),'(e15.6)')a8(1)*rac
        jpos=indexpos(head,"DEC     ",72)
	read(head(jpos)(20:),'(i2,1x,i2,1x,f4.1)')i,j,y
        x=i+j/60.+y/3600.
        if(head(jpos)(19:19).eq.'-')x=-x
	write(head(ipos+4)(11:25),'(e15.6)')x
        write(head(ipos+6)(11:25),'(e15.6)')a8(4)*rac
	x=(a8(2)+a8(3))/(a8(1)-a8(4))*rac
	write(head(ipos+2)(11:25),'(e15.6)')x
	end

	subroutine put_a8(head,a8,match,rms,seeing,sky)      !49-58
	character*80 head(72),ac*12,dc*12
	real*8 a8(8)
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=49
	head(ipos  )(1:46)='SKYADU  =                      / (per pixel)  '
        write(head(ipos)(22:30),'(f9.2)')sky
	head(ipos+1)(1:46)='SEEING  =                      / (arcsec)     '
        write(head(ipos+1)(24:30),'(f7.2)')seeing
	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)          !57-58:59-63
	character*80 head(72)
	character*12 ac,dc
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',72)+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,64
	character*80 head(72)
	character*12 ac,dc
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',72)+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      ',72)
	head(ipos)(11:30)="'                  '"
	head(ipos)(18:29)=ac
	write(head(ipos)(50:),"('(',f6.1,')')")epoch
	ipos=indexpos(head,'EPOCH   ',72)
	if(ipos.ne.72)write(head(ipos)(25:30),'(f6.1)')epoch
	ipos=indexpos(head,'DEC     ',72)
	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      ',72)
	head(ipos)(11:30)="'                  '"
	head(ipos)(18:29)=ac
	xmm=xmass(h*cx,xd*cy)				! middle time
	ipos=indexpos(head,'EXPTIME ',72)
	read(head(ipos)(20:),*)xk				! exp time
	xmb=xmass((h-xk/7200.)*cx,xd*cy)			! begin time
	xme=xmass((h+xk/7200.)*cx,xd*cy)			! end time
	x=(xmb+4.*xmm+xme)/6.
	ipos=indexpos(head,'A81     ',72)+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(1)
	real*8 d1,d
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'DATE-OBS',72)
	read(head(ipos)(12:),'(i4,1x,i2,1x,i2)')iy,im,id	! read obs date
	ipos=indexpos(head,'TIME-OBS',72)
	read(head(ipos)(12:),'(i2,1x,i2,1x,f6.3)')i,j,x
	ipos=indexpos(head,'EXPTIME ',72)
	read(head(ipos)(20:),*)xk				! exp time
	ut=i+j/60.+x/3600.+xk/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=-(111.+36.0/60.+1.6/3600.)/15.       ! 111:36:01.61 lamda of BOK
        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
        call disfithead(head,'OBJNAME ')
        call disfithead(head,'DATE-OBS')
        call disfithead(head,'TIME-OBS')
        call disfithead(head,'EXPTIME ')
        call disfithead(head,'RA      ')
        call disfithead(head,'DEC     ')
c        call disfithead(head,'HA      ')
	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 /bok/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 /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',72)
	if(ipos.eq.109)ipos=indexpos(head,'        ',72)
	if(ipos.gt.51)ipos=indexpos(head,'COMMEMT ',72)
	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 /bok/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,itype,n,nh)
	real a(2048,2048),c(2048,2048)
	integer*2 b(512,512)
	character*60 f1
        leng=n*n*4
	ipos=nh*80
        i=lnblnk(f1)+1
        f1(i:i)=char(0)
	call tvread(f1,a,leng,ipos)
	call swap4(a,leng)
        if(itype.eq.0)then
	  n1=n+1
	  do 10 i=1,n
	  do 10 j=1,n
10	  c(i,j)=a(n1-j,n1-i)
	  do 20 i=1,n
	  do 20 j=1,n
20	  a(i,j)=c(i,j)
	endif	
 	call shrink2(a,n,b)		! 2048-->512
	end

	subroutine shrink2(a,n,b)
	real a(n,n)
	integer*2 b(1)
	k=0
	m=4
	do 10 j=1,n,m
	do 10 i=1,n,m
	k=k+1
	x=a(i,j)
	if(x.gt.32767.)x=32767.
10	b(k)=x+0.5
	end

	function inqh(f1,n1,n2,head,ibit)
	character*(*) f1
	character*80 head(108)
	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=2880*3)
	read(1,rec=1)head
	read(head(2)(22:),*)ibit
	read(head(4)(22:),*)n1          	!NAXIS1
	read(head(5)(22:),*)n2          	!NAXIS2
        close(1)
        end

        subroutine disfithead(head,h8)
        character*80 head(1), h8*8
	ipos=indexpos(head,h8,72)
	write(*,'(1x,a)')head(ipos)(1:70)
	end

	subroutine main_job(head,map,a8,
     c            itype,ma,f1,match,rms,seeing,white,n1)
	character*80 head(72)
	real map(2048,2048)
	integer*2 ma(512,512)	! map buffer
	integer*2 dummy(301*301)
	real*8 a8(8)
	character*12 ac,dc,ac1,dc1		! alpha delta display line
	character*1 ch,old_ch,ch3*3,sg,c3*3	! work char variable
	character*50 ns,dis*60,mns*30		! for display
	byte bell,xyw(99)			! weight of standard posi star
	real x9(99),y9(99),xm(99)		! al_del of starard posi star
	real sx(5),sy(5)			! scale on right_bottom of pic
	real*8 xixn(2,99),xy(2,99),dummy8	! x_y of standard posi star
	real*8 adcoef(2,6),xycoef(2,6)		! 2 set 6 plte coef
	real starx(500),stary(500),starv(500)   ! take ccd stars
	parameter (ngsc=3999)		        ! No. of gsc stars.
	real aagsc(ngsc),ddgsc(ngsc),mmgsc(ngsc)
	character f1*30
	logical logi
	common /bok/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ij=0
	do 1089 i=100,400
	do 1089 j=100,400
	ij=ij+1
1089	dummy(ij)=ma(i,j)
        call whitexblack(dummy,301,301,white,x)
	black=white+20.*x
 	call pgbegin(0,'/xw',1,1 )
	call pgpap(6.5,1.)
	call pgsci(0)
        call pgenv(1.,float(n1),1.,float(n1),1,-1)
	call pgsci(1)
        call pgxgray(ma,512,512,1,512,1,512,white-5*x,black,2)
	call pglabel('R.A. south','Decl. east',f1)

	jy=1
	bell=7
	ir=4
	kexit=0
	down=-120.
	call pgsch(.8)
	call pgsci(2)				! standard star color mark
	n=0					! display star serial number
	call unlink('coord.tmp')
9	if(itype.eq.0)write(*,91)bell
91	format(a1,'use left_button to pick standard position star',
     *            ' (another window, run "gsc")')
c put cursor in center of picture
	x=(n1+1)*.5
	y=x
	if(itype.eq.1)then
	  ipos=indexpos(head,'A81     ',72)
	  read(head(ipos)(11:),*)xycoef(1,1)
	  read(head(ipos+1)(11:),*)xycoef(2,1)
	  read(head(ipos+2)(11:),*)xycoef(1,2)
	  read(head(ipos+3)(11:),*)xycoef(2,2)
	  read(head(ipos+4)(11:),*)xycoef(1,3)
	  read(head(ipos+5)(11:),*)xycoef(2,3)
	  read(head(ipos+6)(21:),*)alpha_c
	  read(head(ipos+7)(21:),*)delta_c
	  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)
	  a8(7)=alpha_c
	  a8(8)=delta_c
88	  epoch=2000.0
	  write(*,89)bell,epoch
89	  format(1x,a,'output epoch, [',f6.1,'] ',$)
	  read(*,'(a)')ns
	  if(ns(1:1).ne.' ')then
	    read(ns(1:),*,err=88)epoch
	    if(epoch.lt.1900.or.epoch.gt.2100.)goto 88
	  endif
	  call astprs(alpha_c/cx,delta_c/cy,2000.,x,y,epoch)	
	  call toms(x,ac,1)
	  call toms(y,dc,0)
	  write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
	  goto 1001
	endif

10      call pgcurse(x,y,ch)
	if(ch.eq.'A')then
	call star(map,n1,x,y,ir,ierr)
c x,y is input & output, x,y is more precision when return	
c "ir" is search radius & star rdius, if ir=3, take out 7*7 pixels
c if no star, and low S/N (<5) then ierr=1
	if(ierr.ne.0)then
	write(*,'(2a)')bell,'No star in this region'
	goto 10
	endif
	n=n+1
	xy(1,n)=x
	xy(2,n)=y
	xyw(n)=1
	write(ch3(1:),'(i3)')n
	call pgtext(x-9,y+9.*jy,ch3)
	endif
	if(ch.eq.'X')goto 20			! pick up standard star end.
	goto 10
20	if(n.lt.3)write(*,*)'too few stars to calaulate 6 plate ceof.'
	if(n.lt.3)goto 9
c use "GSC" results, star order of two system must the same
******* read gsc.tmp begin
21	open(3,file='gsc.tmp',status='old',iostat=ierr)
	if(ierr.ne.0)then
	write(*,*)'I can not find "gsc.tmp" or something wrong !'
	write(*,*)'when ready, enter "return" please'
	read(*,'(a)')ch
	goto 21
	endif
30	read(3,1,end=40,err=430)i,i1,i2,xa,ch,j1,j2,xd,epoch,x
1	format(i3,3x,i2,1x,i2,1x,f5.2,1x,a,i2,1x,i2,1x,f4.1,2x,f6.1,10x,f5.1)
	if(i.gt.n)goto 30
	x9(i)=i1+i2/60.+xa/3600.
	y9(i)=j1+j2/60.+xd/3600.
	if(ch.eq.'-')y9(i)=-y9(i)
	xm(i)=x
	goto 30
430	read(3,'(a)',end=40)ch
	goto 30
c assume center of ccd_frame is fisrt star
40	alpha_c=x9(1)*cx
	delta_c=y9(1)*cy
	close (3)
42	do 41 i=1,n
	xa=x9(i)*cx
	xd=y9(i)*cy
	call standc(alpha_c,delta_c,xa,xd,xi,xn)
	xixn(1,i)=xi
41	xixn(2,i)=xn
*********** read gsc end
c calculate 6 plate coef.
43		write(*,*)'****************'
	call plate(xixn,xy,xyw,n,xycoef,3)
c display 6 coef and ...
	write(*,2)epoch,ac,dc
2	format('Frame center:  (',f6.1,')'/5x,2a)
	write(*,*)(xycoef(1,i),i=1,3)
	write(*,*)(xycoef(2,i),i=1,3)
	sig=0.
	xmaxsig=0.
	m=0.
	do 50 i=1,n
	if(xyw(i).eq.0)goto 50
        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=((x9(i)-xa)*15.*cos(delta_c))**2+(y9(i)-xd)**2
	sigma=sqrt(sigma)*3600.
	if(xyw(i).eq.1)then
	  sig=sig+sigma
	  if(sigma.gt.xmaxsig)then
	    j=i
	    xmaxsig=sigma
	  endif
	  m=m+1
	endif
	call toms(x9(i),ac1,1)
	call toms(y9(i),dc1,0)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	if(i.eq.1)write(*,49)epoch,epoch
49	format(9x,'GSC (',f6.1,')',12x,'Result (',f6.1,')'
     *  /2x,24('-'),3x,24('-'))
	ch=' '
	if(xyw(i).eq.1)ch='*'
50	write(*,'(i3,2a,2x,2a,g12.2,4x,a)')i,ac1,dc1,ac,dc,sigma,ch
	write(*,'(62x,g12.2)')sig/m
	write(*,'(a,i5,g12.2)')'Max_sigma:',j,xmaxsig
***************
c 3 choice here, after rather good result, must adjest plate center, use 0
60	write(*,*)'1,2....n  delete or recover this star to re_cal 6 coef.'
	write(*,*)' 0        calculate center of frame'
	write(*,*)'-1        quit'
	if(kexit.eq.-1)goto 90
	read(*,*)k
	if(k.eq.-2)goto 90
	if(k.eq.-1)then
	  kexit=-1
	  k=0
	endif
	if(k.eq.0)then
	xa=(n1+1)*0.5
	xd=(n1+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
	goto 42
	endif
	if(k.ge.1.and.k.le.n)then
	xyw(k)=-xyw(k)
	goto 43				! re_calculate 6 plate
	endif
	goto 60
******** auto match
90	write(*,*)'wait, auto match stars. Now, pick out CCD stars'
	nstar=500
	call getstar(map,starx,stary,starv,nstar,white,seeing)
c******** take out gsc stars
	call fgsc(alpha_c/cx,delta_c/cy,epoch,40.,aagsc,ddgsc,mmgsc,ng)
	call astprs(alpha_c/cx,delta_c/cy,epoch,xa,xd,2000.)
	epoch=2000.
	alpha_c=xa*cx
	delta_c=xd*cy
	match=kgsc(aagsc,ddgsc,mmgsc,ng,
     *	       starx,stary,starv,nstar,
     *		xycoef,alpha_c,delta_c,rms)
	if(match.eq.0)stop 'match fail, do it again !'
********* 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)
	a8(7)=alpha_c
	a8(8)=delta_c			! 2000 epoch
	call toms(alpha_c/cx,ac,1)
	call toms(delta_c/cy,dc,0)
	write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
********* save a8 end
1001	call xytoad(xycoef,adcoef)	! 2 set of 6 coef convert
	xycoef(1,4)=alpha_c
	xycoef(2,4)=delta_c
c draw scale	
	y=(abs(xycoef(1,1))+abs(xycoef(2,2)))*.5
	x=y*180.*60./pi*n1			! X size in minute
	ii=nint(x/4.)
	k=nint(ii/x*n1)
	k4=k/4
	sx(1)=n1-90
	sx(2)=sx(1)
	sx(3)=sx(1)-k
	sx(4)=sx(3)
	sy(1)=90
	sy(2)=sy(1)-30*jy
	sy(3)=sy(2)
	sy(4)=sy(1)
	write(c3(1:),"(i2)")ii
	c3(3:3)=char(39)
        call pgsci(2)
	call pgline(4,sx,sy)
	call pgtext(n1-280.,sy(1),c3)
	call pgsci(5)
	call pgtext(900.,n1+90.,ns)
********************************************* loop
	num=0
	ip=0				! indicator, show stand star again
	open(99,file='coord.tmp',status='unknown')
	ep_old=epoch
	
100	write(*,*)'  left button -- pick out a star'
	write(*,*)'middle button -- put a false_star'
	write(*,*)' right button -- quit'
	write(*,*)'  "n" or "N"  -- input a false_star  (in alpha delta)'

	if(ip.eq.0)then
c display standard star again in display line position
c another purpose is for serial number
	do 101 i=1,n
	x=xy(1,i)
	y=xy(2,i)
	call star(map,n1,x,y,ir,ierr)
	if(ierr.ne.0)goto 101
        xi=xycoef(1,1)*x+xycoef(1,2)*y+xycoef(1,3)
        xn=xycoef(2,1)*x+xycoef(2,2)*y+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
	num=num+1
	xa=amod(xa,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	sg=' '
	if(xyw(i).eq.1)sg='*'
c in display, * stand for standard star, ? stand for saturated star
	write(*,1335)num,sg,ac,dc,epoch
1335	format(i3,')',2a,1x,a,'(',f6.1,')')
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),1335)num,sg,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
101	continue
	call pgsci(13)
	call pgebuf
	endif
	ip=1			! only do once, display at "help" tail
	x=2032.5
	y=2032.5
***** loop in here many time,
110	old_ch=ch		! for zoom, save last command.
	call pgcurse(x,y,ch)

c left button****************
	if(ch.eq.'A')then
	call pgsci(6)
	oldx=x
	oldy=y
	call star(map,n1,x,y,ir,ierr)	! find star
	call star(map,n1,x,y,ir,ierr)	! find star
	if(ierr.eq.1)goto 200
	xi=xycoef(1,1)*x+xycoef(1,2)*y+xycoef(1,3)
        xn=xycoef(2,1)*x+xycoef(2,2)*y+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
c for display ony
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgtext(x-9,y+9.*jy,ch3)
	xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),1334)num,ac,dc,epoch
1334	format(i3,') ',a,1x,a,'(',f6.1,')')
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis(:40),x,y
	goto 110	
	endif

c right button****************
	if(ch.eq.'X')then
	  close (99)
	  call pgend
	  return
	endif

c middle  button**************
	if(ch.eq.'D')then
	call pgsci(4)
        xi=xycoef(1,1)*x+xycoef(1,2)*y+xycoef(1,3)
        xn=xycoef(2,1)*x+xycoef(2,2)*y+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
	xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	num=num+1
	write(ch3(1:),'(i3)')num
	 call pgpoint(1,x,y,16)
	 call pgtext(x-9,y+9.*jy,ch3)
104	format(i3,') ',a,1x,a,'(',f6.1,') false star')
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),104)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	goto 110
	endif

c file  star******************
	if(ch.eq.'f'.or.ch.eq.'F')then
	ep_new=2000
	write(*,'(a)')' Input file name (epoch=2000): '
	read(*,'(a)')mns
	open(98,file=mns,status='old')
9120	read(98,'(a)',end=91200)mns
	write(*,*)mns
        ii=1
	if(mns(1:1).eq.' ')ii=2 
	if(itohd(mns(ii:),x_a,x_d,x).eq.0)i=i+1
	call astprs(x_a,x_d,ep_new,xa,xd,2000.)
	x_a=xa*cx
	x_d=xd*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)
	if(x.lt.0.or.y.lt.0.or.x.gt.n1.or.y.gt.n2)then
	  write(*,'(2a)')bell,'out of the frame'
	  goto 9120
	endif
	xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgsci(2)
	  call pgpoint(1,x,y,-1)
	  call pgtext(x-9,y+9.*jy,ch3)
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),105)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	goto 9120
91200	close(98)
	endif

c new star******************
	if(ch.eq.'n'.or.ch.eq.'N')then
	ep_new=ep_old
120	write(*,1201)ep_new
1201	format(' Input new star: alpha delta [',f6.1,']')
	read(*,'(a)')mns
	if(mns(1:1).eq.' ')goto 110
	if(itohd(mns,x_a,x_d,x).eq.0)goto 120	
	if(x.ne.0.)ep_new=x
	call astprs(x_a,x_d,ep_new,xa,xd,2000.)
	x_a=xa*cx
	x_d=xd*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)
	if(x.lt.0.or.y.lt.0.or.x.gt.n1.or.y.gt.n2)then
	  write(*,'(2a)')bell,'out of the frame'
	  goto 120
	endif
	xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgsci(2)
	  call pgpoint(1,x,y,18)
	  call pgtext(x-9,y+9.*jy,ch3)
105	format(i3,') ',a,1x,a,'(',f6.1,') input alpha delta')
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),105)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	write(*,"(' Again ? [Y] ',$)")
	read(*,'(a)')sg
	if(sg.ne.'n'.and.sg.ne.'N')goto 120	
	endif

200     if(ch.ne.'A'.and.ch.ne.'D'.and.ch.ne.'X')write(*,'(a)')bell
        goto 110
        close (99)
        call pgiden
        call pgend
        write(*,'(1x,a)')bell
        return
	end
	
	subroutine star(map,n1,xx,yy,ir,ierr)
c when ir=4, 15*15 for center
	real map(2048,2048)
	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.n1-n25)ky=n1-n25
	do 10 j=-n25,n25
	do 10 i=-n25,n25
10	maps(i,j)=map(kx+i,ky+j)	
	n51=n25+n25+1
	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.n1)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.n1)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

        subroutine xy_rad(x,y,a8,alpha,delta)
c  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 kgsc(aa,dd,gg,ng,xx,yy,zz,nc,xycoef,alpha_c,delta_c,x)
	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=2048
	n2=2048
	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)

          call pgsci(2)
	  call pgpoint(ng,wa,wd,24)
	  call pgsci(4)
	  call pgpoint(nc,xx,yy,26)

	do 105 i=1,nc
105	xyw(i)=0
c	m=0
	do 110 i=1,ng
	x=15.
	do 120 j=1,nc
	z1=wa(i)-xx(j)
	z2=wd(i)-yy(j)
c	z=abs(z1)+abs(z2)
	z=sqrt(z1*z1+z2*z2)
	if(z.lt.x)then
	k=j
	x=z
	endif
120	continue
	if(x.ne.15.)then
	xyw(k)=1
	gr(k)=aa(i)
	gd(k)=dd(i)
	endif
110	continue

	j=0
	do 152 i=1,nc
	if(xyw(i).ne.0)then
	  call pgpoint(1,xx(i),yy(i),25)
	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)
        sig=0.
        m=0.
	z=0.
        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
150	continue
	x=sig/m
	if(x.gt.6.0.and.k.gt.9)return	! Fail, too big error
	if(z.gt.1.5)then                ! control precision
	   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
	write(*,*)'GSC:',ng,' CCD:',nc,' Matched:',n
	write(*,3)n-k,x
3	format(9x,'By using',i4,' stars to calculate 8 ceof, RMS:',
     *  f5.2,' (arcsec)')
	kgsc=n-k
	return
	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*1 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)
	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)
	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)
	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)
	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
	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)
	dir1=dir(1:1)
	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
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='/usr/local/regions.bin',status='old',
     *       form='unformatted',iostat=ierr)
	if(ierr.ne.0)then
	close (1)
	open(1,file='/EOD/GSC/regions.bin',
     *             status='old',form='unformatted')
	endif
	read(1)a1,a2,d1,d2
	call swap4(a1,9537*4)
	call swap4(a2,9537*4)
	call swap4(d1,9537*4)
	call swap4(d2,9537*4)
	close (1)
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

        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
