c vlf section

	subroutine vlfsum0(vlf)
	common/trecord/id
c note: sa1(i),sa2(i) E-field (i=1,4), B-field (i=5,8)
	structure /vlfsum/
	   real sa1(8), sa2(8),sfr(8)
	   integer nsa1(8),nsa2(8),nsfr(8)
	   real*8 epoch 
	end structure
	record /vlfsum/ vlf
C**********************************************************************
        INTEGER*4 ID(324),IRAW(64)
        EQUIVALENCE (ID(34),IRAW(1))
C************************************************************************
	INTEGER L_K(4)/2,8,6,4/, M_K(4)/7,5,3,1/
        INTEGER L_K0(2)/7,3/,M_K0(2)/1,5/

	XL = ANT_LENGTH(ISTAT) !ANTENNA LENGTH OF DIPOLE ANTENNA
C
C	SPECTRUM ANALYSER (SA) DATA
C
        K=MOD(ID(100),4)
        L = L_K(K+1)
        M = M_K(K+1)
        VLF.SA1(L) = VLF.SA1(L) + VLF_CAL(IRAW( 6),L,XL,1)
        VLF.SA2(L) = VLF.SA2(L) + VLF_CAL(IRAW( 7),L,XL,2)
        VLF.SA1(M) = VLF.SA1(M) + VLF_CAL(IRAW(35),M,XL,1)
        VLF.SA2(M) = VLF.SA2(M) + VLF_CAL(IRAW(36),M,XL,2)

        VLF.NSA1(L) = VLF.NSA1(L) + 1
        VLF.NSA2(L) = VLF.NSA2(L) + 1
        VLF.NSA1(M) = VLF.NSA1(M) + 1
        VLF.NSA2(M) = VLF.NSA2(M) + 1
C
C       STEP FREQUENCY RECEIVER ( SFR ) DATA
C
        K=MOD(ID(100),2) + 1
	L = L_K0(K)
        M = M_K0(K)
        LL=L+1
        MM=M+1
	VLF.SFR( L) = VLF.SFR( L) + VLF_CAL(IRAW( 9), L,XL,3)
	VLF.SFR(LL) = VLF.SFR(LL) + VLF_CAL(IRAW( 8),LL,XL,3)
	VLF.SFR( M) = VLF.SFR( M) + VLF_CAL(IRAW(38), M,XL,3)
	VLF.SFR(MM) = VLF.SFR(MM) + VLF_CAL(IRAW(37),MM,XL,3)

	VLF.NSFR( L) = VLF.NSFR( L) + 1
	VLF.NSFR(LL) = VLF.NSFR(LL) + 1
	VLF.NSFR( M) = VLF.NSFR( M) + 1
	VLF.NSFR(MM) = VLF.NSFR(MM) + 1
	end   

	subroutine vlfprocess(vlf)
	common/trecord/id
	common/vlf_eng/v(10),iv(7)
c note: sa1(i),sa2(i) E-field (i=1,4), B-field (i=5,8)
	structure /vlf/
	   real sa1(8), sa2(8)
	   real sfr(8,2)
	   real*8 epoch 
	end structure
	record /vlf/ vlf
C**********************************************************************
        INTEGER*4 ID(324),IRAW(64)
        EQUIVALENCE (ID(34),IRAW(1))
C************************************************************************
	INTEGER L_K(4)/2,8,6,4/, M_K(4)/7,5,3,1/
        INTEGER L_K0(2)/7,3/,M_K0(2)/1,5/

	XL = ANT_LENGTH(ISTAT)
C
C	SPECTRUM ANALYSER (SA) DATA
C
        K=MOD(ID(100),4)
        L = L_K(K+1)
        M = M_K(K+1)
        VLF.SA1(L)=VLF_CAL(IRAW( 6),L,XL,1) !td=.45 s
        VLF.SA2(L)=VLF_CAL(IRAW( 7),L,XL,2) !td=.54 s
        VLF.SA1(M)=VLF_CAL(IRAW(35),M,XL,1) !td=3.06 s
        VLF.SA2(M)=VLF_CAL(IRAW(36),M,XL,2) !td=3.15 s
