PROGRAM FRINGE
C
C   Perform a fairly realistic simulation of the performance of the 
C   BOA and AI multi-spectral fringe detectors.
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 NW, ND
      PARAMETER ( ND = 16)
      PARAMETER ( NW = 48)
      REAL *8 PI
      PARAMETER ( PI = 3.1415926535D0 )
      INTEGER *4 I, J, K, NWAVE, NPHOT, IX, IY, COLOR(10), IERR
      REAL *4 X, Y, XMIN, XMAX, YMIN, YMAX, ZMAX, XHI, XLO, YHI, YLO
      REAL *8 WAVE, DELAY, SUM, ARG
      REAL *8 RTRANS(NW,ND), ITRANS(NW,ND)
      REAL *4 IMAGE(NW,ND),  TRANSF(NW,ND), WORK(NW), ALEV(10), TR(6)
      REAL *4 BLEV(10), WAVES(NW)
      REAL *4 SLICE(NW), RSLICE(NW), ISLICE(NW), PHASE(NW), BOZO(NW)
      CHARACTER *50 TITLE, TITLE2
C
C   This common block contains the variables needed to generate the
C   desired photon distribution.  It is needed by PHOTOM
C
      INTEGER*4 NWAVES, NDELAYS, NPHOT
      REAL *8 WAVE0, WAVE1, WAVEA, DELAY0, DELAY1, DELTAD, AIR, VIS
      CHARACTER *64 PLTFILE, RESFILE
      COMMON / CONFIG / NWAVES, NDELAYS, NPHOT, WAVE0, WAVE1, WAVEA,
     $                  DELAY0, DELAY1,  DELTAD, AIR,   VIS
      COMMON / FILES  / PLTFILE, RESFILE
C
C  Contour intervals in fraction of maximum
C
      DATA ALEV / .002, .004, .008, .015, .030,
     $            .060, .125, .250, .500, .950  /
C
C    Color sequence
C
      DATA COLOR / 2, 2, 2, 14, 7, 3, 4, 11, 7, 1 /
C-----------------------------------------------------------------------
      WAVE0  =  0.4D0      ! minimum wavelength in array
      WAVE1  =  0.8D0      ! maximum wavelength in array
      DELAY0 = -1.0D0      ! lower limit to delay range (wavelengths)
      DELAY1 =  1.0D0      ! upper limit to delay range (wavelengths)
      VIS    =  0.8D0      ! visibility (independent of wavelength)
C-----------------------------------------------------------------------
  100 continue
      CALL GETPARM
C
C   Determine the central wavelengths of each filter.
C
      DO I = 1, NW
          WAVES(I) = WAVE0 + (WAVE1-WAVE0)*(FLOAT(I)-.5)/FLOAT(NWAVES)
      END DO
      DO I = 1, NW
          DO J = 1, ND
              IMAGE (I,J) = 0.
              RTRANS(I,J) = 0.
              ITRANS(I,J) = 0.
              TRANSF(I,J) = 0.
          END DO
      END DO
      close (7)
      OPEN (UNIT=7,FILE='RESPLOT',STATUS='UNKNOWN',IOSTAT=IERR)
      IF ( IERR .NE. 0 ) THEN
          WRITE(6,*) ' CANNOT OPE OUTPUT FILE '
      END IF
      CALL PGBEGIN (0, PLTFILE, 1, 1 )
      WRITE(TITLE, 1300) DELTAD, AIR
      WRITE(TITLE2,1305) NPHOT
 1300 FORMAT( ' TRACKING ERROR = ', F5.2, ' AIR MISMATCH = ', F5.2 )
 1305 FORMAT( ' NUMBER OF PHOTONS = ', I5 )
C
C  X axis is wavelength.
C  Y axis is time samples.
      CALL PGSCI   (5)
      CALL PGSCH   ( 1.5 )
      XMIN = WAVE0
      XMAX = WAVE1
      YMIN = DELAY0
      YMAX = DELAY1
