C************************************************************************
      PROGRAM TRAVCALL
C-----------------------------------------------------------------------
C     PROGRAMMER J.D. KOCH 95-10-02
C-----------------------------------------------------------------------
C     MODIFICATIONS
C     96-04-23	S.A.MCLAUGHLIN	ADDED CODE TO ALLOW GIOCALL TO RUN
C		UNDER DOS VERSION 6.2 USING COMPILER:  MS-FORTRAN
C		VERSION 5.0 FOR DOS.  TO RUN UNDER DOS, FIND ALL
C		OCCURRENCES OF THE CHARACTERS 'DOS CHANGES' AND MAKE
C		THE SPECIFIED CODE CHANGES.
C-----------------------------------------------------------------------
************        PROGRAM PURPOSE                         ************
*     THIS PROGRAM IS INTENDED MERELY TO TEST TRAVEL.F AND GIODAT.F    *
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER NAME
      DIMENSION S(3),P(3),G(3)
      OPEN(UNIT=7,FILE='POSITION',ACCESS='SEQUENTIAL',FORM='FORMATTED')
      TUTC = 1000
      PRINT *,'ENTER THE STATION NAME'
      READ *, NAME
      PRINT *,'ENTER THE GEOCENTRIC STATION COORDINATES, X,Y,Z (IN KM)'
      READ *, S(1),S(2),S(3)
      PRINT *,'ENTER SIGNAL RECEPTION TIME (IN EME50)'
      READ *, TUTC
      PRINT *,' '
      CALL TRAVEL(TUTC,S,DTG,P,G,IER)
      IF (IER .EQ. 0 ) THEN 
         PRINT *,'AT TUTC =', TUTC,' :'
	 PRINT *,' SIGNAL TIME FROM GIOTTO TO STATION:', DTG,' SEC'
* DOS CHANGES:	 COMMENT THE FOLLOWING LINE
	 PRINT *,''
* DOS CHANGES:  UNCOMMENT THE FOLLOWING LINE
*	 PRINT *,' '
         PRINT *,' POSITION OF STATION AT TUTC IN HELIOCENTRIC EME50'
         PRINT *,' X',P(1),' KM  Y',P(2),' KM  Z',P(3),' KM'
         PRINT *,' POSITION OF GIOTTO AT TUTC IN HELIOCENTRIC EME50'
         PRINT *,' X',G(1),' KM  Y',G(2),' KM  Z',G(3),' KM'
      ELSE IF (IER .EQ. 1) THEN
         PRINT *,'TUTC TOO EARLY, NO DATA IN POSITION DATA FILE'   
      ELSE IF (IER .EQ. 2) THEN
         PRINT *,'TUTC TOO LATE , NO DATA IN POSITION DATA FILE'   
      ELSE
         PRINT *,'NO CONVERGENCE IN SOLUTION OF LIGHT TIME EQUATION'
      ENDIF
      CLOSE(UNIT=7)
      STOP
      END

C-----------------------------------------------------------------------      
      SUBROUTINE GIODAT(TUTC,PVG,PVH,PVE,PVS,IERR)                            