C
C       STEP FREQUENCY RECEIVER ( SFR ) DATA
C
        K=MOD(ID(100),2) + 1
        N = 1
        IF(MOD(ID(100),4) .GE. 2) N=2
	L = L_K0(K)
        M = M_K0(K)
        LL=L+1
        MM=M+1
	VLF.SFR(L,N)=VLF_CAL(IRAW(9),L,XL,3) !td=-.90
	VLF.SFR(LL,N)=VLF_CAL(IRAW(8),LL,XL,3) !td=.63
	VLF.SFR(M,N)=VLF_CAL(IRAW(38),M,XL,3)  !td=1.98
	VLF.SFR(MM,N)=VLF_CAL(IRAW(37),MM,XL,3) !td=3.06
	end   

	REAL*4 FUNCTION ANT_LENGTH(ISTAT)
	COMMON /VLF_ENG/V(10),IV(7)
C       CALCULATE ANTENNA LENGTH TIP TO TIP IN METERS
        X=IV(3)
	ISTAT = IV(3)
        X=X*0.0196+0.005
        IF (X.GT. 1.4) GO TO 305
        XL=17.0736-1.1774*X
        XL=-9.7273+XL*X
        GO TO 306
  305   XL=16.1819-0.98509*X
        XL=-8.8570+XL*X
  306   CONTINUE
	ANT_LENGTH = XL
	RETURN
	END

	real*4 FUNCTION VLF_CAL(IRAW,K,XL,IJUMP)
C************************************************************************
C       FREQUENCIES FOR THE VLF CHANNELS
        real*4 SFF(8)/13.3,17.8,23.7,31.1,42.2,56.2,100.,178./       !SFR KHZ
        real*4 SAF1(8)/1.78,5.62,17.8,56.2,1.78,5.62,17.8,56.2/      !SA-1 HZ
        real*4 SAF2(8)/178.,562.,1780.,5620.,178.,562.,1780.,5620./  !SA-2 HZ
C************************************************************************
C       BAND WIDTHS FOR VLF CHANNELS
        real*4 SFFBW(8)/878.,987.,1557.,1162.,1559.,2259.,4929.,5609./
        real*4 SA1BW(8)/2.80,1.51,13.25,26.46,3.13,1.98,13.28,24.08/
        real*4 SA2BW(8)/69.98,302.8,830.3,2356.,73.13,309.8,1034.,2532./
        real*4 
	1SA1C(8)/.5,.5,.5,.5,7.571E-05,2.466E-04,8.113E-04,2.266E-03/
        real*4 
	1SA2C(8)/.5,.5,.5,.5,5.508E-03,7.601E-03,7.880E-03,7.0096E-03/
C************************************************************************
	X=IRAW*0.0196+0.005			!compute log comp voltage
	IF (X.GT. 5.03) X=5.03
	IF( IJUMP .EQ. 2) GOTO 2
	IF( IJUMP .EQ. 3) GOTO 3

1	CALL SSA1(X,K,Y)                        !CONVERT TO V rms
	XX=Y/SA1C(k)
	IF( K .LE. 4) XX = XX/XL		!DIVIDE BY DIP ANT LENGTH
	VLF_CAL = XX*XX/SA1BW(K)		!CONVERT TO SPECTRAL DENSITY
	RETURN

2	CALL SSA2(X,K,Y)                        !CONVERT TO V rms
	XX=Y/SA2C(k)
	IF( K .LE. 4) XX = XX/XL		!DIVIDE BY DIP ANT LENGTH
	VLF_CAL = XX*XX/SA2BW(K)		!CONVERT TO SPECTRAL DENSITY
	RETURN

3	CALL SSFR(X,K,Y)                        !CONVERT TO V rms
	XX=Y/(xl/2)
	VLF_CAL = XX*XX/SFFBW(K)		!CONVERT TO SPECTRAL DENSITY
	RETURN
	END

C**********************************************************************
        SUBROUTINE SSA1(X,I,Y)