C      CALL PGRNGE  ( XMIN, XMAX, XLO, XHI )
C      CALL PGRNGE  ( YMIN, YMAX, YLO, YHI )
C      CALL PGENV   ( XLO, XHI, YLO, YHI, 0, 1 )
      CALL PGENV   ( XMIN, XMAX, YMIN, YMAX, 0, 1 )
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, 'SPECTRAL FRINGE DETECTOR' )
      CALL PGSCH   ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, 'WAVELENGTH MICRONS' )
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, 'DELAY (MICRONS)')
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'IMAGE PLANE SIMULATION' )
      CALL PGMTEXT ( 'R', 2.5, 0.5, 0.5, TITLE  )
      CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, TITLE2 )
      CALL PGSCI   ( 2 )
      DO I = 1, NPHOT
C
C     Put each photon in its correct bin.  First the wavelength is
C     quantized.  Then the central wavelength of that channel is 
C     used to determine the delay binning.
C
  134     CONTINUE
          CALL PHOTON  ( WAVE, DELAY )
          X = WAVE
          Y = DELAY
          IX = 1. + FLOAT(NWAVES)*(WAVE-WAVE0)/(WAVE1-WAVE0)
          IY = 1. + FLOAT(NDELAYS)*(DELAY/WAVES(IX)-DELAY0)
     $                            /(DELAY1-DELAY0)
          IF ( (IY.LE.0) .OR. (IY.GT.NDELAYS) ) GO TO 134
          IMAGE(IX,IY) = IMAGE(IX,IY) + 1.
          CALL PGPOINT ( 1, X, Y, 1 )
C
C     Accumulate the fourier transform
C
          ZMAX = 0.
          SUM  = 0.
          DO K = 1, NWAVES
              DO J = 1, NDELAYS
                  ARG = 2.*PI*( FLOAT((K-NWAVES/2)*IX)/FLOAT(NWAVES)
     $                        + FLOAT((J-NDELAYS/2)*IY)/FLOAT(NDELAYS) )
                  RTRANS(K,J) = RTRANS(K,J) + COS( ARG )
                  ITRANS(K,J) = ITRANS(K,J) + SIN( ARG )
                  TRANSF(K,J) = RTRANS(K,J)*RTRANS(K,J)
     $                        + ITRANS(K,J)*ITRANS(K,J)
                  SUM = SUM + TRANSF(K,J)
                  ZMAX = MAX ( ZMAX, TRANSF(K,J) )
              END DO
          END DO
      END DO
C
C  Set up the Transformation matrix for the contour plot
C
      TR(1) = -NWAVES/2.
      TR(2) =  1.0
      TR(3) =  0.0
      TR(4) = -NDELAYS/2.
      TR(5) =  0.0
      TR(6) =  1.0
C
C  Make a contour plot of the power spectrum.
C
      XMIN = -NWAVES/2
      XMAX =  NWAVES/2
      YMIN = -NDELAYS/2
      YMAX =  NDELAYS/2
      CALL PGSCI   (5)
      CALL PGSCH   ( 1.5 )
C      CALL PGRNGE  ( XMIN, XMAX, XLO, XHI )
C      CALL PGRNGE  ( YMIN, YMAX, YLO, YHI )
C      CALL PGENV   ( XLO, XHI, YLO, YHI, 0, 2 )
      CALL PGENV   ( XMIN, XMAX, YMIN, YMAX, 0, 2 )
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, 'SPECTRAL FRINGE DETECTOR' )
      CALL PGSCH   ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, 'WAVELENGTH AXIS (PIXELS)' )
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, 'DELAY AXIS (PIXELS)'      )
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'POWER SPECTRUM SIMULATION')
      CALL PGMTEXT ( 'R', 2.5, 0.5, 0.5, TITLE  )
      CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, TITLE2 )
      CALL PGSCI   ( 2 )
      DO I = 1, 10
          BLEV(I) = ZMAX * ALEV(I)
      END DO
      CALL PGSCI( 2 )
      DO I = 1, 10
          CALL PGSCI(COLOR(I))
          CALL PGCONT( TRANSF, NW, ND, 1, NWAVES, 1, NDELAYS,
     $                 BLEV(I), 1, TR )
      END DO
