diff options
Diffstat (limited to 'src/slalib/pertue.f')
| -rw-r--r-- | src/slalib/pertue.f | 535 |
1 files changed, 535 insertions, 0 deletions
diff --git a/src/slalib/pertue.f b/src/slalib/pertue.f new file mode 100644 index 0000000..d6f6f53 --- /dev/null +++ b/src/slalib/pertue.f @@ -0,0 +1,535 @@ + SUBROUTINE sla_PERTUE (DATE, U, JSTAT) +*+ +* - - - - - - - +* P E R T U E +* - - - - - - - +* +* Update the universal elements of an asteroid or comet by applying +* planetary perturbations. +* +* Given: +* DATE d final epoch (TT MJD) for the updated elements +* +* Given and returned: +* U d(13) universal elements (updated in place) +* +* (1) combined mass (M+m) +* (2) total energy of the orbit (alpha) +* (3) reference (osculating) epoch (t0) +* (4-6) position at reference epoch (r0) +* (7-9) velocity at reference epoch (v0) +* (10) heliocentric distance at reference epoch +* (11) r0.v0 +* (12) date (t) +* (13) universal eccentric anomaly (psi) of date, approx +* +* Returned: +* JSTAT i status: +* +102 = warning, distant epoch +* +101 = warning, large timespan ( > 100 years) +* +1 to +8 = coincident with major planet (Note 5) +* 0 = OK +* -1 = numerical error +* +* Called: sla_PLANET, sla_UE2PV, sla_PV2UE +* +* Notes: +* +* 1 The "universal" elements are those which define the orbit for the +* purposes of the method of universal variables (see reference 2). +* They consist of the combined mass of the two bodies, an epoch, +* and the position and velocity vectors (arbitrary reference frame) +* at that epoch. The parameter set used here includes also various +* quantities that can, in fact, be derived from the other +* information. This approach is taken to avoiding unnecessary +* computation and loss of accuracy. The supplementary quantities +* are (i) alpha, which is proportional to the total energy of the +* orbit, (ii) the heliocentric distance at epoch, (iii) the +* outwards component of the velocity at the given epoch, (iv) an +* estimate of psi, the "universal eccentric anomaly" at a given +* date and (v) that date. +* +* 2 The universal elements are with respect to the J2000 equator and +* equinox. +* +* 3 The epochs DATE, U(3) and U(12) are all Modified Julian Dates +* (JD-2400000.5). +* +* 4 The algorithm is a simplified form of Encke's method. It takes as +* a basis the unperturbed motion of the body, and numerically +* integrates the perturbing accelerations from the major planets. +* The expression used is essentially Sterne's 6.7-2 (reference 1). +* Everhart and Pitkin (reference 2) suggest rectifying the orbit at +* each integration step by propagating the new perturbed position +* and velocity as the new universal variables. In the present +* routine the orbit is rectified less frequently than this, in order +* to gain a slight speed advantage. However, the rectification is +* done directly in terms of position and velocity, as suggested by +* Everhart and Pitkin, bypassing the use of conventional orbital +* elements. +* +* The f(q) part of the full Encke method is not used. The purpose +* of this part is to avoid subtracting two nearly equal quantities +* when calculating the "indirect member", which takes account of the +* small change in the Sun's attraction due to the slightly displaced +* position of the perturbed body. A simpler, direct calculation in +* double precision proves to be faster and not significantly less +* accurate. +* +* Apart from employing a variable timestep, and occasionally +* "rectifying the orbit" to keep the indirect member small, the +* integration is done in a fairly straightforward way. The +* acceleration estimated for the middle of the timestep is assumed +* to apply throughout that timestep; it is also used in the +* extrapolation of the perturbations to the middle of the next +* timestep, to predict the new disturbed position. There is no +* iteration within a timestep. +* +* Measures are taken to reach a compromise between execution time +* and accuracy. The starting-point is the goal of achieving +* arcsecond accuracy for ordinary minor planets over a ten-year +* timespan. This goal dictates how large the timesteps can be, +* which in turn dictates how frequently the unperturbed motion has +* to be recalculated from the osculating elements. +* +* Within predetermined limits, the timestep for the numerical +* integration is varied in length in inverse proportion to the +* magnitude of the net acceleration on the body from the major +* planets. +* +* The numerical integration requires estimates of the major-planet +* motions. Approximate positions for the major planets (Pluto +* alone is omitted) are obtained from the routine sla_PLANET. Two +* levels of interpolation are used, to enhance speed without +* significantly degrading accuracy. At a low frequency, the routine +* sla_PLANET is called to generate updated position+velocity "state +* vectors". The only task remaining to be carried out at the full +* frequency (i.e. at each integration step) is to use the state +* vectors to extrapolate the planetary positions. In place of a +* strictly linear extrapolation, some allowance is made for the +* curvature of the orbit by scaling back the radius vector as the +* linear extrapolation goes off at a tangent. +* +* Various other approximations are made. For example, perturbations +* by Pluto and the minor planets are neglected, relativistic effects +* are not taken into account and the Earth-Moon system is treated as +* a single body. +* +* In the interests of simplicity, the background calculations for +* the major planets are carried out en masse. The mean elements and +* state vectors for all the planets are refreshed at the same time, +* without regard for orbit curvature, mass or proximity. +* +* 5 This routine is not intended to be used for major planets. +* However, if major-planet elements are supplied, sensible results +* will, in fact, be produced. This happens because the routine +* checks the separation between the body and each of the planets and +* interprets a suspiciously small value (0.001 AU) as an attempt to +* apply the routine to the planet concerned. If this condition is +* detected, the contribution from that planet is ignored, and the +* status is set to the planet number (Mercury=1,...,Neptune=8) as a +* warning. +* +* References: +* +* 1 Sterne, Theodore E., "An Introduction to Celestial Mechanics", +* Interscience Publishers Inc., 1960. Section 6.7, p199. +* +* 2 Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983. +* +* P.T.Wallace Starlink 18 March 1999 +* +* Copyright (C) 1999 Rutherford Appleton Laboratory +*- + + IMPLICIT NONE + DOUBLE PRECISION DATE,U(13) + INTEGER JSTAT + +* Coefficient relating timestep to perturbing force + DOUBLE PRECISION TSC + PARAMETER (TSC=1D-4) + +* Minimum and maximum timestep (days) + DOUBLE PRECISION TSMIN,TSMAX + PARAMETER (TSMIN=0.01D0,TSMAX=10D0) + +* Age limit for major-planet state vector (days) + DOUBLE PRECISION AGEPMO + PARAMETER (AGEPMO=5D0) + +* Age limit for major-planet mean elements (days) + DOUBLE PRECISION AGEPEL + PARAMETER (AGEPEL=50D0) + +* Margin for error when deciding whether to renew the planetary data + DOUBLE PRECISION TINY + PARAMETER (TINY=1D-6) + +* Age limit for the body's osculating elements (before rectification) + DOUBLE PRECISION AGEBEL + PARAMETER (AGEBEL=100D0) + +* Gaussian gravitational constant (exact) and square + DOUBLE PRECISION GCON,GCON2 + PARAMETER (GCON=0.01720209895D0,GCON2=GCON*GCON) + +* The final epoch + DOUBLE PRECISION TFINAL + +* The body's current universal elements + DOUBLE PRECISION UL(13) + +* Current reference epoch + DOUBLE PRECISION T0 + +* Timespan from latest orbit rectification to final epoch (days) + DOUBLE PRECISION TSPAN + +* Time left to go before integration is complete + DOUBLE PRECISION TLEFT + +* Time direction flag: +1=forwards, -1=backwards + DOUBLE PRECISION FB + +* First-time flag + LOGICAL FIRST + +* +* The current perturbations +* +* Epoch (days relative to current reference epoch) + DOUBLE PRECISION RTN +* Position (AU) + DOUBLE PRECISION PERP(3) +* Velocity (AU/d) + DOUBLE PRECISION PERV(3) +* Acceleration (AU/d/d) + DOUBLE PRECISION PERA(3) +* + +* Length of current timestep (days), and half that + DOUBLE PRECISION TS,HTS + +* Epoch of middle of timestep + DOUBLE PRECISION T + +* Epoch of planetary mean elements + DOUBLE PRECISION TPEL + +* Planet number (1=Mercury, 2=Venus, 3=EMB...8=Neptune) + INTEGER NP + +* Planetary universal orbital elements + DOUBLE PRECISION UP(13,8) + +* Epoch of planetary state vectors + DOUBLE PRECISION TPMO + +* State vectors for the major planets (AU,AU/s) + DOUBLE PRECISION PVIN(6,8) + +* +* Correction terms for extrapolated major planet vectors +* +* Sun-to-planet distances squared multiplied by 3 + DOUBLE PRECISION R2X3(8) +* Sunward acceleration terms, G/2R^3 + DOUBLE PRECISION GC(8) +* Tangential-to-circular correction factor + DOUBLE PRECISION FC +* Radial correction factor due to Sunwards acceleration + DOUBLE PRECISION FG +* + +* The body's unperturbed and perturbed state vectors (AU,AU/s) + DOUBLE PRECISION PV0(6),PV(6) + +* The body's perturbed and unperturbed heliocentric distances (AU) cubed + DOUBLE PRECISION R03,R3 + +* The perturbating accelerations, indirect and direct + DOUBLE PRECISION FI(3),FD(3) + +* Sun-to-planet vector, and distance cubed + DOUBLE PRECISION RHO(3),RHO3 + +* Body-to-planet vector, and distance cubed + DOUBLE PRECISION DELTA(3),DELTA3 + +* Miscellaneous + INTEGER I,J + DOUBLE PRECISION R2,W,DT,DT2,FT + +* Planetary inverse masses, Mercury through Neptune + DOUBLE PRECISION AMAS(8) + DATA AMAS / 6023600D0,408523.5D0,328900.5D0,3098710D0, + : 1047.355D0,3498.5D0,22869D0,19314D0 / + +*---------------------------------------------------------------------* + + +* Preset the status to OK. + JSTAT = 0 + +* Copy the final epoch. + TFINAL = DATE + +* Copy the elements (which will be periodically updated). + DO I=1,13 + UL(I) = U(I) + END DO + +* Initialize the working reference epoch. + T0=UL(3) + +* Total timespan (days) and hence time left. + TSPAN = TFINAL-T0 + TLEFT = TSPAN + +* Warn if excessive. + IF (ABS(TSPAN).GT.36525D0) JSTAT=101 + +* Time direction: +1 for forwards, -1 for backwards. + FB = SIGN(1D0,TSPAN) + +* Initialize relative epoch for start of current timestep. + RTN = 0D0 + +* Reset the perturbations (position, velocity, acceleration). + DO I=1,3 + PERP(I) = 0D0 + PERV(I) = 0D0 + PERA(I) = 0D0 + END DO + +* Set "first iteration" flag. + FIRST = .TRUE. + +* Step through the time left. + DO WHILE (FB*TLEFT.GT.0D0) + +* Magnitude of current acceleration due to planetary attractions. + IF (FIRST) THEN + TS = TSMIN + ELSE + R2 = 0D0 + DO I=1,3 + W = FD(I) + R2 = R2+W*W + END DO + W = SQRT(R2) + +* Use the acceleration to decide how big a timestep can be tolerated. + IF (W.NE.0D0) THEN + TS = MIN(TSMAX,MAX(TSMIN,TSC/W)) + ELSE + TS = TSMAX + END IF + END IF + TS = TS*FB + +* Override if final epoch is imminent. + TLEFT = TSPAN-RTN + IF (ABS(TS).GT.ABS(TLEFT)) TS=TLEFT + +* Epoch of middle of timestep. + HTS = TS/2D0 + T = T0+RTN+HTS + +* Is it time to recompute the major-planet elements? + IF (FIRST.OR.ABS(T-TPEL)-AGEPEL.GE.TINY) THEN + +* Yes: go forward in time by just under the maximum allowed. + TPEL = T+FB*AGEPEL + +* Compute the state vector for the new epoch. + DO NP=1,8 + CALL sla_PLANET(TPEL,NP,PV,J) + +* Warning if remote epoch, abort if error. + IF (J.EQ.1) THEN + JSTAT = 102 + ELSE IF (J.NE.0) THEN + GO TO 9010 + END IF + +* Transform the vector into universal elements. + CALL sla_PV2UE(PV,TPEL,0D0,UP(1,NP),J) + IF (J.NE.0) GO TO 9010 + END DO + END IF + +* Is it time to recompute the major-planet motions? + IF (FIRST.OR.ABS(T-TPMO)-AGEPMO.GE.TINY) THEN + +* Yes: look ahead. + TPMO = T+FB*AGEPMO + +* Compute the motions of each planet (AU,AU/d). + DO NP=1,8 + +* The planet's position and velocity (AU,AU/s). + CALL sla_UE2PV(TPMO,UP(1,NP),PVIN(1,NP),J) + IF (J.NE.0) GO TO 9010 + +* Scale velocity to AU/d. + DO J=4,6 + PVIN(J,NP) = PVIN(J,NP)*86400D0 + END DO + +* Precompute also the extrapolation correction terms. + R2 = 0D0 + DO I=1,3 + W = PVIN(I,NP) + R2 = R2+W*W + END DO + R2X3(NP) = R2*3D0 + GC(NP) = GCON2/(2D0*R2*SQRT(R2)) + END DO + END IF + +* Reset the first-time flag. + FIRST = .FALSE. + +* Unperturbed motion of the body at middle of timestep (AU,AU/s). + CALL sla_UE2PV(T,UL,PV0,J) + IF (J.NE.0) GO TO 9010 + +* Perturbed position of the body (AU) and heliocentric distance cubed. + R2 = 0D0 + DO I=1,3 + W = PV0(I)+PERP(I)+(PERV(I)+PERA(I)*HTS/2D0)*HTS + PV(I) = W + R2 = R2+W*W + END DO + R3 = R2*SQRT(R2) + +* The body's unperturbed heliocentric distance cubed. + R2 = 0D0 + DO I=1,3 + W = PV0(I) + R2 = R2+W*W + END DO + R03 = R2*SQRT(R2) + +* Compute indirect and initialize direct parts of the perturbation. + DO I=1,3 + FI(I) = PV0(I)/R03-PV(I)/R3 + FD(I) = 0D0 + END DO + +* Ready to compute the direct planetary effects. + +* Interval from state-vector epoch to middle of current timestep. + DT = T-TPMO + DT2 = DT*DT + +* Planet by planet. + DO NP=1,8 + +* First compute the extrapolation in longitude (squared). + R2 = 0D0 + DO J=4,6 + W = PVIN(J,NP)*DT + R2 = R2+W*W + END DO + +* Hence the tangential-to-circular correction factor. + FC = 1D0+R2/R2X3(NP) + +* The radial correction factor due to the inwards acceleration. + FG = 1D0-GC(NP)*DT2 + +* Planet's position, and heliocentric distance cubed. + R2 = 0D0 + DO I=1,3 + W = FG*(PVIN(I,NP)+FC*PVIN(I+3,NP)*DT) + RHO(I) = W + R2 = R2+W*W + END DO + RHO3 = R2*SQRT(R2) + +* Body-to-planet vector, light-time, and distance cubed. + R2 = 0D0 + DO I=1,3 + W = RHO(I)-PV(I) + DELTA(I) = W + R2 = R2+W*W + END DO + W = SQRT(R2) + DELTA3 = R2*SQRT(R2) + +* If too close, ignore this planet and set a warning. + IF (R2.LT.1D-6) THEN + JSTAT = NP + ELSE + +* Accumulate "direct" part of perturbation acceleration. + W = AMAS(NP) + DO I=1,3 + FD(I) = FD(I)+(DELTA(I)/DELTA3-RHO(I)/RHO3)/W + END DO + END IF + END DO + +* Update the perturbations to the end of the timestep. + RTN = RTN+TS + DO I=1,3 + W = (FI(I)+FD(I))*GCON2 + FT = W*TS + PERP(I) = PERP(I)+(PERV(I)+FT/2D0)*TS + PERV(I) = PERV(I)+FT + PERA(I) = W + END DO + +* Time still to go. + TLEFT = TSPAN-RTN + +* Is it either time to rectify the orbit or the last time through? + IF (ABS(RTN).GE.AGEBEL.OR.FB*TLEFT.LE.0D0) THEN + +* Yes: update to the end of the current timestep. + T0 = T0+RTN + RTN = 0D0 + +* The body's unperturbed motion (AU,AU/s). + CALL sla_UE2PV(T0,UL,PV0,J) + IF (J.NE.0) GO TO 9010 + +* Add and re-initialize the perturbations. + DO I=1,3 + J = I+3 + PV(I) = PV0(I)+PERP(I) + PV(J) = PV0(J)+PERV(I)/86400D0 + PERP(I) = 0D0 + PERV(I) = 0D0 + PERA(I) = FD(I)*GCON2 + END DO + +* Use the position and velocity to set up new universal elements. + CALL sla_PV2UE(PV,T0,0D0,UL,J) + IF (J.NE.0) GO TO 9010 + +* Adjust the timespan and time left. + TSPAN = TFINAL-T0 + TLEFT = TSPAN + END IF + +* Next timestep. + END DO + +* Return the updated universal-element set. + DO I=1,13 + U(I) = UL(I) + END DO + +* Finished. + GO TO 9999 + +* Miscellaneous numerical error. + 9010 CONTINUE + JSTAT = -1 + + 9999 CONTINUE + END |