C simple gsc+usno_sa1.0
c if filed >120 sao only
c if filed >=60 sao+gsc only
c if files < 60 sao+gsc+usno

	parameter (nn=2999)			! for max number of gsc
	parameter (nnsao=999)			! for max number of sao
	parameter (nnsa1=9999)			! for max number of sa1

	character cc*25,dis*73,ns1*32,ch3*3	! read & display line	
	character*6 dir		! gsc directory name
	character*8 file1,file2,file3,file4,file5	! gsc file name
	character*1 se,sm,sg                    ! se,sm,sg ~~for display line
	character*12 ac,dc			! normal format of alpha delta
	character*1 color,c5*5,color1	

	real aa(nn),dd(nn)			! alpha delta for gsc
	integer*2 ee(nn),e2			! gsc mag, err .. (last 2 byte)
	integer*2 ingsc(nn)			! index of gsc for mult objects
	character col(nn)
	integer*2 mapx(nn),mapy(nn)		! position of gsc on pic.

	real aasao(nnsao),ddsao(nnsao)		!             for sao
	integer ensao(nnsao)			! sao number on publication
	byte emsao(nnsao),epsao(nnsao)		! magnitude for sao & sp
	character*3 saosp(202)
	integer*2 msaox(nnsao),msaoy(nnsao)	! position of sao on pic.

	real aasa1(nnsa1),ddsa1(nnsa1)		!             for sa1
	byte ebsa1(nnsa1),ersa1(nnsa1)		! magnitude for b & r color
	integer*2 msa1x(nnsa1),msa1y(nnsa1)	! position of sa1 on pic.


	byte w910(2),bell			! (last 2 byte)

	logical*1 lns,lns1			! plate_center in north = .true.
	logical*1 lsky,logi

	character*3 file

        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2

	equivalence (w910,e2)

	open(1,file='/EOD/sao/sao.sp',status='old',iostat=ierr)
	if(ierr.ne.0)then
	  write(*,*)' /EOD  disk not found !'
	  stop
	endif
	read(1,'(20a)')saosp
	close (1)

	pi=4.*atan(1.)
	cx=pi/12.
	cy=pi/180.
	xlimag=23.

        i=iargc()
        if(i.ne.3)then
	write(*,1001)