C
C
C    SKIP THE REST OF THE PLOTS
C
      call pgend
      GO TO 100
C
C
C     Plot one slice of the fourier transform at higher resolution.
C
      ZMAX = 0.
      J = -1
      DO K = 1, NWAVES
          RSLICE(K) = 0.
          ISLICE(K) = 0.
      END DO
      DO K = 1, NWAVES
          DO IX = 1, NWAVES
              DO IY = 1, NDELAYS
                  ARG = 2.*PI*( FLOAT((K-NWAVES/2)*IX)/FLOAT(NWAVES)
     $                        + FLOAT(J*IY)/FLOAT(NDELAYS) )
                  RSLICE(K) = RSLICE(K) + IMAGE(IX,IY)*COS( ARG )
                  ISLICE(K) = ISLICE(K) + IMAGE(IX,IY)*SIN( ARG )
              END DO
          END DO
      END DO
      YMAX = 0.
      YMIN = 0.
      DO K = 1, NWAVES
          SLICE(K) = RSLICE(K)*RSLICE(K) + ISLICE(K)*ISLICE(K)
          ZMAX = MAX ( ZMAX, SLICE(K) )
          PHASE(K) = 180.*ATAN2( RSLICE(K), ISLICE(K) )/PI 
          YMAX = MAX( YMAX, PHASE(K) )
          YMIN = MIN( YMIN, PHASE(K) )
          BOZO(K) = K-NWAVES/2
      END DO
      CALL PGSCI   (5)
      CALL PGSCH   ( 1.5 )
      XMIN = BOZO(1)
      XMAX = BOZO(NWAVES)
      CALL PGENV   ( XMIN, XMAX, 0., ZMAX, 0, 1 )
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, 'SPECTRAL FRINGE DETECTOR' )
      CALL PGSCH   ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, 'WAVELENGTH MICRONS'   )
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, 'POWER' )
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'POWER SPECTRUM SLICE' )
      CALL PGMTEXT ( 'R', 2.5, 0.5, 0.5, TITLE  )
      CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, TITLE2 )
      CALL PGSCI   ( 2 )
      CALL PGSCI   (5)
      CALL PGSCH   ( 1.5 )
      CALL PGLINE ( NWAVES, BOZO, SLICE )
C
      CALL PGENV   ( XMIN, XMAX, -180., 180., 0, 1 )
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, 'SPECTRAL FRINGE DETECTOR' )
      CALL PGSCH   ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, 'WAVELENGTH MICRONS'   )
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, 'PHASE (RADIANS)' )
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'POWER SPECTRUM SLICE' )
      CALL PGMTEXT ( 'R', 2.5, 0.5, 0.5, TITLE  )
      CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, TITLE2 )
      CALL PGSCI   ( 2 )
      CALL PGSCI   (5)
      CALL PGSCH   ( 1.5 )
      CALL PGLINE ( NWAVES, BOZO, PHASE )
C
      CALL PGEND
C
C    Dump the image and power spectrum to disk
C
      WRITE(7,*) '  THE IMAGE '
      DO I = 1, NWAVES
          WRITE(7,1100) ( INT(IMAGE(I,J)),J=1,NDELAYS)
      END DO
      SUM = 999.*SUM/ZMAX
      WRITE(7,*) ' THE POWER SPECTRUM '
      WRITE(7,*) ' MAX ELEMENT IN ARRAY = ', ZMAX
      WRITE(7,*) ' SUM OF ARRAY = ', SUM
      DO I = 1, NWAVES
          WRITE(7,1100) ( INT(999.*TRANSF(I,J)/ZMAX),J=1,NDELAYS )
      END DO
      WRITE(7,*) ' MAX ELEMENT IN ARRAY = ', ZMAX
      WRITE(7,*) ' SUM OF ARRAY = ', SUM
      go to 100
  900 CONTINUE
      STOP
 1100 FORMAT ( 16I5 )
 1200 FORMAT ( 2F7.1, F7.3, F7.1, F10.6 )
      END
      SUBROUTINE PHOTON ( WAVE, DELAY )