C-------------------------------------------------                            
C S.J ISSUE 04 - 14.01.85                                                     
C------------------------                                                     
C                                                                             
************        SUBROUTINE IDENTIFICATION               ************      
*                                                                      *      
*   GIODAT IS AN INTERFACE ROUTINE TO THE EXPERIMENTER SUPPORT DATA    *      
*   FILE ON LOGICAL UNIT 7 WHICH CONTAINS POSITION DATA FOR THE        *      
*   GIOTTO PROBE, HALLEYS COMET, THE EARTH AND THE SUN GIVEN IN THE    *      
*   BARYCENTRIC REFERENCE FRAME OF 1950 ( UNITS = KM )                 *      
*   THE APPROXIMATE VELOCITIES FOR THE FOUR BODIES WILL BE DERIVED     *      
*   FROM THE POSITION DATA AND GIVEN IN KM/S.                          *      
*                                                                      *      
************************************************************************      
*                                                                      *      
************        VARIABLE IDENTIFICATIONS               *************      
*                                                                      *      
*****            INPUT VARIABLES ( ARGUMENT LIST )                   ***      
*                                                                      *      
*   TUTC      -  UNIVERSAL COORDINATED TIME                            *      
*                                                                      *      
*****            OUTPUT VARIABLES ( ARGUMENT LIST )                  ***      
*                                                                      *      
*   PVG(1-3)  -  POSITION OF THE GIOTTO PROBE AT TIME TUTC             *      
*   PVG(4-6)  -  VELOCITY OF THE GIOTTO PROBE AT TIME TUTC             *      
*   PVH(1-3)                                                           *      
*   PVH(4-6)  -  DITTO FOR HALLEYS COMET                               *      
*   PVE(1-3)                                                           *      
*   PVE(4-6)  -  DITTO FOR THE EARTH                                   *      
*   PVS(1-3)                                                           *      
*   PVS(4-6)  -  DITTO FOR THE SUN                                     *      
*   IERR      -  ERROR CODE                                            *      
*                =0 : NO ERROR                                         *      
*                =1 : TUTC TOO EARLY                                   *      
*                =2 : TUTC TOO LATE                                    *      
*                                                                      *      
*****            LOCAL VARIABLES                                     ***      
*                                                                      *      
*   ET        -  MODIFIED JULIAN DAY ( 1950.0 )                        *      
*   PGI(3,3)  -  TO STORE POSITION DATA FOR THREE SAMPLES FOR THE      *      
*                GIOTTO PROBE. USED IN THE QUADRATIC INTERPOLATION     *      
*   PHI(3,3)  -  DITTO FOR HALLEYS COMET                               *      
*   PEI(3,3)  -  DITTO FOR THE EARTH                                   *      
*   PSI(3,3)  -  DITTO FOR THE SUN                                     *      
*   T(3)      -  TO STORE REFERENCE TIME T0 FOR THE INTERPOLATION.     *      
*   SKIP(3)   -  USED FOR SKIPPING RECORDS ON THE DATA FILE            *      
*                                                                      *      
*****            CONSTANT IDENTIFICATION                             ***      
*                                                                      *      
*   HSEC      -  (STEP SIZE)*86400                                     *      
*   LUN       -  LOGICAL UNIT NUMBER OF FILE/TAPE CONTAINING           *      
*                THE POSITION DATA.                                    *      
*   STPUTC    -  STEP IN UTC AT TIME SUTC                              *      
*                                                                      *      
*****            TEMPORARY VARIABLES                                 ***      
*                                                                      *      
*   (1)                                                                *      
*   TMP(1)    -  0.5*(P-1)*(P-2)                                       *      
*   TMP(2)    -  P*(P-2)                                               *      
*   TMP(3)    -  0.5*P*(P-1)                                           *      
*   (2)                                                                *      
*   TMP(1)    -  (P-1.5)/HSEC                                          *      
*   TMP(2)    -  (2*P-2)/HSEC                                          *      
*   TMP(3)    -  (P-0.5)/HSEC                                          *      
*   WHERE ET(P) IS THE EPOCH FOR WHICH THE POSITIONS ARE REQUIRED.     *      
*   ( RANGE OF P :  0:2 )                                              *      
*                                                                      *      
************************************************************************      
*                                                                      *      
************  TYPE DECLARATION AND STORAGE ALLOCATION       ************      
*                                                                      *      
      IMPLICIT REAL*8(A-H,O-Z)                                                
      LOGICAL FIRST                                                           
      CHARACTER*1 SKIP(3)                                                     
      DIMENSION PVG(6),PVH(6),PVE(6),PVS(6),PGI(3,3),                         
     &          PHI(3,3),PEI(3,3),PSI(3,3),T(3),TMP(3)                        
      DATA FIRST/.TRUE./                                                      
      DATA LUN/7/                                                             
*                                                                      *      
************************************************************************      
*                                                                      *      
************  INITIALISE AND VERIFY INPUT DATA              ***BLOCK0000      
*                                                                      *      
      IERR=0                                                                  
C    IF FIRST CALL READ CONTROL RECORD AND SET CONSTANTS                      
      IF(FIRST)THEN                                                           
      READ(LUN,10) IYGEN,MOGEN,DGEN,NS,NR,ET1,ET2,SUTC,ETUTC1,ETUTC2          
   10 FORMAT(2I5,D16.10,2I5,/,5D16.10)                                        
