c  GIF 640*480,  multify 0.77
        parameter (nnsao=1999)			! for max number of sao
	parameter (m695=720)
	parameter (m512=528)                  ! 720*0.734)
	
	character cc*25,dis*73			! read & display line	
	byte bis(73)
	character*1 ch				! ch~~ job indicator
	character*12 ac,dc,tfile		! normal format of alpha delta
	character*17 xc				! line for deal with default

	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.

	integer pgbegin
	integer*2 map(m695,m512),map33(33,33)	! m695*m512 map & star map
	real sx(5),sy(5)
	integer q(13)				! chinese character array

	byte w910(2)		         	! (last 2 byte)

	character*3 s1(72),c6*6,sn(29)
	integer d3(3,72)
	byte b6(6)
	integer*2 rec1(72),rec2(72)

	integer than(30),phan(30),q6(6),q4(4),q5(6)
        integer s3,d18,s4,i3
        character*55 ns
	real xx(4999),yy(4999)
	real xt(4999),yt(4999)
	integer c1,c2,c3,c4

        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512

	equivalence (w910,e2)
	equivalence (c6,b6)
	equivalence (dis,bis)

	data q/1964,3013,308,3647,309,1964,4619,308,2211,309,308,1717,309/
	data q6/0101,2981,4476,2136,4627,5435/
	data q5/0308,3531,4476,5027,2167,0309/
	data q4/5243,4848,2329,4782/

        k=iargc()
	if(k.lt.1)then
	  write(*,1001)
1001	  format(//'             *********** Minor trace (used SAO)
     *    *********'//30x,'640*480 gif               2012.6 Jiang'/)
	  write(*,*)
	  write(*,*)'     Usage: sao8 file [n] [delay,minor_name_len]'
          write(*,*)'   n='
          write(*,*)'  none   screen, all display  1024*768_GIF'
          write(*,*)'     0   screen, trace only   1024*768_GIF'
          write(*,*)'   +-1   screen | gif, minor fixed, 640*480_GIF'
          write(*,*)'   +-2   screen, movie. "delay" works, default 5'
          write(*,*)'   +-3   gif, 1 pic/day             640*480_GIF'
	  write(*,*)'   +-4   gif, 2 pic/day, "-" means, add grid...'
          write(*,*)'   +-5   gif, 4 pic/day'          
          write(*,*)'   if "file" not exites, show 3 data_files struct'
          write(*,*)
          stop
	endif
        s3=-1
        s4=0
        d18=18
        i3=1
        if(k.gt.1)then
          call getarg(2,xc)
          read(xc(1:),*)s3
          if(s3.lt.0.)then
            s3=-s3;
            i3=-1;
          endif
          if(s3.gt.0)d18=1
        endif
        if(s3.lt.0)i3=-1
        if(k.gt.2)then
          call getarg(3,xc)
          read(xc(1:),*)s4
        endif
	call getarg(1,xc)  
	open(11,file=xc,status='old',iostat=ierr)
        if(ierr.ne.0)then
          write(*,*)
          write(*,*)'**ERROR:  ',xc,' input file not found!'
          write(*,*)
          k=lnblnk(xc)
          j=index(xc(1:k),'.')-1
          if(j.lt.0)j=k
        write(*,*)'prepare 2 files, named ',xc(1:k)," & ",xc(1:j),'.coo'
          write(*,*)'------------------'
          write(*,*)xc(1:k),":    include following lines"
        write(*,*)'09:19:38.47 +20:29:11.8 (2000.0)       ! minor coord'
        write(*,*)'08:45:00.00 +00:30:00.0 (2000.0)       ! pic coord'
        write(*,*)'chinese_word,no space,begin with name  ! minor name'
        write(*,*)'ascii_code also use chinese_sytle      ! other mess'
        write(*,*)'260.,56.,210.,52.0, 0, 0       ! last 2: ra,dec posi'
        write(*,*)xc(1:j),'.coo                       ! data filename'
          write(*,*)'------------------'
          write(*,*)xc(1:j),'.coo:   include folowing lines'
        write(*,*)'1997-May-25 16:00  m  04 41 12.07 +21 04 49.3'
        write(*,*)'1997-May-26 16:00  m  04 42 59.63 +21 09 28.9'
          write(*,*)'...'
	  write(*,*)'------------------'
	  write(*,*)'/vega2/pgplot/color shows color_table'
	  write(*,*)'if want to change starname color, edit color.dat'
	  write(*,*)'11           ! minor chinese_code color (0~15)'
	  write(*,*)'5            ! minor mark color         (0~15)'
          write(*,*)'8            ! sao_star color           (0~15)'
	  write(*,*)'7            ! date                     (0~15)'
          write(*,*)
          stop
        endif
	c1=11
	c2=5
	c3=8
	c4=7
	open(1,file='color.dat',status='old',iostat=ierr)
        if(ierr.eq.0)then
	  read(1,*)c1
	  read(1,*)c2
	  read(1,*)c3
	  read(1,*)c4
	endif
 	close(1)
        write(*,*)'color: ',c1,c2,c3
	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)

	open(20,file='/EOD/sao/sao.nam',status='old',access='direct',recl=8)
	
	pi=4.*atan(1.)
	cx=pi/12.
	cy=pi/180.
	w_de=1500.                       ! 1500/900  5/3
	xlimag=12.
           
        s737=0.734                       ! for GIF 640*480 output, 
	s1357=1./s737
        s10152=8.01                    ! pgpap(10.152,0.737)
	x695=m695
	x512=m512

	read(11,'(a)')dis                ! read minor coord
	call itohd(dis,alpha,delta,x)
        a_c1=alpha*cx                    ! monir position
        d_c1=delta*cy                    ! a_c d_c in arc
	read(11,'(a)')dis                ! read pic_center coord
	call itohd(dis,alpha,delta,x)    ! plot center
	
	if(delta.gt. 89.5)delta= 89.5		! near 90 is too large cal.err
	if(delta.lt.-89.5)delta=-89.5
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(cc(1:),'(a,1x,a)')ac,dc
	write(*,'(a)')cc

	i=w_de*0.25/60.+0.5
	sx(3)=640.-i*60/w_de*m512	
	sx(4)=sx(3)
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.*s1357
	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=(x512-1.)/hs*0.5
	xhs=hs*s1357
	
