CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SWPARAMETERS(PM_RATE,SPEED,MACH,ALPHA,SWDENSITY,VTH,CHIMIN,
     &                        NUM_SEXT_IN_AVERAGE, VSW_BEST, VTH_BEST)
C+
C
C FUNCTIONAL DESCRIPTION:
C
C    Reads a file called PM.SEXTUPLETS that contains records that correspond to
C    a variety of speeds and temperatures.  Each record contains 6 simulated PM
C    rates.  The subroutine finds the sextuplet that has the lowest CHI squared
C    when compared to 6 measured rate values.
C
C DESIGN:
C
C    The use of hard-coded values for VOLT, below, prohibits the use of this 
C    program for certain time periods, as shown from the short history below:
C
C       The initial Peq turn on occurred on December 27 1995.
C
C       Stepping was first activated on December 29 1995 at 01:23. The initial 
C       cycle was 180 120 80 60 40 25 (digital values).
C
C       On January 9 1996, the highest step was increased from 180 to 190.
C
C       On January 20 1996, at 17:30, the stepping sequence assumed below 
C       (digital values of 200 124 77 47 29 18) was set.
C
C       Between 19:23 on April 17 1996 and 20:23 on April 18 1996, various mode
C       changes were made.
C
C    Therefore, this program should not be used to try to simulate data from
C    before January 20 1996 17:30, or data from the interval beginning on
C    April 17 1996 at 19:23 and ending on April 18 1996 at 20:23.
C
C    This version fixes Nalpha/Nproton at 0.04; it is meant to read a sextuplet
C    file that only contains sextuplets that have Nalpha/Nproton = 0.04
C
C    This version averages over a range of CHI squared values.  It also passes
C    back the number of CHI squared values used and the VSW corresponding to
C    the best CHI.
C
C FORMAL ARGUMENTS:
C
C    {tbs}
C    PM_RATE:
C      I*4	  I	Array of length 6 containing the measured rate values.
C
C    SPEED:
C      R*4	  O	Upon return will contain 
C
C    MACH:
C      R*4	  O	Upon return will contain 
C
C    ALPHA:
C      R*4	  O	Upon return will contain 
C
C    SWDENSITY:
C      R*4	  O	Upon return will contain 
C
C    VTH:
C      R*4	  O	Upon return will contain 
C
C    CHIMIN:
C      R*4	  O	Upon return will contain 
C
C    NUM_SEXT_IN_AVERAGE:
C      I*4	  O	Upon return will contain
C
C    VSW_BEST:
C      R*4	  O	Upon return will contain 
C
C    VTH_BEST:
C      R*4	  O	Upon return will contain 
C
C INTERFACE:
C
!	INTERFACE
!	  SUBROUTINE SWPARAMETERS(PM_RATE,SPEED,MACH,ALPHA,SWDENSITY,VTH,CHIMIN,
!     &                        NUM_SEXT_IN_AVERAGE, VSW_BEST, VTH_BEST)
!	    INTEGER(KIND=4), INTENT(IN), DIMENSION(6) :: PM_RATE
!	    REAL(KIND=4), INTENT(OUT) :: SPEED,MACH,ALPHA,SWDENSITY,VTH,CHIMIN
!	    INTEGER, INTENT(OUT) :: NUM_SEXT_IN_AVERAGE
!	    REAL(KIND=4), INTENT(OUT) :: VSW_BEST, VTH_BEST
!	  END SUBROUTINE SWPARAMETERS
!	END INTERFACE
C
C IMPLICIT INPUTS:
C
C    The disk file DISKM1:[SOHO.TEST]PM.SEXTUPLETS is read into memory the first
C    time this routine is called, and is used on all subsequent calls.
C
C IMPLICIT OUTPUTS:
C
C    None
C
C SIDE EFFECTS:
C
C    See Implicit Inputs.
C
C COMMON BLOCKS:
C
C    None
C
C INCLUDE FILES:
C
C    None
C
C CREATED:
C
C    22-JAN-1997    John Paquette
C
C This version:   1.17   20-APR-1999
C
C DEVELOPMENT HISTORY:
C
C J. PAQUETTE	22-JAN-1997	get PM.SEXTUPLETS file from DISKM1:[SOHO.TEST]
C J. PAQUETTE	31-JAN-1997	Larger # of VSW values, hence larger # of sextuplets
C J. PAQUETTE	10-FEB-1997     Only 250 Vsw values
C J. PAQUETTE	10-FEB-1997     Different error, differetn CHI sq averaging
C J. PAQUETTE	10-FEB-1997     Last mod allowed error to be 0. Correct by adding 1. inside SQRT
C J. PAQUETTE	10-NOV-1997     New background - to deal with flare periods
C J. PAQUETTE   10-NOV-1997     Derive background from subroutine
C J. PAQUETTE   11-NOV-1997     Convert G_FLOAT to D_FLOAT
C J. PAQUETTE   20-NOV-1997     Background taken care of elsewhere
C J. PAQUETTE    4-DEC-1997     Renamed to SWPARAMETERS
C J. PAQUETTE    6-MAY-1998     Added VTH_BEST to argument list
C L. Bleau	18-MAR-1999	added interface block
C L. Bleau	15-APR-1999	added comments
C L. Bleau	20-APR-1999	added IMPLICIT NONE, declarations
C L. Bleau	10-NOV-1999	added comments
C
C-
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CDEC$ IDENT	'1.17'
      IMPLICIT NONE
      INCLUDE '($CVTDEF)/NOLIST'
      INTEGER NUM_SEXT
      PARAMETER (NUM_SEXT=250*25)
      INTEGER*4 STATUS, CVT$CONVERT_FLOAT, I, J, NUM_SEXT_IN_AVERAGE
      INTEGER*4 BESTINDEX, PM_RATE(6), FLAG
      REAL*4 SPEED,MACH,ALPHA,SWDENSITY, VTH, CHIMIN, VSW_BEST, VTH_BEST
      REAL*8 VOLT(6)/3.149, 1.962, 1.228, 0.758, 0.474, 0.302/, VSW(NUM_SEXT),
     &       MSW(NUM_SEXT), SIM_RATE(6,NUM_SEXT), CHI2, MEASURED_RATE(6), ERROR(6), 
     &       MEASURED_TOTAL/0.0/, DENSITY, BESTCHI, RATE(6), ALPHA_TO_PROTON(NUM_SEXT)
      REAL*8 BACKGND_G, BACKGND
      REAL*8 ChiArray(NUM_SEXT), SimTotArray(NUM_SEXT), DenArray(NUM_SEXT), MaxPMSR, MINPMSR
      REAL*8 AvgCntr, AvgSpeed, AvgMach, AvgAlpha, AvgDensity, AvgVth, Chi2UL
      LOGICAL FIRST/.TRUE./
