PROGRAM VIM5READ
C
C READ VIM-5 DETAIL VOYAGER MAG DATA 
C S. KRAMER  3/29/94
C S. KRAMER  5/24/2006
C
      PARAMETER (IARR=100000)
C
C DATA INPUT PARAMETERS
C
      CHARACTER RECTYPE*4,TELFMT*4,FLTID*4,TIMEFMT*4,RUNMONTH*4,
     &          DSN*50,MODE*4,COORD*2
      INTEGER*2 DATE(6),DETAIL1(6),DETAIL2(6),TBEG(6),TEND(6),
     &          MODCNT(3),DATAID(2)
      INTEGER*4 RUNDAY,RUNYEAR
      LOGICAL*1 MAGSYS(32),SYS2(32),REC(10000)
      REAL*4 DETOUT(2432),DATA(2400),HDR(32),HDR1(100),
     &       PRIDAT(400,3),SECDAT(8,3),SCFDAT(400,3),AMBDAT(400,3),
     &       DATA1(IARR,3),DATA2(IARR,3),
     &       B(IARR),BX(IARR),BY(IARR),BZ(IARR),BDEL(IARR),BLAM(IARR)
      REAL*8 DD,TD,TN,TP,DECDY,DECHR,TDAY(IARR),THOUR(IARR),TIME(IARR)
C
      EQUIVALENCE (HDR(1),RECTYPE),    (HDR(2),TELFMT),
     &            (HDR(3),FLTID),      (HDR(4),DATE(1)),
     &            (HDR(7),DD),         (HDR(9),TD),
     &            (HDR(11),TIMEFMT),   (HDR(12),TIMEPD),
     &            (HDR(13),MODCNT(1)), (HDR(17),DATAID(1))
C
      EQUIVALENCE (HDR1(7),RUNYEAR),   (HDR1(9),RUNMONTH),
     &            (HDR1(10),RUNDAY),   (HDR1(77), SYS2(1))
C
      EQUIVALENCE (DETOUT(1),HDR(1)),     (DETOUT(33),DATA(1)),
     &            (PRIDAT(1,1),DATA(1)),  (SECDAT(1,1),DATA(1201)),
     &            (AMBDAT(1,1),DATA(1)),  (SCFDAT(1,1),DATA(1201)),
     &            (DETOUT(1),HDR1(1))
C
      EQUIVALENCE (DETOUT(1),REC(1))
C
      DATA TITLE/' '/, SUBTITLE/' '/
C
C VIM5 PRIMARY DATA RESOLUTION
C
      DETAIL1(1) = 0 
      DETAIL1(2) = 0
      DETAIL1(3) = 0
      DETAIL1(4) = 0
      DETAIL1(5) = 0
      DETAIL1(6) = 480
C
C VIM5 SECONDARY DATA RESOLUTION
C
      DETAIL2(1) = 0 
      DETAIL2(2) = 0
      DETAIL2(3) = 0
      DETAIL2(4) = 0
      DETAIL2(5) = 24
      DETAIL2(6) = 0
C
      WRITE(6,*) 'ENTER DETAIL DSN'
      READ(5,'(A)') DSN
      OPEN(70,FILE=DSN,STATUS='OLD',FORM='FORMATTED',
     &     RECORDTYPE='VARIABLE',RECL=8191,READONLY)
      WRITE(6,'(1X,A)') DSN
C
      WRITE(6,*)
      WRITE(6,*) 'ENTER START TIME YY/DDD/HH:MM:SS.SSS'
      READ(5,'(I2,1X,I3,3(1X,I2),1X,I3)') TBEG
      WRITE(6,'(1X,I2,1X,I3,3(1X,I2),1X,I3)') TBEG
      WRITE(6,*)
      WRITE(6,*) 'ENTER STOP TIME YY/DDD/HH:MM:SS.SSS'
      READ(5,'(I2,1X,I3,3(1X,I2),1X,I3)') TEND
      WRITE(6,'(1X,I2,1X,I3,3(1X,I2),1X,I3)') TEND