**************************FIND STAR***********************
	nsao=0
	call cat_sao(
     *    aasao,ddsao,ensao,emsao,epsao,nnsao,nsao,msaox,msaoy)
**************FIND STAR END****************

*************PROJECT STAR TO PLANE*********
c for common
	do 49 j=1,m512
	do 49 i=1,m695
49	map(i,j)=2				! for show cursor
	xm_max=-3.				! for pgplot display
	xm_min=20.				! star mag range 

	do 1050 i=1,nsao
	j=emsao(i)
	xmag=j*0.1
	if(xmag.lt.-1.65)xmag=-xmag
	nj1=xmag+5
	xm_min=min(xm_min,xmag)
	xm_max=max(xm_max,xmag)
	call g89(xmag,map33,33,w_de,nj1)
	ix=msaox(i)
	iy=msaoy(i)
	do 1053 j=nj1,34-nj1
	jj=iy-17+j-1
	if(jj.lt.1.or.jj.gt.m512)goto 1053
	do 1052 k=nj1,34-nj1
	ii=ix-17+k-1 
	if(ii.lt.1.or.ii.gt.m695)goto 1052
	map(ii,jj)=map(ii,jj)+map33(k,j)
1052	continue
1053	continue
1050	continue
	write(*,*)'Total SAO',nsao
 
	open(19,file='/EOD/stell.nam',status='old')
	do 910 i=1,72
	read(19,'(a3,1x,a6,2i6)')s1(i),c6,rec1(i),rec2(i)
	call tocap(s1(i))
	k=0
	do 910 j=1,6,2
	k=k+1
	d3(k,i)=(b6(j)+96)*100+b6(j+1)+96
910	if(b6(j).eq.32)d3(k,i)=101
	close(19)

	read(11,'(a)')dis                               ! read minor name
	nk=0
	do 199 i=1,72,2
	if(bis(i).eq.32)goto 201
	nk=nk+1
199	than(nk)=(bis(i)+96)*100+bis(i+1)+96
201	read(11,'(a)')dis                               ! read when find
	nkp=0
	do 299 i=1,72,2
	if(bis(i).eq.32)goto 209
	nkp=nkp+1
