PROGRAM td2uv
c************************************************************************
c
c  PROGRAM NAME: td2uv
c
c  SOURCE FILE:  td2uv.f
c
c  PURPOSE:      Simple interface program between Hipparcos Transit 
c                data (TD) and aperture synthesis imaging software
c                such as difmap (Caltech)
c
c  DESCRIPTION:  This program will, for a given HIP identifier, 
c                extract Hipparcos Transit Data (TD) from CD-ROM #6 
c                of the Hipparcos and Tycho Catalogues (ESA SP-1200) 
c                and create an UV-FITS file that can be used as input
c                for the Caltech DIFMAP program (or other software
c                packages) for creation and deconvolution of images 
c                using aperture synthesis methods.
c
c                The output file is called tdXXXXXX.uvf, where XXXXXX 
c                is the (zero-padded) HIP number.
c
c                The Hipparcos Transit data, though not a real 
c                interferometer data, can be used to create synthesised
c                images by methods similar to those in radio astronomy.  
c                Many software packages for radio astronomy use UV-FITS
c                as the standard input format. This was the motivation 
c                behind this conversion program.
c
c  AUTHORS:      Carl Fredrik Quist and Lennart Lindegren
c                (fredrikq@astro.lu.se, lennart@astro.lu.se)
c                Lund Observatory
c
c  DATE:         1999 January 29 (Version 1.0)
c                2000 October 26 Scisoft ESO version with UNIX as default
c                                file locations.
c
c***********************************************************************
c
c  IMPORTANT NOTES:
c
c  1.  Before compilation it is necessary to set the character
c      variables filedat and fileidx to the correct path+filenames
c      of the TD files on the CD-ROM (see examples below for UNIX
c      and DOS/Windows)
c
c  2.  The program automatically corrects the six records in hipj.dat
c      that contain fields with numerical overfow (asterisks).
c
c***********************************************************************c
c  EXPLICIT DECLARATION OF ALL VARIABLES:
c
      IMPLICIT NONE
c
c  General constants:
c
      REAL*8 pi
      REAL*8 clight
      PARAMETER (pi = 3.14159265358979324d0, clight = 299792458d0)
c
c  Names for the transit data file, transit data index file,
c  and UVFITS file; and associated logical units:
c
      CHARACTER filedat*100, fileidx*100, fileuv*12
      INTEGER ludat, luidx, luuv
c
c  Variables for the transit data index -
c  idx = hip2idx(hip) is the record number in filedat for entry hip:
c
      INTEGER hip, hipmax, idx
      PARAMETER (hipmax = 120416)
      INTEGER hip2idx(hipmax)
c
c  Character string for reading records of the transit data file:
c
      CHARACTER indat*127
c
c  Variables for reading transit data (header record):
c
      INTEGER nhip(3), np, nt
      REAL*8 a0, d0
      REAL*8 par0, pma0, pmd0, col(3)
c
c  Variables for reading transit data (pointing record):
c
      INTEGER pt(3,9)
c
c  Variables for reading transit data (transit records):
c
      INTEGER ip, fx, fy, fp, flag
      REAL*8 time
      REAL*8 lnb1, bnorm(2:5), lnsig(5), s1, s2, sigatt
c
c  Variables for the new reference point, shift in RA, DEC and PAR,
c  and for the decoded and shifted transit data:
c
      REAL*8 a1, d1
      REAL*8 pma1, pmd1, par1
      REAL*8 shfta, shftd, shftpar, phi
      REAL*8 b(5), sig(5), b1(5)
c
c  Buffer to hold all visibility data for one entry
c  and logical array for their inclusion on output file:
c
      INTEGER ntmax
      PARAMETER (ntmax = 583)
      REAL*4 buf4(9,3*ntmax)
      LOGICAL incl(3*ntmax)
c
c  Buffer and counter for output of real values to FITS file:
c
      REAL*4 values(720)
      INTEGER nval
c
c  Counters for transits, visibilities, output records:
c
      INTEGER it, nvis, irec
c
c  General counters and other variables:
c
      INTEGER i, j, k
      REAL*8 m1bar, m2bar, time0, time1, time2, freq
      CHARACTER err1*3,err2*3,err3*5
      LOGICAL mustswap
