PROGRAM SUMMHOURS
C
C PROGRAM MAKES HOUR AVERAGES OF HIGH RESOLUTION SUMMARY DATA
C WRITTEN BY S. KRAMER, CODE 692, 07/15/94
C
      PARAMETER(IARR=1000)
      INTEGER*2 OLDTIME(6),NEWTIME(6),DELTA(6)
      INTEGER*4 MAGCNT(4),SEDRCNT(4),IN(35)
      LOGICAL*1 STASH,FIRST
      REAL*4 MAGSUM(4),MAGAVG(4),MAG(4),
     &       SEDRSUM(4),SEDRAVG(4),SEDR(4),
     &       BR,BT,BN,F2,DEL,LAM,
     &       B(IARR,4),RMS(4),
     &       FIN(35)
      REAL*8 REALTIME
C
      CHARACTER RECTYPE*4,TELFMT*4,FLTID*4,TIMEFMT*4,DSN*50,MODE*4,
     &          COORD*4,CSYS*2,FILETYPE*6
      INTEGER*2 TIME(6),MODCNT(3),DATAID(2)
      INTEGER*4 IB192(25),IB96(5),IB48,SIB48
      LOGICAL*1 MAGSYS(32),REC(10000)
      REAL*4 B192(25,3),DEL192(25),LAM192(25),FMOD192(25),
     &       FMAG192(25),RMS192(3,25),
     &       B96(3,5),DEL96(5),LAM96(5),FMOD96(5),
     &       FMAG96(5),RMS96(3,5),
     &       B48(3),DEL48,LAM48,FMOD48,
     &       FMAG48,RMS48(3)
      REAL*4 SRMS48(3),SX48,SY48,SZ48,SFMOD48
      REAL*4 MTB(3,3),MTB5(3,3),MHG(3,3),POS(3),ANG(2),VEL(3)
      REAL*4 HDR1(100),HDR(32),DATA(341),SCFLD(155),POSN(13),ATTD(27),
     &       SUMOUT(568)
      REAL*8 DD,TD,TN,TP
C
      EQUIVALENCE (REC(1),SUMOUT(1)),
     &            (SUMOUT(1),HDR(1)),
     &            (HDR(1),RECTYPE),       (HDR(2),TELFMT),
     &            (HDR(3),FLTID),         (HDR(4),TIME(1)),
     &            (HDR(7),DD),            (HDR(9),TD),
     &            (HDR(11),TIMEFMT),      (HDR(12),TIMEPD),
     &            (HDR(13),MODCNT(1)),    (HDR(17),DATAID(1)),
     &            (SUMOUT(1),HDR1(1)),
     &            (HDR1(13),COORD),       (HDR1(77), MAGSYS(1)),
     &            (SUMOUT(33),DATA(1)),
     &            (DATA(1),FMAG48),       (DATA(2),FMOD48),
     &            (DATA(3),DEL48),        (DATA(4),LAM48),
     &            (DATA(5),B48(1)),       (DATA(8),RMS48(1)),
     &            (DATA(11),IB48),        (DATA(12),FMAG96(1)),
     &            (DATA(17),FMOD96(1)),   (DATA(22),DEL96(1)),
     &            (DATA(27),LAM96(1)),    (DATA(32),B96(1,1)),
     &            (DATA(47),RMS96(1,1)),  (DATA(62),IB96(1)),
     &            (DATA(67),FMAG192(1)),  (DATA(92),FMOD192(1)),
     &            (DATA(117),DEL192(1)),  (DATA(142),LAM192(1)),
     &            (DATA(167),B192(1,1)),  (DATA(242),RMS192(1,1)),
     &            (DATA(317),IB192(1)),
     &            (SUMOUT(342),SCFLD(1)),
     &            (SCFLD(148),SRMS48(1)), (SCFLD(151),SIB48),
     &            (SCFLD(152),SX48),      (SCFLD(153),SY48),
     &            (SCFLD(154),SZ48),      (SCFLD(155),SFMOD48),
     &            (SUMOUT(529),POSN(1)),
     &            (POSN(1),TN),           (POSN(3),TP),
     &            (POSN(5),POS(1)),       (POSN(8),VEL(1)),
     &            (POSN(11),RANGE),       (POSN(12),ANG(1)),
     &            (SUMOUT(542),ATTD(1)),
     &            (ATTD(1),MTB(1,1)),     (ATTD(10),MTB5(1,1)),
     &            (ATTD(19),MHG(1,1)),
     &            (IN(1),FIN(1))