C
      BESTCHI = 10000000.
      MEASURED_TOTAL = 0.0
      CHI2 = 0.0
C      SIM_TOTAL = 0.0
C
C Open output file
C
      IF(FIRST)THEN
        OPEN(UNIT=16, FORM='UNFORMATTED', TYPE='OLD', READONLY, RECORDTYPE='FIXED',
     &       NAME='DISKM1:[SOHO.TEST]PM.SEXTUPLETS', RECL=18)
        DO 5 I=1,NUM_SEXT
          READ(16) VSW(I),MSW(I),ALPHA_TO_PROTON(I),(SIM_RATE(J,I),J=1,6)
          SimTotArray(I) = SIM_RATE(1,I)
          DO 4 J=2, 6
            SimTotArray(I) = SimTotArray(I) + SIM_RATE(J,I)
   4      CONTINUE
   5    CONTINUE
        FIRST=.FALSE.
        CLOSE(16)          
      ENDIF
C
C Background taken care of elsewhere
C
      BACKGND = 0.0
C
C Subtract off background; find a total of all 6 rates to use for density normalization
C
      MaxPMSR = 0.0
      DO 14 I=1,6
        MEASURED_RATE(I) = PM_RATE(I)
        IF (MEASURED_RATE(I) .GT. MaxPMSR) MaxPMSR = MEASURED_RATE(I)
        MEASURED_RATE(I) = DMAX1(0.0D0, MEASURED_RATE(I) - BACKGND)
