
C************************************************************************
      PROGRAM GIOCALL
C-----------------------------------------------------------------------
C     PROGRAMMER J.D. KOCH 94-05-28
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.  THE CHANGES INVOLVE
C		COMMENTING CODE THAT THE MS-FORTRAN COMPILER DOES
C		NOT ACCEPT AND UNCOMMENTING CORRESPONDING CODE THAT
C		THE MS-FORTRAN COMPILER DOES ACCEPT.
C----------------------------------------------------------------------
************        PROGRAM PURPOSE                         ************
*     THIS PROGRAM IS INTENDED MERELY TO TEST THE SUBROUTINE GIODAT    *
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER I,J,N,M,I1,I2,I4,I5,ONE
      CHARACTER  K*1,B*1,H*1,G*1,C*1
      DIMENSION PVG(6),PVH(6),PVE(6),PVS(6),HVG(6),HVH(6),HVE(6)
      DIMENSION EVG(6),EVH(6),CVG(6)
      OPEN(UNIT=7,FILE='POSITION',ACCESS='SEQUENTIAL',FORM='FORMATTED')
      READ(7,2) I1,I2,D3,I4,I5,TMIN,TMAX,SUTC,ETUTC1,ETUTC2
 2    FORMAT(2I5,D16.10,2I5,/,5D16.10)
      CLOSE(UNIT=7)
      OPEN(UNIT=7,FILE='POSITION',ACCESS='SEQUENTIAL',FORM='FORMATTED')
      TUTC = 1000
      ONE = 1
      PRINT *,'OPTIONS: ENTER t TO FIND THE TIME FOR GIVEN DISTANCE,'
      PRINT *,'OR p TO FIND THE POSITION OF COMET/GIOTTO AT A',
     &' GIVEN TIME.'
      READ *, K
* DOS CHANGES:	COMMENT THE NEXT 2 LINES
       DO 70 I = 1,100
	 I=1
* DOS CHANGES:	UNCOMMENT THE NEXT 3 LINES
*	 I = 1
*	 DO 70 WHILE (I .LE. 100)
*	 I = 1
         IF(K .EQ. 't') THEN
	    E = 0.00001
* DOS CHANGES:	COMMENT THE NEXT 2 LINES
	    DO 15 J = 1,100
	       IF (ONE) THEN
* DOS CHANGES:	UNCOMMENT THE NEXT 3 LINES
*	    J = 1
*	    DO 15 WHILE (J .LE. 100)
*	       IF (ONE .EQ. 1) THEN
                  ONE=0
                  TLOW = TMIN+1
                  THIGH = TMAX-1
                  CALL GIODAT(TMIN+1,PVG,PVH,PVE,PVS,IERR)
                  D1 = (PVG(1)-PVH(1))**2
                  D2 = (PVG(2)-PVH(2))**2
                  D3 = (PVG(3)-PVH(3))**2
                  DMIN = (D1+D2+D3)**0.5
                  CALL GIODAT(TMAX-1,PVG,PVH,PVE,PVS,IERR)
                  D1 = (PVG(1)-PVH(1))**2
                  D2 = (PVG(2)-PVH(2))**2
                  D3 = (PVG(3)-PVH(3))**2
                  DMAX = (D1+D2+D3)**0.5
                  DLOW = DMIN
                  DHIGH = DMAX
                  IF (DMAX .LT. DMIN) THEN
                     TLOW = TMAX-1
                     THIGH = TMIN+1
                     DLOW = DMAX
                     DHIGH = DMIN
		  ENDIF
* DOS CHANGES:	COMMENT THE NEXT LINE
		  DO 5 M = 1,10000
* DOS CHANGES:	UNCOMMENT THE NEXT 2 LINES
*		  M = 1
*		  DO 5 WHILE (M .LE. 10000)
                     TMID = (THIGH + TLOW)/2.0
                     CALL GIODAT(TMID,PVG,PVH,PVE,PVS,IERR)
                     D1 = (PVG(1)-PVH(1))**2
                     D2 = (PVG(2)-PVH(2))**2
                     D3 = (PVG(3)-PVH(3))**2
                     DMID = (D1+D2+D3)**0.5
                     IF (DMID .LT. DLOW) THEN 
                        THIGH = TLOW
                        TLOW = TMID
                        DHIGH = DLOW
                        DLOW = DMID
                     ELSE
                        THIGH = TMID
                        DHIGH = DMID
                     ENDIF
		     IF (ABS(DHIGH-DLOW) .LT. E) M = 10002