c
c  END OF DECLARATIONS
c
c=======================================================================
c
c  Paths to transit data and transit data index file
c  (Hipparcos and Tycho Catalogues CD-ROM #6)
c  - un/comment and edit as appropriate:
c
c  Typical filenames for a UNIX installation:
c
      filedat = '/cdrom/cats/hip_j.dat'
      fileidx = '/cdrom/cats/hip_j.idx'
c
c  Typical filenames for a PC (DOS or Windows) installation:
c
c     filedat = 'E:\cats\hip_j.dat'
c     fileidx = 'E:\cats\hip_j.idx'
c
c=======================================================================
c
c  OTHER GENERAL SPECIFICATIONS:
c
c  Logical units:
c
      ludat = 11
      luidx = 12
      luuv  = 13
c
c  m1bar, m2bar are the theoretical modulation coefficients
c  for a point source:
c
      m1bar = 0.7100d0
      m2bar = 0.2485d0
c
c  UVFITS requires that a frequency (in Hz) is assigned to the
c  observation.  It is computed from the approximate effective
c  wavelength of the Hipparcos instrument, 550 nm:
c
      freq = clight/550d-9
c
c  Test IEEE compliance for binary data storage:
c
      IF (mustswap()) THEN
          WRITE(6,*) 'NOTE: byte swapping will be applied when ',
     +               'writing binary data to UVFITS file'
      ELSE
          WRITE(6,*) 'NOTE: binary storage complies with IEEE - ',
     +               'no byte swapping required'
      ENDIF
c
c  Open and read the transit data index file:
c
      OPEN(UNIT = luidx,
     +     FILE = fileidx,
     +     FORM = 'FORMATTED',
     +     STATUS = 'OLD')
      DO i = 1, hipmax
        READ(luidx,'(i7,1x)') hip2idx(i)
      ENDDO
      CLOSE(luidx)
c
c  Open the transit data file:
c
      OPEN(UNIT = ludat,
     +     FILE = filedat,
     +     ACCESS = 'DIRECT',
     +     RECL = 127,
     +     FORM = 'UNFORMATTED',
     +     STATUS = 'OLD')
c
c  Input a HIP number and construct the output file name (fileuv):
c
  100 WRITE(*,'(/1x,a)') 'Enter a HIP number (or 0 or EOF to quit):'
      READ(*,*,END=200) hip
      IF (hip .LT. 1 .OR. hip .GT. hipmax) GOTO 200
      idx = hip2idx(hip)
      IF (idx .EQ. -1) THEN
          WRITE(6,'(1x,a,i6)') 'no Transit Data available for HIP ',hip
          GOTO 100
      ENDIF
      WRITE(fileuv,'(a2,i6.6,a4)') 'td', hip, '.uvf'
c
c  Read the header record of the transit data:
c
      READ(ludat,REC=idx) indat
      READ(indat,'(3(i6,1x),i2,1x,i3,1x,2(f12.8,1x),f6.2,1x,
     +    2(f8.2,1x),2(f7.3,1x),f7.3)') (nhip(i),i=1,3),
     +    np, nt, a0, d0, par0, pma0, pmd0, (col(i),i=1,3)
c
c  Define the new reference point (centre of map) --- default is
c  equal to the reference point in the transit data header (= Input
c  Catalogue position, parallax, proper motion):
c
      a1 = a0
      d1 = d0
      par1 = par0
      pma1 = pma0
      pmd1 = pmd0
c
      WRITE(*,'(1x,a,i6)') 'Reference point for HIP entry ',hip
      DO  i = 1, 3
          IF (nhip(i) .NE. 0 .AND. nhip(i) .NE. hip) THEN
              WRITE(*,'(1x,a,i6,a)')
     +            ' (and associated entry ',nhip(i),')'
          ENDIF
      ENDDO
      WRITE(*,'(1x,a,f12.8)') 'right ascension (ICRS, deg) = ',a1
      WRITE(*,'(1x,a,f12.8)') 'declination (ICRS, deg)     = ',d1
      WRITE(*,'(1x,a,f7.2)')  'parallax (mas)              = ',par1
      WRITE(*,'(1x,a,f8.2)')  'PM in RA (ICRS, mas/yr)     = ',pma1
      WRITE(*,'(1x,a,f8.2)')  'PM in Dec (ICRS, mas/yr)    = ',pmd1
c
c  Read the pointing record of the transit data
c  (these data are actually not used):
c
      READ(ludat,REC=idx+1) indat
      IF(HIP .eq. 46586)THEN
        READ(indat,'(i1,1x,i3,1x,a3,1x,8(i1,1x,i3,1x,i3,1x))')
     ,           pt(1,1),pt(2,1),err1,((pt(i,j),i=1,3),j=2,9)
        pt(3,1) = 1308
      ELSEIF (HIP .EQ. 99819) THEN
        READ(indat,'(i1,1x,i3,1x,i3,1x,i1,1x,i3,1x,a3,1x,
     ,       7(i1,1x,i3,1x,i3,1x))')pt(1,1),pt(2,1),pt(3,1),
     ,       pt(1,2),pt(2,2),err1,((pt(i,j),i=1,3),j=3,9)
        pt(3,2) = -105
      ELSEIF (HIP .EQ. 101043) THEN
        READ(indat,'(i1,1x,a3,1x,i3,1x,8(i1,1x,i3,1x,i3,1x))')
     ,           pt(1,1),err1,pt(3,1),((pt(i,j),i=1,3),j=2,9)
        pt(2,1) = -157
      ELSEIF (HIP .EQ. 116191) THEN
        READ(indat,'(i1,1x,a3,1x,i3,1x,8(i1,1x,i3,1x,i3,1x))')
     ,           pt(1,1),err1,pt(3,1),((pt(i,j),i=1,3),j=2,9)
        pt(2,1) = -204
      ELSEIF (HIP .EQ. 117011) THEN
        READ(indat,'(i1,1x,i3,1x,a3,1x,i1,1x,i3,1x,a3,1x,
     ,     7(i1,1x,i3,1x,i3,1x))') pt(1,1),pt(2,1),err1,pt(1,2),
     ,     pt(2,2),err2,((pt(i,j),i=1,3),j=3,9)
        pt(3,1) = -234
        pt(3,2) = -234
      ELSE
        READ(indat,'(9(i1,1x,i3,1x,i3,1x))') ((pt(i,j),i=1,3),j=1,9)
      ENDIF
c
c  k = counter for visibilities (all of them, included or not)
c
      k = 0
      DO it = 1, nt
          READ(ludat,REC=idx+1+it) indat
          IF((HIP .EQ. 48036).and.(it .EQ. 22))THEN
            READ(indat,'(i1,1x,f10.7,1x,3(i8,1x),f6.3,1x,4(f7.4,1x),
     +        3(f5.2,1x),a5,1x,f5.2,1x,2(f4.2,1x),f4.1,1x,i1)')
     +        ip, time, fx, fy, fp, lnb1, (bnorm(i),i=2,5),
     +        (lnsig(i),i=1,3),err3,lnsig(5), s1, s2, sigatt, flag
            lnsig(4)=-10.0d0
          ELSE
            READ(indat,'(i1,1x,f10.7,1x,3(i8,1x),f6.3,1x,4(f7.4,1x),
     +        5(f5.2,1x),2(f4.2,1x),f4.1,1x,i1)')
     +        ip, time, fx, fy, fp, lnb1, (bnorm(i),i=2,5),
     +        (lnsig(i),i=1,5), s1, s2, sigatt, flag
          ENDIF
c
c  Decode the signal parameters b() and standard deviations sig():
c
          b(1) = dexp(lnb1)
          DO  i = 2, 5
              b(i) = bnorm(i)*b(1)
          ENDDO
          DO  i = 1, 5
              sig(i) = dexp(lnsig(i))
          ENDDO
c
c  Shift the reference point from (a0,d1,par0,pma0,pmd0) to
c  (a1,d1,par1,pma1,pmd1) --- shfta, shftd, shftpar are the shifts
c  on the sky at the time of observation, expressed in degrees, while
c  phi is the corresponding phase shift in rad; b1() are the signal
c  parameters referred to the new reference point:
c
          shfta = (a1-a0)*dcos(d0*pi/180d0) + time*(pma1-pma0)/3600d3
          shftd = (d1-d0) + time*(pmd1-pmd0)/3600d3
          shftpar = (par1-par0)/3600d3
          phi = (dble(fx)*shfta + dble(fy)*shftd + dble(fp)*shftpar)
     +          *pi/180d0
          b1(1) = b(1)
          b1(2) = b(2)*dcos(phi) - b(3)*dsin(phi)
          b1(3) = b(2)*dsin(phi) + b(3)*dcos(phi)
          b1(4) = b(4)*dcos(2d0*phi) - b(5)*dsin(2d0*phi)
          b1(5) = b(4)*dsin(2d0*phi) + b(5)*dcos(2d0*phi)
c
c  time1, time2 = whole and fractional Julian Date of observation:
c
          time0 = 2447162.0d0 + (3.25d0 + time)*365.25d0
          time1 = dint(time0)
          time2 = time0 - time1
c
c---------------------------------------------------------------------
c  NOTE:  In the subsequent code, all visibility data (first, second
c         and zeroth harmonic for every transit) are included for
c         output.  Specific transits and/or harmonic components may
c         be excluded from output by setting incl(k) = .FALSE. for
c         the corresponding visibility data.
c---------------------------------------------------------------------
c
c  Visibility data for the first harmonic (fundamental frequency):
c
          k = k + 1
          incl(k) = .TRUE.
          buf4(1,k) = sngl(dble(fx)/(2d0*pi))
          buf4(2,k) = sngl(dble(fy)/(2d0*pi))
          buf4(3,k) = 0.0
          buf4(4,k) = 258.0
          buf4(5,k) = sngl(time1)
          buf4(6,k) = sngl(time2)
          buf4(7,k) = sngl(+b1(2)/m1bar)
          buf4(8,k) = sngl(-b1(3)/m1bar)
          buf4(9,k) = 1.0
c
c  Visibility data for the second harmonic (overtone):
c
          k = k + 1
          incl(k) = .TRUE.
          buf4(1,k) = sngl(dble(fx)/pi)
          buf4(2,k) = sngl(dble(fy)/pi)
          buf4(3,k) = 0.0
          buf4(4,k) = 772.0
          buf4(5,k) = sngl(time1)
          buf4(6,k) = sngl(time2)
          buf4(7,k) = sngl(+b1(4)/m2bar)
          buf4(8,k) = sngl(-b1(5)/m2bar)
          buf4(9,k) = 1.0
c
c  Visibility data for the DC component (zero frequency):
c
          k = k + 1
          incl(k) = .TRUE.
          buf4(1,k) = 0.0
          buf4(2,k) = 0.0
          buf4(3,k) = 0.0
          buf4(4,k) = 1286.0
          buf4(5,k) = sngl(time1)
          buf4(6,k) = sngl(time2)
          buf4(7,k) = sngl(b1(1))
          buf4(8,k) = 0.0
          buf4(9,k) = 1.0
      ENDDO
c
c  Count the number of visibilities actually included for output:
c
      nvis = 0
      DO k = 1, 3*nt
          IF (incl(k)) THEN
              nvis = nvis + 1
          ENDIF
      ENDDO
c
c  Open the output file (fileuv):
c
      OPEN(UNIT = luuv,
     +     FILE = fileuv,
     +     FORM = 'UNFORMATTED',
     +     ACCESS = 'DIRECT',
     +     RECL = 80*36,
     +     STATUS = 'UNKNOWN')
c
c  Output UVFITS header:
c
      irec = 0
      CALL uvf_hdr(luuv, hip, a1, d1, par1, pma1, pmd1,
     +             freq, nvis,   irec)
c
c  Output UVFITS data for all visibilities included:
c
      nval = 0
      DO k = 1, 3*nt
          IF (incl(k)) THEN
              DO j = 1, 9
                  nval = nval + 1
                  values(nval) = buf4(j,k)
                  IF (nval .EQ. 720) THEN
                      CALL uvf_dat(luuv, values, nval, irec)
                  ENDIF
              ENDDO
          ENDIF
      ENDDO
      CALL uvf_dat(luuv, values, nval, irec)
c
c  Output UVFITS extension data:
c
      CALL uvf_ext(luuv, freq,   irec)
c
c  Close output file and ask for another HIP number:
c
      CLOSE(luuv)
      WRITE(*,'(1x,a,a)')  'File created: ', fileuv
      WRITE(*,'(1x,a,i4)') 'Number of visibilities: ', nvis
      GOTO 100
c
c  Logical end of program:
c
  200 CONTINUE
      CLOSE(ludat)
        END
c***********************************************************************
      SUBROUTINE uvf_dat(lu, values, nval, irec)
c***********************************************************************
c
c  Write the UVFITS data record on logical unit lu.
c
c  INPUT:  lu     = logical unit opened to UV FITS file
c          values = array of up to 720 real*4 values
c          nval   = number of defined values in values
c                   ** MODIFIED ON OUTPUT **
c          irec   = record number for previous record
c                   ** MODIFIED ON OUTPUT **
c
c  OUTPUT: nval   = 0
c          irec   = record number for the last record written
c
      IMPLICIT NONE
      INTEGER lu, nval, irec, i
      REAL*4 values(720)
      LOGICAL mustswap
      IF (nval .GT. 0) THEN
          DO  i = nval+1, 720
              values(i) = 0.
          ENDDO
          IF (mustswap()) THEN
              CALL swapr4(values, 720)
          ENDIF
          irec = irec + 1
          WRITE(lu, REC=irec) values
          nval = 0
      ENDIF
      RETURN
      END
c***********************************************************************
      SUBROUTINE uvf_hdr(lu, hip, radeg, decdeg, parmas, pmamas,
     +                    pmdmas, freq, nvis,   irec)
c***********************************************************************
c
c  Write the UV FITS header on logical unit lu.
c
c  INPUT:  lu     = logical unit opened to UV FITS file
c          hip    = HIP identifier for the object
c          radeg  = reference Right Ascension in degrees
c          decdeg = reference Declination in degrees
c          freq   = reference frequency in Hz
c          nvis   = number of complex visibilities included
c          irec   = record number for previous record (should be = 0)
c                   ** MODIFIED ON OUTPUT **
c
c  OUTPUT: irec   = record number for the last record in header
c
      IMPLICIT NONE
      INTEGER lu, hip, nvis, irec, j
      REAL*8 radeg, decdeg, parmas, pmamas, pmdmas
      REAL*8 freq, scale, zero
      CHARACTER card(36)*80, key*8, val*20, com*50, lng*70,
     +          q*1, today*9
c
c  The character q is a quote ('):
c
      q = char(39)
c
c  Scale and zero point for real values:
c
      scale = 1.0
      zero = 0.0
c
c  Card formats: 1 = normal, 2 = long value, 3 = blank
c
1     FORMAT(a8,'= ',a20,a50)
2     FORMAT(a8,'= ',a70)
3     FORMAT(80x)
4     FORMAT(a8,72x)
5     FORMAT(a8,'  ',a20,a50)
c---------------------------------------------------------------
      key = 'SIMPLE  '
      val = '                   T'
      com = ' /THIS FOLLOWS STANDARD UVFITS FORMAT'
      WRITE(card(1),1) key, val, com
c---------------------------------------------------------------
      key = 'BITPIX  '
      val = '                 -32'
      com = ' /SINGLE PRECISION FLOATING POINT'
      WRITE(card(2),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS   '
      val = '                   6'
      com = ' /'
      WRITE(card(3),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS1  '
      val = '                   0'
      com = ' /ZERO VALUE IMPLIES UV FITS FILE'
      WRITE(card(4),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS2  '
      val = '                   3'
      com = ' /'
      WRITE(card(5),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS3  '
      val = '                   1'
      com = ' /'
      WRITE(card(6),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS4  '
      val = '                   1'
      com = ' /'
      WRITE(card(7),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS5  '
      val = '                   1'
      com = ' /'
      WRITE(card(8),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS6  '
      val = '                   1'
      com = ' /'
      WRITE(card(9),1) key, val, com
c---------------------------------------------------------------
      key = 'EXTEND  '
      val = '                   T'
      com = ' /THIS MIGHT BE AN ANTENNA FILE'
      WRITE(card(10),1) key, val, com
c---------------------------------------------------------------
      key = 'DATE    '
      WRITE(val,'(a,a,a)') q, today(), q
      com = ' /CREATION DATE FOR UV FITS FILE **NON-STANDARD'
      WRITE(card(11),1) key, val, com
c---------------------------------------------------------------
      key = 'ORIGIN  '
      val = q//'tran2uv.f'//q
      com = ' /NAME OF CONVERSION PROGRAM'
      WRITE(card(12),1) key, val, com
c---------------------------------------------------------------
      key = 'AUTHOR  '
      lng = q//'C.F. Quist & L. Lindegren'//q
      WRITE(card(13),2) key, lng
c---------------------------------------------------------------
      key = 'REFERENC'
      lng = q//'fredrikq@astro.lu.se & lennart@astro.lu.se'//q
      WRITE(card(14),2) key, lng
c---------------------------------------------------------------
      key = 'TELESCOP'
      val = q//'Hipparcos'//q
      com = ' /ESA ASTROMETRY SATELLITE'
      WRITE(card(15),1) key, val, com
c---------------------------------------------------------------
      key = 'OBJECT  '
      WRITE(val,'(a1,i18,a1)') q, hip, q
      com = ' /HIPPARCOS CATALOGUE NUMBER'
      WRITE(card(16),1) key, val, com
c---------------------------------------------------------------
      key = 'EQUINOX '
      val = '                2000'
      com = ' /ICRS'
      WRITE(card(17),1) key, val, com
c---------------------------------------------------------------
      key = 'BSCALE  '
      WRITE(val,'(e20.8)') scale
      com = ' /SCALING FACTOR FOR ARRAY VALUES'
      WRITE(card(18),1) key, val, com
c---------------------------------------------------------------
      key = 'BZERO   '
      WRITE(val,'(e20.8)') zero
      com = ' /OFFSET:  REAL = PIXVAL*BSCALE + BZERO'
      WRITE(card(19),1) key, val, com
c---------------------------------------------------------------
      key = 'BUNIT   '
      val = q//'CU      '//q//'          '
      com = ' /COUNTING UNITS'
      WRITE(card(20),1) key, val, com
c---------------------------------------------------------------
      key = 'OBSRA   '
      WRITE(val,'(e20.8)') radeg
      com = ' /REFERENCE RIGHT ASCENSION IN DEGREES'
      WRITE(card(21),1) key, val, com
c---------------------------------------------------------------
      key = 'OBSDEC  '
      WRITE(val,'(e20.8)') decdeg
      com = ' /REFERENCE DECLINATION IN DEGREES'
      WRITE(card(22),1) key, val, com
c---------------------------------------------------------------
      key = '        '
      WRITE(val,'(e20.8)') parmas
      com = ' /REFERENCE PARALLAX IN mas '
      WRITE(card(23),5) key, val, com
c---------------------------------------------------------------
      key = '          '
      WRITE(val,'(e20.8)') pmamas
      com = ' /REFERENCE PROPER MOTION IN RA IN mas/year'
      WRITE(card(24),5) key, val, com
c---------------------------------------------------------------
      key = '          '
      WRITE(val,'(e20.8)') pmdmas
      com = ' /REFERENCE PROPER MOTION DEC IN mas/year'
      WRITE(card(25),5) key, val, com
c---------------------------------------------------------------
      key = 'CTYPE2  '
      val = q//'COMPLEX '//q//'          '
      com = ' /'
      WRITE(card(26),1) key, val, com
c---------------------------------------------------------------
      key = 'CRVAL2  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(27),1) key, val, com
c---------------------------------------------------------------
      key = 'CDELT2  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(28),1) key, val, com
c---------------------------------------------------------------
      key = 'CRPIX2  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(29),1) key, val, com
c---------------------------------------------------------------
      key = 'CROTA2  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(30),1) key, val, com
c---------------------------------------------------------------
      key = 'CTYPE3  '
      val = q//'STOKES  '//q//'          '
      com = ' /'
      WRITE(card(31),1) key, val, com
c---------------------------------------------------------------
      key = 'CRVAL3  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(32),1) key, val, com
c---------------------------------------------------------------
      key = 'CDELT3  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(33),1) key, val, com
c---------------------------------------------------------------
      key = 'CRPIX3  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(34),1) key, val, com
c---------------------------------------------------------------
      key = 'CROTA3  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(35),1) key, val, com
c---------------------------------------------------------------
      key = 'CTYPE4  '
      val = q//'FREQ    '//q//'          '
      com = ' /'
      WRITE(card(36),1) key, val, com
c
      irec = irec + 1
      WRITE(lu,rec=irec) card
c---------------------------------------------------------------
      key = 'CRVAL4  '
      WRITE(val,'(e20.8)') freq
      com = ' /'
      WRITE(card(1),1) key, val, com
c---------------------------------------------------------------
      key = 'CDELT4  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(2),1) key, val, com
c---------------------------------------------------------------
      key = 'CRPIX4  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(3),1) key, val, com
c---------------------------------------------------------------
      key = 'CROTA4  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(4),1) key, val, com
c---------------------------------------------------------------
      key = 'CTYPE5  '
      val = q//'RA      '//q//'          '
      com = ' '
      WRITE(card(5),1) key, val, com
c---------------------------------------------------------------
      key = 'CRVAL5  '
      WRITE(val,'(e20.8)') radeg
      com = ' /'
      WRITE(card(6),1) key, val, com
c---------------------------------------------------------------
      key = 'CDELT5  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(7),1) key, val, com
c---------------------------------------------------------------
      key = 'CRPIX5  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(8),1) key, val, com
c---------------------------------------------------------------
      key = 'CROTA5  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(9),1) key, val, com
c---------------------------------------------------------------
      key = 'CTYPE6  '
      val = q//'DEC     '//q//'          '
      com = ' '
      WRITE(card(10),1) key, val, com
c---------------------------------------------------------------
      key = 'CRVAL6  '
      WRITE(val,'(e20.8)') decdeg
      com = ' /'
      WRITE(card(11),1) key, val, com
c---------------------------------------------------------------
      key = 'CDELT6  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(12),1) key, val, com
c---------------------------------------------------------------
      key = 'CRPIX6  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(13),1) key, val, com
c---------------------------------------------------------------
      key = 'CROTA6  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(14),1) key, val, com
c---------------------------------------------------------------
      key = 'GROUPS  '
      val = '                   T'
      com = ' /'
      WRITE(card(15),1) key, val, com
c---------------------------------------------------------------
      key = 'GCOUNT  '
      WRITE(val,'(i20)') nvis
      com = ' /NUMBER OF COMPLEX VISIBILITIES'
      WRITE(card(16),1) key, val, com
c---------------------------------------------------------------
      key = 'PCOUNT  '
      val = '                   6'
      com = ' /THERE ARE SIX PARAMETERS NEEDED'
      WRITE(card(17),1) key, val, com
c---------------------------------------------------------------
      key = 'PTYPE1  '
      val = q//'UU      '//q//'          '
      com = ' '
      WRITE(card(18),1) key, val, com
c---------------------------------------------------------------
      key = 'PSCAL1  '
      WRITE(val,'(e20.8)') (1.0d0/freq)
      com = ' /'
      WRITE(card(19),1) key, val, com
c---------------------------------------------------------------
      key = 'PZERO1  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(20),1) key, val, com
c---------------------------------------------------------------
      key = 'PTYPE2  '
      val = q//'VV      '//q//'          '
      com = ' '
      WRITE(card(21),1) key, val, com
c---------------------------------------------------------------
      key = 'PSCAL2  '
      WRITE(val,'(e20.8)') (1.0d0)/freq
      com = ' /'
      WRITE(card(22),1) key, val, com
c---------------------------------------------------------------
      key = 'PZERO2  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(23),1) key, val, com
c---------------------------------------------------------------
      key = 'PTYPE3  '
      val = q//'WW      '//q//'          '
      com = ' '
      WRITE(card(24),1) key, val, com
C---------------------------------------------------------------
      key = 'PSCAL3  '
      WRITE(val,'(e20.8)') (1.0d0)/freq
      com = ' /'
      WRITE(card(25),1) key, val, com
c---------------------------------------------------------------
      key = 'PZERO3  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(26),1) key, val, com
c---------------------------------------------------------------
      key = 'PTYPE4  '
      val = q//'BASELINE'//q//'          '
      com = ' '
      WRITE(card(27),1) key, val, com
c---------------------------------------------------------------
      key = 'PSCAL4  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(28),1) key, val, com
c---------------------------------------------------------------
      key = 'PZERO4  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(29),1) key, val, com
c---------------------------------------------------------------
      key = 'PTYPE5  '
      val = q//'DATE    '//q//'          '
      com = ' '
      WRITE(card(30),1) key, val, com
c---------------------------------------------------------------
      key = 'PSCAL5  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(31),1) key, val, com
c---------------------------------------------------------------
      key = 'PZERO5  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(32),1) key, val, com
c---------------------------------------------------------------
      key = 'PTYPE6  '
      val = q//'DATE    '//q
      com = ' /'
      WRITE(card(33),1) key, val, com
c---------------------------------------------------------------
      key = 'PSCAL6  '
      WRITE(val,'(e20.8)') scale
      com = ' /'
      WRITE(card(34),1) key, val, com
c---------------------------------------------------------------
      key = 'PZERO6  '
      WRITE(val,'(e20.8)') zero
      com = ' /'
      WRITE(card(35),1) key, val, com
c---------------------------------------------------------------
      key = 'END     '
      WRITE(card(36),4) key
c---------------------------------------------------------------
c     DO j = 34, 36
c         WRITE(card(j),3)
c     END DO
c
      irec = irec + 1
      WRITE(lu,rec=irec) card
c
      RETURN
      END
c***********************************************************************
      SUBROUTINE uvf_ext(lu, freq,    irec)
c***********************************************************************
c
c  Write the UV FITS extension header on logical unit lu.
c
c  INPUT:  lu     = logical unit opened to UV FITS file
c          freq   = reference frequency in Hz
c          irec   = record number for previous record
c                   ** MODIFIED ON OUTPUT **
c
c  OUTPUT: irec   = record number for the last record in header
c
      IMPLICIT NONE
      INTEGER lu, irec, j, i
      REAL*8 freq
      CHARACTER card(36)*80, key*8, val*20, com*50, q*1
c
c  Define variables used in extension and extension data:
c
      CHARACTER*8 antnm(6)
      REAL*8 pos(3,6)
      INTEGER*4 stano(6),mntsta
      LOGICAL mustswap
c
c  Initialize values for the extension file:
c
      antnm(1) = 'HHHH    '
      antnm(2) = 'IIII    '
      antnm(3) = 'PPPP    '
      antnm(4) = 'UUUU    '
      antnm(5) = 'VVVV    '
      antnm(6) = 'FFFF    '
      stano(1) = 1
      stano(2) = 2
      stano(3) = 3
      stano(4) = 4
      stano(5) = 5
      stano(6) = 6
      mntsta = 0
      DO  j = 1, 6
          DO  i = 1, 3
              pos(i,j) =  0.0
          ENDDO
      ENDDO
c
c  Swap bytes if necessary:
c
      IF (mustswap()) THEN
          CALL swapr8(pos, 18)
          CALL swapi4(stano, 6)
          CALL swapi4(mntsta, 1)
      ENDIF
c
c  The character q is a quote ('):
c
      q = char(39)
c
c  Card formats: 1 = normal, 3 = blank
c
1     FORMAT(a8,'= ',a20,a50)
3     FORMAT(80x)
4     FORMAT(a8,72x)
c---------------------------------------------------------------
      key = 'XTENSION'
      val = q//'A3DTABLE'//q//'          '
      com = ' /EXTENSION TYPE'
      WRITE(card(1),1) key, val, com
c---------------------------------------------------------------
      key = 'BITPIX  '
      val = '                   8'
      com = ' /BINARY DATA'
      WRITE(card(2),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS   '
      val = '                   2'
      com = ' /TABLE IS A MATRIX'
      WRITE(card(3),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS1  '
      val = '                  40'
      com = ' /WIDTH OF TABLE IN BYTES'
      WRITE(card(4),1) key, val, com
c---------------------------------------------------------------
      key = 'NAXIS2  '
      val = '                   6'
      com = ' /NUMBER OF ENTRIES IN TABLE'
      WRITE(card(5),1) key, val, com
c---------------------------------------------------------------
      key = 'PCOUNT  '
      val = '                   0'
      com = ' /THERE ARE ZERO PARAMETERS NEEDED'
      WRITE(card(6),1) key, val, com
c---------------------------------------------------------------
      key = 'GCOUNT  '
      val = '                   1'
      com = ' /THERE IS ONE GROUP NEEDED'
      WRITE(card(7),1) key, val, com
c---------------------------------------------------------------
      key = 'TFIELDS '
      val = '                   5'
      com = ' /NUMBER OF FIELDS IN EACH ROW'
      WRITE(card(8),1) key, val, com
c---------------------------------------------------------------
      key = 'EXTNAME '
      val = q//'AIPS AN '//q
      com = ' /AIPS TABLE FILE'
      WRITE(card(9),1) key, val, com
c---------------------------------------------------------------
      key = 'EXTVER  '
      val = '                   1'
      com = ' /VERSION NUMBER OF TABLE'
      WRITE(card(10),1) key, val, com
c---------------------------------------------------------------
      key = 'NUMORB  '
      val = '                   0'
      com = '                   '
      WRITE(card(11),1) key, val, com
c---------------------------------------------------------------
      key = 'TFORM1  '
      val = q//'8A      '//q
      com = ' /FORTRAN FORMAT OF FIELD 1'
      WRITE(card(12),1) key, val, com
c---------------------------------------------------------------
      key = 'TTYPE1  '
      val = q//'ANNAME  '//q
      com = ' /TYPE (HEADING) OF FIELD 1'
      WRITE(card(13),1) key, val, com
c---------------------------------------------------------------
      key = 'TUNIT1  '
      val = q//'        '//q
      com = ' /PHYSICAL UNITS OF FIELD 1'
      WRITE(card(14),1) key, val, com
c---------------------------------------------------------------
      key = 'TFORM2  '
      val = q//'3D      '//q
      com = ' /FORTRAN FORMAT OF FIELD 2'
      WRITE(card(15),1) key, val, com
c---------------------------------------------------------------
      key = 'TTYPE2  '
      val = q//'STABXYZ '//q
      com = ' /TYPE (HEADING) OF FIELD 2'
      WRITE(card(16),1) key, val, com
c---------------------------------------------------------------
      key = 'TUNIT2  '
      val = q//'METERS  '//q
      com = ' /PHYSICAL UNITS OF FIELD 2'
      WRITE(card(17),1) key, val, com
c---------------------------------------------------------------
      key = 'TFORM3  '
      val = q//'1J      '//q
      com = ' /FORTRAN FORMAT OF FIELD 3'
      WRITE(card(18),1) key, val, com
c---------------------------------------------------------------
      key = 'TTYPE3  '
      val = q//'NOSTA   '//q
      com = ' /TYPE (HEADING) OF FIELD 3'
      WRITE(card(19),1) key, val, com
c---------------------------------------------------------------
      key = 'TUNIT3  '
      val = q//'        '//q
      com = ' /PHYSICAL UNITS OF FIELD 3'
      WRITE(card(20),1) key, val, com
c---------------------------------------------------------------
      key = 'TFORM4  '
      val = q//'1J      '//q
      com = ' /FORTRAN FORMAT OF FIELD 4'
      WRITE(card(21),1) key, val, com
c---------------------------------------------------------------
      key = 'TTYPE4  '
      val = q//'MNTSTA  '//q
      com = ' /TYPE (HEADING) OF FIELD 4'
      WRITE(card(22),1) key, val, com
c---------------------------------------------------------------
      key = 'TUNIT4  '
      val = q//'        '//q
      com = ' /PHYSICAL UNITS OF FIELD 4'
      WRITE(card(23),1) key, val, com
c---------------------------------------------------------------
      key = 'TFORM5  '
      val = q//'0D      '//q
      com = ' /FORTRAN FORMAT OF FIELD 5'
      WRITE(card(24),1) key, val, com
c---------------------------------------------------------------
      key = 'TTYPE5  '
      val = q//'ORBPARM '//q
      com = ' /TYPE (HEADING) OF FIELD 5'
      WRITE(card(25),1) key, val, com
c---------------------------------------------------------------
      key = 'TUNIT5  '
      val = q//'        '//q
      com = ' /PHYSICAL UNITS OF FIELD 5'
      WRITE(card(26),1) key, val, com
c---------------------------------------------------------------
      key = 'FREQ    '
      WRITE(val,'(e20.8)') freq
      com = ' /'
      WRITE(card(27),1) key, val, com
c---------------------------------------------------------------
      key = 'ARRNAM  '
      val = q//'HIPPARC '//q
      com = ' /'
      WRITE(card(28),1) key, val, com
c---------------------------------------------------------------
      key = 'END     '
      WRITE(card(29),4) key
c---------------------------------------------------------------
      DO j = 30, 36
        WRITE(card(j),3)
      END DO
c
      irec = irec + 1
      WRITE(lu,rec=irec) card
c
      irec = irec + 1
      WRITE(lu,rec=irec)
     +    (antnm(i),pos(1,i),pos(2,i),pos(3,i),stano(i),mntsta,i=1,6),
     +    (' ',i=1,2640)
c
      RETURN
      END
c***********************************************************************
      CHARACTER*9 FUNCTION today()
c***********************************************************************
c
c  User supplied function returning the current date
c  in character form, e.g.: '22-Mar-97'
c
      IMPLICIT NONE
      today = '22-Mar-97'
      RETURN
      END
c***********************************************************************
      SUBROUTINE swapi4(val4, n)
c***********************************************************************
c
c  Swap bytes, in groups of four, in INTEGER*4 array val4(1:n)
c
      IMPLICIT NONE
      INTEGER n, i, j
      INTEGER*4 val4(n), tmp, pmt
      CHARACTER str*4, rts*4
      EQUIVALENCE (tmp, str), (pmt, rts)
      DO  i = 1, n
          tmp = val4(i)
          DO  j = 1, 4
              rts(j:j) = str(5-j:5-j)
          ENDDO
          val4(i) = pmt
      ENDDO
      RETURN
      END
c***********************************************************************
      SUBROUTINE swapr4(val4, n)
c***********************************************************************
c
c  Swap bytes, in groups of four, in REAL*4 array val4(1:n)
c
      IMPLICIT NONE
      INTEGER n, i, j
      REAL*4 val4(n), tmp, pmt
      CHARACTER str*4, rts*4
      EQUIVALENCE (tmp, str), (pmt, rts)
      DO  i = 1, n
          tmp = val4(i)
          DO  j = 1, 4
              rts(j:j) = str(5-j:5-j)
          ENDDO
          val4(i) = pmt
      ENDDO
      RETURN
      END
c***********************************************************************
      SUBROUTINE swapr8(val8, n)
c***********************************************************************
c
c  Swap bytes, in groups of eight, in REAL*8 array val8(1:n)
c
      IMPLICIT NONE
      INTEGER n, i, j
      REAL*8 val8(n), tmp, pmt
      CHARACTER str*8, rts*8
      EQUIVALENCE (tmp, str), (pmt, rts)
      DO  i = 1, n
          tmp = val8(i)
          DO  j = 1, 8
              rts(j:j) = str(9-j:9-j)
          ENDDO
          val8(i) = pmt
      ENDDO
      RETURN
      END
c***********************************************************************
      LOGICAL FUNCTION mustswap()
c***********************************************************************
c
c  Determines whether byte swapping is required for storing binary
c  data according to IEEE standard
c
      IMPLICIT NONE
      CHARACTER str*4
      INTEGER*4 ival
      EQUIVALENCE (str, ival)
      ival = 1
      mustswap = (str .EQ. char(1)//char(0)//char(0)//char(0))
      RETURN
      END