        parameter (nnsao=1999)			! for max number of sao
	character cc*25,dis*73			! read & display line	
	byte bis(73)
	character*1 ch				! ch~~ job indicator
	character ch3*3				! pgplot title,scale,scirpt 
	character*12 ac,dc,tfile		! normal format of alpha delta
	character*17 xc				! line for deal with default
	character*1 c5*5,cv1*8	

	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

	real sx(5),sy(5) 			! scale & cross line at pic.

	integer*2 map(512/8*11,512),map33(33,33)	! 512*512 map & star map
	integer q(13)				! chinese character array

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

	character*3 s1(72),c6*6,sn(19)
	integer d3(3,72),dfx(5),dfy(5)
	byte b6(6)
	character*30 dfile(3)
	integer*2 rec1(72),rec2(72)

	integer than(30),phan(30),q6(6),q4(4),q5(6)

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

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

c	data zx,zy/650,650,0,0,0,15,10,10,15,0/
	data q/1964,3013,308,3647,309,1964,4619,308,2211,309,308,1717,309/
	data q6/2302,4754,4217,2136,4627,5435/
	data q5/0308,3531,4476,5027,2167,0309/
	data q4/5243,4848,2329,4782/

	if(iargc().lt.1)then
	  write(*,1001)
1001	  format(//' *********** FINDING CHART used SAO
     *    *********'/60x,'2000.7 Jiang'/)
	  write(*,*)
	  write(*,*)'     Usage: sao4 file'
          stop
	endif
	call getarg(1,xc)  
	open(11,file=xc,status='old')
	
	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.
	bell=7
	ep_old=2000.0				! always be 1950.0
	w_de=900.
	xlimag=7.
        icall=0
           
	x512=512/8*11.
	n512=nint(x512)
	s512=11./8
	epoch=2000.

	call unlink('sao.tmp')	
	read(11,'(a)')dis
	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
	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
c picture is fixed 512*512 pixels
	sx(3)=650.-i*60/w_de*512.	
	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.*s512
	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
	xhs=hs*s512
	
**************************FIND STAR***********************
1211	nsao=0
	call cat_sao(
     *    aasao,ddsao,ensao,emsao,epsao,nnsao,nsao,msaox,msaoy)
**************FIND STAR END****************

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

ccc for sao
	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.512)goto 1053
	do 1052 k=nj1,34-nj1
	ii=ix-17+k-1 
	if(ii.lt.1.or.ii.gt.n512)goto 1052
	map(ii,jj)=map(ii,jj)+map33(k,j)
1052	continue
1053	continue
1050	continue
	write(*,*)'Total SAO',nsao
 
	if(nsao.eq.0)then
	write(*,'(1x,2a)')bell,'Sorroy, No Star in this region !'
	stop
	endif
 
*************** PICTURE BEGIN ********
	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)

        icall=pgbegin(0,'/xw',1,1)
	call pgpap(13.75,10./13.75)
        if(icall.eq.1)goto 16011
16011   open(99,file='sao.tmp',status='unknown')
16022	call pgsci(5)
	call pgsch(1.)
	call pgslw(1)
	call pgenv(1.,x512,1.,512.,0,-1)
	call pgxgray(map,n512,512,1,n512,1,512,1.,10.,1)
	call pgxchinese(650.,530.,q(11),3,4,0)
	call pgsci(2)
	num=0
	call pgsch(0.7)
	x=256.5*s512
	y=256.5
	inpxy=0
	inpn=0

	call pgrid(0)
	call plb(alpha,delta)
	call pgsci(15)
	call pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3)
	do 9912 i=1,nsn
	do 9912 j=1,72
9912	if(sn(i).eq.s1(j))call pstell(rec1(j),rec2(j))
	if(ch.eq.'g'.or.ch.eq.'G')goto 2009
	read(11,'(a)')dis
	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
	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
