parameter (nn=4999)			! for max number of gsc
	parameter (nnsao=1999)			! for max number of sao

	character cc*25,ns*35,dis*73,ns1*32	! read & display line	
	character*6 dir				! gsc directory name
	character*8 file1,file2,file3,file4,file5	! gsc file name
	character*12 ac,dc			! normal format of alpha delta
	character*1 color,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.


	byte w910(2),bell			! (last 2 byte)
	logical 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.
	bell=7
	ep_old=1950.0				! always be 1950.0
	w_de=60.
	xlimag=17.
        icall=0
           
	epoch=2000.
	ns1='   '
	ns='/EOD/GSC/N0000/0001.BAO'
	isa=1					! switch sao star
	isv=1					! switch gcv star
	icd=0					! if GSC mounted
	inquire(file=ns,exist=logi)
	if(.not.logi)then
	write(*,*)' No "GSC" DISK !'
	write(*,*)
	icd=1
	endif

	n=iargc()
        if(n.lt.2)then
        write(*,*)' 2~3 parameters please ! (2008,10 jiang)'
	write(*,*)' ex. g0 10:20 30 [1950]'
        write(*,*)'    (output epoch 2000, ~ 1 square degree objects) '
        stop
        endif
        call getarg(1,ac)
        call getarg(2,dc)
        call getarg(3,ns1)
        dis=ac(1:lnblnk(ac)+1)//dc(1:lnblnk(dc)+1)
        if(n.eq.3)dis=ac(1:lnblnk(ac)+1)//dc(1:lnblnk(dc)+1)//ns1
        write(*,'(1x,a)')dis
        if(itohd(dis,alpha,delta,ep_old).eq.0)then
          write(*,*)' input error'
          stop
        endif
        if(ep_old.eq.0)ep_old=1950.
        call toms(alpha,ac,1)
        call toms(delta,dc,0)
        write(ns1(1:),'(a,1x,a,f7.1)')ac,dc,ep_old

        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(*,"(9x,a,' (',f6.1,')')")cc,epoch


	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..or.icd.eq.1)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)
	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,mapx,mapy)
	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,mapx,mapy)
	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,mapx,mapy)
	write(*,'(a)')file4
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)
	write(*,'(a)')file5
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'
	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)
**************FIND STAR END****************

ccc for gsc
        open(97,file='posi0.dat',status='unknown')
	num=0 
	do 50 i=1,n
        xa=aa(i)*cx
        xd=dd(i)*cy

	e2=ee(i)
	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
c        if(xmag.gt.lim_mag)goto 50
        write(97,'(2f9.5,f6.2,2x,a)')xa,xd,xmag,col(i)
        num=num+1

50	continue
        write(*,*)'Total GSC',num

ccc for sao
        numsao=0
        do 1050 i=1,nsao

        xa=aasao(i)*cx
        xd=ddsao(i)*cy
        numsao=numsao+1
	j=emsao(i)
	xmag=j*0.1
	if(xmag.eq.0.)xmag=11.
	if(xmag.lt.-1.65)xmag=-xmag
        write(97,"(2f9.5,f6.2,'  S')")xa,xd,xmag
1050	continue
	write(*,*)'Total SAO',numsao

        write(*,*)' ALL ',num+numsao
	write(*,*)' output file is posi0.dat'
        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
	 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

	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

	subroutine pgcircle(x,y,r)
	real xx(64),yy(64)
	pi=3.1416
	i=1
	do 10 z=0.,6.3,0.1
	xx(i)=r*cos(z)+x
	yy(i)=r*sin(z)+y
10	i=i+1
	call pgline(64,xx,yy)
	end