program transform
c
c find standard stars, transform, output file
c requires COLLATEd file with accurate astrometry
c plus a file of standard stars in a specific format.
c This program is quite convoluted because I'm trying
c to automate it as much as possible.
c written 01-April-1997 aah
c 05-may-97 fixed bug in getstds file read
c 28-may-97 fixed xfrm bug for n=1
c
      character*80 file1,file2
c
      print *,'TRANSFORM Version 1.1'
      print *,' '
      print *,'Enter standards filename: '
      read (5,'(a)') file1
      open (unit=2,file=file1,status='old')
      call getstds
      print *,'Enter mean coefficients filename: '
      read (5,'(a)') file2
c
c loop over input files
c
100   continue
      print *,'Enter input filename, xxx ends: '
      read (5,'(a)') file1
      if (file1.ne.'xxx') then
        open (unit=1,file=file1,status='old')
        print *,'Enter output filename: '
        read (5,'(a)') file1
        open (unit=3,file=file1,status='new')
        open (unit=22,file=file2,status='old')
        call getstars
        call calctrans
        call writestars
        goto 100
      endif
      stop
      end

      subroutine getstds
c
c read standards file
c note: maximum of 2000 standard stars allowed
c format is one line header, followed by
c name, gsc#, ra, dec, b, v, r, i
c where ra,dec are J2000 in degrees
c look at sample file to see exact columnar spacing.
c Note: for now, I'm deleting any star that is
c fainter than V=13.
c
      parameter (MAXSTD=2000)
      real*8 ra(MAXSTD),dec(MAXSTD)
      real*4 vstd(MAXSTD),vistd(MAXSTD),a,b,c,d
      integer i,nstd,indx(MAXSTD,2)
      common /blka/ ra,dec,vstd,vistd,indx,nstd
c
      read (2,*)
      i=1
100   continue
      read (2,900,end=200) ra(i),dec(i),a,b,c,d
900   format (20x,2f11.5,2x,4f7.3)
      if (b.gt.13.0) goto 100
      vstd(i) = b
      vistd(i) = b-d
      i=i+1
      if (i.gt.MAXSTD) goto 200
      goto 100
200   continue
      nstd = i-1
      write (6,901) nstd
901   format (i5,' Standards were read.')
      close(2)
      return
      end

      subroutine getstars
c
c read raw collated file, find matches with standards
c matches are 15arcsec radius error box
c can have 2 separate V-I solutions, so save v,vi1,vi2.
c
      parameter (MAXSTD=2000)
      real*8 ra(MAXSTD),dec(MAXSTD),rax,decx
      real*8 hjd1,hjd2,hjd3,err,minerr
      real*4 vstd(MAXSTD),vistd(MAXSTD),dra,ddec
      real*4 fwhmx1,fwhmx2,fwhmx3,fwhmy1,fwhmy2,fwhmy3
      real*4 x1,x2,x3,y1,y2,y3,dmag1,dmag2,dmag3
      real*4 mag1,mag2,mag3,v1(MAXSTD),vi1(MAXSTD),vi2(MAXSTD)
      real*4 v2(MAXSTD),vivi1(MAXSTD),vivi2(MAXSTD)
      integer i,nstd,n1,n2,indx(MAXSTD,2)
      common /blka/ ra,dec,vstd,vistd,indx,nstd
      common /blkb/ v1,v2,vi1,vi2,vivi1,vivi2,n1,n2
c
      minerr = (15./3600.)**2
      n1 = 1
      n2 = 1
c
c loop over stars in collated list
c
100   continue
        read (1,901,end=200) rax,decx,dra,ddec,
     $    hjd1,fwhmx1,fwhmy1,mag1,dmag1,
     $    hjd2,fwhmx2,fwhmy2,mag2,dmag2,
     $    hjd3,fwhmx3,fwhmy3,mag3,dmag3
901     format(4f9.4,f11.4,4f7.3/
     $    5x,f11.4,4f7.3/
     $    5x,f11.4,4f7.3)
c
c brute force check against standards
c for transformation, need V (mag1) + I (mag2 or mag3).  If
c this object does not have V & I, pass it over.
c
        if (mag1.gt.99.0) goto 100
        if (mag2.gt.99.0.and.mag3.gt.99.0) goto 100