C        IF(MEASURED_RATE(I) .LT. 0.0)MEASURED_RATE(I) = 0.0
        MEASURED_TOTAL = MEASURED_TOTAL + MEASURED_RATE(I)
  14  CONTINUE
C
C Calculate errors    
C
      MaxPMSR = 0.2 * MaxPMSR
      DO 15 I=1,6
        ERROR(I) = 10.0*SQRT(MEASURED_RATE(I) + MaxPMSR + 1.0)
  15  CONTINUE
C
      DO 700 J=1,NUM_SEXT
C
C Normalize calculated density and rates using the actual rates.
C
C        DO 100 I =1,6
C          SIM_TOTAL = SIM_TOTAL + SIM_RATE(I,J)
C 100    CONTINUE
C        DENSITY = 10.0 * MEASURED_TOTAL/SIM_TOTAL
        DenArray(J) = 10.0 * MEASURED_TOTAL/SimTotArray(J)
        DO 420 I=1,6
C          RATE(I) = SIM_RATE(I,J) * MEASURED_TOTAL/SIM_TOTAL
          RATE(I) = SIM_RATE(I,J) * MEASURED_TOTAL/SimTotArray(J)
          CHI2 = CHI2 + ((RATE(I) - MEASURED_RATE(I))/ERROR(I))**2
 420    CONTINUE
C
C Divide our sum by the number of points to get a true chi squared.
C
        CHI2 = CHI2/6.0
        ChiArray(J) = CHI2
        IF(CHI2 .LT. BESTCHI)THEN
          BESTINDEX = J
          BESTCHI = CHI2
        ENDIF
C
        CHI2 = 0.0
C        SIM_TOTAL = 0.0
 700  CONTINUE 
C
C Normalize calculated density and rates using the actual rates.
C
C      DO 750 I =1,6
C        SIM_TOTAL = SIM_TOTAL + SIM_RATE(I,BESTINDEX)
C 750  CONTINUE
C      DENSITY = 10.0 * MEASURED_TOTAL/SIM_TOTAL
C      DO 770 I=1,6
C        RATE(I) = SIM_RATE(I,BESTINDEX) * MEASURED_TOTAL/SIM_TOTAL
C 770  CONTINUE
C
C average parameters over a chi^2 range
      AvgCntr = 0
      AvgSpeed = 0
      AvgMach = 0
      AvgAlpha = 0
      AvgDensity = 0
      AvgVth = 0
      Chi2UL = 0.5 + 1.6 * BESTCHI
C
      DO 710 J=1, NUM_SEXT
        IF (ChiArray(J) .LT. Chi2UL) THEN
          AvgCntr = AvgCntr + 1
          AvgSpeed = AvgSpeed + VSW(J)
          AvgMach = AvgMach + MSW(J)
          AvgAlpha = AvgAlpha + ALPHA_TO_PROTON(J)
          AvgDensity = AvgDensity + DenArray(J)
          AvgVth = AvgVth + VSW(J) / MSW(J)
        ENDIF
 710  CONTINUE
C
C Convert to single precision and pass back
C
C      SPEED = VSW(BESTINDEX)
C      MACH = MSW(BESTINDEX)
C      ALPHA = ALPHA_TO_PROTON(BESTINDEX)
C      SWDENSITY = DENSITY
C      VTH = SPEED / MACH
C      CHIMIN = BESTCHI
      SPEED = AvgSpeed / AvgCntr
      MACH = AvgMach / AvgCntr
      ALPHA = AvgAlpha / AvgCntr
      SWDENSITY = AvgDensity / AvgCntr
      VTH = AvgVth / AvgCntr
      NUM_SEXT_IN_AVERAGE = AvgCntr
      VSW_BEST = VSW(BESTINDEX)
      VTH_BEST = VSW(BESTINDEX)/MSW(BESTINDEX)
      CHIMIN = BESTCHI
C
C      WRITE(16, 800) VSW(BESTINDEX), MSW(BESTINDEX), ALPHA_TO_PROTON(BESTINDEX), DENSITY, RATE, BESTCHI
C 800  FORMAT(1X,F5.0,3X,F5.1,2X,F5.3,2X,F10.4,4X,6F7.0,4X,F7.0)
C
      RETURN
      END
