eddy/BTA_lib
1
0
Fork 0
mirror of https://github.com/eddyem/BTA_lib.git synced 2025-12-06 18:55:14 +03:00
Code Issues Packages Projects Releases Wiki Activity
BTA_lib/slalib/pertue.f
eddyem 6b3ca8a34b copy
2015-10-27 19:46:19 +03:00

643 lines
20 KiB
Fortran
Raw Blame History

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 +10 = coincident with major planet (Note 5)
* 0 = OK
* -1 = numerical error
*
* Called: sla_EPJ, sla_PLANET, sla_PV2UE, sla_UE2PV, sla_EPV,
* sla_PREC, sla_DMOON, sla_DMXV
*
* 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 and relativistic
* effects are not taken into account.
*
* 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.
*
* The Earth-Moon system is treated as a single body when the body is
* distant but as separate bodies when closer to the EMB than the
* parameter RNE, which incurs a time penalty but improves accuracy
* for near-Earth objects.
*
* 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 (1-10 = Mercury, Venus, EMB,
* Mars, Jupiter, Saturn, Uranus, Neptune, Earth, Moon) 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.
*
* Last revision: 27 December 2004
*
* Copyright P.T.Wallace. All rights reserved.
*
* License:
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program (see SLA_CONDITIONS); if not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330,
* Boston, MA 02111-1307 USA
*
*-
IMPLICIT NONE
DOUBLE PRECISION DATE,U(13)
INTEGER JSTAT
* Distance from EMB at which Earth and Moon are treated separately
DOUBLE PRECISION RNE
PARAMETER (RNE=1D0)
* Coincidence with major planet distance
DOUBLE PRECISION COINC
PARAMETER (COINC=0.0001D0)
* 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)
* Earth velocity and position vectors (AU,AU/s)
DOUBLE PRECISION VB(3),PB(3),VH(3),PE(3)
* Moon geocentric state vector (AU,AU/s) and position part
DOUBLE PRECISION PVM(6),PM(3)
* Date to J2000 de-precession matrix
DOUBLE PRECISION PMAT(3,3)
*
* 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,R,FT
LOGICAL NE
DOUBLE PRECISION sla_EPJ
* Planetary inverse masses, Mercury through Neptune then Earth and Moon
DOUBLE PRECISION AMAS(10)
DATA AMAS / 6023600D0, 408523.5D0, 328900.5D0, 3098710D0,
: 1047.355D0, 3498.5D0, 22869D0, 19314D0,
: 332946.038D0, 27068709D0 /
*
* License:
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program (see SLA_CONDITIONS); if not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330,
* Boston, MA 02111-1307 USA
*
*----------------------------------------------------------------------*
* 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.
* Reset the "near-Earth" flag.
NE = .FALSE.
* Interval from state-vector epoch to middle of current timestep.
DT = T-TPMO
DT2 = DT*DT
* Planet by planet, including separate Earth and Moon.
DO NP=1,10
* Which perturbing body?
IF (NP.LE.8) THEN
* Planet: 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.
DO I=1,3
RHO(I) = FG*(PVIN(I,NP)+FC*PVIN(I+3,NP)*DT)
END DO
ELSE IF (NE) THEN
* Near-Earth and either Earth or Moon.
IF (NP.EQ.9) THEN
* Earth: position.
CALL sla_EPV(T,PE,VH,PB,VB)
DO I=1,3
RHO(I) = PE(I)
END DO
ELSE
* Moon: position.
CALL sla_PREC(sla_EPJ(T),2000D0,PMAT)
CALL sla_DMOON(T,PVM)
CALL sla_DMXV(PMAT,PVM,PM)
DO I=1,3
RHO(I) = PM(I)+PE(I)
END DO
END IF
END IF
* Proceed unless Earth or Moon and not the near-Earth case.
IF (NP.LE.8.OR.NE) THEN
* Heliocentric distance cubed.
R2 = 0D0
DO I=1,3
W = RHO(I)
R2 = R2+W*W
END DO
R = SQRT(R2)
RHO3 = R2*R
* Body-to-planet vector, and distance.
R2 = 0D0
DO I=1,3
W = RHO(I)-PV(I)
DELTA(I) = W
R2 = R2+W*W
END DO
R = SQRT(R2)
* If this is the EMB, set the near-Earth flag appropriately.
IF (NP.EQ.3.AND.R.LT.RNE) NE = .TRUE.
* Proceed unless EMB and this is the near-Earth case.
IF (.NOT.(NE.AND.NP.EQ.3)) THEN
* If too close, ignore this planet and set a warning.
IF (R.LT.COINC) THEN
JSTAT = NP
ELSE
* Accumulate "direct" part of perturbation acceleration.
DELTA3 = R2*R
W = AMAS(NP)
DO I=1,3
FD(I) = FD(I)+(DELTA(I)/DELTA3-RHO(I)/RHO3)/W
END DO
END IF
END IF
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
Reference in New Issue View Git Blame Copy Permalink