parameter (n=5000)
        real d(n,8),xm(n),dm(n),pho(8)
        character*40 f1,f2,f3,f4,nf4(5),f5,aa
	character*1  ch
	character*3  ch3
	character*12 ac,dc
	character*78 dis,disno,disnold,ph
	byte bell
	real a(2048*2048)
	real b(2048*2048)
	real*8 a8(8),adcoef(2,3)
	real sx(5),sy(5)
	character*80 head(72)	
	logical logi

	if(iargc().lt.3)then
	write(*,*)
	write(*,*)'       ******** DISplay star**********2000.7 '
	write(*,*)
	write(*,*)
	write(*,*)' Usage: no5 no1_out frame1 frame2 [listfile] [dssfit]'
	write(*,*)
	write(*,*)'                  listfile included <=5 fits_file'
	write(*,*)
	stop
	endif
	call getarg(1,f1)
	inquire(file=f1,exist=logi)
	if(.not.logi)then
	  write(*,*)f1,'file not found !'
	  stop
	endif
	call getarg(2,f2)
	k=index(f2,'.')
	if(k.eq.0)f2=f2(1:lnblnk(f2))//'.fit'
	inquire(file=f2,exist=logi)
	if(.not.logi)then
	  write(*,*)f2,'file not found !'
	  stop
	endif
	call getarg(3,f3)
	k=index(f3,'.')
	if(k.eq.0)f3=f3(1:lnblnk(f3))//'.fit'
	inquire(file=f3,exist=logi)
	if(.not.logi)then
	  write(*,*)f3,'file not found !'
	  stop
	endif
	m4=0
	if(iargc().ge.4)then
	  call getarg(4,f4)
	  inquire(file=f4,exist=logi)
	  if(.not.logi)then
	    write(*,*)f4,'file not found !'
	    stop
 	  endif
          call gethead(f2,head)               ! mother file
          call geta8(head,a8)
          call xytoad(a8,adcoef)
          x=a8(7)
          y=a8(8)
	  open(1,file=f4,status='old')
5	  m4=m4+1
	  if(m4.eq.6)goto 6
	  read(1,'(a)',end=6)nf4(m4)
	  inquire(file=nf4(m4),exist=logi)
	  if(.not.logi)then
	    write(*,*)nf4(m4),'file not found !'
	    stop
 	  endif
          call gethead(nf4(m4),head)            
          ipos=indexpos(head,'A87     ')
          read(head(ipos)(20:),*)alpha
          read(head(ipos+1)(20:),*)delta
          call standc(x,y,alpha,delta,xi,xn)
          sx(m4)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)-1024.5
          sy(m4)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)-1024.5
	  goto 5
6	  m4=m4-1
	  close(1)
	endif	
c	write(*,*)(sx(i),sy(i),i=1,m4)
	if5=0
	if(iargc().eq.5)then
	  call getarg(5,f5)
	  inquire(file=f5,exist=logi)
	  if(.not.logi)then
            write(*,*)f1,'file not found !'
          else
	    call dsshead(f5)
	  endif
	endif    
	
	call readfits(f2,a,n1,n2)
	write(*,*)'readin ...',f2
	call readfits(f3,b,n3,n4)
	write(*,*)'readin ...',f3
	
	bell=7

        open(1,file=f1,status='old')
	m=0
	x=0.
10	m=m+1
	read(1,8,end=11)i1,i2,x,j1,j2,y,(d(m,i),i=3,8)
8	format(2(1x,i2,1x,i2,1x,f5.2),2f7.2,4f7.1)
	d(m,1)=i1+i2/60.+x/3600.
	d(m,2)=j1+j2/60.+y/3600.
        xm(m)=(d(m,4)+d(m,3))/2.
        dm(m)=d(m,4)-d(m,3)
	x=x+d(m,2)
        goto 10
11	m=m-1
        write(aa(1:),*)' epoch=2000. total number:',m
	close(1)

	call pgbegin(0,'/xw',1,1)
