program collate
c
c collate TASS starlists, output result
c requirement: all output objects must have been seen in 2 or
c more cameras.  Hardcoded for tm3get format (800x896) with
c 20 overlap rows, and for 10arcsec radius error box, and for
c a maximum of 10000 stars per camera, 3 files per camera.
c These parameters can all be changed.
c I run this program under DEC Ultrix f77.  I create an input
c script and redirect input to that script so that I can run this
c program on the same raw data without having to reenter all the
c filenames.
c written 01-April-1997 aah
c mods:
c  23-April-1997 for latest MG format
c
      parameter (MAXSTARS=10000,MAXFRAMES=3,MAXINDX=20000)
      real*8 ra(MAXSTARS,3), dec(MAXSTARS,3),hjd(MAXSTARS,3)
      real*4 mag(MAXSTARS,3),dmag(MAXSTARS,3),
     $       fwhmx(MAXSTARS,3),fwhmy(MAXSTARS,3)
      integer nstars(3),indx(MAXINDX,3),iflag,nidx
      common /blka/ ra,dec,hjd,mag,dmag,fwhmx,fwhmy,indx,nstars,nidx
c
      print *,'COLLATE Version 1.1'
      print *,' '
c
c MG's STAR has had 3 output formats so far. Sextractor is GG's version.
c
      print *,'Enter input filetype:'
      print *,'    (0=star1,1=star2,2=star3,3=sextractor): '
      read (5,*) iflag
c
c read all input files (hardcoded to max of 3 cameras)
c
      call getstars(1,iflag)
      call getstars(2,iflag)
      call getstars(3,iflag)
c
c now match up objects found on multiple cameras
c
      call match
c
c now create means and output results
c
      call writestars
c
      stop
      end

      subroutine getstars(n,iflag)
c
c read raw extracted files
c while you could figure out filetypes from reading the input
c file itself, I've hardcoded the formatting for now.
c
      parameter (MAXSTARS=10000,MAXFRAMES=3,MAXINDX=20000)
      parameter (XINDEF=999.999,YINDEF=99.999)
      real*8 ra(MAXSTARS,3), dec(MAXSTARS,3),hjd(MAXSTARS,3)
      real*8 hjdx
      real*4 mag(MAXSTARS,3),dmag(MAXSTARS,3),y,
     $       fwhmx(MAXSTARS,3),fwhmy(MAXSTARS,3),x
      integer nstars(3),i,indx(MAXINDX,3),n,i1,i2,iflag,j,nframe
      integer nidx
      character*80 fname(MAXFRAMES)
      common /blka/ ra,dec,hjd,mag,dmag,fwhmx,fwhmy,indx,nstars,nidx
c
c read in the file names as an 'xxx' terminated list
c
      i=1
100   continue
      write (6,900) n-1
900   format ('Enter camera',i1,' frame, xxx ends: ')
      read (5,'(a)') fname(i)
      if (fname(i).eq.'xxx') goto 200
      i = min((i+1),MAXFRAMES)
      goto 100
200   continue
      nframe = i - 1
      nstars(n) = 1
      if (nframe.gt.0) then
c
c read each requested frame
c
      do j=1,nframe
        open (unit=1,file=fname(j),status='old')
c
c first, strip the header (and calculate base JD from filename)
c
        if (iflag.lt.3) then
        do i=1,5
          read (1,*)
        enddo
        endif
        if (iflag.eq.3) read(1,*)
        i = nstars(n)
        read (fname(j),909) hjdx
909     format (4x,f8.3)
        hjdx = hjdx + 50000.0
300     continue
c
c type of read depends on who wrote the file.
c
          if(iflag.eq.0) then
            read(1,*,end=400) x,y,i1,i2,ra(i,n),dec(i,n),
     $      mag(i,n),dmag(i,n)
            hjd(i,n) = hjdx
            fwhmx(i,n) = 0.
            fwhmy(i,n) = 0.
          endif
          if(iflag.eq.1) read(1,901,end=400) x,y,ra(i,n),dec(i,n),
     $      mag(i,n),dmag(i,n),fwhmx(i,n),fwhmy(i,n),hjd(i,n)