1001	format(//' ********* FINDING CHART used by GSC & SAO & USNO_SA1.0
     *  (file) *********'/60x,'2000.7 Jiang'/
     *  ' USNO_SA1.0 is a catalog of 54,787,624 sources over all sky'/
     *  '        data from POSS-I O & E plates  and  UK J & R plates'/
     *  '        astrometric accuracy ~ 0.25 arcseconds'/
     *  '        photometric accuracy   0.25 ~ 0.40 magnitudes'/
     *  ' notice: while finding chart field <10 arcmin, it works'//
     *  '    Usage: gscno_f batch_file output_epoch field'//
     *  '    Ex: gscno_f coord.dat  1997.4  5'//
     *  '    coord.dat include following lines:'/
     *  '    12:34:23.3 -10:30:34.2 1950'/
     *  '    13:00 30:45 1987.2'/
     *  '    3:8:20.34 56:20 1950.0'/
     *  '    ...'//
     *  '  output_file_name is gscno_f.tmp'//) 
        stop
	endif

        call getarg(1,ns1)
	inquire(file=ns1,exist=logi)
	if(.not.logi)then
	  write(*,*)ns1,' file not found !'
	  stop
	endif

	open(1,file=ns1,status='old')
	i=0
2301	i=i+1
	read(1,'(a)',end=2302)dis
	if(itohd(dis,alpha,delta,x).eq.0)goto 2303
	if(x.lt.1900..or.x.gt.2100.)goto 2303
	if(alpha.lt.0..or.alpha.gt.24.)goto 2303
	if(delta.lt.-90.or.delta.gt.90.)goto 2303
	goto 2301

2303	write(*,*)' '
	write(*,*)'something wrong in your batch_file:'
	write(*,'(a,i4,2a)')'line',i,':  ',dis
	stop

2302	close(1)
	open(3,file=ns1,status='old')	  ! chan 3 for read, forever

	call getarg(2,ns1)
	read(ns1(1:),*,err=2304)epoch
	if(epoch.lt.1900.or.epoch.gt.2100.)goto 2304
	goto 2305

2304	write(*,*)' '
	write(*,*)'output_epoch errori: ',ns1
	stop

2305	call getarg(3,ns1)
	read(ns1(1:),*,err=2306)w_de
	if(w_de.lt.1..or.w_de.gt.10.)goto 2306
	goto 2307

2306	write(*,*)' '
	write(*,*)'field error: ',ns1
	stop

2307	open(99,file='gscno_f.tmp',status='unknown')
        call getarg(1,ns1)
	write(99,*)'batch_file: ',ns1
	write(99,*)'field (arcmin): ',w_de
	write(99,*)' '
	igsc=0

10	igsc=igsc+1
	read(3,'(a)',end=2308)ns1
	goto 2309
2308	close(3)
	close(99)
	write(*,*)'gscno_f.tmp is a output file'
	stop
	
2309    if(itohd(ns1,alpha,delta,ep_old).eq.0)goto 10	! char_line  to 3 real
	lsky=.true.
	lns1=.true.
	if(delta.lt.0.)lsky=.false.		! in north sky, lsky=.true.
	call astprs(alpha,delta,ep_old,alph_a,delt_a,epoch)
	call astprs(alph_a,delt_a,epoch,alpha,delta,2000.)
	call toms(alph_a,ac,1)			! chart_epoch alpha in char.
	call toms(delt_a,dc,0)
	write(cc(1:),'(a,1x,a)')ac,dc
	write(*,"(i9,4x,a,' (',f6.1,')')")igsc,cc,epoch
	write(99,"(i9,4x,a,' (',f6.1,')')")igsc,cc,epoch

	lns=.not.lns
c get 4 corner position of pic. unit is hour & degree
	hw_de=w_de/120.
	up_d=delta+hw_de
	dn_d=delta-hw_de

	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.
        a_c=alpha*cx                    ! center of pic.
        d_c=delta*cy                    ! a_c d_c in arc
	hs=hw_de*cy
	hs2=511./hs*0.5
        file(1:1)='N'
        if(delta.lt.0.)file(1:1)='S'
        k=abs(delta)/10.
        k=k*10
        write(file(2:),'(i2.2)')k
	
**************************FIND STAR***********************
********** for gsc
	n=0
c if field > 2 degree, search SAO only.
	if(w_de.gt.120.)goto 1211
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,mapx,mapy)
	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,mapx,mapy)
	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,mapx,mapy)
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,mapx,mapy)
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,mapx,mapy)
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,mapx,mapy)
33	x=(up_a+dn_a)*.5
	if(dn_a.gt.up_a)x=x+12.
	if(x.ge.24.)x=x-24.
	y=(up_d+dn_d)*.5
	if(jindex(x,y,dir,file5).ne.0)goto 332
	if(file5.eq.file1.or.file5.eq.file3)goto 332
	if(file5.eq.file2.or.file5.eq.file4)goto 332
	call find(dir,file5,aa,dd,ee,col,nn,n,mapx,mapy)
332	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.
	if(n.eq.0)goto 1211
	ngsc=n
c 1. determine how many stars
	cosz=15.*cos(delta*cy)
	n=1
	do 34 i=1,ngsc-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(ngsc)=n
c 2. acording color, change positions of alpha and delta
	color1='V'
	if(lns)color1='B'
	i=0
3410	i=i+1
	if(i.ge.ngsc)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
	k=mapx(j)
	mapx(j)=mapx(i)
	mapx(i)=k
	k=mapy(j)
	mapy(j)=mapy(i)
	mapy(i)=k
	endif
