c simple sao
	parameter (nnsao=1999)			! for max number of sao
	character cc*25,ns*35,dis*73,ns1*32	! read & display line	
	character*1 ch				! ch~~ job indicator
	character*1 sg                    ! se,sm,sg ~~for display line
	character label*50,sc*4,ch3*3		! pgplot title,scale,scirpt 
	character*12 ac,dc,ac1,dc1		! normal format of alpha delta
	character*17 xc				! line for deal with default
	character*1 c5*5	

	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(4),sy(4),crox(4),croy(4)	! scale & cross line at pic.
	real xya(99),xyd(99)			! for XY alpha & delta
	real*8 xixn(2,99),xy(2,99),xycoef(2,6)	! for XY ,plate 6 coef

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

	byte w910(2),bell			! (last 2 byte)
	byte xyw(99),ixy(99)			! for XY weight & num

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

	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)

	data sx,sy/500,500,0,0,15,10,10,15/
	data crox,croy/253,260,253,260,253,260,260,253/
	data q/1964,3013,308,3647,309,1964,4619,308,2211,309,308,1717,309/

	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=2000.0				! always be 1950.0
	w_de=240.
	xlimag=7.
        icall=0
           
	epoch=2000.
	ns1='   '
	call unlink('sao.tmp')			! sao.tmp hold line message
	write(*,1001)
