c simple ngsc

	parameter (nn=3999)			! for max number of gsc
	parameter (nnsao=999)			! for max number of sao
	parameter (nnngc=49)			! for max number of ngc
	parameter (nngcv=499)			! for max number of gcv
	parameter (nnpho=99)			! for max number of pho

	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*1 ch				! ch~~ job indicator
	character*1 se,sm,sg                    ! se,sm,sg ~~for display line
	character chgcv1*8,chgcv2*29,chgcv3*8,chgcv4*17
	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 color,c5*5,color1	
	character*3 mon(12)

	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 ch4(4)
	integer*2 m6(6)

	real aangc(nnngc),ddngc(nnngc)		!             for ngc
	integer*2 enngc(nnngc)			! ngc Number
	byte epngc(nnngc),etngc(nnngc)
	integer*2 emngc(nnngc),edngc(nnngc)
	character*3 ngcsp(8)
	integer*2 mngcx(nnngc),mngcy(nnngc)	! position of ngc on pic.

	real aagcv(nngcv),ddgcv(nngcv)		!             for gcv
	integer*4 engcv(nngcv)			! gcv number on publication
	integer*2 emgcv(nngcv),eigcv(nngcv)	! idex of GCV
        integer*2 ibvs(4)
	integer*2 mgcvx(nngcv),mgcvy(nngcv)	! position of gcv on pic.
	integer pgbegin
	real aapho(nnpho),ddpho(nnpho)		!             for pho
	integer*2 enpho(nnpho),empho(nnpho)	! 
	integer*2 mphox(nnpho),mphoy(nnpho)	! position of pho on pic.
	integer*2 m5(5)

	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,log,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)

        data ibvs/2681,3058,3323,3530/
	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/
	data mon/'Jan','Feb','Mar','Apr','May','Jun','Jul',
     *           'Aug','Sep','Oct','Nov','Dec'/
	data ngcsp/'Gx ','OC ','Gb ','Nb ','Pl ','C+N','Ast','Kt '/

	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=30.
	xlimag=17.
        icall=0
	call fdate(cc)
cc      cc=fdate()
cc      character*26 fdate
	read(cc(1:),'(4x,a,i3,9x,i5)')ch3,i,j
           
	do 1009 k=1,12
1009	if(ch3.eq.mon(k))isa=k-1
	j=(j+(isa*30.4+i)/365.)*10.+0.5
	epoch=j*0.1
	ns1='   '
	ns='/EOD/GSC/N0000/0001.BAO'
	call unlink('gsc.tmp')			! gsc.tmp hold line message
	isa=1					! switch sao star
	isv=1					! switch gcv star
	write(*,1001)
1001	format(//' *********** FINDING CHART used by GSC & SAO & NGC
     *  & PHO *********'/60x,'2000.7 Jiang'/)
	icd=0					! if GSC mounted
	inquire(file=ns,exist=logi)
	if(.not.logi)then
	write(*,*)' No "GSC" DISK !'
	write(*,*)
	icd=1
	endif

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**')
	if(icd.eq.0)i=jindex(0.,0.,dir,file1) ! read 'regions.bin' here
	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
	lns=.not.lns
	if(xlimag.lt.5.)goto 12
1000	write(*,9002)nint(w_de)
9002	format('Finding chart size, in arcmin.  [',i3,' ]')
	write(*,*)'if to fit UKST or POSS, choose 200'
	write(*,*)'(4~~~240 min. or >120 for SAO only)'
	read(*,'(a)',err=12)xc
	if(xc(1:2).ne.'  ')read(xc(1:),*,err=1000)w_de
	if(w_de.gt.240.)w_de=240.
	if(w_de.gt.240..and.xlimag.gt.8.)xlimag=8.
	if(w_de.lt.4.)then
	write(*,'(2a)')bell,'too samll field, I change to 4 arcmin'
	w_de=4.
	endif
	write(*,*)'      WAIT !'
	i=w_de*0.25+0.5
	write(sc(1:),'(i3)')i			! display scale on pic.
	sc(4:4)=char(39)
c picture is fixed 512*512 pixels
	sx(3)=500.-i/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***********************
********** 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'
	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)
***********find ngc
	nngc=0	
	call cat_ngc(file,
     *aangc,ddngc,enngc,epngc,etngc,emngc,edngc,nnngc,nngc,mngcx,mngcy)
