c 1995,1,16 modify
c ~40000 --> 512*256 --> 512*512*2 --> display
c       ujcomp     ucomp
c rawdata: 2880 + 17 (table) + ~40000
c x,y gape 2 pixel, 1996,9,2 in xytoad sub
	
	integer*2 b(512*512)
	byte c256(512*256),table(17)
	character*80 f1,head(36)
	real*8 a8(8)

	if(iargc().eq.0)then
	  write(*,*)
	  write(*,*)'         ******** display compress CCD file ***** '
	  write(*,*)' 2008.10'
	  write(*,*)'         position RMS < 1 arcsec'
          write(*,*)'         left button:  pick out a star'
	  write(*,*)'         middle button:  Zoom & Unzoom'
	  write(*,*)'         right button:   Exit'
	  write(*,*)'         f (zoom case): make false star'
	  write(*,*)'         if input -epoch, then hardcopy'
	  write(*,*)'         if input e or q, then interrupt'
	  write(*,*)
	  stop
	endif
	call getarg(1,f1)
		write(*,*)f1
	call getdata(f1,head,b,n,n1,n2,m1,m2,table)
	call geta8(head,a8)
	call ujcomp(b,c256,table,n,m1,m2)
	call ucomp4(c256,b,m1,m2)
	call turnx2(b,m1,m2,-180)
	call dis(b,f1,m1,m2,n1,n2,a8)
	end

	subroutine geta8(head,a8)
	real*8 a8(8)
	character*80 head(36)
	if(head(17)(15:15).eq.' ')return
	do 10 i=1,8
10	read(head(i+16)(12:),*)a8(i)
	end

        subroutine dis(map,f1,n1,n2,m1,m2,a8)
	integer pgbegin
        integer*2 map(n1,n2)
	real*8 a8(8)
	character*50 ns,ns1
        character*80 f1
	character*12 ac,dc
	character*3 ch3
	character*1 ch
	byte bell
	bell=7
	pi=4.*atan(1.0)
        cx=pi/12.
        cy=pi/180.

88      epoch=2000.0
        write(*,89)bell,epoch
89      format(1x,a,'output epoch, [',f6.1,'] ',$)
        read(*,'(a)')ns
        if(ns(1:1).eq.'e'.or.ns(1:1).eq.'q')stop 'abort task'
        if(ns(1:1).ne.' ')then
          read(ns(1:),*,err=88)epoch
	  ihp=0
	  if(epoch.lt.0)then
	    ihp=1
	    epoch=-epoch
	  endif
          if(epoch.lt.1900.or.epoch.gt.2100.)goto 88
	  if(ihp.eq.1)write(*,*)' wait for printing'
        endif

	if(ihp.eq.1)then
	  call unlink('pgplot.hp')
	  icall=pgbegin(0,'/hp',1,1)
	  call pgslw(2)
	  goto 11
	endif
        icall=pgbegin(0,'/xw',1,1 )
	call pgpap(8.5,1.015)
	call pgsci(0)
11      call pgenv(float(m1),1.,float(m2),1.,1,0)
        call pgsci(1)
	do 10 i=1,n2
	map(1,i)=15
10	map(n1,i)=15
	do 20 i=1,n1
	map(i,1)=12
20	map(i,n2)=12
        call pglabel('R.A. north','Decl. west (turn 180 degree)',f1)