C
      IDAYS = TEND(2) - TBEG(2)
      XL = DECHR(TBEG)
      XH = DECHR(TEND) + 24.0*FLOAT(IDAYS)
      WRITE(SUBTITLE,'(I4,'', DAY '',I3.3,A2)') 
     & TBEG(1)+1900,TBEG(2),'\E'
C
      FILL = 999.0
      NCNT = 0
      NPT1 = 0
      NPT2 = 0
   10 CONTINUE
       READ(70,'(Q,<LEN>A1)',END=100,ERR=10) LEN,(REC(I),I=1,LEN)
       NCNT = NCNT + 1
C
C EXTRACT HDR1 INFO
C
       IF ( DATAID(1).EQ.8 .AND. NCNT.EQ.1 ) THEN
C
        DO II = 1,32
         MAGSYS(II) = SYS2(II)
C        WRITE(6,'(1X,''MAGSYS('',I2,'') = '',L1)') II,MAGSYS(II)
        END DO
C
        IF ( MAGSYS(4).OR.MAGSYS(5) ) THEN
         CONTINUE
        ELSE
         WRITE(6,*)
         WRITE(6,*) 'NOT A DETAIL DATA SET'
        END IF
C
        IF ( MAGSYS(7) ) THEN
         IF ( MAGSYS(17) ) THEN
          COORD = 'S3'
         ELSE IF ( MAGSYS(18) ) THEN
          COORD = 'L1'
         ELSE IF ( MAGSYS(19) ) THEN
          COORD = 'U1'
         ELSE IF ( MAGSYS(20) ) THEN
          COORD = 'N1'
         ELSE
          COORD = 'HG'
         END IF
        ELSE
         COORD = 'PL'
        END IF
C
        WRITE(6,*)
        WRITE(6,801) NCNT,RECTYPE,DATE,FLTID,COORD
        READ(FLTID,'(3X,I1)') ID
        WRITE(6,*) '  1 - PRIMARY'
        WRITE(6,*) '  2 - SECONDARY'
        READ(5,*) IDAT
        IF (IDAT.EQ.1) THEN MODE = ' PRI'
        IF (IDAT.EQ.2) THEN MODE = ' SEC'
C
       END IF
C
C ASSUME TIME REQUEST << 1 YR
C
       IF ( DATE(1).LT.TBEG(1) ) GOTO 10
       IF ( DATE(2).LT.TBEG(2) ) GOTO 10
       IF ( DATE(2).EQ.TBEG(2) .AND. 
     &      DECHR(DATE).LT.DECHR(TBEG) ) GOTO 10
       IF ( DATE(2).GE.TEND(2).AND.
     &      DECHR(DATE).GE.DECHR(TEND) ) GOTO 100       
C
       IF ( DATAID(1).NE.1 ) GOTO 10
C
C PROCESS VIM-5 LOW FIELD MAG DATA
C       
       DO I = 1,400
C
        ISEC = (I-1)/50 + 1
C
C GET DETAIL PRI/SEC DATA
C
         WRITE(6,800) 
     &    ID,DATE,I,(PRIDAT(I,J),J=1,3),ISEC,(SECDAT(ISEC,J),J=1,3)
         NPT1 = NPT1 + 1
         DATA1(NPT1,1) = PRIDAT(I,1)
         DATA1(NPT1,2) = PRIDAT(I,2)
         DATA1(NPT1,3) = PRIDAT(I,3)
         IF ( MOD(I-1,50).EQ.0 ) THEN
          NPT2 = NPT2 + 1
          IF (IDAT.EQ.2) THEN
           TDAY(NPT2) = DECDY(DATE)
           THOUR(NPT2) = DECHR(DATE) + 24.0*(DATE(2)-TBEG(2))
C          WRITE(6,*) NPT2,DATE,THOUR(NPT2)
           CALL INC_TIME(DATE,DETAIL2)
          END IF 
          DATA2(NPT2,1) = SECDAT(ISEC,1)
          DATA2(NPT2,2) = SECDAT(ISEC,2)
          DATA2(NPT2,3) = SECDAT(ISEC,3)
         END IF