299	phan(nkp)=(bis(i)+96)*100+bis(i+1)+96
209	read(11,*)sx(1),sy(1),sx(2),sy(2),ira,ide
	read(11,'(a)')tfile
	close(11)

        write(*,*)tfile
	open(11,file=tfile,status='old',iostat=ierr)
        if(ierr.ne.0)then
          write(*,*)tfile,' not found !'
          stop
        endif
        close(11)

        if(s3.ge.3)goto 1000
	if(s3.eq.1)goto 2000                     ! fixed minor
*************** PICTURE BEGIN ********
        icall=pgbegin(0,'/xw',1,1)
	call pgpap(s10152*1.4,s737)
16022	call pgsci(5)
	call pgsch(1.)
	call pgslw(1)
	call pgenv(1.,x695,1.,x512,0,-1)
	call pgxwjpgray(map,m695,m512,1,m695,1,m512,1.,10.,1)
	if(i3.lt.0)call pgxchinese(x695-45.,x512+18.,q(11),3,2,0)      ! north
	x=x695*.5+.5
	y=x512*.5+.5
	call pgsci(2)
	call pgsch(0.7)
	if(i3.lt.0)then
	  call pgrid(0)                              ! net
	  call pgsci(15)
          call pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3,6.)
          call pgsci(7)
	  call plb(alpha,delta)                      ! sun trace
        else
          call pgsci(8)
	  call pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3,5.)  ! star_name
        endif
	do 9912 i=1,nsn
	do 9912 j=1,72
9912	if(sn(i).eq.s1(j))call pstell(rec1(j),rec2(j))   ! stell
        call pgsci(5)
	if(i3.lt.0)call pg4chinese(x695-45.-ira,-20.,q,2,1,0)      ! alpha
	if(i3.lt.0)call pg4chinese(-25.,x512-42.-ide,q(6),2,1,90)  ! delta
        write(*,*)sx(1),sy(1),sx(2),sy(2),ira,ide
        z=x695*0.5-15*nk
        call pg4chinese(z,x512+13.,than,nk,2,0)                ! minor name
	call pgsci(8)
	if(i3.lt.0)call pgxchinese(x695-135.,x512+18.,q5,6,2,0) !  day point
	if(i3.lt.0)call pg4chinese(x695-235.,x512-22.,phan,nkp,1,0) ! find date
	
        if(s3.gt.0)goto 123
        call standc(a_c,d_c,a_c1,d_c1,xi,xn)
        x=hs2*(xhs-xi)+1.5
        y=hs2*(xn-hs)+x512+.5
	call pgsch(2.0)
	call pgpoint(1,x,y,18)
	call pgsch(1.0)
	call pgpoint(1,x,y+13,31)
	call pgsch(0.7)
	call pg4chinese(x-40,y+20,q6,6,1,0)       ! find posi
123     call pgsci(1)
        if(s3.lt.1)call pgsci(11)	
        if(s3.eq.2)d18=18
	call ptrace(tfile,d18,s3,s4,xx,yy)
	if(ch.eq.'g'.or.ch.eq.'G')then
	  call pgend
	  stop
	endif

	write(*,*)'****** On pic, enter "g" or "G"  -- output GIF file'
	call pgcurse(x,y,ch)
	if(ch.eq.'g'.or.ch.eq.'G')then
	  call pgend
	  icall=pgbegin(0,'/jg',1,1)
	  call pgpap(s10152*1.64,s737*1.008)
	  goto 16022
	endif
        stop
c**********************************  movie gif
c wtrace                                   ! s3=+-3,4,5  minor move
1000    jk=0
	open(19,file=tfile,status='old')
1120	read(19,'(a)',end=1121)ns
	read(ns(23:),*)i1,i2,x1,j1,j2,y1
	x_a=i1+i2/60.+x1/3600.	
	if(j1.lt.0)j1=-j1
	x_d=j1+j2/60.+y1/3600.	
	j1=index(ns(23:),'-')
	if(j1.ne.0)x_d=-x_d
	x_a=x_a*cx
	x_d=x_d*cy
	call standc(a_c,d_c,x_a,x_d,xi,xn)
	x=hs2*(xhs-xi)+1.
	y=hs2*(xn-hs)+x512
        if(x.lt.1. .or. x.gt.x695)goto 1120
        if(y.lt.1. .or. y.gt.x512)goto 1120
        jk=jk+1
        xx(jk)=x
	yy(jk)=y
	goto 1120
1121	close(19)
	write(*,*)'points: ',jk
c weak trace end

	open(19,file=tfile,status='old')
        k=0
        kp=0
