      PROGRAM EPHINTRP
      INTEGER*2 IFLAG,INUNIT,OUTUNIT,IOPEN1,IOPEN7,IOPEN9,INEPH,NPOS
      CHARACTER*30 INFILE,OUTFILE,EPHFILE
      CHARACTER*1 ANS
      CHARACTER*3 FPM,EPH
      DATA FPM,EPH/'FPM','EPH'/
      WRITE(*,1008)
1008  FORMAT(' FLOATING-POINT MAGNETIC FIELD DATA  FILE NAME?',/,
     1 ' (MUST HAVE UT,BX,BY,BZ IN FLOATING POINT IN THAT ORDER)',/,
     2 ' SUFFIX MUST BE "FPM"',/)
      READ(*,1004)INFILE
      INUNIT=1
      OUTUNIT=7
      INEPH=9
C
      CALL SUFIX(INFILE,FPM)
      CALL OPENOLD(INFILE,INUNIT,IOPEN1)
      IOPEN1=1
C
      PRINT*,' EPHEMERIS FILE NAME? (MUST HAVE .EPH EXT)'
      READ(*,1004)EPHFILE
      CALL SUFIX(EPHFILE,EPH)
      CALL OPENOLD(EPHFILE,INEPH,IOPEN9)
C
5     WRITE(*,1009)
1009  FORMAT(' OUTPUT FILE NAME?',/)
      READ(*,1004)OUTFILE
1004  FORMAT(A)
      OPEN(UNIT=OUTUNIT,FILE=OUTFILE,STATUS='NEW',FORM='FORMATTED')
      IOPEN7=1
      UT=-1.
      CALL FIRSTEPH(INEPH,UT,APA,IFLAG)
      IF(IFLAG.EQ.1)STOP
      WRITE(*,3002)UT
3002  FORMAT(' FIRST TIME IN EPH FILE IS ',F10.3,/)
      CALL FIRSTPA(INUNIT,UT,APA,IFLAG)
      IF(IFLAG.EQ.1)STOP
      WRITE(*,1002)UT
1002  FORMAT(' FIRST TIME IN MAG FILE IS ',F10.3,/)
      NPOS=24        !MODEL FIELD
      print*,' Start and Stop Times? (2F7.0)'
      READ(*,222)UT1,UTLST
      PRINT*,' TIME INTERVAL IN SECONDS?'
      READ(*,222)DELTEE
222   FORMAT(2F7.0)
      NTEE=INT((DELTEE+0.5)/1.024)
      DELTEE=1.024*FLOAT(NTEE)
      UT1=UT1-DELTEE
      DELTE2=DELTEE*0.5
      N8=INT((DELTEE+.06)/0.128)
      DO 135 JJ=1,370
      DO 134 JK=1,100
      UT1=UT1+DELTEE
      UTE=UT1
      UTB=UT1-DELTE2
      BMAG=0.
      BTMS=0.
      BAVE=0.
      DO 1005 IJ=1,N8
      CALL GETB(INUNIT,UTB,BMAG,IFLAG)
      IF (IFLAG.EQ.1)GOTO140
      BAVE=BAVE+BMAG
      BTMS=BTMS+1.
1005  UTB=UTB+0.128
      IF(BTMS.NE.0.)BAVE=BAVE/BTMS
      CALL GETEPH(INEPH,UTE,BMOD,IFLAG,NPOS)
      IF (IFLAG.EQ.1)GOTO140
      APA=BMAG-BMOD
      WRITE(7,21)UT1,BMOD,BMAG,APA
21    FORMAT(1X,F8.1,3(',',F10.2))
134   CONTINUE
      WRITE(*,19)UT1,BMOD,BMAG,APA
19    FORMAT(1X,F18.1,3(2X,F10.2))
      READ(*,1004)ANS
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GO TO 140
      IF(UTB.GT.UTLST)GO TO 140
135   CONTINUE
140   IF(IOPEN1.EQ.1)CLOSE(INUNIT)
      IF(IOPEN7.EQ.1)CLOSE(OUTUNIT)
      STOP
      END
