program difcal
c
c calculate differential magnitudes
c written 9-April-1997 aah
c this version outputs an unsorted file
c
      PARAMETER (MAXSTARS=10000)
      PARAMETER (XINDEF=99.999,XCHECK=99.0)
      real*8 ra(MAXSTARS),dec(MAXSTARS),hjd(MAXSTARS)
      real*4 smag(MAXSTARS,2),xmag(MAXSTARS,3)
      integer i,n,nmast
      character flist*80,fname*80,fout*80
      common /blka/ ra,dec,hjd,smag,xmag,nmast
c
      print *,'Enter input filelist: '
      read (5,'(a)') flist
      open (unit=1,file=flist,status='old')
      print *,'Enter master star filename: '
      read (5,'(a)') fname
      open (unit=2,file=fname,status='old')
      read (2,*)
      i = 1
100   continue
      read (2,900,end=200) ra(i),dec(i),smag(i,1),n,smag(i,2)
900   format (f8.4,1x,f9.4,16x,f7.3,7x,i4,2x,f7.3)
      if (n.lt.3) goto 100
      i=i+1
      goto 100
200   continue
      nmast = i - 1
      write(6,905) nmast
905   format (i5,' Stars with at least 3 observations.')
      close(2)
      print *,'Enter output filename: '
      read (5,'(a)') fout
      open (unit=3,file=fout,status='new')
300   continue
      read (1,'(a)',end=400) fname
      open (unit=2,file=fname,status='old')
      write (6,901) fname
901   format ('Processing file: ',a)
      call readfile
      call formdif
      close(2)
      goto 300
400   continue
      stop
      end

      subroutine readfile
c
      PARAMETER (MAXSTARS=10000)
      PARAMETER (XINDEF=99.999,XCHECK=99.0)
      real*8 ra(MAXSTARS),dec(MAXSTARS),hjd(MAXSTARS)
      real*8 rax,decx,hjd1,hjd2,hjd3
      real*4 dra,ddec,fwhmx1,fwhmx2,fwhmx3,xmag1,
     $       xmag2,xmag3,dmag1,dmag2,dmag3,x1,x2,x3,
     $       y1,y2,y3,err,minerr
      real*4 smag(MAXSTARS,2),xmag(MAXSTARS,3)
      integer i,j,nmast
      common /blka/ ra,dec,hjd,smag,xmag,nmast
c
      minerr = (15./3600.)**2
      do i=1,nmast
        xmag(i,1) = XINDEF
        xmag(i,2) = XINDEF
        xmag(i,3) = XINDEF
      enddo
      read (2,*)
      read (2,*)
100   continue
         read (2,902,end=200) rax,decx,dra,ddec,
     $  hjd1,fwhmx1,fwhmy1,xmag1,dmag1,x1,y1,
     $  hjd2,fwhmx2,fwhmy2,xmag2,dmag2,x2,y2,
     $  hjd3,fwhmx3,fwhmy3,xmag3,dmag3,x3,y3
902     format (4f9.4,f11.4,6f7.3/
     $    5x,f11.4,6f7.3/
     $    5x,f11.4,6f7.3)
        if (xmag1.gt.XCHECK.or.((xmag2.gt.XCHECK).and.
     $     (xmag3.gt.XCHECK))) goto 100
        do i=1,nmast
           err = (ra(i)-rax)**2 + (dec(i)-decx)**2
           if (err.lt.minerr) then
              xmag(i,1) = xmag1
              xmag(i,2) = xmag2
              xmag(i,3) = xmag3
              hjd(i) = hjd1
              goto 100
           endif
         enddo
         goto 100
200   continue
      return
      end

      subroutine formdif
c
c calculate differential magnitudes for all objects
c
      PARAMETER (MAXSTARS=10000)
      PARAMETER (XINDEF=99.999,XCHECK=99.0)
      real*8 sum1,sum2,ra(MAXSTARS),dec(MAXSTARS),hjd(MAXSTARS)
      real*4 x(100),smag(MAXSTARS,2),xmag(MAXSTARS,3)
      integer indx(100),i,ii,j,k,n1,n2,nmast,j1,j2
      common /blka/ ra,dec,hjd,smag,xmag,nmast
c
      do i=1,nmast
        if (xmag(i,1).lt.XCHECK) then
           j1 = max(1,i-50)
           j2 = min(nmast,j1+100)
           k=0
           do j=j1,j2
             if (j.ne.i) then
               k = k+1
c              x(k) = abs(dec(i)-dec(j))
               x(k) = smag(j,1)
               indx(k) = j
             endif
           enddo
           call sort (k,x,indx)
           n1 = 0
           n2 = 0
           sum1 = 0.
           sum2 = 0.
           do j=1,10
             ii = indx(j)
             if (xmag(ii,1).lt.XCHECK) then
               sum1 = sum1 + (xmag(i,1)-xmag(ii,1))+smag(ii,1)
               n1 = n1 + 1
             endif
             if (xmag(ii,2).lt.XCHECK.and.xmag(i,2).lt.XCHECK
     $         .and.smag(ii,2).lt.XCHECK) then
               sum2 = sum2 + (xmag(i,2)-xmag(ii,2))+smag(ii,2)
               n2 = n2 + 1
             endif
             if (xmag(ii,3).lt.XCHECK.and.xmag(i,3).lt.XCHECK
     $         .and.smag(ii,2).lt.XCHECK) then
               sum2 = sum2 + (xmag(i,3)-xmag(ii,3))+smag(ii,2)
               n2 = n2 + 1
             endif
           enddo
           if (n1.gt.0) then
              sum1 = sum1 / float(n1)
           else
              sum1 = XINDEF
           endif
           if (n2.gt.0) then
              sum2 = sum2 / float(n2)
           else
              sum2 = XINDEF
           endif
           write (3,900) i,ra(i),dec(i),hjd(i),sum1,sum2
900        format (i5,1x,f8.4,2x,f9.4,f11.4,2f7.3)
        endif
      enddo
      return
      end

      SUBROUTINE SORT (n,x,index)
c
c heapsort of array x with corresponding index array index
c from numerical recipies
c note: for real*4 values
c
      DIMENSION x(n),index(n)
      k = n/2+1
      ir = n
10    continue
      IF (k.gt.1) THEN
        k = k-1
        xx = x(k)
        ii = index(k)
      ELSE
        xx = x(ir)
        ii = index(ir)
        x(ir) = x(1)
        index(ir) = index(1)
        ir = ir-1
        IF (ir.eq.1) THEN
          x(1) = xx
          index(1) = ii
          RETURN
        ENDIF
      ENDIF
      i = k
      j = k+k
20    IF (j.le.ir) THEN
        IF (j.lt.ir) THEN
          IF (x(j).lt.x(j+1)) j = j+1
        ENDIF
        IF (xx.lt.x(j)) THEN
          x(i) = x(j)
          index(i) = index(j)
          i = j
          j = j+j
        ELSE
          j = ir+1
        ENDIF
      GOTO 20
      ENDIF
      x(i) = xx
      index(i) = ii
      GOTO 10
      END