c
c now check vs. standard star list to see if have coordinate match
c vx = vstd - vins
c vix = (V-I)std
c vivix = (V-I)std - (V-I)ins
c
        do i=1,nstd
          err = (ra(i)-rax)**2 + (dec(i)-decx)**2
          if (err.lt.minerr) then
            if (mag2.lt.99.0) then
              v1(n1) = vstd(i) - mag1
              vivi1(n1) = vistd(i) - (mag1-mag2)
              vi1(n1) = vistd(i)
              indx(i,1) = i
              n1 = n1 + 1
            endif
            if (mag3.lt.99.0) then
              v2(n2) = vstd(i) - mag1
              vivi2(n2) = vistd(i) - (mag1-mag3)
              vi2(n2) = vistd(i)
              indx(i,2) = i
              n2 = n2 + 1
            endif
            goto 100
          endif
        enddo
        goto 100
200   continue
      n1 = n1-1
      n2 = n2-1
      write (6,904) n1,n2
904   format (i5,' Stars with I1 and ',i5,
     $  ' Stars with I2 were read.')
      return
      end

      subroutine calctrans
c
c calculate the transformation coefficients
c this is the heart of the program.  I usually include
c interactive graphics here so that you can delete
c discrepant points, but we need to be as automated as
c possible for TASS.  Therefore, I just iterate once after
c removing the worst points.
c
      parameter (MAXSTD=2000)
      real*4 v1(MAXSTD),vi1(MAXSTD),vi2(MAXSTD)
      real*4 vivi1(MAXSTD),vivi2(MAXSTD),v2(MAXSTD)
      real*4 a1,a2,a3,a4,b1,b2,b3,b4,s1,s2,s3,s4
      integer n1,n2,n3,n4
      common /blkb/ v1,v2,vi1,vi2,vivi1,vivi2,n1,n2
      common /blkc/ a1,a2,a3,a4,b1,b2,b3,b4
c
      read(22,902) a3,b3,a1,b1
      read(22,902) a4,b4,a2,b2
902   format (4f7.3)
      close(22)
c
      print *,'Zerocalc(0), Fullcalc(1), or use Means(2)? '
      read (5,*) iflag
      if (iflag.eq.1) then
c
c perform full transformation calculation w/one iteration
c
        if (n1.ge.2) then
          call dofulltran(v1,vi1,vivi1,n1,a3,b3,a1,b1)
          call sigma(1,s3,s1)
          call dofulltran(v1,vi1,vivi1,n1,a3,b3,a1,b1)
        endif
        n3 = n1
        call sigma(1,s3,s1)
        if (n2.ge.2) then
          call dofulltran(v2,vi2,vivi2,n2,a4,b4,a2,b2)
          call sigma(2,s4,s2)
          call dofulltran(v2,vi2,vivi2,n2,a4,b4,a2,b2)
        endif
        n4 = n2
        call sigma(2,s4,s2)
        write (3,900) a3,b3,a1,b1,s3,s1,n3
        write (3,901) a4,b4,a2,b2,s4,s2,n4
900     format (' Full transform using I1: ',6f8.3,i5)
901     format (' Full transform using I2: ',6f8.3,i5)
      else
        if (iflag.eq.0) then
c
c just calculate zero points w/one iteration
c
          if (n1.ge.1) then
            call dozerotran(v1,vi1,vivi1,n1,a3,b3,a1,b1)
            call sigma(1,s3,s1)
            call dozerotran(v1,vi1,vivi1,n1,a3,b3,a1,b1)
          endif
          n3 = n1
          call sigma(1,s3,s1)
          if (n2.ge.1) then
            call dozerotran(v2,vi2,vivi2,n2,a4,b4,a2,b2)
            call sigma(2,s4,s2)
            call dozerotran(v2,vi2,vivi2,n2,a4,b4,a2,b2)
          endif
          n4 = n2
          call sigma(2,s4,s2)
          write (3,903) a3,b3,a1,b1,s3,s1,n3
          write (3,904) a4,b4,a2,b2,s4,s2,n4
903       format (' Zero transform using I1: ',6f8.3,i5)
904       format (' Zero transform using I2: ',6f8.3,i5)
        else
          write (3,905) a3,b3,a1,b1,s3,s1,n1
          write (3,906) a4,b4,a2,b2,s4,s2,n2
905       format (' Mean transform using I1: ',6f8.3,i5)
906       format (' Mean transform using I2: ',6f8.3,i5)
        endif
      endif
      return
      end

      subroutine dofulltran (x,y,z,n,a,b,c,d)
