CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	SUBROUTINE HEFFALUMP_H(PM_RATE, CURRENTTIME, SPEED, MACH, DEN)
C
C Hybrid of Hefti of Scott
C 17-Jan-97 modified to match Hefti's 1996/11/11 code
C
C This version:   1.9    4-DEC-1997
C
C  Author       Date            Description of Change 
C  --------     -----------     --------------------- 
C  J. PAQUETTE  17-FEB-1997     Corrected FORTRAN density to match original density in 
C                               Hefti's C code
C  J. PAQUETTE  18-FEB-1997     Removed deadtime correction(ours will be used instead)
C  J. PAQUETTE  18-FEB-1997     Change slope and offset for V
C  J. PAQUETTE  18-FEB-1997     Apply our time corrections in place of Hefti
C  J. PAQUETTE  18-FEB-1997     New derivation of Vth
C  J. PAQUETTE  24-FEB-1997     Name change
C  J. PAQUETTE	10-NOV-1997	Modified to use variable background
C  J. PAQUETTE  10-NOV-1997     Derive background from subroutine
C  J. PAQUETTE  20-NOV-1997     Background taken care of elsewhere
C  J. PAQUETTE   4-DEC-1997     Renamed to HEFFALUMP_H
C 
C
CDEC$   IDENT	'1.9'
	PARAMETER EOQRES_R = 6.0D-1
	PARAMETER SLP_V = 4.6901D2
	PARAMETER OFF_V = 2.2249D1
	PARAMETER SLP_N = 4.3041D-2
	PARAMETER OFF_N = -8.965D-1
	PARAMETER SLP_T = 3.219D2
	PARAMETER OFF_T = -1.04D1
	PARAMETER SEC06MAR962110 = 1204837800.D0
	PARAMETER SEC10JUN961855 = 1213124100.D0
	PARAMETER SEC01JAN960000 = 1199145600.D0
        INTEGER*4 PM_RATE(6), FLAG
	INTEGER I, M
	REAL*8 BACKGND
	REAL*8 PMRATE(6), SQRTPMEQ(6), SQRPMEQ1, MEASURED_TOTAL
	REAL*8 PMEQ(6)/0.302,0.474,0.758,1.228,1.962,3.149/
	REAL*8 SPEED, MACH, DEN, MOMENTS(2), CURRENTTIME
	REAL*8 Instr_Vel(6) ! velocity associated with PMEQ
	REAL*8 Sig_Instr(6) ! Delta_PMEQ expressed as velocity
	REAL*8 Sig_Beam     ! Vth^2
	REAL*8 Sig_Tmp2     ! Sig_Instr^2 + Sig_Beam^2
	REAL*8 Den_Cor      ! correction of density due to Delta_PMEQ
C
	MEASURED_TOTAL = 0.0
	Sig_Beam = 0.0
	Sig_Tmp2 = 0.0
	Den_Cor = 0.0
	MOMENTS(1) = 0.0
	MOMENTS(2) = 0.0
C
C Convert rates to REAL*8; Calculate square root of E/Q voltages
C
	DO 10 I=1, 6
          PMRATE(7-I) = PM_RATE(I)
	  SQRTPMEQ(I) = SQRT(PMEQ(I))
	  Instr_Vel(I) = SLP_V * SQRTPMEQ(I)
	  Sig_Instr(I) = 0.5 * EOQRES_R * Instr_Vel(I)
10	CONTINUE
C
C Background taken care of elsewhere
C
        BACKGND = 0.0
C
C Get measured rates, correct for E/Q window
C	 
	SQRPMEQ1 = PMEQ(1)*PMEQ(1)
	DO 20 I=1, 6
C
C Background (Ipavich)
C
	  PMRATE(I) = DMAX1(0.0D0, PMRATE(I) - BACKGND)
C
C Total for density
C
	  MEASURED_TOTAL = MEASURED_TOTAL + PMRATE(I)
C
C Account for E/Q window for moments
C
!	  PMRATE(I) = PMRATE(I)*SQRPMEQ1/(PMEQ(I)*PMEQ(I))
20	CONTINUE
C
	CALL HEFMOM(SQRTPMEQ, PMRATE, MOMENTS)
	SPEED = MOMENTS(1)
	MACH  = SQRT(MOMENTS(2))
	DEN   = MEASURED_TOTAL
C
C Apply corrections (in any case)
C	
	SPEED = SPEED * SLP_V + OFF_V
        CALL CORRECT_VTH_H(MACH,MOMENTS(2))
	IF (SPEED .GT. 0.0) THEN
	  DEN = DEN/SPEED
	ELSE
	  DEN = 0.0
	ENDIF
	
	IF (EOQRES_R .GT. 0.0) THEN
	  Sig_Beam = MACH
	  Den_Cor = 0.0
	  
	    DO 323 M = 1, 6
	     Sig_Tmp2 = Sig_Instr(M) * Sig_Instr(M) + Sig_Beam * Sig_Beam
	     Den_Cor = Den_Cor + Sig_Instr(M) / SQRT(Sig_Tmp2) * EXP(-0.5 * (Instr_Vel(M) - SPEED)**2 / Sig_Tmp2)
323         CONTINUE
          DEN = DEN/Den_Cor
        ENDIF
C
        DEN = DEN * SLP_N + OFF_N
C
C Renormalize density if time is before Mar 6 1996 at 2110, when bias changed (PM MCP bias=level 150
C before 6 MAR 96, 155 thereafter). Since the bias increased, and we are normalizing the early data 
C to the data after Mar 6, the density from before Mar 6 must be increased
C
	IF (CURRENTTIME .LT. SEC06MAR962110) DEN = DEN * 1.28897
C
C On Jun 10, the threshold was changed. Since we are normalizing the data after Jun 10 to the data 
C in the interval from Mar 6 to Jun 10, after the threshold change, density after Jun 10 must be 
C increased.
C
	IF (CURRENTTIME .GE. SEC10JUN961855) DEN = DEN * 0.620065
C
C Remove decay in density caused by loss of MCP gain. Since the decay was at first exponential, 
C but then leveled off, the function chosen to fit the decay was of the form A+B*EXP(-T/C)
C
        DEN = DEN/(0.42935 + 1.6509*EXP(-(CURRENTTIME - SEC01JAN960000)/5.9173E6))
C
        IF (DEN .LT. 0.0) DEN = 0.0
        IF (SPEED .LT. 100.0) SPEED = 100.0
        IF (MACH .LT. 0.0) MACH = 0.0
        RETURN
        END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	SUBROUTINE HEFMOM(X, Y, MOMENTS)
C HEFMOM ROUTINE TO CALCULATE MOMENTS
	REAL*8 X(6), Y(6), MOMENTS(2)
	INTEGER*4  I, J
	REAL*8 TOTY, MEAN, VAR

	TOTY = Y(1)
	DO 30 I=2, 6
	  TOTY = TOTY + Y(I)
30	CONTINUE

	MOMENTS(1) = 0.0
	MOMENTS(2) = 0.0
	IF (TOTY .GT. 0.0) THEN
	  MEAN = X(1)*Y(1)
	  DO 40 J = 2, 6
	    MEAN = MEAN + X(J)*Y(J)
40	  CONTINUE
	  MEAN = MEAN / TOTY
	  VAR  = (X(1)-MEAN)*(X(1)-MEAN)*Y(1)
	  DO 50 J = 2, 6
	    VAR = VAR + (X(J)-MEAN)*(X(J)-MEAN)*Y(J)
50	  CONTINUE
	  MOMENTS(1) = MEAN
	  MOMENTS(2) = VAR / TOTY
	ENDIF
	RETURN
	END
