program difcal
c
c calculate differential magnitudes
c written 9-April-1997 aah
c this version outputs an unsorted file
c mod 20-May-1997  put filename on line
c mod 23-May-1997  improve comparison selection
c mod 24-May-1997  add error calc to dif phot + rejection
c
      PARAMETER (MAXSTARS=10000,MAXCOMP=100)
      PARAMETER (XINDEF=99.999,XCHECK=99.0)
      real*8 ra(MAXSTARS),dec(MAXSTARS),hjd(MAXSTARS)
      real*4 smag(MAXSTARS,2),xmag(MAXSTARS,3)
      integer indx(MAXCOMP),i,n,nmast
      character flist*80,fname*80,fout*80
      common /blka/ ra,dec,hjd,smag,xmag,indx,nmast,fname
c
      print *,'DIFCAL Version 1.3'
      print *
      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')
      write (3,901)
901   format('# ID',5x,'RA',8x,'DEC',8x,'JD',7x,'V',6x,'I',
     $  5x,'sigV',3x,'sigI',3x,'nc',5x,'fname')
300   continue
      read (1,'(a)',end=400) fname
      open (unit=2,file=fname,status='old')
      write (6,902) fname
902   format ('Processing file: ',a)
      call readfile
      call formdif
      close(2)
      goto 300
400   continue
      stop
      end

      subroutine readfile
c
c read the input collated/transformed star lists
c
      PARAMETER (MAXSTARS=10000,MAXCOMP=100)
      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,err,minerr
      real*4 smag(MAXSTARS,2),xmag(MAXSTARS,3)
      integer indx(MAXCOMP),i,j,nmast
      character*80 fname
      common /blka/ ra,dec,hjd,smag,xmag,indx,nmast,fname
c
c minerr is the coordinate error radius for matching
c
      minerr = (15./3600.)**2
c
c use magnitudes as flags to indicate if data is avaliable
c
      do i=1,nmast
        xmag(i,1) = XINDEF
        xmag(i,2) = XINDEF
        xmag(i,3) = XINDEF
      enddo
c skip header
      read (2,*)
      read (2,*)
c
c loop over number of objects in list
c
100   continue
         read (2,902,end=200) rax,decx,dra,ddec,
     $  hjd1,fwhmx1,fwhmy1,xmag1,dmag1,
     $  hjd2,fwhmx2,fwhmy2,xmag2,dmag2,
     $  hjd3,fwhmx3,fwhmy3,xmag3,dmag3
902     format (4f9.4,f11.4,4f7.3/
     $    5x,f11.4,4f7.3/
     $    5x,f11.4,4f7.3)
c
c note: for now, only accept data that has V plus at least one I point
c
        if (xmag1.gt.XCHECK.or.((xmag2.gt.XCHECK).and.
     $     (xmag3.gt.XCHECK))) goto 100
c
c for each object, find its corresponding entry in the master phot file
c
        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
c note: using V-frames JD as the date for now
              hjd(i) = hjd1
              goto 100
           endif
         enddo
         goto 100
200   continue
      return
      end

      subroutine formdif
c
c calculate differential magnitudes for all objects
c heart of program.  Uses up to 100 nearby stars as comparisons
c and does ensemble differential photometry.
c One iteration pass to remove discrepant points.
c
      PARAMETER (MAXSTARS=10000,MAXCOMP=100)
      PARAMETER (XINDEF=99.999,XCHECK=99.0)
      real*8 ra(MAXSTARS),dec(MAXSTARS),hjd(MAXSTARS)
      real*4 smag(MAXSTARS,2),xmag(MAXSTARS,3)
      real*4 sig1,sig2,ave1,ave2,s1err,s2err
      integer indx(MAXCOMP),i,nmast,nc,iflag
      character*80 fname
      common /blka/ ra,dec,hjd,smag,xmag,indx,nmast,fname
c
      do i=1,nmast
        if (xmag(i,1).lt.XCHECK) then
         call findcomp(i,nc)
         if (nc.ge.2) then
           s1err = XINDEF
           s2err = XINDEF
           ave1 = 0.
           ave2 = 0.
           iflag = 0
           call calccomp(i,nc,s1err,s2err,ave1,sig1,ave2,sig2,iflag)
           s1err = sig1*1.05
           s2err = sig2*1.05
           iflag = 1
           call calccomp(i,nc,s1err,s2err,ave1,sig1,ave2,sig2,iflag)
           write (3,900) i,ra(i),dec(i),hjd(i),ave1,ave2,sig1,
     $       sig2,nc,fname(1:12)
900        format (i5,1x,f8.4,2x,f9.4,f11.4,2f7.3,2f7.3,i4,2x,a12)
         endif
        endif
      enddo
      return
      end

      subroutine calccomp(i,nc,s1err,s2err,ave1,sig1,ave2,sig2,iflag)
c
c do actual comparison calculation
c
      PARAMETER (MAXSTARS=10000,MAXCOMP=100)
      PARAMETER (XINDEF=99.999,XCHECK=99.0)
      real*8 sum1,sum2,ra(MAXSTARS),dec(MAXSTARS),hjd(MAXSTARS)
      real*8 sumx1,sumx2,sumn
      real*4 smag(MAXSTARS,2),xmag(MAXSTARS,3)
      real*4 sig1,sig2,ave1,ave2,s1err,s2err,yy,xx
      integer indx(MAXCOMP),i,ii,j,k,n1,n2,nmast,nc,iflag
      character*80 fname
      common /blka/ ra,dec,hjd,smag,xmag,indx,nmast,fname