120     k=k+1
	read(19,'(a)',end=121)ns
	read(ns(23:),*)i1,i2,x1,j1,j2,y1
	x_a=i1+i2/60.+x1/3600.	
	if(j1.lt.0)j1=-j1
	x_d=j1+j2/60.+y1/3600.	
	j1=index(ns(23:),'-')
	if(j1.ne.0)x_d=-x_d
	x_a=x_a*cx
	x_d=x_d*cy
	call standc(a_c,d_c,x_a,x_d,xi,xn)
        ox=zx
	oy=zy
	x=hs2*(xhs-xi)+1.
	y=hs2*(xn-hs)+x512
        zx=x
        zy=y
        if(ox.lt.1. .or. ox.gt.x695)goto 120
        if(oy.lt.1. .or. oy.gt.x512)goto 120
        nz=1
	if(s3.eq.4)nz=2
	if(s3.gt.4)nz=4
	dx=(zx-ox)/nz
        dy=(zy-oy)/nz
        goto 1235	
1234    x=x+dx
	y=y+dy
	nz=nz-1
	if(nz.eq.0)goto 120
	k=k+1
1235    icall=pgbegin(0,'/jg',1,1)
        call pgpap(s10152,s737)
	call pgsch(0.)
	call pgenv(1.,x695,1.,x512,0,-1)
	call pgsch(1.)
	call pgxwjpgray(map,m695,m512,1,m695,1,m512,1.,10.,1)
        if(i3.eq.-1)then
        call pgsci(2)
	call pgrid(0)                              ! net
	call pgsci(15)
	call pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3,6.)  ! star_name
        call pgsci(7)
	call plb(alpha,delta)                      ! sun trace
        else
        call pgsci(8)
	call pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3,4.5)
        endif
	call pgsci(55)
        call pgline(jk,xx,yy)
        call pgsci(3)
	call pgpoint(1,x-1,y-1,18)
	call pgsci(7)
        call pgtext(540.,20.,ns(1:12))
        j=3
	if(s4.ne.0)j=s4
	call pgsci(11)
	call pg4chinese(x,y+5.,than,j,2,0)      ! minor name
        call pgend
        j=index(tfile,".")-1
        write(xc(1:),'(a,i5,a)')tfile(1:j),k,'.gif'
        do i=1,j+5
        if(xc(i:i).eq.' ')xc(i:i)='0'
	enddo
        if(kp.eq.0)then
          write(dis(1:),*)"mv pgplot.gif ",xc
        else
          write(dis(1:),*)"mv pgplot?.gif ",xc
          dis(11:11)=char(mod(kp,10)+48) 
        endif
        write(*,*)dis
        call system(dis)
        kp=kp+1
	goto 1234
121	close(19)
	return
c--------------------------
2000	continue                      ! S3=+-1,  SAO move
c trace_data
        jk=0
	open(19,file=tfile,status='old')
3120	read(19,'(a)',end=3121)ns
	read(ns(23:),*)i1,i2,x1,j1,j2,y1
	x_a=i1+i2/60.+x1/3600.	
	if(j1.lt.0)j1=-j1
	x_d=j1+j2/60.+y1/3600.	
	j1=index(ns(23:),'-')
	if(j1.ne.0)x_d=-x_d
        jk=jk+1
	xt(jk)=x_a*cx
	yt(jk)=x_d*cy
	goto 3120
3121	close(19)

        if(i3.eq.1)then
          icall=pgbegin(0,'/xw',1,1)
	  call pgpap(s10152,s737)
	  call pgsch(1.)
	  call pgenv(1.,x695,1.,x512,0,-1)
	  call pgsci(c1)
          call delay(5)                              ! otherwise, no display
          call pg4chinese(90.,x512+13.,than,nk,2,0)  ! minor name
        endif
	open(19,file=tfile,status='old')
        kp=0
        kk=0
        ij=1
c**************** loop
2120    kk=kk+1                             ! loop here
        read(19,'(a)',end=2121)ns
	read(ns(23:),*)i1,i2,x1,j1,j2,y1
	x_a=i1+i2/60.+x1/3600.	
	if(j1.lt.0)j1=-j1
	x_d=j1+j2/60.+y1/3600.	
	j1=index(ns(23:),'-')
	if(j1.ne.0)x_d=-x_d
	alpha=x_a
	delta=x_d

	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.*s1357
	if(hw_al.gt.12..or.x.lt.0.)hw_al=12.
	up_a=alpha+hw_al
	dn_a=alpha-hw_al	
	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=(x512-1.)/hs*0.5
	xhs=hs*s1357
	
	nsao=0
	call cat_sao(
     *    aasao,ddsao,ensao,emsao,epsao,nnsao,nsao,msaox,msaoy)
