	character*30 f1,cc*1
	integer*2 b(2048*2048)                  ! limit 2048,2048 CCD
 	integer*2 a(2048*2048*2)		! -32 change to 16 when cal.
	character*80 head(180)
	real*8 a8(8)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	pi=4.*atan(1.0)
	cx=pi/12.
	cy=pi/180.
  	phi=(40.+23/60.+36/3600.)*pi/180.	! XL latitude
  	sphi=sin(phi)
	cphi=cos(phi)
c ! means: force to do a8(8) again 
	k=iargc()
        if(k.eq.0)stop 'usage:  tv_world filename [!]'
c check file
        call getarg(1,f1)
        if(inqh(f1,n1,n2,head,ibit,ihead,lhead).eq.0)
     *    stop 'DATA_file not found'
	call getpar(head,lhead)				! get s0,iy,im,id,ut
c get 6 plate coef. and plate center & XL local sidereal time of OBS.
	ip=0
	if(k.ge.2)then
	  call getarg(2,cc)
	  if(cc.eq.'!')goto 10
	endif
	if(indexpos(head,'A81     ',lhead).ne.lhead)then
	  call dis_5(head,a8,lhead,n1,n2)
	  ip=1
	endif
c read data
10	call readata(a,b,f1,ibit,sky,n1,n2)
        m1=512
        m2=512
        ipos=indexpos(head,'INSTRUME',lhead)
        if(ipos.ne.73)then
          do i=12,20
          if(head(ipos)(i:i).eq.'3')iturn=0
          if(head(ipos)(i:i).eq.'4')iturn=4
          if(head(ipos)(i:i).eq.'6')iturn=6
          if(head(ipos)(i:i).eq.'7')iturn=7
          enddo
        endif
        if(iturn.eq.0)call turnx2(b,m1,m2,-180)
        if(iturn.eq.1)call turnx2(b,m1,m2,-360)
        if(iturn.eq.4)call turnx2(b,m1,m2,180)
        if(iturn.ge.6)call turnx2(b,m1,m2,-180)

c main_job
	call main_job(head,a8,ip,b,a,f1,match,rms,seeing,
     *  n1,n2,lhead)
	if(ip.eq.1)stop
c write something to fits head
	call put_a8(head,a8,lhead) 	        ! 45-52	 write down a8
20	call put_gr(head,lhead)		! 53-57  write gal_co. redden
	call put_xmass(head,lhead)	! 15-17, 58 write center,h_angle,airmass
	call put_aa(head,lhead)		! 59,60  write Azimuth, Altitude
	call put_moon(head,a8,lhead)	! 61-65  phase, position, pos_angle
c	ipos=indexpos(head,'COMMENT ',lhead)
c	head(ipos)(1:32)='CRVAL1  =                    0 /'  
c	ipos=indexpos(head,'COMMENT ',lhead)
c	head(ipos)(1:32)='CRVAL2  =                    0 /'  
c	ipos=indexpos(head,'CTYPE1  ',lhead)
c	if(ipos.ne.lhead)head(ipos)(1:8)='COMMENT '
	close (1)
	write(*,*) ' re_write fits head !'
	i=lnblnk(f1)+1
	f1(i:i)=char(0)
	call tvcoordc(f1,head)
	end

	subroutine putdata2(a,b)
	integer*2 a(1),b(1)
	do 10 i=1,2048*2048
10	b(i)=a(i)
	end

        function see(map,n1,n2,x,y,white)
c (a81+a84)/2 *180*3600/pi = arcsec / pixel
c schmidt is 1.709"/pixel
c seeing= 2.355*sigma*1.709
        integer*2 map(n1,n2)
        integer b(113)        ! (15-1)*8+1
        ix=x+0.5 
        iy=y+0.5
        j=1
        do 10 i=ix-7,ix+7
        b(j)=map(i,iy)-white
10      j=j+8    
        do 20 j=1,112,8
20      b(j+4)=(b(j)+b(j+8))*0.5+0.5
        do 30 j=1,112,4
30      b(j+2)=(b(j)+b(j+4))*0.5+0.5
        do 40 j=1,112,2
40      b(j+1)=(b(j)+b(j+2))*0.5+0.5
        call histat(b,mode,maxh,xmean,xpeak,fsigma1,gsigma,113)
        j=1
        do 110 i=iy-7,iy+7
        b(j)=map(ix,i)-white
110     j=j+8
        do 120 j=1,112,8
120     b(j+4)=(b(j)+b(j+8))*0.5+0.5
        do 130 j=1,112,4
130     b(j+2)=(b(j)+b(j+4))*0.5+0.5
        do 140 j=1,112,2
140     b(j+1)=(b(j)+b(j+2))*0.5+0.5
        call histat(b,mode,maxh,xmean,xpeak,fsigma2,gsigma,113)
c       write(*,'(a,f9.2,a)')'seeing:',(fsigma1+fsigma2)/16.*2.355*1.709,'"'
        see=(fsigma1+fsigma2)/16.*2.355*1.709
        end

	subroutine getmxy(a8,x,y,z,phase,pl,pb,xx)
	real*8 a8(8),b8(12),jd
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	call xytoad(a8,b8)
	j=(iy+((im-1)*30.4+id)/365.)*10.+0.5
	epoch_g=j*0.1
	x=a8(7)/cx
	y=a8(8)/cy
	call astprs(x,y,2000.,alpha_g,delta_g,epoch_g)
	call ad_azl(alpha_g,delta_g,a1,d1)
	alpha_g=alpha_g*cx
	delta_g=delta_g*cy
	call azl_ad(a1,d1-0.1,x,y)
	call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	a3=b8(1)*xi+b8(3)*xn+b8(5)
	d3=b8(2)*xi+b8(4)*xn+b8(6)
	call azl_ad(a1,d1+0.1,x,y)
	call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	a2=b8(1)*xi+b8(3)*xn+b8(5)
	d2=b8(2)*xi+b8(4)*xn+b8(6)
	z=atan2(a2-a3,d2-d3)/cy
	call fjd(jd,iy,im,id)
	jd=jd+ut/24.
	call moon2(jd,a2,d2)				! get monn l,b
	call lb_ad(a2,d2,x,y)				! get moon alpha delta
	phase=amod(a2-sunlamda(jd)+360.,360.)/360.*29.53  ! get moon phase
	call ad_azl(x,y,a2,d2)				! get moon A, H
	pl=a2
	pb=d2
	x=sin(d1*cy)*sin(d2*cy)+cos(d1*cy)*cos(d2*cy)*cos(a1*cy-a2*cy)
        y=acos(x)                                    
	xx=y/cy
	cosA=(sin(d2*cy)-sin(d1*cy)*x)/cos(d1*cy)/sin(y)
	if(cosa.gt.0.999999)cosa=0.999999
	if(cosa.lt.-.999999)cosa=-.999999
	y=acos(cosA)/cy
	x=a2-a1
	if(abs(x).gt.180.)x=a1-a2
	if(x.lt.0.)y=-y
	z=270.+z-y
	z=amod(z,360.)
	if(z.lt.0.)z=z+360.
	x=cos(z*cy)
	y=sin(z*cy)
	end

	subroutine dis_5(head,a8,lhead,n1,n2)
	character*80 head(180),ac*12,dc*12
	real*8 a8(8)
	ipos=indexpos(head,'A81     ',lhead)-1
        
	do 10 i=1,8
10	read(head(ipos+i)(11:),*)a8(i)
        call xy_rad(1.,1.,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,1)
1	format(5x,'Letf_up corner:  ',35x,'Right_up corner:')
	write(*,2)ac,dc
2	format(a,1x,a,'                          ',$)
        call xy_rad(float(n1),1.,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,3)ac,dc
3	format(a,1x,a/28x,'plate center: (2000.0)')
        call xy_rad(n1/2+.5,n2/2+.5,a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,4)ac,dc