2009    call pgsci(5)
	call pg4chinese(660.-ira,-20.,q,2,2,0)
	call pg4chinese(-25.,470.-ide,q(6),2,2,90)
            write(*,*)sx(1),sy(1),sx(2),sy(2)
	call pg4chinese(50.,525.,than,nk,4,0)
	call pgsci(8)
        call pg4chinese(380.,530.,q5,6,2,0)
	call pg4chinese(400.,485.,phan,nkp,2,0)
        call standc(a_c,d_c,a_c1,d_c1,xi,xn)
        x=hs2*(xhs-xi)+1.5
        y=hs2*(xn-hs)+512.5
	call pgsch(2.0)
	call pgpoint(1,x,y-1,18)
	call pgsch(1.0)
	call pgpoint(1,x,y+13,31)
	call pgsch(0.7)
	call pg4chinese(x-40,y+20,q6,6,2,0)

	x=atan2(sy(2)-sy(1),sx(2)-sx(1))+0.785
	sx(3)=sx(2)-6*cos(x)
	sy(3)=sy(2)-6*sin(x)
	call pgline(3,sx,sy)
	x=atan2(sy(2)-sy(1),sx(2)-sx(1))-0.785
	sx(3)=sx(2)-6*cos(x)
	sy(3)=sy(2)-6*sin(x)
	call pgline(2,sx(2),sy(2))
	call pgxchinese(sx(1),sy(1)-15,q4,4,2,0)
        call pgsci(11)	
	if(ch.eq.'g'.or.ch.eq.'G')then
	  call ptrace(tfile)
	  call display(dfile,dfx,dfy,1)
	  call pgend
	  stop
	endif

	read(11,'(a)')tfile
	call ptrace(tfile)
	do 202 i=1,3
	read(11,'(a)')dis
	k=index(dis,' ')
	if(k.eq.0)then
          write(*,*)'no cursor position'
	  stop
	endif
        read(dis(k:),*)dfx(i),dfy(i)
	dfile(i)=dis(1:k)
202	write(*,'(a,2i3)')dfile(i),dfx(i),dfy(i)
	close(11)
	call display(dfile,dfx,dfy,-1)

********************************************* loop
100	write(*,*)'  left button -- pick out a star'
	write(*,*)'middle button -- put a false_star'
	write(*,*)' right button -- quit & next turn'
	write(*,*)'  "p" or "P"  -- print finding chart'
	write(*,*)'  "g" or "G"  -- output GIF file'
	
110	call pgcurse(x,y,ch)
	ii=ichar(ch)

c GIF*****************
	if(ch.eq.'g'.or.ch.eq.'G')then
	call pgend

	icall=pgbegin(0,'/jg',1,1)
	call pgpap(13.75,10./13.75)
	goto 16022
	endif

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

c middle  button**************
	if(ch.eq.'D')then
	call pgsci(4)
	x=x+1
	y=y+1
	call pgpoint(1,x-1,y-1,16)
	xi=x
	xn=y
c accordint cursor position to cal alpha delta
	xi=xhs-(xi-1.)/hs2
	xn=hs+(xn-512.)/hs2
	call astand(a_c,d_c,xi,xn,xa,xd)
	x_a=xa/cx
	x_d=xd/cy
