*** SOLSYS VERSION 1 PACKAGE: SOLSYS, FILDEF, IDSS, BLOCK DATA *** SUBROUTINE SOLSYS (TJD,M,K,POS,VEL,IERR) * * SUBROUTINE SOLSYS VERSION 1. * THIS SUBROUTINE READS A COORDINATE FILE CONTAINING BARYCENTRIC * POSITIONS OF SOLAR SYSTEM BODIES AT DAILY INTERVALS AND PROVIDES * THE POSITION AND VELOCITY OF BODY M AT EPOCH TJD. * * TJD = TDB JULIAN DATE OF DESIRED EPOCH (IN) * M = BODY IDENTIFICATION NUMBER (IN) * K = ORIGIN SELECTION CODE (IN) * SET K=0 FOR ORIGIN AT SOLAR SYSTEM BARYCENTER * SET K=1 FOR ORIGIN AT CENTER OF SUN * POS = POSITION VECTOR, EQUATORIAL RECTANGULAR * COORDINATES, REFERRED TO ICRS, COMPONENTS IN AU (OUT) * VEL = VELOCITY VECTOR, EQUATORIAL RECTANGULAR * COORDINATES, REFERRED TO ICRS, COMPONENTS IN AU/DAY (OUT) * IERR = ERROR INDICATOR (OUT) * IERR=0 MEANS EVERYTHING OK * IERR=1 MEANS TJD BEFORE FIRST USABLE DATE IN FILE * IERR=2 MEANS TJD AFTER LAST USABLE DATE IN FILE * IERR=3 MEANS BAD VALUE OF M * IERR=4 MEANS PROBLEM OPENING FILE * * NOTE 1: INFORMATION ON THE COORDINATE FILE READ IN: * - PATH/NAME OF FILE SPECIFIED IN COMMON /SSFILE/ * - ALL RECORDS ASCII (FORMATTED) * - FIRST RECORD: HEADER OR IDENTIFYING INFORMATION, IGNORED * HERE * - OTHER RECORDS: TDB JULIAN DATE, X,Y,Z COORDINATES OF SUN, * X,Y,Z COORDINATES OF MERCURY, X,Y,Z COORDINATES OF VENUS, * X,Y,Z COORDINATES OF EARTH, X,Y,Z COORDINATES OF MARS, ... * READ IN AS PER FORMAT IN COMMON /SSFILE/ * - RECORDS AT FIXED INTERVALS OF +1 DAY OF TDB * - X,Y,Z VALUES IN AU WITH RESPECT TO BCRS (ICRS AXES) * - FILE MUST CONTAIN AT LEAST THE COORDINATES OF THE SUN AND * THE EARTH * - EARTH REFERS TO GEOCENTER, NOT EARTH-MOON BARYCENTER * MANY ASPECTS OF THE FILE CAN BE CONTROLLED AT EXECUTION TIME VIA * SUBROUTINE FILDEF. DEFAULTS: FILE PATH/NAME 'SS_EPHEM.TXT' * READ ON LOGICAL UNIT 20, CONTAINING, IN EACH RECORD, A TDB * JULIAN DATE AND COORDINATES OF 11 BODIES (SUN, MERCURY, VENUS, * EARTH, ..., PLUTO, MOON), READ IN WITH FORMAT (F10.2,11(3F16.12)). * * NOTE 2: IN SUCCESSIVE CALLS TO THIS SUBROUTINE, INPUT * JULIAN DATES (TJD) SHOULD GENERALLY BE IN ASCENDING ORDER * TO AVOID MULTIPLE SEARCHES STARTING AT BEGINNING OF FILE. * * NOTE 3: THIS SUBROUTINE USES A 7-POINT LAGRANGIAN INTERPOLATION * SCHEME ON FIXED INTERVAL DATA, WHICH PRODUCES INTERPOLATION ERRORS * THAT VARY WITH BODY. * * DOUBLE PRECISION TJD,POS,VEL,XJD,XYZ,BPOS,BVEL,T, . TBEG,TEND,TLAST,ASTART,ORIGIN,TI,AK,AI,P, . DABS,DFLOAT CHARACTER FILNAM*80, FORMT*80 DIMENSION POS(3), VEL(3), XJD(13), XYZ(13,50,3), . BPOS(13,3), BVEL(13,3) * COMMON BLOCK SSFILE CONTAINS INFORMATION ON THE COORDINATE FILE COMMON /SSFILE/ LU,N,FILNAM,FORMT SAVE DATA TBEG, TEND, TLAST, MLAST, KLAST / 0.D0, 1.D10, 0.D0, 0, 0 / 1 FORMAT ( A ) 3 FORMAT ( ' SOLSYS: ERROR ',I1,' AT JD ', F10.1, ', BODY ID ', . I2 ) 4 FORMAT ( ' SOLSYS: ERROR 4 TRYING TO OPEN FILE ', A, ' ON UNIT ', . I2 ) IF ( TLAST .LE. 0.D0 ) THEN OPEN ( UNIT=LU, FILE=FILNAM, STATUS='UNKNOWN', ERR=84 ) READ ( LU, 1 ) READ ( LU, FORMT ) T REWIND LU READ ( LU, 1 ) MSUN = IDSS ( 'SUN' ) TLAST = 1.D0 TBEG = T + 8.D0 NPTS = 13 INTPTS = 7 LMIDDL = NPTS / 2 + 1 LSTART = LMIDDL - INTPTS/2 - 1 ASTART = LSTART END IF * CHECK FOR COMMON ERROR CONDITIONS IERR = 0 IF ( M .LT. 0 .OR. M .GT. N - 1 ) GO TO 73 IF ( TJD .LT. TBEG ) GO TO 71 IF ( TJD .GT. TEND ) GO TO 78 * LOGIC TO DETERMINE BEST WAY TO SEARCH COORDINATE FILE * CHECK IF NEEDED DATA ALREADY IN ARRAYS IF ( DABS ( TJD - TLAST ) .LE. 0.8D0 ) THEN IF ( M .NE. MLAST .OR. K .NE. KLAST ) GO TO 30 GO TO 60 END IF * IF INPUT JD LESS THAN LAST JD, START FROM BEGINNING OF FILE IF ( TJD .LT. TLAST ) THEN REWIND LU READ ( LU, 1 ) TLAST = 1.D0 END IF * DECIDE ON COURSE OR FINE SEARCH IF ( TJD - TLAST .LE. 15.D0 ) GO TO 20 * COURSE SEARCH THROUGH COORDINATE FILE 15 READ ( LU, FORMT, END=77 ) T IF ( TJD - T .GT. 10.D0 ) GO TO 15 DO 18 I = 1, NPTS READ ( LU, FORMT, END=77 ) T, ((XYZ(I,L,J), J=1,3), L=1,N) XJD(I) = T 18 CONTINUE * FINE SEARCH THROUGH COORDINATE FILE 20 DO 22 I = 1, NPTS - 1 IOLD = I + 1 XJD(I) = XJD(IOLD) DO 21 L = 1, N DO 21 J = 1, 3 XYZ(I,L,J) = XYZ(IOLD,L,J) 21 CONTINUE 22 CONTINUE READ ( LU, FORMT, END=77 ) T, ((XYZ(NPTS,L,J), J=1,3), L=1,N) XJD(NPTS) = T IF ( DABS ( TJD - XJD(LMIDDL) ) .GT. 0.5D0 ) GO TO 20 TLAST = XJD(LMIDDL) * AT THIS POINT, THE FILE IS POSITIONED CORRECTLY AND ARRAYS * XJD AND XYZ ARE FILLED IN 30 CONTINUE * FILL ARRAY BPOS WITH DAILY POSITIONS OF BODY M (WITH INDEX = M+1) * IF K=1, MOVE ORIGIN TO SUN DO 40 I = 1, NPTS DO 40 J = 1, 3 ORIGIN = 0.D0 IF ( K .GE. 1 ) ORIGIN = XYZ(I,MSUN+1,J) BPOS(I,J) = XYZ(I,M+1,J) - ORIGIN 40 CONTINUE * FILL ARRAY BVEL WITH DAILY VELOCITIES OF BODY M * COMPUTED FROM NUMERICAL DIFFERENTIATION OF POSITIONS IN ARRAY BPOS DO 50 I = 1, NPTS DO 50 J = 1, 3 BVEL(I,J) = 0.D0 IF ( I .GE. 4 .AND. I .LE. 10 ) . BVEL(I,J) = ( BPOS(I+3,J) - 9.D0 * BPOS(I+2,J) . + 45.D0 * BPOS(I+1,J) - 45.D0 * BPOS(I-1,J) . + 9.D0 * BPOS(I-2,J) - BPOS(I-3,J) ) . / 60.D0 50 CONTINUE MLAST = M KLAST = K * PERFORM LAGRANGIAN INTERPOLATION FOR POSITION AND VELOCITY AT * EPOCH TJD 60 TI = TJD - XJD(LMIDDL) + LMIDDL DO 63 J = 1, 3 POS(J) = 0.D0 VEL(J) = 0.D0 DO 62 L = 1, INTPTS AK = ASTART + DFLOAT(L) P = 1.D0 DO 61 I = 1, INTPTS IF ( I .EQ. L ) GO TO 61 AI = ASTART + DFLOAT(I) P = P * (TI-AI) / (AK-AI) 61 CONTINUE POS(J) = POS(J) + P * BPOS(LSTART+L,J) VEL(J) = VEL(J) + P * BVEL(LSTART+L,J) 62 CONTINUE 63 CONTINUE GO TO 99 71 IERR = 1 GO TO 80 73 IERR = 1 GO TO 80 77 TEND = T - 6.D0 REWIND LU READ ( LU, 1 ) TLAST = 1.D0 78 IERR = 2 80 WRITE ( *, 3 ) IERR, TJD, M GO TO 99 84 IERR = 4 WRITE ( *, 4 ) FILNAM, LU 99 RETURN END SUBROUTINE FILDEF (LUN,NBOD,FILNM,FMT) * * FOR USE WITH SUBROUTINE SOLSYS VERSION 1. * THIS SUBROUTINE MAY BE CALLED TO CHANGE THE VALUES IN * COMMON BLOCK SSFILE, WHICH CONTAINS INFORMATION ON THE * COORDINATE FILE USED BY SUBROUTINE SOLSYS. * * LUN = FORTRAN LOGICAL UNIT NUMBER TO BE USED FOR * COORDINATE FILE (IN) * NBOD = NUMBER OF BODIES WITH COORDINATES IN FILE (IN) * FILNM = CHARACTER VARIABLE CONTAINING PATH AND FILE NAME * OF COORDINATE FILE (IN) * FMT = CHARACTER VARIABLE CONTAINING FORMAT STATEMENT, * INCLUDING PARENTHESES AND EVERYTHING IN BETWEEN (IN) * * NOTE: IF LUN OR NBOD IS ZERO OR LESS, OR FILNM OR FMT IS BLANK, * THE CORRESPONDING DEFAULT VALUE IS NOT CHANGED. DEFAULT VALUES * ARE SET IN BLOCK DATA FOR COMMON /SSFILE/. * * CHARACTER FILNM*(*), FMT*(*), FILNAM*80, FORMT*80 COMMON /SSFILE/ LU,N,FILNAM,FORMT SAVE IF ( LUN .GE. 1 ) LU = LUN IF ( NBOD .GE. 1 ) N = NBOD IF ( FILNM .NE. ' ' ) FILNAM = FILNM IF ( FMT .NE. ' ' ) FORMT = FMT RETURN END INTEGER FUNCTION IDSS ( NAME ) * * FOR USE WITH SOLSYS VERSION 1. * THIS FUNCTION RETURNS THE ID NUMBER OF A SOLAR SYSTEM BODY * FOR THE VERSION OF SOLSYS (OR SOLSYS-AUXPOS COMBINATION) IN USE. * FOR SOLSYS VERSION 1, THE ID NUMBER OF A BODY REFERS TO ITS * ORDER WITHIN EACH RECORD OF THE COORDINATE FILE, WITH ID NUMBERS * BEGINNING AT 0 FOR THE FIRST BODY (NORMALLY THE SUN). * * NAME = NAME OF BODY WHOSE ID NUMBER IS DESIRED, E.G., * 'SUN', 'MOON, 'MERCURY', ETC., EXPRESSED AS ALL * UPPER-CASE LETTERS (IN) * IDSS = ID NUMBER OF BODY, FOR USE IN CALLS TO SOLSYS * (FUNCTION VALUE RETURNED) * * NOTE 1: IN THIS VERSION, ONLY THE FIRST THREE LETTERS OF THE * BODY'S NAME ARE USED FOR IDENTIFICATION. ALTERNATIVE VERSIONS * MIGHT USE MORE LETTERS. * * NOTE 2: IF NAME IS 'JD', IDSS RETURNS IDSS=1, SINCE SOLSYS * VERSION 1 DOES NOT PROCESS SPLIT JULIAN DATES. * * NOTE 3: ALL VERSIONS OF IDSS MUST RETURN IDSS=-9999 FOR OBJECTS * THAT IT CANNOT IDENTIFY OR ARE UNSUPPORTED BY SOLSYS. * * CHARACTER NAME*(*), NAMEIN*3, NAMES*3 DIMENSION NAMES(50), IDS(50) DATA NAMES / 'SUN', 'MER', 'VEN', 'EAR', 'MAR', 'JUP', 'SAT', . 'URA', 'NEP', 'PLU', 'MOO', '---', '---', '---', . '---', '---', '---', '---', '---', '---', '---', . '---', '---', '---', '---', '---', '---', '---', . '---', '---', '---', '---', '---', '---', '---', . '---', '---', '---', '---', '---', '---', '---', . '---', '---', '---', '---', '---', '---', '---', . '---'/ DATA IDS / 0, 1, 2, 3, 4, 5, 6, . 7, 8, 9, 10, 0, 0, 0, . 0, 0, 0, 0, 0, 0, 0, . 0, 0, 0, 0, 0, 0, 0, . 0, 0, 0, 0, 0, 0, 0, . 0, 0, 0, 0, 0, 0, 0, . 0, 0, 0, 0, 0, 0, 0, . 0/ DATA NUM / 11 / 3 FORMAT ( ' IDSS ERROR: NO BODY ID NUMBER FOUND FOR ', A ) IDSS = -9999 NAMEIN = NAME * LOOK THROUGH LIST OF BODY NAMES TO FIND MATCH DO 20 I = 1, NUM IF ( NAMEIN .EQ. NAMES(I) ) THEN IDSS = IDS(I) GO TO 30 END IF 20 CONTINUE * IF NO MATCH, CHECK FOR INQUIRY ABOUT SPLIT JULIAN DATES IF ( NAMEIN .EQ. 'JD ' ) THEN * IN THIS CASE, SET IDSS=1, SINCE SOLSYS VERSION 1 DOES NOT * PROCESSES SPLIT JULIAN DATES (IN SUCCESSIVE CALLS) IDSS = 1 GO TO 30 END IF WRITE ( *, 3 ) NAME 30 RETURN END BLOCK DATA * * FOR USE WITH SUBROUTINE SOLSYS VERSION 1. * COMMON BLOCK /SSFILE/ CONTAINS INFORMATION ON THE COORDINATE FILE * USED BY SUBROUTINE SOLSYS. THIS BLOCK DATA SEGMENT SETS UP THE * DEFAULT VALUES FOR THE PARAMETERS IN /SSFILE/. THESE DEFAULTS CAN * BE ALTERED BY EXECUTABLE STATEMENTS IN THE MAIN PROGRAM OR ANY * SUBROUTINE, OR BY A CALL TO SUBROUTINE FILDEF. * * LU = FORTRAN LOGICAL UNIT NUMBER OF COORDINATE FILE * N = NUMBER OF BODIES WITH COORDINATES IN FILE * FILNAM = CHARACTER VARIABLE CONTAINING PATH AND FILE NAME * OF COORDINATE FILE * FORMT = CHARACTER VARIABLE CONTAINING FORMAT STATEMENT, * INCLUDING PARENTHESES AND EVERYTHING BETWEEN * * CHARACTER FILNAM*80, FORMT*80 COMMON /SSFILE/ LU,N,FILNAM,FORMT DATA LU / 20 / DATA N / 11 / DATA FILNAM / 'SS_EPHEM.TXT' / DATA FORMT / '(F10.2,33F16.12)' / END