program master
c
c make master photometric starlist from transformed files
c This only handles 10K stars, so not entire stripe
c input is a list of files to process.  Only stars with V & I
c  are taken from each input file as those are the only stars
c  for which proper photometric transformations are available.
c output is the master photometric file list, where
c  each star has at least 3 measures
c written 23-March-1997 aah
c
      parameter (MAXSTARS=10000)
      parameter (XINDEF=99.999,XCHECK=99.0)
      real*8 ra(MAXSTARS), dec(MAXSTARS),rasq(MAXSTARS),
     $   decsq(MAXSTARS), mag(MAXSTARS,2), magsq(MAXSTARS,2),
     $   magsum(MAXSTARS,2),
     $   minerr,err,hjd1,hjd2,hjd3,rax,decx,xnr,xnp,
     $   ramean,decmean,sum(20)
      real*4  fwhmx1,fwhmx2,fwhmx3,fwhmy1,fwhmy2,fwhmy3
      real*4  xmag1,xmag2,xmag3,dmag1,dmag2,dmag3
      real*4 x1,x2,x3,y1,y2,y3,dra,ddec,dr(MAXSTARS),dd(MAXSTARS)
      real*4 magmean1,magmean2,dmag(MAXSTARS,2),x(20)
      real*4 rasig(MAXSTARS),decsig(MAXSTARS),errfac
      real*4 smag(20),serr(20)
      integer nr(MAXSTARS),np(MAXSTARS),i,j,n,indx(MAXSTARS),
     $     nused(MAXSTARS,2),nsum(20),nmags,k
      character*80 flist,fname,fout
      common /blka/ ra,dec,rasq,decsq,magsum,magsq,dr,dd,nr,np
      common /blkb/ rax,decx,xmag1,xmag2,xmag3,dra,ddec,n
c
      minerr = (15./3600.)**2
      errfac = 1.0
      print *,'MASTER Version 1.0'
      print *,'Enter input filelist: '
      read (5,'(a)') flist
      open (unit=1,file=flist,status='old')
      print *,'Enter output filename: '
      read (5,'(a)') fout
      open (unit=3,file=fout,status='new')
      write (3,920)
920   format ('#  RA',7x,'DEC',6x,'dra',5x,'ddec',4x,'V',5x,'dV',
     $   5x,'nV',5x,'I',5x,'dI',5x,'nI')
      print *,'Enter statistics filename: '
      read (5,'(a)') fout
      open (unit=22,file=fout,status='new')
      open (unit=11,file='sigma.txt',status='old')
      read (11,941) nmags
941   format (i5)
      do i=1,nmags
        read (11,942) smag(i),serr(i)
942     format (f5.2,5x,f7.3)
        serr(i) = serr(i)*4.0
      enddo
      close(11)
c
c read in initial set of stars
c
      read (1,'(a)') fname
      open (unit=2,file=fname,status='old')
      write (6,912) fname
      read (2,*)
      read (2,*)
      n=0
100   continue
        read (2,901,end=200) rax,decx,dra,ddec,
     $    hjd1,fwhmx1,fwhmy1,xmag1,dmag1,
     $    hjd2,fwhmx2,fwhmy2,xmag2,dmag2,
     $    hjd3,fwhmx3,fwhmy3,xmag3,dmag3
901     format(4f9.4,f11.4,4f7.3/
     $    5x,f11.4,4f7.3/
     $    5x,f11.4,4f7.3)
        call addstar
        goto 100
200   continue
      close(2)
c
c now read in set of stars to match with first
c
300   continue
      read (1,'(a)',end=600) fname
      open (unit=2,file=fname,status='old')
      write (6,912) fname
912   format ('Processing file ',a)
      read (2,*)
      read (2,*)
400   continue
        read (2,901,end=500) rax,decx,dra,ddec,
     $    hjd1,fwhmx1,fwhmy1,xmag1,dmag1,
     $    hjd2,fwhmx2,fwhmy2,xmag2,dmag2,
     $    hjd3,fwhmx3,fwhmy3,xmag3,dmag3
c if no V & I magnitudes, ignore this star
        if (xmag1.gt.XCHECK.or.((xmag2.gt.XCHECK).and.
     $      (xmag3.gt.XCHECK))) goto 400
c else (brute force) see if the star is already in the list,
c by comparing its coordinates.  Error is minerr (15arcsec radius)
        do i=1,n
          xnr = float(nr(i))
          err = (ra(i)/xnr-rax)**2 + (dec(i)/xnr-decx)**2
          if (err.lt.minerr) then
c we've found a match. add current star's data to sums
            nr(i) = nr(i) + 1
            ra(i) = ra(i) + rax
            rasq(i) = rasq(i) + rax*rax
            dec(i) = dec(i) + decx
            decsq(i) = decsq(i) + decx*decx
            magsum(i,1) = magsum(i,1) + xmag1
            magsq(i,1) = magsq(i,1) + xmag1*xmag1