C    DETERMINE STEP IN UTC AT TIME SUTC                                       
      STPUTC=ETUTC2-ETUTC1                                                    
C    DETERMINE STEP SIZE                                                      
      H=(ET2-ET1)/(DFLOAT(NR)-1.D0)                                           
      HSEC=H*86400.D0                                                         
      FIRST=.FALSE.                                                           
      ENDIF                                                                   
C    DETERMINE ET ( MODIFIED JULIAN DAY )                                     
      ET=TUTC+ETUTC1                                                          
      IF(ET.GE.SUTC)ET=ET+STPUTC                                              
C    VERIFY EPOCH ON FILE                                                     
      IF(ET.LT.ET1-1.D-4)GOTO 910                                             
      IF(ET.GT.ET2+1.D-4)GOTO 920                                             
      DO 20 I=1,6                                                             
      PVG(I)=0.D0                                                             
      PVH(I)=0.D0                                                             
      PVE(I)=0.D0                                                             
      PVS(I)=0.D0                                                             
   20 CONTINUE                                                                
*                                                                      *      
************************************************************************      
*                                                                      *      
*******    DETERMINE THE POSITION VECTORS FOR EPOCH ET       **BLOCK0100      
*                                                                      *      
C    FIND RECORD NUMBER                                                       
      IREC=DINT((ET-ET1)/H)+1                                                 
  110 CONTINUE                                                                
C    POSITION THE T0 SAMPLE TO ENSURE T1 AND T2 AVAILABLE                     
      IF(NR-IREC.LT.2)THEN                                                    
      IREC=IREC-1                                                             
      GOTO 110                                                                
      ENDIF                                                                   
      IF(FIRST)THEN                                                           
      IF(IREC.EQ.1)GOTO 135                                                   
      DO 130 I=1,IREC-1                                                       
      READ(LUN,120)SKIP                                                       
  120 FORMAT(A1)                                                              
  130 CONTINUE                                                                
      GOTO 135                                                                
      ENDIF                                                                   
      IDIF=IREC-IRECOL-1                                                      
      IF(IDIF.EQ.0)GOTO 135                                                   
      ISKP=ABS(IDIF)                                                          
      IF(IDIF.LT.0)THEN                                                       
      DO 131 I=1,ISKP*3                                                       
      BACKSPACE LUN                                                           
  131 CONTINUE                                                                
      ELSE                                                                    
      DO 132 I=1,ISKP                                                         
      READ(LUN,120)SKIP                                                       
  132 CONTINUE                                                                
      ENDIF                                                                   
  135 CONTINUE                                                                
      IRECOL=IREC+2                                                           
      DO 150 I=1,3                                                            
      READ(LUN,140)T(I),(PGI(I,J),J=1,3),(PHI(I,J),J=1,3),                    
     &                  (PEI(I,J),J=1,3),(PSI(I,J),J=1,3)                     
  140 FORMAT(5D16.10)                                                         
  150 CONTINUE                                                                
*                                                                      *      
*                                                                      *      
************************************************************************      
*                                                                      *      
*******  DETERMINE POSITIONS USING QUADRATIC INTERPOLATION   **BLOCK0200      
*                                                                      *      
      P=(ET-T(1))/H                                                           
      TMP(1)=0.5*(P-1.D0)*(P-2.D0)                                            
      TMP(2)=P*(2.D0-P)                                                       
      TMP(3)=0.5D0*P*(P-1.D0)                                                 
      DO 220 I=1,3                                                            
      DO 210 J=1,3                                                            
      PVG(I)=PVG(I)+PGI(J,I)*TMP(J)                                           
      PVH(I)=PVH(I)+PHI(J,I)*TMP(J)                                           
      PVE(I)=PVE(I)+PEI(J,I)*TMP(J)                                           
      PVS(I)=PVS(I)+PSI(J,I)*TMP(J)                                           
  210 CONTINUE                                                                
  220 CONTINUE                                                                