4	format(26x,a,1x,a/
     * 4x,'Left_down corner:',35x,'Right_down corner:')
        call xy_rad(1.,float(n2),a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,2)ac,dc
        call xy_rad(float(n2),float(n1),a8,alpha,delta)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,'(a,1x,a)')ac,dc
	end

	subroutine azl_ad(a,h,alpha,delta)	!dd->hd
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	aa=a*cy
	hh=h*cy
	sina=sin(aa)
	cosa=cos(aa)
	sinh=sin(hh)
	cosh=cos(hh)
	sinde=sphi*sinh-cphi*cosa*cosh
	cosde=sqrt(1.-sinde*sinde)
	sint=cosh*sina/cosde
	cost=(cosa*cosh+cphi*sinde)/(sphi*cosde)
        delta=atan(sinde/cosde)/cy
        t=atan(sint/cost)
        if(cost.lt.0.)t=t+pi
        alpha=s0-t/cx	
	if(alpha.lt.0.)alpha=alpha+24.
	if(alpha.gt.24.)alpha=alpha-24.
	end

	subroutine ad_azl(alpha,delta,a,h)     ! hd,dd
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        hour=(s0-alpha)*cx
        delt=delta*cy
        s_delta=sin(delt)
	c_delta=cos(delt)
        s_hour=sin(hour)
	c_hour=cos(hour)
        sinh=sphi*s_delta+cphi*c_delta*c_hour
        cosh=sqrt(1.-sinh*sinh)
        sinA=c_delta*s_hour/cosh
        cosA=(sphi*c_delta*c_hour-cphi*s_delta)/cosh
        h=atan(sinh/cosh)/cy
        a=atan(sinA/cosA)
        if(cosA.lt.0.)a=a+pi
	if(a.lt.0.)a=a+pi+pi
	if(a.gt.pi+pi)a=a-pi-pi
	a=a/cy
	end

	subroutine put_aa(head,lhead)     ! 59,60
	character*80 head(180)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        ipos=indexpos(head,'HA      ',lhead)
	read(head(ipos)(19:),'(i2,1x,i2,1x,f5.2)')i,j,x
	t=i+j/60.+x/3600.
	ipos=indexpos(head,'DEC     ',lhead)
	read(head(ipos)(20:),'(i2,1x,i2,1x,f4.1)')i,j,x
	d=i+j/60.+x/3600.
	if(head(ipos)(19:19).eq.'-')d=-d
	call ad_azl(s0-t,d,a,h) 		    ! hd,dd
	ipos=indexpos(head,'A81     ',lhead)+14           ! 59 row
        head(ipos)(1:32)='AZIMUTH =                      /'
        head(ipos+1)(1:32)='ALTITUDE=                      /'
 	head(ipos)(34:76)='Azimuth is measured from south through west'
 	head(ipos+1)(34:54)='both units in degrees'
	write(head(ipos)(21:30),'(f10.2)')a        ! in degree
	write(head(ipos+1)(21:30),'(f10.2)')h
	end

	subroutine put_moon(head,a8,lhead)		! 61-65
	character*80 head(180)
	real*8 a8(8)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',lhead)+16           ! 61 row
        head(ipos)(1:32)='MPHASE  =                      /'
        head(ipos+1)(1:32)='MAZIMUTH=                      /'
        head(ipos+2)(1:32)='MALITIUD=                      /'
        head(ipos+3)(1:32)='MANGLE  =                      /'
        head(ipos+4)(1:32)='MDIRECT =                      /'
	head(ipos)(34:71)='Define the length of Moon_Month 29.53d'
 	head(ipos+1)(34:78)='Moon Azimuth is measured from south thr. west'
 	head(ipos+2)(34:54)='both units in degrees'
	head(ipos+3)(34:75)='Position angle of Moon to CCD_Field center'
	head(ipos+4)(34:77)='Moon direction in the X_Y plane of CCD frame'
	call getmxy(a8,x,y,z,phase,pl,pb,xx)
	write(head(ipos)(21:30),'(f10.1)')phase
	write(head(ipos+1)(21:30),'(f10.2)')pl
	write(head(ipos+2)(21:30),'(f10.2)')pb
	write(head(ipos+3)(21:30),'(f10.2)')xx		
        write(head(ipos+4)(21:30),'(f10.2)')z	
	end
	
	function sunlamda(jd)			! in degree
c sun lamda in degree, beta==0
	real*8 jd,n,l,g,x
	n=jd-2451545.d0
	l=280.460d0+0.9856474d0*n
	g=357.528d0+0.9856993d0*n
	x=dmod(g,360d0)*datan(1d0)/45d0
	x=l+1.915d0*dsin(x)+0.020d0*dsin(x+x)
	x=dmod(x,360d0)
	if(x.lt.0d0)x=x+360d0
	sunlamda=x
	end

	subroutine moon2(jd,pl,pb)		!in degree
	real*8 jd,t,g,r,x1,x2,x3,x4,x5,x6
	g=datan(1d0)/45d0
	t=(jd-2451545.d0)/36525.d0
	r=218.32d0+dmod(481267.833d0*t,360d0)
	x1=(134.9d0+dmod(477198.85d0*t,360d0))*g
	x2=(259.2d0-dmod(413335.38d0*t,360d0))*g
	x3=(235.7d0+dmod(890534.23d0*t,360d0))*g
	x4=(269.9d0+dmod(954397.70d0*t,360d0))*g
	x5=(357.5d0+dmod( 35999.05d0*t,360d0))*g
	x6=(186.6d0+dmod(966404.05d0*t,360d0))*g
	r=r+6.29d0*dsin(x1)-1.27d0*dsin(x2)+0.66d0*dsin(x3)
     *     +0.21d0*dsin(x4)-0.19d0*dsin(x5)-0.11d0*dsin(x6)
	r=dmod(r,360d0)
	if(r.lt.0d0)r=r+360d0
	pl=r
	x1=( 93.3d0+dmod(483202.03d0*t,360d0))*g
	x2=(228.2d0+dmod(960400.87d0*t,360d0))*g
	x3=(318.3d0+dmod(  6003.18d0*t,360d0))*g
	x4=(217.6d0-dmod(407332.20d0*t,360d0))*g
	r=5.13d0*dsin(x1)+0.28d0*dsin(x2)-0.28d0*sin(x3)-0.17d0*dsin(x4)
	r=dmod(r,360d0)
	if(r.gt.90d0)r=r-360d0
	pb=r
	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 ad_lb(al,de,pl,pb)		!hd-->dd
	cy=atan(1.)/45.
	cx=cy*15.
	xl=cos(de*cy)*cos(al*cx)
	xm=0.9175*cos(de*cy)*sin(al*cx)+0.3978*sin(de*cy)
	xn=-.3978*cos(de*cy)*sin(al*cx)+0.9175*sin(de*cy)
	pl=atan2(xm,xl)/cy
	if(pl.lt.0.)pl=pl+360.
	pb=asin(xn)/cy
	end

	subroutine put_a8(head,a8,lhead)      !45-52
	character*80 head(180),ac*12,dc*12
	real*8 a8(8)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',lhead)
	if(ipos.eq.lhead)ipos=indexpos(head,'        ',lhead)
	if(ipos.gt.51)ipos=indexpos(head,'COMMEMT ',lhead)
	if(ipos.gt.51)ipos=51
	if(ipos.eq.51)then
	  do 10 i=51,71
	  head(i)(30:60)='                              '
10	  head(i)(61:80)='                    '
	endif
	head(ipos)(1:46)='A81     =                      / 6 plate coef.'
	head(ipos+1)(1:46)='A82     =                      / matched star:'
	head(ipos+2)(1:46)='A83     =                      / RMS (arcsec):'
	head(ipos+3)(1:32)='A84     =                      /'
	head(ipos+4)(1:32)='A85     =                      /'
	head(ipos+5)(1:32)='A86     =                      /'
	head(ipos+6)(1:46)='A87     =                      / p_c (2000.0) '
	head(ipos+7)(1:46)='A88     =                      / in rad.      '
	do 11 i=1,6