c now look at I0 and I2, and add them in if present
            if (xmag2.lt.XCHECK) then
              np(i) = np(i) + 1
              magsum(i,2) = magsum(i,2) + xmag2
              magsq(i,2) = magsq(i,2) + xmag2*xmag2
            endif
            if (xmag3.lt.XCHECK) then
              np(i) = np(i) + 1
              magsum(i,2) = magsum(i,2) + xmag3
              magsq(i,2) = magsq(i,2) + xmag3*xmag3
            endif
            goto 400
          endif
        enddo
c star not found in current list...add it to list
        call addstar
        goto 400
500   continue
      close(2)
      goto 300
600   continue
c
c now have all stars, compute means
c
      do i=1,n
        ramean = ra(i)
        decmean = dec(i)
        magmean1 = magsum(i,1)
        if (nr(i).ge.2) then
          xnr = float(nr(i))
          ramean = ramean / xnr
          decmean = decmean / xnr
          magmean1 = magmean1 / xnr
          dra = dsqrt(abs((xnr*rasq(i)-ra(i)*ra(i))/
     $           (xnr*(xnr-1.))))
          ddec = dsqrt(abs((xnr*decsq(i)-dec(i)*dec(i))/
     $           (xnr*(xnr-1.))))
          dmag1 = dsqrt(abs((xnr*magsq(i,1)-magsum(i,1)*
     $           magsum(i,1))/(xnr*(xnr-1.))))
        else
          dra = dr(i)
          ddec = dd(i)
          dmag1 = XINDEF
        endif
        magmean2 = magsum(i,2)
        if (np(i).ge.2) then
          xnp = float(np(i))
          magmean2 = magmean2 / xnp
          dmag2 = dsqrt(abs((xnp*magsq(i,2)-magsum(i,2)*
     $       magsum(i,2))/(xnp*(xnp-1.))))
        else
          dmag2 = XINDEF
        endif
c save values for second pass
        ra(i) = ramean
        dec(i) = decmean
        rasig(i) = dra
        decsig(i) = ddec
        mag(i,1) = magmean1
        mag(i,2) = magmean2
        magsum(i,1) = 0.
        magsum(i,2) = 0.
        magsq(i,1) = 0.
        magsq(i,2) = 0.
        dmag(i,1) = dmag1*errfac
        dmag(i,2) = dmag2*errfac
        nused(i,1) = 0
        nused(i,2) = 0
        do k=1,nmags-1
          if (magmean1.ge.smag(k).and.magmean1.lt.smag(k+1)
     $      .and.dmag1.gt.serr(k).and.dmag1.lt.XCHECK) then
             write (22,934) ramean,decmean,magmean1,dmag1,
     $           magmean2,dmag2
934          format (f8.4,1x,f9.4,4f7.3)
             goto 888
          endif
        enddo
888     continue
      enddo
c
c now do second iteration.  We have means, but they may be
c contaminated by bad points.  Here we remove points that
c are more than 2*sigma from the mean.
c
      write (6,914)
914   format (' Starting iteration pass.')
c
c first, sort by ra
c
      do i=1,n
        indx(i) = i
      enddo
      call sort (n,ra,indx)
      rewind(1)
700   continue
      read (1,'(a)',end=870) fname
      open (unit=2,file=fname,status='old')
      write (6,912) fname
      read (2,*)
      read (2,*)
800   continue
        read (2,901,end=850) rax,decx,dra,ddec,
     $    hjd1,fwhmx1,fwhmy1,xmag1,dmag1,
     $    hjd2,fwhmx2,fwhmy2,xmag2,dmag2,
     $    hjd3,fwhmx3,fwhmy3,xmag3,dmag3
c if no V & I magnitudes, ignore this star
        if (xmag1.gt.XCHECK.or.((xmag2.gt.XCHECK).and.
     $      (xmag3.gt.XCHECK))) goto 800
c else (brute force) see if the star is already in the list,
c by comparing its coordinates.  Error is minerr (15arcsec radius)
         do j=1,n
          i = indx(j)
          err = (ra(j)-rax)**2 + (dec(i)-decx)**2
          if (err.lt.minerr) then
c we've found a match. add current star's data to sums
            if (abs(mag(i,1)-xmag1).lt.dmag(i,1)) then
              magsum(i,1) = magsum(i,1) + xmag1
              magsq(i,1) = magsq(i,1) + xmag1*xmag1
              nused(i,1) = nused(i,1) + 1
            endif