c	call pgslw(2)
	call pgsch(0.8)
	call pgenv(11.5,20.5,-4.,2.,0,0)
	call pglabel(' ','delta mag',aa)
	do 20 i=1,m
20	call pgpoint(1,xm(i),dm(i),21)

	num=0
	kk=0
110     call pgcurse(x,y,ch)
        ii=ichar(ch)

c right button****************
        if(ii.eq.132)then
          call pgend
	  stop
        endif


c left button****************
        if(ii.eq.128)then
	  r=1e10
	  do 30 i=1,m
	  rr=(xm(i)-x)**2+(dm(i)-y)**2
          if(rr.lt.r)then
            k=i
            r=rr
          endif
30	  continue
	  if(r.gt.0.01)then
	    write(*,'(a)')bell
	    goto 110
	  endif
c	  write(*,*)(d(k,i),i=1,8)
          num=num+1
	  call pgsci(9)
	  if(kk.eq.k)then
	    num=num-1
	    call pgsci(0)
 	    call pgpoint(1,xm(k),dm(k),21)
	  endif
          write(ch3(1:),'(i3)')num
	  call pgsch(0.6)
          call pgtext(xm(k)-0.02,dm(k)+0.03,ch3)
	  call pgsch(0.8)
          call toms(d(k,1),ac,1)
          call toms(d(k,2),dc,0)
	  write(*,12)num,ac,dc,(d(k,i),i=3,8)
12	  format(i3,') ',a,1x,a,2f7.2,4f7.1)
          call pgsci(0)
          call pgtext(11.6,-1.5,dis)
          write(dis(1:),12)num,ac,dc,(d(k,i),i=3,8)
	  if(kk.eq.k)then
	    k=0
	    num=num-1
	  else 
            call pgsci(9)
            call pgtext(11.6,-1.5,dis)
	    pho(2)=disp(a,n1,n2,d(k,5),d(k,6),80,400,1)
	    pho(3)=disp(b,n3,n4,d(k,7),d(k,8),190,400,0)
	    if(if5.ne.0)call dsstar(f5,d(k,1),d(k,2),520,400)
	    call usnosa(num,d(k,1),d(k,2),disno)
	    call pgsci(0)
	    call pgtext(11.6,-4.4,disnold)
            call pgsci(9)
	    call pgtext(11.6,-4.4,disno)
	    write(*,'(a)')disno
	    disnold=disno
	    do 13 i=1,m4
	    i1=d(k,5)+sx(i)+0.5
	    i2=d(k,6)+sy(i)+0.5
13	    pho(i+3)=disp4(nf4(i),i1,i2,i*110-30,525)
	    kk=k
	    pho(1)=d(k,3)-pho(2)
	    do 14 i=2,m4+3
14	    pho(i)=pho(i)+pho(1)
	    call pgsci(0)
	    call pgtext(11.6,-3.0,ph)
	    call pgsci(9)
	    write(ph(1:),4)(pho(i),i=1,m4+3)
	    write(*,4)(pho(i),i=1,m4+3)
4	    format('  Radius photometry: ',8f7.2)
	    call pgtext(11.6,-3.0,ph)
	  endif
	  goto 110
	endif

c display data
        if(ch.eq.'d')then
	  r=1e10
	  do 130 i=1,m
	  rr=(xm(i)-x)**2+(dm(i)-y)**2
          if(rr.lt.r)then
            k=i
            r=rr
          endif
130	  continue
	  if(r.gt.0.01)then
	    write(*,'(a)')bell
	    goto 110
	  endif

	  call disd(a,n1,n2,d(k,5),d(k,6))
	  call disd(b,n3,n4,d(k,7),d(k,8))

	  goto 110
	endif

	goto 110
	end

cccccccccc
	subroutine dsstar(f1,xx,yy,ix,iy)
	character*40 f1
	integer*2 map(101*101),ma(101*101)
	common /dss/nn1,nn2,a1,a2,d1,d2,npix1,npix2
        xc=xx
	yc=yy
	ip=0