C
       END DO
C
      GOTO 10
  100 CONTINUE
C
      WRITE(6,*)
      WRITE(6,*) NCNT, ' RECORDS READ'
      WRITE(6,*) NPT1, ' PRI/AMB/SCF POINTS'
      WRITE(6,*) NPT2, ' SEC POINTS'
C
      IF ( NCNT.LE.0.OR.NCNT.GT.IARR ) GOTO 999
      IF ( NPT1.LE.0.OR.NPT1.GT.IARR ) GOTO 999
      IF ( IDAT.EQ.2 .AND. (NPT2.LE.0.OR.NPT2.GT.IARR) ) GOTO 999
C
C PRI/AMB/SCF PLOT DATA
C
      IF ( IDAT.NE.2 ) THEN
       NMAX = NPT1
       DO I = 1,NMAX
C
C ASSIGN DECIMAL HOUR TO TIME ARRAY
C
        TIME(I) = THOUR(I)
C
C ASSIGN DECIMAL DAY TO TIME ARRAY
C
C       TIME(I) = TDAY(I)
C       WRITE(6,805) 
C    &   ID,TIME(I),I,(DATA1(I,J),J=1,3)
        BX(I) = DATA1(I,1)
        BY(I) = DATA1(I,2)
        BZ(I) = DATA1(I,3)
        IF ( BX(I).NE.FILL.AND.
     &       BY(I).NE.FILL.AND.
     &       BZ(I).NE.FILL) THEN
C
C COMPUTE FIELD MODULUS
C
         B2 = SQRT(BX(I)**2+BY(I)**2+BZ(I)**2)
         B(I) = B2
C
C COMPUTE LATITUDINAL ANGLE
C
         IF (B2.NE.0.0) THEN
          BDEL(I) = ASIN(BZ(I)/B2) * 180.0/3.14159
         ELSE
          BDEL(I) = FILL
         END IF
C
C COMPUTE LONGITUDINAL ANGLE
C
         IF (BX(I).NE.0.0.OR.BY(I).NE.0.0) THEN
          BLAM(I) = 180.0 - ATAN2(BY(I),-BX(I))  * 180.0/3.14159
         END IF
C
        ELSE
C
         B(I) = FILL
         BDEL(I) = FILL
         BLAM(I) = FILL
C
        END IF
       END DO
C
C SEC PLOT DATA
C
      ELSE IF ( IDAT.EQ.2 ) THEN
       NMAX = NPT2
       DO I = 1,NMAX
        TIME(I) = THOUR(I)
C       TIME(I) = TDAY(I)
C       WRITE(6,805) 
C    &   ID,TIME(I),I,(DATA2(I,J),J=1,3)       
        BX(I) = DATA2(I,1)
        BY(I) = DATA2(I,2)
        BZ(I) = DATA2(I,3)
        IF ( BX(I).NE.FILL.AND.
     &       BY(I).NE.FILL.AND.
     &       BZ(I).NE.FILL) THEN
C
C COMPUTE FIELD MODULUS
C
         B2 = SQRT(BX(I)**2+BY(I)**2+BZ(I)**2)
         B(I) = B2
C
C COMPUTE LATITUDINAL ANGLE
C
         IF (B2.NE.0.0) THEN
          BDEL(I) = ASIN(BZ(I)/B2) * 180.0/3.14159
         ELSE
          BDEL(I) = FILL
         END IF
C
C COMPUTE LONGITUDINAL ANGLE
C
         IF (BX(I).NE.0.0.OR.BY(I).NE.0.0) THEN
          BLAM(I) = 180.0 - ATAN2(BY(I),-BX(I))  * 180.0/3.14159
         END IF
C
        ELSE
C
         B(I) = FILL
         BDEL(I) = FILL
         BLAM(I) = FILL
C
        END IF
       END DO
      END IF