c now look at I0 and I2, and add them in if present
            if (xmag2.lt.XCHECK) then
              if (abs(mag(i,2)-xmag2).lt.dmag(i,2)) then
                magsum(i,2) = magsum(i,2) + xmag2
                magsq(i,2) = magsq(i,2) + xmag2*xmag2
                nused(i,2) = nused(i,2) + 1
              endif
            endif
            if (xmag3.lt.XCHECK) then
              if (abs(mag(i,2)-xmag3).lt.dmag(i,2)) then
                magsum(i,2) = magsum(i,2) + xmag3
                magsq(i,2) = magsq(i,2) + xmag3*xmag3
                nused(i,2) = nused(i,2) + 1
              endif
            endif
            goto 800
          endif
        enddo
850   continue
      close(2)
      goto 700
870   continue
c
c now have all stars, compute means
c
      write (6,913)
913   format ('Now form means and standard deviations')
      do j=1,n
        i = indx(j)
        if (nused(i,1).ge.3) then
          xnr = nused(i,1)
          nr(i) = nused(i,1)
          mag(i,1) = magsum(i,1) / xnr
          dmag(i,1) = dsqrt(abs((xnr*magsq(i,1)-magsum(i,1)*
     $           magsum(i,1))/(xnr*(xnr-1.))))
        else
          dmag(i,1) = dmag(i,1)/errfac
        endif
        if (nused(i,2).ge.3) then
          xnr = nused(i,2)
          np(i) = nused(i,2)
          mag(i,2) = magsum(i,2) / xnr
          dmag(i,2) = dsqrt(abs((xnr*magsq(i,2)-magsum(i,2)*
     $           magsum(i,2))/(xnr*(xnr-1.))))
        else
          dmag(i,2) = dmag(i,2)/errfac
        endif
        if (nr(i).ge.3.and.np(i).ge.3) 
     $  write (3,908) ra(j),dec(i),rasig(i),decsig(i),mag(i,1),
     $    dmag(i,1),nr(i),mag(i,2),dmag(i,2),np(i)
908     format (f8.4,1x,f9.4,2f8.4,2f7.3,i4,2x,2f7.3,i4)
      enddo
      close(1)
      close(2)
      close(3)
      x(1) = 6.
      do i=2,16
        x(i) = x(i-1) + 0.5
        sum(i) = 0.
        nsum(i) = 0
      enddo
      do j=1,n
        do i=2,16
          if (mag(j,1).ge.x(i-1).and.mag(j,1).le.x(i).and.
     $        dmag(j,1).lt.XCHECK) then
             sum(i-1) = sum(i-1) + dmag(j,1)
             nsum(i-1) = nsum(i-1) + 1
             goto 880
          endif
        enddo
880     continue
      enddo
      do i=2,16
        if (nsum(i).eq.0) then
          write (22,924) x(i)
        else
           sum(i) = sum(i) / float(nsum(i))
           write (22,925) x(i),nsum(i),sum(i)
        endif
       enddo
924   format (f5.2,'    0','  0.000')
925   format (f5.2,i5,f7.3)
      close(22)
      stop
      end

      subroutine addstar
c
c if Vmag present, add current star to end of starlist
c
      parameter (MAXSTARS=10000)
      parameter (XINDEF=99.999,XCHECK=99.0)
      real*8 ra(MAXSTARS), dec(MAXSTARS),rasq(MAXSTARS),
     $   decsq(MAXSTARS), magsum(MAXSTARS,2), magsq(MAXSTARS,2),
     $   rax,decx
      real*4  xmag1,xmag2,xmag3,dra,ddec,dr(MAXSTARS),dd(MAXSTARS)
      integer nr(MAXSTARS),np(MAXSTARS),n
      common /blka/ ra,dec,rasq,decsq,magsum,magsq,dr,dd,nr,np
      common /blkb/ rax,decx,xmag1,xmag2,xmag3,dra,ddec,n
c
        if (xmag1.gt.XCHECK) return
        n=n+1
        ra(n) = rax
        dec(n) = decx
        dd(n) = ddec
        dr(n) = dra
        rasq(n) = ra(n)*ra(n)
        decsq(n) = dec(n)*dec(n)
        magsum(n,1) = xmag1
        magsq(n,1) = xmag1*xmag1
        nr(n) = 1
        if (xmag2.lt.XCHECK) then
          np(n) = 1
          magsum(n,2) = xmag2
          magsq(n,2) = xmag2*xmag2
        else
          np(n) = 0
          magsum(n,2) = 0.
          magsq(n,2) = 0.
        endif
        if (xmag3.lt.XCHECK) then
          np(n) = np(n) + 1
          magsum(n,2) = magsum(n,2) + xmag3
          magsq(n,2) = magsq(n,2) + xmag3*xmag3
        endif
        return
        end

      SUBROUTINE SORT (n,x,index)
c
c heapsort of array x with corresponding index array index
c from numerical recipies
c
      REAL*8 x(n),xx
      INTEGER 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