*                                                                      *      
************************************************************************      
*                                                                      *      
******* DETERMINE VELOCITIES FROM DERIVATIVES OF POSITIONS   **BLOCK0300      
*                                                                      *      
      TMP(1)=(P-1.5D0)/HSEC                                                   
      TMP(2)=(2.D0-2.D0*P)/HSEC                                               
      TMP(3)=(P-0.5D0)/HSEC                                                   
      DO 320 I=1,3                                                            
      I1=I+3                                                                  
      DO 310 J=1,3                                                            
      PVG(I1)=PVG(I1)+PGI(J,I)*TMP(J)                                         
      PVH(I1)=PVH(I1)+PHI(J,I)*TMP(J)                                         
      PVE(I1)=PVE(I1)+PEI(J,I)*TMP(J)                                         
      PVS(I1)=PVS(I1)+PSI(J,I)*TMP(J)                                         
  310 CONTINUE                                                                
  320 CONTINUE                                                                
*                                                                      *      
************************************************************************      
*                                                                      *      
************             ERROR MESSAGE BLOCK                ***BLOCK0900      
*                                                                      *      
      GOTO 1000                                                               
  910 CONTINUE                                                                
      IERR=1                                                                  
      WRITE(6,911)                                                            
  911 FORMAT(/' GIOPOS :  ET TOO EARLY ')                                     
      GOTO 1000                                                               
  920 CONTINUE                                                                
      IERR=2                                                                  
      WRITE(6,921)                                                            
  921 FORMAT(/' GIOPOS :  ET TOO LATE ')                                      
 1000 CONTINUE                                                                
      RETURN                                                                  
      END                                                                     
*************************************************************************

C-----------------------------------------------------------------------       
      SUBROUTINE TRAVEL(TUTC,S,DTG,P,G,IER)                                    
C-------------------------------------------------                             
C   ISSUE 01 - 18.01.85                                                        
C----------------------                                                        
C                                                                              
************        SUBROUTINE IDENTIFICATION               ************       
*                                                                      *       
*   TRAVEL PROVIDES THE FOLLOWING GIOTTO EXPERIMENTER SUPPORT DATA:    *       
*   - SIGNAL TRAVEL TIME BETWEEN THE GIOTTO S/C AND A STATION ON EARTH *       
*   - POSITIONS OF STATION AND S/C IN THE HELIOCENTRIC EME50           *       
*     AT SIGNAL RECEPTION AND TRANSMISSION TIMES INVOLVED.             *       
*   STATION COORDINATES ARE COMPUTED WITHOUT PRECESSION, I.E.          *       
*   THE TRAVEL TIMES GIVEN ARE ACCURATE TO ABOUT 1 MSEC.               *       
*   TRAVEL RETRIEVES THE STATES OF THE PARTICIPATING BODIES FROM       *       
*   THE POSITION DATA FILE BY MEANS OF THE INTERFACE ROUTINE GIODAT.   *       
*                                                                      *       
************************************************************************       
*                                                                      *       
************        VARIABLE IDENTIFICATIONS               *************       
*                                                                      *       
*****            INPUT VARIABLES ( ARGUMENT LIST )                   ***       
*                                                                      *       
*   TUTC      -  SIGNAL RECEPTION TIME (UTC DAYS SINCE 1950.0)         *       
*   S(1-3)    -  CARTESIAN COORDINATES OF STATION IN BODY-FIXED SYSTEM *       
*                                                                      *       
*****            OUTPUT VARIABLES ( ARGUMENT LIST )                  ***       
*                                                                      *       
*   DTG       -  SIGNAL TRAVEL TIME(SEC) GIOTTO S/C - STATION          *       
*   P(1-3)    -  POSITION OF STATION AT TUTC IN HELIOCENTRIC EME50     *       
*   G(1-3)    -  POSITION OF GIOTTO AT TUTC-DTG IN HELIOCENTRIC EME50  *       
*   IER       -  ERROR CODE                                            *       
*                =0 : NO ERROR                                         *       
*                =1 : TUTC TOO EARLY, NO DATA ON POSITION DATA FILE    *       
*                =2 : TUTC TOO LATE , NO DATA ON POSITION DATA FILE    *       
*                =3 : NO CONVERGENCE IN SOLUTION OF LIGHT TIME EQUATION*       
*                                                                      *       
*****            LOCAL VARIABLES OF SPECIAL INTEREST                 ***       
*                                                                      *       
*   SG(6)        GIOTTO STATE IN BARYCENTRIC EME50                     *       
*   SH(6)     -  DITTO FOR HALLEYS COMET                               *       
*   SU(6)     -  DITTO FOR SUN                                         *       
*                                                                      *       
*****            CONSTANT IDENTIFICATION                             ***       
*                                                                      *       
*   RJ,RK,RL  -  COEFFICIENTS IN FORMULA FOR SIDERIAL ANGLE            *       
*   C         -  VELOCITY OF LIGHT                                     *       
*   TWOPI     -  TWO PI                                                *       
*                                                                      *       
************************************************************************       
*                                                                      *       
************  TYPE DECLARATION AND STORAGE ALLOCATION       ************       
*                                                                      *       
      IMPLICIT REAL*8(A-H,O-Z)                                                 
      DIMENSION S(3),P(3),G(3),SS(6),SG(6),SH(6),SU(6)                         