10	x=(nn1-1)/(a2-a1)*(xc-a1)+1.
	y=(nn2-1)/(d2-d1)*(yc-d1)+1.
	if(x.lt.1.or.y.lt.1.or.x.gt.nn1.or.y.gt.nn2)then
	  do 20 i=1,101*101
20	  map(i)=0
	  goto 209
	endif
        x1=x+npix1
        y1=y+npix2
        call amdpos(x1,y1,x2,y2)
        x=(x-1)/(x2-a1)*(xc-a1)+1.
        y=(y-1)/(y2-d1)*(yc-d1)+1.
cccccccccccc
        x1=x+npix1
        y1=y+npix2
        call amdpos(x1,y1,x2,y2)
        x1=xx-x2
	y1=yy-y2
c	   write(*,'(3f10.6)')xx,x2,x1
c	   write(*,'(3f10.6)')yy,y2,y1
        if(abs(x1).gt.0.0001.or.abs(y1).gt.0.0015)then
	   xc=xc+x1
	   yc=yc+y1
 	   ip=ip+1
	   if(ip.lt.10)goto 10
	   do 30 i=1,101*101
30	   map(i)=0
	   goto 209
	endif
cccccccccccc
	i1=x+0.5
	i2=y+0.5
	call readd(f1,map,i1,i2,nn1,nn2)
	ip=0
        do 203 i=1,101*101
	if(map(i).eq.32767)goto 203
	ip=ip+1
	ma(ip)=map(i)
203	continue
	call whitexblack(ma,1,ip,white,sigma)
209     call ximage2(map,101,ix,iy,white,white+20.*sigma)
	end

	function disp4(f1,i1,i2,ix,iy)
	character*40 f1
	integer*2 map(101*101),ma(101*101)
	real xmap(101*101)
	call readcc(f1,xmap,i1,i2)
	ip=0
        do 203 i=1,101*101
	if(xmap(i).lt.0.)then
	  map(i)=32767
	else
   	  ip=ip+1
	  j=xmap(i)*.25
	  if(j.gt.32766)j=32766
	  ma(ip)=j
          map(i)=j
	endif
203	continue
	call whitexblack(ma,1,ip,white,sigma)
	disp4=photo(xmap,white*4.)
        call ximage2(map,101,ix,iy,white*0.95,white+20.*sigma)
	end

	function photo(a,white)
	real a(101,101)
	xx=0
        do 206 j1=-5,+5
        do 206 j2=-5,+5
	k=j1*j1+j2*j2
	if(k.gt.25)goto 206
	x=a(51+j1,51+j2)-white
	if(x.gt.0)xx=xx+x
206	continue
	photo=25.-2.5*alog10(xx)
	end
	

	subroutine disd(a,n1,n2,zx,zy)
	real a(n1,n2)
	integer map(15,15)
	i1=nint(zx)
	i2=nint(zy)
        do 203 j1=i1-7,i1+7
        do 203 j2=i2-7,i2+7
	x=a(j1,j2)
	if(x.gt.99999.)x=99999.
	if(x.lt.-9999.)x=-9999.
        map(j1-i1+8,j2-i2+8)=x		
203	continue
	write(*,*)
	do 204 j=1,15
204	write(*,'(15i5)')(map(i,j),i=1,15)
	write(*,*)
	end

