character*40 f1,f2*5,f3*5
 	integer*2 a(2048*2048*2)
	character*80 head(72)
	real*8 a8(8)
	k=iargc()
        if(k.lt.3)stop 'usage:  sel file.fit ixc iyc'
        call getarg(1,f1)
	call getarg(2,f2)
	read(f2(1:),*)ixc
	call getarg(3,f3)
	read(f3(1:),*)iyc
        if(inqh(f1,n1,n2,head,ibit).eq.0)stop 'DATA not found'
	iturn=1
        ipix=211
c read data
10	call readata(a,f1,ibit)
	k=2048*2048
	do 20 i=1,k
20	a(k+i)=a(i)
	m1=2048
	m2=2048
        if(iturn.eq.0)call turnx2(a(k+1),m1,m2,-180)
        if(iturn.eq.1)call turnx2(a(k+1),m1,m2,-360)
        if(iturn.eq.4)call turnx2(a(k+1),m1,m2, 180)
	f1=f1(1:lnblnk(f1))//'  '//f2(1:lnblnk(f2))//'  '//f3
	  write(*,*)f1
	call main_job(head,a,a(k+1),a8,f1,iturn,ipix,ixc,iyc)
	end

	subroutine main_job(head,map,pam,a8,f1,it,ipix,ixc,iyc)
	integer pgbegin
	character*80 head(72)
	integer*2 map(2048,2048)		! map buffer
	integer*2 pam(2048,2048)
	integer*2 ma(512*512)
	real*8 a8(8)
	character*12 ac,dc			! alpha delta display line
	character*1 ch,ch3*3,sg,c3*3		! work char variable
	character*50 ns,dis*60,mns*30		! for display
	byte bell				! weight of standard posi star
	real sx(5),sy(5)			! scale on right_bottom of pic
	real*8 dummy8				! x_y of standard posi star
	real*8 adcoef(2,6),xycoef(2,6)		! 2 set 6 plte coef
	character f1*40

	pi=4.*atan(1.0)
	cx=pi/12.
	cy=pi/180.

	iu1=ixc-ipix                              ! 211
	iu2=ixc+ipix
	iv1=iyc-ipix
	iv2=iyc+ipix
	u1=iu1
	u2=iu2
	v1=iv1
	v2=iv2
	xc=ixc
	yc=iyc

	m=0
	do 10 i=iu1,iu2
	do 10 j=iv1,iv2
	k=map(i,j)
	if(k.le.0.or.k.gt.30000)goto 10
	m=m+1
	ma(m)=k
10	continue
        call whitexblack(ma,1,m,white,x)
	white=white-5
	black=white+20.*x

 	icall=pgbegin(0,'/xw',1,1 )
        if(it.eq.1)then
	call pgenv(u2,u1,v1,v2,0,0)
  	call pgxgray(pam,2048,2048,2049-iu2,2049-iu1,iv1,iv2,white,black,1)
	endif
	call pglabel('R.A. south','Decl. east',f1)

	bell=7
	ir=4
	down=v1-28.
	n1=2048
	n2=2048

	call pgsch(.8)
	n=0					! display star serial number
	call unlink('coord.tmp')
c put cursor in center of picture
	x=xc
	y=yc
	  ipos=indexpos(head,'A81     ')
	  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
	xi=xycoef(1,1)*(2049-x)+xycoef(1,2)*(2049-y)+xycoef(1,3)
        xn=xycoef(2,1)*(2049-x)+xycoef(2,2)*(2049-y)+xycoef(2,3)
        call astand(alpha_c,delta_c,xi,xn,x_a,x_d)
	xa=x_a/cx
	xd=x_d/cy
	  call astprs(xa,xd,2000.,x,y,epoch)	
	  call toms(x,ac,1)
	  call toms(y,dc,0)
	  write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch

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))+abs(xycoef(2,2)))*.5
	x=y*180.*60./pi*ipix*2			! X size in minute
	ii=nint(x/4.)
	k=nint(ii/x*ipix*2)
	k4=k/4
	sx(1)=u2-20.
	sx(2)=sx(1)
	sx(3)=sx(1)-k
	sx(4)=sx(3)
	sy(1)=v1+10.
	sy(2)=sy(1)-3.
	sy(3)=sy(2)
	sy(4)=sy(1)
	write(c3(1:),"(i2)")ii
	c3(3:3)=char(39)
	call pgsci(2)
	call pgline(4,sx,sy,1)
	call pgtext(u2-50.,sy(1),c3)
	call pgsci(3)
	call pgtext(u2-20.,v2+10.,ns)
********************************************* loop
	num=0
	open(99,file='coord.tmp',status='unknown')
	ep_old=epoch
	
	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(*,*)'  "c" or "C"  -- Coordination'

	x=xc
	y=yx
***** loop in here many time,
110	call pgcurse(x,y,ch)
        ich=ichar(ch)		