11	write(head(ipos-1+i)(11:30),'(E20.12)')a8(i)
	write(head(ipos+6)(11:30),'(f20.7)')a8(7)
	write(head(ipos+7)(11:30),'(f20.7)')a8(8)
	x=a8(7)/cx				! put lamda beta
	y=a8(8)/cy
	call ad_lb(x,y,alpha,delta)
	call toms(alpha,ac,2)
	call toms(delta,dc,0)
	head(ipos+3)(57:63)='l_sun ='
	head(ipos+4)(57:63)='b_sun ='
	head(ipos+3)(65:73)=ac(1:9)
	head(ipos+4)(65:73)=dc(1:9)
	end

	subroutine put_gr(head,lhead)          !51:53-57
	character*80 head(180)
	character*12 ac,dc
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',lhead)+6
	read(head(ipos)(21:),*)xa			! center plate in rad.
	read(head(ipos+1)(21:),*)xd
	xa=xa/cx
	xd=xd/cy				! to h,d
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	head(ipos+2)(1:46)="RA2000  = '                  ' / plate center "
	head(ipos+3)(1:32)="DEC2000 = '                  ' /"
	head(ipos+2)(18:29)=ac
	head(ipos+3)(19:29)=dc(1:11)
        call astprs(xa,xd,2000.,x,y,1950.)
        call galactic(x*cx,y*cy,xa,xd)
        head(ipos+4)(1:32)='GALLONG =                      /'
        head(ipos+4)(34:62)='galactic coordinates (1950.0)'
        head(ipos+5)(1:32)='GALLATI =                      /'
        head(ipos+5)(34:54)='both units in degrees'
	write(head(ipos+4)(21:30),'(f10.2)')xa        ! in degree
	write(head(ipos+5)(21:30),'(f10.2)')xd
	head(ipos+6)(1:32)='EXTINCT =                      /'
	head(ipos+6)(34:63)='BH Bmag, EXTINCTION = 4*E(B-V)'
	end

	subroutine put_xmass(head,lhead)          !15,16,17,58
	character*80 head(180)
	character*12 ac,dc
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',lhead)+6
	read(head(ipos)(21:),*)xa			! center plate in rad.
	read(head(ipos+1)(21:),*)xd
	xa=xa/cx
	xd=xd/cy					! to h,d
        j=(iy+((im-1)*30.4+id)/365.)*10.+0.5
	epoch=j*0.1
        call astprs(xa,xd,2000.,x,y,epoch)
	call toms(x,ac,1)
	call toms(y,dc,0)
	ipos=indexpos(head,'RA      ',lhead)
	head(ipos)(11:30)="'                  '"
	head(ipos)(18:29)=ac
	write(head(ipos)(50:),"('(',f6.1,')')")epoch
	ipos=indexpos(head,'DEC     ',lhead)
	head(ipos)(11:30)="'                  '"
	head(ipos)(19:29)=dc(1:11)
	write(head(ipos)(50:),"('(',f6.1,')')")epoch
	h=amod(s0-x,24.)				! get HA, here I wast
        if(h.lt.0.)h=h+24.				! more time, x->xa
	call toms(h,ac,1)
	ipos=indexpos(head,'HA      ',lhead)
	head(ipos)(11:30)="'                  '"
	head(ipos)(18:29)=ac
	head(ipos+1)(1:23)="EPOCH   =              "
	write(head(ipos+1)(24:30),'(f7.1)')epoch
	xmm=xmass(h*cx,xd*cy)				! middle time
	ipos=indexpos(head,'EXPOSURE',lhead)
	read(head(ipos)(20:),*)k				! exp time
	xmb=xmass((h-k/7200.)*cx,xd*cy)			! begin time
	xme=xmass((h+k/7200.)*cx,xd*cy)			! end time
	xme=xme             !!!!
	x=(xmb+4.*xmm+xme)/6.
	ipos=indexpos(head,'A81     ',lhead)+13
	head(ipos)(1:32)='AIRMASS =                      /'
	write(head(ipos)(21:30),'(f10.3)')x
	end

	subroutine getpar(head,lhead)
c get s0 iy im id ut
	character*80 head(180)
	real*8 d1,d
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'DATE-OBS',lhead)
	read(head(ipos)(22:),'(i2,1x,i2,1x,i2)',err=101)id,im,iy ! read obs date
	iy=iy+1900
	if(iy.lt.1950)iy=iy+100
102	ipos=indexpos(head,'TIME    ',lhead)
	read(head(ipos)(20:),'(i2,1x,i2,1x,f4.1)',err=101)i,j,x
        ipos=indexpos(head,'EXPOSURE',lhead)
	read(head(ipos)(20:),*)k				! exp time
ccc 010329 due to some expure in milli_sec unit
	k=0
	ut=i+j/60.+x/3600.+k/7200.			! center UT time, in h
        call fjd(d1,1992,12,31)			! 1992,12,31
        d2=(6+38/60.+40.1954/3600.)/24.         ! 6 38 40.1954
        x1=7+50/60.+18.344/3600.                ! 7:50:18.344 lamda of XL
        call fjd(d,iy,im,id)
        d=(d-d1)*1.0027379093d0+d2
        d=dmod(d,1d0)
        if(d.lt.0d0)d=d+1d0
        s0=d*24d0
        s0=s0+x1+ut*1.002738
	return
101     write(*,*)'***** unexpect data format *******'
        stop
	write(*,*)head(ipos)(1:30),'  please input yyyy,mm,dd'
	read(*,*)iy,im,id
	head(ipos)(10:32)='                      /'
	write(head(ipos)(22:30),"(i2.2,'/',i2.2,'/',i2.2)")id,im,mod(iy,100)
	goto 102
	end

        subroutine fjd(d,iy,im,id)
        real*8 d
        i=(im+9)*0.09
        i=(i+iy-1900)*1.75+0.01
        j=30.56*im
        j=(iy-1950)*367+id+j-i
        d=2433338.5d0+j
        end

	function xmass(angle,delta)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
  	x=sin(delta)*sphi+cos(delta)*cos(angle)*cphi
c  	z=90.-atan2(x,sqrt(1.-x*x))*180./pi      ! Zenth
  	x=1./x-1.
  	xmass=1.+x-x*(0.0018167+x*(0.002875+x*0.0008083))
	end

	subroutine put_a2(head,lhead)
	character*80 head(180)
	character*35 dis
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
	ipos=indexpos(head,'A81     ',lhead)
	if(ipos.eq.lhead)ipos=indexpos(head,'        ',lhead)
	if(ipos.gt.51)ipos=indexpos(head,'COMMEMT ',lhead)
	if(ipos.gt.51)ipos=51
	head(ipos)(1:32)='A81     =                      /'
        head(ipos)(34:53)='6 plate coefficients'
	head(ipos+1)(1:32)='A82     =                      /'
	head(ipos+2)(1:32)='A83     =                      /'
	head(ipos+3)(1:32)='A84     =                      /'
	head(ipos+4)(1:32)='A85     =                      /'
	head(ipos+5)(1:32)='A86     =                      /'
	head(ipos+6)(1:46)='A87     =                      / p_c (2000.0) '
	head(ipos+7)(1:46)='A88     =                      / in radians   '
10	write(*,*)'input center R.A, Decl, [1950.0]: ex. 3:46:02.1 23:0'
	read(*,'(a)')dis
	if(itohd(dis,x,y,epoch).eq.0)goto 10
	if(epoch.eq.0.)epoch=1950.
	call astprs(x,y,epoch,x,y,2000.)
	write(head(ipos+6)(11:30),'(f20.7)')x*cx
	write(head(ipos+7)(11:30),'(f20.7)')y*cy
	end

        subroutine galactic(alpha,delta,xl,xb)
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        a0=282.25*cy
        theta=62.6*cy
        Xl0=33.0*cy
        x1=cos(delta)*cos(alpha-a0)
        x2=cos(delta)*sin(alpha-a0)*cos(theta)+sin(delta)*sin(theta)
        xl=atan2(x2,x1)+xl0
        if(xl.lt.0.)xl=xl+2.0*pi
        xb=asin(sin(delta)*cos(theta)-
     *          cos(delta)*sin(alpha-a0)*sin(theta))
        xl=xl/cy
        xb=xb/cy
        end

	subroutine readata(a,b,f1,ibit,sky,n1,n2)
	integer c(720),a(1),b(1)
	character*30 f1,ch*2880
	equivalence (c,ch)
	open(1,file=f1,status='old',access='direct',recl=2880)
	read(1,rec=1)ch
	n=1
	nn=720
55	k=index(ch,'END     ')
	if(k.eq.0)then
	  n=n+1
	  read(1,rec=n)ch
	  goto 55
	endif
	write(*,"(' reading data ')")
	i=0
	if(ibit.lt.0)ibit=-ibit
	k1=n+1
	k2=n+(n1*n2*ibit/32-1)/nn
	do 60 k=k1,k2
	if(mod(k,120).eq.0)call dispdot()
	read(1,rec=k)ch
	do 60 j=1,nn
	i=i+1
60	a(i)=c(j)
	close (1)
	if(ibit.eq.32)then
	  call swap4(a,n1*n2*4)
	  call tob4(a,b,n1*n2)
	else
	  call swap2(a,n1*n2*2)
	  do 10 i=1,n1*n2/2
10	  b(i)=a(i)
	endif
	call whitexblack(b,n1,n2,sky,sigma)
	write(*,*)
	end

	subroutine tob4(a,b,n)
	real a(1)
	integer*2 b(1)
	do 10 i=1,n
	x=a(i)+0.5
	if(x.lt.0)x=0
	if(x.gt.32767.)x=32767.
10	b(i)=x
	end

	function inqh(f1,n1,n2,head,ibit,ihead,lhead)
	character*(*) f1
	character*80 head(180)
	logical logi
	k=index(f1,'.')
	if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
	inquire(file=f1,exist=logi)
	if(.not.logi)return
	open(1,file=f1,status='old',access='direct',recl=80)
	lhead=73
	do 6 i=1,180
	read(1,rec=i)head(i)
	if(lhead.eq.0.and.head(i)(1:8).eq.'        ')lhead=i
	if(head(i)(1:8).eq.'END     ')goto 7
6	continue
	write(*,*)' not a fits file'
	inqh=0
	return
7	ihead=(i-1)/36+1


	read(head(2)(22:),*)ibit
	read(head(4)(22:),*)n1	!NAXIS1
	read(head(5)(22:),*)n2	!NAXIS2

	ipos=indexpos(head,'OBJECT  ',lhead)
	write(*,'(1x,a)')head(ipos)(1:70)       	! object
	ipos=indexpos(head,'DATE-OBS',lhead)
        head(ipos)(63:72)='(dd/mm/yy)' 
	do 10 i=ipos,ipos+5				! da.ti.ex.ra.de.ha