c because of frame in 2000_epoch, so change to chart_epoch for display
	call astprs(x_a,x_d,2000.,xa,xd,epoch)
	xa=amod(xa+24.,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgtext(x-3.,y+4.,ch3)
	write(*,104)num,ac,dc,epoch
104	format(i3,') ',a,1x,a,'(',f6.1,') middle button')
	call pgsci(0)
	call pgtext(0.,-35.,dis)
	write(dis(1:),104)num,ac,dc,epoch
	call pgsci(9)
	call pgtext(0.,-35.,dis)
	write(99,'(a)')dis
	endif

c left button****************
	if(ch.eq.'A')then
        ip=0
	call pgsci(6)
c pick out the nestest star of cursor

c search SAO (3)
3010	if(nsao.eq.0)goto 1012
	r=1e10
	do 1011 i=1,nsao
	rr=(msaox(i)-x)**2+(msaoy(i)-y)**2
	if(rr.lt.r)then
	m=i
	r=rr
	endif
1011	continue
	if(r.gt.32.)goto 1012		! if fail, search gsc next.
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgtext(msaox(m)-3.,msaoy(m)+4.,ch3)
	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
	cv1=' '
        if(j.ne.0)read(20,rec=j)cv1
	
	k=mod(ensao(m),260000)	
	
	call pgsci(0)
	call pgtext(0.,-35.,dis)
	write(*,1333)num,ac,dc,epoch,k,c5,ch3,cv1
1333	format(i3,') ',a,1x,a,'(',f6.1,')      SAO',i7,3(' ',a))
	call greek(cv1)
	write(dis(1:),1333)num,ac,dc,epoch,k,c5,ch3,cv1
	call pgsci(9)
	call pgtext(0.,-35.,dis)
	write(99,'(a)')dis
	goto 110
	endif

1012    write(*,'(1x,a)')bell
	goto 110	
	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,
     *  epoch,w_de,xlimag,hs2,xhs,a_c1,d_c1
	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)+512.5

	call pgslw(1)
	call pgsci(3)
	call pgline(j,x,y)
	end	

	subroutine pname(ensao,msaox,msaoy,emsao,nsao,sn,nsn,s1,d3)
	integer ensao(nsao)
	integer*2 msaox(nsao),msaoy(nsao)
	byte emsao(nsao)
	character*8 cv1
	character*3 sn(19),ch3,s1(72)
        integer d3(3,72)

	integer*2 mx(99),my(99),mj(99)
	real xm(99)
	jp=0

	do 1011 i=1,nsao
	j=ensao(i)/260000
	if(j.eq.0)goto 1011
	if(msaox(i).lt.20.or.msaox(i).gt.680)goto 1011
	if(msaoy(i).lt.20)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
	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
        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    k=3
        if(d3(3,j).eq.101)k=2
        call pgxchinese(mx(i)-5.,my(i)-20.,d3(1,j),k,2,0)
	call pgtext(mx(i)-19.,my(i)-18.,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
	write(*,*)'nsn',nsn,' ',ch3
	sn(nsn)=ch3
1015	continue

	end

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

	subroutine pgxcross(x,y,h)
c draw cross in Hi_pre pic
c h is center hollow size
	real xx(2),yy(2)
	xx(1)=x-3.
	xx(2)=x-h
	yy(1)=y
	yy(2)=y
	call pgline(2,xx,yy)
	xx(1)=x+h
	xx(2)=x+3.
	call pgline(2,xx,yy)
	xx(1)=x
	xx(2)=x
	yy(1)=y-3.
	yy(2)=y-h
	call pgline(2,xx,yy)
	yy(1)=y+h
	yy(2)=y+3.
	call pgline(2,xx,yy)
	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+20.-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,
     *  epoch,w_de,xlimag,hs2,xhs,a_c1,d_c1
	equivalence (ee,bb)
	dy=2000.-epoch
	open(1,file='/EOD/sao.dat',status='old',access='direct',recl=512)
	read(1,rec=1)ind
c	 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
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)
        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 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,xhs,a_c1,d_c1
        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)+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	
	
	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 peq
	real x(2),y(2)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2,xhs,a_c1,d_c1
        call standc(a_c,d_c,a_c,0.,xi,xn)
        y(1)=hs2*(xn-hs)+512.5
	y(2)=y(1)
	x(1)=0.
	x(2)=512.*11/8.
	call pgline(2,x,y)	
	end

	subroutine plb(alpha,delta)             ! draw sun trace
	real x(30),y(30)
	integer lbhan(2)
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2,xhs,a_c1,d_c1
	data lbhan/2738,2132/
	if(delta.gt.30.or.delta.lt.-30.)return
	a1=(alpha-0.9)*15.
	a2=(alpha+0.9)*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)+512.5
10	continue
	call pgsci(7)
	call pgline(k,x,y)
	do 20 i=1,30
	if((x(i).gt.0.and.x(i).lt.100.)
     *  .and.(y(i).gt.0.and.y(i).lt.512))goto 30
20	continue
	return