c
c main part of linear transformation routine
c if n=1, just return old values
c
      real*4 x(n),y(n),z(n)
      real*4 a,b,c,d
      integer n
c
      if (n.gt.1) then
c transform for v as function of (v-i)
        call lin (y,x,n,a,b)
c transform for (v-i) as function of (v-i)
        call lin (y,z,n,c,d)
        c = 1./(1.-c)
        d = c*d
      else
        c = 1. - (1./c)
        d = z(1)-c*y(1)
        b = x(1)-a*y(1)
        c = 1./(1.-c)
        d = c*d
      endif
      return
      end

      subroutine dozerotran (x,y,z,n,a,b,c,d)
c
c main part of zero transformation routine
c this uses mean slope and just determines intercepts b,d
c
      real*8 sum1,sum2
      real*4 x(n),y(n),z(n)
      real*4 a,b,c,d
      integer n
c
      if (n.gt.1) then
        sum1 = 0.
        sum2 = 0.
        c = 1. - (1./c)
        do i=1,n
          sum2 = sum2 + (x(i)-a*y(i))
          sum1 = sum1 + (z(i)-c*y(i))
        enddo
        d = sum1 / float(n)
        b = sum2 / float(n)
        c = 1./(1.-c)
        d = c*d
      else
        c = 1. - (1./c)
        d = z(1)-c*y(1)
        b = x(1)-a*y(1)
        c = 1./(1.-c)
        d = c*d
      endif
      return
      end

      subroutine sigma (m,s1,s2)
c
c calculate rms deviation for standards with this fit
c
      parameter (MAXSTD=2000,XINDEF=9.999)
      real*8 ra(MAXSTD),dec(MAXSTD),sum1,sum2
      real*4 vstd(MAXSTD),vistd(MAXSTD),s1,s2
      real*4 v1(MAXSTD),vi1(MAXSTD),vi2(MAXSTD)
      real*4 v2(MAXSTD),vivi1(MAXSTD),vivi2(MAXSTD)
      real*4 a1,a2,a3,a4,b1,b2,b3,b4,xx,yy,x1,xi1
      real*4 dx(MAXSTD)
      integer i,nstd,n1,n2,m,indx(MAXSTD,2)
      common /blka/ ra,dec,vstd,vistd,indx,nstd
      common /blkb/ v1,v2,vi1,vi2,vivi1,vivi2,n1,n2
      common /blkc/ a1,a2,a3,a4,b1,b2,b3,b4
c
      sum1 = 0.0
      sum2 = 0.0
c check whether we are working with camera0-1 or camera1-2
      if (m.eq.1) then
c if more than 1 point in fit, then calculate deviations
        if (n1.gt.1) then
c loop over total number of standards used
          do i=1,n1
            j = indx(i,m)
            xx = vistd(j) - vivi1(i)
            yy = vstd(j) - v1(i)
            xi1 = a1*(xx)+b1
            x1 = yy + a3*xi1 + b3
            dx(i) = abs(vstd(j)-x1)
c sum all errors
            sum1 = sum1 + dx(i)
            sum2 = sum2 + abs(vistd(j)-xi1)
          enddo
c find mean errors
          s1 = sum1 / float(n1)
          s2 = sum2 / float(n1)
c and remove the 2(max) most discrepant points for next iteration
          call remove(dx,v1,vi1,vivi1,indx(1,1),n1)
        else
          s1 = XINDEF
          s2 = XINDEF
        endif
c using camera1-2
      else
c if more than 1 point in fit, then calculate deviations
        if (n2.gt.1) then
c loop over total number of standards used
          do i=1,n2
            j = indx(i,m)
            xx = vistd(j) - vivi2(i)
            yy = vstd(j) - v2(i)
            xi1 = a2*(xx)+b2
            x1 = yy + a4*xi1 + b4
            dx(i) = abs(vstd(j)-x1)
c sum all errors
            sum1 = sum1 + dx(i)
            sum2 = sum2 + abs(vistd(j)-xi1)
           enddo
c find mean errors
          s1 = sum1 / float(n2)
          s2 = sum2 / float(n2)
c and remove the 1 or 2 most discrepant points for next iteration
          call remove(dx,v2,vi2,vivi2,indx(1,2),n2)
        else
          s1 = XINDEF
          s2 = XINDEF
        endif
      endif
      return
      end

      subroutine remove (dx,v,vi,vivi,indx,n)
