PROGRAM saocfit
C
C	This program is designed to be used with the optical constants
C  from Tisdale et al., 1998.  The paper reports optical constants for
C  sulfuric acid solutions between 45 and 80 wt% H2SO4 at 215 K.  This
C  program allows one to obtain optical constants for other values of
C  wt% H2SO4 using fits of the optical constants for each frequency, 
C  which were also reported (in abbreviated form) in the paper.
C	There are two ways of using the program.  The program will ask
C  you to decide between the two.  One way is to have a file prepared
C  with the frequencies for which you want optical constants.  The 
C  program will prompt you for the filename of this file if you choose
C  this option.  The other is to specify a starting and stopping 
C  frequency values, as well as a step value.  For example, you could
C  ask for optical constants between 1000 and 2000 cm-1, every 25 cm-1.
C  IMPORTANT:  You must choose values between 500 and 7000 cm-1 for
C  the starting and stopping frequencies.  If you have frequencies
C  above 3650 cm-1, there will be zeros for the k values, as we were
C  unable to measure k above this frequency.
C	The number of frequencies in the data is 1685.  The maximum
C  number of frequency points to fit is 5000, though this is arbitrary
C  and can be set higher (edit the DIMENSION statement) if desired.
C	The citation for the source of the data is Tisdale et al.,
C  "Infrared optical constants of low temperature H2SO4 solutions
C  representative of stratospheric sulfate aerosols," J. Geophys.
C  Res., ?, ?-?, 1998.
C
C	Variable list:
C 
C  freqfile     filename where frequencies are stored (User Input)
C  outfile      filename of output file (UI) 
C  fstart	      start frequency (UI)
C  fstop        stop frequency (UI)
C  fstep        step for frequency (UI)
C  nfreqs       number of frequencies for output
C  mfreq        frequencies of measured optical constant fit parameters
C  cfreq        frequencies optical constants are desired for
C  p3,p2,p1,p0  parameters for fitting n
C  q3,q2,q1,q0  parameters for fitting k
C  nh,nl        n values for bounding frequencies in interpolation
C  kh,kl        k values for bounding frequencies in interpolation
C  ksa          output values of k
C  nsa          output values of n
C
      IMPLICIT REAL (A-H,O-Z)
      REAL ksa,nsa,mfreq,nslope,ninter,kslope,kinter,nh,nl,kh,kl
      CHARACTER*12 outfile, freqfile, dummy
      DIMENSION mfreq(1685), p3(1685), p2(1685), p1(1685), p0(1685),
     $  q3(1685), q2(1685), q1(1685), q0(1685), ksa(20000), nsa(20000),
     $  cfreq(20000)
C
C  Read in the fit parameters from the file fitparam.csv
C
	OPEN(10, file='fitparam.txt', status='old')
      DO 90 k=1,1685
	  IF  (k.lt.818) THEN
	    READ(10,*) mfreq(k),p3(k),p2(k),p1(k),p0(k),
     $      q3(k),q2(k),q1(k),q0(k)
        END IF
	  IF (k.ge.818) THEN
          READ(10,*) mfreq(k),p3(k),p2(k),p1(k),p0(k)
        END IF
90    CONTINUE
      CLOSE(10)
C
C  Prompt the user for information
C
4     CONTINUE
      PRINT 5
5     FORMAT (1X, 'What weight % H2SO4? (Ex: type 50 for 50 wt%)')
      READ *, wtp
	IF (wtp.gt.80. .or. wtp.lt.45.) THEN
	  PRINT *, 'Must be between 45 and 80 wt% H2SO4'
        GOTO 4
      END IF
	wtp=wtp/100.
      PRINT 10
10    FORMAT (1X, 'Do you have a file with frequencies set up?')
      READ *, dummy
	jflag=0
	IF (dummy(:1).eq.'y'.or.dummy(:1).eq.'Y') THEN
        jflag = 1
	  PRINT 20
20      FORMAT (1X, 'What is the filename with the frequencies?')
          READ *, freqfile
      ELSE
29      CONTINUE
	  PRINT 30
30      FORMAT (1X, 'Enter the starting frequency in wavenumbers')
        READ *, fstart
        IF (fstart.lt.499.4792 .or. fstart.gt.6996.565) THEN
          PRINT *, 'Must be between 499.4792 and 6996.565 cm-1'
          GOTO 29
        END IF
39      CONTINUE
	  PRINT 40