901         format (10x,f6.2,2x,f6.2,13x,2f8.4,1x,2f6.3,10x,
     $       2f6.2,3x,f11.5)
          if(iflag.eq.2) read(1,902,end=400) x,y,ra(i,n),dec(i,n),
     $      mag(i,n),dmag(i,n),fwhmx(i,n),fwhmy(i,n),hjd(i,n)
902         format (10x,f6.2,2x,f6.2,14x,2f8.4,1x,2f6.3,10x,
     $       2f6.2,3x,f11.5)
          if (iflag.eq.3) then
            read(1,903,end=400) y,x,ra(i,n),dec(i,n),mag(i,n)
903         format(f8.2,1x,f8.2,e14.6,3x,e14.6,3x,f7.4)
            if (mag(i,n).gt.90.0) goto 300
            hjd(i,n) = hjdx
            dmag(i,n) = YINDEF
            fwhmx(i,n) = YINDEF
            fwhmy(i,n) = YINDEF
          endif
c two checks to remove overlap
          if (j.ne.1.and.y.lt.10.0) goto 300
          if (j.lt.nframe.and.nframe.gt.1.and.y.ge.886.0) goto 300
          i = min((i+1),MAXSTARS)
          goto 300
400     continue
        nstars(n) = i-1
        close(1)
      enddo
      write (6,904) n,nstars(n)
904   format('nstars(',i1,')= ',i5)
      endif
      return
      end

      subroutine match
c
c Heart of COLLATE.
c Match star lists using 10arcsec positional error.
c This subroutine is a kludge. Surely there is a better way!
c Output is indx(n,m) where m=1,2,3 indicating the index of
c  a corresponding star for cameras 0,1,2 respectively, and
c  n = running index of matched pairs and triplets
c
c A better way would be to sort camera1 & 2 by RA, then binary
c search in the loops so that you don't have to do N x M x 3
c searches.
c
      parameter (MAXSTARS=10000,MAXFRAMES=3,MAXINDX=20000)
      real*8 ra(MAXSTARS,3), dec(MAXSTARS,3),hjd(MAXSTARS,3),
     $       minerr,err
      real*4 mag(MAXSTARS,3),dmag(MAXSTARS,3),
     $       fwhmx(MAXSTARS,3),fwhmy(MAXSTARS,3)
      integer nstars(3),i,j,k,indx(MAXINDX,3),nidx
      integer ix(MAXSTARS,2)
      common /blka/ ra,dec,hjd,mag,dmag,fwhmx,fwhmy,indx,nstars,nidx
c
c minerr is error radius (10arcsec) in degrees, squared
c
      minerr = (10./3600.)**2
      do i=1,2
        do j=1,MAXSTARS
          ix(j,i) = 1
        enddo
      enddo
c
c first, include all camera0 objects
c  that are matched with a camera1 and/or camera2 object
c
      nidx = 1
      if (nstars(1).ne.0) then
        do i=1,nstars(1)
          indx(nidx,1) = i
          indx(nidx,2) = 0
          indx(nidx,3) = 0
c compare against camera1
          if (nstars(2).ne.0) then
            do j=1,nstars(2)
              if (ix(j,1).ne.0) then
                err = (ra(i,1)-ra(j,2))**2 + (dec(i,1)-dec(j,2))**2
                if (err.lt.minerr) then
                  indx(nidx,2) = j
                  ix (j,1) = 0
                  goto 10
                endif
              endif
            enddo
          endif
10        continue
c compare against camera2
          if (nstars(3).ne.0) then
            do j=1,nstars(3)
              if (ix(j,2).ne.0) then
                err = (ra(i,1)-ra(j,3))**2 + (dec(i,1)-dec(j,3))**2
                if (err.lt.minerr) then
                  indx(nidx,3) = j
                  ix (j,2) = 0
                  goto 20
                endif
              endif
            enddo
          endif
20        continue
          if (indx(nidx,2).ne.0.or.indx(nidx,3).ne.0) nidx=nidx+1
        enddo
      endif
c
c now, find all unmatched camera2 objects that are
c matched with camera3
c
      if (nstars(2).ne.0) then
        do i=1,nstars(2)
          if (ix(i,1).ne.0) then
            indx(nidx,1) = 0
            indx(nidx,2) = i
            indx(nidx,3) = 0
            if (nstars(3).ne.0) then
              do j=1,nstars(3)
                if(ix(j,2).ne.0) then
                  err = (ra(i,2)-ra(j,3))**2 + (dec(i,2)-dec(j,3))**2
                  if (err.lt.minerr) then
                    indx(nidx,3) = j
                    ix(j,2) = 0
                    nidx = nidx + 1
                    goto 30
                  endif
                endif
              enddo
            endif
          endif