***********find gcv
	ngcv=0	
	if(isv.eq.1)call cat_gcv(file,
     *    aagcv,ddgcv,engcv,emgcv,eigcv,nngcv,ngcv,mgcvx,mgcvy)
***********find pho
	npho=0	
	if(isv.eq.1)call cat_pho(file,
     *    aapho,ddpho,enpho,empho,nnpho,npho,mphox,mphoy)
**************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 ngc
	do 2050 i=1,nngc
	x=edngc(i)*0.1
	if(x.le.0.)x=1.
	nr=x*256./w_de
	if(nr.lt.4)nr=4
	if(nr.gt.255)nr=255
	nr2=nr*nr
	ix=mngcx(i)
	iy=mngcy(i)
	do 2054 k=0,nr
	do 2054 j=0,nr
	if(k*k+j*j.gt.nr2)goto 2054
	if(inpic(ix+k,iy+j).ne.0)map(ix+k,iy+j)=map(ix+k,iy+j)+1
	if(inpic(ix+k,iy-j).ne.0)map(ix+k,iy-j)=map(ix+k,iy-j)+1
	if(inpic(ix-k,iy+j).ne.0)map(ix-k,iy+j)=map(ix-k,iy+j)+1
	if(inpic(ix-k,iy-j).ne.0)map(ix-k,iy-j)=map(ix-k,iy-j)+1
2054	continue	
2050	continue
c delte overlap
	do 2055 j=1,512
	do 2055 i=1,512
2055	if(map(i,j).gt.4)map(i,j)=4
	write(*,*)'Total NGC',nngc

ccc for gsc
1212	n33=35-20*w_de/120			! according w_de to choice
	if(n33.gt.33)n33=33			! star_map size for gsc
	if(n33.lt.15)n33=15
	nj1=(35-n33)/2
	nj2=34-nj1
	do 50 i=1,n
c release last 2 byte message
	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 star, istar=-1          if non_star, istar=1      as APM defined
	istar=-1
	if(mag.lt.0)istar=1
c to make star map
	call g89(xmag,istar,map33,33,w_de,nj1)
	xm_max=max(xm_max,xmag)
	xm_min=min(xm_min,xmag)
c merge star map to big map
	ix=mapx(i)
	iy=mapy(i)
	do 53 j=nj1,nj2
	jj=iy-17+j-1
	if(jj.lt.1.or.jj.gt.512)goto 53
	do 52 k=nj1,nj2
	ii=ix-17+k-1 
	if(ii.lt.1.or.ii.gt.512)goto 52
	map(ii,jj)=map(ii,jj)+map33(k,j)
52	continue
53	continue
c merge end
50	continue
	write(*,*)'Total gsc',n,'  /',ngsc

ccc for sao
	do 1050 i=1,nsao
	j=emsao(i)
	xmag=j*0.1
	if(xmag.eq.0.)xmag=11.
	if(xmag.lt.-1.65)xmag=-xmag
c if use SAO only, to make star_shape samll, use samll star map for fast cal.
	if(w_de.gt.120.)nj1=xmag+3
	call g89(xmag,-1,map33,33,w_de,nj1)
	xm_min=min(xm_min,xmag)
	xm_max=max(xm_max,xmag)
	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

ccc for gcv
        do 3050 i=1,ngcv
        xmag=emgcv(i)*0.01
	if(xmag.eq.0.)xmag=16.
        call g88(xmag,w_de,nj1)
        xm_min=min(xm_min,xmag)
        if(xmag.lt.25.)xm_max=max(xm_max,xmag)
	i1=0
	i2=0
	ix=mgcvx(i)
	iy=mgcvy(i)
        do 3053 j=-nj1,nj1
	log=.true.
	if(j.eq.nj1.or.j.eq.-nj1)log=.false.
        jj=iy+j-1
        if(jj.lt.1.or.jj.gt.512)goto 3053
        do 3052 k=-nj1,nj1
	if(k.ne.nj1.and.k.ne.-nj1.and.log)goto 3052
        ii=ix+k-1 
        if(ii.lt.1.or.ii.gt.512)goto 3052
        if(map(ii,jj).lt.8.)then
	i1=i1+1
			else
	i2=i2+1
			endif
