C$Procedure LSRGP_TEST C PROGRAM LSRGPJ_TEST C C$ Log C C 19-JUN-1990 - Eric Cannell - creation of LSRGP_TEST C C$ Purpose C C LSRGP_TEST exercises the LSRGP library. Bosically, LSRGP_TEST propagates C an orbit for a specified number of days and then plots the mean C elements during that time. C C$ Namelist_Input C C Namelist: $INPUT C Name Type Dim Default Units Description C ----------------------------------------------------------------------------- C DATE C*25 1 - --> epoch of O0 in TIMETRANS C character format C DAYS I 1 20 days the approximate number of days to C propagate O0 C DT DP 1 1000.0 sec time step size C LSFLAG L 1 T - if true, luni-solar effects are on C LTOP I 1 0 - LTOP in LSRGP call C O0 DP 6 Tpx Ref #5 --> classical orbital elements, a, e, C i, LAN, w, M in km and deg. C OFILE C*20 1 'OUTPUT.LIS' - output file name C PFILE C*13 1 'PGPLOT.IMPLOT' - PGPLOT file name C C$ Library_Links C C LSRGP C PGPLOT C TPXUTIL C C$ Files C C * - - standard I/O C 7 - IN_LSRGP_TEST - namelist_input_file C 8 - OFILE - output file C 9 - GRAVITY_FIELD - earth field C ? - - PGPLOT plot file C C$ Parameters C INTEGER MAXSTP PARAMETER ( MAXSTP = 20000 ) C C$ Declarations_of_Input_and_Output_Arguments C CHARACTER*25 DATE INTEGER DAYS DOUBLE PRECISION DT LOGICAL LSFLAG INTEGER LTOP DOUBLE PRECISION O0(6) CHARACTER*20 OFILE CHARACTER*13 PFILE C Via $ZFRAME. C CHARACTER*12 DEVICE LOGICAL SUMMRY CHARACTER*60 TITLE REAL WINDO ( 4 ) CHARACTER*23 XTITLE CHARACTER*35 YTITLE C C Via $ZLINE. C INTEGER NPTS INTEGER STYLE REAL X ( MAXSTP ) REAL Y ( MAXSTP ) C via $GRAVITY double precision jearth (2:29) C C$ Declarations_of_Local_Variables C REAL CPU0 REAL CPU1 REAL DE ( MAXSTP ) DOUBLE PRECISION DI ( MAXSTP ) REAL DO ( MAXSTP ) DOUBLE PRECISION DTDAY REAL DW ( MAXSTP ) LOGICAL ERROR INTEGER I DOUBLE PRECISION JD INTEGER NSTEPS DOUBLE PRECISION O1(6) DOUBLE PRECISION O2(6) CHARACTER*21 PGFILE REAL TIME ( MAXSTP ) REAL XHI REAL XLO REAL YHI REAL YLO C C$ External_Statements C DOUBLE PRECISION CH2JD EXTERNAL CH2JD REAL GETCPU EXTERNAL GETCPU CHARACTER*25 JD2CH EXTERNAL JD2CH INTEGER RMNI EXTERNAL RMNI INTEGER RMXI EXTERNAL RMXI REAL SETCPU EXTERNAL SETCPU C C$ Namelists C NAMELIST / INPUT / DATE , DAYS , DT , LSFLAG , LTOP , & O0 , OFILE , PFILE NAMELIST / ZFRAME / DEVICE , SUMMRY , TITLE , WINDO , XTITLE , & YTITLE NAMELIST / ZLINE / NPTS , STYLE , X , Y namelist /gravity/ jearth C C$ Data_Statements C DATA DAYS / 20 / DATA DT / 1000.000000 / DATA LSFLAG / .TRUE. / DATA LTOP / 0 / DATA O0 / 7714.408000 , & .000080 , & 66.018000 , & 155.101000 , & 90.000000 , & 0.000000 / DATA OFILE / 'OUTPUT.LIS' / DATA PFILE / 'PGPLOT.IMPLOT' / C C$ Method C-& C1> Use OPSFOR to open the input file. CALL OPSFOR ( 7 , 'IN_LSRGP_TEST' , ERROR ) IF ( ERROR ) THEN WRITE(*,'(/1X,''LSRGP_TEST: cannot open IN_LSRGP_TEST.'')') STOP END IF C1 Read the input namelist $INPUT and then close it. READ ( 7 , INPUT ) CLOSE ( 7 ) C1> Use OPSFOR to open the input file. CALL OPSFOR ( 9, 'GRAVITY_FIELD', ERROR ) IF ( ERROR ) THEN WRITE(*,'(/1X,''LSRGP_TEST: cannot open GRAVITY_FIELD.'')') STOP END IF C1 Read the input namelist $INPUT and then close it. READ ( 9 , GRAVITY) CLOSE ( 9 ) write (6, gravity ) C1> Use OPSFN to open the output file. CALL OPSFN ( 8 , OFILE , ERROR ) IF ( ERROR ) THEN WRITE(*,'(/1X,''LSRGP_TEST: cannot open '',A,''.'')') OFILE STOP END IF C1 Check that number of steps is less than MAXSTP. DTDAY = DT / 86400D0 NSTEPS = DFLOAT( DAYS ) / DTDAY + 1 IF ( NSTEPS .GT. MAXSTP ) THEN WRITE(8,'(//,1X,''LSRGP_TEST: Too many steps.'')') STOP END IF C1 Setup first step at time zero. JD = CH2JD( DATE , .TRUE. , ERROR ) TIME(1) = 0.0 DE (1) = O0(2) DI (1) = O0(3) DO (1) = O0(4) DW (1) = O0(5) C1 Copy initial orbit into working space. DO I = 1 , 6 O1(I) = O0(I) END DO C1 Set CPU timer. CPU0 = SETCPU() C1 Loop it. DO I = 2 , NSTEPS C2 Propagate the orbit. CALL LSRGPJ ( O1 , DATE , DT , JEARTH, LTOP , LSFLAG , O2 ) C2 Copy orbits. O1(1) = O2(1) O1(2) = O2(2) O1(3) = O2(3) O1(4) = O2(4) O1(5) = O2(5) O1(6) = O2(6) C2 Update epoch. JD = JD + DTDAY DATE = JD2CH( JD ) C2 Save the propagated orbit. TIME(I) = TIME(I-1) + DTDAY DE(I) = O2(2) DI(I) = O2(3) DO(I) = O2(4) DW(I) = O2(5) END DO C1 Get CPU time. CPU1 = GETCPU() PRINT * PRINT * PRINT *,'CPU time for LSRGP = ',CPU1-CPU0 PRINT * PRINT * C1 Print the results. 301 FORMAT('1',/,1X,' days |', & ' e', & ' i', & ' lan', & ' w', & /,1X,132('-')) DO I = 1 , NSTEPS IF ( MOD(I-1,50) .EQ. 0 ) WRITE(8,301) WRITE(8,302) TIME(I) , DE(I) , DI(I) , DO(I) , DW(I) 302 FORMAT(1X,F14.6,' |',4F29.15) END DO C1 Plot the results. DEVICE = '/TEK' SUMMRY = .FALSE. TITLE = 'LSRGP Test' WINDO(3) = DE( RMNI( DE , NSTEPS ) ) WINDO(4) = DE( RMXI( DE , NSTEPS ) ) WINDO(1) = 0 WINDO(2) = TIME(NSTEPS) XTITLE = 'DAYS' YTITLE = 'eccentricity' style = 1 NPTS = NSTEPS do j=1,NSTEPS Y(j) = de(j) x(j) = time(j) end do open ( unit = 11, name = 'eccentricity.plot', & status = 'NEW' ) WRITE(11,ZFRAME) WRITE(11,ZLINE) close ( 11 ) windo(3) = DO( RMNI( DO , NSTEPS ) ) windo(4) = DO( RMXI( DO , NSTEPS ) ) ytitle = 'LONGITUDE OF THE ASCENDING NODE via LSRGP Library' do j=1,NSTEPS Y(j) = do(j) end do open ( unit = 11, name = 'lan.plot', & status = 'NEW' ) WRITE(11,ZFRAME) WRITE(11,ZLINE) close ( 11 ) windo(3) = DW( RMNI( DW , NSTEPS ) ) windo(4) = DW( RMXI( DW , NSTEPS ) ) ytitle = 'ARGUMENT OF PERIAPSIS via LSRGP Library' do j=1,NSTEPS Y(j) = dw(j) end do open ( unit = 11, name = 'perigee.plot', & status = 'NEW' ) WRITE(11,ZFRAME) WRITE(11,ZLINE) close ( 11 ) windo(1) = DW( RMNI( DW , NSTEPS ) ) windo(2) = DW( RMXI( DW , NSTEPS ) ) windo(3) = DE( RMNI( DE , NSTEPS ) ) windo(4) = DE( RMXI( DE , NSTEPS ) ) xtitle = 'argument of periapsis, deg' ytitle = 'eccentricity' do j=1,NSTEPS x(j) = dw(j) y(j) = de(j) end do open ( unit = 11, name = 'evsw.plot', & status = 'NEW' ) WRITE(11,ZFRAME) WRITE(11,ZLINE) close ( 11 ) END