30        continue
        enddo
      endif
      nidx = nidx - 1
      print *,'Found ',nidx,' matches'
      return
      end

      subroutine writestars
c
c create means for coordinates and magnitudes, output
c Note: I do not sort the output by RA.  The basic ordering
c is the same as in camera0's datafile, plus some additional
c stars from the camera1-camera2 matching.
c
      parameter (MAXSTARS=10000,MAXFRAMES=3,MAXINDX=20000)
      parameter (XINDEF=999.999,YINDEF=99.999)
      real*8 ra(MAXSTARS,3), dec(MAXSTARS,3),hjd(MAXSTARS,3),
     $   sumn,sumra,sumra2,sumdec,sumdec2,dra,ddec,ramean,decmean
      real*4 mag(MAXSTARS,3),dmag(MAXSTARS,3),
     $       fwhmx(MAXSTARS,3),fwhmy(MAXSTARS,3),x(5)
      integer nstars(3),indx(MAXINDX,3),i,j,k,m,np,nidx
      character*80 fout
      common /blka/ ra,dec,hjd,mag,dmag,fwhmx,fwhmy,indx,nstars,nidx
c
c perform some initialization stuff
c
      do i=1,5
        x(i) = YINDEF
      enddo
c
c get output file name
c
      print *,'Enter output filename: '
      read (5,'(a)') fout
      open (unit=2,file=fout,status='new')
c
c loop over number of matched pairs and triplets
c
      do i=1,nidx
        j = indx(i,1)
        k = indx(i,2)
        m = indx(i,3)
c
c form standard deviation partial sums
c
        sumra = 0.0
        sumdec = 0.0
        sumra2 = 0.0
        sumdec2 = 0.0
        np = 0
        if (j.ne.0) then
          sumra = sumra + ra(j,1)
          sumdec = sumdec + dec(j,1)
          sumra2 = sumra2 + ra(j,1)**2
          sumdec2 = sumdec2 + dec(j,1)**2
          np = np + 1
        endif
        if (k.ne.0) then
          sumra = sumra + ra(k,2)
          sumdec = sumdec + dec(k,2)
          sumra2 = sumra2 + ra(k,2)**2
          sumdec2 = sumdec2 + dec(k,2)**2
          np = np + 1
        endif
        if (m.ne.0) then
          sumra = sumra + ra(m,3)
          sumdec = sumdec + dec(m,3)
          sumra2 = sumra2 + ra(m,3)**2
          sumdec2 = sumdec2 + dec(m,3)**2
          np = np + 1
        endif
c
c calculate means and standard deviations
c
        ramean = sumra / float(np)
        decmean = sumdec / float(np)
        sumn = float(np)
        dra = dsqrt(abs((sumn*sumra2-sumra*sumra)/
     $        (sumn*(sumn-1.))))
        ddec = dsqrt(abs((sumn*sumdec2-sumdec*sumdec)/
     $         (sumn*(sumn-1.))))
c
c now write output file
c note: first line is camera1 (assumed = V),
c  second line is camera0, third line is camera2.
c
        if (k.ne.0) then
          write(2,901) ramean,decmean,dra,ddec,hjd(k,2),
     $    fwhmx(k,2),fwhmy(k,2),mag(k,2),dmag(k,2)
        else
          write(2,901) ramean,decmean,dra,ddec,x
        endif
        if (j.ne.0) then
          write (2,902) hjd(j,1),
     $    fwhmx(j,1),fwhmy(j,1),mag(j,1),dmag(j,1)
        else
          write (2,902) x
        endif
        if (m.ne.0) then
          write (2,902) hjd(m,3),
     $    fwhmx(m,3),fwhmy(m,3),mag(m,3),dmag(m,3)
        else
          write (2,902) x
        endif
901     format(4f9.4,f11.4,4f7.3)
902     format(5x,f11.4,4f7.3)
100     continue
      enddo
      close(2)
      return
      end