10	write(*,'(1x,a)')head(i)(1:78)
	ipos=indexpos(head,'EPOCH   ',lhead)
        if(ipos.ne.lhead)write(*,'(1x,a)')head(ipos)(1:78)
	close (1)
	end

	subroutine main_job(head,a8,jp,ma,map,f1,match,rms,seeing,
     *  n1,n2,lhead)
	integer pgbegin
	character*80 head(180)
	integer*2 ma(n1,n2)         	! map buffer
	integer*2 map(n1,n2)
	real*8 a8(8)
	character*12 ac,dc,ac1,dc1		! alpha delta display line
	character*1 ch,old_ch,ch3*3,sg,c3*3	! work char variable
	character*50 ns,dis*60,mns*30		! for display
	byte bell,xyw(99)			! weight of standard posi star
	real x9(99),y9(99),xm(99)		! al_del of starard posi star
	real sx(5),sy(5)			! scale on right_bottom of pic
	real*8 xixn(2,99),xy(2,99)		! x_y of standard posi star
	real*8 adcoef(2,6),xycoef(2,6)		! 2 set 6 plte coef
	character f1*30
	common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut

        call whitexblack(ma,n1/2,n2/2,white,x)
	black=white+20.*x
 	icall=pgbegin(0,'/xw',1,1 )
	call	pgpap(8.5,1.015)
	call pgsci(0)
	x048=n1
	y048=n2
	i9=(n1+n2)/400
              write(*,*)x048,y048
        call pgenv(1.,x048,1.,y048,0.,0)
	call pgsci(1)

  	call pgxgray(ma,n1,n2,1,n1,1,n2,white,black,1)
	call pglabel(' ',' ',f1)

	bell=7
	ir=4
	kexit=0
	down=-i9*8
	call pgsch(.8)
	call pgsci(2)				! standard star color mark
	n=0					! display star serial number
	call unlink('coord.tmp')
9	if(jp.eq.0)write(*,91)bell
91	format(a1,'use left_button to pick standard position star',
     *            ' (another window, run "gsc,gscno" or "dss")')
c put cursor in center of picture
	x=(n1+1)*.5
	y=(n2+1)*.5
	if(jp.eq.1)then
	  ipos=indexpos(head,'A81     ',lhead)
	  read(head(ipos)(11:),*)xycoef(1,1)
	  read(head(ipos+1)(11:),*)xycoef(2,1)
	  read(head(ipos+2)(11:),*)xycoef(1,2)
	  read(head(ipos+3)(11:),*)xycoef(2,2)
	  read(head(ipos+4)(11:),*)xycoef(1,3)
	  read(head(ipos+5)(11:),*)xycoef(2,3)
	  read(head(ipos+6)(21:),*)alpha_c
	  read(head(ipos+7)(21:),*)delta_c
	  a8(1)=xycoef(1,1)
	  a8(2)=xycoef(2,1)
	  a8(3)=xycoef(1,2)
	  a8(4)=xycoef(2,2)
	  a8(5)=xycoef(1,3)
	  a8(6)=xycoef(2,3)
	  a8(7)=alpha_c
	  a8(8)=delta_c
88	  epoch=2000.0
	  write(*,89)bell,epoch
89	  format(1x,a,'output epoch, [',f6.1,'] ',$)
	  read(*,'(a)')ns
	  if(ns(1:1).ne.' ')then
	    read(ns(1:),*,err=88)epoch
	    if(epoch.lt.1900.or.epoch.gt.2100.)goto 88
	  endif
	  call astprs(alpha_c/cx,delta_c/cy,2000.,x,y,epoch)	
	  call toms(x,ac,1)
	  call toms(y,dc,0)
	  write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
	  goto 1001
	endif

10	call pgcurse(x,y,ch)
	if(ch.eq.'A')then
 	  call star(ma,n1,n2,x,y,ir,ierr)
	  if(ierr.ne.0)then
	    write(*,'(2a)')bell,'No star in this region'
	    goto 10
	  endif
	  n=n+1
	  xy(1,n)=x
	  xy(2,n)=513.-y
	  xyw(n)=1
	  write(ch3(1:),'(i3)')n
	  call pgtext(x-i9,y-i9,ch3)
	endif
	if(ch.eq.'D')goto 20			! pick up standard star end.
	goto 10
20	if(n.lt.3)write(*,*)'too few stars to calaulate 6 plate ceof.'
	if(n.lt.3)goto 9
c use "GSC" results, star order of two system must the same
******* read gsc.tmp or dss.tmp begin
21	gd=1                                    ! gd=0 meams read dss
        open(3,file='gsc.tmp',status='old',iostat=ierr)
	if(ierr.ne.0)then
	  gd=0
	  close(3)
          open(3,file='dss.tmp',status='old',iostat=ierr)
	  if(ierr.eq.0)goto 30
	  write(*,*)'I can not find "gsc.tmp" or "dss.tmp"'
	  write(*,*)'when ready, enter "return" please'
	  read(*,'(a)')ch
	  goto 21
	endif
30	if(gd.eq.1)read(3,1,end=40,err=430)i,i1,i2,xa,ch,j1,j2,xd,epoch,x
1	format(i3,3x,i2,1x,i2,1x,f5.2,1x,a,i2,1x,i2,1x,f4.1,2x,f6.1,10x,f5.1)
	if(gd.eq.0)read(3,3,end=40,err=430)i,i1,i2,xa,ch,j1,j2,xd,epoch
3	format(i3,4x,i2,1x,i2,1x,f5.2,1x,a,i2,1x,i2,1x,f4.1,3x,f6.1)
	if(i.gt.n)goto 30
	x9(i)=i1+i2/60.+xa/3600.
	y9(i)=j1+j2/60.+xd/3600.
	if(ch.eq.'-')y9(i)=-y9(i)
	xm(i)=x
	goto 30
430	read(3,'(a)',end=40)ch
	goto 30
c assume center of ccd_frame is fisrt star
40	alpha_c=x9(1)*cx
	delta_c=y9(1)*cy
	close (3)
42	do 41 i=1,n
	xa=x9(i)*cx
	xd=y9(i)*cy
	call standc(alpha_c,delta_c,xa,xd,xi,xn)
	xixn(1,i)=xi
41	xixn(2,i)=xn
*********** read gsc end
c calculate 6 plate coef.
43		write(*,*)'****************'
	call plate(xixn,xy,xyw,n,xycoef,3)
c display 6 coef and ...
	write(*,2)epoch,ac,dc
2	format('Frame center:  (',f6.1,')'/5x,2a)
	write(*,*)(xycoef(1,i),i=1,3)
	write(*,*)(xycoef(2,i),i=1,3)
	sig=0.
	xmaxsig=0.
	m=0.
	do 50 i=1,n
	if(xyw(i).eq.0)goto 50
        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(alpha_c,delta_c,xi,xn,xa,xd)
	xa=xa/cx
	xd=xd/cy
	sigma=((x9(i)-xa)*15.*cos(delta_c))**2+(y9(i)-xd)**2
	sigma=sqrt(sigma)*3600.
	if(xyw(i).eq.1)then
	  sig=sig+sigma
	  if(sigma.gt.xmaxsig)then
	    j=i
	    xmaxsig=sigma
	  endif
	  m=m+1
	endif
	call toms(x9(i),ac1,1)
	call toms(y9(i),dc1,0)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	if(i.eq.1)write(*,49)epoch,epoch
49	format(9x,'GSC (',f6.1,')',12x,'Result (',f6.1,')'
     *  /2x,24('-'),3x,24('-'))
	ch=' '
	if(xyw(i).eq.1)ch='*'
50	write(*,'(i3,2a,2x,2a,g12.2,4x,a)')i,ac1,dc1,ac,dc,sigma,ch
	write(*,'(62x,g12.2)')sig/m
	write(*,'(a,i5,g12.2)')'Max_sigma:',j,xmaxsig
***************
c 3 choice here, after rather good result, must adjest plate center, use 0
60	write(*,*)'1,2....n  delete or recover this star to re_cal 6 coef.'
	write(*,*)' 0        calculate center of frame'
	write(*,*)'-1        quit'
	if(kexit.eq.-1)goto 90
	read(*,*)k
	if(k.eq.-2)goto 90
	if(k.eq.-1)then
	  kexit=-1
	  k=0
	endif
	if(k.eq.0)then
	xa=(n1+1)*0.5
	xd=(n2+1)*0.5
        xi=xycoef(1,1)*xa+xycoef(1,2)*xd+xycoef(1,3)
        xn=xycoef(2,1)*xa+xycoef(2,2)*xd+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,xa,xd)
	alpha_c=xa
	delta_c=xd
        call toms(xa/cx,ac,1)           ! 010329 last 2 line for many tv_?
        call toms(xd/cy,dc,0)
	goto 42
	endif
	if(k.ge.1.and.k.le.n)then
	xyw(k)=-xyw(k)
	goto 43				! re_calculate 6 plate
	endif
	goto 60
