Talk:Help with Fortran90
From Stellarium
Summary of Program Function: Part1
Summary of Program Function: Part1 <sdwizard at users.sourceforge.net>
Summary from the internal comments and code:
Ephemeride files Nsa8-fs.dat (sequential) and Nsat8-fd.dat (direct) are available for the time span January 1, 1800 to January 1, 2050. The data is provided at 2 day intervals. The data file Nsat8-fd.dat provides the initial conditions for ILGS8_D/ILGS8_Q.
The times of interest are coded into the main routine (ILGS8_D.f/ILGS8_Q.f). Up to 150 times are allowed. The code locates the selected times within the intervals of the data file and integrates the system from the base of the interval to the desired point.
ILGS8_D/ILGS8_Q calls INS8_D/INS8_Q. INS8_D/INS8_Q computes initial time, step-size, final time, the initial positions and velocities of the system of satellites from the ephemeride "Nsat-fd.dat", and those of the perturbing bodies from the ephemeride JPL DE403 about the Neptunian center. INS8_D/INS8_Q calls ISAT8_D/ISAT8_Q, which handle the numerical integration of the three natural satellites, the other major planets of the solar system, and the Sun.
ILGS8_D
The user needs to enter the times of interest, as follows:
t(i) = value
where i = ranges from 1 to 150; represents the sequence of time values
nti = first sequence number of t(i)
ntf = last sequence number of t(i)
This main program calls
Call INSD8(nc1,nc2,nc3,nc4,nc5,nc6,nc7,nc8,
1 ipa,ico,nt,t,nti,ntf)
where nc1 ... nc8 are body identifiers for the moons of Neptune
ipa and ico are parameters affecting the output
nt = 150
t is a vector of length nt
nti and ntf are the initial and final time values in vector t
SUBROUTINE INSD8
This subroutine integrates 17 bodies (sun, 8 planets, 8 moons of Neptune) and displays the results for Triton, Nereid, Naiad, Thalassa, Despina, Galatea, Larissa, and Proteus. Integration can progress in a forward with respect to time or backward with respect to time fashion depending on the sequential values in the t vector.
A loop begins by extracting the first desired value of vector t. The code appears to skip each unrealistic value in vector t (outside of the range January 1, 1800 to January 1, 2050). Warnings are issued if this occurs.
File Nsat8-fd.dat is opened and the 2 day period that begins just at or before the selected time is located.
The data for that 2 day period is read. The initial time for the integration is set to the time value. The integration step size is set (positive for forward integration, negative for reverse integration, and zero when the desired time matches the data file time value). Matrix v0 is charged from Nsat8-fd.dat as follows: v0(1, body, 0) = record( 6*(body-1)+2) v0(2, body, 0) = record( 6*(body-1)+3) v0(3, body, 0) = record( 6*(body-1)+4) v0(1, body, 1) = record( 6*(body-1)+5) v0(2, body, 1) = record( 6*(body-1)+6) v0(3, body, 1) = record( 6*(body-1)+7)
My guess is that each record in Nsat8-fd.dat is set up as follows: Time; x position satellite 1; y position satellite 1; z position satellite 1;
x velocity satellite 1; y velocity satellite 1; z velocity satellite 1;
x position satellite 2; y position satellite 2; z position satellite 2;
x velocity satellite 2; y velocity satellite 2; z velocity satellite 2;
x position satellite 3; y position satellite 3; z position satellite 3;
x velocity satellite 3; y velocity satellite 3; z velocity satellite 3;
x position satellite 4; y position satellite 4; z position satellite 4;
x velocity satellite 4; y velocity satellite 4; z velocity satellite 4;
x position satellite 5; y position satellite 5; z position satellite 5;
x velocity satellite 5; y velocity satellite 5; z velocity satellite 5;
x position satellite 6; y position satellite 6; z position satellite 6;
x velocity satellite 6; y velocity satellite 6; z velocity satellite 6;
x position satellite 7; y position satellite 7; z position satellite 7;
x velocity satellite 7; y velocity satellite 7; z velocity satellite 7;
x position satellite 8; y position satellite 8; z position satellite 8;
x velocity satellite 8; y velocity satellite 8; z velocity satellite 8
where the satellite sequence is:
Triton; Nereid; Naiad; Thalassa; Despina; Galatea; Larissa; Proteus
The values in v0 are part of the initial conditions for the integration. Next, the subroutine JPL is called.
Call JPL(nc,nce,ET,v0)
where nc = 17 (total bodies)
nce = 8 (bodies to be integrated)
ET = initial time
If integration is needed:
Call ISATD8(tinit,painit,tfin,v0)
where tinit = the initial time
painit = the initial integration step size
tfin = the final time
SUBROUTINE JPL
JPL fills matrix v0. The regions filled are: v0(1, 9, 0) .. v0(1, 17, 0) v0(2, 9, 0) .. v0(2, 17, 0) v0(3, 9, 0) .. v0(3, 17, 0) v0(1, 9, 1) .. v0(1, 17, 1) v0(2, 9, 1) .. v0(2, 17, 1) v0(3, 9, 1) .. v0(3, 17, 1)
JPL gets the data from PLEPH_De403, which in turn calls STATEde403, which accesses the 403 data file.
Call PLEPH_De403 ( ET, NTARG, NCENT, RRD )
where NCENT = 8
RRD is a vector of data from the 403 file
NTARG = 1, 2, 13, 4, 5, 6, 7, 11, or 9
v0 element NTARG Body
9 1 = MERCURY
10 2 = VENUS
3 = EARTH
12 4 = MARS
13 5 = JUPITER
14 6 = SATURN
15 7 = URANUS
8 = NEPTUNE
17 9 = PLUTO
10 = MOON
16 11 = SUN
12 = SOLAR-SYSTEM BARYCENTER
11 13 = EARTH-MOON BARYCENTER
14 = NUTATIONS (LONGITUDE AND OBLIQ)
15 = LIBRATIONS, IF ON EPH FILE
At this time, I am not following the trail from PLEPH_De403.
SUBROUTINE INSATD8 description to follow
Summary of Program Function - Part 2
Summary of Program Function - Part 2 <sdwizard at users.sourceforge.net>
SUBROUTINE INSATD8
The code given below is modified to reflect parameter settings and internal logic. For example, code structured to handle output to files has been deleted because id and is are set to zero. So, for that matter, are code segments between triggered GOTO - CONTINUE structures. What remains is documented. Please do not extract and use this code for testing - It is my intention to supply code separately later if there is any interest.
c setting parameters
parameter(ndat=1) ! obtain a file "Ephemeride"
parameter(id=0) ! not recording directly on a file
parameter(is=0) ! not recording sequentially on a file
parameter(nc=17) ! number of bodies integrated
parameter(nd=21) ! number of derivatives
parameter(np=nc*(nc+1)/2) ! number of auxilliary variables
parameter(lb=6) ! inverse power of distance
parameter(nap=1) ! take flatness into account
parameter(npi=1) !
parameter(npf=8) ! first and last bodies with flatness
parameter(nid=1) !
parameter(nif=8) ! first and last bodies to print on screen
parameter(ned=1) !
parameter(nef=8) ! number of the first and last body to record
parameter(nce=nef-ned+1) ! number of successive bodies to record
parameter(nch=40) ! number of authorized step changes (?)
parameter(no=5) ! no>=1: an exact number of particular dates
parameter(nda=0) ! do not print the result at each date
c dimensioning matrices
dimension v0(3,nc,0:1),v(3,np,0:nd)
c vo = initial values built in earlier subroutines c v = derivative storage c coefficient and position storage start
dimension cr(nd,nd),cg(lb,nd,nd)
dimension r(np,0:nd),g(lb,np,0:nd)
dimension fx(lb,np,0:nd),fy(lb,np,0:nd),fz(lb,np,0:nd)
c coefficient and position storage end
dimension mas(0:nc),gmas(0:nc)
c mas = inverse of mass c gmas = ( gauss constant ) * ( gauss constant ) / mas
dimension gmu(nc,np)
c Le Guyader 2.2 Matrix
dimension aia(4)
c aia - 3 surfaces, energy
dimension t(0:nd)
c subroutine generated time steps = ( step size ) / (i!)
dimension ck(0:nd,0:nd)
dimension vx(no,nc,0:1),vy(no,nc,0:1),vz(no,nc,0:1)
c vx, vy, vz = position, speed by date c flatness matrices start
dimension pn(0:nd),ap(0:nd),dp(0:nd)
dimension cap(0:nd),cdp(0:nd),sap(0:nd),sdp(0:nd)
dimension cc(0:nd),cs(0:nd),sd(0:nd)
dimension acc(nc,0:nd),acs(nc,0:nd),asd(nc,0:nd)
dimension sa(4,nc,0:nd),psa(4,nc,0:nd),ppa(4,nc,0:nd)
dimension gcj(4,nc,0:nd)
dimension pcc(4,nc,0:nd),pcs(4,nc,0:nd),psd(4,nc,0:nd)
dimension plx(nc,0:nd),ply(nc,0:nd),plz(nc,0:nd)
c flatness matrices end c set constants
unan=1.d0
daou=2.d0
pi=3.141592653589793d0
dr=pi/180.d0 ! degrees in radians
mr=pi/10800.d0 ! minutes in radians
sr=pi/648000.d0 ! seconds in radians
rd=unan/dr ! radians in degrees
rm=unan/dmr ! radians in minutes
rs=unan/sr ! radians in seconds
dpir=daou*pi ! 2*pi
dpis=dpir*rs ! 2*pi in seconds
hj=1.d0/24.d0 ! 1 / ( hours in a day)
mj=1.d0/1440.d0 ! 1 / ( minutes in a day )
sj=1.d0/86400.d0 ! 1 / ( seconds in a day )
gk=0.01720209895d0 ! (IERS 1992) Gauss constant
gk2=gk*gk
UA=149597870.691d0 ! km (DE403 1995)
precr = 1.d-16 ! precision
ting = tfin-tinit ! integration interval
nbp = int((tfin-tinit)/painit) ! number of initial steps
api=abs(painit/1000.d0) ! if the step becomes too small
c flatness parameters
pn0=357.85d0*dr
ap0=299.36d0*dr
dp0= 43.46d0*dr
pn1=52.316d0*dr
ap1=0.7d0*dr
dp1= -0.51d0*dr
pn1=pn1/36525.d0 days are the unit of time
pj2= 3410.47359149002860d-6 ! J2 (Jacobson 1991)
pj4=-34.7030857830163906d-6 ! J4 (Jacobson 1991)
Rek=25225.d0 ! km (Jacobson 1991)
c inverses of masses (?), scaled masses, sum of scaled masses
mas(0)= 19416.29264728873d+00
mas(1)= 92944475.76880403d+00
mas(2)= 1.d+34
mas(3)= 1.d+34
mas(4)= 1.d+34
mas(5)= 1.d+34
mas(6)= 1.d+34
mas(7)= 1.d+34
mas(8)= 1.d+34
mas(9)= 6023600.000000000d+00
mas(10)= 408523.7100000000d+00
mas(11)= 328900.5603917904d+00
mas(12)= 3098708.000000000d+00
mas(13)= 1047.348600000000d+00
mas(14)= 3497.898000000000d+00
mas(15)= 22902.9800000000d+00
mas(16)= 1.00000000000000d+00
mas(17)= 135200000.000000d+00
do n=0,nc
gmas(n)=gk2/mas(n)
enddo
sum=0.d0
do n=0,nc
sum=sum+gmas(n)
enddo
c flatness of Neptune
Rua=Rek/UA
cj2=gmas(0)*pj2*(Rua**2)
cj4=gmas(0)*pj4*(Rua**4)
c characteristics of the files
if(ndat.eq.1) then
if(nimp.eq.0) then
timp=tfin-tinit ! interval between 2 display
nbk=0
endif
if(nimp.ne.0) then
timp = nimp*painit ! interval between 2 display
nbk = (nbp/nimp) ! number of intervals
endif
nbi=nbk+1 ! number of display
tps=abs(tfin-tinit)
tps=abs(tps-nbk*abs(timp))
if(tps.gt.1.d-9) nbi=nbi+1 ! number of display
Endif ! (ndat=1)
go to 1633
1633 Continue c Coefficients for r
cr(1,1)=1.d0
cr(1,2)=1.d0
do i=2,nd-1
ih=i+1
j=i-1
cr(i,1)=1.d0
cr(i,ih)=1.d0
do k=2,i
cr(i,k)=cr(j,k)+cr(j,k-1)
enddo
enddo
do i=2,nd-2,2
j=i/2+1
cr(i,j)=cr(i,j)/2.d0
enddo
c Coefficients for g
do l=1,lb
cg(l,1,1)=l
cg(l,2,1)=l
cg(l,2,2)=l+1
do i=3,nd-1
cg(l,i,1)=l
cg(l,i,i)=i+(l-1)
j=i-1
do k=2,j
cg(l,i,k)=cg(l,j,k)+cg(l,j,k-1)
enddo
enddo
enddo
c Paper by Le Guyader 2.2 Matrix begin
ic = 1
do i=1,nc
do j=1,nc
gmu(i,j)=gmas(j)
enddo
enddo
do i=1,nc
gmu(i,i)=gmas(0)+gmas(i)
enddo
jmin=nc+1
jd=nc-1
do i=1,nc-1,ic
l=jmin-i-1
jd=jd-1
jmax=jmin+jd
do j=jmin,jmax
gmu(i,j)=-gmas(j-l)
enddo
jmin=jmax+1
enddo
n=nc
do j=1,nc-1,ic
l=nc-j
do i=1,l
gmu(j+i,n+i)=gmas(j)
enddo
n=n+l
enddo
c Paper by Le Guyader 2.2 Matrix end c Binomial coefficients
Call cik(ck,nd)
c time and time steps
tpi= tinit
tpf= tfin
Time=tinit
pas=painit
pat=0.d0
if(ndat.eq.1) ki=no+1
jch=0
ipr=0
iep=0
tip=abs(tinit)
inr=0
ier=0
tep=abs(tinit)
5000 continue c initialize matrices
do j=1,np
do k=0,nd
r(j,k)=0.d0
do l=1,3
v(l,j,k)=0.d0
enddo
do l=1,lb
g(l,j,k)=0.d0
fx(l,j,k)=0.d0
fy(l,j,k)=0.d0
fz(l,j,k)=0.d0
enddo
enddo
enddo
c Auxilliary variables
do k=0,1
do n=1,nc
do l=1,3
v(l,n,k)=v0(l,n,k)
enddo
enddo
enddo
if(nc.gt.1) then
do k=0,1
call vd(k,nc,nd,np,v)
enddo
endif
rmin=1.d+38
c Paper Le Guyader Equation 7
do j=1,np
r(j,0)=sqrt(v(1,j,0)**2+v(2,j,0)**2+v(3,j,0)**2)
g(3,j,0)=1.d0/(r(j,0)**3)
fx(3,j,0)=v(1,j,0)*g(3,j,0)
fy(3,j,0)=v(2,j,0)*g(3,j,0)
fz(3,j,0)=v(3,j,0)*g(3,j,0)
rmin=min(rmin,r(j,0))
enddo ! end do loop j
pmin=rmin*precr
c second derivatives - Le Guyader Equation 8
do n=1,nc
do j=1,np
v(1,n,2)=v(1,n,2)-gmu(n,j)*fx(3,j,0)
v(2,n,2)=v(2,n,2)-gmu(n,j)*fy(3,j,0)
v(3,n,2)=v(3,n,2)-gmu(n,j)*fz(3,j,0)
enddo
enddo
c flatness region
if(nap.eq.0) go to 340
c alpha and delta of the peel (?) surface (?) of Neptune at time t0
i=0
tj0=(2451545.0d0-Time)
pn(0)=pn0+pn1*tj0
ap(0)=ap0+ap1*sin(pn(0))
dp(0)=dp0+dp1*cos(pn(0))
cap(0)=cos(ap(0))
sap(0)=sin(ap(0))
cdp(0)=cos(dp(0))
sdp(0)=sin(dp(0))
cc(0)=cdp(0)*cap(0)
cs(0)=cdp(0)*sap(0)
sd(0)=sdp(0)
do jp=npi,npf
g(1,jp,0) =1.d0/r(jp,0)
fx(1,jp,0)=v(1,jp,0)*g(1,jp,0)
fy(1,jp,0)=v(2,jp,0)*g(1,jp,0)
fz(1,jp,0)=v(3,jp,0)*g(1,jp,0)
sa(1,jp,0)=cc(0)*fx(1,jp,0)+cs(0)*fy(1,jp,0)+sd(0)*fz(1,jp,0)
sa(2,jp,0)=sa(1,jp,0)*sa(1,jp,0)
sa(3,jp,0)=sa(2,jp,0)*sa(1,jp,0)
sa(4,jp,0)=sa(3,jp,0)*sa(1,jp,0)
acc(jp,0)=sa(1,jp,0)*cc(0)
acs(jp,0)=sa(1,jp,0)*cs(0)
asd(jp,0)=sa(1,jp,0)*sd(0)
psa(2,jp,0)=15.d0/2.d0*sa(2,jp,0)-3.d0/2.d0
psa(4,jp,0)=315d0/8.d0*sa(4,jp,0)-210.d0/8.d0*sa(2,jp,0)
1 +15.d0/8.d0
ppa(4,jp,0)=35.d0/2.d0*sa(2,jp,0)-30.d0/4.d0
do lp=4,6,2
g(lp,jp,0)=1.d0/(r(jp,0)**lp)
enddo
gcj(2,jp,0)=cj2*g(4,jp,0)
gcj(4,jp,0)=cj4*g(6,jp,0)
pcc(2,jp,0)=psa(2,jp,0)*fx(1,jp,0)-3.d0*acc(jp,0)
pcs(2,jp,0)=psa(2,jp,0)*fy(1,jp,0)-3.d0*acs(jp,0)
psd(2,jp,0)=psa(2,jp,0)*fz(1,jp,0)-3.d0*asd(jp,0)
pcc(4,jp,0)=psa(4,jp,0)*fx(1,jp,0)-ppa(4,jp,0)*acc(jp,0)
pcs(4,jp,0)=psa(4,jp,0)*fy(1,jp,0)-ppa(4,jp,0)*acs(jp,0)
psd(4,jp,0)=psa(4,jp,0)*fz(1,jp,0)-ppa(4,jp,0)*asd(jp,0)
plx(jp,0)=gcj(2,jp,0)*pcc(2,jp,0)+gcj(4,jp,0)*pcc(4,jp,0)
ply(jp,0)=gcj(2,jp,0)*pcs(2,jp,0)+gcj(4,jp,0)*pcs(4,jp,0)
plz(jp,0)=gcj(2,jp,0)*psd(2,jp,0)+gcj(4,jp,0)*psd(4,jp,0)
v(1,jp,2)=v(1,jp,2)+plx(jp,0)
v(2,jp,2)=v(2,jp,2)+ply(jp,0)
v(3,jp,2)=v(3,jp,2)+plz(jp,0)
enddo ! end do loop jp (npi,npf)
c end of flatness region (not in Le Guyader) 340 continue
if(nc.gt.1) then
call vd(2,nc,nd,np,v)
endif
tps=abs(tfin-Time)
if(tps.le.1.d-9) then
go to 2255
inr=inr+1
if(id.eq.0) go to 201
201 Continue
if(is.eq.0) go to 301
301 Continue
ipr=ipr+1
write(6,210) pas
210 format(8x,' Last step:',d23.16)
write(6,211) ipr,Time
211 format(8x,'(',i3,')',' Time final:',d23.16)
do n=1,nc
if(n.lt.10) then
do l=1,3
write(6,212) l,n,v0(l,n,0)
212 format(8x,'v0(',i1,',',i1,',0) =',d23.16)
enddo
do l=1,3
write(6,213) l,n,v0(l,n,1)
213 format(8x,'v0(',i1,',',i1,',1) =',d23.16)
enddo
endif
if(n.ge.10) then
do l=1,3
write(6,214) l,n,v0(l,n,0)
214 format(8x,'v0(',i1,',',i2,',0)=',d23.16)
enddo
do l=1,3
write(6,215) l,n,v0(l,n,1)
215 format(8x,'v0(',i1,',',i2,',1)=',d23.16)
enddo
endif
Enddo ! end do loop n
Call intp(nc,nd,np,gmas,v,aia)
if(nap.eq.1) then
write(*,*) ' Values close to the first integrals '
endif
write(6,216) aia(4)
216 format(8x,'Energy:',d23.16)
write(6,217) aia(1),aia(2),aia(3)
217 format(8x,'Surfaces :',d23.16,d23.16,d23.16)
if(ndat.eq.1) then
write(6,218)
218 format(27x,'End " File Ephemeride "')
endif
write(6,223) jch
223 format(20x,' number of changes of steps :',i3)
write(6,224)
224 format(28x,'Stop: end program ') 2255 Continue
Return
Endif
go to 2399
2399 Continue
if(pas.eq.painit) go to 450
pat=pat+pas
if(abs(painit-pat).gt.1.d-12) go to 450
pat=0.d0
pas=painit
450 Time=Time+pas
if(Time*pas.gt.tfin*pas) then
Time=Time-pas
pas=tfin-Time
Time=Time+pas
endif
c GUESS higher order derivatives
do i=1,nd-2
k=i+2
ip=i/2
iq=i-1
do j=1,np
wx=v(1,j,0)
wy=v(2,j,0)
wz=v(3,j,0)
do l=1,ip
ih=i-l
im=l+1
r(j,i)=r(j,i)+cr(i,im)*(v(1,j,l)*v(1,j,ih)+
1 v(2,j,l)*v(2,j,ih)+v(3,j,l)*v(3,j,ih)-r(j,l)*r(j,ih))
enddo
r(j,i)=(r(j,i)+wx*v(1,j,i)+wy*v(2,j,i)+wz*v(3,j,i))/r(j,0)
do l=0,iq
ih=i-l
im=l+1
g(3,j,i)=g(3,j,i)+cg(3,i,im)*g(3,j,l)*r(j,ih)
enddo
g(3,j,i)=-g(3,j,i)/r(j,0)
do l=0,ip
ih=i-l
im=l+1
w1=cr(i,im)*g(3,j,ih)
w2=cr(i,im)*g(3,j,l)
fx(3,j,i)=fx(3,j,i)+w1*v(1,j,l)+w2*v(1,j,ih)
fy(3,j,i)=fy(3,j,i)+w1*v(2,j,l)+w2*v(2,j,ih)
fz(3,j,i)=fz(3,j,i)+w1*v(3,j,l)+w2*v(3,j,ih)
enddo
enddo ! end do loop j
c kth order derivatives - Le Guyader Equation 8
do n=1,nc
do j=1,np
v(1,n,k)=v(1,n,k)-gmu(n,j)*fx(3,j,i)
v(2,n,k)=v(2,n,k)-gmu(n,j)*fy(3,j,i)
v(3,n,k)=v(3,n,k)-gmu(n,j)*fz(3,j,i)
enddo
enddo
c flatness region
if(nap.eq.0) go to 350
c derivatives (?) of sine cosine of the peel (?) surface (?)
if(i.eq.1) then
ap(1)=-ap1*pn1*cos(pn(0))
dp(1)= dp1*pn1*sin(pn(0))
cap(1)=-ap(1)*sap(0)
cdp(1)=-dp(1)*sdp(0)
sap(1)= ap(1)*cap(0)
sdp(1)= dp(1)*cdp(0)
go to 799
endif
if(i.eq.2) then
ap(2)=-ap1*(pn1**2)*sin(pn(0))
dp(2)=-dp1*(pn1**2)*cos(pn(0))
cap(2)=-ap(2)*sap(0)-ap(1)*sap(1)
cdp(2)=-dp(2)*sdp(0)-dp(1)*sdp(1)
sap(2)= ap(2)*cap(0)+ap(1)*cap(1)
sdp(2)= dp(2)*cdp(0)+dp(1)*cdp(1)
go to 799
endif
ap(i)=-(pn1**2)*ap(i-2)
dp(i)=-(pn1**2)*dp(i-2)
cap(i)=0.d0
cdp(i)=0.d0
sap(i)=0.d0
sdp(i)=0.d0
do kp=0,i-1
cap(i)=cap(i)-ck(i-1,kp)*ap(kp+1)*sap(i-1-kp)
cdp(i)=cdp(i)-ck(i-1,kp)*dp(kp+1)*sdp(i-1-kp)
sap(i)=sap(i)+ck(i-1,kp)*ap(kp+1)*cap(i-1-kp)
sdp(i)=sdp(i)+ck(i-1,kp)*dp(kp+1)*cdp(i-1-kp)
enddo
799 continue
cc(i)=0.d0
cs(i)=0.d0
sd(i)=0.d0
do kp=0,i
cc(i)=cc(i)+ck(i,kp)*cdp(kp)*cap(i-kp)
cs(i)=cs(i)+ck(i,kp)*cdp(kp)*sap(i-kp)
enddo
sd(i)=sdp(i)
do jp=npi,npf
g(1,jp,i)=0.d0
do l=0,iq
ih=i-l
im=l+1
g(1,jp,i)=g(1,jp,i)+cg(1,i,im)*g(1,jp,l)*r(jp,ih)
enddo
g(1,jp,i)=-g(1,jp,i)/r(jp,0)
fx(1,jp,i)=0.d0
fy(1,jp,i)=0.d0
fz(1,jp,i)=0.d0
do l=0,ip
ih=i-l
im=l+1
w1=cr(i,im)*g(1,jp,ih)
w2=cr(i,im)*g(1,jp,l)
fx(1,jp,i)=fx(1,jp,i)+w1*v(1,jp,l)+w2*v(1,jp,ih)
fy(1,jp,i)=fy(1,jp,i)+w1*v(2,jp,l)+w2*v(2,jp,ih)
fz(1,jp,i)=fz(1,jp,i)+w1*v(3,jp,l)+w2*v(3,jp,ih)
enddo
sa(1,jp,i)=0.d0
sa(2,jp,i)=0.d0
sa(3,jp,i)=0.d0
sa(4,jp,i)=0.d0
do l=0,ip
ih=i-l
im=l+1
w11=cr(i,im)*cc(ih)
w21=cr(i,im)*cc(l)
w12=cr(i,im)*cs(ih)
w22=cr(i,im)*cs(l)
w13=cr(i,im)*sd(ih)
w23=cr(i,im)*sd(l)
sa(1,jp,i)=sa(1,jp,i)+w11*fx(1,jp,l)+w21*fx(1,jp,ih)
1 +w12*fy(1,jp,l)+w22*fy(1,jp,ih)
2 +w13*fz(1,jp,l)+w23*fz(1,jp,ih)
enddo
do lp=2,4
do l=0,ip
ih=i-l
im=l+1
w11=cr(i,im)*sa(1,jp,ih)
w21=cr(i,im)*sa(1,jp,l)
sa(lp,jp,i)=sa(lp,jp,i)+w11*sa(lp-1,jp,l)+w21*sa(lp-1,jp,ih)
enddo
enddo
acc(jp,i)=0.d0
acs(jp,i)=0.d0
asd(jp,i)=0.d0
do l=0,ip
ih=i-l
im=l+1
w11=cr(i,im)*cc(ih)
w21=cr(i,im)*cc(l)
w12=cr(i,im)*cs(ih)
w22=cr(i,im)*cs(l)
w13=cr(i,im)*sd(ih)
w23=cr(i,im)*sd(l)
acc(jp,i)=acc(jp,i)+w11*sa(1,jp,l)+w21*sa(1,jp,ih)
acs(jp,i)=acs(jp,i)+w12*sa(1,jp,l)+w22*sa(1,jp,ih)
asd(jp,i)=asd(jp,i)+w13*sa(1,jp,l)+w23*sa(1,jp,ih)
enddo
c derivatives (?) of the Legendre polynomials
psa(2,jp,i)=15.d0/2.d0*sa(2,jp,i)
psa(4,jp,i)=315.d0/8.d0*sa(4,jp,i)-210.d0/8.d0*sa(2,jp,i)
ppa(4,jp,i)=35.d0/2.d0*sa(2,jp,i)
do lp=4,6,2
g(lp,jp,i)=0.d0
do l=0,iq
ih=i-l
im=l+1
g(lp,jp,i)=g(lp,jp,i)+cg(lp,i,im)*g(lp,jp,l)*r(jp,ih)
enddo
g(lp,jp,i)=-g(lp,jp,i)/r(jp,0)
enddo ! end do loop lp(4,6)
gcj(2,jp,i)=cj2*g(4,jp,i)
gcj(4,jp,i)=cj4*g(6,jp,i)
do lp=2,4,2
pcc(lp,jp,i)=0.d0
pcs(lp,jp,i)=0.d0
psd(lp,jp,i)=0.d0
enddo
do kp=0,i
pcc(2,jp,i)=pcc(2,jp,i)+ck(i,kp)*psa(2,jp,kp)*fx(1,jp,i-kp)
pcs(2,jp,i)=pcs(2,jp,i)+ck(i,kp)*psa(2,jp,kp)*fy(1,jp,i-kp)
psd(2,jp,i)=psd(2,jp,i)+ck(i,kp)*psa(2,jp,kp)*fz(1,jp,i-kp)
enddo
pcc(2,jp,i)=pcc(2,jp,i)-3.d0*acc(jp,i)
pcs(2,jp,i)=pcs(2,jp,i)-3.d0*acs(jp,i)
psd(2,jp,i)=psd(2,jp,i)-3.d0*asd(jp,i)
do kp=0,i
pcc(4,jp,i)=pcc(4,jp,i)+ck(i,kp)*(psa(4,jp,kp)*fx(1,jp,i-kp)
1 -ppa(4,jp,kp)*acc(jp,i-kp))
pcs(4,jp,i)=pcs(4,jp,i)+ck(i,kp)*(psa(4,jp,kp)*fy(1,jp,i-kp)
1 -ppa(4,jp,kp)*acs(jp,i-kp))
psd(4,jp,i)=psd(4,jp,i)+ck(i,kp)*(psa(4,jp,kp)*fz(1,jp,i-kp)
1 -ppa(4,jp,kp)*asd(jp,i-kp))
enddo
plx(jp,i)=0.d0
ply(jp,i)=0.d0
plz(jp,i)=0.d0
do kp=0,i
plx(jp,i)=plx(jp,i)+ck(i,kp)*(gcj(2,jp,kp)*pcc(2,jp,i-kp)
1 +gcj(4,jp,kp)*pcc(4,jp,i-kp))
ply(jp,i)=ply(jp,i)+ck(i,kp)*(gcj(2,jp,kp)*pcs(2,jp,i-kp)
1 +gcj(4,jp,kp)*pcs(4,jp,i-kp))
plz(jp,i)=plz(jp,i)+ck(i,kp)*(gcj(2,jp,kp)*psd(2,jp,i-kp)
1 +gcj(4,jp,kp)*psd(4,jp,i-kp))
enddo
v(1,jp,k)=v(1,jp,k)+plx(jp,i)
v(2,jp,k)=v(2,jp,k)+ply(jp,i)
v(3,jp,k)=v(3,jp,k)+plz(jp,i)
enddo ! end do loop jp (npi,npf)
350 continue c end of flatness region (not in Le Guyader)
if(nc.gt.1) then
call vd(k,nc,nd,np,v)
endif
enddo ! end do loop i
! Test
trt=0.d0
do n=1,nc
do l=1,3
trt=trt+abs(v(l,n,nd))
enddo
enddo
410 call tn(nd,pas,t)
tst=trt*abs(t(k))
if(tst.lt.pmin) go to 420
pas=pas/2.d0
Time=Time-pas
apa=abs(pas)
if(apa.lt.api) then
write(6,246) jch,painit,pas
246 format(8x,'Stop: steps have become too small:',i4,d23.16,d23.16)
stop
endif
jch=jch+1
if(jch.gt.nch) then
write(6,249) jch
249 format(18x,'Stop: too many changes of steps:',i4)
stop
endif
go to 410
420 if(ki.gt.no) go to 430 430 Continue ! Reprise of normal steps c Integrate to get coordinates and speeds
do k=0,1
do n=1,nc
do l=1,3
v0(l,n,k)=0.d0
do kp=nd-k,0,-1
v0(l,n,k)=v0(l,n,k)+v(l,n,kp+k)*t(kp)
enddo
enddo
enddo
enddo
go to 5000
Return
End
subroutine cik(ck,nd)
implicit double precision (a-h,m,o-z)
dimension ck(0:nd,0:nd)
ck(0,0)=1.d0
ck(1,0)=1.d0
ck(1,1)=1.d0
do i=2,nd
j=i-1
ck(i,0)= 1.d0
ck(i,i)= 1.d0
do k=1,j
ck(i,k)=ck(j,k-1)+ck(j,k)
enddo
enddo
return
end
c time steps = ( step size ) / (i!)
subroutine tn(nd,pas,t)
implicit double precision (a-h,m,o-z)
dimension t(0:nd)
t(0)=1.d0
do i=1,nd
t(i)=t(i-1)*pas/i
enddo
return
end
subroutine vd(k,nc,nd,np,v)
implicit double precision (a-h,m,o-z)
dimension v(3,np,0:nd),w(3,np,0:nd)
jmin=nc+1
jd=nc-1
do i=1,nc-1
jm=jmin-i-1
jd=jd-1
jmax=jmin+jd
do j=jmin,jmax
do l=1,3
v(l,j,k)=v(l,j-jm,k)-v(l,i,k)
enddo
enddo
jmin=jmax+1
enddo
return
end
subroutine wd(k,nc,nd,np,w)
implicit double precision (a-h,m,o-z)
dimension w(3,np,0:nd)
jmin=nc+1
jd=nc-1
do i=1,nc-1
jm=jmin-i-1
jd=jd-1
jmax=jmin+jd
do j=jmin,jmax
do l=1,3
w(l,j,k)=w(l,j-jm,k)-w(l,i,k)
enddo
enddo
jmin=jmax+1
enddo
return
end
Usage documentation supplied at end of code
First, after an integration step-size has been fixed, the optimal number of derivatives is determined, which will be used for every solution of the Taylor series.
It is important to obtain a numerical drift as smallest as possible for a minimal computation time.
Indeed if the number of derivatives used in the series is not sufficient, the program changes the step and this increases the round-off errors. If the number of derivatives is large, then the integration time will also be large with no changes in the precision.
In order to avoid these disadvantages, the following procedures are established:
One chooses one of the fastest bodies, (for example, Triton) and computes the position and velocity vectors at the end of one or several revolutions, beginning with an insufficient number of derivatives.
This integration is reinitiated many times but after at every integration a new derivative is added to the series until the difference between two successive integrations has no effect on the required accuracy. Thus, the optimal number of derivatives is chosen. However, if this number is too large, it will be necessary to choose a smaller integration size.