c
       n1 = 0
       n2 = 0
       sum1 = 0.
       sum2 = 0.
       sumx1 = 0.
       sumx2 = 0.
       do j=1,nc
         ii = indx(j)
         if (xmag(ii,1).lt.XCHECK) then
           xx = (xmag(i,1)-xmag(ii,1))+smag(ii,1)
           if(iflag.eq.0) then
             sum1 = sum1 + xx
             sumx1 = sumx1 + xx*xx
             n1 = n1 + 1
           else
             yy = abs(ave1-xx)
             if (yy.le.s1err) then
               sum1 = sum1 + xx
               sumx1 = sumx1 + xx*xx
               n1 = n1 + 1
             endif
           endif
         endif
         if (xmag(ii,2).lt.XCHECK.and.xmag(i,2).lt.XCHECK
     $     .and.smag(ii,2).lt.XCHECK) then
           xx = (xmag(i,2)-xmag(ii,2))+smag(ii,2)
           if(iflag.eq.0) then
             sum2 = sum2 + xx
             sumx2 = sumx2 + xx*xx
             n2 = n2 + 1
           else
             yy = abs(ave2-xx)
             if (yy.le.s2err) then
               sum2 = sum2 + xx
               sumx2 = sumx2 + xx*xx
               n2 = n2 + 1
             endif
           endif
         endif
         if (xmag(ii,3).lt.XCHECK.and.xmag(i,3).lt.XCHECK
     $     .and.smag(ii,2).lt.XCHECK) then
           xx = (xmag(i,3)-xmag(ii,3))+smag(ii,2)
           if(iflag.eq.0) then
             sum2 = sum2 + xx
             sumx2 = sumx2 + xx*xx
             n2 = n2 + 1
           else
             yy = abs(ave2-xx)
             if (yy.le.s2err) then
               sum2 = sum2 + xx
               sumx2 = sumx2 + xx*xx
               n2 = n2 + 1
             endif
           endif
         endif
       enddo
       if (n1.ge.2) then
         sumn = float(n1)
         sig1 = dsqrt(abs((sumn*sumx1-sum1*sum1)/
     $           (sumn*(sumn-1.))))
       else
         sig1 = XINDEF
       endif
       if (n2.ge.2) then
         sumn = float(n2)
         sig2 = dsqrt(abs((sumn*sumx2-sum2*sum2)/
     $           (sumn*(sumn-1.))))
       else
         sig2 = XINDEF
       endif
       if (n1.gt.0) then
          ave1 = sum1 / float(n1)
       else
          ave1 = XINDEF
       endif
       if (n2.gt.0) then
         ave2 = sum2 / float(n2)
       else
         ave2 = XINDEF
       endif
      return
      end

      subroutine findcomp(i,k)
c
c find a reasonable set of comparison stars for this object
c
      PARAMETER (MAXSTARS=10000,MAXCOMP=100)
      PARAMETER (XINDEF=99.999,XCHECK=99.0)
      real*8 ra(MAXSTARS),dec(MAXSTARS),hjd(MAXSTARS)
      real*4 smag(MAXSTARS,2),xmag(MAXSTARS,3)
      real*4 dx1,dx2,dy,maxerr
      integer indx(MAXCOMP),i,k,nmast,j1,j2
      character*80 fname
      common /blka/ ra,dec,hjd,smag,xmag,indx,nmast,fname
c
c maxerr is 20arcmin in any direction
c
      maxerr = (20./60.)
      j1 = max(1,i-1)
      j2 = min(nmast,i+1)
      dx1 = 0.0
      dx2 = 0.0
      k = 0
100   continue
      if (j1.ne.1.and.j2.ne.nmast.and.
     $   dx1.le.maxerr.and.dx2.le.maxerr) then
        if (j1.ne.1) then
          dx1 = abs(ra(i)-ra(j1))
          dy = abs(dec(i)-dec(j1))
          if (dx1.le.maxerr.and.dy.le.maxerr.and.
     $      xmag(j1,1).lt.12.0) then
            k = min((k+1),100)
            indx(k) = j1
          endif
        endif
        if (j2.ne.nmast) then
          dx2 = abs(ra(i)-ra(j2))
          dy = abs(dec(i)-dec(j2))
          if (dx2.le.maxerr.and.dy.le.maxerr.and.
     $      xmag(j2,1).lt.12.0) then
            k = min((k+1),100)
            indx(k) = j2
          endif
        endif
        if (k.eq.100) return
        j1 = max(1,j1-1)
        j2 = min(nmast,j2+1)
        goto 100
      endif
      if (k.ge.10) return
c
c fewer than 10 comp stars; adjust limits once
c
      maxerr = (30./60.)
      j1 = max(1,i-1)
      j2 = min(nmast,i+1)
      dx1 = 0.0
      dx2 = 0.0
      k = 0
200   continue
      if (j1.ne.1.and.j2.ne.nmast.and.
     $   dx1.le.maxerr.and.dx2.le.maxerr) then
        if (j1.ne.1) then
          dx1 = abs(ra(i)-ra(j1))
          dy = abs(dec(i)-dec(j1))
          if (dx1.le.maxerr.and.dy.le.maxerr.and.
     $      xmag(j1,1).lt.12.0) then
            k = min((k+1),100)
            indx(k) = j1
          endif
        endif
        if (j2.ne.nmast) then
          dx2 = abs(ra(i)-ra(j2))
          dy = abs(dec(i)-dec(j2))
          if (dx2.le.maxerr.and.dy.le.maxerr.and.
     $      xmag(j2,1).lt.12.0) then
            k = min((k+1),100)
            indx(k) = j2
          endif
        endif
        if (k.eq.100) return
        j1 = max(1,j1-1)
        j2 = min(nmast,j2+1)
        goto 200
      endif
      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