30	call pg4chinese(x(i),y(i),lbhan,2,2,0)
	write(*,*)i,x(i),y(i)
	call pgsci(8)
	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(9),y9(9),z(9)
	character*12 ns
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2,xhs,a_c1,d_c1
	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-4,i+4
	x=j*cx/3.
	k=k+1
	z(k)=x
	y=(d_c/cy+8)*cy
        call standc(a_c,d_c,x,y,xi,xn)
        x9(1)=hs2*(xhs-xi)+1.5
        y9(1)=hs2*(xn-hs)+512.5
	y=(d_c/cy-8)*cy
        call standc(a_c,d_c,x,y,xi,xn)
        x9(2)=hs2*(xhs-xi)+1.5
        y9(2)=hs2*(xn-hs)+512.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.704)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)+512.5
11      x9(i)=hs2*(xhs-xi)+1.5
	if(y9(1).lt.1..or.y9(1).gt.512.)goto 10
	call pgline(k,x9,y9)	
	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)
	character*12 tfile,ns*35
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2,xhs,a_c1,d_c1
	open(19,file=tfile,status='old')
120	read(19,'(23x,a)',end=121)ns
	read(ns(1:),*)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,'-')
	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)+512.
	call pgpoint(1,x-1,y-1,18)
	goto 120
121	close(19)
	end

	subroutine display(nf4,dfx,dfy,ip)
        character*30 nf4(3)
	integer dfx(5),dfy(5)
	character*8 dft(5)
	character*12 ac,dc
	real*8 a8(8),adcoef(2,3)
	real sx(5),sy(5), pho(5),xw(5)
	character*80 head(72)	
        common /fid/a_c,d_c,up_d,dn_d,up_a,dn_a,iaf,hs,cx,cy,
     *  epoch,w_de,xlimag,hs2,xhs,a_c1,d_c1

	do 10 m4=1,3
        if(m4.eq.1)then
          call gethead(nf4(1),head)               ! mother file
          call geta8(head,a8)
          call xytoad(a8,adcoef)
          x=a8(7)
          y=a8(8)
          sx(1)=0
          sy(1)=0
	  xx=x/cx
	  yy=y/cy
	  call toms(xx,ac,1)
	  call toms(yy,dc,0)
	  write(*,*)' Center  ',ac,dc,'(2000.0)'
          ipos=indexpos(head,'VOLT3   ')
	  read(head(ipos)(25:),*)pho(m4)
        else 
          call gethead(nf4(m4),head)            
          ipos=indexpos(head,'A87     ')
          read(head(ipos)(20:),*)alpha
          read(head(ipos+1)(20:),*)delta
          ipos=indexpos(head,'VOLT3   ')
	  read(head(ipos)(25:),*)pho(m4)
          call standc(x,y,alpha,delta,xi,xn)
          sx(m4)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)-1024.5
          sy(m4)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)-1024.5
        endif
        ipos=indexpos(head,'TIME    ')
	dft(m4)=head(ipos)(20:27)
10	continue
	m4=3
        do 7 i=1,m4