C
C OPEN MAG SUMMARY DATA SET
C
      WRITE(6,*) 
      WRITE(6,*) 'ENTER INPUT SUMMARY DSN'
      READ(5,'(A)') DSN
      OPEN(70,FILE=DSN,STATUS='OLD',FORM='FORMATTED',
     &     RECORDTYPE='VARIABLE',RECL=8191,READONLY)
      WRITE(6,'(1X,A)') DSN
C
C DETERMINE DESIRED OUTPUT FORM (ASCII/BINARY)  SBK 12/07/95
C
      READ(5,'(A)') FILETYPE
C
C OPEN ASCII MAG HOUR AVERAGES DATA SET (PARKER PROCESSING)
C
      IF ( FILETYPE.EQ.'ASCII' ) THEN
       WRITE(6,*) 
       WRITE(6,*) 'ENTER OUTPUT ASCII HOUR AVERAGES DSN'
       READ(5,'(A)') DSN
       OPEN(11,FILE=DSN,STATUS='NEW',FORM='FORMATTED',
     &      CARRIAGECONTROL='LIST',RECORDTYPE='VARIABLE',RECL=8191)
       WRITE(6,'(1X,A)') DSN
      ELSE IF ( FILETYPE.EQ.'BINARY' ) THEN
C
C OPEN BINARY MAG HOUR AVERAGES DATA SET (MASTER HOUR AVERAGE DATABASE)
C
       WRITE(6,*) 
       WRITE(6,*) 'ENTER OUTPUT BINARY HOUR AVERAGES DSN'
       READ(5,'(A)') DSN
       OPEN(12,FILE=DSN,STATUS='NEW',FORM='UNFORMATTED',
     &      CONVERT='VAXG',RECORDTYPE='FIXED',RECL=35)
       WRITE(6,'(1X,A)') DSN
      ELSE
       WRITE(6,*)
       WRITE(6,*) 'INVALID OUTPUT FILE TYPE SPECIFIED: ',FILETYPE
       STOP
      END IF
C
C AVERAGING PERIOD  YEARS/DAYS/HOURS/MINS/SECS/MSECS
C
      DELTA(1) = 0
      DELTA(2) = 0
      DELTA(3) = 1  ! DELTA SET TO 1 HOUR
      DELTA(4) = 0
      DELTA(5) = 0
      DELTA(6) = 0
C
C INITIALIZE VARIABLES
C
      DO I = 1,4
       MAGCNT(I) = 0
       MAGSUM(I) = 0.0
       SEDRCNT(I) = 0
       SEDRSUM(I) = 0.0
      END DO
C
C READ SUMMARY RECORDS
C
      FIRST = .TRUE.
      STASH = .FALSE.
      IREAD = 0
      IWRITE = 0
      INREC = 0
      FILL = 999.0
      WRITE(6,*)
   10 CONTINUE
       READ(70,'(Q,<LEN>A1)',END=100,ERR=10) LEN,(REC(I),I=1,LEN)
       IREAD = IREAD + 1
C
C READ HDR1 DATA
C
       IF (DATAID(1).EQ.8) THEN
        CSYS = COORD(1:2)
        IF (MAGSYS(21).AND.MAGSYS(22)) THEN
         MODE = 'ERR '
        ELSE IF (MAGSYS(21)) THEN
         MODE = 'PRI '
        ELSE IF (MAGSYS(22)) THEN
         MODE = 'SEC '
        ELSE
         MODE = 'DUAL'
        END IF
        GOTO 10
       END IF