cccccccccccccccccccccc

	function disp(a,n1,n2,zx,zy,ix,iy,jp)
	real a(n1,n2),uy(101*101),ux(101*101),xmap(101*101)
	integer*2 map(101,101),ma(101*101)
	ip=0
	i1=nint(zx)
	i2=nint(zy)
	j=0
        do 203 j1=i1-50,i1+50
        do 203 j2=i2-50,i2+50
	j=j+1
	if(j1.lt.1 .or. j1.gt.n1 .or. j2.lt.1 .or. j2.gt.n2)then
	  map(j1-i1+51,j2-i2+51)=32767
	  xmap(j)=-1.
	else 
          white=a(j1,j2)
	  if(white.lt.0.)white=0.
	  xmap(j)=white
	  ip=ip+1
          uy(ip)=white*0.25
          ux(ip)=sqrt((j1-zx)**2+(j2-zy)**2)
          k=white*.25
          if(k.gt.32766)k=32766
          map(i1-j1+51,i2-j2+51)=k
	  ma(ip)=k
	endif
203	continue
	call whitexblack(ma,1,ip,white,sigma)
	disp=photo(xmap,white*4.)
	black=white+20.*sigma
	do 204 i=1,101
	map(i,1)=0
	map(i,101)=0
	map(1,i)=0
	map(101,i)=0
204	continue

        call ximage2(map,101,ix,iy,white*.95,black)

	do 205 i=1,101
	do 205 j=1,101
205	map(i,j)=11
	x=0.
	do 206 i=1,101*101
	i1=ux(i)*16.+1
	if(i1.gt.101)goto 206
	if(uy(i).gt.x)x=uy(i)
206	continue

	if(jp.eq.1)scale=101./x

	do 207 i=1,101*101
	i1=ux(i)*16.+1
	if(i1.gt.101)goto 207
	i2=uy(i)*scale
	if(i2.gt.101)i2=101
	if(i2.lt.1)i2=1
	map(i1,102-i2)=1		
207	continue		
        call ximage3(map,101,ix+220,iy)
	end

        subroutine ximage2(map,n,ix,iy,black,white)
        integer*2 map(n,n)
        character dummy
        real rbuf(101*101+2)
        n16=86
        d=n16/(black-white)
        rbuf(1)=ix
        rbuf(2)=iy
        i6=2
        do 30 j=1,n
        do 30 i=1,n
        k=nint(d*(map(i,j)-white))
        if(k.lt.1)k=1
        if(k.gt.n16)k=n16
        i6=i6+1
30      rbuf(i6)=k+15
        i6=n+2
        call xwdriv(27,rbuf,i6,dummy)
        call pgebuf
        end

        subroutine ximage3(map,n,ix,iy)
        integer*2 map(n,n)
        character dummy
        real rbuf(101*101+2)
        rbuf(1)=ix
        rbuf(2)=iy
	m=2
        do 30 j=1,n
        do 30 i=1,n
	m=m+1
30      rbuf(m)=map(i,j)
        m=n+2
        call xwdriv(27,rbuf,m,dummy)
        call pgebuf
        end

        subroutine gethead(f1,head)
        character*80 head(72),f1*40
        open(52,file=f1,status='old',access='direct',recl=5760)
        read(52,rec=1)head
        close (52)
        end

        subroutine geta8(head,a8)
        real*8 a8(8)
        character*80 head(72)
        ipos=indexpos(head,'A81     ')-1
        do 10 i=1,8
10      read(head(i+ipos)(12:),*)a8(i)
        end

        function indexpos(head,f1)
        character*80 head(72),f1*8
        do 10 indexpos=1,72
10      if(head(indexpos)(1:8).eq.f1)return
        end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine dsshead(f1)
	character*(*) f1
	real*8 pse,xz,yz,px,py,a_c,d_c,ax(13),ay(13),cx,cy
	common /head/pse,xz,yz,a_c,d_c,px,py,ax,ay,cx,cy
	common /dss/nn1,nn2,a1,a2,d1,d2,npix1,npix2
	call getdsshead(f1)
	x1=1+npix1
	y1=1+npix2
	call amdpos(x1,y1,a1,d1)
	x1=nn1+npix1
	y1=nn2+npix2
	call amdpos(x1,y1,a2,d2)