C**********************************************************************
C       ROUTINE TO CONVERT THE VLF LOG COMPRESSOR OUTPUT
C       TO THE CORRECTED INPUT VOLTAGES FOR THE SA-1 SYSTEM
C**********************************************************************
C       FREQUENCIES 1.78,5.62,17.8,56.2 HZ ELECTRIC AND MAGNETIC
C**********************************************************************
        DIMENSION B(8)
        DATA B/1.447,0.388,-0.006,-0.004,-0.133,-0.712,-0.746,-0.716/
        IF (X.LT. 4.601) GO TO 10
        Y=-5.03435+0.96335*X
        GO TO 20
 10     IF (X.LT. 4.192) GO TO 11
        Y=-6.32831+1.24456*X
        GO TO 20
 11     IF (X.LT. 3.599) GO TO 12
        Y=-4.18081+0.73229*X
        GO TO 20
 12     IF (X.LT. 3.081) GO TO 13
        Y=-5.62028+1.13227*X
        GO TO 20
 13     IF (X.LT. 2.524) GO TO 14
        Y=-4.33788+0.71605*X
        GO TO 20
 14     IF (X.LT. 1.988) GO TO 15
        Y=-5.35324+1.11837*X
        GO TO 20
 15     IF (X.LT. 1.401) GO TO 16
        Y=-4.57466+0.72671*X
        GO TO 20
 16     IF (X.LT. 1.004) GO TO 17
        Y=-5.30468+1.24766*X
        GO TO 20
 17     IF (X.LT. 0.414) GO TO 18
        Y=-4.88784+0.83264*X
        GO TO 20
 18     Y=-5.23333+1.66667*X
 20     Y=Y+B(I)
        Y=10.0**Y
        RETURN
        END

C**********************************************************************
        SUBROUTINE SSA2(X,I,Y)
C**********************************************************************
C       ROUTINE TO CONVERT THE VLF LOG COMPRESSOR OUTPUT
C       TO THE CORRECTED INPUT VOLTAGES FOR THE SA-2 SYSTEM
C**********************************************************************
C       FREQUENCIES 178,562,1.78K,5.62K HZ ELECTRIC AND MAGNETIC
C**********************************************************************
        DIMENSION B(8)
        DATA B/0.042,-0.004,-0.0055,0.093,-0.682,-0.75,-0.75,-0.6365/
        IF (X.LT. 4.641) GO TO 10
        Y=-5.27+X
        IF (I.EQ.4) Y=Y+0.083
        GO TO 20
 10     IF (X.LT. 4.155) GO TO 11
        Y=-6.12912+1.18510*X
        IF (I.EQ.4) Y=Y+0.048
        GO TO 20
 11     IF (X.LT. 3.623) GO TO 12
        Y=-4.09281+0.695018*X
        IF (I.EQ.4) Y=Y-0.022
        GO TO 20
 12     IF (X.LT. 3.087) GO TO 13
        Y=-5.67735+1.13235*X
        GO TO 20
 13     IF (X.LT. 2.527) GO TO 14
        Y=-4.36360+0.706838*X
        GO TO 20
 14     IF (X.LT. 2.007) GO TO 15
        Y=-5.44755+1.13578*X
        GO TO 20
 15     IF (X.LT. 1.391) GO TO 16
        Y=-4.60093+0.713992*X
        GO TO 20
 16     IF (X.LT. 0.943) GO TO 17
        Y=-5.32231+1.23250*X
        GO TO 20
 17     IF (X.LT. 0.486) GO TO 18
        Y=-4.84802+0.729694*X
        GO TO 20
 18     IF (X.LT. 0.264) GO TO 19
        Y=-5.03333+1.11111*X
        GO TO 20
 19     Y=-5.40+2.5*X
 20     Y=Y+B(I)
        Y=10.0**Y
        RETURN
        END

C************************************************************************
        SUBROUTINE SSFR(X,I,Y)