C   
      SUBROUTINE GETB(INUNIT,UT,B,IFLAG)
      DIMENSION ARRAY(4,33),BCMP(3)
      INTEGER*2 INDX,IFLAG,NDX,INUNIT
1     NDX=INDX
      DO 100 J=NDX,33
      IF(UT.GT.ARRAY(1,J))GOTO100
      TDEL=ARRAY(1,J)-ARRAY(1,J-1)
      TPART1=(UT-ARRAY(1,J-1))/TDEL
      BCMP(1)=TPART1*(ARRAY(2,J)-ARRAY(2,J-1))+ARRAY(2,J-1)
      BCMP(2)=TPART1*(ARRAY(3,J)-ARRAY(3,J-1))+ARRAY(3,J-1)
      BCMP(3)=TPART1*(ARRAY(4,J)-ARRAY(4,J-1))+ARRAY(4,J-1)
      B=SQRT(DOT(BCMP,BCMP))
      RETURN
100   INDX=INDX+1
      DO 110 IJ=1,4
110   ARRAY(IJ,1)=ARRAY(IJ,33)
      CALL GETSECT(INUNIT,ARRAY(1,2),IFLAG)
      IF(IFLAG.EQ.1)THEN
         PRINT*,' MAGNETIC FIELD FILE EXHAUSTED'
         RETURN
      ENDIF
      INDX=2
      GO TO 1
      ENTRY FIRSTPA
      CALL GETSECT(INUNIT,ARRAY(1,2),IFLAG)
        IF(IFLAG.EQ.1)THEN
           PRINT*,' MAGNETIC FIELD FILE CAN NOT BE READ'
           RETURN
        ENDIF
      DO 120 IJ=1,4
120   ARRAY(IJ,1)=ARRAY(IJ,2)
      INDX=2
      UT=ARRAY(1,1)
      RETURN
      END


      SUBROUTINE GETSECT(INUNIT,BREC,IFLAG)
C     A subroutine to read sectors from disk.
C     IFLAG = 1 if reached end of file.
      LOGICAL*1 BREC(512)
      INTEGER*2 INUNIT,IFLAG
      IFLAG=0
      READ(UNIT=INUNIT,END=1000)BREC
      RETURN
1000  IFLAG=1
      PRINT *,' END OF DATA ON INPUT FILE'
      CLOSE (UNIT=INUNIT)
      RETURN
      END


      SUBROUTINE SUFIX(NAME,SUF)
      CHARACTER*1 NAME(51),SUF(3),SUF1
      DATA SUF1/'.'/
      DO 1 I=1,51
      IF(NAME(I).EQ.'.')GOTO4
1     IF(NAME(I).EQ.' ')GOTO3
      WRITE(*,2)NAME