3052    continue
3053    continue
	i3=0
	if(i1.ge.i2)i3=9
        do 3055 j=-nj1,nj1
	log=.true.
	if(j.eq.nj1.or.j.eq.-nj1)log=.false.
        jj=iy+j-1
        if(jj.lt.1.or.jj.gt.512)goto 3055
        do 3054 k=-nj1,nj1
	if(k.ne.nj1.and.k.ne.-nj1.and.log)goto 3054
        ii=ix+k-1 
        if(ii.lt.1.or.ii.gt.512)goto 3054
	if(iabs(j).eq.nj1.and.iabs(k).eq.nj1)goto 3054
		map(ii,jj)=i3
3054    continue
3055    continue
3050    continue
        write(*,*)'Total GCV',ngcv

ccc for pho
        do 5050 i=1,npho
        xmag=empho(i)*0.001
	if(xmag.lt.0.)xmag=-xmag*10.
        call g88(xmag,w_de,nj1)
        xm_min=min(xm_min,xmag)
        if(xmag.lt.25.)xm_max=max(xm_max,xmag)
	ix=mphox(i)
	iy=mphoy(i)
        do 5055 j=-nj1,nj1
	log=.true.
	if(j.eq.nj1.or.j.eq.-nj1)log=.false.
        jj=iy+j-1
        if(jj.lt.1.or.jj.gt.512)goto 5055
        do 5054 k=-nj1,nj1
	if(k.ne.nj1.and.k.ne.-nj1.and.log)goto 5054
        ii=ix+k-1 
        if(ii.lt.1.or.ii.gt.512)goto 5054
		map(ii,jj)=20
5054    continue
5055    continue
5050    continue
        write(*,*)'Total PHO',npho

	if(n+nsao+nngc+ngcv+npho.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.0)
        open(99,file='gsc.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.,20.,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
	if(isv.eq.-1)call pgtext(518.,270.,' OFF')
********************************************* 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'
	write(*,*)' classification of NGC & IC'
	write(*,*)'  1. Gx  Galaxy'
	write(*,*)'  2. OC  Open star cluster'
	write(*,*)'  3. Gb  Globular star cluster'
	write(*,*)'  4. Nb  Bright emission or reflection nebula'
	write(*,*)'  5. Pl  Planetary nebula'
	write(*,*)'  6. C+N Cluster associated with nebulosity'
	write(*,*)'  7. Ast Asterism or group of a new star'
	write(*,*)'  8. Kt  Knot or nebulous region (external galaxy)'
	
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
	isv=-isv
	if(isv.eq.1)call pgsci(0)
	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(8)
	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]'
	if(w_de.ge.120.)write(*,1544)w_de*6./67.14
1544	format(' for UKST or Palomar plate, 67.14/mm, suggest:',f6.1,'cm')
	write(*,*)'(Ex. 17.06~30*30 for 216),  -value: goto next'
	write(*,*)' 17.06~S2HO4 paper, 17.04~~normal paper'
	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='gsc.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(8)
	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(8)
	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 GCVS (1)
	if(ngcv.eq.0.or.isv.eq.-1)goto 3009
        r=1e10
        do 3011 i=1,ngcv
        rr=(mgcvx(i)-x)**2+(mgcvy(i)-y)**2
        if(rr.lt.r)then
        m=i  
        r=rr
        endif
3011    continue
	if(r.gt.32.)goto 3009
