c Filename: mietable.f
c Author:   Mike Schubmehl
c Date(s):  07/11/2002 01/03/2002 06/25/2001 07/14/2000
c Purpose:  Calls Wiscombe's MIEV0 routine to generate tables of
c           intensity as a function of angle for three different
c           polarizations of incident light. It outputs tables for
c           angles ranging from 0 to 180 degrees, for size parameters
c           at a given interval between two user-specified endpoints.
c           Command line syntax is as follows:
c
c        mietable xmin xmax xstep nr [maxangle] [anglestep] [> outputfile]
c
c           xmin  - A real value greater than zero, and in practice
c                   greater than around 0.1. This is the first size
c                   parameter in the table.
c           xmax  - A real value greater than or equal to xmin, and
c                   in practice less than around 10,000 or so. This is
c                   the largest size parameter to appear in the table.
c           xstep - A positive real value giving the step size between
c                   consecutive size parameters in the table. 0.05
c                   should be the smallest needed for the experiment.
c           nr    - Relative index of refraction of scatterers to
c                   surrounding medium (e.g., latex spheres in water
c                   are 1.59/1.33 = 1.2).
C           maxangle (optional)
C                 - Maximum angle to output. Must be in the range 5--180.
C           anglestep (optional)
C                 - Spacing between angles. Must be in the range 0.1--5.
c
C           If the maxangle is specified, then the program outputs
C           from 0 to maxangle instead of from 0 to 180 degrees.
c           The output of the program goes to standart output, and can
c           thus be redirected with the command > outputfile to be
c           written to an output file.
c
c           Output is in the form:
c
c    *   0.500
c        0.00      0.26902789E-03      0.26902789E-03      0.26902789E-03
c        1.00      0.26902361E-03      0.26894236E-03      0.26898299E-03
c        2.00      0.26901066E-03      0.26868581E-03      0.26884823E-03
c       <angle>    <perp intens>       <par intens>        <unp intens>
c          .             .                   .                    .
c          .             .                   .                    .
c          .             .                   .                    .
c    *  <size parameter>
c          .             .                   .                    .
c          .             .                   .                    .
c          .             .                   .                    .
c
c Notes:    Requires miev0.f and ErrPack.f, Wiscombe's publicly
c           available routines, to generate mie tables. Note that the
c           include command is not in ANSI Fortran 77, but is supported
c           by the g77 compiler. If things aren't working, and you're
c           using a different compiler, see if it supports INCLUDE.

      INCLUDE 'ErrPack.f'
      INCLUDE 'miev0.f'

      PROGRAM MieTable

c     Declarations ****************************************************
        IMPLICIT NONE

        INTEGER i
        REAL theta
        CHARACTER *20 BUFFER
        REAL TEMP
        REAL XMIN, XMAX, XSTEP, NR, MAXOUTANGLE, ANGLESTEP
        REAL I1, I2

        REAL XX
        REAL XSAVE

        REAL XMU(400)
        COMPLEX CREFIN
        LOGICAL PERFCT
        REAL MIMCUT
        LOGICAL ANYANG
        INTEGER NUMANG
        INTEGER NMOM
        INTEGER IPOLZN
        INTEGER MOMDIM
        LOGICAL PRNT(2)

        REAL QEXT
        REAL QSCA
        REAL GQSC
        REAL PMOM(2,2)

        COMPLEX S1(181)
        COMPLEX S2(181)
        COMPLEX SFORW
        COMPLEX SBACK

        COMPLEX TFORW(2)
        COMPLEX TBACK(2)

        REAL SPIKE
c     END Declarations *************************************************

c     Main Program *****************************************************

c       First, check for proper number of arguments and print a usage
c         message if necessary:

        IF ((IARGC() .LT. 4) .OR. (IARGC() .GT. 6)) THEN
          WRITE(*, *)
     &     'Usage: mietable xmin xmax xstep nr [maxangle] [anglestep]'
          WRITE(*, *)
          WRITE(*, *) ' xmin  - Starting size parameter for table.'
          WRITE(*, *) ' xmax  - Largest size parameter for table.'
          WRITE(*, *) ' xstep - Spacing between size parameters.'
          WRITE(*, *)
     &     ' nr    - Relative index of scatterer to surrounding medium'
          WRITE(*, *)
     &     '         (e.g., latex spheres in water = 1.59/1.33 = 1.2).'
          WRITE(*, *)
     &     ' [maxangle]  - Largest angle to write (6 to 180).'
          WRITE(*, *)
     &     ' [anglestep] - Spacing between angles (0.5 to 5.0).'
          STOP
        END IF

