Talk:Help with Fortran90

From Stellarium Wiki
Jump to: navigation, search

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.

Personal tools
Namespaces
Variants
Actions
in this wiki
other languages
Toolbox