c by using a1,a2; d1,d2 calculate new star
	if(a1.lt.a2)a1=a1+24.
	end

	subroutine getdsshead(f1)
	character*(*) f1
	character*80 head(108)
	logical logi
	real*8 pse,xz,yz,px,py,a_c,d_c,ax(13),ay(13),cx,cy
	common /head/pse,xz,yz,a_c,d_c,px,py,ax,ay,cx,cy
	common /dss/nn1,nn2,a1,a2,d1,d2,npix1,npix2
	cy=atan(1d0)/45d0
	cx=cy*15.
	if(index(f1,'.').eq.0)f1=f1(1:lnblnk(f1))//'.fits'
	inquire(file=f1,exist=logi)
	if(.not.logi)stop "File not found !"
	open(41,file=f1,status='old',access='direct',recl=8640)
	read(41,rec=1)head
	close (41)
	call ival(head,'NAXIS1  =',nn1) 
	call ival(head,'NAXIS2  =',nn2) 
	call ival(head,'CNPIX1  =',npix1) 
	call ival(head,'CNPIX2  =',npix2) 
        call fval(head,'PLTSCALE=',pse)
	call fval(head,'XPIXELSZ=',xz) 
	call fval(head,'YPIXELSZ=',yz) 
	call cval(head,'PLTRAH  =')
	call fval(head,'PPO3    =',px)
	call fval(head,'PPO6    =',py)
	call aval(head,'AMDX1   =',ax)
	call aval(head,'AMDY1   =',ay)
	end

	subroutine ival(head,word,k)
	character*80 head(110)
	character*9 word
	do 10 i=1,110
	if(head(i)(1:9).eq.word)goto 20
10	continue
	stop "not DSS fits file !"
20	read(head(i)(21:),*)k
	end

	subroutine rval(head,word,tmp)
	character*80 head(110)
	character*9 word
	do 10 i=1,110
	if(head(i)(1:9).eq.word)goto 20
10	continue
	stop "not DSS fits file !"
20	read(head(i)(11:),*)tmp
	end

	subroutine fval(head,word,tmp)
	character*80 head(110)
	character*9 word
	real*8 tmp
	do 10 i=1,110
	if(head(i)(1:9).eq.word)goto 20
10	continue
	stop "not DSS fits file !"
20	read(head(i)(11:),*)tmp
	end

	subroutine cval(head,word)
	character*80 head(110)
	character*9 word
	real*8 a,d
	real*8 pse,xz,yz,px,py,a_c,d_c,ax(13),ay(13),cx,cy
	common /head/pse,xz,yz,a_c,d_c,px,py,ax,ay,cx,cy
	do 10 i=1,110
	if(head(i)(1:9).eq.word)goto 20
10	continue
	stop "not DSS fits file !"
20	read(head(i)(22:),*)j
	read(head(i+1)(22:),*)k
	read(head(i+2)(11:),*)x    ! 95,3  modify for north plate 
	a=j+k/60.+x/3600.
	read(head(i+4)(22:),*)j
	read(head(i+5)(22:),*)k
	read(head(i+6)(11:),*)x
	d=j+k/60.+x/3600.
	if(head(i+3)(12:12).eq.'-')d=-d
	a_c=a*cx
	d_c=d*cy
	end

	subroutine aval(head,word,a)
	character*80 head(110)
	character*9 word
	real*8 a(1)
	do 10 i=1,110
	if(head(i)(1:9).eq.word)goto 20
10	continue
	stop "not DSS fits file !"
20	do 30 k=1,13
	read(head(i)(11:),*)a(k)