c left button****************
	if(ich.eq.128)then
        oldx=x
        oldy=y
	call star(map,n1,n2,x,y,ir,ierr)	! find star
	call star(map,n1,n2,x,y,ir,ierr)	! find star
	if(ierr.eq.1)goto 200
        if((oldx-x)**2+(oldy-y)**2.gt.6)then
          x=oldx
          y=oldy
          call starpeak(map,n1,n2,x,y,white)
          call starpeak(map,n1,n2,x,y,white)
        endif
	xi=xycoef(1,1)*(2049-x)+xycoef(1,2)*(2049-y)+xycoef(1,3)
        xn=xycoef(2,1)*(2049-x)+xycoef(2,2)*(2049-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
	call pgsci(6)
	call pgtext(x,y+8.,ch3)
	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(u2,down,dis)
	write(dis(1:),1334)num,ac,dc,epoch
1334	format(i3,') ',a,1x,a,'(',f6.1,')')
	call pgsci(9)
	call pgtext(u2,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	goto 110	
	endif

c right button****************
	if(ich.eq.132)then
	  close (99)
	  call pgend
	  return
	endif

c middle  button**************
	if(ich.eq.130)then
	call pgsci(4)
	xi=xycoef(1,1)*(2049-x)+xycoef(1,2)*(2049-y)+xycoef(1,3)
        xn=xycoef(2,1)*(2049-x)+xycoef(2,2)*(2049-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
	 call pgpoint(1,x,y,16)
	 call pgtext(x,y+8.,ch3)
104	format(i3,') ',a,1x,a,'(',f6.1,') false star')
	call pgsci(0)
	call pgtext(u2,down,dis)
	write(dis(1:),104)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(u2,down,dis)
	write(99,'(a,2f7.1)')dis,x,y
	write(*,'(a,2f7.1)')dis,x,y
	goto 110
	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)
	x=2049-x
	y=2049-y
	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
	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)
	  call pgpoint(1,x,y,18)
	  call pgtext(x,y+8.,ch3)
105	format(i3,') ',a,1x,a,'(',f6.1,') input alpha delta')
	call pgsci(0)
	call pgtext(u2,down,dis)
	write(dis(1:),105)num,ac,dc,epoch
	call pgsci(8)
	call pgtext(u2,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 hardcopy ****************
	if(ch.eq.'p'.or.ch.eq.'P')then
	  write(*,'(a)')bell
	  call pgend
	  close (99)
	
	  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=15.85
          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(u2,u1,v1,v2)
          call pgbox('bc',0,0,'bc',0,0)
  	call pgxgray(pam,2048,2048,2049-iu2,2049-iu1,iv1,iv2,white,black,2)
	  call pgsch(0.8)
          call pglabel(' ','Decl.  East',f1//ns)
	call pgline(4,sx,sy,1)
	call pgtext(u2-50.,sy(1),c3)
          x1=u2
	  dy=15.
          call pgsch(.6)
          i=0
149       y1=v1-30.
150       i=i+1
          read(99,'(a,2f7.1)',end=151)dis,x,y
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
          write(dis(31:),131)x,y
131	format(' ',2f8.1) 
          call pgtext(x1,y1,dis)
          y1=y1-dy
c write No. on picture
          call pgtext(x,y+8,dis(1:3))
c see if second column for writting
          if(i.eq.15)then
            x1=xc-20
            goto 149
          endif
          goto 150
151       close (99)
          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(ich.lt.128)write(*,'(a)')bell
	goto 110
	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	  do 5 j=ky-4,ky+4
c5	  write(*,'(9i6)')(map(i,j),i=kx+4,kx-4,-1)
	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-25+center_1(xmap,sky,51,1)
	y=iy-25+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
        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)
	character*80 head(72),f1*8
	byte bell
	data bell/7/
	do 10 indexpos=1,72
10	if(head(indexpos)(1:8).eq.f1)return
        if(indexpos.eq.73)write(*,'(1x,a)')bell
	end

	subroutine readata(a,f1,ibit)
	integer c(720),a(1)
	character*30 f1,ch*2880
	equivalence (c,ch)
	nn=720						
	open(1,file=f1,status='old',access='direct',recl=nn*4)
	read(1,rec=1)ch
	n1=2048
	n2=2048
	n=1
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 shrink4(a,n1,n2,a)
	else
	  call swap2(a,n1*n2*2)
	endif
	write(*,*)
	end

	subroutine shrink4(a,n1,n2,b)
	real a(n1*n2)
	integer*2 b(n1*n2)
	do 10 i=1,n1*n2
	x=a(i)
        if(x.gt.32767.)x=32767.
        if(x.lt.0.)x=0.
10	b(i)=x
	end

	function inqh(f1,n1,n2,head,ibit)
	character*(*) f1
	character*80 head(72)
	logical inqh
	k=index(f1,'.')
	if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
	inquire(file=f1,exist=inqh)
	if(.not.inqh)return
	open(1,file=f1,status='old',access='direct',recl=5760)
	read(1,rec=1)head
	read(head(2)(22:),*)ibit
	read(head(4)(22:),*)n1	!NAXIS1
	read(head(5)(22:),*)n2	!NAXIS2
	ipos=indexpos(head,'OBJECT  ')
	write(*,'(1x,a)')head(ipos)(1:70)       	! object
	ipos=indexpos(head,'DATE-OBS')
        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   ')
        if(ipos.ne.73)write(*,'(1x,a)')head(ipos)(1:78)
        ipos=indexpos(head,'TELESCOP')
        if(ipos.ne.73)head(ipos)(12:28)='BAO 60-90 SCHIMDT'
        ipos=indexpos(head,'BZERO   ')
        if(ipos.ne.73)head(ipos)(1:1)='P'
        ipos=indexpos(head,'BSCALE  ')
        if(ipos.ne.73)head(ipos)(1:1)='P'
        ipos=indexpos(head,'TIME    ') 
        if(ipos.ne.73)head(ipos)(42:57)='OF START OF OBS.' 
        ipos=indexpos(head,'HA      ')
        if(ipos.ne.73)head(ipos)(45:61)='OF MIDDLE OF OBS.' 
	close (1)
	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