******** auto match
90	match=m
	call astprs(alpha_c/cx,delta_c/cy,epoch,xa,xd,2000.)
	epoch=2000.
	alpha_c=xa*cx
	delta_c=xd*cy
	rms=sig/m
********* save a8
	a8(1)=xycoef(1,1)
	a8(2)=xycoef(2,1)
	a8(3)=xycoef(1,2)
	a8(4)=xycoef(2,2)
	a8(5)=xycoef(1,3)
	a8(6)=xycoef(2,3)
	a8(7)=alpha_c
	a8(8)=delta_c			! 2000 epoch
	call toms(alpha_c/cx,ac,1)
	call toms(delta_c/cy,dc,0)
	write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
********* save a8 end
1001	call xytoad(xycoef,adcoef)	! 2 set of 6 coef convert
	xycoef(1,4)=alpha_c
	xycoef(2,4)=delta_c
c draw scale	
	y=abs(xycoef(1,1))
	x=abs(xycoef(1,2))
	if(x.gt.y)y=x
	x=y*180.*60./pi*n1			! X size in minute
	y=4.
111	ii=nint(x/y)
	k=nint(ii/x*n1)
	if(k.lt.n1*0.2.or.ii.lt.4)then
	  y=y*0.9
	  goto 111
	endif
	k4=k/4
	sx(1)=n1-i9*10
	sx(2)=sx(1)
	sx(3)=sx(1)-k
	sx(4)=sx(3)
	sy(1)=i9*10
	sy(2)=sy(1)-i9*3
	sy(3)=sy(2)
	sy(4)=sy(1)
	write(c3(1:),"(i2)")ii
	c3(3:3)=char(39)
	call pgline(4,sx,sy)
	call pgtext(n1-120.,sy(1),c3)
	call pgsci(5)
	call pgtext(n1*0.4,n2+i9*3.,ns)
********************************************* loop
	num=0
	ip=0				! indicator, show stand star again
	open(99,file='coord.tmp',status='unknown')
	ep_old=epoch
	
100	write(*,*)'  left button -- pick out a star'
	write(*,*)'middle button -- put a false_star'
	write(*,*)' right button -- quit'
	write(*,*)'  "n" or "N"  -- input a false_star  (in alpha delta)'
	write(*,*)'  "p" or "P"  -- hard copy'
        write(*,*)'    if "printer.ps" in your defalt directory, then'
        write(*,*)'        produces "pgplot.jp" with postscript style'
	write(*,*)'  "c" or "C"  -- Coordination'
	write(*,*)'  "z" or "Z"  -- Zoom'

	if(ip.eq.0)then
c display standard star again in display line position
c another purpose is for serial number
	do 101 i=1,n
	x=xy(1,i)
	y=xy(2,i)
	call star(ma,n1,n2,x,y,ir,ierr)
	if(ierr.ne.0)goto 101
        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(alpha_c,delta_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
	num=num+1
	xa=amod(xa,24.)
	call toms(xa,ac,1)
	call toms(xd,dc,0)
	sg=' '
	if(xyw(i).eq.1)sg='*'
c in display, * stand for standard star, ? stand for saturated star
	write(*,1335)num,sg,ac,dc,epoch
1335	format(i3,')',2a,1x,a,'(',f6.1,')')
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),1335)num,sg,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
101	continue
	call pgsci(13)
	call pgebuf
	endif
	ip=1			! only do once, display at "help" tail
	x=n1/2+0.5
	y=n2/2+0.5
	izz=0
	ixx=0
	iyy=0
***** loop in here many time,
110	old_ch=ch		! for zoom, save last command.
	x=x-ixx
	y=y-iyy
	call pgcurse(x,y,ch)

****** ZOOM
	if(ch.eq.'z'.or.ch.eq.'Z')then
	  if(izz.eq.1)then
	    x=x*4.
	    y=y*4.
	  endif
	  ixx=n1+1-x
	  iyy=n2+1-y
	  if(ixx.lt.256)ixx=256
	  if(iyy.lt.256)iyy=256
	  k=n1-256
	  if(ixx.gt.k)ixx=k
	  if(iyy.gt.k)iyy=k
	  x=n1+1-ixx
	  y=n2+1-iyy
	  if(izz.eq.0)then
	    call pgsci(3)
	    sx(1)=x-256
	    sx(2)=x+256
	    sx(3)=x+256
	    sx(4)=x-256
	    sy(1)=y-256
	    sy(2)=y-256
	    sy(3)=y+256
	    sy(4)=y+256
	    sx(5)=sx(1)
	    sy(5)=sy(1)
	    call pgline(5,sx,sy)
	  endif    
	  ii=0
	  do 25364 i=ixx-255,ixx+256
	  ii=ii+1
	  jj=0
	  do 25364 j=iyy-255,iyy+256
	  jj=jj+1
25364	  map(ii,jj)=ma(i,j)
	  m1=512
	  m2=512

          call whitexblack(map,512,512,white,z)
	    write(*,*)white,z
	    if(white.le.1.)then
		white=whiteold
	        z=zold
	    endif
	    whiteold=white
	    zold=z
	  black=white+20.*z
	  call pgsch(1.)
          call pgsci(0)
          call pgenv(1.,513.,-1.,511.,0,0)
          call pgsci(1)
          call pgximage(map,512,512,white-5*z,black)
          call pglabel('R.A. south','Decl. east',f1)
	  call pgsch(2.)
	  call pgpoint(1,256.5,256.5,5)
	  call pgsch(.8)
	  
           
        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(alpha_c,delta_c,xi,xn,x_a,x_d)
        xa=x_a/cx
        xd=x_d/cy
        xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
        call toms(x_a,ac,1)
        call toms(x_d,dc,0)
        write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
	call pgsci(5)
        call pgtext(0.,522.,ns)
	write(*,*)'-------'
	izz=1
	ixx=x-256.
	iyy=y-256.
	down=m2+20.
	ir=3
	endif

c left button****************
	if(ch.eq.'A')then
	call pgsci(6)
	if(izz.eq.1)then
	  x=x+ixx
	  y=y+iyy
	endif
	oldx=x
	oldy=y
	call star(ma,n1,n2,x,y,ir,ierr)	! find star
	call star(ma,n1,n2,x,y,ir,ierr)	! find star
	if(ierr.eq.1.and.izz.eq.0)goto 200
	if(izz.eq.1)then
	  if((oldx-x)**2+(oldy-y)**2.gt.6)then
	    x=oldx
	    y=oldy
	    call starpeak(ma,n1,n2,x,y,white)
	    call starpeak(ma,n1,n2,x,y,white)
	  endif
	endif   
	xi=xycoef(1,1)*x+xycoef(1,2)*(513.-y)+xycoef(1,3)
        xn=xycoef(2,1)*x+xycoef(2,2)*(513.-y)+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
c for display ony
	num=num+1
	write(ch3(1:),'(i3)')num
	if(izz.eq.0)then
	call pgtext(x-i9,y-i9,ch3)
	else
	call pgtext(x-ixx-6.,y-iyy-6.,ch3)
	endif
	xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),1334)num,ac,dc,epoch
1334	format(i3,') ',a,1x,a,'(',f6.1,')')
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1,i8)')dis(1:50),x,y,ma(nint(x),nint(y))
	goto 110	
	endif

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

c middle  button**************
	if(ch.eq.'D')then
	call pgsci(4)
	if(izz.eq.1)then
	  x=x+ixx
	  y=y+iyy
	endif
        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(alpha_c,delta_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
	Xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	num=num+1
	write(ch3(1:),'(i3)')num
	if(izz.eq.0)then
	 call pgpoint(1,x,y,16)
	 call pgtext(x-i9,y-i9,ch3)
	else
	 call pgpoint(1,x-ixx+2.,y-iyy,16)
	 call pgtext(x-ixx-6.,y-iyy-6,ch3)
	endif
104	format(i3,') ',a,1x,a,'(',f6.1,') false star')
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),104)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	goto 110
	endif

c file  star******************
	if(ch.eq.'f'.or.ch.eq.'F')then
	ep_new=2000
	write(*,'(a)')' Input file name (epoch=2000): '
	read(*,'(a)')mns
	open(98,file=mns,status='old')