c99      call pg_image4(map,n1,n2)
99      call pgxgray(map,n1,n2,1,n1,1,n2,0.,16.,2)
	call pgsch(2.)
	call pgsci(2)
	x=a8(7)/cx
	y=a8(8)/cy
        call astprs(x,y,2000.,xx,yy,epoch)
        call toms(xx,ac,1)
        call toms(yy,dc,0)
        write(ns1(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
	x=m1*0.5
	y=m2*0.5
	call pgpoint(1,x,y,5)
	call pgsch(0.8)
	call pgsci(3)
        call pgtext(2049*2.-900*2.,2049*2.-2080*2.,ns1)
	if(ihp.eq.1)then
	  call pgiden
	  call pgend
	  call system('print pgplot.hp')
	  stop 'normal finished'
	endif	
	iz=m1/n1
	ip=1
100	x=4097.-x
	y=4097.-y
	call pgcurse(x,y,ch)
	x=4097.-x
	y=4097.-y
	if(ch.eq.'A')ii=128
	if(ch.eq.'D')ii=130
	if(ch.eq.'X')ii=132
	if(ch.eq.'A')then
c	  call box(map,n1,n2,x,y)
	  call star(map,n1,n2,x,y,1)
c	  call box(map,n1,n2,x,y)
	  call star(map,n1,n2,x,y,1)
c	  write(*,*)x,y
	  call coord(a8,x,4097-y,ac,dc,epoch)
	  call pgsci(0)
	  call pgtext(4097.-20.,4097.+199.,ns)
	  write(ns(1:),"(i3,': ',a,1x,a,' (',f6.1,')')")ip,ac,dc,epoch
	  call pgsci(3)
	  call pgtext(4097.-20.,4097.+199.,ns)
	  write(*,*)ns
	  write(ch3(1:),'(i3)')ip
	  call pgsci(2)
	  call pgtext(4097.-x,4097.-(y+60.),ch3)
	  ip=ip+1
	endif
	if(ch.eq.'D')then
	  call pgsci(0)
          call pgtext(2049*2.-900*2.,2049*2.-2080*2.,ns1)
	  call pgtext(2049*2.-10*2.,2049*2.+99*2.,ns)
	  call zoom(map,n1,n2,x,y,iz,ii,m1,m2,a8,epoch)
	endif
	if(ch.eq.'X')then
	  call pgend
	  stop 'normal finished !'
	endif
	if(ii.eq.134)goto 99
	goto 100
        end
	
	subroutine zoom(map,n1,n2,x1,y1,iz,ii,m1,m2,a8,epoch)
	integer*2 map(n1,n2)
	real*8 a8(8)
	character ch*1,ch3*3,ac*12,dc*12,ns*50,ns1*50
	write(*,*)
        call pgsci(3)
c*****************************
c	xx=x1-2.
c	yy=y1-2.5
	xx=x1
	yy=y1
	call coord(a8,2049*2.-xx,yy,ac,dc,epoch)
        write(ns1(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch
	call pgtext(2049*2.-900*2.,2049*2.-2080*2.,ns1)
	ip=1
	if(iz.eq.2)call zoom2(map,n1,n2,x1,y1,iz,m1,m2)
	if(iz.eq.4)call zoom4(map,n1,n2,x1,y1,iz,m1,m2)
	xx=1024*2.
	yy=1024*2.
10	xx=2049*2.-xx
	yy=2049*2.-yy
	call pgcurse(xx,yy,ch)
	xx=2049*2.-xx
	yy=2049*2.-yy
	if(ch.eq.'A')ii=128
	if(ch.eq.'D')ii=130
	if(ch.eq.'X')ii=132
	if(ch.eq.'f')then
	  x=(xx-1024*2.)/4.+x1
	  y=(yy-1024*2.)/4.+y1
	  call coord(a8,2049*2.-x,y,ac,dc,epoch)
	  call pgsci(0)
	  call pgtext(2049*2.-10*2.,2049*2+99*2.,ns)
	  write(ns(1:),"(i3,': ',a,1x,a,' (',f6.1,')   f')")
     *            ip,ac,dc,epoch
	  call pgsci(3)
	  call pgtext(2049*2-10*2.,2049*2+99*2.,ns)
	  write(*,*)ns
	  write(ch3(1:),'(i3)')ip
	  call pgsci(2)
	  call pgpoint(1,xx+8.,yy,18)
	  call pgtext(2049*2.-(xx-15*2.),2049*2.-(yy+25*2.),ch3)
	  ip=ip+1
	endif
	if(ii.eq.128)then
	  x=(xx-1024*2.)/4.+x1
	  y=(yy-1024*2.)/4.+y1
	  call star(map,n1,n2,x,y,0)
	  call star(map,n1,n2,x,y,0)
	  xx=(x-x1)*4.+1024*2.
	  yy=(y-y1)*4.+1024*2. 
	  call coord(a8,2049*2.-x,y,ac,dc,epoch)
	  call pgsci(0)
	  call pgtext(2049*2.-10*2.,2049*2.+99*2.,ns)
	  write(ns(1:),"(i3,': ',a,1x,a,' (',f6.1,')')")ip,ac,dc,epoch
	  call pgsci(3)
	  call pgtext(2049*2.-10*2.,2049*2.+99*2.,ns)
	  write(*,*)ns
	  write(ch3(1:),'(i3)')ip
	  call pgsci(2)
	  call pgtext(2049*2.-(xx-15*2.),2049*2.-(yy+25*2.),ch3)
	  ip=ip+1
	endif
	if(ii.eq.130)ii=134
	if(ii.eq.132)then
	        call pgend
		stop 'normal finished'
	endif
	if(ii.ne.134)goto 10
	write(*,*)
	call pgsci(0)
	call pgtext(2049*2-10*2.,2049*2+99*2.,ns)
	end

	subroutine coord(a8,x,y,ac,dc,epoch)
	character*12 ac,dc
	real*8 a8(8)
	call xy_rad(x,y,a8,a,d)
        call astprs(a,d,2000.,aa,dd,epoch)
	call toms(aa,ac,1)
	call toms(dd,dc,0)
	end

	subroutine xy_rad(x11,y11,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.)
	x=x11+2.
	y=y11-2.
        xi=a8(1)*x+a8(3)*y+a8(5)
        xn=a8(2)*x+a8(4)*y+a8(6)
c becauese of 4-->1pixel, diff 2 pixel in coord.
        rc=a8(7)                    ! center of ccd (in rad.)o
        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

	subroutine box(map,n1,n2,x,y)
	integer*2 map(n1,n2)
	write(*,*)x,y
	ix=x/4+0.5
	iy=y/4+0.5
	do 20 j=iy+4,iy-4,-1
20	write(*,'(9i4)')(map(i,j),i=ix-4,ix+4)
	end

	subroutine zoom4(map,n1,n2,x,y,iz,m1,m2)
	integer*2 map(n1,n2),map1(512,512)
	  ix=x/iz-0.5
	  iy=y/iz-0.5
	  id=512/8-1
	  if(ix+id.gt.510)ix=510-id
	  if(iy+id.gt.510)iy=510-id
	  if(ix-id.lt.1)ix=1+id
	  if(iy-id.lt.1)iy=1+id
	
	  ii=-3
	  do 110 i=ix-id,ix+id
	  ii=ii+iz
	  jj=-3
	  do 110 j=iy-id,iy+id
	  jj=jj+iz
110	  map1(ii,jj)=map(i,j)	
	  do 120 i=3,512,iz
	  do 120 j=3,512,iz
120	  map1(i,j)=(map1(i-2,j-2)+map1(i-2,j+2)+
     *               map1(i+2,j-2)+map1(i+2,j+2))/4
	  k=0
	  do 130 i=3,512,2
	  k=2-k
	  do 130 j=3+k,512,iz
130	  map1(i,j)=(map1(i,j-2)+map1(i,j+2)+
     *               map1(i+2,j)+map1(i-2,j))/4
	  do 220 i=2,512,2
	  do 220 j=2,512,2
220	  map1(i,j)=(map1(i-1,j-1)+map1(i-1,j+1)+
     *               map1(i+1,j-1)+map1(i+1,j+1))/4
	  k=0
	  do 230 i=2,512
	  k=1-k
	  do 230 j=2+k,512,2
230	  map1(i,j)=(map1(i,j-1)+map1(i,j+1)+
     *               map1(i+1,j)+map1(i-1,j))/4
	  do 240 i=1,512
	  do 240 j=1,5
	  map1(j,i)=0
240	  map1(i,j)=0
c          call pg_image4(map1,n1,n2)
	  call pgxgray(map1,n1,n2,1,n1,1,n2,0.,16.,2)
	call pgsch(2.)
	call pgsci(2)
	xx=m1*0.5
	yy=m2*0.5
	call pgpoint(1,xx,yy,5)
	call pgsch(0.8)
	end

	subroutine zoom2(map,n1,n2,x,y,iz,m1,m2)
	integer*2 map(n1,n2),map1(512,512)
	  ix=x/iz
	  iy=y/iz
	  ii=-1
	  do 110 i=ix-127,ix+127
	  ii=ii+2
	  jj=-1
	  do 110 j=iy-127,iy+127
	  jj=jj+2
110	  map1(ii,jj)=map(i,j)	
	  do 120 i=2,511,2
	  do 120 j=2,512,2
120	  map1(i,j)=(map1(i-1,j-1)+map1(i-1,j+1)+
     *               map1(i+1,j-1)+map1(i+1,j+1))/4
	  k=0
	  do 130 i=2,510
	  k=1-k
	  do 130 j=2+k,512,2
130	  map1(i,j)=(map1(i,j-1)+map1(i,j+1)+
     *               map1(i+1,j)+map1(i-1,j))/4
c          call pg_image4(map1,n1,n2)
	 call pgxgray(map1,n1,n2,1,n1,1,n2,0.,16.,2)
	end

	subroutine ujcomp(a,b,table,n,m1,m2)
	byte a(n),b(512*256),table(17)
	integer pb
	mm=m1*m2/2
	pb=0
10	pb=pb+1
	if(pb.gt.mm)return
	k=kgetbit(1,a)
	if(k.eq.0)then
	  b(pb)=0              ! 0
	  goto 10
	endif	
	k=kgetbit(2,a)
	if(k.lt.3)then
	  b(pb)=table(k+1) 	! 3 3_data
	  goto 10
	endif	
  	k=kgetbit(1,a)
  	if(k.eq.0)then
    	  b(pb)=kgetbit(8,a)  	! 4 for char 
	  goto 10
	endif	
  	k=kgetbit(2,a)
  	j=3
  	if(k.lt.2)then
    	  b(pb)=table(j+k+1)     ! 6 2_data 
	  goto 10
	endif	
20	j=j+2
  	if(j.eq.15)then
	  b(pb)=table(j+m+1)
	  goto 10
	endif	
  	m=kgetbit(1,a)
  	if(k.eq.2)then
	  b(pb)=table(j+m+1)
	  goto 10
	endif	
	if(m.eq.0)k=2
	goto 20
	end

	function kgetbit(m,a)
	byte a(1)
	integer pc,pa
  	integer*2 i,k
	byte c
	data pc,pa/0,0/
  	k=0
	do 10 i=1,m
    	if(pc.eq.0)then
	  pa=pa+1
	  c=a(pa)
	endif
	k=k+k
        if(c.lt.0)k=k+1
	c=c+c
	pc=pc+1
10	if(pc.eq.8)pc=0
	kgetbit=k
	end

	subroutine getdata(f1,head,b,n,n1,n2,m1,m2,table)
	character*80 f1,head(36)
	byte b(1),table(17)
	logical logi
	k=index(f1,'.')
	if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
	inquire(file=f1,exist=logi)
	if(.not.logi)stop 'file not found !'
	open(3,file=f1,status='old',access='direct',recl=2880)
	read(3,rec=1)head
	close (3)
	read(head(3)(67:),*,err=10)n
	read(head(4)(20:),*)n1
	read(heaD(4)(69:),*)m1
	read(head(5)(20:),*)n2
	read(head(5)(69:),*)m2
	open(3,file=f1,status='old',access='direct',recl=2897+n)
	read(3,rec=1)head,table,(b(i),i=1,n)
	close(3)
	return
10	write(*,*)' not my format'
	stop
	end
	
	subroutine ucomp4(c,b,m1,m2)
	integer*2 b(512*512)
	byte c(512*256)
	i=1
	do 10 j=1,m1*m2/2
	k=c(j)
	if(k.lt.0)k=k+256
	j1=k/16
	j2=mod(k,16)
	b(i)=j1
	b(i+1)=j2
10	i=i+2
	end

	subroutine pg_image4(map,n1,n2)
	integer*2 map(n1,n2)
	real rbuf(520)
	character*16 c
	character*1 cdummy

	i1=1
	call pgqinf('TYPE',c,k)
	if(c(1:2).eq.'XW')then
	  ii=(641-n1)*.5
	  jj=(641-n2)*.5
	  rbuf(1)=ii
	  do 40 j=1,n2
	  rbuf(2)=j+jj
	  i6=2
	  do 30 i=1,n1
	  i6=i6+1
30	  rbuf(i6)=90-map(i,j)*5
40	  call xwdriv(26,rbuf,i6)
	  call pgebuf
	endif	
	if(c(1:2).eq.'HP')then
	call pgqwin(x1,x2,y1,y2)
	call grtxy0(.false.,x1,y1,x1,y1)
        ix1=nint(x1)
        iy1=nint(y1)
	call grtxy0(.false.,x2,y2,x2,y2)
	ix2=nint(x2)
	iy2=nint(y2)
	sx=(n1-1)/(x2-x1)
	sy=(n2-1)/(y2-y1)
	ip=3
	scale=18.
	ip2=ip+ip
	ip05=ip/2
        if(ix2.gt.ix1)then
                ix11=ix1+ip
                ix22=ix2-ip05
                ixstep=ip
        else
                ix11=ix1-ip
                ix22=ix2+ip05
                ixstep=-ip
        endif
        if(iy2.gt.iy1)then
                iy11=iy1+ip
                iy22=iy2-ip05
                iystep=ip2
                iyd=ip
        else
                iy11=iy1-ip
                iy22=iy2+ip05
                iystep=-ip2
                iyd=-ip
        endif
        ratio=scale/15.
        k=-1
        do 50 ix=ix11,ix22,ixstep
        k=-k
        kd=0
        if(k.eq.1)kd=iyd
	dx=sx*(ix-ix1)
	i1=dx
	i1=i1+1
	i2=i1+1
	x2=mod(dx,1.)
	x1=1.-x2
	rbuf(1)=ix
        do 50 iy=iy1+kd,iy22,iystep
	rbuf(2)=iy
	dy=sy*(iy-iy1)
	j1=dy
	j1=j1+1
	j2=j1+1
	y2=mod(dy,1.)
	y1=1.-y2
	z=map(i1,j1)*x1*y1+map(i1,j2)*x1*y2+
     *    map(i2,j1)*x2*y1+map(i2,j2)*x2*y2
        rbuf(3)=ratio*z
	rbuf(4)=scale
	call hpdriv(18,rbuf,6,cdummy,1)
50	continue
	endif
        end

	subroutine star(map,n1,n2,x,y,m)
	integer*2 map(n1,n2),sky
	byte bell
	real b(-15:15)
	bell=7
	oldx=x
	oldy=y
	sky=3
	i5=5
	if(m.eq.0)i5=4
        ix=nint(x*.25)
        iy=nint(y*.25)
cccccc x center    weight_center method
        do 30 i=ix-i5,ix+i5
        z=0.
        do 40 j=iy-i5,iy+i5
40      if(map(i,j).gt.sky)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
	if(z.ne.0.)then
        s=s/z
        x=(ix+s)*4.
	endif
cccccc y center
        do 130 j=iy-i5,iy+i5
        z=0.
        do 140 i=ix-i5,ix+i5
140     if(map(i,j).gt.sky)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
	if(z.ne.0.)then
        s=s/z
        y=(iy+s)*4.
	endif
	if(m.eq.1)return
	if((oldx-x)**2+(oldy-y)**2.gt.8)then
	  x=oldx
	  y=oldy
	  write(*,'(1x,a)')bell
	endif
	end