C
C READ LFM SUMMARY DATA
C
       IF (DATAID(1).EQ.1) THEN
        READ(FLTID,'(3X,I1)') ID
        MAG(1) = DATA(2)  ! FIELD MODULUS
        MAG(2) = DATA(5)  ! X COMPONENT FIELD
        MAG(3) = DATA(6)  ! Y COMPONENT FIELD
        MAG(4) = DATA(7)  ! Z COMPONENT FIELD
        SEDR(1) = POSN(11) ! S/C POSITION RANGE
        SEDR(2) = POSN(5)  ! X COMPONENT S/C POSITION
        SEDR(3) = POSN(6)  ! Y COMPONENT S/C POSITION
        SEDR(4) = POSN(7)  ! Z COMPONENT S/C POSITION
C
C SAVE START TIME FROM FIRST LFM RECORD - ROUND BACK TO START OF HOUR
C
        IF (FIRST) THEN
         TIME(4) = 0
         TIME(5) = 0
         TIME(6) = 0
C
C ESTABLISH STARTING AVERAGING WINDOW
C
         CALL MOVETIME(TIME,OLDTIME)  ! WINDOW START TIME
         CALL INCRTIME(TIME,DELTA)
         CALL MOVETIME(TIME,NEWTIME)  ! WINDOW END TIME
         FIRST = .FALSE.
        END IF
       ELSE
        GOTO 10
       END IF
C
C SKIP RECORDS WITH TIME TAG LESS THAN START TIME
C
       IF (REALTIME(TIME).LT.REALTIME(OLDTIME)) GOTO 10
C
   15  CONTINUE
C
C SUM VALUES IN AVERAGING PERIOD
C
       IF (REALTIME(TIME).GE.REALTIME(OLDTIME).AND.
     &     REALTIME(TIME).LT.REALTIME(NEWTIME)) THEN
C
        INREC = INREC + 1
C
C REJECT FIELD MAGNITUDE > 3 nT
C
        IF ( MAG(1).EQ.999.0 .OR. MAG(1).GT.3.0 ) GOTO 33
C
C ACCUMULATE FIELD MAGNITUDE AND COMPONENTS IF MODULUS IS NON-FILL
C
        DO I = 1,4
         MAGSUM(I) = MAGSUM(I) + MAG(I)
         MAGCNT(I) = MAGCNT(I) + 1
         B(MAGCNT(I),I) = MAG(I)
        END DO
C
   33   CONTINUE
C
C ACCUMULATE NAVIGATION DATA IF RANGE IS NON-FILL
C
        IF ( SEDR(1).LT.1.0 .OR. SEDR(1).GE.999.0 ) GOTO 66
C
        DO I = 1,4
         SEDRSUM(I) = SEDRSUM(I) + SEDR(I)
         SEDRCNT(I) = SEDRCNT(I) + 1
        END DO

   66   CONTINUE
C
C FLAG ACCUMULATORS AS LOADED
C
        STASH = .TRUE.
C
C COMPUTE AVERAGES WHEN INCREMENTING AVERAGING WINDOW
C
       ELSE IF (REALTIME(TIME).GE.REALTIME(NEWTIME)) THEN
C
        DO I = 1,4
C
         IF (MAGCNT(I).GE.1) THEN
          MAGAVG(I) = MAGSUM(I) / FLOAT(MAGCNT(I))
C
C COMPUTE FIELD MAGNITUDE RMS AND COMPONENT RMS'S
C
          DO J = 1,MAGCNT(I)
           RMS(I) = ABS(B(J,I)-MAGAVG(I))**2.0
          END DO
          RMS(I) = SQRT(RMS(I)/FLOAT(MAGCNT(I)))
C
         ELSE
          MAGAVG(I) = 999.0
          RMS(I) = 999.0
         END IF
C
        END DO