c trace
3458    nj=0
        k=0
        j=0
        do 3456 i=ij,jk
	call standc(a_c,d_c,xt(i),yt(i),xi,xn)
	x=hs2*(xhs-xi)+1.
	y=hs2*(xn-hs)+x512
        if(x.lt.1. .or. x.gt.x695)goto 3456
        if(y.lt.1. .or. y.gt.x512)goto 3456
        if(k.ne.0 .and. i-k.ne.1)goto 3457
        k=i
        nj=nj+1
        xx(nj)=x
	yy(nj)=y
        if(j.eq.0)then
          z=(x-x695*.5)**2+(y-x512*.5)**2
          if(z.lt.4.)j=1
        endif
3456    continue
3457    if(j.eq.0)then
          ij=ij+30
          goto 3458
        endif
c trace put to array end
	do 2049 j=1,m512
	do 2049 i=1,m695
2049	map(i,j)=2				! for show cursor
	do 2050 i=1,nsao
	j=emsao(i)
	xmag=j*0.1
	if(xmag.lt.-1.65)xmag=-xmag
	nj1=xmag+5
	call g89(xmag,map33,33,w_de,nj1)
	ix=msaox(i)
	iy=msaoy(i)
	do 2053 j=nj1,34-nj1
	jj=iy-17+j-1
	if(jj.lt.1.or.jj.gt.m512)goto 2053
	do 2052 k=nj1,34-nj1
	ii=ix-17+k-1 
	if(ii.lt.1.or.ii.gt.m695)goto 2052
	map(ii,jj)=map(ii,jj)+map33(k,j)
2052	continue
2053	continue
2050	continue
        if(i3.eq.-1)then
          icall=pgbegin(0,'/jg',1,1) 
          call pgpap(s10152,s737)
	  call pgsch(0.)
	  call pgenv(1.,x695,1.,x512,0,-1)
	  call pgxwjpgray(map,m695,m512,1,m695,1,m512,1.,10.,1)
          call pgsci(c3)
          call pgsch(1.)
	  call pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3,4.)  ! star_name
	  call pgsci(c2)
	  call pgsch(2.)
	  call pgpoint(1,x695*0.5,x512*.5,18)
	  call pgsci(c4)
          call pgsch(1.)
          call pgtext(x695-155.,20.,ns(1:12))
	  call pgsci(55)
          call pgline(nj,xx,yy)
          j=3
	  if(s4.ne.0)j=s4
	  call pgsci(c1)
	  call pg4chinese(x695*0.5,x512*.5+5.,than,j,2,0)      ! minor name
          call pgend
          j=index(tfile,".")-1
          write(xc(1:),'(a,i5,a)')tfile(1:j),kk,'.gif'
          do i=1,j+5
          if(xc(i:i).eq.' ')xc(i:i)='0'
	  enddo
          if(kp.eq.0)then
            write(dis(1:),*)"mv pgplot.gif ",xc
          else
            write(dis(1:),*)"mv pgplot?.gif ",xc
            dis(11:11)=char(mod(kp,10)+48) 
          endif
          write(*,*)dis
          call system(dis)
          kp=kp+1
       else
        do  ii=-2,2
        do  jj=-2,2
        if(ii*ii+jj*jj.ne.8)map(m695/2+ii,jj+m512/2)=99
        enddo
        enddo
	call pgxwjpgray(map,m695,m512,1,m695,1,m512,1.,10.,1)
        call pgsci(c3)
        call pgsch(1.0)
	call pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3,4.)  ! star_name
	call pgsci(c2)
	call pgsch(2.)
	call pgpoint(1,x695*0.5,x512*.5,18)
	call pgsch(1.)
	call pgsci(55)
        call pgline(nj,xx,yy)
        call pgsci(0)
        call pglabel(cc,"","")
        call pgsci(c4)
        call pglabel(ns(1:12),"","")
        cc=ns(1:12)
       endif
        goto 2120