1001	format(//' *********** FINDING CHART used SAO
     *  *********'/60x,'1996.12 Jiang'/)

10	if(ns1(1:2).eq.'  ')write(*,9001)ep_old
9001	format('Input: alpha,delta,epoch
     *  (Ex.3:46:02.11 23:35:47.8 [',f6.1,'])  "q"~quit **gs**')
	if(ns1(1:2).ne.'  ')write(*,9009)ns1
9009	format('Input: alpha,delta,epoch [',a,']  "q"~quit **gs**')
	read(*,'(a)')dis
	if(dis(1:2).eq.'  ')dis=ns1(2:)
	if(dis(1:1).eq.'q'.or.dis(1:1).eq.'Q')stop
	if(itohd(dis,alpha,delta,x).eq.0)goto 10	! char_line  to 3 real
	if(delta.gt. 89.5)delta= 89.5		! near 90 is too large cal.err
	if(delta.lt.-89.5)delta=-89.5
	if(x.ne.0.)ep_old=x
	goto 112				! if char_line include 3 item
11	write(*,*)'epoch (ex. 1950.0)'
	read(*,*,err=11)ep_old
112	if(ep_old.lt.1900..or.ep_old.gt.2100.)goto 11
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	lsky=.true.
	lns1=.true.
	if(delta.lt.0.)lsky=.false.		! in north sky, lsky=.true.
	write(ns1(1:),'(a,1x,a,f7.1)')ac,dc,ep_old

111	write(*,9004)epoch
9004	format('finding chart epoch [',f6.1,']')
	read(*,'(a)',err=111)xc
	if(xc(1:2).ne.'  ')read(xc(1:),*,err=111)epoch
	if(epoch.lt.1900..or.epoch.gt.2100.)goto 111
c first change epoch to finding chart_epoch, use to display
c secondly, change to 2000_epoch, because of gsc is 2000_epoch
c so, in any where of this program, when calculate: use 2000_epoch
c				    when display: use chart_epoch 
c sao is 2000 epoch
c ngc is the same epoch as sao
c alpha delta in main program: in hour & degree express
	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

*********** "Center" loop
	dis='   '
12	write(*,9003)xlimag
9003	format('magnitude limited, ( -value: change color) [',f4.1,'] ')
	read(*,'(a)',err=12)xc
	if(xc(1:2).ne.'  ')then
	read(xc(1:),*,err=12)xlimag
	if(xlimag.lt.0.)then
	xlimag=-xlimag
	lns1=.false.
	endif
	endif
ccccccc	lns=lns1.xor.lsky
	w_de=w_de/60.
	if(xlimag.lt.5.)goto 12
1000	write(*,9002)nint(w_de)
9002	format('Finding chart size, in degree.  [',i3,' ]')
	read(*,'(a)',err=12)xc
	if(xc(1:2).ne.'  ')read(xc(1:),*,err=1000)w_de
	w_de=w_de*60
	if(w_de.gt.750.)w_de=750.
	write(*,*)'      WAIT !'
	i=w_de*0.25/60.+0.5
	write(sc(1:),'(i2)')i			! display scale on pic.
	sc(4:4)='d'
c picture is fixed 512*512 pixels
	sx(3)=500.-i*60/w_de*512.			! how long is scale line
	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.
	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***********************
1211	nsao=0
	call cat_sao(file,
     *    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,512
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)
        xmag=xmag*1.8
	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.512)goto 1052
	map(ii,jj)=map(ii,jj)+map33(k,j)
1052	continue
1053	continue
1050	continue
	write(*,*)'Total SAO',nsao

	if(n+nsao.eq.0)then
	write(*,'(1x,2a)')bell,'Sorroy, No Star in this region !'
	goto 10
	endif
 
*************** PICTURE BEGIN ********
        icall=pgbegin(0,'/xw',1,1)
	call pgpap(8.5,1.)
        if(icall.eq.1)goto 16011
16011   open(99,file='sao.tmp',status='unknown')
	call pgsci(1)
	call pgsch(1.)
	call pgslw(1)
	call pgenv(1.,512.,1.,512.,1,-1)
	call pgxgray(map,512,512,1,512,1,512,1.,10.,1)
	write(label(1:),"(' (',f6.1,')')")epoch
	label=cc(1:lnblnk(cc))//label(1:9)
	il=lnblnk(label)+1
	write(label(il:),"(f6.1,'~~',f4.1)")xm_min,xm_max
	call pglabel('R.A. South','Decl, East',label(1:il+11)//'\um\d')
	call pgxchinese(400.,-50.,q,5,0,0)
	call pgxchinese(-45.,400.,q(6),5,0,90)
	call pgxchinese(450.,540.,q(11),3,0,0)
	call pgsci(2)
	call pgline(4,sx,sy)
	call pgtext(450.,15.,sc)
	call pgline(2,crox,croy)
	call pgline(2,crox(3),croy(3))
	num=0
	call pgsch(0.7)
	x=256.5
	y=256.5
	inpxy=0
	inpn=0
	do 99 i=1,99
99	xyw(i)=0
********************************************* loop
100	write(*,*)'  left button -- pick out a star'
	write(*,*)'middle button -- put a false_star'
	write(*,*)' right button -- quit & next turn'
	write(*,*)'  "n" or "N"  -- input a false_star  (in alpha delta)'
	if(inpxy.eq.1)write(*,*)'  "m" or "M"  -- input a false_star
     *  (in X & Y)'
	write(*,*)'  "p" or "P"  -- print finding chart'
	write(*,*)'  "c" or "C"  -- new center'
	write(*,*)'  "x"         -- input X,Y coordinates'
	write(*,*)'  "h" or "H"  -- help'
	write(*,*)'  "s" or "S"  -- switch'
	
110	call pgcurse(x,y,ch)
	ii=ichar(ch)
	if(ch.eq.'h'.or.ch.eq.'H')goto 100
	if(ch.eq.'s'.or.ch.eq.'S')then
	call pgtext(518.,270.,' OFF')
	call pgsci(1)
	goto 110
	endif

c input X,Y as a false_star*******
	if((ch.eq.'m'.or.ch.eq.'M').and.inpxy.eq.1)then
170	write(*,*)'Input X , Y  (ex. 12.234, -123.230)'
	read(*,*,err=170)x,y
	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(a_c,d_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
	xa=amod(xa+24.,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	call astprs(xa,xd,epoch,x_a,x_d,2000.)
	xa=x_a*cx
	xd=x_d*cy
	call standc(a_c,d_c,xa,xd,xi,xn)
	x=hs2*(hs-xi)+1.
	y=hs2*(xn-hs)+512.
	call pgsci(13)
	call pgpoint(1,x-1,y-1,18)
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgtext(x-3.,y+4.,ch3)
	write(*,171)num,ac,dc,epoch
171	format(i3,') ',a,1x,a,'(',f6.1,') input X Y')
	call pgsci(0)
	call pgtext(0.,-20.,dis)
	write(dis(1:),171)num,ac,dc,epoch
	call pgsci(9)
	call pgtext(0.,-20.,dis)
	write(99,'(a)')dis
	write(*,"(' Again ? [Y] ',$)")
	read(*,'(a)')sg
	if(sg.ne.'n'.and.sg.ne.'N')goto 170	
	endif

c input X,Y coordinates********
	if(ch.eq.'x')then
	if(inpxy.eq.0)then
	write(*,*)'Notice: use left_button first, pick out 4 stars at least'
	else
	write(*,*)'You can modify your data which inputed last time OR add XY'
	endif
	if(inpn.lt.3)goto 190
	write(*,*)'if you want to quit this sub_task, turn to next stage,
     *  then input: -1'
	write(*,*)'input X,Y example: number,X_posi,Yposi ( 6,12.125,-23.230)'
	do 180 i=1,inpn
	if(xyw(i).ne.0)then
	sg=' '
	if(xyw(i).eq.1)sg='*'
	write(*,'(i3,2f8.3,3x,a)')ixy(i),xy(1,i),xy(2,i),sg
	endif
180	continue
181	write(*,"(1x,' * ',$)")
	read(*,'(a)',err=181)xc
	if(xc(1:2).eq.'-1')goto 185
	read(xc(1:),*,err=181)k,x,y
	do 182 i=1,inpn
	if(k.eq.ixy(i))goto 183
182	continue
	write(*,'(2a)')bell,
     *  'sorry, cannot find star picked out by left button.'
	goto 181
c read keyborad finished
183	xy(1,i)=x
	xy(2,i)=y
	xa=xya(i)*cx
	xd=xyd(i)*cy
	call standc(a_c,d_c,xa,xd,xi,xn)
	xixn(1,i)=xi
	xixn(2,i)=xn
	xyw(i)=1		! weight, adjust it in next stage
	goto 181
185	k=0
	do 186 i=1,inpn
186	if(xyw(i).eq.1)k=k+1
	if(k.lt.4)then
	write(*,*)k,' X,Y not enough, at least 4'
	goto 190
	endif
c get plate 6 coef.
c here, plate center is fixed already, so it is not an avariable.
	call plate(xixn,xy,xyw,inpn,xycoef,3) 
	write(*,*)(xycoef(1,i),i=1,3)
	write(*,*)(xycoef(2,i),i=1,3)
	sig=0.
	m=0
	do 187 i=1,inpn
	if(xyw(i).eq.0)goto 187
	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(a_c,d_c,xi,xn,xa,xd)
	xa=xa/cx
	xd=xd/cy
	sigma=((xya(i)-xa)*15.*cos(d_c))**2+(xyd(i)-xd)**2
	sigma=sqrt(sigma)*3600.
	if(xyw(i).eq.1)then
	sig=sig+sigma
	m=m+1
	endif
	call toms(xya(i),ac1,1)
	call toms(xyd(i),dc1,0)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	if(i.eq.1)write(*,184)epoch,epoch
184	format(12x,'old (',f6.1,')',14x,'new (',f6.1,')'/2(2x,24('-')))
	sg=' '
	if(xyw(i).eq.1)sg='*'
187	write(*,'(i3,2a,2x,2a,g12.2,,4x,a)')ixy(i),ac1,dc1,ac,dc,sigma,sg
	write(*,'(62x,g12.2)')sig/m
	inpxy=1				! indicate: 6 coef existed now !
188	write(*,*)'input number: re_calculate plate constants, -1 for quit'
	read(*,*)k
	if(k.eq.-1)goto 190
	do 189 i=1,inpn
	if(k.ne.ixy(i))goto 189
	if(xyw(i).eq.0)goto 188
	xyw(i)=-xyw(i)			! delete or recover a star
	goto 185
189	continue
	goto 188
190	write(*,'(a)')bell
	goto 110
	endif	

c new center***************
	if(ch.eq.'c'.or.ch.eq.'C')then
	if(dis(1:3).eq.'   ')then	! "dis" is  for a dis_line on pic, 
	write(*,*)'Please use left button first !'
	write(*,'(1x,a)')bell
	goto 110
	endif
	cc=dis(7:30)			! R.A. begin at no. 7 in dis)line
c use last display position as new center	! take out R.A. Decl only
c it is chart_epoch, so change to 2000_epoch 
c repeat principle: display in current epoch, cal in 2000 epoch
	k=itohd(cc,alph_a,delt_a,dummy)
	call astprs(alph_a,delt_a,epoch,alpha,delta,2000.)
	close (99)
	write(*,'(a)')bell
c notice, here must call pgend to end a pic frame
c at first time of this program, no this statment here, so when call pgbegin
c again, prog requir end old_frame, then create new pic_frame immidately
c in this case, hardware of this SUN machine cannot serve, dead sometime. 
	call pgend
	goto 12
	endif

c print********************
	if(ch.eq.'p'.or.ch.eq.'P')then
	write(*,'(2a)')bell
	call pgend
	write(*,*)'Mirror reflect ? [n]'
	im=0
	read(*,'(a)')xc
	if(xc(1:1).eq.'y'.or.xc(1:1).eq.'Y')im=1
153	write(*,*)'HARDCOPY:  frame size in "cm",    [17.06]'
	read(*,'(a)',err=153)xc
	xp=17.06
	if(xc(1:2).ne.'  ')read(xc(1:),*,err=153)xp
	if(xp.lt.0.)xp=-xp
	if(xp.gt.19.5)then
	write(*,*)'Max size is 19.5cm, try again !'
	goto 153
	endif
	if(xp.lt.5.)then
	write(*,*)'I wonder the "Size" is too small, try again!'
	goto 153
	endif
	close (99)
	open(99,file='sao.tmp',status='old',iostat=i)
	call unlink('pgplot.hp')
	icall=pgbegin(0,'/hp',1,1)
	call pgslw(2)
	call pgsch(.8)
	write(*,*)'wait for 1 minute.'
c physical size of laser_printer, 300 dots/inch, 2336 dots in width
c                                 3200 dots in length
	dy=300./xp
	x1=xp/2.54*300./2336.
	y1=xp/2.54*300./3200.
c create pg window with input scale
	x=(1.-x1)*0.5
	x1=(1.+x1)*0.5
	y=0.95
	y1=y-y1
	call pgadvance
	call pgvport(x,x1,y1,y)
	if(im.eq.0)call pgwindow(1.,512.,1.,512.)
	if(im.eq.1)call pgwindow(512.,1.,1.,512.)
	call pgbox('bc',0,0,'bc',0,0)
	i=1				! use 2 sort gray_scale according
	if(w_de.le.30.)i=2		! to chart field size
c 18.+i for i=1, it can print out NGC gray, if choice 20. it cannot print
	n1=512
	n2=512
	if(im.ne.0)call turnx2(map,n1,n2,-360)
	call pgxgray(map,512,512,1,512,1,512,3.,18.+i,i)
	if(im.eq.0)then
	call pglabel('R.A.  South','Decl,  East',label(1:il+11)//'\um\d')
	else
	call pglabel('R.A.  South','Decl,  West',label(1:il+11)//'\um\d')
	endif	
	call pgxchinese(450.,-dy*1.2,q(3),3,0,0)
	call pgxchinese(220.,-dy*1.2,q(1),2,0,0)
	call pgxchinese(-dy*1.2,450.,q(8),3,0,90)
	call pgxchinese(-dy*1.2,50.,q(6),2,0,90)
	call pgxchinese(450.,512.+dy*0.5,q(11),3,0,0)
	call pgline(4,sx,sy)
	call pgtext(450.,15.,sc)
	call pgline(2,crox,croy)
	call pgline(2,crox(3),croy(3))
	call pgsch(0.5)
c here, is too complicate, grtext can make job easy, but lazzy to change
	y1=-4.*dy			! for print display line on paper
	x1=(1.-17./xp)*256.		! in 2 coloumn, cal x position first
	if(im.ne.0)x1=512.-x1
	i=0
150	i=i+1
	read(99,'(a)',end=151)dis
	call dis_pm(dis)
	call pgtext(x1,y1,dis)
	y1=y1-dy
	ns=dis(7:30)
c cal position om star map for mark number above the star
	k=itohd(ns,xa,xd,dummy)
	call astprs(xa,xd,epoch,x_a,x_d,2000.)
	x_a=x_a*cx
	x_d=x_d*cy
	call standc(a_c,d_c,x_a,x_d,xi,xn)
	x=hs2*(hs-xi)+1.
	y=hs2*(xn-hs)+512.
	sg=dis(40:40)
c mark false star
	if(sg.eq.'i')then
	call pgslw(1) 
	call pgpoint(1,x-1,y-1,23)
	call pgslw(2)
	endif
	if(sg.eq.'m')call pgpoint(1,x-1,y-1,6)
	call pgtext(x-3.,y+4.,dis(1:3))
	if(i.eq.10)then
	x1=260.
	if(im.ne.0)x1=250.
	y1=-4.*dy
	endif
	goto 150
151	close (99)
	call pgiden
	call pgend
c submit "print map job" to system, in queeue, print map automatically.
	if(xc(1:1).ne.'-')call system('print pgplot.hp')
	write(*,'(1x,a)')bell
	goto 10
	endif

c right button****************
	if(ch.eq.'X')then
	call pgend
	close (99)
	goto 10
	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=hs-(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.,-20.,dis)
	write(dis(1:),104)num,ac,dc,epoch
	call pgsci(9)
	call pgtext(0.,-20.,dis)
	write(99,'(a)')dis
	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)')ns
	if(itohd(ns,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,epoch)
	call astprs(xa,xd,epoch,x_a,x_d,2000.)
	x_a=x_a*cx
	x_d=x_d*cy
	call standc(a_c,d_c,x_a,x_d,xi,xn)
	x=hs2*(hs-xi)+1.
	y=hs2*(xn-hs)+512.
	call pgsci(2)
	call pgpoint(1,x-1,y-1,18)
	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(*,105)num,ac,dc,epoch
105	format(i3,') ',a,1x,a,'(',f6.1,') input alpha delta')
	call pgsci(0)
	call pgtext(0.,-20.,dis)
	write(dis(1:),105)num,ac,dc,epoch
	call pgsci(9)
	call pgtext(0.,-20.,dis)
	write(99,'(a)')dis
	write(*,"(' Again ? [Y] ',$)")
	read(*,'(a)')sg
	if(sg.ne.'n'.and.sg.ne.'N')goto 120	
	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.)
	inpn=inpn+1
	if(inpn.gt.99)inpn=99
	xya(inpn)=xa
	xyd(inpn)=xd
	ixy(inpn)=num
	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)	
	
	call pgsci(0)
	call pgtext(0.,-20.,dis)
	write(*,1333)num,ac,dc,epoch,k,c5,ch3
1333	format(i3,') ',a,1x,a,'(',f6.1,')      SAO',i7,2(' ',a))
	write(dis(1:),1333)num,ac,dc,epoch,k,c5,ch3
	call pgsci(9)
	call pgtext(0.,-20.,dis)
	write(99,'(a)')dis
	goto 110
	endif

1012    write(*,'(1x,a)')bell
	goto 110	
	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 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+5
	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=.5+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(f1,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)
	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.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
        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