C
C COMPUTE FIELD MODULUS, ANGLES AND COMPONENT RMS NORM
C
        F1 = MAGAVG(1)
        BR = MAGAVG(2)
        BT = MAGAVG(3)
        BN = MAGAVG(4)
        IF (BR.NE.999.0.AND.
     &      BT.NE.999.0.AND.
     &      BN.NE.999.0) THEN
         CALL ANGLES(BR,BT,BN,F2,DEL,LAM)
         RMSC = SQRT(RMS(2)**2+RMS(3)**2+RMS(4)**2)
        ELSE
         F2 = 999.0
         DEL = 90.0
         LAM = 45.0
         RMSC = 999.0
        END IF
C
C COMPUTE RANGE AND POSITION COMPONENT AVERAGES
C
        DO I = 1,4
         IF (SEDRCNT(I).GE.1) THEN
          SEDRAVG(I) = SEDRSUM(I) / FLOAT(SEDRCNT(I))
         ELSE
          SEDRAVG(I) = 999.0
         END IF
        END DO
C
        IF ( MAGAVG(1).NE.999.0 ) THEN
C
C WRITE TO SYS$OUTPUT
C
         WRITE(6,'(1X,I1,1X,A2,1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,
     &         4(1X,F7.3),1X,I2,1X,F8.4,1X,I2)')
     &         ID,CSYS,OLDTIME,(MAGAVG(I),I=1,4),MAGCNT(1),
     &         SEDRAVG(1),SEDRCNT(1)
C
C WRITE HOUR AVERAGE ASCII FILE
C S/C ID, COORD SYS, TIME, F1, B1, B2, B3,
C RMSF1, RMSB1, RMSB2, RMSB3, NUMF1,
C RMAG, X, Y, Z, NUMR
C
         IF ( FILETYPE.EQ.'ASCII' ) 
     &    WRITE(11,890) ID,CSYS,(OLDTIME(I),I=1,3),
     &                  (MAGAVG(I),I=1,4),(RMS(I),I=1,4),MAGCNT(1),
     &                  (SEDRAVG(I),I=1,4),SEDRCNT(1)
C
C MASTER HOUR AVERAGE OUTPUT
C
         IN(1) = ID
         IN(2) = OLDTIME(1)
         IN(3) = OLDTIME(2)
         IN(4) = OLDTIME(3)
         FIN(5) = SEDRAVG(2)
         FIN(6) = SEDRAVG(3)
         FIN(7) = SEDRAVG(4)
         FIN(8) = SEDRAVG(1)
         FIN(9) = F1
         FIN(10) = F2
         FIN(11) = DEL
         FIN(12) = LAM
         FIN(22) = RMSC
         FIN(23) = RMS(1)
         FIN(24) = RMS(2)
         FIN(25) = RMS(3)
         FIN(26) = RMS(4)
         IF ( FILETYPE.EQ.'BINARY' ) 
     &    WRITE(12) (IN(I),I=1,4),(FIN(I),I=5,35)
C
         IWRITE = IWRITE + 1
C
        END IF
C
        DO I = 1,4
         MAGCNT(I) = 0
         MAGSUM(I) = 0.0
         MAGAVG(I) = 999.0
         RMS(I) = 999.0
         SEDRCNT(I) = 0
         SEDRSUM(I) = 0.0
         SEDRAVG(I) = 0.0
        END DO
        RMSC = 999.0
C
C INCREMENT AVERAGING WINDOW
C
        CALL MOVETIME(NEWTIME,OLDTIME)
        CALL INCRTIME(NEWTIME,DELTA)
        STASH = .FALSE.
C
       END IF
C
       IF (.NOT.STASH) GOTO 15
C
      GOTO 10
  100 CONTINUE
C
C COMPUTE AND OUTPUT FINAL AVERAGING PERIOD
C
      IF (STASH) THEN
C
       DO I = 1,4
C
        IF (MAGCNT(I).GE.1) THEN
         MAGAVG(I) = MAGSUM(I) / FLOAT(MAGCNT(I))
C
C COMPUTE FIELD MAGNITUDE RMS AND COMPONENT RMS'S
C
         DO J = 1,MAGCNT(I)
          RMS(I) = ABS(B(J,I)-MAGAVG(I))**2.0
         END DO
         RMS(I) = SQRT(RMS(I)/FLOAT(MAGCNT(I)))
