C San Marco Project -- Attitude Software C NSSDC Technical Reference File #B55111 (Entered 26 April 2006) C Provided in conjunction with archived DBI data. C---------------------------------------------------------------------- SUBROUTINE ULPAOR (IT) ULP00010 C ULP00020 C----------------------------------------------------------------------+ULP00030 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - 0012 |ULP00040 C----------------------------------------------------------------------|ULP00050 C Title : ULPAOR |ULP00060 C----------------------------------------------------------------------|ULP00070 C Function : ULPAOR gives the NORAD orbital parameters closest to |ULP00080 C : the 'IT' time included in the NORAD STORY data file. |ULP00090 C----------------------------------------------------------------------|ULP00100 C Author : Daniele Mortari |ULP00110 C----------------------------------------------------------------------|ULP00120 C Programmer : Daniele Mortari |ULP00130 C----------------------------------------------------------------------|ULP00140 C Date-coded : 23 April 1988 |ULP00150 C----------------------------------------------------------------------|ULP00160 C Last update : 02 December 1988 |ULP00170 C----------------------------------------------------------------------|ULP00180 C Notes : None |ULP00190 C----------------------------------------------------------------------|ULP00200 C Input data : |ULP00210 C IT = Reference time |ULP00220 C----------------------------------------------------------------------|ULP00230 C Output data : COMMON /ORB/ |ULP00240 C ITE = Epoch time |ULP00250 C EN0P2 = 1-st mean motion derivative |ULP00260 C EN0PP6 = 2-nd mean motion derivative |ULP00270 C BSTI = Drag term |ULP00280 C CLII = Inclination (Deg) |ULP00290 C RAI = Right ascension of the ascending node (Deg) |ULP00300 C ECI = Eccentricity |ULP00310 C API = Perigee argument (Deg) |ULP00320 C AMI = Mean anomaly (Deg) |ULP00330 C ENI = Mean motion (Rev/Day) |ULP00340 C----------------------------------------------------------------------+ULP00350 C ULP00360 IMPLICIT REAL*8(A-H,O-Z) ULP00370 INTEGER IT(6) ULP00380 COMMON /ORB/ ITE(6),EN0P2,EN0PP6,BSTI,CLII,RAI,ECI,API,AMI,ENI ULP00390 C ULP00400 ITE(1)=IT(1) ULP00410 ITE(2)=1 ULP00420 ITE(3)=1 ULP00430 ITE(4)=0 ULP00440 ITE(5)=0 ULP00450 ITE(6)=0 ULP00460 CALL T (ITE,FDGE) ULP00470 IY=ITE(1)-1900 ULP00480 C ULP00490 C Reference data ULP00500 C ULP00510 CALL T (IT,FDG) ULP00520 DDAY=FDG-FDGE+1.D0 ULP00530 C ULP00540 C Read in the file NORAD STORY (Filedef=99) ULP00550 C ULP00560 REWIND 99 ULP00570 READ(99,1000,END=9999) IA1,A1,B1,C1,IC1,D1,ID1,E1,F1,G1,H1,O1,P1 ULP00580 1000 FORMAT(18X,I2,F12.8,1X,F10.8,2(1X,F6.5,I2),/,7X,2(1X,F8.4),1X, ULP00590 $ F7.7,2(1X,F8.4),1X,F11.8) ULP00600 C ULP00610 IF(IY.LT.IA1) GO TO 2 ULP00620 IF(IY.GT.IA1.OR.DDAY.GE.A1) GO TO 1 ULP00630 C ULP00640 2 DATE=1000.D0*IA1+A1 ULP00650 CALL EPOCH (DATE,ITE,DITE) ULP00660 C ULP00670 EN0P2=B1 ULP00680 EN0PP6=C1*(10.D0**IC1) ULP00690 BSTI=D1*(10.D0**ID1) ULP00700 CLII=E1 ULP00710 RAI=F1 ULP00720 ECI=G1 ULP00730 API=H1 ULP00740 AMI=O1 ULP00750 ENI=P1 ULP00760 C ULP00770 RETURN ULP00800 C ULP00810 C Read the NORAD 2-Lines ULP00820 C ULP00830 1 READ(99,1000,END=2) IA2,A2,B2,C2,IC2,D2,ID2,E2,F2,G2,H2,O2,P2 ULP00840 C ULP00850 IF(IY.LT.IA2) GO TO 2 ULP00860 IF(IY.GT.IA2) GO TO 3 ULP00870 IF(DDAY.LT.A2) GO TO 2 ULP00880 C ULP00890 3 IA1=IA2 ULP00900 IC1=IC2 ULP00910 ID1=ID2 ULP00920 A1=A2 ULP00930 B1=B2 ULP00940 C1=C2 ULP00950 D1=D2 ULP00960 E1=E2 ULP00970 F1=F2 ULP00980 G1=G2 ULP00990 H1=H2 ULP01000 O1=O2 ULP01010 P1=P2 ULP01020 GO TO 1 ULP01030 C ULP01040 C NO DATA in the NORAD STORY file ULP01050 C ULP01060 9999 WRITE(3,1001) ULP01070 1001 FORMAT(/,' From ULPAOR : NO NORAD data from file NORAD STORY ???')ULP01080 STOP ULP01090 C ULP01100 END ULP01110 SUBROUTINE NOI (IT,X,Y,Z,R,XD,YD,ZD,V,CHI,FLAT,H,P,ENER,A,ECO,RA, NOI00010 $ RP,PSI,ANVE,E,AMO,ENO,PER,CLIO,FMU,AP,RAO,FI,IFLAG,IEPT)NOI00020 C NOI00030 C----------------------------------------------------------------------+NOI00040 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - 0010 |NOI00050 C----------------------------------------------------------------------|NOI00060 C Title : NOI |NOI00070 C----------------------------------------------------------------------|NOI00080 C Function : Norad Orbital Integrator (Obsculating elements) |NOI00090 C----------------------------------------------------------------------|NOI00100 C Author : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI00110 C----------------------------------------------------------------------|NOI00120 C Programmer : Mortari Daniele |NOI00130 C----------------------------------------------------------------------|NOI00140 C Date-coded : 11 October 1988 |NOI00150 C----------------------------------------------------------------------|NOI00160 C Last update : 16 December 1988 |NOI00170 C----------------------------------------------------------------------|NOI00180 C Notes : None |NOI00190 C----------------------------------------------------------------------|NOI00200 C Input data : |NOI00210 C IT = Time at which the orbit computations are whished |NOI00220 C IFLAG = Flag index (IFLAG must be = 1 for initializations) |NOI00230 C IEPT = Selects which orbit integrator is choosen |NOI00240 C (1=SGP, 2=SGP4, 3=SGP8, 4=SDP4, 5=SDP8) |NOI00250 C----------------------------------------------------------------------|NOI00260 C Input data : COMMON /ORB/ |NOI00270 C ITE = Epoch time vector (Y-M-D-H-M-S) |NOI00280 C EN0P2 = 1-st time derivative of the mean motion (Rev/Day**2) |NOI00290 C EN0PP6 = 2-nd time derivative of the mean motion (Rev/Day**3) |NOI00300 C BSTI = BSTAR drag term |NOI00310 C CLII = Orbital inclination (Deg) |NOI00320 C RAI = Right ascension of ascending node (Deg) |NOI00330 C ECI = Eccentricity |NOI00340 C API = Perigee argument (Deg) |NOI00350 C AMI = Mean anomaly (Deg) |NOI00360 C ENI = Mean motion (Rev/Day) |NOI00370 C----------------------------------------------------------------------|NOI00380 C Input data : COMMON /CONST/ |NOI00390 C PIG = Pi greek |NOI00400 C CW = Degrees/Radiant |NOI00410 C REQ = Earth equatorial radius (Km) |NOI00420 C----------------------------------------------------------------------|NOI00430 C Output data : COMMON /CM2/ |NOI00440 C OMPP = Perigee argument derivative (Deg/Day) |NOI00450 C OMGP = Right ascension of ascending node derivative (Deg/Day) |NOI00460 C----------------------------------------------------------------------|NOI00470 C Output data : |NOI00480 C X, Y, Z = Vernal cartesian components (Km) |NOI00490 C R = Zenith radius (Km) |NOI00500 C XD,YD,ZD = Vernal velocity components (Km/Sec) |NOI00510 C V = Velocity (Km/Sec) |NOI00520 C CHI = Longitude from vernal axis (Deg) |NOI00530 C FLAT = Latitude from equatorial plane (Deg) |NOI00540 C H = Altitude from earth surface (Km) |NOI00550 C P = Semilatum rectum (Km) |NOI00560 C ENER = Energy (Km**2/Sec**2) |NOI00570 C A = Semimajor axis (Km) |NOI00580 C ECO = Eccentricity |NOI00590 C RA = Apogee radius (Km) |NOI00600 C RP = Perigee radius (Km) |NOI00610 C PSI = Longitude from Greenwitch meridian (Deg) |NOI00620 C ANVE = True anomaly (Deg) |NOI00630 C E = Eccentric anomaly (Deg) |NOI00640 C AMO = Mean anomaly (Deg) |NOI00650 C ENO = Mean motion (Rev/Day) |NOI00660 C PER = Orbit period (Min) |NOI00670 C CLIO = Orbital inclination (Deg) |NOI00680 C FMU = Complessive anomaly (Deg) |NOI00690 C AP = Perigee argument (Deg) |NOI00700 C RAO = Right ascension of ascending node (Deg) |NOI00710 C FI = Longitude from the ascending node (Deg) |NOI00720 C----------------------------------------------------------------------+NOI00730 C NOI00740 IMPLICIT REAL*8(A-H,O-Z) NOI00750 INTEGER IT(6) NOI00760 COMMON /CONST/ PIG,CW,REQ NOI00770 COMMON /ORB/ ITE(6),EN0P2,EN0PP6,BSTI,CLII,RAI,ECI,API,AMI,ENI NOI00780 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI00790 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI00800 C NOI00810 IF(IFLAG.EQ.0) GO TO 2 NOI00820 C NOI00830 C Input check for selected case NOI00840 C NOI00850 IF(IEPT.LE.0.OR.IEPT.GE.6) GO TO 2222 NOI00860 C NOI00870 C General constants NOI00880 C NOI00890 TWOPI=2.D0*PIG NOI00900 ZERO=1.D-6 NOI00910 TOTHRD=2.D0/3.D0 NOI00920 QO=120.D0 NOI00930 SO=78.D0 NOI00940 XJ2=1.082616D-3 NOI00950 XJ3=-0.253881D-5 NOI00960 XJ4=-1.65597D-6 NOI00970 XKE=0.743669161D-1 NOI00980 AE=1.D0 NOI00990 C NOI01000 C AMU = Earth gravitational constant (Km**3/Sec**2) NOI01010 C NOI01020 AMU=3.9860064D+5 NOI01030 C NOI01040 C Select ephemeris type and output times NOI01050 C WGS-72 Physical and geopotential constants are CK2 and CK4 NOI01060 C NOI01070 CK2=0.5D0*XJ2*AE**2 NOI01080 CK4=-0.375D0*XJ4*AE**4 NOI01090 QOMS2T=((QO-SO)*AE/REQ)**4 NOI01100 S=AE*(1.D0+SO/REQ) NOI01110 C NOI01120 BSTAR=BSTI/AE NOI01130 ECC=ECI NOI01140 RAAN=RAI/CW NOI01150 OMPE=API/CW NOI01160 TETAE=AMI/CW NOI01170 CLIN=CLII/CW NOI01180 TEMP=TWOPI/1440.D0/1440.D0 NOI01190 XNO=ENI*TEMP*1440.D0 NOI01200 IF(XNO.LE.0.D0) GO TO 1111 NOI01210 XNDT2O=EN0P2*TEMP NOI01220 XNDD6O=EN0PP6*TEMP/1440.D0 NOI01230 C NOI01240 C Input check for period (Period GE 225 minutes is deep space) NOI01250 C NOI01260 IF(IEPT.LE.3.AND.(TWOPI/XNO/1440.D0).LT.(0.15625D0)) GO TO 1 NOI01270 IF(IEPT.GE.4.AND.(TWOPI/XNO/1440.D0).GE.(0.15625D0)) GO TO 1 NOI01280 IF(IEPT.LE.3) WRITE(3,1000) NOI01290 IF(IEPT.LE.3) WRITE(6,1000) NOI01300 1000 FORMAT(/,' Should use deep-space ephemeris',/) NOI01310 IF(IEPT.GE.4) WRITE(3,1001) NOI01320 IF(IEPT.GE.4) WRITE(6,1001) NOI01330 1001 FORMAT(/,' Should use near-earth ephemeris',/) NOI01340 STOP NOI01350 C NOI01360 1 A1=(XKE/XNO)**TOTHRD NOI01370 TEMP=1.5D0*CK2*(3.D0*DCOS(CLIN)**2-1.D0)/(1.D0-ECC**2)**1.5D0 NOI01380 DEL1=TEMP/(A1*A1) NOI01390 AO=A1*(1.D0-DEL1*(0.5D0*TOTHRD+DEL1*(1.D0+134.D0/81.D0*DEL1))) NOI01400 DELO=TEMP/(AO*AO) NOI01410 XNODP=XNO/(1.D0+DELO) NOI01420 CALL T (ITE,FDE) NOI01430 C NOI01440 2 CALL T (IT,FD) NOI01450 TSINCE=(FD-FDE)*1440.D0 NOI01460 C NOI01470 C Orbital integrators NOI01480 C NOI01490 IF(IEPT.EQ.1) CALL SGP (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI01500 IF(IEPT.EQ.2) CALL SGP4 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI01510 IF(IEPT.EQ.3) CALL SGP8 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI01520 IF(IEPT.EQ.4) CALL SDP4 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI01530 IF(IEPT.EQ.5) CALL SDP8 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI01540 C NOI01550 C Position NOI01560 C NOI01570 X=X*REQ/AE NOI01580 Y=Y*REQ/AE NOI01590 Z=Z*REQ/AE NOI01600 R=DSQRT(X**2+Y**2+Z**2) NOI01610 C NOI01620 C Velocity (Km/Sec) NOI01630 C NOI01640 XD=XD*REQ/AE/60.D0 NOI01650 YD=YD*REQ/AE/60.D0 NOI01660 ZD=ZD*REQ/AE/60.D0 NOI01670 V2=XD**2+YD**2+ZD**2 NOI01680 V=DSQRT(V2) NOI01690 C NOI01700 C CHI = Longitude from vernal axis (Deg) NOI01710 C NOI01720 SCHI=Y/DSQRT(X**2+Y**2) NOI01730 CCHI=X/DSQRT(X**2+Y**2) NOI01740 CALL FTAN (SCHI,CCHI,CHI) NOI01750 C NOI01760 C FLAT = Latitude (Deg) NOI01770 C NOI01780 SLAT=Z/R NOI01790 FLAT=CW*DASIN(SLAT) NOI01800 CLAT=DCOS(FLAT/CW) NOI01810 C NOI01820 C H = Altitude from earth surface (Km) NOI01830 C NOI01840 H=R-REQ NOI01850 IF(H.GT.100.D0) GO TO 3 NOI01860 WRITE(3,1002) NOI01870 WRITE(6,1002) NOI01880 1002 FORMAT(//,' From NOI: Satellite altitude less than 100 Km |||',//)NOI01890 STOP NOI01900 C NOI01910 C P = Semilatum rectum (Km) NOI01920 C NOI01930 3 W2=(Y*ZD-Z*YD)**2+(Z*XD-X*ZD)**2+(X*YD-Y*XD)**2 NOI01940 P=W2/AMU NOI01950 C NOI01960 C ENER = Energy (Km**2/Sec**2) NOI01970 C NOI01980 ENER=V2/2.D0-AMU/R NOI01990 C NOI02000 C A = Semi-major axis (Km) NOI02010 C NOI02020 A=-AMU/(2.D0*ENER) NOI02030 C NOI02040 C ECC = Eccentricity NOI02050 C NOI02060 ECO=DSQRT(1.D0-P/A) NOI02070 C NOI02080 C RA = Apogee radius (Km) NOI02090 C NOI02100 RA=A*(1.D0+ECO) NOI02110 C NOI02120 C RP = Perigee radius (Km) NOI02130 C NOI02140 RP=A*(1.D0-ECO) NOI02150 C NOI02160 C PSI = Longitude from Greenwich (Deg) NOI02170 C NOI02180 CALL SUN (IT,CHIS,DECL) NOI02190 PSI=PRIMO(CHI-CHIS+15.D0*(12.D0-IT(4)-(IT(5)+IT(6)/60.D0)/60.D0)) NOI02200 C NOI02210 C ANVE = True anomaly (Deg) NOI02220 C NOI02230 VRAD=ZD*SLAT+(XD*CCHI+YD*SCHI)*CLAT NOI02240 SANVE=DSQRT(W2)*VRAD/(ECO*AMU) NOI02250 CANVE=(P/R-1.D0)/ECO NOI02260 CALL FTAN (SANVE,CANVE,ANVE) NOI02270 C NOI02280 C E = Eccentric anomaly (Deg) NOI02290 C NOI02300 SE=SANVE*DSQRT(1.D0-ECO**2)/(1.D0+ECO*CANVE) NOI02310 CE=(CANVE+ECO)/(1.D0+ECO*CANVE) NOI02320 CALL FTAN (SE,CE,E) NOI02330 C NOI02340 C EM = Mean anomaly (Deg) NOI02350 C NOI02360 AMO=E-ECO*SE*CW NOI02370 C NOI02380 C EN = Mean motion (Rev/Day) NOI02390 C NOI02400 ENO=86400.D0*DSQRT(AMU/A)/A/TWOPI NOI02410 C NOI02420 C PER = Orbital period (Min) NOI02430 C NOI02440 PER=1440.D0/ENO NOI02450 C NOI02460 C CLIO = Inclination (Deg) NOI02470 C NOI02480 VPSI=ZD*CLAT-(XD*CCHI+YD*SCHI)*SLAT NOI02490 VLAM=YD*CCHI-XD*SCHI NOI02500 VTET=DSQRT(VLAM**2+VPSI**2) NOI02510 SANG=VPSI/VTET NOI02520 CANG=VLAM/VTET NOI02530 CALL FTAN (SANG,CANG,ANG) NOI02540 CLIO=CW*DACOS(CLAT*CANG) NOI02550 C NOI02560 C FMU = Complessive anomaly (Deg) NOI02570 C NOI02580 SFMU=SLAT/DSIN(CLIO/CW) NOI02590 CFMU=DSQRT(1.D0-SFMU**2) NOI02600 IF(SANG.LE.0.D0) CFMU=-CFMU NOI02610 CALL FTAN (SFMU,CFMU,FMU) NOI02620 C NOI02630 C AP = Perigee argument (Deg) NOI02640 C NOI02650 AP=PRIMO(FMU-ANVE) NOI02660 C NOI02670 C RAO = Right ascension of ascending node (Deg) NOI02680 C NOI02690 CANG=CFMU/CLAT NOI02700 SANG=CANG*DCOS(CLIO/CW)*SFMU/CFMU NOI02710 CALL FTAN (SANG,CANG,ANG) NOI02720 RAO=PRIMO(CHI-ANG) NOI02730 C NOI02740 C FI = Longitude from the ascending node (Deg) NOI02750 C NOI02760 FI=PRIMO(CHI-RAO) NOI02770 RETURN NOI02780 C NOI02790 C Wrong returns NOI02800 C NOI02810 1111 WRITE(3,1003) XNO NOI02820 WRITE(6,1003) XNO NOI02830 1003 FORMAT(/,'From NOI : Mean motion = ',D15.7,' < 0.D0 (???)',/) NOI02840 STOP NOI02850 C NOI02860 2222 WRITE(3,1004) IEPT NOI02870 WRITE(6,1004) IEPT NOI02880 1004 FORMAT(//,'Ephemeris number ',I2,' NOT legal (IEPT=1,2,3,4,5)',//)NOI02890 STOP NOI02900 C NOI02910 END NOI02920 SUBROUTINE SGP (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI02930 C NOI02940 C----------------------------------------------------------------------+NOI02950 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI02960 C----------------------------------------------------------------------|NOI02970 C Title : SGP |NOI02980 C----------------------------------------------------------------------|NOI02990 C Function : Orbital integrator |NOI03000 C----------------------------------------------------------------------|NOI03010 C Author : NORAD |NOI03020 C----------------------------------------------------------------------|NOI03030 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI03040 C----------------------------------------------------------------------|NOI03050 C Date-coded : 11 October 1988 |NOI03060 C----------------------------------------------------------------------|NOI03070 C Last update : 14 December 1988 |NOI03080 C----------------------------------------------------------------------|NOI03090 C Notes : Notes |NOI03100 C----------------------------------------------------------------------|NOI03110 C Input data : |NOI03120 C----------------------------------------------------------------------|NOI03130 C Output data : |NOI03140 C----------------------------------------------------------------------+NOI03150 C NOI03160 IMPLICIT REAL*8(A-H,O-Z) NOI03170 COMMON /CONST/ DUM1,CW,DUM2 NOI03180 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI03190 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI03200 COMMON /CM2/ OMPP,OMGP NOI03210 C NOI03220 IF(IFLAG.EQ.0) GO TO 1 NOI03230 C NOI03240 C Initialization NOI03250 C NOI03260 C1=CK2*1.5D0 NOI03270 C2=CK2/4.D0 NOI03280 C3=CK2/2.D0 NOI03290 C4=XJ3*AE**3/(4.D0*CK2) NOI03300 COSIO=DCOS(CLIN) NOI03310 SINIO=DSIN(CLIN) NOI03320 A1=(XKE/XNO)**TOTHRD NOI03330 D1=C1/A1/A1*(3.D0*COSIO**2-1.D0)/(1.D0-ECC**2)**1.5D0 NOI03340 AO=A1*(1.D0-1.D0/3.D0*D1-D1**2-134.D0/81.D0*D1**3) NOI03350 PO=AO*(1.D0-ECC**2) NOI03360 QO=AO*(1.D0-ECC) NOI03370 XLO=TETAE+OMPE+RAAN NOI03380 D10=C3*SINIO**2 NOI03390 D20=C2*(7.D0*COSIO**2-1.D0) NOI03400 D30=C1*COSIO NOI03410 D40=D30*SINIO NOI03420 PO2NO=XNO/PO**2 NOI03430 OMGDT=C1*PO2NO*(5.D0*COSIO**2-1.D0) NOI03440 OMPP= OMGDT*1440.D0*CW NOI03450 XNODOT=-2.D0*D30*PO2NO NOI03460 OMGP=XNODOT*1440.D0*CW NOI03470 C5=0.5D0*C4*SINIO*(3.D0+5.D0*COSIO)/(1.D0+COSIO) NOI03480 C6=C4*SINIO NOI03490 IFLAG=0 NOI03500 C NOI03510 C Update for secular gravity and atmospheric drag NOI03520 C NOI03530 1 A=XNO+(2.D0*XNDT2O+3.D0*XNDD6O*TSINCE)*TSINCE NOI03540 A=AO*(XNO/A)**TOTHRD NOI03550 E=ZERO NOI03560 IF(A.GT.QO) E=1.D0-QO/A NOI03570 P=A*(1.D0-E**2) NOI03580 XNODES=RAAN+XNODOT*TSINCE NOI03590 OMGAS=OMPE+OMGDT*TSINCE NOI03600 XLS=FMOD2P(XLO+(XNO+OMGDT+XNODOT+(XNDT2O+XNDD6O*TSINCE)* NOI03610 $ TSINCE)*TSINCE) NOI03620 C NOI03630 C Long period periodics NOI03640 C NOI03650 AXNSL=E*DCOS(OMGAS) NOI03660 AYNSL=E*DSIN(OMGAS)-C6/P NOI03670 XL=FMOD2P(XLS-C5/P*AXNSL) NOI03680 C NOI03690 C Solve Keplers equation NOI03700 C NOI03710 U=FMOD2P(XL-XNODES) NOI03720 ITEM3=0 NOI03730 EO1=U NOI03740 TEM5=1.D0 NOI03750 C NOI03760 2 SINEO1=DSIN(EO1) NOI03770 COSEO1=DCOS(EO1) NOI03780 IF(DABS(TEM5).LT.ZERO.OR.ITEM3.GE.10) GO TO 3 NOI03790 ITEM3=ITEM3+1 NOI03800 TEM5=1.D0-COSEO1*AXNSL-SINEO1*AYNSL NOI03810 TEM5=(U-AYNSL*COSEO1+AXNSL*SINEO1-EO1)/TEM5 NOI03820 TEM2=DABS(TEM5) NOI03830 IF(TEM2.GT.1) TEM5=TEM2/TEM5 NOI03840 EO1=EO1+TEM5 NOI03850 GO TO 2 NOI03860 C NOI03870 C Short period preliminary quantities NOI03880 C NOI03890 3 ECOSE=AXNSL*COSEO1+AYNSL*SINEO1 NOI03900 ESINE=AXNSL*SINEO1-AYNSL*COSEO1 NOI03910 EL2=AXNSL**2+AYNSL**2 NOI03920 PL=A*(1.D0-EL2) NOI03930 PL2=PL**2 NOI03940 R=A*(1.D0-ECOSE) NOI03950 RDOT=XKE*DSQRT(A)/R*ESINE NOI03960 RVDOT=XKE*DSQRT(PL)/R NOI03970 TEMP=ESINE/(1.D0+DSQRT(1.D0-EL2)) NOI03980 SINU=A/R*(SINEO1-AYNSL-AXNSL*TEMP) NOI03990 COSU=A/R*(COSEO1-AXNSL+AYNSL*TEMP) NOI04000 CALL FTAN (SINU,COSU,SU) NOI04010 SU=SU/CW NOI04020 SIN2U=2.D0*COSU*SINU NOI04030 COS2U=1.D0-2.D0*SINU**2 NOI04040 C NOI04050 C Update for short periodics NOI04060 C NOI04070 RK=R+D10/PL*COS2U NOI04080 UK=SU-D20/PL2*SIN2U NOI04090 XNODEK=XNODES+D30*SIN2U/PL2 NOI04100 XINCK=CLIN+D40/PL2*COS2U NOI04110 C NOI04120 C Orientation vectors NOI04130 C NOI04140 SINUK=DSIN(UK) NOI04150 COSUK=DCOS(UK) NOI04160 SINNOK=DSIN(XNODEK) NOI04170 COSNOK=DCOS(XNODEK) NOI04180 SINIK=DSIN(XINCK) NOI04190 COSIK=DCOS(XINCK) NOI04200 C NOI04210 XMX=-SINNOK*COSIK NOI04220 XMY=COSNOK*COSIK NOI04230 C NOI04240 UX=XMX*SINUK+COSNOK*COSUK NOI04250 UY=XMY*SINUK+SINNOK*COSUK NOI04260 UZ=SINIK*SINUK NOI04270 C NOI04280 VX=XMX*COSUK-COSNOK*SINUK NOI04290 VY=XMY*COSUK-SINNOK*SINUK NOI04300 VZ=SINIK*COSUK NOI04310 C NOI04320 C Position and velocity NOI04330 C NOI04340 X=RK*UX NOI04350 Y=RK*UY NOI04360 Z=RK*UZ NOI04370 C NOI04380 XD=RVDOT*VX+RDOT*UX NOI04390 YD=RVDOT*VY+RDOT*UY NOI04400 ZD=RVDOT*VZ+RDOT*UZ NOI04410 RETURN NOI04420 C NOI04430 END NOI04440 SUBROUTINE SGP4 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI04450 C NOI04460 C----------------------------------------------------------------------+NOI04470 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI04480 C----------------------------------------------------------------------|NOI04490 C Title : SGP4 |NOI04500 C----------------------------------------------------------------------|NOI04510 C Function : Orbital integrator |NOI04520 C----------------------------------------------------------------------|NOI04530 C Author : NORAD |NOI04540 C----------------------------------------------------------------------|NOI04550 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI04560 C----------------------------------------------------------------------|NOI04570 C Date-coded : 11 October 1988 |NOI04580 C----------------------------------------------------------------------|NOI04590 C Last update : 14 December 1988 |NOI04600 C----------------------------------------------------------------------|NOI04610 C Notes : None |NOI04620 C----------------------------------------------------------------------|NOI04630 C Input data : |NOI04640 C----------------------------------------------------------------------|NOI04650 C Output data : |NOI04660 C----------------------------------------------------------------------+NOI04670 C NOI04680 IMPLICIT REAL*8(A-H,O-Z) NOI04690 COMMON /CONST/ DUM1,CW,REQ NOI04700 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI04710 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI04720 COMMON /CM2/ OMPP,OMGP NOI04730 C NOI04740 IF(IFLAG.EQ.0) GO TO 2 NOI04750 IFLAG=0 NOI04760 C NOI04770 C Recover original mean motion (XNODP) and semimajor axis (AODP) NOI04780 C from input elements NOI04790 C NOI04800 A1=(XKE/XNO)**TOTHRD NOI04810 COSIO=DCOS(CLIN) NOI04820 THETA2=COSIO**2 NOI04830 X3THM1=3.D0*THETA2-1.D0 NOI04840 ECCSQ=ECC**2 NOI04850 BETAO2=1.D0-ECCSQ NOI04860 BETAO=DSQRT(BETAO2) NOI04870 DEL1=1.5D0*CK2*X3THM1/(A1**2*BETAO*BETAO2) NOI04880 AO=A1*(1.D0-DEL1*(0.5D0*TOTHRD+DEL1*(1.D0+134.D0/81.D0*DEL1))) NOI04890 DELO=1.5D0*CK2*X3THM1/(AO**2*BETAO*BETAO2) NOI04900 XNODP=XNO/(1.D0+DELO) NOI04910 AODP=AO/(1.D0-DELO) NOI04920 C ---> INITIALIZATIONNOI04930 C NOI04940 C For perigee less than 220 Km, the ISIMP flag is set and the equations NOI04950 C are truncated to linear variation in DSQRT and quadratic variation NOI04960 C in mean anomaly. Also, the C3 term, the delta omega term, and the NOI04970 C DELTAM term are dropped. NOI04980 C NOI04990 ISIMP=0 NOI05000 IF((AODP*(1.D0-ECC)/AE).LT.(220.D0/REQ+AE)) ISIMP=1 NOI05010 C NOI05020 C For perigee below 156 Km, the values of S and QOMS2T are altered. NOI05030 C NOI05040 S4=S NOI05050 QOMS24=QOMS2T NOI05060 PERIGE=(AODP*(1.D0-ECC)-AE)*REQ NOI05070 IF(PERIGE.GE.(156.D0)) GO TO 1 NOI05080 S4=PERIGE-78.D0 NOI05090 IF(PERIGE.LE.(98.D0)) S4=20.D0 NOI05100 QOMS24=((120.D0-S4)*AE/REQ)**4 NOI05110 S4=S4/REQ+AE NOI05120 C NOI05130 1 PINVSQ=1.D0/(AODP**2*BETAO2**2) NOI05140 TSI=1.D0/(AODP-S4) NOI05150 ETA=AODP*ECC*TSI NOI05160 ETASQ=ETA**2 NOI05170 EETA=ECC*ETA NOI05180 PSISQ=DABS(1.D0-ETASQ) NOI05190 COEF=QOMS24*TSI**4 NOI05200 COEF1=COEF/PSISQ**3.5D0 NOI05210 C2=COEF1*XNODP*(AODP*(1.D0+1.5D0*ETASQ+EETA*(4.D0+ETASQ)) NOI05220 $ +0.75D0*CK2*TSI/PSISQ*X3THM1*(8.D0+3.D0*ETASQ*(8.D0+ETASQ))) NOI05230 C1=BSTAR*C2 NOI05240 SINIO=DSIN(CLIN) NOI05250 A3OVK2=-XJ3/CK2*AE**3 NOI05260 C3=COEF*TSI*A3OVK2*XNODP*AE*SINIO/ECC NOI05270 X1MTH2=1.D0-THETA2 NOI05280 C4=2.D0*XNODP*COEF1*AODP*BETAO2*(ETA*(2.D0+0.5D0*ETASQ) NOI05290 $ +ECC*(0.5D0+2.D0*ETASQ)-2.D0*CK2*TSI/(AODP*PSISQ)*( NOI05300 $ -3.D0*X3THM1*(1.D0-2.D0*EETA+ETASQ*(1.5D0-0.5D0*EETA)) NOI05310 $ +0.75D0*X1MTH2*(2.D0*ETASQ-EETA*(1.D0+ETASQ))*DCOS(2.D0*OMPE))) NOI05320 C5=2.D0*COEF1*AODP*BETAO2*(1.D0+2.75D0*(ETASQ+EETA)+EETA*ETASQ) NOI05330 THETA4=THETA2**2 NOI05340 TEMP1=3.D0*CK2*PINVSQ*XNODP NOI05350 TEMP2=TEMP1*CK2*PINVSQ NOI05360 TEMP3=1.25D0*CK4*PINVSQ**2*XNODP NOI05370 XMDOT=XNODP+0.5D0*TEMP1*BETAO*X3THM1+0.0625D0*TEMP2*BETAO* NOI05380 $ (13.D0-78.D0*THETA2+137.D0*THETA4) NOI05390 X1M5TH=1.D0-5.D0*THETA2 NOI05400 OMGDOT=-0.5D0*TEMP1*X1M5TH+0.0625D0*TEMP2*(7.D0-114.D0*THETA2+ NOI05410 $ 395.D0*THETA4)+TEMP3*(3.D0-36.D0*THETA2+49.D0*THETA4) NOI05420 OMPP=OMGDOT*1440.D0*CW NOI05430 XHDOT1=-TEMP1*COSIO NOI05440 XNODOT=XHDOT1+(0.5D0*TEMP2*(4.D0-19.D0*THETA2)+2.D0*TEMP3*(3.D0 NOI05450 $ -7.D0*THETA2))*COSIO NOI05460 OMGP=XNODOT*1440.D0*CW NOI05470 OMGCOF=BSTAR*C3*DCOS(OMPE) NOI05480 XMCOF=-TOTHRD*COEF*BSTAR*AE/EETA NOI05490 XNODCF=3.5D0*BETAO2*XHDOT1*C1 NOI05500 T2COF=1.5D0*C1 NOI05510 XLCOF=0.125D0*A3OVK2*SINIO*(3.D0+5.D0*COSIO)/(1.D0+COSIO) NOI05520 AYCOF=0.25D0*A3OVK2*SINIO NOI05530 DELMO=(1.D0+ETA*DCOS(TETAE))**3 NOI05540 SINMO=DSIN(TETAE) NOI05550 X7THM1=7.D0*THETA2-1.D0 NOI05560 IF(ISIMP.EQ.1) GO TO 2 NOI05570 C1SQ=C1**2 NOI05580 D2=4.D0*AODP*TSI*C1SQ NOI05590 TEMP=D2*TSI*C1/3.D0 NOI05600 D3=(17.D0*AODP+S4)*TEMP NOI05610 D4=0.5D0*TEMP*AODP*TSI*(221.D0*AODP+31.D0*S4)*C1 NOI05620 T3COF=D2+2.D0*C1SQ NOI05630 T4COF=0.25D0*(3.D0*D3+C1*(12.D0*D2+10.D0*C1SQ)) NOI05640 T5COF=0.2D0*(3.D0*D4+12.D0*C1*D3+6.D0*D2**2+15.D0*C1SQ*( NOI05650 $ 2.D0*D2+C1SQ)) NOI05660 C NOI05670 C Update for secular gravity and atmospheric drag NOI05680 C NOI05690 2 XMDF=TETAE+XMDOT*TSINCE NOI05700 OMGADF=OMPE+OMGDOT*TSINCE NOI05710 XNODDF=RAAN+XNODOT*TSINCE NOI05720 OMEGA=OMGADF NOI05730 XMP=XMDF NOI05740 TSQ=TSINCE**2 NOI05750 XNODE=XNODDF+XNODCF*TSQ NOI05760 TEMPA=1.D0-C1*TSINCE NOI05770 TEMPE=BSTAR*C4*TSINCE NOI05780 TEMPL=T2COF*TSQ NOI05790 IF(ISIMP.EQ.1) GO TO 3 NOI05800 DELOMG=OMGCOF*TSINCE NOI05810 DELM=XMCOF*((1.D0+ETA*DCOS(XMDF))**3-DELMO) NOI05820 TEMP=DELOMG+DELM NOI05830 XMP=XMDF+TEMP NOI05840 OMEGA=OMGADF-TEMP NOI05850 TCUBE=TSQ*TSINCE NOI05860 TFOUR=TSINCE*TCUBE NOI05870 TEMPA=TEMPA-D2*TSQ-D3*TCUBE-D4*TFOUR NOI05880 TEMPE=TEMPE+BSTAR*C5*(DSIN(XMP)-SINMO) NOI05890 TEMPL=TEMPL+T3COF*TCUBE+TFOUR*(T4COF+TSINCE*T5COF) NOI05900 C NOI05910 3 A=AODP*TEMPA**2 NOI05920 E=ECC-TEMPE NOI05930 XL=XMP+OMEGA+XNODE+XNODP*TEMPL NOI05940 BETA=DSQRT(1.D0-E**2) NOI05950 XN=XKE/A**1.5D0 NOI05960 C NOI05970 C Long period periodics NOI05980 C NOI05990 AXN=E*DCOS(OMEGA) NOI06000 TEMP=1.D0/(A*BETA**2) NOI06010 XLL=TEMP*XLCOF*AXN NOI06020 AYNL=TEMP*AYCOF NOI06030 XLT=XL+XLL NOI06040 AYN=E*DSIN(OMEGA)+AYNL NOI06050 C NOI06060 C Solve Keplers equation NOI06070 C NOI06080 CAPU=FMOD2P(XLT-XNODE) NOI06090 TEMP2=CAPU NOI06100 DO 4 I=1,10 NOI06110 SINEPW=DSIN(TEMP2) NOI06120 COSEPW=DCOS(TEMP2) NOI06130 TEMP3=AXN*SINEPW NOI06140 TEMP4=AYN*COSEPW NOI06150 TEMP5=AXN*COSEPW NOI06160 TEMP6=AYN*SINEPW NOI06170 EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.D0-TEMP5-TEMP6)+TEMP2 NOI06180 IF(DABS(EPW-TEMP2).LE.ZERO) GO TO 5 NOI06190 4 TEMP2=EPW NOI06200 C NOI06210 C Short period preliminary quantities NOI06220 C NOI06230 5 ECOSE=TEMP5+TEMP6 NOI06240 ESINE=TEMP3-TEMP4 NOI06250 ELSQ=AXN**2+AYN**2 NOI06260 TEMP=1.D0-ELSQ NOI06270 PL=A*TEMP NOI06280 R=A*(1.D0-ECOSE) NOI06290 TEMP1=1.D0/R NOI06300 RDOT=XKE*DSQRT(A)*ESINE*TEMP1 NOI06310 RFDOT=XKE*DSQRT(PL)*TEMP1 NOI06320 TEMP2=A*TEMP1 NOI06330 BETAL=DSQRT(TEMP) NOI06340 TEMP3=1.D0/(1.D0+BETAL) NOI06350 COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3) NOI06360 SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3) NOI06370 CALL FTAN (SINU,COSU,U) NOI06380 U=U/CW NOI06390 SIN2U=2.D0*SINU*COSU NOI06400 COS2U=2.D0*COSU**2-1.D0 NOI06410 TEMP=1.D0/PL NOI06420 TEMP1=CK2*TEMP NOI06430 TEMP2=TEMP1*TEMP NOI06440 C NOI06450 C Update for short periodics NOI06460 C NOI06470 RK=R*(1.D0-1.5D0*TEMP2*BETAL*X3THM1)+0.5D0*TEMP1*X1MTH2*COS2U NOI06480 UK=U-0.25D0*TEMP2*X7THM1*SIN2U NOI06490 XNODEK=XNODE+1.5D0*TEMP2*COSIO*SIN2U NOI06500 XINCK=CLIN+1.5D0*TEMP2*COSIO*SINIO*COS2U NOI06510 RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U NOI06520 RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5D0*X3THM1) NOI06530 C NOI06540 C Orientation vectors NOI06550 C NOI06560 SINUK=DSIN(UK) NOI06570 COSUK=DCOS(UK) NOI06580 SINIK=DSIN(XINCK) NOI06590 COSIK=DCOS(XINCK) NOI06600 SINNOK=DSIN(XNODEK) NOI06610 COSNOK=DCOS(XNODEK) NOI06620 XMX=-SINNOK*COSIK NOI06630 XMY=COSNOK*COSIK NOI06640 UX=XMX*SINUK+COSNOK*COSUK NOI06650 UY=XMY*SINUK+SINNOK*COSUK NOI06660 UZ=SINIK*SINUK NOI06670 VX=XMX*COSUK-COSNOK*SINUK NOI06680 VY=XMY*COSUK-SINNOK*SINUK NOI06690 VZ=SINIK*COSUK NOI06700 C NOI06710 C Position and velocity NOI06720 C NOI06730 X=RK*UX NOI06740 Y=RK*UY NOI06750 Z=RK*UZ NOI06760 C NOI06770 XD=RDOTK*UX+RFDOTK*VX NOI06780 YD=RDOTK*UY+RFDOTK*VY NOI06790 ZD=RDOTK*UZ+RFDOTK*VZ NOI06800 RETURN NOI06810 C NOI06820 END NOI06830 SUBROUTINE SGP8 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI06840 C NOI06850 C----------------------------------------------------------------------+NOI06860 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI06870 C----------------------------------------------------------------------|NOI06880 C Title : SGP8 |NOI06890 C----------------------------------------------------------------------|NOI06900 C Function : Orbital integrator |NOI06910 C----------------------------------------------------------------------|NOI06920 C Author : NORAD |NOI06930 C----------------------------------------------------------------------|NOI06940 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI06950 C----------------------------------------------------------------------|NOI06960 C Date-coded : 11 October 1988 |NOI06970 C----------------------------------------------------------------------|NOI06980 C Last update : 14 December 1988 |NOI06990 C----------------------------------------------------------------------|NOI07000 C Notes : None |NOI07010 C----------------------------------------------------------------------|NOI07020 C Input data : |NOI07030 C----------------------------------------------------------------------|NOI07040 C Output data : |NOI07050 C----------------------------------------------------------------------+NOI07060 C NOI07070 IMPLICIT REAL*8(A-H,O-Z) NOI07080 COMMON /CONST/ DUM1,CW,DUM2 NOI07090 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI07100 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI07110 COMMON /CM2/ OMPP,OMGP NOI07120 DATA RHO/.15696615D0/ NOI07130 C NOI07140 IF(IFLAG.EQ.0) GO TO 2 NOI07150 IFLAG=0 NOI07160 C NOI07170 C Recover original mean motion (XNODP) and semimajor axis (AODP) from NOI07180 C input elements ------ calculate ballistic coefficient (B term) from NOI07190 C input B* drag term NOI07200 C NOI07210 A1=(XKE/XNO)**TOTHRD NOI07220 COSI=DCOS(CLIN) NOI07230 THETA2=COSI**2 NOI07240 TTHMUN=3.D0*THETA2-1.D0 NOI07250 ECCSQ=ECC**2 NOI07260 BETAO2=1.D0-ECCSQ NOI07270 BETAO=DSQRT(BETAO2) NOI07280 DEL1=1.5D0*CK2*TTHMUN/(A1**2*BETAO*BETAO2) NOI07290 AO=A1*(1.D0-DEL1*(0.5D0*TOTHRD+DEL1*(1.D0+134.D0/81.D0*DEL1))) NOI07300 DELO=1.5*CK2*TTHMUN/(AO**2*BETAO*BETAO2) NOI07310 AODP=AO/(1.D0-DELO) NOI07320 XNODP=XNO/(1.D0+DELO) NOI07330 B=2.D0*BSTAR/RHO NOI07340 C NOI07350 C Initialization NOI07360 C NOI07370 ISIMP=0 NOI07380 PO=AODP*BETAO2 NOI07390 POM2=1.D0/PO**2 NOI07400 SINI=DSIN(CLIN) NOI07410 COSI=DCOS(CLIN) NOI07420 SING=DSIN(OMPE) NOI07430 COSG=DCOS(OMPE) NOI07440 TEMP=0.5D0*CLIN NOI07450 SINIO2=DSIN(TEMP) NOI07460 COSIO2=DCOS(TEMP) NOI07470 THETA4=THETA2**2 NOI07480 UNM5TH=1.D0-5.D0*THETA2 NOI07490 UNMTH2=1.D0-THETA2 NOI07500 A3COF=-XJ3/CK2*AE**3 NOI07510 PARDT1=3.D0*CK2*POM2*XNODP NOI07520 PARDT2=PARDT1*CK2*POM2 NOI07530 PARDT4=1.25D0*CK4*POM2**2*XNODP NOI07540 XMDT1=0.5D0*PARDT1*BETAO*TTHMUN NOI07550 XGDT1=-0.5D0*PARDT1*UNM5TH NOI07560 XHDT1=-PARDT1*COSI NOI07570 XLLDOT=XNODP+XMDT1+0.0625D0*PARDT2*BETAO*(13.D0-78.D0*THETA2 NOI07580 $ +137.D0*THETA4) NOI07590 OMGDT=XGDT1+0.0625D0*PARDT2*(7.D0-114.D0*THETA2+395.D0*THETA4) NOI07600 $ +PARDT4*(3.D0-36.D0*THETA2+49.D0*THETA4) NOI07610 OMPP= OMGDT*1440.D0*CW NOI07620 XNODOT=XHDT1+(0.5D0*PARDT2*(4.D0-19.D0*THETA2)+2.D0*PARDT4*(3.D0 NOI07630 $ -7.D0*THETA2))*COSI NOI07640 OMGP=XNODOT*1440.D0*CW NOI07650 TSI=1.D0/(PO-S) NOI07660 ETA=ECC*S*TSI NOI07670 ETA2=ETA**2 NOI07680 PSIM2=DABS(1.D0/(1.D0-ETA2)) NOI07690 ALPHA2=1.D0+ECCSQ NOI07700 EETA=ECC*ETA NOI07710 COS2G=2.D0*COSG**2-1.D0 NOI07720 D5=TSI*PSIM2 NOI07730 D1=D5/PO NOI07740 D2=12.D0+ETA2*(36.D0+4.5D0*ETA2) NOI07750 D3=ETA2*(15.D0+2.5D0*ETA2) NOI07760 D4=ETA*(5.D0+3.75D0*ETA2) NOI07770 B1=CK2*TTHMUN NOI07780 B2=-CK2*UNMTH2 NOI07790 B3=A3COF*SINI NOI07800 C0=0.5D0*B*RHO*QOMS2T*XNODP*AODP*TSI**4*PSIM2**3.5D0/DSQRT(ALPHA2)NOI07810 C1=1.5D0*XNODP*ALPHA2**2*C0 NOI07820 C4=D1*D3*B2 NOI07830 C5=D5*D4*B3 NOI07840 XNDT=C1*((2.D0+ETA2*(3.D0+34.D0*ECCSQ)+5.D0*EETA*(4.D0+ETA2) NOI07850 $ +8.5D0*ECCSQ)+D1*D2*B1+C4*COS2G+C5*SING) NOI07860 XNDTN=XNDT/XNODP NOI07870 C NOI07880 C If drag is very small, the ISIMP flag is set and the equations are NOI07890 C truncated to linear variation in mean motion and quadratic variationNOI07900 C in mean anomaly NOI07910 C NOI07920 IF(DABS(XNDTN*1440.D0).LT.(2.16D-3)) GO TO 1 NOI07930 D6=ETA*(30.D0+22.5D0*ETA2) NOI07940 D7=ETA*(5.D0+12.5D0*ETA2) NOI07950 D8=1.D0+ETA2*(6.75D0+ETA2) NOI07960 C8=D1*D7*B2 NOI07970 C9=D5*D8*B3 NOI07980 EDOT=-C0*(ETA*(4.D0+ETA2+ECCSQ*(15.5D0+7.D0*ETA2))+ECC*(5.D0 NOI07990 $ +15.D0*ETA2)+D1*D6*B1+C8*COS2G+C9*SING) NOI08000 D2O=0.5D0*TOTHRD*XNDTN NOI08010 ALDTAL=ECC*EDOT/ALPHA2 NOI08020 TSDTTS=2.D0*AODP*TSI*(D2O*BETAO2+ECC*EDOT) NOI08030 ETDT=(EDOT+ECC*TSDTTS)*TSI*S NOI08040 PSDTPS=-ETA*ETDT*PSIM2 NOI08050 SIN2G=2.D0*SING*COSG NOI08060 CODTCO=D2O+4.D0*TSDTTS-ALDTAL-7.D0*PSDTPS NOI08070 C1DTC1=XNDTN+4.D0*ALDTAL+CODTCO NOI08080 D9=ETA*(6.D0+68.D0*ECCSQ)+ECC*(20.D0+15.D0*ETA2) NOI08090 D10=5.D0*ETA*(4.D0+ETA2)+ECC*(17.D0+68.D0*ETA2) NOI08100 D11=ETA*(72.D0+18.D0*ETA2) NOI08110 D12=ETA*(30.D0+10.D0*ETA2) NOI08120 D13=5.D0+11.25*ETA2 NOI08130 D14=TSDTTS-2.D0*PSDTPS NOI08140 D15=2.D0*(D2O+ECC*EDOT/BETAO2) NOI08150 D1DT=D1*(D14+D15) NOI08160 D2DT=ETDT*D11 NOI08170 D3DT=ETDT*D12 NOI08180 D4DT=ETDT*D13 NOI08190 D5DT=D5*D14 NOI08200 C4DT=B2*(D1DT*D3+D1*D3DT) NOI08210 C5DT=B2*(D5DT*D4+D5*D4DT) NOI08220 D16=D9*ETDT+D10*EDOT+B1*(D1DT*D2+D1*D2DT)+ NOI08230 $ C4DT*COS2G+C5DT*SING+XGDT1*(C5*COSG-2.D0*C4*SIN2G) NOI08240 XNDDT=C1DTC1*XNDT+C1*D16 NOI08250 EDDOT=C0DTC0*EDOT-C0*((4.D0+3.D0*ETA2+30.D0*EETA+ECCSQ*(15.5D0 NOI08260 $ +21.D0*ETA2))*ETDT+(5.D0+15.D0*ETA2+EETA*(31.D0 NOI08270 $ +14.D0*ETA2))*EDOT+B1*(D1DT*D6+D1*ETDT*(30.D0+67.5D0*ETA2)) NOI08280 $ +B2*(D1DT*D7+D1*ETDT*(5.D0+37.5D0*ETA2))*COS2O+B3*(D5DT*D8 NOI08290 $ +D5*ETDT*ETA*(13.5D0+4.D0*ETA2))*SING+XGDT1*(C9*COSG NOI08300 $ -2.D0*C8*SIN2G)) NOI08310 D25=EDOT**2 NOI08320 D17=XNDDT/XNODP-XNDTN**2 NOI08330 TSDDTS=2.D0*TSDTTS*(TSDTTS-D2O)+AODP*TSI*(TOTHRD*BETAO2*D17 NOI08340 $ -4.D0*D20*ECC*EDOT+2.D0*(D25+ECC*EDDOT)) NOI08350 ETDDT=(EDDOT+2.D0*EDOT*TSDTTS)*TSI*S+TSDDTS*ETA NOI08360 D18=TSDDTS-TSDTTS**2 NOI08370 D19=-PSDTPS**2/ETA2-ETA*ETDDT*PSIM2-PSDTPS**2 NOI08380 D23=ETDT**2 NOI08390 D1DDT=D1DT*(D14+D15)+D1*(D18-2.D0*D19+TOTHRD*D17+2.D0*(ALPHA2*D25 NOI08400 $ /BETAO2+ECC*EDDOT)/BETAO2) NOI08410 XNTRDT=XNDT*(2.D0*TOTHRD*D17+3.D0*(D25+ECC*EDDOT)/ALPHA2-6.D0* NOI08420 $ ALDTAL**2+4.D0*D18-7.D0*D19)+C1DTC1*XNDDT+C1*(C1DTC1*D16+D9*ETDDTNOI08430 $ +D10*EDDOT+D23*(6.D0+30.D0*EETA+68.D0*ECCSQ)+ETDT*EDOT*(40.D0+ NOI08440 $ 30.D0*ETA2+272.D0*EETA)+D25*(17.D0+68.D0*ETA2)+B1*(D1DDT*D2+2.D0*NOI08450 $ D1DT*D2DT+D1*(ETDDT*D11+D23*(72.D0+54.D0*ETA2)))+B2*(D1DDT*D3+ NOI08460 $ 2.D0*D1DT*D3DT+D1*(ETDDT*D12+D23*(30.D0+30.D0*ETA2)))+COS2G+B3* NOI08470 $ ((D5DT*D14+D5*(D18-2.D0*D19))*D4+2.D0*D4DT*D5DT+D5*(ETDDT*D13+ NOI08480 $ 22.5D0*ETA*D23))*SING+XGDT1*((7.D0*D20+4.D0*ECC*EDOT/BETAO2)*(C5 NOI08490 $ *COSG-2.D0*C4*SIN2G)+((2.D0*C5DT*COSG-4.D0*C4DT*SIN2G)-XGDT1*(C5 NOI08500 $ *SING+4.D0*C4*COS2G)))) NOI08510 TMNDDT=XNDDT*1.D9 NOI08520 TEMP=TMNDDT**2-XNDT*1.D18*XNTRDT NOI08530 PP=(TEMP+TMNDDT**2)/TEMP NOI08540 GAMMA=-XNTRDT/(XNDDT*(PP-2.D0)) NOI08550 XND=XNDT/(PP*GAMMA) NOI08560 QQ=1.D0-EDDOT/(EDOT*GAMMA) NOI08570 ED=EDOT/(QQ*GAMMA) NOI08580 OVGPP=1.D0/(GAMMA*(PP+1.D0)) NOI08590 GO TO 2 NOI08600 C NOI08610 1 ISIMP=1 NOI08620 EDOT=-TOTHRD*XNDTN*(1.D0-ECC) NOI08630 C NOI08640 C Update for secular gravity and atmospheric drag NOI08650 C NOI08660 2 XMAM=FMOD2P(TETAE+XLLDOT*TSINCE) NOI08670 OMGASM=OMPE+OMGDT*TSINCE NOI08680 XNODES=RAAN+XNODOT*TSINCE NOI08690 IF(ISIMP.EQ.1) GO TO 3 NOI08700 TEMP=1.D0-GAMMA*TSINCE NOI08710 TEMP1=TEMP**PP NOI08720 XN=XNODP+XND*(1.D0-TEMP1) NOI08730 EM=ECC+ED*(1.D0-TEMP**QQ) NOI08740 Z1=XND*(TSINCE+OVGPP*(TEMP*TEMP1-1.D0)) NOI08750 GO TO 4 NOI08760 C NOI08770 3 XN=XNODP+XNDT*TSINCE NOI08780 EM=ECC+EDOT*TSINCE NOI08790 Z1=0.5D0*XNDT*TSINCE**2 NOI08800 C NOI08810 4 Z7=3.5D0*TOTHRD*Z1/XNODP NOI08820 XMAM=FMOD2P(XMAM+Z1+Z7*XMDT1) NOI08830 OMGASM=OMGASM+Z7*XGDT1 NOI08840 XNODES=XNODES+Z7*XHDT1 NOI08850 C NOI08860 C Solve Keplers equation NOI08870 C NOI08880 ZC2=XMAM+EM*DSIN(XMAM)*(1.D0+EM*DCOS(XMAM)) NOI08890 DO 5 I=1,10 NOI08900 SINE=DSIN(ZC2) NOI08910 COSE=DCOS(ZC2) NOI08920 ZC5=1.D0/(1.D0-EM*COSE) NOI08930 CAPE=(XMAM+EM*SINE-ZC2)*ZC5+ZC2 NOI08940 IF(DABS(CAPE-ZC2).LE.ZERO) GO TO 6 NOI08950 5 ZC2=CAPE NOI08960 C NOI08970 C Short period preliminary quantities NOI08980 C NOI08990 6 AM=(XKE/XN)**TOTHRD NOI09000 BETA2M=1.D0-EM**2 NOI09010 SINOS=DSIN(OMGASM) NOI09020 COSOS=DCOS(OMGASM) NOI09030 AXNM=EM*COSOS NOI09040 AYNM=EM*SINOS NOI09050 PM=AM*BETA2M NOI09060 G1=1.D0/PM NOI09070 G2=0.5D0*CK2*G1 NOI09080 G3=G2*G1 NOI09090 BETA=DSQRT(BETA2M) NOI09100 G4=0.25D0*A3COF*SINI NOI09110 G5=0.25D0*A3COF*G1 NOI09120 SNF=BETA*SINE*ZC5 NOI09130 CSF=(COSE-EM)*ZC5 NOI09140 CALL FTAN (SNF,CSF,FM) NOI09150 FM=FM/CW NOI09160 SNFG=SNF*COSOS+CSF*SINOS NOI09170 CSFG=CSF*COSOS-SNF*SINOS NOI09180 SN2F2G=2.D0*SNFG*CSFG NOI09190 CS2F2G=2.D0*CSFG**2-1.D0 NOI09200 ECOSF=EM*CSF NOI09210 G10=FM-XMAM+EM*SNF NOI09220 RM=PM/(1.D0+ECOSF) NOI09230 AOVR=AM/RM NOI09240 G13=XN*AOVR NOI09250 G14=-G13*AOVR NOI09260 DR=G2*(UNMTH2*CS2F2G-3.D0*TTHMUN)-G4*SNFG NOI09270 DIWC=3.D0*G3*SINI*CS2F2G-G5*AYNM NOI09280 DI=DIWC*COSI NOI09290 C NOI09300 C Update for short periodics NOI09310 C NOI09320 SNI2DU=SINIO2*(G3*(0.5D0*(1.D0-7.D0*THETA2)*SN2F2G-3.D0*UNM5TH* NOI09330 $ G10)-G5*SINI*CSFG*(2.D0+ECOSF))-0.5D0*G5*THETA2*AXNM/COSIO2 NOI09340 XLAMB=FM+OMGASM+XNODES+G3*(0.5D0*(1.D0+6.D0*COSI-7.D0*THETA2)* NOI09350 $ SN2F2G-3.D0*(UNM5TH+2.D0*COSI)*G10)+G5*SINI*(COSI*AXNM/(1.D0+ NOI09360 $ COSI)-(2.D0+ECOSF)*CSFG) NOI09370 Y4=SINIO2*SNFG+CSFG*SNI2DU+0.5D0*SNFG*COSIO2*DI NOI09380 Y5=SINIO2*CSFG-SNFG*SNI2DU+0.5D0*CSFG*COSIO2*DI NOI09390 R=RM+DR NOI09400 RDOT=XN*AM*EM*SNF/BETA+G14*(2.D0*G2*UNMTH2*SN2F2G+G4*CSFG) NOI09410 RVDOT=XN*AM**2*BETA/RM+G14*DR+AM*G13*SINI*DIWC NOI09420 C NOI09430 C Orientation vectors NOI09440 C NOI09450 SNLAMB=DSIN(XLAMB) NOI09460 CSLAMB=DCOS(XLAMB) NOI09470 TEMP=2.D0*(Y5*SNLAMB-Y4*CSLAMB) NOI09480 UX=Y4*TEMP+CSLAMB NOI09490 VX=Y5*TEMP-SNLAMB NOI09500 TEMP=2.D0*(Y5*CSLAMB+Y4*SNLAMB) NOI09510 UY=-Y4*TEMP+SNLAMB NOI09520 VY=-Y5*TEMP+CSLAMB NOI09530 TEMP=2.D0*DSQRT(1.D0-Y4**2-Y5**2) NOI09540 UZ=Y4*TEMP NOI09550 VZ=Y5*TEMP NOI09560 C NOI09570 C Position and velocity NOI09580 C NOI09590 X=R*UX NOI09600 Y=R*UY NOI09610 Z=R*UZ NOI09620 C NOI09630 XD=RDOT*UX+RVDOT*VX NOI09640 YD=RDOT*UY+RVDOT*VY NOI09650 ZD=RDOT*UZ+RVDOT*VZ NOI09660 RETURN NOI09670 C NOI09680 END NOI09690 SUBROUTINE SDP4 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI09700 C NOI09710 C----------------------------------------------------------------------+NOI09720 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI09730 C----------------------------------------------------------------------|NOI09740 C Title : SDP4 |NOI09750 C----------------------------------------------------------------------|NOI09760 C Function : Orbital integrator |NOI09770 C----------------------------------------------------------------------|NOI09780 C Author : NORAD |NOI09790 C----------------------------------------------------------------------|NOI09800 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI09810 C----------------------------------------------------------------------|NOI09820 C Date-coded : 11 October 1988 |NOI09830 C----------------------------------------------------------------------|NOI09840 C Last update : 14 December 1988 |NOI09850 C----------------------------------------------------------------------|NOI09860 C Notes : None |NOI09870 C----------------------------------------------------------------------|NOI09880 C Input data : |NOI09890 C----------------------------------------------------------------------|NOI09900 C Output data : |NOI09910 C----------------------------------------------------------------------+NOI09920 C NOI09930 IMPLICIT REAL*8(A-H,O-Z) NOI09940 COMMON /CONST/ DUM1,CW,REQ NOI09950 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI09960 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI09970 COMMON /CM2/ OMPP,OMGP NOI09980 C NOI09990 IF(IFLAG.EQ.0) GO TO 2 NOI10000 C NOI10010 C Recover original mean motion (XNODP) and semimajor axis (AODP) NOI10020 C from input elements NOI10030 C NOI10040 A1=(XKE/XNO)**TOTHRD NOI10050 COSIO=DCOS(CLIN) NOI10060 THETA2=COSIO**2 NOI10070 X3THM1=3.D0*THETA2-1.D0 NOI10080 EOSQ=ECC**2 NOI10090 BETAO2=1.D0-EOSQ NOI10100 BETAO=DSQRT(BETAO2) NOI10110 DEL1=1.5D0*CK2*X3THM1/(A1**2*BETAO*BETAO2) NOI10120 AO=A1*(1.D0-DEL1*(0.5D0*TOTHRD+DEL1*(1.D0+134.D0/81.D0*DEL1))) NOI10130 DELO=1.5*CK2*X3THM1/(AO**2*BETAO*BETAO2) NOI10140 XNODP=XNO/(1.D0+DELO) NOI10150 AODP=AO/(1.D0-DELO) NOI10160 C NOI10170 C Initialization : NOI10180 C for perigee below 156 km, the values of S and QOMS2T are altered. NOI10190 C NOI10200 S4=S NOI10210 QOMS24=QOMS2T NOI10220 PERIGE=(AODP*(1.D0-ECC)-AE)*REQ NOI10230 IF(PERIGE.GE.156) GO TO 1 NOI10240 S4=PERIGE-78.D0 NOI10250 IF(PERIGE.LE.98) S4=20.D0 NOI10260 QOMS24=((120.D0-S4)*AE/REQ)**4 NOI10270 S4=S4/REQ+AE NOI10280 C NOI10290 1 PINVSQ=1.D0/(AODP**2*BETAO2**2) NOI10300 SING=DSIN(OMPE) NOI10310 COSG=DCOS(OMPE) NOI10320 TSI=1.D0/(AODP-S4) NOI10330 ETA=AODP*ECC*TSI NOI10340 ETASQ=ETA**2 NOI10350 EETA=ECC*ETA NOI10360 PSISQ=DABS(1.D0-ETASQ) NOI10370 COEF=QOMS24*TSI**4 NOI10380 COEF1=COEF/PSISQ**3.5D0 NOI10390 C2=COEF1*XNODP*(AODP*(1.D0+1.5D0*ETASQ+EETA*(4.D0+ETASQ)) NOI10400 $ +0.75D0*CK2*TSI/PSISQ*X3THM1*(8.D0+3.D0*ETASQ*(8.D0+ETASQ))) NOI10410 C1=BSTAR*C2 NOI10420 SINIO=DSIN(CLIN) NOI10430 A3OVK2=-XJ3/CK2*AE**3 NOI10440 X1MTH2=1.D0-THETA2 NOI10450 C4=2.D0*XNODP*COEF1*AODP*BETAO2*(ETA* NOI10460 $ (2.D0+0.5D0*ETASQ)+ECC*(0.5D0+2.D0*ETASQ)-2.D0*CK2*TSI/ NOI10470 $ (AODP*PSISQ)*(-3.D0*X3THM1*(1.D0-2.D0*EETA+ETASQ* NOI10480 $ (1.5D0-0.5D0*EETA))+0.75D0*X1MTH2*(2.D0*ETASQ-EETA* NOI10490 $ (1.D0+ETASQ))*DCOS(2.D0*OMPE))) NOI10500 THETA4=THETA2**2 NOI10510 TEMP1=3.D0*CK2*PINVSQ*XNODP NOI10520 TEMP2=TEMP1*CK2*PINVSQ NOI10530 TEMP3=1.25D0*CK4*PINVSQ**2*XNODP NOI10540 XMDOT=XNODP+0.5D0*TEMP1*BETAO*X3THM1+0.0625D0*TEMP2*BETAO* NOI10550 $ (13.D0-78.D0*THETA2+137.D0*THETA4) NOI10560 X1M5TH=1.D0-5.D0*THETA2 NOI10570 OMGDOT=-0.5D0*TEMP1*X1M5TH+0.0625D0*TEMP2*(7.D0-114.D0*THETA2+ NOI10580 $ 395.D0*THETA4)+TEMP3*(3.D0-36.D0*THETA2+49.D0*THETA4) NOI10590 OMPP= OMGDOT*1440.D0*CW NOI10600 XHDOT1=-TEMP1*COSIO NOI10610 XNODOT=XHDOT1+(.5D0*TEMP2*(4.D0-19.D0*THETA2)+2.D0*TEMP3*(3.D0- NOI10620 $ 7.D0*THETA2))*COSIO NOI10630 OMGP=XNODOT*1440.D0*CW NOI10640 XNODCF=3.5D0*BETAO2*XHDOT1*C1 NOI10650 T2COF=1.5D0*C1 NOI10660 XLCOF=0.125D0*A3OVK2*SINIO*(3.D0+5.D0*COSIO)/(1.D0+COSIO) NOI10670 AYCOF=0.25D0*A3OVK2*SINIO NOI10680 X7THM1=7.D0*THETA2-1.D0 NOI10690 IFLAG=0 NOI10700 CALL DPINIT (EOSQ,SINIO,COSIO,BETAO,AODP,THETA2, NOI10710 $ SING,COSG,BETAO2,XMDOT,OMGDOT,XNODOT,XNODP) NOI10720 C NOI10730 C Update for secular gravity and atmospheric drag NOI10740 C NOI10750 2 XMDF=TETAE+XMDOT*TSINCE NOI10760 OMGADF=OMPE+OMGDOT*TSINCE NOI10770 XNODDF=RAAN+XNODOT*TSINCE NOI10780 TSQ=TSINCE**2 NOI10790 XNODE=XNODDF+XNODCF*TSQ NOI10800 TEMPA=1.D0-C1*TSINCE NOI10810 TEMPE=BSTAR*C4*TSINCE NOI10820 TEMPL=T2COF*TSQ NOI10830 XN=XNODP NOI10840 CALL DPSEC (XMDF,OMGADF,XNODE,EM,XINC,XN,TSINCE) NOI10850 A=(XKE/XN)**TOTHRD*TEMPA**2 NOI10860 E=EM-TEMPE NOI10870 XMAM=XMDF+XNODP*TEMPL NOI10880 CALL DPPER (E,XINC,OMGADF,XNODE,XMAM) NOI10890 XL=XMAM+OMGADF+XNODE NOI10900 BETA=DSQRT(1.D0-E**2) NOI10910 XN=XKE/A**1.5D0 NOI10920 C NOI10930 C Long period periodics NOI10940 C NOI10950 AXN=E*DCOS(OMGADF) NOI10960 TEMP=1.D0/(A*BETA**2) NOI10970 XLL=TEMP*XLCOF*AXN NOI10980 AYNL=TEMP*AYCOF NOI10990 XLT=XL+XLL NOI11000 AYN=E*DSIN(OMGADF)+AYNL NOI11010 C NOI11020 C Solve Keplers equation NOI11030 C NOI11040 CAPU=FMOD2P(XLT-XNODE) NOI11050 TEMP2=CAPU NOI11060 DO 3 I=1,10 NOI11070 SINEPW=DSIN(TEMP2) NOI11080 COSEPW=DCOS(TEMP2) NOI11090 TEMP3=AXN*SINEPW NOI11100 TEMP4=AYN*COSEPW NOI11110 TEMP5=AXN*COSEPW NOI11120 TEMP6=AYN*SINEPW NOI11130 EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.D0-TEMP5-TEMP6)+TEMP2 NOI11140 IF(DABS(EPW-TEMP2).LE.ZERO) GO TO 4 NOI11150 3 TEMP2=EPW NOI11160 C NOI11170 C Short period preliminary quantities NOI11180 C NOI11190 4 ECOSE=TEMP5+TEMP6 NOI11200 ESINE=TEMP3-TEMP4 NOI11210 ELSQ=AXN**2+AYN**2 NOI11220 TEMP=1.D0-ELSQ NOI11230 PL=A*TEMP NOI11240 R=A*(1.D0-ECOSE) NOI11250 TEMP1=1.D0/R NOI11260 RDOT=XKE*DSQRT(A)*ESINE*TEMP1 NOI11270 RFDOT=XKE*DSQRT(PL)*TEMP1 NOI11280 TEMP2=A*TEMP1 NOI11290 BETAL=DSQRT(TEMP) NOI11300 TEMP3=1.D0/(1.D0+BETAL) NOI11310 COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3) NOI11320 SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3) NOI11330 CALL FTAN (SINU,COSU,U) NOI11340 U=U/CW NOI11350 SIN2U=2.D0*SINU*COSU NOI11360 COS2U=2.D0*COSU**2-1.D0 NOI11370 TEMP=1.D0/PL NOI11380 TEMP1=CK2*TEMP NOI11390 TEMP2=TEMP1*TEMP NOI11400 C NOI11410 C Update for short periodics NOI11420 C NOI11430 RK=R*(1.D0-1.5D0*TEMP2*BETAL*X3THM1)+0.5D0*TEMP1*X1MTH2*COS2U NOI11440 UK=U-0.25D0*TEMP2*X7THM1*SIN2U NOI11450 XNODEK=XNODE+1.5D0*TEMP2*COSIO*SIN2U NOI11460 XINCK=XINC+1.5D0*TEMP2*COSIO*SINIO*COS2U NOI11470 RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U NOI11480 RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5D0*X3THM1) NOI11490 C NOI11500 C Orientation vectors NOI11510 C NOI11520 SINUK=DSIN(UK) NOI11530 COSUK=DCOS(UK) NOI11540 SINIK=DSIN(XINCK) NOI11550 COSIK=DCOS(XINCK) NOI11560 SINNOK=DSIN(XNODEK) NOI11570 COSNOK=DCOS(XNODEK) NOI11580 XMX=-SINNOK*COSIK NOI11590 XMY=COSNOK*COSIK NOI11600 C NOI11610 UX=XMX*SINUK+COSNOK*COSUK NOI11620 UY=XMY*SINUK+SINNOK*COSUK NOI11630 UZ=SINIK*SINUK NOI11640 C NOI11650 VX=XMX*COSUK-COSNOK*SINUK NOI11660 VY=XMY*COSUK-SINNOK*SINUK NOI11670 VZ=SINIK*COSUK NOI11680 C NOI11690 C Position and velocity NOI11700 C NOI11710 X=RK*UX NOI11720 Y=RK*UY NOI11730 Z=RK*UZ NOI11740 C NOI11750 XD=RDOTK*UX+RFDOTK*VX NOI11760 YD=RDOTK*UY+RFDOTK*VY NOI11770 ZD=RDOTK*UZ+RFDOTK*VZ NOI11780 RETURN NOI11790 C NOI11800 END NOI11810 SUBROUTINE SDP8 (IFLAG,TSINCE,X,Y,Z,XD,YD,ZD) NOI11820 C NOI11830 C----------------------------------------------------------------------+NOI11840 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI11850 C----------------------------------------------------------------------|NOI11860 C Title : SDP8 |NOI11870 C----------------------------------------------------------------------|NOI11880 C Function : Orbital integrator |NOI11890 C----------------------------------------------------------------------|NOI11900 C Author : NORAD |NOI11910 C----------------------------------------------------------------------|NOI11920 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI11930 C----------------------------------------------------------------------|NOI11940 C Date-coded : 11 October 1988 |NOI11950 C----------------------------------------------------------------------|NOI11960 C Last update : 14 December 1988 |NOI11970 C----------------------------------------------------------------------|NOI11980 C Notes : None |NOI11990 C----------------------------------------------------------------------|NOI12000 C Input data : |NOI12010 C----------------------------------------------------------------------|NOI12020 C Output data : |NOI12030 C----------------------------------------------------------------------+NOI12040 C NOI12050 IMPLICIT REAL*8 (A-H,O-Z) NOI12060 COMMON /CONST/ DUM1,CW,DUM2 NOI12070 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI12080 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI12090 COMMON /CM2/ OMPP,OMGP NOI12100 DATA RHO/.15696615/ NOI12110 C NOI12120 IF(IFLAG.EQ.0) GO TO 1 NOI12130 C NOI12140 C Recover original mean motion (XNODP) and semimajor axis (AODP) from NOI12150 C input elements ------ calculate ballistic coefficient (B term) from NOI12160 C input B* drag term NOI12170 C NOI12180 A1=(XKE/XNO)**TOTHRD NOI12190 COSI=DCOS(CLIN) NOI12200 THETA2=COSI*COSI NOI12210 TTHMUN=3.D0*THETA2-1.D0 NOI12220 EQSQ=ECC**2 NOI12230 BETAO2=1.D0-EQSQ NOI12240 BETAO=DSQRT(BETAO2) NOI12250 DEL1=1.5D0*CK2*TTHMUN/(A1*A1*BETAO*BETAO2) NOI12260 AO=A1*(1.D0-DEL1*(.5D0*TOTHRD+DEL1*(1.D0+134.D0/81.D0*DEL1))) NOI12270 DELO=1.5D0*CK2*TTHMUN/(AO*AO*BETAO*BETAO2) NOI12280 AODP=AO/(1.D0-DELO) NOI12290 XNODP=XNO/(1.D0+DELO) NOI12300 B=2.D0*BSTAR/RHO NOI12310 C NOI12320 C Initialization NOI12330 C NOI12340 PO=AODP*BETAO2 NOI12350 POM2=1.D0/(PO*PO) NOI12360 SINI=DSIN(CLIN) NOI12370 SING=DSIN(OMPE) NOI12380 COSG=DCOS(OMPE) NOI12390 TEMP=0.5D0*CLIN NOI12400 SINIO2=DSIN(TEMP) NOI12410 COSIO2=DCOS(TEMP) NOI12420 THETA4=THETA2**2 NOI12430 UNM5TH=1.D0-5.D0*THETA2 NOI12440 UNMTH2=1.D0-THETA2 NOI12450 A3COF=-XJ3/CK2*AE**3 NOI12460 PARDT1=3.D0*CK2*POM2*XNODP NOI12470 PARDT2=PARDT1*CK2*POM2 NOI12480 PARDT4=1.25D0*CK4*POM2*POM2*XNODP NOI12490 XMDT1=0.5D0*PARDT1*BETAO*TTHMUN NOI12500 XGDT1=-0.5D0*PARDT1*UNM5TH NOI12510 XHDT1=-PARDT1*COSI NOI12520 XLLDOT=XNODP+XMDT1+0.0625D0*PARDT2*BETAO*(13.D0-78.D0*THETA2+ NOI12530 $ 137.D0*THETA4) NOI12540 OMGDT=XGDT1+0.0625D0*PARDT2*(7.D0-114.D0*THETA2+395.D0*THETA4)+ NOI12550 $ PARDT4*(3.D0-36.D0*THETA2+49.D0*THETA4) NOI12560 OMPP= OMGDT*1440.D0*CW NOI12570 XNODOT=XHDT1+(.5D0*PARDT2*(4.D0-19.D0*THETA2)+2.D0*PARDT4*( NOI12580 $ 3.D0-7.D0*THETA2))*COSI NOI12590 OMGP=XNODOT*1440.D0*CW NOI12600 TSI=1.D0/(PO-S) NOI12610 ETA=ECC*S*TSI NOI12620 ETA2=ETA**2 NOI12630 PSIM2=DABS(1.D0/(1.D0-ETA2)) NOI12640 ALPHA2=1.D0+EQSQ NOI12650 EETA=ECC*ETA NOI12660 COS2G=2.D0*COSG**2-1.D0 NOI12670 D5=TSI*PSIM2 NOI12680 D1=D5/PO NOI12690 D2=12.D0+ETA2*(36.D0+4.5D0*ETA2) NOI12700 D3=ETA2*(15.D0+2.5D0*ETA2) NOI12710 D4=ETA*(5.D0+3.75D0*ETA2) NOI12720 B1=CK2*TTHMUN NOI12730 B2=-CK2*UNMTH2 NOI12740 B3=A3COF*SINI NOI12750 C0=0.5D0*B*RHO*QOMS2T*XNODP*AODP*TSI**4*PSIM2**(14.D0/4.D0)/ NOI12760 $ DSQRT(ALPHA2) NOI12770 C1=1.5D0*XNODP*ALPHA2**2*C0 NOI12780 C4=D1*D3*B2 NOI12790 C5=D5*D4*B3 NOI12800 XNDT=C1*((2.D0+ETA2*(3.D0+34.D0*EQSQ)+5.D0*EETA*(4.D0+ETA2)+ NOI12810 $ 8.5D0*EQSQ)+D1*D2*B1+C4*COS2G+C5*SING) NOI12820 XNDTN=XNDT/XNODP NOI12830 EDOT=-TOTHRD*XNDTN*(1.D0-ECC) NOI12840 IFLAG=0 NOI12850 CALL DPINIT (EQSQ,SINI,COSI,BETAO,AODP,THETA2,SING,COSG, NOI12860 $ BETAO2,XLLDOT,OMGDT,XNODOT,XNODP) NOI12870 C NOI12880 C Update for secular gravity and atmospheric drag NOI12890 C NOI12900 1 Z1=0.5D0*XNDT*TSINCE*TSINCE NOI12910 Z7=3.5D0*TOTHRD*Z1/XNODP NOI12920 XMAMDF=TETAE+XLLDOT*TSINCE NOI12930 OMGASM=OMPE+OMGDT*TSINCE+Z7*XGDT1 NOI12940 XNODES=RAAN+XNODOT*TSINCE+Z7*XHDT1 NOI12950 XN=XNODP NOI12960 CALL DPSEC (XMAMDF,OMGASM,XNODES,EM,XINC,XN,TSINCE) NOI12970 XN=XN+XNDT*TSINCE NOI12980 EM=EM+EDOT*TSINCE NOI12990 XMAM=XMAMDF+Z1+Z7*XMDT1 NOI13000 CALL DPPER (EM,XINC,OMGASM,XNODES,XMAM) NOI13010 XMAM=FMOD2P(XMAM) NOI13020 C NOI13030 C Solve Keplers equation NOI13040 C NOI13050 ZC2=XMAM+EM*DSIN(XMAM)*(1.D0+EM*DCOS(XMAM)) NOI13060 DO 2 I=1,10 NOI13070 SINE=DSIN(ZC2) NOI13080 COSE=DCOS(ZC2) NOI13090 ZC5=1.D0/(1.D0-EM*COSE) NOI13100 CAPE=(XMAM+EM*SINE-ZC2)*ZC5+ZC2 NOI13110 IF(DABS(CAPE-ZC2).LE.ZERO) GO TO 3 NOI13120 2 ZC2=CAPE NOI13130 C NOI13140 C Short period preliminary quantities NOI13150 C NOI13160 3 AM=(XKE/XN)**TOTHRD NOI13170 BETA2M=1.D0-EM*EM NOI13180 SINOS=DSIN(OMGASM) NOI13190 COSOS=DCOS(OMGASM) NOI13200 AXNM=EM*COSOS NOI13210 AYNM=EM*SINOS NOI13220 PM=AM*BETA2M NOI13230 G1=1.D0/PM NOI13240 G2=0.5D0*CK2*G1 NOI13250 G3=G2*G1 NOI13260 BETA=DSQRT(BETA2M) NOI13270 G4=0.25D0*A3COF*SINI NOI13280 G5=0.25D0*A3COF*G1 NOI13290 SNF=BETA*SINE*ZC5 NOI13300 CSF=(COSE-EM)*ZC5 NOI13310 CALL FTAN (SNF,CSF,FM) NOI13320 FM=FM/CW NOI13330 SNFG=SNF*COSOS+CSF*SINOS NOI13340 CSFG=CSF*COSOS-SNF*SINOS NOI13350 SN2F2G=2.D0*SNFG*CSFG NOI13360 CS2F2G=2.D0*CSFG**2-1.D0 NOI13370 ECOSF=EM*CSF NOI13380 G10=FM-XMAM+EM*SNF NOI13390 RM=PM/(1.D0+ECOSF) NOI13400 AOVR=AM/RM NOI13410 G13=XN*AOVR NOI13420 G14=-G13*AOVR NOI13430 DR=G2*(UNMTH2*CS2F2G-3.D0*TTHMUN)-G4*SNFG NOI13440 DIWC=3.D0*G3*SINI*CS2F2G-G5*AYNM NOI13450 DI=DIWC*COSI NOI13460 SINI2=DSIN(.5D0*XINC) NOI13470 C NOI13480 C Update for short period periodics NOI13490 C NOI13500 SNI2DU=SINIO2*(G3*(.5D0*(1.D0-7.D0*THETA2)*SN2F2G-3.D0*UNM5TH* NOI13510 $ G10)-G5*SINI*CSFG*(2.D0+ECOSF))-0.5D0*G5*THETA2*AXNM/COSIO2 NOI13520 XLAMB=FM+OMGASM+XNODES+G3*(.5D0*(1.D0+6.D0*COSI-7.D0*THETA2)* NOI13530 $ SN2F2G-3.D0*(UNM5TH+2.D0*COSI)*G10)+G5*SINI*(COSI*AXNM/ NOI13540 $ (1.D0+COSI)-(2.D0+ECOSF)*CSFG) NOI13550 Y4=SINI2*SNFG+CSFG*SNI2DU+0.5D0*SNFG*COSIO2*DI NOI13560 Y5=SINI2*CSFG-SNFG*SNI2DU+0.5D0*CSFG*COSIO2*DI NOI13570 R=RM+DR NOI13580 RDOT=XN*AM*EM*SNF/BETA+G14*(2.D0*G2*UNMTH2*SN2F2G+G4*CSFG) NOI13590 RVDOT=XN*AM**2*BETA/RM+G14*DR+AM*G13*SINI*DIWC NOI13600 C NOI13610 C Orientation vectors NOI13620 C NOI13630 SNLAMB=DSIN(XLAMB) NOI13640 CSLAMB=DCOS(XLAMB) NOI13650 TEMP=2.D0*(Y5*SNLAMB-Y4*CSLAMB) NOI13660 UX=Y4*TEMP+CSLAMB NOI13670 VX=Y5*TEMP-SNLAMB NOI13680 TEMP=2.D0*(Y5*CSLAMB+Y4*SNLAMB) NOI13690 UY=-Y4*TEMP+SNLAMB NOI13700 VY=-Y5*TEMP+CSLAMB NOI13710 TEMP=2.D0*DSQRT(1.D0-Y4*Y4-Y5*Y5) NOI13720 UZ=Y4*TEMP NOI13730 VZ=Y5*TEMP NOI13740 C NOI13750 C Position and velocity NOI13760 C NOI13770 X=R*UX NOI13780 Y=R*UY NOI13790 Z=R*UZ NOI13800 XD=RDOT*UX+RVDOT*VX NOI13810 YD=RDOT*UY+RVDOT*VY NOI13820 ZD=RDOT*UZ+RVDOT*VZ NOI13830 RETURN NOI13840 C NOI13850 END NOI13860 SUBROUTINE DPINIT (EQSQ,SINIQ,COSIQ,RTEQSQ,AO,COSQ2,SINOMO,COSOMO,NOI13870 $ BSQ,XLLDOT,OMGDT,XNODOT,XNODP) NOI13880 C NOI13890 C----------------------------------------------------------------------+NOI13900 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI13910 C----------------------------------------------------------------------|NOI13920 C Title : DPINIT |NOI13930 C----------------------------------------------------------------------|NOI13940 C Function : Entrance for deep space initialization NOI13950 C----------------------------------------------------------------------|NOI13960 C Author : NORAD |NOI13970 C----------------------------------------------------------------------|NOI13980 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI13990 C----------------------------------------------------------------------|NOI14000 C Date-coded : 11 October 1988 |NOI14010 C----------------------------------------------------------------------|NOI14020 C Last update : 14 December 1988 |NOI14030 C----------------------------------------------------------------------|NOI14040 C Notes : None |NOI14050 C----------------------------------------------------------------------|NOI14060 C Input data : |NOI14070 C----------------------------------------------------------------------|NOI14080 C Output data : |NOI14090 C----------------------------------------------------------------------+NOI14100 C NOI14110 IMPLICIT REAL*8 (A-H,O-Z) NOI14120 INTEGER ITX(6) NOI14130 COMMON /CONST/ PIG,CW,REQ NOI14140 COMMON /ORB/ ITE(6),DUM1(9) NOI14150 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI14160 C NOI14170 DATA ZNS,C1SS,ZES/1.19459D-5,2.9364797D-6,0.01675D0/ NOI14180 DATA ZNL, C1L,ZEL/1.5835218D-4,4.7968035D-7,0.05490D0/ NOI14190 DATA ZCOSIS,ZSINIS,ZSINGS/0.91744867D0,0.39785416D0,-0.98088458D0/NOI14200 DATA ZCOSGS/0.1945905D0/ NOI14210 DATA Q22,Q31,Q33/1.7891679D-6,2.1460748D-6,2.2123015D-7/ NOI14220 DATA ROOT22,ROOT32/1.7891679D-6,3.7393792D-7/ NOI14230 DATA ROOT44,ROOT52/7.3636953D-9,1.1428639D-7/ NOI14240 DATA ROOT54/2.1765803D-9/ NOI14250 DATA THDT/4.3752691D-3/ NOI14260 C NOI14270 DATA ITX/1988, 1, 1, 0, 0, 0/ NOI14280 ITX(1)=ITE(1) NOI14290 CALL DT (ITE,ITX,FDG) NOI14300 EPOCH=(ITE(1)-1900.D0)*1000.D0+FDG NOI14310 CALL CICCIO (THETAG,DS50,EPOCH) NOI14320 C NOI14330 EQ=ECC NOI14340 XNQ=XNODP NOI14350 AQNV=1.D0/AO NOI14360 XQNCL=CLIN NOI14370 XMAO=TETAE NOI14380 XPIDOT=OMGDT+XNODOT NOI14390 SINQ=DSIN(RAAN) NOI14400 COSQ=DCOS(RAAN) NOI14410 OMEGAQ=OMPE NOI14420 C NOI14430 C Initialize lunar solar terms NOI14440 C NOI14450 5 DAY=DS50+18261.5D0 NOI14460 IF(DAY.EQ.PREEP) GO TO 10 NOI14470 PREEP=DAY NOI14480 XNODCE=4.5236020D0-9.2422029D-4*DAY NOI14490 STEM=DSIN(XNODCE) NOI14500 CTEM=DCOS(XNODCE) NOI14510 ZCOSIL=0.91375164D0-0.03568096D0*CTEM NOI14520 ZSINIL=DSQRT(1.D0-ZCOSIL*ZCOSIL) NOI14530 ZSINHL=0.089683511D0*STEM/ZSINIL NOI14540 ZCOSHL=DSQRT(1.D0-ZSINHL*ZSINHL) NOI14550 C=4.7199672D0+0.22997150D0+DAY NOI14560 GAM=5.8351514D0+0.001944368D0*DAY NOI14570 ZMOL=FMOD2P(C-GAM) NOI14580 ZX=0.39785416D0*STEM/ZSINIL NOI14590 ZY=ZCOSHL*CTEM+0.91744367D0*ZSINHL*STEM NOI14600 CALL FTAN (ZX,ZY,QWER) NOI14610 ZX=QWER/CW NOI14620 ZX=GAM+ZX-XNODCE NOI14630 ZCOSGL=DCOS(ZX) NOI14640 ZSINGL=DSIN(ZX) NOI14650 ZMOS=6.2565837D0+0.017201977D0*DAY NOI14660 ZMOS=FMOD2P(ZMOS) NOI14670 C NOI14680 C Do solar terms NOI14690 C NOI14700 10 LS=0 NOI14710 SAVTSN=1.D20 NOI14720 ZCOSG=ZCOSGS NOI14730 ZSING=ZSINGS NOI14740 ZCOSI=ZCOSIS NOI14750 ZSINI=ZSINIS NOI14760 ZCOSH=COSQ NOI14770 ZSINH=SINQ NOI14780 CC=C1SS NOI14790 ZN=ZNS NOI14800 ZE=ZES NOI14810 ZMO=ZMOS NOI14820 XNOI=1.D0/XNQ NOI14830 ASSIGN 30 TO LS NOI14840 20 A1=ZCOSG*ZCOSH+ZSING*ZCOSI*ZSINH NOI14850 A3=-ZSING*ZCOSH+ZCOSG*ZCOSI*ZSINH NOI14860 A7=-ZCOSG*ZSINH+ZSING*ZCOSI*ZCOSH NOI14870 A8=ZSING*ZSINI NOI14880 A9=ZSING*ZSINH+ZCOSG*ZCOSI*ZCOSH NOI14890 A10=ZCOSG*ZSINI NOI14900 A2=COSIQ*A7+SINIQ*A8 NOI14910 A4=COSIQ*A9+SINIQ*A10 NOI14920 A5=-SINIQ*A7+COSIQ*A8 NOI14930 A6=-SINIQ*A9+COSIQ*A10 NOI14940 C NOI14950 X1=A1*COSOMO+A2*SINOMO NOI14960 X2=A3*COSOMO+A4*SINOMO NOI14970 X3=-A1*SINOMO+A2*COSOMO NOI14980 X4=-A3*SINOMO+A4*COSOMO NOI14990 X5=A5*SINOMO NOI15000 X6=A6*SINOMO NOI15010 X7=A5*COSOMO NOI15020 X8=A6*COSOMO NOI15030 C NOI15040 Z31=12.D0*X1*X1-3.D0*X3*X3 NOI15050 Z32=24.D0*X1*X2-6.D0*X3*X4 NOI15060 Z33=12.D0*X2*X2-3.D0*X4*X4 NOI15070 Z1=3.D0*(A1*A1+A2*A2)+Z31*EQSQ NOI15080 Z2=6.D0*(A1*A3+A2*A4)+Z32*EQSQ NOI15090 Z3=3.D0*(A3*A3+A4*A4)+Z33*EQSQ NOI15100 Z11=-6.D0*A1*A5+EQSQ*(-24.D0*X1*X7-6.D0*X3*X5) NOI15110 Z12=-6.D0*(A1*A6+A3*A5)+EQSQ*(-24.D0*(X2*X7+X1*X8)-6.D0* NOI15120 $ (X3*X6+X4*X5)) NOI15130 Z13=-6.D0*A3*A6+EQSQ*(-24.D0*X2*X8-6.D0*X4*X6) NOI15140 Z21=6.D0*A2*A5+EQSQ*(24.D0*X1*X5-6.D0*X3*X7) NOI15150 Z22=6.D0*(A4*A5+A2*A6)+EQSQ*(24.D0*(X2*X5+X1*X6)-8.D0*(X4*X7 NOI15160 $ +X3*X8)) NOI15170 Z23=6.D0*A4*A6+EQSQ*(24.D0*X2*X5-6.D0*X4*X8) NOI15180 Z1=Z1+Z1+BSQ+Z31 NOI15190 Z2=Z2+Z2+BSQ+Z32 NOI15200 Z3=Z3+Z3+BSQ+Z33 NOI15210 S3=CC*XNOI NOI15220 S2=-0.5D0*S3/RTEQSQ NOI15230 S4=S3*RTEQSQ NOI15240 S1=-15.D0*EQ*S4 NOI15250 S5=X1*X3+X2*X4 NOI15260 S6=X2*X3+X1*X4 NOI15270 S7=X2*X4-X1*X3 NOI15280 SE=S1*ZN*S5 NOI15290 SI=S2*ZN*(Z11+Z13) NOI15300 SL=-ZN*S3*(Z1+Z3-14.D0-6.D0*EQSQ) NOI15310 SHG=S4*ZN*(Z31+Z33-6.D0) NOI15320 SH=-ZN*S2*(Z21+Z23) NOI15330 IF(XQNCL.LT.5.2359877D-2) SH=0.D0 NOI15340 EE2=2.D0*S1*S6 NOI15350 E3=2.D0*S1*S7 NOI15360 XI2=2.D0*S2*Z12 NOI15370 XI3=2.D0*S2*(Z13-Z11) NOI15380 XL2=-2.D0*S3*Z2 NOI15390 XL3=-2.D0*S3*(Z3-Z1) NOI15400 XL4=-2.D0*S3*(-21.D0-9.D0*EQSQ)*ZE NOI15410 XGH2=2.D0*S4*Z32 NOI15420 XGH3=2.D0*S4*(Z33-Z31) NOI15430 XGH4=-18.D0*S4*ZE NOI15440 XH2=-2.D0*S2*Z22 NOI15450 XH3=-2.D0*S2*(Z23-Z21) NOI15460 GO TO LS,(30,40) NOI15470 C NOI15480 C Do lunar terms NOI15490 C NOI15500 30 SSE=SE NOI15510 SSI=SI NOI15520 SSL=SL NOI15530 SSH=SH/SINIQ NOI15540 SSG=SGH-COSIQ*SSH NOI15550 SE2=EE2 NOI15560 SI2=XI2 NOI15570 SL2=XL2 NOI15580 SGH2=XGH2 NOI15590 SH2=XH2 NOI15600 SE3=E3 NOI15610 SI3=XI3 NOI15620 SL3=XL3 NOI15630 SGH3=XGH3 NOI15640 SH3=XH3 NOI15650 SL4=XL4 NOI15660 SGH4=XGH4 NOI15670 LS=1 NOI15680 ZCOSG=ZCOSGL NOI15690 ZSING=ZSINGL NOI15700 ZCOSI=ZCOSIL NOI15710 ZSINI=ZSINIL NOI15720 ZCOSH=ZCOSHL*COSQ+ZSINHL*SINQ NOI15730 ZSINH=SINQ*ZCOSHL-COSQ*ZSINHL NOI15740 ZN=ZNL NOI15750 CC=C1L NOI15760 ZE=ZEL NOI15770 ZMO=ZMOL NOI15780 ASSIGN 40 TO LS NOI15790 GO TO 20 NOI15800 40 SSE=SSE+SE NOI15810 SSI=SSI+SI NOI15820 SSL=SSL+SL NOI15830 SSG=SSG+SHG-COSIQ/SINIQ*SH NOI15840 SSH=SSH+SH/SINIQ NOI15850 C NOI15860 C Geopotential resonance initialization for 12 hour orbits NOI15870 C NOI15880 IRESFL=0 NOI15890 ISYNFL=0 NOI15900 IF(XNQ.LT.(.0052359877D0).AND.XNQ.GT.(.0034906585D0)) GO TO 70 NOI15910 IF(XNQ.LT.(8.26D-3).OR.XNQ.GT.(9.24D-3)) RETURN NOI15920 IF(EQ.LT.0.5D0) RETURN NOI15930 IRESFL=1 NOI15940 EOC=EQ*EQSQ NOI15950 G201=-0.306D0-(EQ-0.64D0)*0.44D0 NOI15960 IF(EQ.GT.(.65D0)) GO TO 45 NOI15970 G211=3.616D0-13.247D0*EQ+16.290*EQSQ NOI15980 G310=-19.302D0+117.390D0*EQ-228.419D0*EQSQ+156.591D0*EOC NOI15990 G322=-18.9068D0+109.7927D0*EQ-214.6334D0*EQSQ+146.5816D0*EOC NOI16000 G410=-41.122D0+242.694D0*EQ-471.094D0*EQSQ+313.953D0*EOC NOI16010 G422=-146.407D0+841.880D0*EQ-1629.014D0*EQSQ+1083.435D0*EOC NOI16020 G520=-532.114D0+3017.977D0*EQ-5740.D0*EQSQ+3708.276D0*EOC NOI16030 GO TO 55 NOI16040 C NOI16050 45 G211=-72.099D0+331.819D0*EQ-508.378D0*EQSQ+266.724D0*EOC NOI16060 G310=-346.844D0+1582.851D0*EQ-2415.925D0*EQSQ+1246.113D0*EOC NOI16070 G322=-342.585D0+1554.908D0*EQ-2366.899D0*EQSQ+1215.972D0*EOC NOI16080 G410=-1052.797D0+4758.686D0*EQ-7193.992D0*EQSQ+3651.957D0*EOC NOI16090 G422=-3581.69D0+16178.11D0*EQ-24462.77D0*EQSQ+12422.52D0*EOC NOI16100 IF(EQ.GT.(.715D0)) GO TO 50 NOI16110 G520=1464.74D0-4664.75D0*EQ+3763.64D0*EQSQ NOI16120 GO TO 55 NOI16130 C NOI16140 50 G520=-5149.66D0+29936.92D0*EQ-54087.36D0*EQSQ+31324.56D0*EOC NOI16150 C NOI16160 55 IF(EQ.GE.(.7D0)) GO TO 60 NOI16170 G533=-919.2277D0+4988.61D0*EQ-9064.77D0*EQSQ+5542.21D0*EOC NOI16180 G521=-822.71072D0+4568.6173D0*EQ-8491.4146D0*EQSQ+5337.524D0*EOC NOI16190 G532=-853.666D0+4690.25D0*EQ-8624.77D0*EQSQ+5341.4D0*EOC NOI16200 GO TO 65 NOI16210 C NOI16220 60 G533=-37995.78D0+161616.52D0*EQ-229838.2D0*EQSQ+109377.94D0*EOC NOI16230 G521=-51752.104D0+218913.95D0*EQ-309468.16D0*EQSQ+146349.42D0*EOC NOI16240 G532=-40023.88D0+170470.89D0*EQ-242699.48D0*EQSQ+115605.82D0*EOC NOI16250 65 SINI2=SINIQ*SINIQ NOI16260 F220=0.75D0*(1.D0+2.D0*COSIQ*COSQ2) NOI16270 F221=1.5D0*SINI2 NOI16280 F321=1.875D0*SINIQ*(1.D0-2.D0*COSIQ-3.D0*COSQ2) NOI16290 F322=-1.875D0*SINIQ*(1.D0+2.D0*COSIQ-3.D0*COSQ2) NOI16300 F441=35.D0*SINI2*F220 NOI16310 F442=39.3750D0*SINI2*SINI2 NOI16320 F522=9.84375D0*SINIQ*(SINI2*(1.D0-2.D0*COSIQ-5.D0*COSQ2) NOI16330 $ +(-2.D0+4.D0*COSIQ+6.D0*COSQ2)/3.D0) NOI16340 F523=SINIQ*(4.92187512*SINI2*(-2.D0-4.D0*COSIQ+10.D0*COSQ2) NOI16350 $ +6.56250012D0*(1.D0+2.D0*COSIQ-3.D0*COSQ2)) NOI16360 F542=29.53125D0*SINIQ*(2.D0-8.D0*COSIQ+COSQ2*(-12.D0+8.D0*COSIQ NOI16370 $ +10.D0*COSQ2)) NOI16380 F543=29.53125D0*SINIQ*(-2.D0-8.D0*COSIQ+COSQ2*(12.D0+8.D0*COSIQ NOI16390 $ -10.D0*COSQ2)) NOI16400 XNO2=XNQ*XNQ NOI16410 AINV2=AQNV*AQNV NOI16420 TEMP1=3.D0*XNO2*AINV2 NOI16430 TEMP=TEMP1*ROOT22 NOI16440 D2201=TEMP*F220*G201 NOI16450 D2211=TEMP*F221*G221 NOI16460 TEMP1=TEMP1*AQNV NOI16470 TEMP=TEMP1*ROOT32 NOI16480 D3210=TEMP*F321*G310 NOI16490 D3222=TEMP*F322*G322 NOI16500 TEMP1=TEMP1*AQNV NOI16510 TEMP=2.D0*TEMP1*ROOT44 NOI16520 D4410=TEMP*F441*G410 NOI16530 D4422=TEMP*F442*G442 NOI16540 TEMP1=TEMP1*AQNV NOI16550 TEMP=TEMP1*ROOT52 NOI16560 D5220=TEMP*F522*G520 NOI16570 D5232=TEMP*F523*G532 NOI16580 TEMP=2.D0*TEMP1*ROOT54 NOI16590 D5421=TEMP*F542*G521 NOI16600 D5433=TEMP*F543*G533 NOI16610 XLAMO=XMAO+RAAN+RAAN-THETAG-THETAG NOI16620 BFACT=XLLDOT+XNODOT+XNODOT-THDT-THDT NOI16630 BFACT=BFACT+SSL+SSH+SSH NOI16640 GO TO 80 NOI16650 C NOI16660 C Syncronous resonance terms initialization NOI16670 C NOI16680 70 IRESFL=1 NOI16690 ISYNFL=1 NOI16700 G200=1.D0+EQSQ*(-2.5D0+0.8125D0*EQSQ) NOI16710 G310=1.D0+2.D0*EQSQ NOI16720 G300=1.D0+EQSQ*(-6.D0+6.60937D0*EQSQ) NOI16730 F220=0.75*(1.D0+COSIQ)*(1.D0+COSIQ) NOI16740 F311=0.9375D0*SINIQ*SINIQ*(1.D0+3.D0*COSIQ)-0.75D0*(1.D0+COSIQ) NOI16750 F330=1.D0+COSIQ NOI16760 F330=1.875D0+F330*F330*F330 NOI16770 DEL1=3.D0*XNQ*XNQ*AQNV*AQNV NOI16780 DEL2=2.D0*DEL1*F220*G200*Q22 NOI16790 DEL3=3.D0*DEL1*F330*G300*Q33*AQNV NOI16800 DEL1=DEL1*F311*G310*Q31*AQNV NOI16810 FASX2=0.13130908D0 NOI16820 FASX4=2.8843198D0 NOI16830 FASX6=0.37448087D0 NOI16840 XLAMO=XMAO+RAAN+OMPE-THETAG NOI16850 BFACT=XLLDOT+XPIDOT-THDT NOI16860 BFACT=BFACT+SSL+SSG+SSH NOI16870 80 XFACT=BFACT-XNQ NOI16880 C NOI16890 C Initialization integrator NOI16900 C NOI16910 XLI=XLAMO NOI16920 XNI=XNQ NOI16930 ATIME=0.D0 NOI16940 STEPP=720.D0 NOI16950 STEPN=-720.D0 NOI16960 STEP2=259200.D0 NOI16970 RETURN NOI16980 C NOI16990 END NOI17000 SUBROUTINE DPSEC (XLL,OMGASM,XNODES,EM,XINC,XN,TSINCE) NOI17010 C NOI17020 C----------------------------------------------------------------------+NOI17030 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI17040 C----------------------------------------------------------------------|NOI17050 C Title : DPSEC |NOI17060 C----------------------------------------------------------------------|NOI17070 C Function : Entrance for deep space secular effects NOI17080 C----------------------------------------------------------------------|NOI17090 C Author : NORAD |NOI17100 C----------------------------------------------------------------------|NOI17110 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI17120 C----------------------------------------------------------------------|NOI17130 C Date-coded : 11 October 1988 |NOI17140 C----------------------------------------------------------------------|NOI17150 C Last update : 14 December 1988 |NOI17160 C----------------------------------------------------------------------|NOI17170 C Notes : None |NOI17180 C----------------------------------------------------------------------|NOI17190 C Input data : |NOI17200 C----------------------------------------------------------------------|NOI17210 C Output data : |NOI17220 C----------------------------------------------------------------------+NOI17230 C NOI17240 IMPLICIT REAL*8 (A-H,O-Z) NOI17250 COMMON /CM0/ XNDT2O,XNDD6O,BSTAR,CLIN,RAAN,ECC,OMPE,TETAE,XNO NOI17260 DATA G22,G32/5.7636396D0,.95240898D0/,THDT/4.3752691D-3/ NOI17270 DATA G44,G52/1.8014998D0,1.0508330D0/,G54/4.4108898/ NOI17280 C NOI17290 XLL=XLL+SSL*TSINCE NOI17300 OMGASM=OMGASM+SSG*TSINCE NOI17310 XNODES=XNODES+SSH*TSINCE NOI17320 EM=ECC+SSE*TSINCE NOI17330 XINC=CLIN+SSI*TSINCE NOI17340 IF(IRESFL.EQ.0) RETURN NOI17350 C NOI17360 1 IF(ATIME.EQ.0.D0) GO TO 11 NOI17370 IF(TSINCE.GE.0.D0.AND.ATIME.LT.0.D0) GO TO 11 NOI17380 IF(TSINCE.LT.0.D0.AND.ATIME.GE.0.D0) GO TO 11 NOI17390 IF(DABS(TSINCE).GE.DABS(ATIME)) GO TO 2 NOI17400 DELT=STEPP NOI17410 IF(TSINCE.GE.0.D0) DELT=STEPN NOI17420 ASSIGN 1 TO IRET NOI17430 GO TO 9 NOI17440 C NOI17450 2 DELT=STEPN NOI17460 IF(TSINCE.GT.0.D0) DELT=STEPP NOI17470 3 IF(DABS(TSINCE-ATIME).LT.STEPP) GO TO 4 NOI17480 ASSIGN 3 TO IRET NOI17490 GO TO 9 NOI17500 C NOI17510 4 FT=TSINCE-ATIME NOI17520 ASSIGN 5 TO IRETN NOI17530 GO TO 6 NOI17540 C NOI17550 5 XN=XNI+XNDOT*FT+XNDDT*FT*FT*0.5D0 NOI17560 XL=XLI+XLDOT*FT+XNDOT*FT*FT*0.5D0 NOI17570 TEMP=-XNODES+THETAG+TSINCE*THDT NOI17580 XLL=XL-OMGASM+TEMP NOI17590 IF(ISYNFL.EQ.0) XLL=XL+TEMP+TEMP NOI17600 RETURN NOI17610 C NOI17620 C Dot terms calculated NOI17630 C NOI17640 6 IF(ISYNFL.EQ.0) GO TO 7 NOI17650 XNDOT=DEL1*DSIN(XLI-FASX2)+DEL2*DSIN(2.D0*(XLI-FASX4)) NOI17660 $ +DEL3*DSIN(3.D0*(XLI-FASX6)) NOI17670 XNDDT=DEL1*DCOS(XLI-FASX2)+2.D0*DEL2*DCOS(2.D0*(XLI-FASX4)) NOI17680 $ +3.D0*DEL3*DCOS(3.D0*(XLI-FASX6)) NOI17690 GO TO 8 NOI17700 C NOI17710 7 XOMI=OMEGAQ+OMGDT*ATIME NOI17720 X2OMI=XOMI+XOMI NOI17730 X2LI=XLI+XLI NOI17740 XNDOT=D2201*DSIN(X2OMI+XLI-G22)+D2211*DSIN(XLI-G22) NOI17750 $ +D3210*DSIN(XOMI+XLI-G32)+D3222*DSIN(-XOMI+XLI-G52) NOI17760 $ +D4410*DSIN(X2OMI+X2LI-G44)+D4422*DSIN(X2LI-G44) NOI17770 $ +D5220*DSIN(XOMI+XLI-G52)+D5232*DSIN(-XOMI+XLI-G52) NOI17780 $ +D5421*DSIN(XOMI+X2LI-G54)+D5433*DSIN(-XOMI+X2LI-G54) NOI17790 XNDDT=D2201*DCOS(X2OMI+XLI-G22)+D2211*DCOS(XLI-G22) NOI17800 $ +D3210*DCOS(XOMI+XLI-G32)+D3222*DCOS(-XOMI+XLI-G32) NOI17810 $ +D5220*DCOS(XOMI+XLI-G52)+D5232*DCOS(-XOMI+XLI-G52) NOI17820 $ +2.D0*(D4410*DCOS(X2OMI+X2LI-G44)+D4422*DCOS(X2LI-G44) NOI17830 $ +D5421*DCOS(XOMI+X2LI-G54)+D5433*DCOS(-XOMI+X2LI-G54)) NOI17840 C NOI17850 8 XLDOT=XNI+XFACT NOI17860 XNDDT=XNDDT+XLDOT NOI17870 GO TO IRETN,(5,10) NOI17880 C NOI17890 C Integrator NOI17900 C NOI17910 9 ASSIGN 10 TO IRETN NOI17920 GO TO 6 NOI17930 C NOI17940 10 XLI=XLI+XLDOT*DELT+XNDOT*STEP2 NOI17950 XNI=XNI+XNDOT*DELT+XNDDT*STEP2 NOI17960 ATIME=ATIME+DELT NOI17970 GO TO IRET,(1,3) NOI17980 C NOI17990 C Epoch restart NOI18000 C NOI18010 11 DELT=STEPN NOI18020 IF(TSINCE.GE.0.D0) DELT=STEPP NOI18030 ATIME=0.D0 NOI18040 XNI=XNQ NOI18050 XLI=XLAMO NOI18060 GO TO 3 NOI18070 C NOI18080 END NOI18090 SUBROUTINE DPPER (EM,XINC,OMGASM,XNODES,XLL) NOI18100 C NOI18110 C----------------------------------------------------------------------+NOI18120 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI18130 C----------------------------------------------------------------------|NOI18140 C Title : DPPER |NOI18150 C----------------------------------------------------------------------|NOI18160 C Function : Entrances for lunar-solar periodics NOI18170 C----------------------------------------------------------------------|NOI18180 C Author : NORAD |NOI18190 C----------------------------------------------------------------------|NOI18200 C Programmer : NORAD (Modified by Daniele Mortari for SAN MARCO D/M) |NOI18210 C----------------------------------------------------------------------|NOI18220 C Date-coded : 11 October 1988 |NOI18230 C----------------------------------------------------------------------|NOI18240 C Last update : 14 December 1988 |NOI18250 C----------------------------------------------------------------------|NOI18260 C Notes : None |NOI18270 C----------------------------------------------------------------------|NOI18280 C Input data : |NOI18290 C----------------------------------------------------------------------|NOI18300 C Output data : |NOI18310 C----------------------------------------------------------------------+NOI18320 C NOI18330 IMPLICIT REAL*8 (A-H,O-Z) NOI18340 COMMON /CONST/ PIG,CW,REQ NOI18350 DATA ZNS,ZES/1.19459D-5,0.01675D0/,ZNL,ZEL/1.5835218D-4,0.05490D0/NOI18360 C NOI18370 SINIS=DSIN(XINC) NOI18380 COSIS=DCOS(XINC) NOI18390 IF(DABS(SAVTSN-T).LT.30.D0) GO TO 1 NOI18400 SAVTSN=T NOI18410 ZM=ZMOS+ZNS*T NOI18420 ZF=ZM+2.D0*ZES*DSIN(ZM) NOI18430 SINZF=DSIN(ZF) NOI18440 F2=0.5D0*SINZF*SINZF-0.25D0 NOI18450 F3=-0.5D0*SINZF*DCOS(ZF) NOI18460 SES=SE2*F2+SE3*F3 NOI18470 SIS=SI2*F2+SI3*F3 NOI18480 SLS=SL2*F2+SL3*F3+SL4*SINZF NOI18490 SGHS=SGH2*F2+SGH3*F3+SGH4*SINZF NOI18500 SHS=SH2*F2+SH3*F3 NOI18510 ZM=ZMOL+ZNL*T NOI18520 ZF=ZM+2.D0*ZEL*DSIN(ZM) NOI18530 SINZF=DSIN(ZF) NOI18540 F2=0.5D0*SINZF*SINZF-0.25D0 NOI18550 F3=-0.5D0*SINZF*DCOS(ZF) NOI18560 SEL=EE2*F2+E3*F3 NOI18570 SIL=XI2*F2+XI3*F3 NOI18580 SLL=XL2*F2+XL3*F3+XL4*SINZF NOI18590 SGHL=XGH2*F2+XGH3*F3+XGH4*SINZF NOI18600 SHL=XH2*F2+XH3*F3 NOI18610 PE=SES+SEL NOI18620 PINC=SIS+SIL NOI18630 PL=SLS+SLL NOI18640 C NOI18650 1 PGH=SGHS+SGHL NOI18660 PH=SHS+SHL NOI18670 XINC=XINC+PINC NOI18680 EM=EM+PE NOI18690 IF(XQNCL.LT.0.2D0) GO TO 2 NOI18700 C NOI18710 C Apply periodics directly NOI18720 C NOI18730 PH=PH/SINIQ NOI18740 PGH=PGH-COSIQ*PH NOI18750 OMGASM=OMGASM+PGH NOI18760 XNODES=XNODES+PH NOI18770 XLL=XLL+PL NOI18780 RETURN NOI18790 C NOI18800 C Apply periodics with lyddane modification NOI18810 C NOI18820 2 SINOK=DSIN(XNODES) NOI18830 COSOK=DCOS(XNODES) NOI18840 ALFDP=SINIS*SINOK NOI18850 BETDP=SINIS*COSOK NOI18860 DALF=PH*COSOK+PINC*COSIS*SINOK NOI18870 DBET=-PH*SINOK+PINC*COSIS*COSOK NOI18880 ALFDP=ALFDP+DALF NOI18890 BETDP=BETDP+DBET NOI18900 XLS=XLL+OMGASM+COSIS*XNODES NOI18910 DLS=PL+PGH-PINC*XNODES*SINIS NOI18920 XLS=XLS+DLS NOI18930 CALL FTAN (ALFDP,BETDP,XNODES) NOI18940 XNODES=XNODES/CW NOI18950 XLL=XLL+PL NOI18960 OMGASM=XLS-XLL-DCOS(XINC)*XNODES NOI18970 RETURN NOI18980 C NOI18990 END NOI19000 DOUBLE PRECISION FUNCTION FMOD2P (ANGLE) NOI19010 C NOI19020 C----------------------------------------------------------------------+NOI19030 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI19040 C----------------------------------------------------------------------|NOI19050 C Title : FMOD2P |NOI19060 C----------------------------------------------------------------------|NOI19070 C Function : Brings the angle between 0 and 2*PI (Radiants) |NOI19080 C----------------------------------------------------------------------|NOI19090 C Author : Daniele Mortari |NOI19100 C----------------------------------------------------------------------|NOI19110 C Programmer : Daniele Mortari |NOI19120 C----------------------------------------------------------------------|NOI19130 C Date-coded : 11 October 1988 |NOI19140 C----------------------------------------------------------------------|NOI19150 C Last update : 14 December 1988 |NOI19160 C----------------------------------------------------------------------|NOI19170 C Notes : None |NOI19180 C----------------------------------------------------------------------|NOI19190 C Input data : |NOI19200 C----------------------------------------------------------------------|NOI19210 C Output data : |NOI19220 C----------------------------------------------------------------------+NOI19230 C NOI19240 IMPLICIT REAL*8(A-H,O-Z) NOI19250 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI19260 C NOI19270 FMOD2P=ANGLE-IDINT(ANGLE/TWOPI)*TWOPI NOI19280 IF(FMOD2P.LT.0.D0) FMOD2P=FMOD2P+TWOPI NOI19290 RETURN NOI19300 C NOI19310 END NOI19320 SUBROUTINE CICCIO (THETAG,DS50,EPOCH) NOI19330 C NOI19340 C----------------------------------------------------------------------+NOI19350 C Identification ===> SAN MARCO PROJECT - ATTITUDE SOFTWARE - X |NOI19360 C----------------------------------------------------------------------|NOI19370 C Title : THETAG |NOI19380 C----------------------------------------------------------------------|NOI19390 C Function : |NOI19400 C----------------------------------------------------------------------|NOI19410 C Author : Daniele Mortari |NOI19420 C----------------------------------------------------------------------|NOI19430 C Programmer : Daniele Mortari |NOI19440 C----------------------------------------------------------------------|NOI19450 C Date-coded : 11 October 1988 |NOI19460 C----------------------------------------------------------------------|NOI19470 C Last update : 14 December 1988 |NOI19480 C----------------------------------------------------------------------|NOI19490 C Notes : None |NOI19500 C----------------------------------------------------------------------|NOI19510 C Input data : |NOI19520 C----------------------------------------------------------------------|NOI19530 C Output data : |NOI19540 C----------------------------------------------------------------------+NOI19550 C NOI19560 IMPLICIT REAL*8(A-H,O-Z) NOI19570 COMMON /CM1/ CK2,CK4,ZERO,QOMS2T,S,TOTHRD,XJ3,XKE,TWOPI,AE NOI19580 C NOI19590 YR=(EPOCH+2.D-7)*1.D-3 NOI19600 JY=YR NOI19610 YR=JY NOI19620 D=EPOCH-YR*1.D3 NOI19630 IF(JY.LT.10) JY=JY+80 NOI19640 N=(JY-69)/4 NOI19650 IF(JY.LT.70) N=(JY-72)/4 NOI19660 DS50=7305.D0+365.D0*(JY-70)+N+D NOI19670 THETA=1.72944494D0+6.3003880987D0*DS50 NOI19680 THETAG=FMOD2P(THETA) NOI19690 RETURN NOI19700 C NOI19710 END NOI19720