9120	read(98,'(a)',end=91200)mns
	write(*,*)mns
        ii=1
	if(mns(1:1).eq.' ')ii=2 
	if(itohd(mns(ii:),x_a,x_d,x).eq.0)i=i+1
	call astprs(x_a,x_d,ep_new,xa,xd,2000.)
	x_a=xa*cx
	x_d=xd*cy
	call standc(alpha_c,delta_c,x_a,x_d,xi,xn)
	x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
	y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	if(x.lt.0.or.y.lt.0.or.x.gt.n1.or.y.gt.n2)then
	  write(*,'(2a)')bell,'out of the frame'
	  goto 9120
	endif
	if(izz.eq.1)then
	  xx=x-ixx
	  yy=y-iyy
	  if(xx.lt.0.or.yy.lt.0.or.xx.gt.m1.or.yy.gt.m2)then
	  write(*,'(2a)')bell,'out of the frame'
	  goto 9120
	endif
	endif
	xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgsci(2)
	if(izz.eq.0)then
	  call pgpoint(1,x,y,-1)
	  call pgtext(x-i9,y-i9,ch3)
	else
	 call pgpoint(1,xx+2,yy,-1)
	 call pgtext(xx-6.,yy-6,ch3)
	endif
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),105)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	goto 9120
91200	close(98)
	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)')mns
	if(mns(1:1).eq.' ')goto 110
	if(itohd(mns,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,2000.)
	x_a=xa*cx
	x_d=xd*cy
	call standc(alpha_c,delta_c,x_a,x_d,xi,xn)
	x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
	y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	if(x.lt.0.or.y.lt.0.or.x.gt.n1.or.y.gt.n2)then
	  write(*,'(2a)')bell,'out of the frame'
	  goto 120
	endif
	if(izz.eq.1)then
	  xx=x-ixx
	  yy=y-iyy
	  if(xx.lt.0.or.yy.lt.0.or.xx.gt.m1.or.yy.gt.m2)then
	  write(*,'(2a)')bell,'out of the frame'
	  goto 120
	endif
	endif
	xa=amod(xa,24.)
	call astprs(xa,xd,2000.,x_a,x_d,epoch)
	call toms(x_a,ac,1)
	call toms(x_d,dc,0)
	num=num+1
	write(ch3(1:),'(i3)')num
	call pgsci(2)
	if(izz.eq.0)then
	  call pgpoint(1,x,y,18)
	  call pgtext(x-i9,y-i9,ch3)
	else
	 call pgpoint(1,xx+2,yy,18)
	 call pgtext(xx-6.,yy-6,ch3)
	endif
105	format(i3,') ',a,1x,a,'(',f6.1,') input alpha delta')
	call pgsci(0)
	call pgtext(20.,down,dis)
	write(dis(1:),105)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(20.,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	write(*,"(' Again ? [Y] ',$)")
	read(*,'(a)')sg
	if(sg.ne.'n'.and.sg.ne.'N')goto 120	
	endif


c display coordination
	if(ch.eq.'c'.or.ch.eq.'C')then
	  alpha_save=alpha_c
	  delta_save=delta_c
	  call astprs(alpha_save/cx,delta_save/cy,2000.,
     *              alpha_c,delta_c,epoch)
	  alpha_c=alpha_c*cx
	  delta_c=delta_c*cy
	xycoef(1,4)=alpha_c
	xycoef(2,4)=delta_c
	  x=delta_c/cy
	  x60=60.
	  if(x.gt.60.)x60=30.
	  if(x.gt.80.)x60=12.
	  if(x.gt.86.)x60=6.
	  if(x.gt.88.)x60=3.
	  if(x.gt.89.)x60=2.
          call xy_rad(1.,1.,xycoef,a1,d1)
          call xy_rad(1.,2048.,xycoef,a2,d2)
          call xy_rad(2048.,1.,xycoef,a3,d3)
          call xy_rad(2048.,2048.,xycoef,a4,d4)
	  x=amin1(a1,a2,a3,a4)
	  y=amax1(a1,a2,a3,a4)
	  i1=x*x60
	  i2=y*x60+1
	  x=amin1(d1,d2,d3,d4)
	  y=amax1(d1,d2,d3,d4)
	  j1=x*6.
	  if(x.lt.0)j1=j1-1
	  j2=y*6.
	  if(y.gt.0)j2=j2+1
	  do 183 i=i1,i2
	  xa=i/x60*cx
	  k=0
	  do 182 j=j1,j2
	  k=k+1
	  xd=j/6.*cy
	  call standc(alpha_c,delta_c,xa,xd,xi,xn)
	  x9(k)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
182	  y9(k)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
183	  call pgline(k,x9,y9)
	  do 184 j=j1,j2
	  xd=j/6.*cy
	  k=0
	  do 185 i=i1,i2
	  k=k+1
	  xa=i/x60*cx
	  call standc(alpha_c,delta_c,xa,xd,xi,xn)
	  x9(k)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
185	  y9(k)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
184	  call pgline(k,x9,y9)
c draw AZIMUTH & ALTITUDE grid, if to near pole, not draw
189	  j=(iy+((im-1)*30.4+id)/365.)*10.+0.5
	  epoch_g=j*0.1
          call astprs(a1,d1,epoch,a1,d1,epoch_g)
          call astprs(a2,d2,epoch,a2,d2,epoch_g)
          call astprs(a3,d3,epoch,a3,d3,epoch_g)
          call astprs(a4,d4,epoch,a4,d4,epoch_g)
	  call astprs(alpha_c/cx,delta_c/cy,epoch,alpha_g,delta_g,epoch_g)
	  call ad_azl(a1,d1,az1,al1)
          call ad_azl(a2,d2,az2,al2)
          call ad_azl(a3,d3,az3,al3)
          call ad_azl(a4,d4,az4,al4)
	  call ad_azl(alpha_g,delta_g,a1,d1)
	  alpha_g=alpha_g*cx
	  delta_g=delta_g*cy
	  x60=6.
	  if(d1.gt.30.)x60=4.
	  if(d1.gt.50.)x60=3.
	  if(d1.gt.70.)x60=2.
	  if(d1.gt.80.)x60=1.
	  if(d1.gt.84.)x60=0.5
	  if(d1.gt.86.)x60=0.4
	  if(d1.gt.88.)x60=0.2
	  if(d1.gt.89.)x60=0.1
	  call pgsci(3)
	  call pgsls(3)
	  x=amin1(az1,az2,az3,az4)
	  y=amax1(az1,az2,az3,az4)
          if(y-x.gt.200.)then
            xa=x
            x=y-360-5./x60
            y=xa+5./x60
          endif
	  i1=x*x60
	  i2=y*x60+1
	  x=amin1(al1,al2,al3,al4)
	  y=amax1(al1,al2,al3,al4)
	  j1=x*6.
	  j2=y*6.+1
	  kk=0
	  do 193 i=i1,i2
	  xa=i/x60
	  k=0
	  do 192 j=j1,j2
	  k=k+1
	  xd=j/6.
	  call azl_ad(xa,xd,x,y)
	  call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	  x9(k)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
192	  y9(k)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
193	  call pgline(k,x9,y9)
	  do 194 j=j1,j2
	  xd=j/6.
	  k=0
	  do 195 i=i1,i2
	  k=k+1
	  xa=i/x60
	  call azl_ad(xa,xd,x,y)
	  call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	  x9(k)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
195	  y9(k)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
194	  call pgline(k,x9,y9)
	  call pgsls(1)
	  alpha_c=alpha_save
	  delta_c=delta_save
	xycoef(1,4)=alpha_c
	xycoef(2,4)=delta_c
	endif

c hardcopy ****************
	if(ch.eq.'p'.or.ch.eq.'P')then
	  write(*,'(a)')bell
	  call pgend
	  close (99)
	  if(izz.eq.1)goto 130
	
	  alpha_save=alpha_c
	  delta_save=delta_c
	  call astprs(alpha_save/cx,delta_save/cy,2000.,
     *              alpha_c,delta_c,epoch)
	  alpha_c=alpha_c*cx
	  delta_c=delta_c*cy
	xycoef(1,4)=alpha_c
	xycoef(2,4)=delta_c
	  open(99,file='coord.tmp',status='old',iostat=i)
	  call unlink('pgplot.hp')
          icall=pgbegin(0,'/hp',1,1)
          call pgslw(2)
          write(*,*)'wait for 1 minute.'
c physical size of laser_printer, 300 dots/inch, 2336 dots in width
c                                 3200 dots in length
	  xp=17.06
          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)
          call pgwindow(1.,512.,1.,512.)
          call pgbox('bc',0,0,'bc',0,0)
  	  call pgxgray(ma,512,512,1,512,1,512,white,black,2)
	  call pgsch(0.8)
          call pglabel('R.A  South','Decl.  East',f1//ns)
	  sx(1)=500
	  sx(2)=sx(1)
	  sx(3)=sx(1)-k4
	  sx(4)=sx(3)
	  sy(1)=15
	  sy(2)=10
	  sy(3)=sy(2)
	  sy(4)=sy(1)
	  call pgline(4,sx,sy)
	  call pgtext(450.,15.,c3)
          y1=-4.*dy                       ! for print display line on paper
          x1=(1.-17./xp)*256.             ! in 2 coloumn, cal x position first
          call pgsch(.6)
          i=0
150       i=i+1
          read(99,'(a,2f7.1)',end=151)dis,x,y
	  x=x*.25
	  y=(n2+1-y)*0.25
          call pgtext(x1,y1,dis)
          y1=y1-dy
c mark false star
          if(dis(40:40).eq.'f'.or.dis(40:40).eq.'i')then
            call pgslw(1)
            call pgpoint(1,x,y,23)
            call pgslw(2)
          endif
c write No. on picture
          call pgtext(x-3.,y+4.,dis(1:3))
c see if second column for writting
          if(i.eq.10)then
            x1=260.
            y1=-4.*dy
          endif
          goto 150
151       close (99)
c draw grid, if to near pole, not draw
          x=delta_c/cy
          x60=60.
          if(x.gt.60.)x60=30.
          if(x.gt.80.)x60=12.
          if(x.gt.86.)x60=6.
          if(x.gt.88.)x60=3.
          if(x.gt.89.)x60=2.
	  call pgslw(1)
          call xy_rad(1.,1.,xycoef,a1,d1)
          call xy_rad(1.,2048.,xycoef,a2,d2)
          call xy_rad(2048.,1.,xycoef,a3,d3)
          call xy_rad(2048.,2048.,xycoef,a4,d4)
	  x=amin1(a1,a2,a3,a4)
	  y=amax1(a1,a2,a3,a4)
	  i1=x*x60
	  i2=y*x60+1
	  x=amin1(d1,d2,d3,d4)
	  y=amax1(d1,d2,d3,d4)
	  j1=x*6.
	  if(x.lt.0)j1=j1-1
	  j2=y*6.
	  if(y.gt.0)j2=j2+1
	  do 153 i=i1,i2
	  xa=i/x60*cx
	  k=0
	  do 152 j=j1,j2
	  k=k+1
	  xd=j/6.*cy
	  call standc(alpha_c,delta_c,xa,xd,xi,xn)
	  x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
	  y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	  x9(k)=x*0.25
152	  y9(k)=(n2+1-y)*0.25
	  call pgline(k,x9,y9)
	  call toms(i/x60,ac,1)
	  ns=ac(2:3)//'\\uh\\d'//ac(5:6)//'\\um\\d'
153	  if(x9(1).gt.0..and.x9(1).lt.512.)call pgtext(x9(1)-15.,-9.,ns)
	  do 154 j=j1,j2
	  xd=j/6.*cy
	  k=0
	  do 155 i=i1,i2
	  k=k+1
	  xa=i/x60*cx
	  call standc(alpha_c,delta_c,xa,xd,xi,xn)
	  x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
	  y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	  x9(k)=x*0.25
155	  y9(k)=(n2+1-y)*0.25
	  call pgline(k,x9,y9)
	  call toms(j/6.,dc,0)
	  ns=dc(1:3)//'\\uo\\d'//dc(5:6)//char(39)
154	  if(y9(k).gt.0..and.y9(k).lt.512.)call pgtext(-33.,y9(k),ns)
c draw AZIMUTH & ALTITUDE grid, if to near pole, not draw
159       j=(iy+((im-1)*30.4+id)/365.)*10.+0.5
          epoch_g=j*0.1
          call astprs(a1,d1,epoch,a1,d1,epoch_g)
          call astprs(a2,d2,epoch,a2,d2,epoch_g)
          call astprs(a3,d3,epoch,a3,d3,epoch_g)
          call astprs(a4,d4,epoch,a4,d4,epoch_g)
	  call astprs(alpha_c/cx,delta_c/cy,epoch,alpha_g,delta_g,epoch_g)
          call ad_azl(a1,d1,az1,al1)
          call ad_azl(a2,d2,az2,al2)
          call ad_azl(a3,d3,az3,al3)
          call ad_azl(a4,d4,az4,al4)
          call ad_azl(alpha_g,delta_g,a1,d1)
          alpha_g=alpha_g*cx
          delta_g=delta_g*cy
          x60=6.
          if(d1.gt.30.)x60=4.
          if(d1.gt.50.)x60=3.
          if(d1.gt.70.)x60=2.
          if(d1.gt.80.)x60=1.
	  if(d1.gt.84.)x60=0.5
          if(d1.gt.86.)x60=0.4
          if(d1.gt.88.)x60=0.2
          if(d1.gt.89.)x60=0.1
	  call pgsci(3)
	  call pgsls(3)
c judge direction of A,H
c center (1024.5,1024.5)--a1,l1
c another point, a2=a1, l2=l1-1
	  call azl_ad(a1,d1-1.,x,y)
	  call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	  x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
	  y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	  a3=x-1024.5
	  d3=y-1024.5
	  x=a3/sqrt(a3*a3+d3*d3)
	  x=asin(x)/cy
	  if(x.lt.0.)x=-x
	  iah=0
	  if(x.gt.48.)iah=1
c iah=0:     x direction is (A); y direction is (H)
c iah=1:     x direction is (H); y direction is (A)
	  x=amin1(az1,az2,az3,az4)
	  y=amax1(az1,az2,az3,az4)
          if(y-x.gt.200.)then
            xa=x
            x=y-360-5./x60
            y=xa+5./x60
          endif
	  i1=x*x60
	  i2=y*x60+1
	  x=amin1(al1,al2,al3,al4)
	  y=amax1(al1,al2,al3,al4)
	  j1=x*6.
	  j2=y*6.+1
	  kk=0
	  do 163 i=i1,i2
	  xa=i/x60
	  k=0
	  do 162 j=j1,j2
	  k=k+1
	  xd=j/6.
	  call azl_ad(xa,xd,x,y)
	  call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	  x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
	  y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	  x9(k)=x*0.25
	  y9(k)=(n2+1-y)*0.25
162	continue
	  call pgline(k,x9,y9)
	  xa=xa+0.0001
	  j=xa
	  k=(xa-j)*60.
	  write(ac(1:),'(i3,1x,i2.2)')j,k
	  ns=ac(1:3)//'\\uo\\d'//ac(5:6)//char(39)
	  if(iah.eq.0)then
	    x=x9(1)-(x9(1)-x9(2))*(y9(1)-512.)/(y9(1)-y9(2))
	    if(x.gt.0..and.x.lt.512.)then
	      call pgtext(x-15.,516.,ns)
	    endif
	  else
	    y=y9(1)-(y9(1)-y9(2))*(x9(1)-512.)/(x9(1)-x9(2))
	    if(y.gt.0.and.y.lt.512)call pgtext(512.,y-3.,ns)
  	    if(y.gt.200.and.y.lt.350.and.kk.eq.0)then
	      call pgtext(520.,y+40.,'(A)')
	      kk=1
	    endif
	  endif
163	  continue
	   xd=525.
	  if(iah.eq.1)call pgtext(490.,xd,'(H)')
	  if(iah.eq.0)call pgtext(490.,xd,'(A)')
	  do 164 j=j1,j2
	  xd=j/6.
	  k=0
	  do 165 i=i1,i2
	  k=k+1
	  xa=i/x60
	  call azl_ad(xa,xd,x,y)
	  call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn)
	  x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
	  y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
	  x9(k)=x*0.25
165	  y9(k)=(n2+1-y)*0.25
	  call pgline(k,x9,y9)
	  call toms(xd,dc,0)
	  ns=dc(1:3)//'\\uo\\d'//dc(5:6)//char(39)
	  if(iah.eq.0)then
	    if(x9(1).gt.x9(2))then
	      do 166 i=1,k
	      if(x9(i).lt.512.)goto 167
166	      continue
	    else
	      do 168 i=1,k
	      if(x9(i).gt.512.)goto 167
168	      continue
	    endif
167	    if(i.gt.k)goto 164
	    y=y9(i)-(y9(i)-y9(i-1))*(x9(i)-512.)/(x9(i)-x9(i-1))
	    if(y.gt.0..and.y.lt.512.)then
	      call pgtext(512.,y-3.,ns)
  	      if(y.gt.200.and.y.lt.350.and.kk.eq.0)then
	        call pgtext(520.,y+50.,'(H)')
	        kk=1
	      endif
	    endif
	  else
	    if(y9(1).gt.y9(2))then
	      do 266 i=1,k
	      if(y9(i).lt.0.)goto 267
266	      continue
	    else
	      do 268 i=1,k
	      if(y9(i).gt.0.)goto 267
268	      continue
	    endif
267	    if(i.gt.k)goto 164
	    x=x9(i)-(x9(i)-x9(i-1))*(y9(i)-512.)/(y9(i)-y9(i-1))
	    if(x.gt.0..and.x.lt.512.)then
	      call pgtext(x-15.,516.,ns)
	    endif
	  endif
164	  continue
c/* draw somthing about moon */
 	  call getmxy(a8,x,y,z,phase,pl,pb,xx)
	  xa=2100.
	  xd=250.
	  a1=xa+60.
	  d1=xd+60.
	  a2=xa+60.
	  d2=xd-60.
	  a3=xa-60.
	  d3=xd-60.
	  a4=xa-60.
	  d4=xd+60.
	  call grlin2(a1,d1,a2,d2)
	  call grlin2(a2,d2,a3,d3)
	  call grlin2(a3,d3,a4,d4)
	  call grlin2(a4,d4,a1,d1)
	  if(pb.gt.0)then
	    a1=x*150.+xa
	    d1=y*150.+xd
	    call grlin2(a1,d1,xa,xd)
	    dis='moon(??.?) to ccd '
	    write(dis(6:9),'(f4.1)')phase
	    i=xx+0.5
	    write(dis(19:),'(i3)')i
	    dis=dis(1:21)//'\\uo'
	    call pgtext(440.,-230.,dis)
	  else
	    call pgtext(500.,-230.,"no moon")
	  endif
          call pgiden
          call pgend
c submit "print map job" to system, in queeue, print map automatically.
          call system('print pgplot.hp')
	  write(*,'(a)')bell
	  alpha_c=alpha_save
	  delta_c=delta_save
	xycoef(1,4)=alpha_c
	xycoef(2,4)=delta_c
	  return
	endif
	
c *********
200	if(ch.lt.'A')write(*,'(a)')bell
	goto 110

130	  open(99,file='coord.tmp',status='old',iostat=i)
	  call unlink('pgplot.hp')
          icall=pgbegin(0,'/hp',1,1)
          call pgslw(2)
          write(*,*)'wait for 1 minute.'
	  xp=17.06
          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)
          call pgwindow(1.,512.,1.,512.)
          call pgbox('bc',0,0,'bc',0,0)
  	  call pgxgray(ma,512,512,1,512,1,512,white,black,2)
	  call pgsch(0.8)
          call pglabel('R.A  South','Decl.  East',f1//ns)
	  sx(1)=500
	  sx(2)=sx(1)
	  sx(3)=sx(1)-141
	  sx(4)=sx(3)
	  sy(1)=15
	  sy(2)=10
	  sy(3)=sy(2)
	  sy(4)=sy(1)
	  call pgline(4,sx,sy)
	  c3(1:2)=' 4'
	  call pgtext(450.,15.,c3)
          y1=-4.*dy                       ! for print display line on paper
          x1=(1.-17./xp)*256.             ! in 2 coloumn, cal x position first
          call pgsch(.6)
          i=0
1150       read(99,'(a,2f7.1)',end=1151)dis,x,y
	  x=x-ixx
	  y=513.-(y-iyy)
	  if(x.lt.1.or.x.gt.511.or.y.lt.0.or.y.gt.511)goto 1150
          i=i+1
          call pgtext(x1,y1,dis)
          y1=y1-dy
c mark false star
          if(dis(40:40).eq.'f'.or.dis(40:40).eq.'i')then
            call pgslw(1)
            call pgpoint(1,x,y,23)
            call pgslw(2)
          endif
c write No. on picture
          call pgtext(x-3.,y+4.,dis(1:3))
c see if second column for writting
          if(i.eq.10)then
            x1=260.
            y1=-4.*dy
          endif
          goto 1150
1151       close (99)
	
        call pgiden
        call pgend
c print map job to system
        call system('print pgplot.hp')
        write(*,'(1x,a)')bell
        return

	end
	
	subroutine starpeak(map,n1,n2,x,y,sky)
	integer*2 map(n1,n2)
	real b(-2:2)
	kx=nint(x)
	ky=nint(y)
	iz=0
	do 10 i=kx-2,kx+2
	do 10 j=ky-2,ky+2
	if(map(i,j).gt.iz)then
	  iz=map(i,j)
	  ix=i
	  iy=j
	endif
	i5=2
10	continue
cccccc x center    weight_center method
        do 30 i=ix-i5,ix+i5
        z=0.
        do 40 j=iy-i5,iy+i5
40      z=z+map(i,j)-sky
30      b(i-ix)=z
        z=0.
        s=0.
        do 50 i=-i5,i5
        z=z+b(i)
50      s=s+b(i)*i
	x=ix
        if(z.ne.0.) x=x+s/z
cccccc y center
        do 130 j=iy-i5,iy+i5
        z=0.
        do 140 i=ix-i5,ix+i5
140     z=z+map(i,j)-sky
130     b(j-iy)=z
        z=0.
        s=0.
        do 150 i=-i5,i5
        z=z+b(i)
150     s=s+b(i)*i
	y=iy
        if(z.ne.0.) y=y+s/z
	end

	subroutine star(map,n1,n2,xx,yy,ir,ierr)
c when ir=4, 15*15 for center
	integer*2 map(n1,n2),maps(-25:25,-25:25)
	real b(-25:25),xmap(-25:25,-25:25)
c sky first
	n25=25
	x=xx
	y=yy
	kx=nint(x)
	ky=nint(y)
	if(kx.le.n25)kx=n25+1
	if(kx.gt.n1-n25)kx=n1-n25
	if(ky.le.n25)ky=n25+1
	if(ky.gt.n2-n25)ky=n2-n25
c	  write(*,*)'xy: ',kx,ky
c	  do 5 j=ky-4,ky+4
c5	  write(*,'(9i6)')(map(i,j),i=kx+4,kx-4,-1)
c	  write(*,*)
	do 10 j=-n25,n25
	do 10 i=-n25,n25
10	maps(i,j)=map(kx+i,ky+j)	
	n51=n25+n25+1
	call whitexblack(maps,n51,n51,sky,sigma)
	ierr=1					! error =1
	i5=ir*1.1+2
	do 1000 ll=1,7
	 if(x.lt.1..or.x.gt.n1)return
	 if(y.lt.1..or.y.gt.n2)return
	ix=nint(x)
	iy=nint(y)
cccccc x center    weight_center method
	do 30 i=ix-i5,ix+i5
	z=0.
	do 40 j=iy-i5,iy+i5
40	z=z+map(i,j)-sky		! must minus sky, other wise wrong
30	b(i-ix)=z			! for weight_center cal.
	z=0.
	s=0.
	do 50 i=-i5,i5
	z=z+b(i)
50	s=s+b(i)*i
	s=s/z
	x=ix+s
cccccc y center
	do 130 j=iy-i5,iy+i5
	z=0.
	do 140 i=ix-i5,ix+i5
140	z=z+map(i,j)-sky
130	b(j-iy)=z
	z=0.
	s=0.
	do 150 i=-i5,i5
	z=z+b(i)
150	s=s+b(i)*i
	s=s/z
	y=iy+s
1000	continue
	ix=nint(x)
	iy=nint(y)
	do 101 j=-n25,n25
	do 101 i=-n25,n25
101	xmap(i,j)=map(ix+i,iy+j)	
	x=ix-n25+center_1(xmap,sky,51,1)
	y=iy-n25+center_1(xmap,sky,51,0)
	 if(x.lt.1..or.x.gt.n1)return
	 if(y.lt.1..or.y.gt.n2)return
	ierr=0
          xx=x
	  yy=y
	return
	end

        function center_1(a,sky,n,ip)
        real a(51,51)
        integer b(51)
        skyn=sky*n
        if(ip.eq.1)then
        do 10 i=1,n
        x=0.
        do 20 j=1,n
20      x=x+a(i,j)
10      b(i)=x-skyn
        else
        do 110 i=1,n
        x=0.
        do 120 j=1,n
120     x=x+a(j,i)
110     b(i)=x-skyn
        endif
        call histat(b,mode,maxh,xmean,xpeak,fsigma,gsigma,n)
        center_1=xpeak
        if(abs(xpeak-25.).gt.3.)center_1=25.
        end

        subroutine xy_rad(x,y,a8,alpha,delta)
c 2048*2048 ccd's X Y (in pixel) -------------> alpha delta (in h,d)
	real x,y,alpha,delta
	real*8 a8(8)				! from fits_head "A81"--"A88"
	pi=4.*atan(1.)
	xi=a8(1)*x+a8(3)*y+a8(5)
	xn=a8(2)*x+a8(4)*y+a8(6)
	rc=a8(7)				! center of ccd (in rad.)
	dc=a8(8)
        cd=cos(dc)
        td=tan(dc)
        alpha=atan(xi/cd/(1.-xn*td))
        delta=atan((xn+td)*cos(alpha)/(1.-xn*td))
        alpha=alpha+rc
        if(abs(dc-delta).gt.2.)then
          delta=-delta
          alpha=alpha+pi
        endif
	alpha=alpha*12./pi
	delta=delta*180./pi
        end

	function indexpos(head,f1,lhead)
	character*80 head(180),f1*8
	byte bell
	data bell/7/
	do 10 indexpos=1,lhead-1
10	if(head(indexpos)(1:8).eq.f1)return
        if(indexpos.eq.lhead)write(*,'(1x,a)')bell
	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

        function xmedian(a,n)
        real a(n)
        do 10 i=1,n-1
        do 10 j=i+1,n
        if(a(i).lt.a(j))goto 10
        x=a(i)   
        a(i)=a(j)
        a(j)=x  
10      continue
        m=n/2
        xmedian=a(m+1)
        if(m+m.ne.n)return
        xmedian=(a(m)+a(m+1))/2
        end
