      PROGRAM AIR
C
C   Make a plot of air dispersion versus wavelength.
C
      IMPLICIT UNDEFINED (A-Z)
      REAL *8 PI, RAD
      PARAMETER ( PI = 3.1415926535D0 )
      PARAMETER ( RAD = PI/180. )
      INTEGER *4 I, J, K, NPTS, IERR
      REAL *4 XMIN, XMAX, YMIN, YMAX, XHI, XLO, YHI, YLO
      REAL *8 WAVE0, WAVE1, T, P, E(3), WAVE, INDEX
      REAL *4 XPLOT(300), YPLOT(300,3)
C
      P = 1013.
      T =  285.
      E(1) = 0.
      E(2) = 2.
      E(3) = 4.
      NPTS = 300
      WAVE0 = .450
      WAVE1 = .800
      YMIN = 1.E6
      YMAX = 0.
      DO I = 1, NPTS
          WAVE = WAVE0 + (WAVE1-WAVE0)*FLOAT(I-1)/FLOAT(NPTS)
          DO J = 1, 3
              CALL AINDEX ( P, T, E(J), WAVE, INDEX )
              YPLOT(I,J) = 1.E6*INDEX
              YMAX = MAX ( YMAX, YPLOT(I,J) )
              YMIN = MIN ( YMIN, YPLOT(I,J) )
          END DO
          XPLOT(I) = WAVE
      END DO
      XMIN = XPLOT(1)
      XMAX = XPLOT(NPTS)
      CALL PGBEGIN (0, '?', 1, 1 )
      CALL PGSCI   (5)
      CALL PGSCH   ( 1.5 )
      CALL PGENV   ( XMIN, XMAX, YMIN, YMAX, 0, 2 )
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, 'DISPERSION PLOT')
      CALL PGSCH   ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, 'WAVELENGTH (MICRONS)')
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, 'CHANGE IN DELAY (MICRONS)' )
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, 'FOR 1 METER OF AIR MISMATCH')
      CALL PGMTEXT ( 'R', 2.5, 0.5, 0.5, ' ' )
      CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, 'LINES ARE 0, 4 mb OF H2O')
      CALL PGSCI   ( 2 )
      CALL PGLINE ( NPTS, XPLOT, YPLOT(1,1) )
      CALL PGLINE ( NPTS, XPLOT, YPLOT(1,3) )
      CALL PGEND
  900 CONTINUE
      STOP
 1100 FORMAT ( 16I5 )
 1200 FORMAT ( 2F7.1, F7.3, F7.1, F10.6 )
      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