3420	if(i.ge.ngsc-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,ngsc
	if(ingsc(j).eq.i)goto 343
342	continue
343	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)
	ingsc(j)=ingsc(i)
	aa(i)=x
	dd(i)=y
	ee(i)=k
	col(i)=color
	ingsc(i)=i
	k=mapx(j)
	mapx(j)=mapx(i)
	mapx(i)=k
	k=mapy(j)
	mapy(j)=mapy(i)
	mapy(i)=k
	endif
341	continue

********* for sao
1211	nsao=0
	call cat_sao(file,
     *    aasao,ddsao,ensao,emsao,epsao,nnsao,nsao,msaox,msaoy)

********* for sa1
	nsa1=0
	i=(up_d+90.)*10.
	i=i/75*75
	write(file1(1:),"('zone',i4.4)")i	
	i=(dn_d+90.)*10.
	i=i/75*75
	write(file2(1:),"('zone',i4.4)")i	
	call cat_sa1(file1,
     *    aasa1,ddsa1,ebsa1,ersa1,nnsa1,nsa1,msa1x,msa1y)
	if(file2.ne.file1)call cat_sa1(file2,
     *    aasa1,ddsa1,ebsa1,ersa1,nnsa1,nsa1,msa1x,msa1y)
**************FIND STAR END****************

	if(n+nsao+nsa1.eq.0)then
	write(99,'(a)')'Sorroy, No Star in this region !'
	goto 10
	endif
 
        num=0