C
      WRITE(6,*)
      WRITE(6,'(1X,''NUMBER OF RECORDS READ: '',I5)') NCNT
      WRITE(6,'(1X,''NUMBER OF POINTS EXTRACTED: '',I5)') NMAX
      CLOSE(70)
C
  800 FORMAT(1X,I1,1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,
     &       2(1X,I6,3(1X,F7.3)))
  801 FORMAT(1X,I5,1X,A4,1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,1X,A4,1X,A2)
  802 FORMAT(1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,2(1X,I4,3(1X,F7.3)))
  805 FORMAT(1X,I1,1X,F12.8,1X,I6,3(1X,F7.3))
  999 CONTINUE
C
      STOP
      END
      REAL*8 FUNCTION DECDY(DATE)
C
C CONVERT HIGH RESOLUTION CALENDAR DATA TO DECIMAL DAY
C
      INTEGER*2 DATE(6)
      REAL*8 DAY,HOUR,MIN,SEC,MSEC
      MSEC = DFLOAT(DATE(6))
      SEC = DFLOAT(DATE(5)) + MSEC/1000.0D0
      MIN = DFLOAT(DATE(4)) + SEC/60.0D0
      HOUR = DFLOAT(DATE(3)) + MIN/60.0D0
      DAY = DFLOAT(DATE(2)) + HOUR/24.0D0
      DECDY = DAY
      RETURN
      END
      REAL*8 FUNCTION DECHR(DATE)
C
C CONVERT HIGH RESOLUTION CALENDAR DATA TO DECIMAL HOUR
C
      INTEGER*2 DATE(6)
      REAL*8 DAY,HOUR,MIN,SEC,MSEC
      MSEC = DFLOAT(DATE(6))
      SEC = DFLOAT(DATE(5)) + MSEC/1000.0D0
      MIN = DFLOAT(DATE(4)) + SEC/60.0D0
      HOUR = DFLOAT(DATE(3)) + MIN/60.0D0
      DAY = DFLOAT(DATE(2))
      DECHR = HOUR 
      RETURN
      END
      SUBROUTINE INC_TIME(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) 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)
C
      RETURN
      END
      REAL*8 FUNCTION DECYR(DATE)
C
C CONVERT HIGH RESOLUTION CALENDAR TIME INTO DECIMAL TIME
C
      INTEGER*2 DATE(6)
      REAL*8 YEAR,DAY,HOUR,MIN,SEC,MSEC,LEAP
      LEAP = 365.0D0
      IF (MOD(DATE(1),4).EQ.0) LEAP = 366.0D0
      MSEC = DFLOAT(DATE(6))
      SEC = DFLOAT(DATE(5)) + MSEC/1000.0D0
      MIN = DFLOAT(DATE(4)) + SEC/60.0D0
      HOUR = DFLOAT(DATE(3)) + MIN/60.0D0
      DAY = DFLOAT(DATE(2)-1) + HOUR/24.0D0
      YEAR = DFLOAT(DATE(1)) + DAY/LEAP
      DECYR = YEAR
      RETURN
      END
      REAL*4 FUNCTION DVAL(VAL)
C
C RETURN REAL*4 DECIMAL DAY FROM CALENDAR INPUT (YYDDDHH)
C
      CHARACTER*10 VALUE
      IVAL = INT(VAL)
      WRITE(VALUE,'(I7)') IVAL
      READ(VALUE,'(I2,I3.3,I2.2)') IY,ID,IH
      DVAL = FLOAT(ID) + FLOAT(IH)/24.0
      RETURN
      END
      SUBROUTINE GETDATE(SYSDATE)
C
C CALL DEC FORTRAN SUBROUTINE TO GET SYSTEM DATE 
C
      CHARACTER SYSDATE*9
C
      CALL DATE(SYSDATE)
C
      RETURN 
      END
      SUBROUTINE GETTIME(SYSTIME)
C
C CALL DEC FORTRAN SUBROUTINE TO GET SYSTEM TIME
C
      CHARACTER SYSTIME*8
C
      CALL TIME(SYSTIME)
C
      RETURN
      END