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 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 (microns)
      DELAY1 =  1.0D0      ! upper limit to delay range (microns)
      VIS    =  0.8D0      ! visibility (independent of wavelength)
C-----------------------------------------------------------------------
      DO I = 1, NWAVES
          DO J = 1, NDELAYS
              IMAGE (I,J) = 0.
              RTRANS(I,J) = 0.
              ITRANS(I,J) = 0.
              TRANSF(I,J) = 0.
          END DO
      END DO
      CALL GETPARM
      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 (WAVELENGTHS)')
      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
          CALL PHOTON  ( WAVE, DELAY )
          X = WAVE
          Y = DELAY
          CALL PGPOINT ( 1, X, Y, 4 )
C
C     Put each photon in its correct bin
C
          IX = 1. + FLOAT(NWAVES)*(WAVE-WAVE0)/(WAVE1-WAVE0)
          IY = 1. + FLOAT(NDELAYS)*(DELAY-DELAY0)/(DELAY1-DELAY0)
          IMAGE(IX,IY) = IMAGE(IX,IY) + 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
          ALEV(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,
     $                 ALEV(I), 1, TR )
      END DO
C
C
C     Plot one slice of the fourier transform at higher resolution.
C
      ZMAX = 0.
      J = -2
      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(-2*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
  900 CONTINUE
      STOP
 1100 FORMAT ( 16I5 )
 1200 FORMAT ( 2F7.1, F7.3, F7.1, F10.6 )
      END