C
        ELSE
         MAGAVG(I) = 999.0
         RMS(I) = 999.0
        END IF
C
       END DO
C
C COMPUTE FIELD MODULUS, ANGLES AND COMPONENT RMS NORM
C
       F1 = MAGAVG(1)
       BR = MAGAVG(2)
       BT = MAGAVG(3)
       BN = MAGAVG(4)
       IF (BR.NE.999.0.AND.
     &     BT.NE.999.0.AND.
     &     BN.NE.999.0) THEN
        CALL ANGLES(BR,BT,BN,F2,DEL,LAM)
        RMSC = SQRT(RMS(2)**2+RMS(3)**2+RMS(4)**2)
       ELSE
        F2 = 999.0
        DEL = 90.0
        LAM = 45.0
        RMSC = 999.0
       END IF
C
C COMPUTE RANGE AND POSITION COMPONENT AVERAGES
C
       DO I = 1,4
        IF (SEDRCNT(I).GE.1) THEN
         SEDRAVG(I) = SEDRSUM(I) / FLOAT(SEDRCNT(I))
        ELSE
         SEDRAVG(I) = 999.0
        END IF
       END DO
C
       IF ( MAGAVG(1).NE.999.0 ) THEN
C
C WRITE TO SYS$OUTPUT
C
        WRITE(6,'(1X,I1,1X,A2,1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,
     &        4(1X,F7.3),1X,I2,1X,F8.4,1X,I2)')
     &        ID,CSYS,OLDTIME,(MAGAVG(I),I=1,4),MAGCNT(1),
     &        SEDRAVG(1),SEDRCNT(1)
C
C WRITE HOUR AVERAGE ASCII FILE
C S/C ID, COORD SYS, TIME, F1, B1, B2, B3,
C RMSF1, RMSB1, RMSB2, RMSB3, NUMF1,
C RMAG, X, Y, Z, NUMR
C
        IF ( FILETYPE.EQ.'ASCII' ) 
     &   WRITE(11,890) ID,CSYS,(OLDTIME(I),I=1,3),
     &                 (MAGAVG(I),I=1,4),(RMS(I),I=1,4),MAGCNT(1),
     &                 (SEDRAVG(I),I=1,4),SEDRCNT(1)
C
C MASTER HOUR AVERAGE OUTPUT
C
        IN(1) = ID
        IN(2) = OLDTIME(1)
        IN(3) = OLDTIME(2)
        IN(4) = OLDTIME(3)
        FIN(5) = SEDRAVG(2)
        FIN(6) = SEDRAVG(3)
        FIN(7) = SEDRAVG(4)
        FIN(8) = SEDRAVG(1)
        FIN(9) = F1
        FIN(10) = F2
        FIN(11) = DEL
        FIN(12) = LAM
        FIN(22) = RMSC
        FIN(23) = RMS(1)
        FIN(24) = RMS(2)
        FIN(25) = RMS(3)
        FIN(26) = RMS(4)
        IF ( FILETYPE.EQ.'BINARY' ) 
     &   WRITE(12) (IN(I),I=1,4),(FIN(I),I=5,35)
C
        IWRITE = IWRITE + 1
C
       END IF
C
      END IF
C
      WRITE(6,*)
      WRITE(6,'(1X,''RECORDS READ:       '',I6.6)') IREAD
      WRITE(6,'(1X,''RECORDS PROCESSED:  '',I6.6)') INREC
      WRITE(6,'(1X,''RECORDS WRITTEN:    '',I6.6)') IWRITE
C
      STOP
  890 FORMAT(I1,1X,A2,1X,I2,1X,I3.3,1X,I2.2,
     &       8(1X,F7.3),1X,I2,4(1X,F8.4),1X,I2)
      END
      SUBROUTINE MOVETIME(INTIME,OUTTIME)
C
C ASSIGN INTIME ARRAY VALUES TO OUTTIME ARRAY
C
      INTEGER*2 INTIME(6),OUTTIME(6)
C
      DO I = 1,6
       OUTTIME(I) = INTIME(I)
      END DO
