character*80 a,f1*70
        real*8 a8(8),adcoef(2,3)
        character*80 head(72)
	logical logi

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

	if(iargc().lt.1)then
	  write(*,*)' ****** Usage: pho dat_file    (2000.7)**'
	  write(*,*)
	  write(*,*)'data_file example: default RADIUS OF APERTURE is 5'
          write(*,*)'p1021129TA03c032.fit   '
	  write(*,*)'07:07:23.35  63:38:07.0'
	  write(*,*)'07:06:53.30  63:36:16.2'
	  write(*,*)'07:07:22.43  63:28:17.2'
	  write(*,*)'07:07:48.53  63:32:15.5'
	  write(*,*)'p1021129TA03c033.fit  4'
	  write(*,*)'07:07:23.35  63:38:07.0'
	  write(*,*)'07:06:53.30  63:36:16.2'
	  write(*,*)'07:06:38.36  63:29:01.8'
	  write(*,*)'07:07:22.43  63:28:17.2'
	  stop
	endif
	call getarg(1,f1)
        inquire(file=f1,exist=logi)
	if(.not.logi)then
	  write(*,*)' file not found!'
	  stop
	endif
	open(1,file=f1,status='old')
10	read(1,'(a)',end=100)a
	i1=lnblnk(a)
        i2=index(a,' ')
	if(ichar(a(1:1)).gt.57)then
	  i3=5
	  if(i1.gt.i2)read(a(i2:),*)i3
	  write(*,*)a(1:i2),i3
	  f1=a(1:i2)
          call gethead(f1,head)             
          call geta8(head,a8)
          call xytoad(a8,adcoef)
	else
	  call itohd(a,alpha,delta,dummy)
	  alp=alpha*cx
          del=delta*cy
	a87=a8(7)
	a88=a8(8)
          call standc(a87,a88,alp,del,xi,xn)
          alp=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
          del=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
          ia=2049.-alp+.5
          id=2049.-del+.5
	  if(ia.gt.2048.or.ia.lt.0.or.id.gt.2048.or.id.lt.0)then
	    write(*,'(a,2i5,a)')a(1:i1),ia,id,' out of frame'
	    goto 10
	  endif
	  xmag=disp4(f1,ia,id,i3)
	  write(*,'(a,2i5,f6.2)')a(1:i1),ia,id,xmag
	endif
	goto 10
100	close(1)
	end

        subroutine gethead(f1,head)
        character*80 head(72),f1*70
        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

        function disp4(f1,i1,i2,i3)
        character*70 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).le.0.)then
          map(i)=32767
        else
          ip=ip+1
          j=xmap(i)
          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,i3)
        end

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