30	i=i+1
	end

        subroutine amdpos(x,y,ra,dec)
	real*8 ox,oy,ox2,oy2,ox3,oy3,xi,yi,n,d,r
	real*8 pse,xz,yz,px,py,a_c,d_c,ax(13),ay(13),cx,cy
	common /head/pse,xz,yz,a_c,d_c,px,py,ax,ay,cx,cy
        ox = (px - x*xz)/1000.0
        oy = (y*yz - py)/1000.0
        ox2 = ox*ox
        oy2 = oy*oy
        ox3 = ox*ox2
        oy3 = oy*oy2
        xi = ax(1)*ox + ax(2)*oy + ax(3) + ax(4)*ox2 + ax(5)*ox*oy +
     *       ax(6)*oy2 + ax(7)*(ox2+oy2) + ax(8)*ox3 + ax(9)*ox2*oy +
     *       ax(10)*ox*oy2 + ax(11)*oy3  + ax(12)*ox*(ox2+oy2) +
     *       ax(13)*ox*(ox2+oy2)*(ox2*oy2)
        yi = ay(1)*oy + ay(2)*ox + ay(3) + ay(4)*oy2 + ay(5)*oy*ox +
     *       ay(6)*ox2 + ay(7)*(oy2+ox2) + ay(8)*oy3 + ay(9)*oy2*ox +
     *       ay(10)*oy*ox2 + ay(11)*ox3  + ay(12)*oy*(oy2+ox2) +
     *       ay(13)*oy*(oy2+ox2)*(oy2*ox2)
	xi=xi*cy/3600.
	yi=yi*cy/3600.
	n=xi/dcos(d_c)
	d=1-yi*dtan(d_c)
	r=datan2(n,d)+a_c
	n=dcos(r-a_c)
	d=d/(yi+dtan(d_c))
	dec=datan(n/d)/cy
	ra=r/cx
	if(ra.lt.0.)ra=ra+24.
	end

	subroutine usnosa(num,alpha,delta,disno)
	character*(*)disno
	character*8 f1
	character*12 ac,dc
        integer*4 a(3)
	real*8 xx
	byte eb,er
        i=(delta+90.)*10.
        i=i/75*75
        write(f1(1:),"('zone',i4.4)")i
        open(1,file='/USNOSA10/'//f1//'.cat',
     *     status='old',access='direct',recl=12)
        open(2,file='/USNOSA10/'//f1//'.acc',status='old')
          read(f1(5:),*)ifsa1
          ifsa1=ifsa1/75*100000000
        do 20 k=1,96
        read(2,*)x1,i
        if(x1+0.25.ge.alpha)goto 21
20      continue
21      close(2)
        x1=alpha-0.0004
        x2=alpha+0.0004
	d1=delta-0.005
	d2=delta+0.005
	cosx=cos(delta*3.141592654/180.)*15.
	rr=10.
110     read(1,rec=i,err=100)a
	  i=i+1
          call swap4(a,12)
          xx=a(1)
          ra=xx/5400000.
	if(ra.lt.x1)goto 110
        if(ra.gt.x2)goto 100
          xx=a(2)
          de=xx/360000.-90.
        if(de.gt.d2.or.de.lt.d1)goto 110
	r=sqrt(((alpha-ra)*cosx)**2+(delta-de)**2)
	if(r.lt.rr)then
	  rr=r
          k=a(3)
          if(k.lt.0)k=-k
          j=mod(k,1000)
          if(j.lt.0)j=j+256
          aa=ra
          dd=de
          er=mod(k,1000)
          k=k/1000
          eb=mod(k,1000)
          numsa1=ifsa1+i-1
        endif
        goto 110
100     close (1)
	if(rr.ge.10.)then
	  write(disno(1:),'(a)')' Cannot find USNOSA10'
	else
        call toms(aa,ac,1)
        call toms(dd,dc,0)
        j=eb
        if(j.lt.0)j=j+256
        xbbb=j*0.1
        j=er
        if(j.lt.0)j=j+256
        xrrr=j*0.1
        i=numsa1
        j=i
        i=i/100000000*75
        j=mod(j,100000000)
        write(disno(1:),203)num,ac,dc,i,j,xbbb,xrrr
203     format(i3,') ',a,1x,a,'(2000)  USNO',i4.4,'.',i8.8,
     *  f5.1,'B ',f4.1,'R')
	endif
        end