2121	close(19)
        if(i3.eq.1)call pgend
	end

	subroutine pstell(rec1,rec2)
	integer*2 rec1,rec2
	real x(544),y(544)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512
	open(19,file='/EOD/stell.bin',status='old',
     *  access='direct',recl=8)
	j=1
	do 10 i=rec1,rec2
	read(19,rec=i)x(j),y(j)
10	j=j+1
	close(19)
	x(j)=x(1)
	y(j)=y(1)
	do 20 i=1,j
        xa=x(i)*cx
        xd=y(i)*cy
        call standc(a_c,d_c,xa,xd,xi,xn)
        x(i)=hs2*(xhs-xi)+1.5
20      y(i)=hs2*(xn-hs)+x512+.5

	call pgsls(4)
	call pgslw(1)
	call pgsci(3)
	if(j.le.300)call pgline(j,x,y)
	call pgsls(1)
	end	

	subroutine pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3,xp)
	integer ensao(nsao)
	integer*2 msaox(nsao),msaoy(nsao)
	byte emsao(nsao)
	character*8 cv1
	character*3 sn(1),ch3,s1(1)
        integer d3(3,72)
	integer*2 mx(999),my(999),mj(999)
	real xm(999)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512
	jp=0
	do 1011 i=1,nsao
	j=ensao(i)/260000
	if(j.eq.0)goto 1011
	if(msaox(i).lt.22.or.msaox(i).gt.x695-30.)goto 1011
	if(msaoy(i).lt.25)goto 1011
	jp=jp+1
	mj(jp)=j
	mx(jp)=msaox(i)
	my(jp)=msaoy(i)
	j=emsao(i)
	xmag=j*0.1
	if(xmag.lt.-1.65)xmag=-xmag
	xm(jp)=xmag
1011	continue
c	write(*,*)'jp ',jp
	do 2011 i=1,jp
	do 2012 j=1,jp
	if(i.eq.j)goto 2012
	k=(mx(i)-mx(j))**2/2+(my(i)-my(j))**2
	if(k.lt.399)then
	  if(xm(j).gt.xm(i))then
	    mj(j)=0
	  else
	    mj(i)=0
	    goto 2011 
	  endif
	endif
2012	continue 
2011	continue
	nsn=0
	do 1015 i=1,jp
	if(mj(i).eq.0)goto 1015
        if(xm(i).gt.xp)goto 1015
        read(20,rec=mj(i))cv1
	call greek(cv1)
	if(cv1(4:4).eq.' ') cv1(1:3)=cv1(2:4)
	ch3=cv1(6:8)
	call tocap(ch3)
	do 1020 j=1,72
	if(ch3.eq.s1(j))goto 1021
1020	continue
1021	call pgxchinese(mx(i)-5.,my(i)-25.,d3(1,j),3,2,0)
	call pgtext(mx(i)-23.,my(i)-21.,cv1(1:5))
	if(nsn.eq.0)goto 1012
	do 1013 j=1,nsn
1013	if(sn(j).eq.ch3)goto 1015
1012	nsn=nsn+1
c	write(*,*)'nsn',nsn,' ',ch3
	sn(nsn)=ch3
1015	continue
	end

	subroutine g89(xxmag,a,n,ww,nj1)
	integer*2 a(n,n)
	real y(33)
	xmag=xxmag*2.75
	w=ww
	nj2=34-nj1
	do 5 i=1,33		! for guass curve, y(i) > y(i+1)
5	y(i)=0.
	high=16.		! gray scale in unit
	x=600./w+19.-xmag	! for sigma, accoring mag & w_de & sao only
	sigma=0.3+x*.25
	call gauss(y,n,sigma)			! get guass curve
	rr=max((sigma*0.5+0.5)/sqrt(640./w),1.)	! for spark
	do 30 j=nj1,17
	jj=(j-17)**2
	do 30 i=nj1,17
	r=jj+(i-17)*(i-17)
	kr=sqrt(r)+1.5
	r=y(kr)*high+0.5
	if(j.eq.17.or.i.eq.17)r=r*rr		! now spark get value
	a(i,j)=r				! cal 1/4, other give
	a(i,34-j)=r				! value here
	a(34-i,j)=r
30	a(34-i,34-j)=r
	end
	
	subroutine gauss(y,n,sigma)
c cal guass curve
	real y(n)
	do 10 i=1,n
	x=((i-1.)/sigma)**2
	if(x.gt.80.)return