c for display only
        num=num+1
        write(ch3(1:),'(i3)')num
        call pgtext(mgcvx(m)-3.,mgcvy(m)+4.,ch3)
	call astprs(aagcv(m),ddgcv(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)
	xmag=emgcv(m)*0.01
        write(*,3333)num,ac,dc,epoch,engcv(m),xmag
3333    format(i3,') ',a,1x,a,'(',f6.1,') GCV',i7,f7.2,'m')
        call pgsci(0)
        call pgtext(0.,-20.,dis)
        write(dis(1:),3333)num,ac,dc,epoch,engcv(m),xmag
        call pgsci(8)
        call pgtext(0.,-20.,dis)
        write(99,'(a)')dis
c for gcvs.dex
	i=eigcv(m)
	if(i.gt.0)then
        open(19,file='/EOD/gcvs.dex',status='old',
     *       access='direct',recl=62)
	read(19,rec=i)chgcv1,chgcv2,chgcv3,chgcv4
	close (19)
	write(*,'(9x,4(2x,a))')chgcv1,chgcv2,chgcv3,chgcv4
			else
	j=mod(engcv(m),1000)
	i=-i
	write(*,*)'      look supliment vol.',ibvs(i),'   number:',j
			endif
        ip=1

c search pho (2)
3009	if(npho.eq.0.or.isv.eq.-1)goto 3010
        r=1e10
        do 4011 i=1,npho
        rr=(mphox(i)-x)**2+(mphoy(i)-y)**2
        if(rr.lt.r)then
        m=i  
        r=rr
        endif
4011    continue
	if(r.gt.32.)goto 3010
c for display only
        num=num+1
        write(ch3(1:),'(i3)')num
        call pgtext(mphox(m)-3.,mphoy(m)+4.,ch3)
	call astprs(aapho(m),ddpho(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)
	xmag=empho(m)*0.001
	if(xmag.lt.0.)xmag=-xmag*10.
        write(*,3334)num,ac,dc,epoch,xmag
3334    format(i3,') ',a,1x,a,'(',f6.1,') PHO',f7.2,'v')
        call pgsci(0)
        call pgtext(0.,-20.,dis)
        write(dis(1:),3334)num,ac,dc,epoch,xmag
        call pgsci(8)
        call pgtext(0.,-20.,dis)
        write(99,'(a)')dis
c for land.dex,lanz.dex
	i=0
	if(enpho(m).gt.20000)i=2
	if(enpho(m).gt.30000)i=3
	if(i.eq.3)then
         open(19,file='/EOD/sph/land.dex',status='old',
     *       access='direct',recl=10)
	 read(19,rec=enpho(m)-30000)m5
	 close (19)
	endif
	if(i.eq.2)then
         open(19,file='/EOD/sph/lanz.dex',status='old',
     *       access='direct',recl=10)
	 read(19,rec=enpho(m)-20000)m5
	 close (19)
	endif
        call swap2(m5,10)
	if(i.eq.0)then
	 write(*,"(9x,'GSPC: V',f8.2,6x,'B-V',f8.2)")xmag,enpho(m)*0.01
	else
	 write(*,*)'             V    B-V    U-B    V-R    R-I    V-I'
	 if(i.eq.3)write(*,'(a,6f7.3)')'LANDOLT:',xmag,(m5(j)*0.001,j=1,5)
	 if(i.eq.2)then
	  write(*,'(a,f7.3,$)')'  LANZ: ',xmag
	  do j=1,5
	  if(m5(j).eq.9999)write(*,'(a,$)')'       '
	  if(m5(j).ne.9999)write(*,'(f7.3,$)')m5(j)*0.001
	  enddo
	  write(*,*)
	 endif
	endif	
	ip=1

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
        	if(isv.eq.-1)goto 110
        write(*,"('      ',$)")
	if(j.ne.0)then
         open(19,file='/EOD/sao/sao.nam',status='old',
     *       access='direct',recl=8)
	 read(19,rec=j)chgcv1
	 close (19)
	 write(*,"(a,'   ',$)")chgcv1
1334	 format(i3,') ',a,1x,a,'(',f6.1,') ',a,i7,2(' ',a))
	endif
        open(19,file='/EOD/CROSS/hd.dat',status='old',
     *       access='direct',recl=3)
	read(19,rec=k)ch3
	close (19)
	j=ichar(ch3(1:1))*10000+ichar(ch3(2:2))*100+ichar(ch3(3:3))
	if(j.ne.0)write(*,"('HD',i6,'   ',$)")j
        open(19,file='/EOD/CROSS/dm.dat',status='old',
     *       access='direct',recl=4)
	read(19,rec=k)ch4
	close (19)
        j=ch4(1)
        if(j.ne.0)then
         ch3(3:3)='+'
         if(j.lt.0)ch3(3:3)='-'
         if(j.lt.0)j=-j
         if(j.eq.1)ch3(1:2)='BD'
         if(j.eq.2)ch3(1:2)='CD'
         if(j.eq.3)ch3(1:2)='CP'
         j=ch4(3)
         if(j.lt.0)j=256+j
         j=100*j+ch4(4)
	 write(*,'(a,i2.2,i5.4,$)')ch3,ch4(2),j
	endif
        write(*,*)
        i=muvby(k,m6)
        if(i.ne.0)then
         if(i.gt.9)i=9
         write(*,'(4(a,f6.3,3x),a,$)')'      v:',m6(1)*0.001,
     * 'b-y:',m6(2)*0.001,'m1:',m6(3)*0.001,'c1:',m6(4)*0.001,'beta:'
         if(m6(5).ne.0)write(*,'(f6.3,$)')m6(5)*0.001
         write(*,"('   Ref:',i4,'-',i1)")mod(m6(6),1000),i
        endif
	call pgsci(8)
	call pgtext(0.,-20.,dis)
	write(99,'(a)')dis
	ip=1			! if success, not search gsc, so some	

c search gsc next (4)
1012    r=1e10
	do 101 i=1,n
	rr=(mapx(i)-x)**2+(mapy(i)-y)**2
	if(rr.lt.r)then
	m=i
	r=rr
	endif
101	continue
	if(r.gt.32.)goto 1013
	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
	write(ch3(1:),'(i3)')num
	call pgtext(mapx(m)-3.,mapy(m)+4.,ch3)
	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.)
	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)
	call pgsci(0)				! earse existed words
	call pgtext(0.,-20.,dis)
	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)
	call pgsci(2)
	call pgpoint(1,mapx(m)-.9,mapy(m)-.9,-1)
	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
	call pgsci(8)
	call pgtext(0.,-20.,dis)
	write(99,'(a)')dis
	write(*,'(a)')dis
	endif
	ip=1