C
C  Each time PHOTON is called it returns the wavelength and 
C  delay of a detected photon, both in microns.
C
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      REAL *8 PI
      PARAMETER ( PI = 3.1415926653589793D0 )
      INTEGER *4 ISEED
      REAL *8  P, T, E, INDEX, WAVE, DELAY, FLUX, PHASE, TEMP
      REAL *8  A, B, C, F
      REAL *4 RAN0
      EXTERNAL RAN0, FLUX
C
      INTEGER*4 NWAVES, NDELAYS, NPHOT
      REAL *8 WAVE0, WAVE1, WAVEA, DELAY0, DELAY1, DELTAD, AIR, VIS
      CHARACTER *64 PLTFILE, RESFILE
      COMMON / CONFIG / NWAVES, NDELAYS, NPHOT, WAVE0, WAVE1, WAVEA,
     $                  DELAY0, DELAY1,  DELTAD, AIR,   VIS
      COMMON / FILES  / PLTFILE, RESFILE
C
C  Atmospheric model for index of refraction calculation.
C  P = pressure, T = temperature, E= partial pressure of water vapor
C
      DATA  P, T, E /  1013., 285.,  0. /
C
      DATA ISEED / 12534 /
C
  100 CONTINUE
      A = RAN0(ISEED)
      B = RAN0(ISEED)
      C = RAN0(ISEED)
      WAVE  =  WAVE0 + (WAVE1 -WAVE0 )*A
      DELAY = DELAY0 + (DELAY1-DELAY0)*B
      CALL AINDEX ( P, T, E, WAVEA, TEMP  )
      CALL AINDEX ( P, T, E, WAVE,  INDEX )
      INDEX = INDEX - TEMP
C
C  The phase has contributions from the stroke (DELAY), an error in
C  fringe tracking (DELTAD), and an error in the air path length (AIR).
C  Note that INDEX is the difference in index of refraction between
C  air and a vaccuum.
C
      PHASE = 2.*PI*( DELAY + DELTAD + 1.D6*AIR*INDEX ) / WAVE
      F = FLUX(WAVE)*( 1 + VIS*COS(PHASE) ) / (1.+VIS)
      IF ( C .GT. F ) THEN
          GO TO 100
      END IF
      RETURN
      END
      REAL *8 FUNCTION FLUX(WAVE)
C
C   Flux distribution of star.
C   WAVE is the wavelength is microns
C   FLUX must not exceed unity.
C
      IMPLICIT UNDEFINED (A-Z)
      REAL *8 WAVE
      IF ( WAVE .LE. 0.D0 ) THEN
          WRITE(6,*) ' WAVELENGTH IS TOO SMALL '
      END IF
      FLUX = 1.0D0
      RETURN
      END
      SUBROUTINE AINDEX ( P, T, E, WAVE, INDEX )
C
C   Calculates the index of refraction of air for pressure P (mb)
C   temperature T (kelvins) paritial pressure of water E (mb) and
C   wavelength WAVE (microns).
C
      IMPLICIT UNDEFINED (A-Z)
      REAL *8 P, T, E, INDEX, A, B, K, WAVE
C
      K = 1./(WAVE*WAVE)
      A = 1.E-8*( 2372.434 + 684255.24/(130.-K) + 4549.40/(38.9-K) )
      B = A - 1.E-8*(6487.31 +K*( 58.058 +K*(-0.71150 +K*0.08851 ) ) )
      INDEX = ( A*P - B*E ) / T
      RETURN
      END
      REAL *4 FUNCTION RAN0(IDUM)