c
c remove discrepant points from standards
c up to two points are removed.
c This is convoluted logic.  I'm sure there is a better way
c to remove discrepant points from arrays.
c
      parameter (MAXSTD=2000)
      real*4 dx(*),v(*),vi(*),vivi(*)
      integer indx(*),id(MAXSTD),idx(MAXSTD),n
c
c first, check to see if sufficient points to make remove effective
c
      if (n.ge.5) then
c
c sort dx
c
        do i=1,n
          id(i) = i
        enddo
        call sort (n,dx,id)
c
c rearrange data arrays according to sort
c note: reuse dx since no longer needed
c
        do i=1,n
          dx(i) = v(id(i))
        enddo
        do i=1,n
          v(i) = dx(i)
          dx(i) = vi(id(i))
        enddo
        do i=1,n
          vi(i) = dx(i)
          dx(i) = vivi(id(i))
        enddo
        do i=1,n
          vivi(i) = dx(i)
          idx(i) = indx(id(i))
        enddo
        do i=1,n
          indx(i) = idx(i)
        enddo
c this drops last two values
        n = n - 2
      endif
      return
      end

      subroutine writestars
c
c this routine calculates V & I and then writes updated array to file
c
      parameter (XINDEF=99.999)
      real*8 rax,decx,hjd1,hjd2,hjd3
      real*4 dra,ddec
      real*4 fwhmx1,fwhmx2,fwhmx3,fwhmy1,fwhmy2,fwhmy3
      real*4 x1,x2,x3,y1,y2,y3,dmag1,dmag2,dmag3
      real*4 mag1,mag2,mag3,xx,v1,v2,vi1,vi2
      real*4 a1,a2,a3,a4,b1,b2,b3,b4
      common /blkc/ a1,a2,a3,a4,b1,b2,b3,b4
c
c rewind input file and loop over number of entries
c
      rewind(1)
100   continue
        read (1,901,end=200) rax,decx,dra,ddec,
     $    hjd1,fwhmx1,fwhmy1,mag1,dmag1,
     $    hjd2,fwhmx2,fwhmy2,mag2,dmag2,
     $    hjd3,fwhmx3,fwhmy3,mag3,dmag3
901     format(4f9.4,f11.4,4f7.3/
     $    5x,f11.4,4f7.3/
     $    5x,f11.4,4f7.3)
        if (mag1.gt.99.0) then
          vi1 = mag2
          vi2 = mag3
          xx = mag1
        else
          if (mag2.gt.99.0) then
            vi1 = XINDEF
            v1 = XINDEF
          else
            vi1 = a1*(mag1-mag2)+b1
            v1 = mag1 + a3*vi1 + b3
            vi1 = v1 - vi1
          endif
          if (mag3.gt.99.0) then
            vi2 = XINDEF
            v2 = XINDEF
          else
            vi2 = a2*(mag1-mag3)+b2
            v2 = mag1 + a4*vi2 + b4
            vi2 = v2 - vi2
          endif
          if (v1.eq.XINDEF) xx = v2
          if (v2.eq.XINDEF) xx = v1
          if (v1.ne.XINDEF.and.v2.ne.XINDEF) xx=(v1+v2)/2.
        endif
        write (3,901) rax,decx,dra,ddec,
     $    hjd1,fwhmx1,fwhmy1,xx,dmag1,
     $    hjd2,fwhmx2,fwhmy2,vi1,dmag2,
     $    hjd3,fwhmx3,fwhmy3,vi2,dmag3
      goto 100
200   continue
      close(1)
      close(3)
      return
      end

	subroutine lin (x,y,m,a,b)
c
c  *** Linear least squares routine from Nielson ***
c  Simple solution with no weighting factors
c
c  Inputs:
c		x,y   data arrays
c		m     number of points in arrays
c  Outputs:
c		a,b   where y=ax+b
c
c  *** written by A. Henden 1973 ***
c
	dimension x(m),y(m)
	double precision a1,a2,a3,c1,c2,det
c  initialize summing parameters
	a2=0.
	a3=0.
	c1=0.
	c2=0.
	a1=m
c  loop to set up matrix coefficients
	do 10 i=1,m
	a2=a2+x(i)
	a3=a3+x(i)*x(i)
	c1=c1+y(i)
	c2=c2+y(i)*x(i)
10	continue
c  solve matrix - simple since only 2x2
	det=1./(a1*a3-a2*a2)
	b=-(a2*c2-c1*a3)*det
	a=(a1*c2-c1*a2)*det
	return
	end

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