C************************************************************************
C       ROUTINE TO CONVERT THE VLF LOG COMPRESSOR OUTPUT
C       TO THE CORRECTED INPUT VOLTAGES FOR THE SFR SYSTEM
C************************************************************************
        DIMENSION B(8),A(8),C(8)
        DATA B/0.102,0.0188,-0.0069,-0.0315,0.0519,0.0169,0.0288,0.1383/
        DATA A/0.1452,0.0545,0.01815,0.0,0.0726,0.0363,0.0363,0.2087/
        DATA C/0.304,0.0967,0.0975,-0.0022,0.0553,-0.0029,0.0471,0.4526/
        IF (X.LT. 4.108) GO TO 10
        Y=-4.62046+0.90759*X
        Y=Y+A(I)
        IF (Y.LT. -0.345) GO TO 21
        Y=-11.95+2.5*X+C(I)
        GO TO 21
 10     IF (X.LT. 3.820) GO TO 11
        Y=-7.41225+1.58727*X
        GO TO 20
 11     IF (X.LT. 3.254) GO TO 12
        Y=-4.50012+0.82484*X
        GO TO 20
 12     IF (X.LT. 2.917) GO TO 13
        Y=-6.88878+1.55885*X
        GO TO 20
 13     IF (X.LT. 2.335) GO TO 14
        Y=-4.67732+0.80070*X
        GO TO 20
 14     IF (X.LT. 1.972) GO TO 15
        Y=-6.24591+1.47253*X
        GO TO 20
 15     IF (X.LT. 1.406) GO TO 16
        Y=-4.95012+0.81568*X
        GO TO 20
 16     IF (X.LT. 1.038) GO TO 17
        Y=-5.85304+1.45795*X
        GO TO 20
 17     IF (X.LT. 0.503) GO TO 18
        Y=-5.19802+0.82691*X
        GO TO 20
 18     Y=-7.7513E-06+4.82696E-05*X
c        Y=LOG10(Y)
	return
 20     Y=Y+B(I)
 21     Y=10.**Y
        RETURN
        END

c magnetometer section
	SUBROUTINE MAGPROCESS(mag,IREC)
	structure /mag/
	   real*8 epoch
	   real*4 b(3),colat(3),long(3)	!gsm coord (magnitude,colat,long)
	   real*4 bx(3),by(3),bz(3)	!Sch cal., s/c coord.
	   byte range
	end structure
	record /mag/mag
C**********************************************************************
        INTEGER*4 ID(324),IRAW(64)
        REAL*8 CD(21),EW(2),EM(9),EI(9),CD0(21),cd1(21)
C**********************************************************************
        real*4 BX(3),BY(3),BZ(3),BB(3)
C**********************************************************************
        EQUIVALENCE (ID(34),IRAW(1)),(ID(104),EW(1)),(ID(130),CD(1))
        EQUIVALENCE (ID(242),EM(1)),(ID(260),EI(1))
	common/trecord/id
C************************************************************************
        IMAG = (ID(12) .AND. 2**16 )/2**16      !ON/OFF
c	PRINT *,' IMAG=',IMAG
        IF (IMAG.EQ. 0)goto 9898               !MAGNETOMETER NOT ON
        MAGRANGE =ID(12)-(ID(12)/4)*4                   !RANGE 0 1 2 3
	mag.range = magrange
        JC=MOD(ID(101),1024)                    !S/C SEQUENCE NUMBER
	if( (jc .eq. 0) .or. (jc .eq. 1) )Then
	   print *,' Magnetometer calib tone encountered'
	    goto 9898 !remove calib tone
	endif
c        J=MOD (ID(100),32)+1                    !FRAME NUMBER
c        IF (J.EQ. 6) SENTEMP=EW(2)
c        IF (J.EQ. 7) ELETEMP=EW(2)

	do i = 1,21
	   CD0(i)=cd(i)
	enddo
C
C       DESPIN CONVERTED MAGNETOMETER DATA
C
        XI=360.D0/CD(10)
        IF (CD(12).GE. 0.D0) GO TO 205          !FOR DESPUN DATA
        YI=irec-1                             !FOR NON-DESPUN DATA
        CD0(12)=XI*YI*5.76D0                     !ASSUME START = 0.
 205    CONTINUE
       CD0(1)=MOD(CD(1)+CD0(12)+1.62D0*XI,360.D0)
       CD0(4)=MOD(CD(4)+CD0(12)+3.51D0*XI,360.D0)
       CD0(7)=MOD(CD(7)+CD0(12)+5.40D0*XI,360.D0)
C
       CALL ANGLE (EM,EI,CD0)                   ! TRANSFORM TO GSM
C
C       PULL OUT RAW MAGNETOMETER DATA
C
        BX(1)=IRAW(19)
        BY(1)=IRAW(20)
        BZ(1)=IRAW(21)
        BX(2)=IRAW(40)
        BY(2)=IRAW(41)
        BZ(2)=IRAW(42)
        BX(3)=IRAW(61)
        BY(3)=IRAW(62)
        BZ(3)=IRAW(63)