c search NGC (5)
1013	r=1e10
	do 2014 i=1,nngc
	rr=(mngcx(i)-x)**2+(mngcy(i)-y)**2
	if(rr.lt.r)then
	m=i
	r=rr
	endif
2014	continue
	if(r.gt.32.)goto 4012
c for display only
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgtext(mngcx(m)-3.,mngcy(m)+4.,ch3)
	call astprs(aangc(m),ddngc(m),2000.,xa,xd,epoch)
	xa=amod(xa+24.,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	xmag=emngc(m)*0.1
	if(xmag.lt.0.)xmag=-xmag
	write(c5(1:),"(f4.1,'v')")xmag
	if(emngc(m).lt.0)c5(5:5)='p'
	if(emngc(m).eq.0)c5='NoMag'
	call pgsci(0)
	call pgtext(0.,-20.,dis)
	ch3='NGC'
	k=enngc(m)
	if(k.lt.0)then
	 k=-k
	 ch3=' IC'
	endif
	j=epngc(m)
	xc=ngcsp(j)
	if(etngc(m).ne.0)write(xc(4:),"(' M',i3)")etngc(m)
	write(*,2335)num,ac,dc,epoch,ch3,k,c5,xc
2335	format(i3,') ',a,1x,a,'(',f6.1,') ',a,i4,1x,a,1x,a)
	write(dis(1:),1334)num,ac,dc,epoch,ch3,k,c5,xc
	call pgsci(8)
	call pgtext(0.,-20.,dis)
	write(99,'(a)')dis
	ip=1
4012    if(ip.eq.0)write(*,'(1x,a)')bell
	goto 110	
	end

********** uvby
        function muvby(k,m6)
        integer*4 a(4,32),a3(3),k4,k
        integer*2 a256(256),a6(6),m6(6)
        equivalence (a3,a6)
        data ip/0/
	data a256/256*0/
        open(19,file='/EOD/sph/uvby.bin',status='old',
     *       access='direct',recl=512)
        if(ip.eq.0)then
	  read(19,rec=1)a256
	  call swap2(a256,512)
	  ip=1
	endif
        muvby=0
        j=k/1000
        if(j.lt.1)j=1
        if(j.gt.256)j=256
        j=a256(j)
20      read(19,rec=j)a
	call swap4(a,512)
        do 30 i=1,32
        k4=a(1,i)
        if(k.lt.k4)goto 40
        if(k.gt.k4)goto 30
        a3(1)=a(2,i)
        a3(2)=a(3,i)
        a3(3)=a(4,i)
        do 10 m=1,6
10      m6(m)=a6(m)
        muvby=muvby+1
30      continue
        j=j+1
        goto 20
40	close(19)
        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 g88(xmag,ww,nj1)
c according xmag & ww to determine box size
c ww:  chart field width, in arcmin
c nj1: start number of star_map
	x=xmag
	if(x.gt.16.)x=16.
	w=ww
	if(w.gt.120)w=120.
	nj1=12-w/30-x/2
	if(nj1.lt.2)nj1=2
	end
	
	subroutine g89(xxmag,istar,a,n,ww,nj1)
c make false star and non_star, output is star_map "a"
c xxmag: input, star mag
c istar: =1 means non_star, =-1 means star
c a(n,n) n=33 star_map, in many case, not use all in order to save time
c n = 33 forever
c ww:  chart field width, in arcmin
c nj1: start number of star_map
c c:   hollow in center or not
	integer*2 a(n,n)
	real y(33)
	xmag=xxmag
	w=ww
	if(w.gt.120.)then	! when in > 120 arcmin chart, draw SAO
	xmag=xmag+4.		! only, to make star smaller and less
	w=w/4.			! cross spark in center lines
	endif
	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=60./w+15.-xmag	! for sigma, accoring mag & w_de & sao only
	if(w.le.8.)x=60./w+(15.-xmag)*2.	! make big gap in samll field
c for star
	if(istar.eq.-1)then
	if(xmag.lt.6.)x=((6.-xmag)*0.3+1.)*x	! make < 6m star more bright
	sigma=.5+x*.25
	call gauss(y,n,sigma)			! get guass curve
	rr=max((sigma*0.5+0.5)/sqrt(40./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
c for galaxic objects
	else	
	high=10.			! not dark in center for galaxy.
	sigma=1.+x*.3
	call gauss(y,n,sigma)
	do 130 j=nj1,17
	jj=(j-17)**2
	do 130 i=nj1,nj2
	r=jj+(i-17)**2+(i-17)*(17-j)		! (i-17)(17-j) make ecllipe
	kr=sqrt(r)+1.5				! shape of non_star
	r=y(kr)*high			
	if(r.gt.0.5.and.r.lt.1.)r=1.		! make big profile of 
	if(r.gt.9.)r=9.
	a(i,j)=r
130	a(34-i,34-j)=r
	endif
	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 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

	subroutine cat_gcv(f1,aa,dd,en,em,ei,nn,n,mx,my)
        real aa(nn),dd(nn),en(nn)
	integer*2 em(nn),ei(nn),mx(nn),my(nn)
        character*3 f1 
        real a(4,32)
        integer*2 ind(256),bb(2)
        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),(ee,mg)
        open(1,file='/EOD/gcv/'//f1//'.GCV',
     *       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)
	 call swap2(a(4,j),4)
        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)
	if(bb(1).gt.10000)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)
        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)
	 ee=a(4,j)
	 em(n)=bb(1)
	 ei(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_pho(f1,aa,dd,en,em,nn,n,mx,my)
        real aa(nn),dd(nn)
	integer*2 en(nn),em(nn),mx(nn),my(nn)
        character*3 f1 
        real a(3,32)
        integer*2 ind(192),bb(2)
        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)
        open(1,file='/EOD/sph/'//f1//'.PHO',
     *       status='old',access='direct',recl=384)
        read(1,rec=1)ind
	 call swap2(ind,384)
        x1=dn_a
        x2=up_a
        k1=x1*5.
	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(192)
        read(1,rec=i)a
        do 20 j=1,32
	 call swap2(a(1,j),4)
	 call swap4(a(2,j),8)
        if(a(2,j).gt.x2.and.i87.eq.0)goto 100
        if(i.eq.ind(192).and.a(1,j).eq.0)goto 30
        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)
        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=a(1,j)
	 en(n)=bb(1)
	 em(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_ngc(f1,aa,dd,en,ep,et,em,ed,nn,n,mx,my)
	real aa(nn),dd(nn)
	integer*2 en(nn),em(nn),ed(nn),mx(nn),my(nn)
	byte ep(nn),et(nn),b(4)
	character*3 f1
        real a(4,32)
        integer*2 ind(256),bb(2)
        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),(ee,mg),(bb,b)
        mag=xlimag*100.
        open(1,file='/EOD/ngc/'//f1//'.NGC',
     *       status='old',access='direct',recl=512)
        read(1,rec=1)ind
	 call swap2(ind,512)
	if(ind(256).eq.2)goto 100
        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 swap2(a(1,j),2)
	 call swap4(a(2,j),8)
	 call swap2(a(4,j),4)
        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
        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)
        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=a(1,j)
	 en(n)=bb(1)
	 et(n)=b(3)
	 ep(n)=b(4)
	 ee=a(4,j)
	 em(n)=bb(1)
	 ed(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	