7       write(*,8)i,nf4(i)(1:25),sx(i),sy(i),pho(i),dft(i)
8       format(i2,':  ',a,2f6.2,f9.2,2x,a)
        alpha=a_c1
	delta=d_c1
          call standc(x,y,alpha,delta,xi,xn)
          alpha=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
          delta=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	  alpha=2049.-alpha
	  delta=2049.-delta
        j1=600
	do 13 i=1,m4
        i1=alpha+sx(i)+0.5
        i2=delta+sy(i)+0.5
        pho(i)=disp4(nf4(i),i1,i2,j1,630,xw(i),dfx(i),dfy(i),ip)
	call pgsch(0.6)
	call pgtext(float(j1-50)*0.785,382.,dft(i)//' (UT)')
13      j1=j1+120
	end	

	function disp4(f1,i1,i2,ix,iy,xw,dfx,dfy,jp)
	character*30 f1
	integer*2 map(101*101),ma(101*101)
	integer dfx,dfy
	real xmap(101*101)
	call readcc(f1,xmap,i1,i2)
	ip=0
        do 203 i=1,101*101
	if(xmap(i).le.0.)then
	  map(i)=32767
	else
   	  ip=ip+1
	  j=xmap(i)*.25
	  if(j.gt.32766)j=32766
	  ma(ip)=j
          map(i)=j
	endif
203	continue
	call whitexblack(ma,1,ip,white,sigma)
	  write(*,*)white,sigma
	if(white.gt.0.)then
	  xw=white
	else
	  white=xw
	endif
	dfy=dfy*jp
	i=dfx+51
	do 204 j=dfy+55,dfy+63
204	map(j*101+i)=32767
	i=dfx+50
	j=dfy+56
	map(j*101+i)=32767
	map(j*101+i+2)=32767
	j=dfy+57
	map(j*101+i)=32767
	map(j*101+i+2)=32767
	i=dfx+49
	map(j*101+i)=32767
	map(j*101+i+4)=32767
	iix=ix
	iiy=iy
	if(jp.eq.1)iiy=iiy+5
        call ximage2(map,101,iix,iiy,white*0.95,white+20.*sigma)
	end

        subroutine ximage2(map,n,iix,iiy,black,white)
        integer*2 map(n,n)
        character dummy,c*16
        real rbuf(101*101+2)
        call pgqinf('TYPE',c,k)
        if(c(1:2).eq.'XW'.or.c(1:2).eq.'JG')then
        n16=84
        d=n16/(black-white)
        rbuf(1)=iix
	do 40 j=1,n
        rbuf(2)=iiy-j
        i6=2
        do 30 i=1,n
        k=nint(d*(map(i,j)-white))
        if(k.lt.1)k=1
        if(k.gt.n16)k=n16
        i6=i6+1
30      rbuf(i6)=k+15
	if(c(1:2).eq.'XW')then
        call xwdriv(26,rbuf,i6,dummy)
        call pgebuf
	else
        call jgdriv(26,rbuf,i6,dummy)
	endif	
40	continue
	endif

        if(c(1:2).eq.'LH')then
	i=n
	call turn2(map,n)
        ix1=iix*4
	  ix1=ix1-245
	  iy1=1612

        ix2=ix1+(n-1)*4
        iy2=iy1+(n-1)*4
        sx=.25
        sy=.25
        ip=2
        scale=8.
        ip2=ip+ip
        ip05=ip/2
                ix11=ix1+ip
                ix22=ix2-ip05
                ixstep=ip

                iy11=iy1+ip
                iy22=iy2-ip05
                iystep=ip2
                iyd=ip
        ratio=scale/(white-black)
        k=-1
        do 50 ix=ix11,ix22,ixstep
        k=-k
        kd=0
        if(k.eq.1)kd=iyd
        dx=sx*(ix-ix1)
        i1=dx
        i1=i1+1
        i2=i1+1
        x2=mod(dx,1.)
        x1=1.-x2
        rbuf(1)=ix
        do 50 iy=iy1+kd,iy22,iystep
        rbuf(2)=iy
        dy=sy*(iy-iy1)
        j1=dy
        j1=j1+1
        j2=j1+1
        y2=mod(dy,1.)
        y1=1.-y2
        z=map(i1,j1)*x1*y1+map(i1,j2)*x1*y2+
     *    map(i2,j1)*x2*y1+map(i2,j2)*x2*y2
        rbuf(3)=ratio*(z-black)
        rbuf(4)=scale
        call lhdriv(28,rbuf,6,cdummy,1)
50      continue
	endif
        end

	subroutine turn2(map,n)
	integer*2 map(n,n)
	do 10 j=1,n/2
	do 10 i=1,n
	k=map(i,j)
	map(i,j)=map(i,n-j+1)
10	map(i,n-j+1)=k
	end

        subroutine gethead(f1,head)
        character*80 head(72),f1*30
        open(52,file=f1,status='old',access='direct',recl=5760)
        read(52,rec=1)head
        close (52)
        end

        subroutine geta8(head,a8)
        real*8 a8(8)
        character*80 head(72)
        ipos=indexpos(head,'A81     ')-1
        do 10 i=1,8
10      read(head(i+ipos)(12:),*)a8(i)
        end

        function indexpos(head,f1)
        character*80 head(72),f1*8
        do 10 indexpos=1,72
10      if(head(indexpos)(1:8).eq.f1)return
        end