* DOS CHANGES:	UNCOMMENT THE NEXT LINE
*		     M = M + 1
 5                CONTINUE
                  TCA = TLOW
                  DCA = DLOW
                  PRINT *,'AT TIME OF CLOSEST APPROACH',TCA,' DAYS ',
     &'EME50'
                  PRINT *,'    GIOTTO WAS ',DCA,' KM FROM',' COMET'
                  PRINT *,' '
               ENDIF
               J = 1
               K = 'r'
               PRINT *,'ENTER a IF GIOTTO IS APPROACHING COMET, r IF',
     &' RECEEDING'
               READ *, K
               PRINT *,' '
               PRINT *,'ENTER DISTANCE FROM GIOTTO TO COMET IN KM'
               PRINT *,' '
               READ *,DIST
               PRINT *,' '
               TLOW = TCA
               THIGH = TMAX-1
	       IF (K .EQ. 'a') THIGH = TMIN+1
* DOS CHANGES:	COMMENT THE NEXT LINE
	       DO 10 M = 1,1000
* DOS CHANGES:	UNCOMMENT THE NEXT 2 LINES
*	       M = 1
*	       DO 10 WHILE (M .LE. 1000)
                  TMID = (THIGH + TLOW)/2.0
                  CALL GIODAT(TMID,PVG,PVH,PVE,PVS,IERR)
                  D1 = (PVG(1)-PVH(1))**2
                  D2 = (PVG(2)-PVH(2))**2
                  D3 = (PVG(3)-PVH(3))**2
                  D = (D1+D2+D3)**0.5
                  IF (M .EQ. 1000) THEN 
                     PRINT *,'No time for distance =',DIST,'km found.'
                     PRINT *,'Approximate time =',TMID,
     &' at a distance of',D,'km'
                  ELSE IF (DIST .LT. (D-E)) THEN 
                     THIGH = TMID
                  ELSE IF (DIST .GT. (D+E)) THEN 
                     TLOW = TMID
                  ELSE 
                     PRINT *,'AT TUTC',TMID,', GIOTTO WAS ',D,' KM FROM',
     &' COMET'
                     M = 1001
		  ENDIF
* DOS CHANGES:	UNCOMMENT THE NEXT LINE
*		  M = M + 1
 10	       CONTINUE
	       PRINT *, ' '
               PRINT *,'SEARCH FOR ANOTHER TIME OR POSITION?(t/p/q)'
               READ *, K
               IF (K .EQ. 'q') I=101
	       J = 101
* DOS CHANGES:	UNCOMMENT THE NEXT LINE
*	       J = J + 1
 15         CONTINUE
         ELSE IF(K .EQ. 'p') THEN
            B = 'y'
            H = 'y'
            G = 'y'
            C = 'y'
            PRINT *,'BARYCENTRIC, HELIOCENTRIC, GEOCENTRIC, AND/OR', 
     &' COMETOCENTRIC COORDINATES? ','(y/n,y/n,y/n,y/n)'
	    READ *,B,H,G,C
* DOS CHANGES:	COMMENT THE NEXT 2 LINES
	    DO 65 J = 1,100
	       J=2