c left button****************
c search SAO (3)
	do 1011 m=1,nsao

	num=num+1
	call astprs(aasao(m),ddsao(m),2000.,xa,xd,epoch)
	xa=amod(xa+24.,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	j=emsao(m)
	xmag=j*0.1
	if(xmag.lt.-1.65)xmag=-xmag
	write(c5(1:),"(f4.1,'v')")xmag
	if(j*.1.lt.-1.65)c5(5:5)='p'
	if(xmag.eq.0.)c5='NoMag'
	j=epsao(m)
	if(j.lt.0)j=j+256
	ch3=saosp(j)
	j=ensao(m)/260000
	k=mod(ensao(m),260000)	
	
1333	format(i3,') ',a,1x,a,'(',f6.1,')  SAO',i7,2(' ',a))
	write(dis(1:),1333)num,ac,dc,epoch,k,c5,ch3
	write(99,'(a)')dis

1011	continue

c search usnosa10
	do 201 m=1,nsa1

	num=num+1
	call astprs(aasa1(m),ddsa1(m),2000.,xa,xd,epoch)
	xa=amod(xa+24.,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	j=ebsa1(m)
        if(j.lt.0)j=j+256
	xbbb=j*0.1
	j=ersa1(m)
        if(j.lt.0)j=j+256
	xrrr=j*0.1
	write(dis(1:),203)num,ac,dc,epoch,xbbb,xrrr
203	format(i3,') ',a,1x,a,'(',f6.1,') USNO_SA1 ',f4.1,'B ',f4.1,'R')
	write(99,'(a)')dis

201	continue

c search gsc next (4)
	do 101 m=1,n

	e2=ee(m)
	mag=w910(1)
	if(mag.lt.0)xmag=mag*0.1-5.
	if(mag.gt.0)xmag=mag*0.1+5.
	if(mag.lt.0)xmag=-xmag
	istar=-1
	if(mag.lt.0)istar=1
	mag=w910(2)
	if(mag.lt.0)mag=256+mag
c position & mag error for use here, only for display
c mark +- in one symbol, is more complicate in write down to file
c so make thing easy, use + only.
	pe=mag/16*0.1
	pm=mod(mag,16)*0.1
	num=num+1
	sg='S'					! star
	if(istar.eq.1)sg='G'			! galaxy
	se='+'
	if(pe.eq.0.)then
	se='>'
	pe=1.5
	endif
	sm='+'
	if(pm.eq.0.)then
	sm='>'
	pm=1.5
	endif
	call astprs(aa(m),dd(m),2000.,xa,xd,epoch)
	xa=amod(xa+24.,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	write(dis(1:),103)num,ac,dc,epoch,se,pe,sg,xmag,col(m),sm,pm
103	format(i3,') ',a,1x,a,'(',f6.1,') ',a,f3.1,'"  ',a,f5.1,a,1x,a,f3.1,'m')
c for n+1 to ngsc
	ipgsc=0
	do i=n+1,ngsc
	if(ingsc(m).eq.ingsc(i))then
	e2=ee(i)
	mag=w910(1)
	if(mag.lt.0)xmag1=mag*0.1-5.
	if(mag.gt.0)xmag1=mag*0.1+5.
	if(mag.lt.0)xmag1=-xmag1
	if(ipgsc.eq.0)then
	write(dis(54:),'(f5.1,a)')xmag1,col(i)
	ipgsc=1
	else
	if(col(m).ne.col(i))write(dis(54:),'(f5.1,a)')xmag1,col(i)
	endif	
	endif
	enddo
	if(ipgsc.eq.1)then
	sm=dis(59:59)
	if((lns.and.sm.eq.'V').or.(.not.lns.and.sm.eq.'B'))then
	c5=dis(49:53)
	dis(49:53)=dis(55:59)
	dis(55:59)=c5
	endif
	endif
	write(99,'(a)')dis

101	continue

	goto 10	
	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 find(dir,f1,aa,dd,ee,col,nn,n,mx,my)
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),mx(nn),my(nn)
	character*1 col(nn),color,dir1
	character*6 dir
	character*8 f1
	integer*2 k512,mag
	byte u(512),v(10,51),w(10)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2
	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,iostat=i)
	if(i.ne.0)goto 100
	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(i_star(alpha,delta,ix,iy).eq.1)then
	j=w(9)
	if(j.lt.0)xmag=j*0.1-5.
	if(j.gt.0)xmag=j*0.1+5.
	if(j.lt.0)xmag=-xmag
        if(xmag.gt.xlimag)goto 20
	 if(n+1.gt.nn)goto 100
	 n=n+1
	 aa(n)=alpha
	 dd(n)=delta
	 ee(n)=mag
	 col(n)=color
	 mx(n)=ix
	 my(n)=iy
	endif
20	continue
	goto 10
100	close (1)
	end

	subroutine cat_sao(f1,aa,dd,en,em,ep,nn,n,mx,my)
c same as sub find, but here for SAO
c sao has no index file like regions.bin of gsc
c epoch here for correct sao star self moving year by year
c w_de > 120, pick out all sao star
	real aa(nn),dd(nn),en(nn),a(4,32)
	byte em(nn),ep(nn),bb(4)
	integer*2 mx(nn),my(nn),ind(256)
	character*3 f1
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2
	equivalence (ee,bb)
	dy=2000.-epoch
	open(1,file='/EOD/sao/'//f1//'.SAO',
     *     status='old',access='direct',recl=512)
	read(1,rec=1)ind
	 call swap2(ind,512)
        x1=dn_a
        x2=up_a
        k1=x1*10.
	if(k1.lt.1)k1=1
        k1=ind(k1)
        i87=0
        if(amax1(abs(up_d),abs(dn_d)).gt.89.)i87=1
        if(i87.eq.1)k1=2
        if(x2.lt.x1)x2=x2+24.
110     do 10 i=k1,ind(256)
        read(1,rec=i)a
        do 20 j=1,32
	 call swap4(a(1,j),12)
        if(a(2,j).gt.x2.and.i87.eq.0)goto 100
        if(i.eq.ind(256).and.a(1,j).eq.0)goto 30
        ee=a(4,j)
	k=bb(1)
	emag=k*0.1
	if(emag.lt.-1.65)emag=-emag
	if(emag.eq.0.)emag=11.
        if(emag.gt.xlimag)goto 20
c        if(w_de.lt.120..and.emag.gt.10.)goto 20
        if(a(3,j).gt.up_d.or.a(3,j).lt.dn_d)goto 20
        if(a(2,j).lt.x1.and.i87.eq.0)goto 20
        if(iaf.eq.1)a(2,j)=a(2,j)+24.
        if(iaf.eq.-1)a(2,j)=a(2,j)-24.
        alpha=a(2,j)
        delta=a(3,j)
	k=bb(3)
        alpha=alpha-dy*k/15.*1e-6
	k=bb(4)
        delta=delta-dy*k*1e-6
        if(i_star(alpha,delta,ix,iy).eq.1)then
	 if(n+1.gt.nn)goto 100
	 n=n+1
	 aa(n)=alpha
	 dd(n)=delta
	 en(n)=a(1,j)
	 em(n)=bb(1)
	 ep(n)=bb(2)
	 mx(n)=ix
	 my(n)=iy
	endif
20      continue
10      continue
30      x1=x1-24
        x2=x2-24
        k1=2
        if(i87.eq.0)goto 110
100     close (1)
        end

	subroutine cat_sa1(f1,aa,dd,eb,er,nn,n,mx,my)
	real aa(nn),dd(nn)
	byte eb(nn),er(nn)
	integer*2 mx(nn),my(nn)
	character*8 f1
	integer*4 a(3)
	real*8 xx
	logical logi
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2
	inquire(file='/USNOSA10/'//f1//'.cat',exist=logi)
	if(.not.logi)return
	open(1,file='/USNOSA10/'//f1//'.cat',
     *     status='old',access='direct',recl=12)
	open(2,file='/USNOSA10/'//f1//'.acc',status='old')

	do 20 k=1,96
	read(2,*)x1,i
	if(x1+0.25.ge.dn_a)goto 21
20	continue
21	close(2)
	x1=dn_a
        x2=up_a
        i87=0
        if(amax1(abs(up_d),abs(dn_d)).gt.89.)i87=1
        if(i87.eq.1)i=1
        if(x2.lt.x1)x2=x2+24.
110     read(1,rec=i,err=30)a
	i=i+1
	  call swap4(a,12)
	  xx=a(1)
	  ra=xx/5400000.
	  xx=a(2)
	  de=xx/360000.-90.
        if(ra.gt.x2.and.i87.eq.0)goto 100
        if(de.gt.up_d.or.de.lt.dn_d)goto 110
        if(ra.lt.x1.and.i87.eq.0)goto 110
        if(i_star(ra,de,ix,iy).eq.1)then
	  if(n+1.gt.nn)goto 100
	  k=a(3)
	  if(k.lt.0)k=-k
  	  j=mod(k,1000)
	  if(j.lt.0)j=j+256
	  if(j*0.1.gt.xlimag)goto 110
	  n=n+1
	  aa(n)=ra
	  dd(n)=de
	  er(n)=mod(k,1000)
  	  k=k/1000
  	  eb(n)=mod(k,1000)
	  mx(n)=ix
	  my(n)=iy
	endif
        goto 110
30      x1=x1-24.
        x2=x2-24.
        i=1
        if(i87.eq.0)goto 110
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
	  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 i_star(alpha,delta,ix,iy)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2
        xa=alpha*cx
        xd=delta*cy
        call standc(a_c,d_c,xa,xd,xi,xn)
        ix=hs2*(hs-xi)+1.5
        iy=hs2*(xn-hs)+512.5
	i_star=inpic(ix,iy)
	end

	subroutine dis_pm(dis)
	character*(*) dis
	character c1*1,c2*1,c3*9,dum*73
	c1=dis(40:40)
	c2=dis(55:55)
	c3='\(2233)'
	dum=dis
	if(c1.eq.'+'.and.c2.eq.'+')then
	dis=dis(:39)//c3(1:8)//dum(41:54)//c3(1:8)//dum(56:)
	else if(c1.eq.'+')then
	dis=dis(:39)//c3(1:8)//dum(41:)
	else if(c2.eq.'+')then
	dis=dis(:54)//c3(1:8)//dum(56:)
	endif
	end
