PROGRAM FRINGE C C Perform a fairly realistic simulation of the performance of the C BOA and AI multi-spectral fringe detectors. C Makes 1). a plot of transform phase versus fringe phase C 2). a plot of transform amplitude versus fringe phase. C IMPLICIT UNDEFINED (A-Z) INTEGER *4 NW, ND PARAMETER ( ND = 16) PARAMETER ( NW = 48) REAL *8 PI, RAD PARAMETER ( PI = 3.1415926535D0 ) PARAMETER ( RAD = PI/180. ) INTEGER *4 I, J, K, NWAVE, NPHOT, IX, IY, COLOR(10), IERR, JJ REAL *4 X, Y, XMIN, XMAX, YMIN, YMAX, ZMAX, XHI, XLO, YHI, YLO REAL *8 WAVE, DELAY, SUM, ARG REAL *8 RTRANS(NW), ITRANS(NW) REAL *4 TRANSF(NW), WORK(NW), TR(6) REAL *4 WAVES(NW), XPLOT(NW), YPLOT(NW) REAL *4 PHASE(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, NSAMPLE REAL *8 WAVE0, WAVE1, WAVEA, DELAY0, DELAY1, DELTAD, AIR, VIS, $ DELTADEL CHARACTER *64 PLTFILE, RESFILE COMMON / CONFIG / NWAVES, NDELAYS, NPHOT, WAVE0, WAVE1, WAVEA, $ DELAY0, DELAY1, DELTAD, AIR, VIS, $ NSAMPLE, DELTADEL COMMON / FILES / PLTFILE, RESFILE C C Contour intervals in fraction of maximum C 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 OPEN (UNIT=7,FILE=RESFILE,STATUS='UNKNOWN',IOSTAT=IERR) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' CANNOT OPE OUTPUT FILE ' END IF DELTAD = DELTAD - DELTADEL CALL PGBEGIN (0, PLTFILE, 1, 1 ) C C DO JJ = 1, NSAMPLE DO I = 1, NW RTRANS(I) = 0. ITRANS(I) = 0. TRANSF(I) = 0. END DO DELTAD = DELTAD + DELTADEL 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. C 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 C C Accumulate the fourier transform. C All we need is one row C ZMAX = 0. SUM = 0. DO K = 1, NWAVES J = 1 ARG = 2.*PI*( FLOAT((K-NWAVES/2)*IX)/FLOAT(NWAVES) $ + FLOAT(-IY)/FLOAT(NDELAYS) ) RTRANS(K) = RTRANS(K) + COS( ARG ) ITRANS(K) = ITRANS(K) + SIN( ARG ) TRANSF(K) = RTRANS(K)*RTRANS(K) $ + ITRANS(K)*ITRANS(K) SUM = SUM + TRANSF(K) ZMAX = MAX ( ZMAX, TRANSF(K) ) END DO END DO C C Make plots of the power spectrum. C XMIN = -NWAVES/2 XMAX = NWAVES/2 YMIN = 0. YMAX = ZMAX CALL PGSCI (5) CALL PGSCH ( 1.5 ) 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 ') 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, NWAVES YPLOT(I) = TRANSF(I) XPLOT(I) = I - NWAVES/2 END DO CALL PGLINE ( NWAVES, XPLOT, YPLOT ) C C Make a plot of the phase of the transform. C XMIN = -NWAVES/2 XMAX = NWAVES/2 YMIN = -200. YMAX = 200. CALL PGSCI (5) CALL PGSCH ( 1.5 ) 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,'TRANSFORM PHASE') 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, NWAVES YPLOT(I) = ATAN2( ITRANS(I), RTRANS(I) )/RAD END DO CALL PGLINE ( NWAVES, XPLOT, YPLOT ) WRITE(7,1400) DELTAD, NPHOT, (TRANSF(NWAVES/2+K),K=1,5), $ ( YPLOT(NWAVES/2+K),K=1,5) IF ( PLTFILE .NE. '/EGA' ) THEN WRITE(6,1400) DELTAD, NPHOT, (TRANSF(NWAVES/2+K),K=1,5), $ ( YPLOT(NWAVES/2+K),K=1,5) END IF END DO 1400 FORMAT ( F7.3, I6, 2(4X,5F12.2) ) C C C SKIP THE REST OF THE PLOTS C CALL PGEND close (7) GO TO 100 900 CONTINUE STOP 1100 FORMAT ( 16I5 ) 1200 FORMAT ( 2F7.1, F7.3, F7.1, F10.6 ) END SUBROUTINE POISSON ( N, MEAN ) C C On input, MEAN is the mean of a poisson distribution C On output, N is a poisson distributed random variable. C IMPLICIT UNDEFINED (A-Z) REAL *4 MEAN, RAN0, CUTOFF, A, B INTEGER *4 N, IA, ISEED EXTERNAL RAN0 DATA ISEED / 171717 / CUTOFF = 5.0*MEAN 100 CONTINUE A = RAN0(ISEED) B = RAN0(ISEED) IA = CUTOFF*A P = MEAN**(IA)*EXP(-MEAN)/FACT(IA) IF ( REAL *4 FUNCTION FACT(K) IMPLICIT UNDEFINED ( A-Z) SAVE INTEGER *4 K, I REAL*4 F(20) LOGICAL FIRST DATA FIRST / .TRUE. / IF ( FIRST ) THEN F(1) = 1. DO I = 2, 20 F(I) = FLOAT(I)*F(I-1) END DO END IF IF ( K .LE. 0 ) THEN FACT = 1. ELSE IF ( K .GT. 20 ) THEN WRITE(6,*) ' FACTORIAL IS TOO BIG = ', K FACT = 1. ELSE FACT = F(K) END IF RETURN 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, NSAMPLE REAL *8 WAVE0, WAVE1, WAVEA, DELAY0, DELAY1, DELTAD, AIR, VIS, $ DELTADEL CHARACTER *64 PLTFILE, RESFILE COMMON / CONFIG / NWAVES, NDELAYS, NPHOT, WAVE0, WAVE1, WAVEA, $ DELAY0, DELAY1, DELTAD, AIR, VIS, $ NSAMPLE, DELTADEL 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 = 29 ) REAL*8 INPARS(NPARS), INVALS(NPARS), ENDMRK CHARACTER *1 BACKSLSH CHARACTER *8 LABELS(NPARS), CVALS(1:3) INTEGER*4 NWAVES, NDELAYS, NPHOT, NSAMPLE REAL *8 WAVE0, WAVE1, WAVEA, DELAY0, DELAY1, DELAY, AIR, VIS, $ DELTADEL CHARACTER *64 PLTFILE, RESFILE COMMON / CONFIG / NWAVES, NDELAYS, NPHOT, WAVE0, WAVE1, WAVEA, $ DELAY0, DELAY1, DELAY, AIR, VIS, $ NSAMPLE, DELTADEL COMMON / FILES / PLTFILE, RESFILE DATA BACKSLSH / '\\' / DATA LABELS / 'PLTFILE ', 7*' ', 'RESFILE ', 7*' ', $ 'WAVE0 ', 'WAVE1 ', 'WAVEA ', 'DELAY0 ', 'DELAY1 ', $ 'NWAVES ', 'NDELAYS ', 'DELAY ', 'AIR ', 'NPHOT ', $ 'VIS ', 'NSAMPLE ', 'DELTADEL' / 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, 10.00D0, 0.2D0 / 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) NSAMPLE = INVALS(28) DELTADEL = INVALS(29) C---------------------------------------------------------------------- RETURN 1100 FORMAT ( A8 ) END