c	print *,' mag raw =',bx(1),by(1),bz(1)
	btt = 257
	do i = 1,3
	if( bx(i) .ge. btt) goto 9898
	if( by(i) .ge. btt) goto 9898
	if( bz(i) .ge. btt) goto 9898
	enddo

 	DO K=1,3
	   BX(k)=BX(k)*0.0196+0.005            !CONVERT TO VOLTS
	   BY(k)=BY(k)*0.0196+0.005
	   BZ(k)=BZ(k)*0.0196+0.005
	ENDDO

        IF (MAGRANGE.EQ. 3) GO TO 260
        IF (MAGRANGE.EQ. 2) GO TO 250
        IF (MAGRANGE.EQ. 1) GO TO 240

	DO K=1,3
   	   BX(k)=BX(k)*59.996-149.962+3.685    !RANGE 1
	   BY(k)=BY(k)*60.200-150.438-0.908    !SCHONSTEDT
	   BZ(k)=BZ(k)*60.186-150.363          !CALIBRATIONS
	ENDDO
        GO TO 264
  240   CONTINUE
        DO K=1,3
	   BX(k)=BX(k)*180.565-451.345+3.685   !RANGE 2
	   BY(k)=BY(k)*179.854-449.488-0.908   !SCHONSTEDT
	   BZ(k)=BZ(k)*180.883-452.133         !CALIBRATIONS
	ENDDO
        GO TO 264
  250   CONTINUE
	DO K=1,3
	   BX(k)=BX(k)*597.133-1493.403+3.685  !RANGE 3
	   BY(k)=BY(k)*599.672-1498.909-0.908  !SCHONSTEDT
	   BZ(k)=BZ(k)*601.870-1503.170        !CALIBRATIONS
	ENDDO
        GO TO 264
  260   CONTINUE
	DO K=1,3
	   BX(k)=BX(k)*09977.3-24937.8+3.685   !RANGE 4
	   BY(k)=BY(k)*10003.7-24994.3-0.908   !SCHONSTEDT
	   BZ(k)=BZ(k)*10042.4-25100.9         !CALIBRATIONS
	ENDDO
  264   CONTINUE
	do k = 1,3
	   mag.bx(k)=bx(k)
	   mag.by(k)=by(k)
	   mag.bz(k)=bz(k)
	   mag.colat(k) =cd0((k-1)*3 + 1)
	   mag.long(k)=cd0((k-1)*3 + 2)
	   mag.b(k)    =cd0((k-1)*3 + 3)
	enddo
C 	TRANFORM SCH. CAL. B FIELD INTO GSM COORD.
c	pi = 3.1415927
c        arad = 180./pi
c        DO K=1,3
c	   brho =BX(k)*BX(k) + BY(k)*BY(k) 
c	   X = brho + BZ(k)*BZ(k)
c	   brho = sqrt(brho)
c	   BB(k)=SQRT(X)
c	   theta = acos(bz(k)/bb(K))*arad
c	   phi = acos(bx(k)/brho) 
c	   if( by(k) .lt. 0) phi = 2*pi - phi
c	   phi = phi*arad
c	   CD1((k-1)*3 + 1) = phi
c	   CD1((k-1)*3 + 2) = theta
c	   CD1((k-1)*3 + 3) = bb(k)
c	ENDDO
c	CD1(1)=MOD( CD(1)+CD0(12)+1.62D0*XI, 360.D0)
c	CD1(4)=MOD( CD(4)+CD0(12)+3.51D0*XI, 360.D0)
c	CD1(7)=MOD( CD(7)+CD0(12)+5.40D0*XI, 360.D0)
c	CALL ANGLE (EM,EI,CD1)                   ! TRANSFORM TO GSM
	RETURN
9898	do i = 1,3
	   mag.bx(i)=-1
	   mag.by(i)=-1
	   mag.bz(i)=-1
	   mag.b(i)=-1
	   mag.colat(i)=-1
	   mag.long(i)=-1
	enddo
	   mag.epoch =-1
	print *,' magnetometer record skipped, irec=',irec
	END