40      FORMAT (1X, 'Enter the ending frequency in wavenumbers')
	  READ *, fstop
        IF (fstop.lt.499.4792 .or. fstop.gt.6996.565) THEN
          PRINT *, 'Must be between 499.4792 and 6996.565 cm-1'
          GOTO 39
        END IF
	  PRINT 50
50      FORMAT (1X, 'Enter the step for frequency in wavenumbers')
	  READ *, fstep
      END IF
      PRINT 60
60    FORMAT (1X, 'What output filename do you want to use?')
      READ *, outfile
C
C  Set up the frequency array.  If jflag is 1, read from freqfile.
C  If jflag is 0, set it up by calculation. There are provisions
C  set up to account for whether or not the frequencies in the
C  calculation method are input in increasing or decreasing order.
C  These are the lines "fstep=-1*fstep" and the two if then loops
C  which check to see if the final frequency has been reached
C  (these loops count the number of frequencies and then send 
C  the program to label 150).
C
      IF (jflag .eq. 1) THEN
        OPEN(20, file=freqfile, status='old')
        DO 100 j=1,20000
99        CONTINUE
          READ(20,*,END=101) cfreq(j)
          IF (cfreq(j).lt.499.4792 .or. cfreq(j).gt.6996.565) THEN
	      GOTO 99
          END IF
          nfreqs=j
100     CONTINUE
101     CONTINUE
        CLOSE(20)
      END IF
      IF (jflag .eq. 0) THEN
        IF (fstart.gt.fstop.and.fstep.gt.0.) THEN
          fstep=-1.*fstep
        END IF
        j=0
125     CONTINUE
        j=j+1
        cfreq(j) = fstart + (j-1) * fstep
        IF (cfreq(j) .gt. fstop .and. fstep .gt. 0.) THEN
          nfreqs = j-1
          GOTO 150
        END IF
        IF (cfreq(j) .lt. fstop .and. fstep .lt. 0.) THEN
          nfreqs = j-1
          GOTO 150
        END IF
        GOTO 125
      END IF
150	CONTINUE
      PRINT*, nfreqs
C
C  Calculate indices.  The first IF..THEN checks to see if a frequency
C  exactly matches a measured frequency.  Calculation is then simple.
C  the second IF..THEN checks to make sure that the value of m
C  is high enough so that mfreq(m) is greater than cfreq(i).
C  This forces mfreq(m) and mfreq(m-1) to bound cfreq(i).  If not, m
C  is incremented until the condition is met.  Then, the indices can
C  be calculated using a linear interpolation.  The value of k is 
C  set to zero in the final IF..THEN if cfreq(i) is greater than the
C  last point in the dataset for k, which is 3647.711 cm-1.
C
      DO 200 i=1,nfreqs
        m=1
        IF (cfreq(i) .eq. mfreq(m)) THEN
          nsa(i) = p3(m)*wtp**3 + p2(m)*wtp**2 + p1(m)*wtp + p0(m)
          ksa(i) = q3(m)*wtp**3 + q2(m)*wtp**2 + q1(m)*wtp + q0(m)
        END IF
300     IF (cfreq(i) .gt. mfreq(m)) THEN
          m=m+1
          GO TO 300
        END IF
        nh = p3(m)*wtp**3 + p2(m)*wtp**2 + p1(m)*wtp + p0(m)
        nl = p3(m-1)*wtp**3 + p2(m-1)*wtp**2 + p1(m-1)*wtp + p0(m-1)
        nslope = (nh-nl)/(mfreq(m)-mfreq(m-1))
        ninter = nh-nslope*mfreq(m)
        nsa(i) = nslope * cfreq(i) + ninter
        IF (cfreq(i) .lt. 3647.711) THEN
          kh = q3(m)*wtp**3 + q2(m)*wtp**2 + q1(m)*wtp + q0(m)
          kl = q3(m-1)*wtp**3 + q2(m-1)*wtp**2 + q1(m-1)*wtp + q0(m-1)
          kslope = (kh-kl)/(mfreq(m)-mfreq(m-1))
          kinter = kh-kslope*mfreq(m)
          ksa(i) = kslope * cfreq(i) + kinter
        ELSE
          ksa(i) = 0
        END IF
200   CONTINUE
C
C  Save indices into outfile.	 If you want to eliminate the title
C  lines of output, comment out the first two write statements.
C
      OPEN(30, file=outfile, status='new')
      WRITE(30,*) 'Indices for ', wtp*100
	WRITE(30,*) 'Frequency ', 'n ', 'k '
      DO 400 i=1,nfreqs
        WRITE(30,*) cfreq(i), nsa(i), ksa(i)
400   CONTINUE
      CLOSE(30)
      END