C
      RETURN
      END
      SUBROUTINE INCRTIME(TIME,DELTA)
C
C INCREMENT CALENDAR TIME (YY,DDD,HH,MM,SS,FFF) BY DELTA(6)
C
      INTEGER*2 TIME(6),DELTA(6)
C
      IYR = TIME(1)
      LEAP = 365
      IF (MOD(IYR,4).EQ.0.AND.IYR.NE.0) LEAP = 366
C
      TIME(6) = TIME(6) + DELTA(6)
      IF (TIME(6).GT.999) THEN
       TIME(5) = TIME(5) + 1
       TIME(6) = TIME(6) - 1000
      END IF
C
      TIME(5) = TIME(5) + DELTA(5)
      IF (TIME(5).GT.59) THEN
       TIME(4) = TIME (4) + 1
       TIME(5) = TIME (5) - 60
      END IF
C
      TIME(4) = TIME(4) + DELTA(4)
      IF (TIME(4).GT.59) THEN
       TIME(3) = TIME(3) + 1
       TIME(4) = TIME(4) - 60
      END IF
C
      TIME(3) = TIME(3) + DELTA(3)
      IF (TIME(3).GT.23) THEN
       TIME(2) = TIME(2) + 1
       TIME(3) = TIME(3) - 24
      END IF
C
      TIME(2) = TIME(2) + DELTA(2)
      IF (TIME(2).GT.LEAP) THEN
       TIME(1) = TIME(1) + 1
       TIME(2) = TIME(2) - LEAP
       IF (TIME(1).GT.99) TIME(1) = TIME(1) - 100
      END IF
C
      TIME(1) = TIME(1) + DELTA(1)
      IF (TIME(1).GT.99) TIME(1) = TIME(1) - 100
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION REALTIME(TIME)
C
C CONVERT CALENDAR TIME INTO DECIMAL YEAR REAL TIME.
C
      INTEGER*2 TIME(6)
      REAL*8 DAYS
C
      IYEAR = TIME(1)
      IDAY = TIME(2)
      IHOUR = TIME(3)
      IMIN = TIME(4)
      ISEC = TIME(5)
      MS = TIME(6)
C
      IF (IYEAR.GE.0.AND.IYEAR.LT.77) IYEAR = IYEAR + 2000
      IF (IYEAR.GE.77.AND.IYEAR.LE.99) IYEAR = IYEAR + 1900
C
      DAYS = 365.0D0
      IF (MOD(IYEAR,4).EQ.0.AND.IYEAR.NE.2000) DAYS = 366.0D0
      REALTIME = DBLE(IYEAR) +
     &           DBLE(IDAY-1)/DAYS +
     &           DBLE(IHOUR)/24.0D0/DAYS +
     &           DBLE(IMIN)/60.0D0/24.0D0/DAYS +
     &           DBLE(ISEC)/60.0D0/60.0D0/24.0D0/DAYS +
     &           DBLE(MS)/1000.0D0/60.0D0/60.0D0/24.0D0/DAYS
C
      RETURN
      END
      SUBROUTINE ANGLES(BX,BY,BZ,F2,DEL,LAM)
C
      REAL*4 LAM
C
C CALCULATE FIELD ANGLES AND MODULUS FROM COMPONENTS.
C BAD COMPONENTS ARE FILTERED BEFORE CALLING THIS ROUTINE.
C DELTA FILL VALUE IS 90.0
C LAMBDA FILL VALUE IS 45.0
C
      F2 = SQRT(BX**2+BY**2+BZ**2)
C
      IF (F2.NE.0.0) THEN
       DEL = ASIN(BZ/F2) * 57.2957795D0
      ELSE
       DEL = 90.0
      END IF
C
      IF (BX.NE.0) THEN
       LAM = 180.0 - ATAN2(BY,-BX) * 57.2957795D0
      ELSE IF (BX.EQ.0.0.AND.BY.NE.0.0) THEN
       LAM = 90.0
      ELSE
       LAM = 45.0
      END IF
C
      RETURN
      END