c       Get arguments from command line, and put them into variables:

        CALL GETARG(1, BUFFER)
        READ(BUFFER, *) XMIN
        CALL GETARG(2, BUFFER)
        READ(BUFFER, *) XMAX
        CALL GETARG(3, BUFFER)
        READ(BUFFER, *) XSTEP
        CALL GETARG(4, BUFFER)
        READ(BUFFER, *) NR

C       Now, if there was a fifth/sixth argument, read max output angle
C         and/or angle step arguments. These are differentiated by range.

        MAXOUTANGLE = 180.0
        ANGLESTEP   = 1.0

        IF (IARGC() .GE. 5) THEN
          CALL GETARG(5, BUFFER)
          READ(BUFFER, *) TEMP
          IF ((TEMP .GE. 0.5) .AND. (TEMP .LE. 5.0)) THEN
            ANGLESTEP   = TEMP
            MAXOUTANGLE = 180
          ELSE
            MAXOUTANGLE = TEMP
            ANGLESTEP   = 1.0
          END IF
        END IF
        IF (IARGC() .EQ. 6) THEN
          CALL GETARG(6, BUFFER)
          READ(BUFFER, *) TEMP
          IF ((TEMP .GE. 0.5) .AND. (TEMP .LE. 5.0)) THEN
            ANGLESTEP   = TEMP
          ELSE
            MAXOUTANGLE = TEMP
          END IF
        END IF

C       Make sure the max output angle and angle step are valid.

        IF ((MAXOUTANGLE .LT. 6) .OR. (MAXOUTANGLE .GT. 180)) THEN
          WRITE(*,*)
     &     'Error: Invalid maximum output angle. Must be 6 to 180.'
          STOP
        END IF
        IF ((ANGLESTEP .LT. 0.5) .OR. (ANGLESTEP .GT. 5.0)) THEN
          WRITE(*,*)
     &     'Error: Invalid angle step. Must be 0.5 to 5.0.'
          STOP
        END IF

c       Setup other values for the call to MIEV0:

        CREFIN = NR
        PERFCT = .FALSE.
        MIMCUT = 1.0
        IF ((ABS(ANGLESTEP - 1.0) .GT. 0.01) .OR.
     &      (MAXOUTANGLE .NE. 180)) THEN
          ANYANG = .TRUE.
          NUMANG = INT(MAXOUTANGLE/ANGLESTEP) + 1
        ELSE
          ANYANG = .FALSE.
          NUMANG = 181
        ENDIF

        PRNT(1) = .FALSE.
        PRNT(2) = .FALSE.

        NMOM = 0
        IPOLZN = 0
        MOMDIM = 1

c       Generate array of angle cosines to be passed into MIEV0 (see
c       miev.doc for further details):

        theta = 0
        i = 1
        DO WHILE (theta <= MAXOUTANGLE)
          XMU(i) = COS(theta * ACOS(-1.0) / 180.0)
          theta = theta + ANGLESTEP
          i = i + 1
        END DO

c       Step over the range of size parameters specified by the user
c       and generate tables by calling MIEV0:

        XX = XMIN
        DO WHILE (XX <= XMAX)
          XSAVE = XX

          WRITE( *, '(''* '',F7.3)' ) XX

          CALL MIEV0(XX, CREFIN, PERFCT, MIMCUT, ANYANG, NUMANG, XMU,
     &               NMOM, IPOLZN, MOMDIM, PRNT, QEXT, QSCA, GQSC, PMOM,
     &               SFORW, SBACK, S1, S2, TFORW, TBACK, SPIKE)
          XX = XSAVE

c         Output tables to file in format described above:

          theta = 0
          i = 1
          DO WHILE (theta <= MAXOUTANGLE)
            I1 = REAL( S1( i ) )**2 + AIMAG( S1( i ) )**2
            I2 = REAL( S2( i ) )**2 + AIMAG( S2( i ) )**2
            WRITE( *, '( F6.2, 3E20.8)' )
     &              ACOS(XMU(i))*180./ACOS(-1.0),
     &              I1, I2, 0.5*(I1+I2)
            theta = theta + ANGLESTEP
            i = i + 1
          END DO

          XX = XX + XSTEP
        END DO

        STOP
      END PROGRAM MieTable
c     END Main Program *************************************************