rk78cn.f

Go to the documentation of this file.

      SUBROUTINE rk78cn
C  VERSION OF 4/1/85
C  PURPOSE
C    COMPUTES THE FEHLBERG COEFFICIENTS FOR A RUNGE-KUTTA 78 INTEGRATOR
C  INPUT
C    NONE
C  COMMON BLOCK OUTPUT
C    CH, ALPH, BETA, NORDER, ORDRCP, NTIMES
C  REFERENCES
C    JPL EM 312/85-153, 20 APRIL 1987
C    NASA TR R-287, OCTOBER 1968
C  ANALYSIS
C    J. H. KWOK - JPL
C  PROGRAMMER
C    J. H. KWOK - JPL
C  PROGRAM MODIFICATIONS
C    NONE
C  COMMENTS
C    THIS ROUTINE MUST BE CALLED BEFORE CALLING ROUTINE RK78
C
      IMPLICIT DOUBLE PRECISION (a-h,o-z)   
      common/felcon/ch(13),alph(13),beta(13,12),ordrcp,norder,ntimes
      DATA zero,one,two,three/0.d0,1.d0,2.d0,3.d0/
      norder=8
      ordrcp=one/norder
      ntimes=13 
      DO 30 i=1,13  
      DO 40 j=1,12  
   40 beta(i,j)=zero
      alph(i)=zero
   30 ch(i)=zero
      ch(6)=34.d0/105.d0
      ch(7)=9.d0/35.d0  
      ch(8)=ch(7)   
      ch(9)=9.d0/280.d0 
      ch(10)=ch(9)  
      ch(12)=41.d0/840.d0   
      ch(13)=ch(12) 
      alph(2)=two/27.d0
      alph(3)=one/9.d0 
      alph(4)=one/6.d0 
      alph(5)=5.d0/12.d0
      alph(6)=.5d0  
      alph(7)=5.d0/6.d0 
      alph(8)=one/6.d0 
      alph(9)=two/three 
      alph(10)=one/three
      alph(11)=one
      alph(13)=one 
      beta(2,1)=two/27.d0  
      beta(3,1)=one/36.d0  
      beta(4,1)=one/24.d0  
      beta(5,1)=5.d0/12.d0  
      beta(6,1)=.5d-1   
      beta(7,1)=-25.d0/108.d0   
      beta(8,1)=31.d0/300.d0
      beta(9,1)=two
      beta(10,1)=-91.d0/108.d0  
      beta(11,1)=2383.d0/4100.d0
      beta(12,1)=three/205.d0
      beta(13,1)=-1777.d0/4100.d0   
      beta(3,2)=one/12.d0  
      beta(4,3)=one/8.d0   
      beta(5,3)=-25.d0/16.d0
      beta(5,4)=-beta(5,3)  
      beta(6,4)=.25d0   
      beta(7,4)=125.d0/108.d0   
      beta(9,4)=-53.d0/6.d0 
      beta(10,4)=23.d0/108.d0   
      beta(11,4)=-341.d0/164.d0 
      beta(13,4)=beta(11,4) 
      beta(6,5)=.2d0
      beta(7,5)=-65.d0/27.d0
      beta(8,5)=61.d0/225.d0
      beta(9,5)=704.d0/45.d0
      beta(10,5)=-976.d0/135.d0 
      beta(11,5)=4496.d0/1025.d0
      beta(13,5)=beta(11,5) 
      beta(7,6)=125.d0/54.d0
      beta(8,6)=-two/9.d0  
      beta(9,6)=-107.d0/9.d0
      beta(10,6)=311.d0/54.d0   
      beta(11,6)=-301.d0/82.d0  
      beta(12,6)=-6.d0/41.d0
      beta(13,6)=-289.d0/82.d0  
      beta(8,7)=13.d0/900.d0
      beta(9,7)=67.d0/90.d0 
      beta(10,7)=-19.d0/60.d0   
      beta(11,7)=2133.d0/4100.d0
      beta(12,7)=-three/205.d0   
      beta(13,7)=2193.d0/4100.d0
      beta(9,8)=three
      beta(10,8)=17.d0/6.d0 
      beta(11,8)=45.d0/82.d0
      beta(12,8)=-three/41.d0
      beta(13,8)=51.d0/82.d0
      beta(10,9)=-one/12.d0
      beta(11,9)=45.d0/164.d0   
      beta(12,9)=three/41.d0 
      beta(13,9)=33.d0/164.d0   
      beta(11,10)=18.d0/41.d0   
      beta(12,10)=6.d0/41.d0
      beta(13,10)=12.d0/41.d0   
      beta(13,12)=one  
      RETURN
      END