C
C  Improver for the canned radnom number generator.
C  see Numerical Recipes Press, Flannery, Teukolsky, Vetterling
C  for details.
C
      IMPLICIT UNDEFINED (A-Z)
      REAL *4 DUM, V(97), Y
      REAL *4 RAN
      INTEGER *4 IFF, J, IDUM, ISEED
      EXTERNAL RAN
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
          IFF = 1
          ISEED = ABS(IDUM)
          IDUM = 1
          DO J = 1,97
              DUM = RAN(ISEED)
          END DO
          DO J=1,97
              V(J) = RAN(ISEED)
          END DO
          Y = RAN(ISEED)
      END IF
      J = 1 + INT(97.*Y)
      IF (J.GT.97 .OR. J.LT.1 ) THEN
          WRITE(6,*) ' RAN SCREWED UP '
          STOP
      END IF
      Y = V(J)
      RAN0 = Y
      IF ( (RAN0 .LT. 0.) .OR. (RAN0 .GT. 1.0) ) THEN
          WRITE(6,*) ' RAN0 IS OUT OF RANGE '
          STOP
      END IF
      V(J) = RAN(ISEED)
      RETURN
      END
      SUBROUTINE GETPARM
C
C            by David Mozurkewich         Sept 20, 1989
C
C      This subroutine gets the parameters input from the terminal
C      for the spectral fringe detection simulator.
C      
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      INTEGER *4 MODE, NPARS, I
      PARAMETER (MODE   = 2  )
      PARAMETER ( NPARS = 27 )
      REAL*8  INPARS(NPARS), INVALS(NPARS), ENDMRK
      CHARACTER *1 BACKSLSH
      CHARACTER *8 LABELS(NPARS), CVALS(1:3)
      INTEGER*4 NWAVES, NDELAYS, NPHOT
      REAL *8 WAVE0, WAVE1, WAVEA, DELAY0, DELAY1, DELAY, AIR, VIS
      CHARACTER *64 PLTFILE, RESFILE
      COMMON / CONFIG / NWAVES, NDELAYS, NPHOT, WAVE0, WAVE1, WAVEA,
     $                  DELAY0, DELAY1,  DELAY, AIR,   VIS
      COMMON / FILES  / PLTFILE, RESFILE
      DATA BACKSLSH / '\\' /
      DATA LABELS   /   'PLTFILE ', 7*' ',  'RESFILE ', 7*' ',
     $     'WAVE0   ', 'WAVE1   ', 'WAVEA   ', 'DELAY0  ', 'DELAY1  ',
     $     'NWAVES  ', 'NDELAYS ', 'DELAY   ', 'AIR     ', 'NPHOT   ',
     $     'VIS     '  /
      DATA CVALS  /   '/', '/EGA', ' ' /
      DATA INVALS  / 8*0.D0,  8*0.D0,
     $     0.45D0,     0.80D0,      0.8D0,     -1.0D0,       1.0D0,
     $    16.00D0,     8.00D0,      0.0D0,      0.0D0,     500.0D0,
     $     0.80D0     /
C-----------------------------------------------------------------------
C     Set up parameter arrays for KEYIN
      DO I = 1, NPARS
          READ ( LABELS(I), 1100 ) INPARS(I)
      END DO
      READ (CVALS(1), 1100 ) ENDMRK
      READ (CVALS(2), 1100 ) INVALS(1)
      DO I = 2, 16
          READ (CVALS(3), 1100 ) INVALS(I)
      END DO
      READ (LABELS(9), 1100 ) INVALS(9)
C----------------------------------------------------------------------
C               Read parameters from input file using KEYIN
      WRITE (6, *) ' input configuration.   / when finished '
      CALL KEYIN (INPARS, INVALS, NPARS, ENDMRK, MODE, 5, 6 )
C----------------------------------------------------------------------
      WRITE ( PLTFILE, '(8A8)' ) (INVALS(I), I=1, 8)
      WRITE ( RESFILE, '(8A8)' ) (INVALS(I), I=9,16)
      WAVE0  = INVALS(17)
      WAVE1  = INVALS(18)
      WAVEA  = INVALS(19)
      DELAY0 = INVALS(20)
      DELAY1 = INVALS(21)
      NWAVES = INVALS(22)
      NDELAYS= INVALS(23)
      DELAY  = INVALS(24)
      AIR    = INVALS(25)
      NPHOT  = INVALS(26)
      VIS    = INVALS(27)
C----------------------------------------------------------------------
      RETURN
 1100 FORMAT ( A8 )
      END