10	y(i)=exp(-.5*x)
	end

	subroutine cat_sao(aa,dd,en,em,ep,nn,n,mx,my)
	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)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512
	equivalence (ee,bb)
	open(1,file='/EOD/sao.dat',status='old',access='direct',recl=512)
	read(1,rec=1)ind
        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
c	 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
        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)
	k=bb(4)
        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 i_star(alpha,delta,ix,iy)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512
        xa=alpha*cx
        xd=delta*cy
        call standc(a_c,d_c,xa,xd,xi,xn)
        ix=hs2*(xhs-xi)+1.5
        iy=hs2*(xn-hs)+x512+.5
	i_star=1
	if(ix.lt.0 .or. ix.gt.x695)i_star=0
	if(iy.lt.0 .or. iy.gt.x512)i_star=0
	end

	subroutine greek(c)
	character*8 c
	character*24 cg
	character*3 dg(24)
	data cg/'abgdezyhiklmncoprstufxqw'/
	data dg/'Alp','Bet','Gam','Del','Eps','Zet','Eta','The',
     *          'Iot','Kap','Lam','Mu ','Nu ','Xi ','Omi','Pi ',
     *          'Rho','Sig','Tau','Ups','Phi','Chi','Psi','Ome'/
	if(c(1:1).lt.'A')return
	do 20 j=1,24
	if(c(1:3).eq.dg(j))goto 30
20	continue
30	if(c(4:4).ne.' ')then
	  c(1:3)=char(92)//'g'//cg(j:j)
	else
	  c(1:4)=' '//char(92)//'g'//cg(j:j)
	endif	
	end

	subroutine plb(alpha,delta)             ! draw sun trace
	real x(50),y(50),z
	integer lbhan(2)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512
	data lbhan/2738,2132/
	if(delta.gt.30.or.delta.lt.-30.)return
	a1=(alpha-1.5)*15.
	a2=(alpha+1.5)*15.
	k=0
	do 10 z=a1,a2
	call lb_ad(z,0.,al,de)
        call standc(a_c,d_c,al*cx,de*cy,xi,xn)
	k=k+1
        x(k)=hs2*(xhs-xi)+1.5
        y(k)=hs2*(xn-hs)+x512+.5
10	continue
        j=0
        do 15 i=1,k
        if(x(i).lt.1. .or. x(i).gt.x695)goto 15
        if(y(i).lt.1. .or. y(i).gt.x512)goto 15
        j=j+1
        if(j.eq.1)jj=i
15      continue
        if(j.lt.3)return
	call pgline(k,x,y)
        z=x(jj)
        if(z.lt.10.)z=10.
        if(z.gt.x695-50.)z=x695-50.
	call pg4chinese(z,y(jj)+5.,lbhan,2,1,0)
	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 pgrid(lw)
	real x9(19),y9(19),z(19)
	character*12 ns
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512
	if(d_c.gt.70*cy.or.d_c.lt.-70*cy)return
	i=nint(a_c/cx*3)
	k=0
	do 20 j=i-6,i+6
	x=j*cx/3.
	k=k+1
	z(k)=x
	y=(d_c/cy+13)*cy
        call standc(a_c,d_c,x,y,xi,xn)
        x9(1)=hs2*(xhs-xi)+1.5
        y9(1)=hs2*(xn-hs)+x512+.5
	y=(d_c/cy-13)*cy
        call standc(a_c,d_c,x,y,xi,xn)
        x9(2)=hs2*(xhs-xi)+1.5
        y9(2)=hs2*(xn-hs)+x512+.5
	call pgline(2,x9,y9)
	call toms(j/3.,ns,1)
        ns=ns(2:3)//'\uh\d'//ns(5:6)//'\um\d'
	if(lw.eq.1)call pgslw(4)
        if(x9(2).gt.0.and.x9(2).le.x695)call pgtext(x9(2)-15.,-19.,ns)
	if(lw.eq.1)call pgslw(2)
20	continue
	
	x=d_c/cy
	if(x.lt.0.)x=-x
	j=3
	if(x.gt.40.)j=2
	if(x.gt.50.)j=1
	do 10 x=-80.,80,5.
	do 11 i=1,k
        call standc(a_c,d_c,z(i),x*cy,xi,xn)
        y9(i)=hs2*(xn-hs)+x512+.5