*                                                                      *       
************************************************************************       
*                                                                      *       
************  SET CONSTANTS                                 ***BLOCK0000       
      DATA RJ,RK,RL /0.2769194D0,100.0021359028D0,0.107523D-5/                 
      DATA C,TWOPI /299792.458D0,6.283185307179587D0/                          
*                                                                      *       
************************************************************************       
*                                                                      *       
*******    STATION STATE AT TUTC                             **BLOCK0100       
*                                                                      *       
      CALL GIODAT(TUTC,SG,SH,SS,SU,IER)                                        
      IF(IER.NE.0) RETURN                                                      
      Y=0.5D0+TUTC/36525.D0                                                    
      IT=TUTC                                                                  
      X=TUTC-IT                                                                
      X=(X+RJ+(RK+RL*Y)*Y)*TWOPI                                               
      Y=DCOS(X)                                                                
      Z=DSIN(X)                                                                
      P(1)=S(1)*Y-S(2)*Z+SS(1)                                                 
      P(2)=S(1)*Z+S(2)*Y+SS(2)                                                 
      P(3)=S(3)+SS(3)                                                          
      G(1)=SU(1)                                                               
      G(2)=SU(2)                                                               
      G(3)=SU(3)                                                               
*                                                                      *       
************************************************************************       
*                                                                      *       
*******  SOLUTION OF LIGHT TIME EQUATION                     **BLOCK0200       
*                                                                      *       
      DTG=0.D0                                                                 
      DO 210 I=1,5                                                             
      X=SG(1)-P(1)                                                             
      Y=SG(2)-P(2)                                                             
      Z=SG(3)-P(3)                                                             
      DGB=DSQRT(X*X+Y*Y+Z*Z)                                                   
      Z=(X*SG(4)+Y*SG(5)+Z*SG(6))/DGB                                          
      Y=(DTG-DGB/C)/(1.D0-Z/C)                                                 
      DTG=DTG-Y                                                                
      IF(DABS(Y).LT.1.D-4) GOTO 1000                                           
      X=TUTC-DTG/864.D2                                                        
      CALL GIODAT(X,SG,SH,SS,SU,IER)                                           
      IF(IER.NE.0) RETURN                                                      
  210 CONTINUE                                                                 
*                                                                      *       
************************************************************************       
*                                                                      *       
*******  ERROR MESSAGES                                      **BLOCK0900       
*                                                                      *       
      IER=3                                                                    
      WRITE(6,910) IER                                                         
  910 FORMAT(' TRAVEL ERROR MESSAGE',I3/                                       
     &' = NO CONVERGENCE AFTER 5 ITERATIONS'/                                  
     &'   IN SOLUTION OF LIGHT TIME EQUATION')                                 
*                                                                      *       
************************************************************************       
*                                                                      *       
*******  SET OUTPUT VARIABLES AND RETURN                     **BLOCK1000       
*                                                                      *       
 1000 P(1)=P(1)-G(1)                                                           
      P(2)=P(2)-G(2)                                                           
      P(3)=P(3)-G(3)                                                           
      G(1)=SG(1)-SU(1)                                                         
      G(2)=SG(2)-SU(2)                                                         
      G(3)=SG(3)-SU(3)                                                         
      RETURN                                                                   
      END                                                                      