c engineering routines
	SUBROUTINE ENG_LOAD(IN_UNIT,ILOAD)
	COMMON /VLF_ENG/V(10),IV(7) !SPECIFIC TO VLF INSTRUMENT
        DIMENSION IW24(16),IW25(32)
	INTEGER*2 DUMMY
        INTEGER*4 ID(324),IRAW(64)
        REAL*8 TRJ(77),CD(21),EW(2),WORD24(16),WORD25(32)
	real*8 Z,ZZ,XI,YY,YI,Y1,Y2,Y3,EM(9),EI(9)
C**********************************************************************
        EQUIVALENCE (ID(34),IRAW(1)),(ID(104),EW(1)),(ID(130),CD(1))
        EQUIVALENCE (ID(180),TRJ(1)),(ID(242),EM(1)),(ID(260),EI(1))
C************************************************************************
        J=MOD(ID(100),64)+1
        V(1)= (ID(12) .AND. 2**8 )/2**8	        !BAND WIDTH
        V(2)= (ID(12) .AND. 2**4 )/2**4         !MODE
        V(3)= (ID(12) .AND. 4)/4                !V DIFF GAIN
        IF (J.EQ.  9) V(4)=EW(2)                !TEMPERATURE
        IF (J.EQ.  9) IV(1)=IRAW(26)
        IF (J.EQ. 22) V(5)=EW(2)                !BOOM MOTOR TEMPERATURE
        IF (J.EQ. 22) IV(2)=IRAW(26)
        IF (J.EQ. 25) V(6)=EW(2)                !ANTENNA LENGTH WRONG
        IF (J.EQ. 25) IV(3)=IRAW(26)
        IF (J.EQ. 26) V(7)=EW(2)                !LOW VOLTAGE MONITOR
        IF (J.EQ. 26) IV(4)=IRAW(26)
        IF (J.EQ. 27) V(8)=EW(2)                !POWER SUPPLY STATUS
        IF (J.EQ. 27) IV(5)=IRAW(26)
        IF (J.EQ. 28) V(9)=EW(2)                !ELECTRIC ANTENNA LIMIT
        IF (J.EQ. 28) IV(6)=IRAW(26)
        IF (J.EQ. 29) V(10)=EW(2)               !TELESCOPIC BOOM LIMIT
        IF (J.EQ. 29) IV(7)=IRAW(26)

	IF(ILOAD .EQ. 0)THEN
	   PRINT *,'INITIALIZING ENGINEERING DATA'
	   DO ILOAD = 1,64
10	      READ(IN_UNIT,END=500) dummy,ID
	      if( cd(10) .le. 0.)then
	          print *,' cd(10)=',cd(10),'skip record'
	          goto 10
	      endif	
	        J=MOD(ID(100),64)+1
	        IF (J.EQ.  9) V(4)=EW(2)                !TEMPERATURE
	        IF (J.EQ.  9) IV(1)=IRAW(26)
	        IF (J.EQ. 22) V(5)=EW(2)                !BOOM MOTOR TEMPERATURE
	        IF (J.EQ. 22) IV(2)=IRAW(26)
	        IF (J.EQ. 25) V(6)=EW(2)                !ANTENNA LENGTH WRONG
	        IF (J.EQ. 25) IV(3)=IRAW(26)
	        IF (J.EQ. 26) V(7)=EW(2)                !LOW VOLTAGE MONITOR
	        IF (J.EQ. 26) IV(4)=IRAW(26)
	        IF (J.EQ. 27) V(8)=EW(2)                !POWER SUPPLY STATUS
	        IF (J.EQ. 27) IV(5)=IRAW(26)
	        IF (J.EQ. 28) V(9)=EW(2)                !ELECTRIC ANTENNA LIMIT
	        IF (J.EQ. 28) IV(6)=IRAW(26)
	        IF (J.EQ. 29) V(10)=EW(2)               !TELESCOPIC BOOM LIMIT
	        IF (J.EQ. 29) IV(7)=IRAW(26)
	    ENDDO
	    REWIND IN_UNIT
	    READ(IN_UNIT,END=500) dummy,ID !RESET TO FIRST RECORD
	    PRINT *,' ENGINEERING DATA INITIALIZED'
	ENDIF
	RETURN
500     PRINT *,'EOF ENCOUNTERED'
	STOP
	END

c time routines