2     FORMAT(' ***** ERROR ***** ERROR ***** ERROR *****',//,
     1 A,//,' exceeds the limit on file names (50 characters)')
      STOP
3     NAME(I)=SUF1
4     NAME(I+1)=SUF(1)
      NAME(I+2)=SUF(2)
      NAME(I+3)=SUF(3)
      RETURN
      END

      FUNCTION DOT(A,B)
      DIMENSION A(3),B(3)
      DOT=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
      RETURN
      END
C   
      SUBROUTINE GETEPH(INUNIT,UT,APA,IFLAG,NPOS)
      DIMENSION DATA(60,10)
      INTEGER*2 INDX,IFLAG,NDX,NPOS,INUNIT,JDEX,KFLAG
      KFLAG=0
      IFLAG=0
      UTI=UT
1     NDX=INDX
      DO 100 J=NDX,JDEX
      IF(UT.GT.DATA(2,J))GO TO 100
      IF(UT.LT.DATA(2,J-1))GO TO 131
      TDEL=DATA(2,J)-DATA(2,J-1)
      TPART1=(UT-DATA(2,J-1))/TDEL
      APA=TPART1*(DATA(NPOS,J)-DATA(NPOS,J-1))+DATA(NPOS,J-1)
        RETURN
100   INDX=INDX+1
      IF(JFLAG.EQ.1)GO TO 200
      IF((UT-DATA(2,10)).GT.6000.)GO TO 400
      DO 110 IJ=1,60
110   DATA(IJ,1)=DATA(IJ,10)
      JDEX=1
      DO 111 I=2,10
      CALL REEDEPH(INUNIT,UTI,DATA(1,I),JFLAG)
      IF(JFLAG.EQ.1)GO TO 410
      JDEX=JDEX+1
111   CONTINUE
      INDX=2
      GO TO 1
131   IF(KFLAG.EQ.1)GO TO 141
      KFLAG=1
      INDX=2
      GO TO 1
141   PRINT*,' UT=',UT,' LESS THAN ARRAY TIME',DATA(2,1)
      IFLAG=1
      RETURN
400   JDEX=0
      INDX=2
      DO 401 I=1,10
      CALL REEDEPH(INUNIT,UTI,DATA(1,I),JFLAG)
      IF(JFLAG.EQ.1)GO TO 410
      JDEX=JDEX+1
401   CONTINUE
      GO TO 1
410   PRINT*,' EOF ON EPHEMERIS FILE'
      IF(JDEX.GT.1)GO TO 1
      RETURN
C
      ENTRY FIRSTEPH
      JDEX=0
      DO 112 I=1,10
      CALL REEDEPH(INUNIT,UT,DATA(1,I),JFLAG)
      IF(JFLAG.EQ.1)GO TO 200
      JDEX=JDEX+1
112   CONTINUE
      UT=DATA(2,1)
      INDX=2
      RETURN
200   PRINT*,' EPHEMERIS FILE CAN NOT BE READ'
      IFLAG=1
      RETURN
      END

      SUBROUTINE REEDEPH(INUNIT,T1,DATA,IFLAG)
      INTEGER*4 JATA(60),ICON,ITMP
      INTEGER*2 IFLAG,INDFYL,IEOR,INUNIT
      DIMENSION DATA(60),FACTOR(60)
      DATA INDFYL/0/
      DATA FACTOR/1.,.001,3*.0001,3*1.E-7,2*.0001,2*1.0E-6,
     1 3*1.E-7,2*1.E-6,1.E-7,1.E-6,2*1.E-7,1.E-6,1.E-7,4*.0001,
     2 1.E-7,6*1.E-6,1.E-7,.0001,2*1.E-6,.0001,2*1.E-6,.0001,6*1.,
     3 1.E-6,8*.0001,3*1./
C     1 JULIAN DATE (DAYS)
C     2 UT, Milliseconds
C     3 X, ECI, km
C     4 Y, ECI, km
C     5 Z, ECI, km
C     6 VX, ECI, km/SEC
C     7 VY, ECI, km/SEC
C     8 VZ, ECI, km/SEC
C     9 RADIUS, EARTH CENTER TO SATELLITE, km
C    10 Altitude, km
C    11 Latitude, Deg
C    12 Longitude, Deg
C    13 Velocity, km/sec
C    14 Local Time, hr
C    15 Radius, Mag, EMR (6371.2 km)
C    16 Latitude. Mag, Deg
C    17 Longitude, Mag, Deg
C    18 Radius, SM, EMR
C    19 Latitude, SM, Deg
C    20 Local Time, SM, hr
C    21 Radius, GSM, EMR
C    22 Latitude, GSM, Deg
C    23 Local Time, GSM, hr
C    24 B, nT
C    25 BX, ECI, nT
C    26 BY, ECI, nT
C    27 BZ, ECI, nT
C    28 Mag Local Time, hr
C    29 Solar Zenith Angle, Deg
C    30 Invariant Latitude, Deg
C    31 B100N latitude, Deg
C    32 B100N longitude, Deg
C    33 B100S latitude, Deg
C    34 B100S longitude, Deg
C    35 L-shell, EMR
C    36 Bmin, nT
C    37 Bmin Latitude, Deg
C    38 Bmin Longitude, Deg
C    39 Bmin Altitude, km
C    40 Bconj Latitude, Deg
C    41 Bconj Longitude, Deg
C    42 Bconj Altitude, km
C    43 X Sun Position, ECI, km
C    44 Y Sun Position, ECI, km
C    45 Z Sun Position, ECI, km
C    46 X Moon Position, ECI, km
C    47 Y Moon Position, ECI, km
C    48 Z Moon Position, ECI, km
C    49 Right Ascension of Greenwich
C    50 B100N, nT
C    51 B100S, nT
C    52 Mx dipole moment, ECI, nT
C    53 My dipole moment, ECI, nT
C    54 Mz dipole moment, ECI, nT
C    55 Dx dipole offset, ECI, nT
C    56 Dy dipole offset, ECI, nT
C    57 Dz dipole offset, ECI, nT
C    58 Vacant
C    59 Vacant
C    60 Vacant
C    

      IF(INDFYL.EQ.1.AND.IFLAG.EQ.1)RETURN
      ICON=2**30
      IFLAG=0
1     READ(UNIT=INUNIT,END=1000)JATA
      READ(UNIT=INUNIT,END=1000)IEOR
      CALL SWAP4(JATA(2),ITMP)
      TIME=.001*FLOAT(ITMP-ICON)
      IF(TIME.LT.(T1-300.))GO TO 1
      DO 2 I=1,57
      CALL SWAP4(JATA(I),ITMP)
2     DATA(I)=FLOAT(ITMP-ICON)*FACTOR(I)
      T1=TIME
      RETURN
1000  WRITE(*,1001)
1001  FORMAT(' HIT THE END OF THE EPHEMERIS FILE',/,
     1 ' PROCESSING OF EPHEMERIS STOPPED')
      IFLAG=1
      INDFYL=1
      RETURN
      END

      SUBROUTINE OPENOLD(FILNAM,NUNIT,NOPEN)
      CHARACTER*30 FILNAM
      CHARACTER*1 ANS
      INTEGER*2 NUNIT,NOPEN
      INTEGER*4 LEN,NLEN
      LOGICAL OPEND,EXST
1     INQUIRE(FILE=FILNAM,EXIST=EXST,OPENED=OPEND)
      IF(EXST)GOTO5
      NOPEN=0
      WRITE(*,2)CHARNB(FILNAM)
2     FORMAT(' FILE ',A,' DOES NOT EXIST',/,' DO YOU NEED IT? (Y,N)')
      CALL PAUZ(ANS)
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')RETURN
      WRITE(*,3)
3     FORMAT(' ENTER FILE NAME, INCLUDING FULL PATH',/)
      READ(*,4)FILNAM
4     FORMAT(A)
      GO TO 1
5     IF(OPEND)GO TO 6
      OPEN (UNIT=NUNIT,FILE=FILNAM,STATUS='OLD',ACCESS='TRANSPARENT',
     1 FORM='UNFORMATTED')
      NOPEN=1
C     GET THE FILE LENGTH IN BYTES
6     INQUIRE(FILE=FILNAM,FLEN=LEN)
      NLEN=LEN/512
      WRITE(*,7)CHARNB(FILNAM),LEN,NLEN
7     FORMAT(' FILE ',A,' IS ',I8,' BYTES (',I5,'+ SECTORS) LONG')
      RETURN
      END
      SUBROUTINE SWAP4(IN,IOUT)
      LOGICAL*1 IN(4),IOUT(4)
      IOUT(1)=IN(4)
      IOUT(2)=IN(3)
      IOUT(3)=IN(2)
      IOUT(4)=IN(1)
      RETURN
      ENTRY SWAP2
      IOUT(1)=IN(2)
      IOUT(2)=IN(1)
      RETURN
      END


      SUBROUTINE PAUZ(ANS)
      CHARACTER*1 ANS
      READ(*,1)ANS
1     FORMAT(A)
      RETURN
      END