11      x9(i)=hs2*(xhs-xi)+1.5
	if(y9(1).lt.-29. .or.y9(1).gt.x512+29.)goto 10
	call pgline(k,x9,y9)	                ! dec line
	write(ns(1:),'(i3)')nint(x)
        ns=ns(1:3)//'\uo'
	if(lw.eq.1)call pgslw(4)
        call pgtext(-33.,y9(j),ns)
	if(lw.eq.1)call pgslw(2)
10	continue
	end

	subroutine tocap(s)
	character*3 s
	do 10 i=1,3
	j=ichar(s(i:i))
	if(j.gt.96)j=j-32
10	s(i:i)=char(j)
	end

	subroutine ptrace(tfile,d18,s3,s4,xx,yy)
	real xx(1),yy(1)
	character*12 tfile,ns*55,cc*35,ns1
        integer s3,s4
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  w_de,xlimag,hs2,xhs,a_c1,d_c1,s373,s10152,x695,x512
        k=0
	open(19,file=tfile,status='old')
120	read(19,'(a)',end=121)ns
	read(ns(23:),*)i1,i2,x1,j1,j2,y1
	x_a=i1+i2/60.+x1/3600.	
	if(j1.lt.0)j1=-j1
	x_d=j1+j2/60.+y1/3600.	
	j1=index(ns(23:),'-')
	if(j1.ne.0)x_d=-x_d
	x_a=x_a*cx
	x_d=x_d*cy
	call standc(a_c,d_c,x_a,x_d,xi,xn)
	x=hs2*(xhs-xi)+1.
	y=hs2*(xn-hs)+x512
        if(x.lt.1. .or. x.gt.x695)goto 120
        if(y.lt.1. .or. y.gt.x512)goto 120
	if(s3.ne.1)then
	  call pgpoint(1,x-1,y-1,d18)
        else
          k=k+1
          xx(k)=x
	  yy(k)=y
	endif 
        if(s3.eq.2)then
          call pgsci(0)
          call pglabel(ns1,"","")
          call pgsci(11)
          call pglabel(ns(1:12),"","")
          ns1=ns(1:12)
          if(s4.le.0)s4=5
          call delay(s4)
        endif
	goto 120
121	close(19)
	if(s3.eq.1)then
	  call pgsci(55)
          call pgline(k,xx,yy,-1)
	  write(*,*)'points: ',k
	endif
	end

        subroutine delay(s4)
        integer s4
        real x
        if(s4.le.0)return
        do 20 j=1,s4
        k=0
        do 10 i=1,100000
        x=i-99
        k=k+amod(x,123.)
10      continue
20      continue
        end

	subroutine pgxwjpgray(map,n1,n2,nx1,nx2,ny1,ny2,white,black)
	integer*2 map(n1,n2)
	character dummy
	character*16 c
	real rbuf(1900)
	call pgqinf('TYPE',c,k)
	n16=84
	i1=1
	call pgqwin(x1,x2,y1,y2)
	call grtxy0(.false.,x1,y1,x1,y1)
	ibx1=nint(x1)
	iby1=nint(y1)
	call grtxy0(.false.,x2,y2,x2,y2)
	ibx2=nint(x2)
	iby2=nint(y2)
	sx=(nx2-nx1)/(x2-x1)
	sy=(ny2-ny1)/(y2-y1)
	d=n16/(black-white)
	rbuf(1)=ibx1
	do 20 j=iby1,iby2
	dy=sy*(j-iby1)
	j1=dy				! cannot be j1=dy+ny1
	j1=j1+ny1
	j2=j1+1
	if(j2.gt.n2)j2=n2
	y2=mod(dy,1.)
	y1=1.-y2	
	rbuf(2)=j
	i6=2
	do 10 i=ibx1,ibx2
	dx=sx*(i-ibx1)
	i1=dx
	i1=i1+nx1
	i2=i1+1
	if(i2.gt.n1)i2=n1
	x2=mod(dx,1.)
	x1=1.-x2
	z=map(i1,j1)*x1*y1+map(i1,j2)*x1*y2+
     *    map(i2,j1)*x2*y1+map(i2,j2)*x2*y2
	k=nint(d*(z-white))
	if(k.lt.1)k=1
	if(k.gt.n16)k=n16
	i6=i6+1
10	rbuf(i6)=k+15
	if(c(1:2).eq.'XW')then
	call xwdriv(26,rbuf,i6,dummy)
	call pgebuf
	endif
20	if(c(1:2).eq.'JG')call jgdriv(26,rbuf,i6,dummy)
	end