* DOS CHANGES:	UNCOMMENT THE NEXT 3 LINES
*	    J = 1
*	    DO 65 WHILE (J .LE. 100)
*	       J = 2
               PRINT *,' '
               PRINT *,'ENTER SIGNAL RECEPTION TIME(-999 T0 QUIT)'
               READ *, TUTC
               IF(TUTC .NE. -999)THEN
                  CALL GIODAT(TUTC,PVG,PVH,PVE,PVS,IERR)
                  IF(IER .EQ. 0) THEN
                     DO 60 N = 1,6
                        HVG(N) = PVG(N)-PVS(N)
                        HVH(N) = PVH(N)-PVS(N)
                        HVE(N) = PVE(N)-PVS(N)
                        EVG(N) = PVG(N)-PVE(N)
                        EVH(N) = PVH(N)-PVE(N)
                        CVG(N) = PVG(N)-PVH(N)
 60                  CONTINUE
                     PRINT *,' '
                     IF(B .NE. 'n') THEN
                       PRINT *,'BARYCENTRIC GIOTTO POSITION (X,Y,Z)'
                       PRINT *,PVG(1),' KM ',PVG(2),' KM ',PVG(3),' KM'
                       PRINT *,'BARYCENTRIC GIOTTO VELOCITY (Vx,Vy,Vz)'
                       PRINT *,PVG(4),' KM/S ',PVG(5),' KM/S ',PVG(6),
     &' KM/S'
                       PRINT *,'BARYCENTRIC COMET POSITION (X,Y,Z)'
                       PRINT *,PVH(1),' KM ',PVH(2),' KM ',PVH(3),' KM'
                       PRINT *,'BARYCENTRIC COMET VELOCITY (Vx,Vy,Vz)'
                       PRINT *,PVH(4),' KM/S ',PVH(5),' KM/S ',PVH(6),
     &' KM/S'
                       PRINT *,'BARYCENTRIC EARTH POSITION (X,Y,Z)'
                       PRINT *,PVE(1),' KM ',PVE(2),' KM ',PVE(3),' KM'
                       PRINT *,'BARYCENTRIC EARTH VELOCITY (Vx,Vy,Vz)'
                       PRINT *,PVE(4),' KM/S ',PVE(5),' KM/S ',PVE(6),
     &' KM/S'
                       PRINT *,' '    
                     ENDIF
                     IF (H .NE. 'n') THEN
                       PRINT *,'HELIOCENTRIC GIOTTO POSITION (X,Y,Z)'
                       PRINT *,HVG(1),' KM ',HVG(2),' KM ',HVG(3),' KM'
                       PRINT *,'HELIOCENTRIC GIOTTO VELOCITY (Vx,Vy,Vz)'
                       PRINT *,HVG(4),' KM/S ',HVG(5),' KM/S ',HVG(6),
     &' KM/S'
                       PRINT *,'HELIOCENTRIC COMET POSITION (X,Y,Z)'
                       PRINT *,HVH(1),' KM ',HVH(2),' KM ',HVH(3),' KM'
                       PRINT *,'HELIOCENTRIC COMET VELOCITY (Vx,Vy,Vz)'
                       PRINT *,HVH(4),' KM/S ',HVH(5),' KM/S ',HVH(6),
     &' KM/S'
                       PRINT *,'HELIOCENTRIC EARTH POSITION (X,Y,Z)'
                       PRINT *,HVE(1),' KM ',HVE(2),' KM ',HVE(3),' KM'
                       PRINT *,'HELIOCENTRIC EARTH VELOCITY (Vx,Vy,Vz)'
                       PRINT *,HVE(4),' KM/S ',HVE(5),' KM/S ',HVE(6),
     &' KM/S'
                       PRINT *,' '
                     ENDIF
                     IF (G .NE. 'n') THEN
                       PRINT *,'GEOCENTRIC GIOTTO POSITION (X,Y,Z)'
                       PRINT *,EVG(1),' KM ',EVG(2),' KM ',EVG(3),' KM'
                       PRINT *,'GEOCENTRIC GIOTTO VELOCITY (Vx,Vy,Vz)'
                       PRINT *,EVG(4),' KM/S ',EVG(5),' KM/S ',EVG(6),
     &' KM/S'
                       PRINT *,'GEOCENTRIC COMET POSITION (X,Y,Z)'
                       PRINT *,EVH(1),' KM ',EVH(2),' KM ',EVH(3),' KM'
                       PRINT *,'GEOCENTRIC COMET VELOCITY (Vx,Vy,Vz)'
                       PRINT *,EVH(4),' KM/S ',EVH(5),' KM/S ',EVH(6),
     &' KM/S'
                       PRINT *,' '
                     ENDIF  
                     IF (C .NE. 'n') THEN
                       PRINT *,'COMETOCENTRIC GIOTTO POSITION (X,Y,Z)'
                       PRINT *,CVG(1),' KM ',CVG(2),' KM ',CVG(3),' KM'
                       PRINT *,'COMETOCENTRIC GIOTTO VELOCITY(Vx,Vy,Vz)'
                       PRINT *,CVG(4),' KM/S ',CVG(5),' KM/S ',CVG(6),
     &' KM/S'
                       PRINT *,' '
                     ENDIF  
                     D1 = (PVG(1)-PVH(1))**2
                     D2 = (PVG(2)-PVH(2))**2
                     D3 = (PVG(3)-PVH(3))**2
                     D = (D1+D2+D3)**0.5
                     PRINT *,'DISTANCE FROM GIOTTO TO COMET'
                     PRINT *,D,' KM'
                     PRINT *,' '
                  ENDIF
               ELSE
                  PRINT *,'SEARCH FOR ANOTHER TIME OR POSITION?(t/p/q)'
                  READ *, K
                  IF (K .EQ. 'q') I=101
                  J = 101
	       ENDIF
* DOS CHANGES:	UNCOMMENT THE NEXT LINE
*	       J = J + 1
 65         CONTINUE
         ELSE 
            I=101 
	 ENDIF
* DOS CHANGES:	UNCOMMENT THE NEXT LINE
*	 I = I + 1
 70   CONTINUE
      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, THE TARGET 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 THE TARGET 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 THE TARGET 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                                                                     
*************************************************************************
