Hi everyone, I can not compile a program here (RADAU-INSTRUCOES.for), which uses the program radau ra15, (ra27.for)
When I compile the program RADAU-INSTRUCOES.for, the following errors appear:
reading the errors, you can see that one of the problems is the dimension of the variable x, which was defined in x (1).
However after defining the variable x as x (3), the following errors appear in ra27.for:
RADAU-INSTRUCOES.FOR (WITHOUT CHANGES IN X):
ra27.for:
.
When I compile the program RADAU-INSTRUCOES.for, the following errors appear:
Code:
./RADAU-INSTRUCOES.for:229:38:
f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! f (i) are the first derivatives of x, y, z
1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:229:25:
f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! f (i) are the first derivatives of x, y, z
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:230:30:
f(2)=-bb*x(2)+ff*x(3)-x(2)
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:230:25:
f(2)=-bb*x(2)+ff*x(3)-x(2)
1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:230:17:
f(2)=-bb*x(2)+ff*x(3)-x(2)
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:231:35:
f(3)=5.d-4*bb*x(1)-ff*x(3)-x(2)
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:231:30:
f(3)=5.d-4*bb*x(1)-ff*x(3)-x(2)
1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:285:29:
WRITE(42,123)tm,X(1),X(2),X(3) ! each output_step (tm current) outputs the results
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:285:34:
WRITE(42,123)tm,X(1),X(2),X(3) ! each output_step (tm current) outputs the results
1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:286:28:
WRITE(*,123)tm,X(1),X(2),X(3) ! save to file and display on screen
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:286:33:
WRITE(*,123)tm,X(1),X(2),X(3) ! save to file and display on screen
1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
ra27.for:131:72:
+OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
1
Warning: More actual than formal arguments in procedure call at (1)
ra27.for:138:72:
call OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
1
Warning: More actual than formal arguments in procedure call at (1)
ra27.for:264:72:
CALL OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,NPER,iemax)
1
Warning: More actual than formal arguments in procedure call at (1)
/tmp/cc9kC96U.o: Na função "MAIN__":
RADAU-INSTRUCOES.for:(.text+0x36e9): referência não definida para "ra27_"
collect2: error: ld returned 1 exit statusreading the errors, you can see that one of the problems is the dimension of the variable x, which was defined in x (1).
However after defining the variable x as x (3), the following errors appear in ra27.for:
Code:
ra27.for:107:18:
CALL FORCE (X, V, ZERO, F1)
1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:131:19:
+OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:138:24:
call OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:264:24:
CALL OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,NPER,iemax)
1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:268:18:
78 CALL FORCE (X,V,TM,F1)
1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
/tmp/ccfQadBO.o: Na função "MAIN__":
RADAU-INSTRUCOES.for:(text+0x36e9): referência não definida para "ra27_"
collect2: error: ld returned 1 exit statusRADAU-INSTRUCOES.FOR (WITHOUT CHANGES IN X):
Code:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c: here is this trivial prog that uses the package: RA27.for (next to
c: Ra15 which is in fact an integration routine that uses here)
c: Still this black box is modified slightly
c: View original: Everhart article data output.
c: User must give Force.for and as an example goes the OUTPUT.for (routine
c: data output)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
INCLUDE 'ra27.for' ! include the package RA27.for
C: Just copy the RA27.for file to the same folder you are working on (same as this main program)
IMPLICIT REAL*8(A-H,O-Z)
logical fixed
CHARACTER FName1*30
common/ctes/aa,bb,cc,dd,ee,ff
data aa,bb,cc,dd,ee,ff/1.d0,0.2d0,0.1d0,0.6d0,-2.d0,-1.d0/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C C
C C
C
C C
C C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C dimension x: coordinate and v: velocity
c: in its case x = vector (x, y, z)
c: v = vector velocity of x, y, z: you do not need, your units are from 1st order
DIMENSION X(3),v(3) ! Assuming 3 first-order equations
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C NV
C
C
C
c: Nv: number of simultaneous events (maximum is 18, but
c: user can move (F1, FJ)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
NV=3
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C NCLASS describes the type of system to be integrated
C y'=F(y,t) NCLASS=1
C y"=F(y,t) NCLASS= -2
C y"=F(y',y,t) NCLASS=2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
NCLASS=1 ! your system is the first order
c: eLagrange equations always sound first order, so nclass = 1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C NOR order of the integrator ie (15 or 27)
C In case here, always use NOR = 15 (Radau15)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
NOR=15
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C LL controls to a certain extent accuracy (sequence size: see article).
C SS=10**(-LL) .
C A good value of LL e ': half of NOR (maximum I think it is 12).
C LL large-> high precision
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
LL=10
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C XL if LL <= 0 then XL forms a constant-length sequence
C The smallest XL toleravel: XL = 0.01
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
XL=0.01D0
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TF integration interval TF = t (final) - t (initial)
C Logic that may be negative: integrate to Back
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
tm=0.d0 ! time starts
TF=10.D0 ! final integration time
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C X e V : position and initial velocity given by the user
C
C
C Entre aqui X(t=To) e V(t=To)=X'
c In case of first order, you do not need to give V (t = to)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C PI with 32 digits. (German things)
C
PI=3.1415926535897932384626433832795D0
C:INICIAL CONDITIONS of the variables x (1) x (2), x (3) (speed you do not have, do not need)
x(1)=0.2D0 ! would be x (1) at t = t0
x(2)=0.3d0 ! x(2) em t0
x(3)=-0.05d0 ! x(3) em t0
c: ! in your case you do not need speeds
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C OUTPUT_STEP: tells the integrator: exit result. to each (APROX) OUTPUT_step (time t)
C OUTPUT and called OUTPUT_STEP
C
C fixed: true if I want to exit in exactly OUTPUT_STEP
C
C: false if I do not ask too much, this is completed
c: sequence releases an output
C: Whenever possible use FIXED = .false. (much faster)
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
OUTPUT_STEP=1.5D0 !
c: in fact it will not be exactly 1.5 (eh approx), if you want exactly 1.5, then change fixed to .true,
c: but alas, it gets a little slower, of course.
c: fixed=.true. ! here every exact output_step output the results, fixed pitch
fixed=.false. ! in this case each output_step output the results (approximately)
c: use fixed = false, the integrator walks much faster
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Filename for Output file: from the file where your outputs will be saved
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c: FName1='/tmp/bla2.xyz'
FName1='apg.dat' !name of the file, to save the outputs
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
open (42,file=FName1, status='unknown')
write(42,*) ' Ordem da Integracao = ',NOR,
+' SequenzSize = ', LL
write(*,*) ' Ordem da Integracao = ',NOR,
+' LL (precisao) = ', LL
123 FORMAT(7F20.10)
5 FORMAT (A)
WRITE(*,5) ' Favor esperar: estou calculando ........'
IF (NOR.eq.15)
+ CALL RA15 (X,V,TF,XL,LL,NV,NCLASS,NOR,OUTPUT_STEP,fixed)
IF (NOR.gt.15)
+ CALL RA27 (X,V,TF,XL,LL,NV,NCLASS,NOR,OUTPUT_STEP,fixed)
Close(42)
write(*,5) ' Pronto, FIM '
STOP
END
c:ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C C Now provide the Force.for subroutine, where it must conform to
C Nclass that you have defined before.
c
c
C SUBROUTINE FORCE (X, V, TM, F)
C X, V: see previous definition
C TM: Current time
C F: forca: vector of NV components
c: in general, Lagrange has no time = tm in the second limb
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE FORCE(X,V,TM,F)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION X(1),V(1),F(3)
common/ctes/aa,bb,cc,dd,ee,ff
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Here you should give the components of the velocities of vector X (be careful
c: there must be agreement with NCLASS from the Main prog).
c:
c: F (1) = f1 (x, t)
c: F (2) = f2 (x, t)
c: ..............
c: F (nv) = fnv (x, t)
c Let X = [x (1), x (2), x (3)] be the easiest thing to do here:
c (1) = dx (1) / dt, f (2) = dx (2) / dt and f (3) = dx (3 / dt:
c: aa, bb, ff are constants given in prog. main
f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! f (i) are the first derivatives of x, y, z
f(2)=-bb*x(2)+ff*x(3)-x(2)
f(3)=5.d-4*bb*x(1)-ff*x(3)-x(2)
RETURN
END
C ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C This is the routine that prints the outputs on the screen and in the file apg.apg
c
C SUBROUTINE OUTPUT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
SUBROUTINE OUTPUT(NF,NS,X,V,F,T,TM,TF,FIRST,LAST)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C NF :Number of iterations made by Radau
C
C NS : number of steps in time
C
C X,V X and V : said before
C
C F power vector
C
C T size of the next step (time t)
C
C TM current time
C
C TF TF=t(final) - t (initial) may be negative
C
C FIRST logical variable: truth only at the beginning
C The integrator can test several times at the beginning until it finds a good
c: to start
c:
C
cc All this information can be ignored in a first read, Just remember tm = current time
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
IMPLICIT REAL*8(A-H,O-Z)
INTEGER NF,NS,N
LOGICAL LAST,FIRST
DIMENSION X(1),V(1),F(1) ! as this is a fortran 77, in general the vectors of
c: subroutines, are scaled to 1 only. If you are going to use Gfortran, then you should set the maximum, like
c: was defined in the main programl (3)
5 FORMAT (A)
123 FORMAT(F16.8, 3d21.10)
C Here I am writing the outputs: screen and file in unit 42 (attention
c: this 42 must agree with the OPEN: prog. Main)
C: if there is an integral integral.teste a preciso com several LL !!!
WRITE(42,123)tm,X(1),X(2),X(3) ! each output_step (tm current) outputs the results
WRITE(*,123)tm,X(1),X(2),X(3) ! save to file and display on screen
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccc Fim ccccccccccccccccccccccccccccccccccccra27.for:
Code:
C $$$$ RA15T.FOR
C
SUBROUTINE RA15(X,V,TF,XL,LL,NV,NCLASS,NOR,OS,FIXED)
C Integrator RADAU by E. Everhart, Physics Department, University of Denver
C This 15th-order version, called RA15, is written out for faster execution.
C y'=F(y,t) is NCLASS=1, y"=F(y,t) is NCLASS= -2, y"=F(y',y,t) is NCLASS=2
C TF is t(final) - t(initial). It may be negative for backward integration.
C NV = the number of simultaneous differential equations.
C The dimensioning below assumes NV will not be larger than 18.
C LL controls sequence size. Thus SS=10**(-LL) controls the size of a term.
C A typical LL-value is in the range 6 to 12 for this order 11 program.
C However, if LL.LT.0 then XL is the constant sequence size used.
C X and V enter as the starting position-velocity vector, and are output as
C the final position-velocity vector.
C Integration is in double precision. A 64-bit double-word is assumed.
IMPLICIT REAL*8 (A-H,O-Z)
REAL *4 TVAL,PW
REAL*8 NEXTOUTPUT
DIMENSION X(1),V(1),F1(18),FJ(18),C(21),D(21),R(21),Y(18),Z(18),
A B(7,18),G(7,18),E(7,18),BD(7,18),H(8),W(7),U(7),NW(8)
c: 25/11/98 introduzi iemax: introduce iemax: stop integration if eccentricity reaches 0.7
c: iemax goes into the OUTPUT_STEP subroutine
LOGICAL NPQ,NSF,NPER,NCL,NES,fixed
DATA NW/0,0,1,3,6,10,15,21/
DATA ZERO, HALF, ONE,SR/0.0D0, 0.5D0, 1.0D0,1.4D0/
C These H values are the Gauss-Radau spacings, scaled to the range 0 to 1,
C for integrating to order 15.
DATA H/ 0.D0, .05626256053692215D0, .18024069173689236D0,
A.35262471711316964D0, .54715362633055538D0, .73421017721541053D0,
B.88532094683909577D0, .97752061356128750D0/
C The sum of the H-values should be 3.73333333333333333
c: 26/11/98 : large eccentricity control with IEMAX
iemax=0
NPER=.FALSE.
NSF=.FALSE.
NCL=NCLASS.EQ.1
NPQ=NCLASS.LT.2
Out=OS
C y'=F(y,t) NCL=.TRUE. y"=F(y,t) NCL=.FALSE. y"=F(y',y,t) NCL=.FALSE.
C NCLASS=1 NPQ=.TRUE. NCLASS= -2 NPQ=.TRUE. NCLASS= 2 NPQ=.FALSE.
C NSF is .FALSE. on starting sequence, otherwise .TRUE.
C NPER is .TRUE. only on last sequence of the integration.
C NES is .TRUE. only if LL is negative. Then the sequence size is XL.
DIR=ONE
IF(TF.LT.ZERO) DIR=-ONE
NES=LL.LT.0
XL=DIR*DABS(XL)
PW=1./9.
C Evaluate the constants in the W-, U-, C-, D-, and R-vectors
DO 14 N=2,8
WW=N+N*N
IF(NCL) WW=N
W(N-1)=ONE/WW
WW=N
14 U(N-1)=ONE/WW
DO 22 K=1,NV
IF(NCL) V(K)=ZERO
DO 22 L=1,7
BD(L,K)=ZERO
22 B(L,K)=ZERO
W1=HALF
IF(NCL) W1=ONE
C(1)=-H(2)
D(1)=H(2)
R(1)=ONE/(H(3)-H(2))
LA=1
LC=1
DO 73 K=3,7
LB=LA
LA=LC+1
LC=NW(K+1)
C(LA)=-H(K)*C(LB)
C(LC)=C(LA-1)-H(K)
D(LA)=H(2)*D(LB)
D(LC)=-C(LC)
R(LA)=ONE/(H(K+1)-H(2))
R(LC)=ONE/(H(K+1)-H(K))
IF(K.EQ.3) GO TO 73
DO 72 L=4,K
LD=LA+L-3
LE=LB+L-4
C(LD)=C(LE)-H(K)*C(LE+1)
D(LD)=D(LE)+H(L-1)*D(LE+1)
72 R(LD)=ONE/(H(K+1)-H(L-1))
73 CONTINUE
SS=10.**(-LL)
C The statements above are used only once in an integration to set up the
C constants. They use less than a second of execution time. Next set in
C a reasonable estimate to TP based on experience. Same sign as DIR.
C An initial first sequence size can be set with XL even with LL positive.
TP=0.1D0*DIR
IF(XL.NE.ZERO) TP=XL
IF(NES) TP=XL
IF(TP/TF.GT.HALF) TP=HALF*TF
NCOUNT=0
C An * is the symbol for writing on the monitor. The printer is unit 4.
C Line 4000 is the starting place of the first sequence.
4000 NS=0
NF=0
NI=6
TM=ZERO
CALL FORCE (X, V, ZERO, F1)
NF=NF+1
C Line 722 is begins every sequence after the first. First find new beta-
C values from the predicted B-values, following Eq. (2.7) in text.
722 DO 58 K=1,NV
G(1,K)=B(1,K)+D(1)*B(2,K)+
X D(2)*B(3,K)+D(4)*B(4,K)+D( 7)*B(5,K)+D(11)*B(6,K)+D(16)*B(7,K)
G(2,K)= B(2,K)+
X D(3)*B(3,K)+D(5)*B(4,K)+D( 8)*B(5,K)+D(12)*B(6,K)+D(17)*B(7,K)
G(3,K)=B(3,K)+D(6)*B(4,K)+D( 9)*B(5,K)+D(13)*B(6,K)+D(18)*B(7,K)
G(4,K)= B(4,K)+D(10)*B(5,K)+D(14)*B(6,K)+D(19)*B(7,K)
G(5,K)= B(5,K)+D(15)*B(6,K)+D(20)*B(7,K)
G(6,K)= B(6,K)+D(21)*B(7,K)
58 G(7,K)= B(7,K)
T=TP
T2=T*T
IF(NCL) T2=T
TVAL=DABS(T)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c: 25/11/98
c: here are the OUTPUT_STEP calls: Then, as in this BERND3, I interrupt the integration
c: if eccentricity reaches 0.7 (in this case iemax = 222), make RA15 return
IF(fixed.and.DABS(DIR*TM-OUT+OS).le.1.d-8) call
+OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
if(iemax.eq.222) return
c: if iemax = 222, I can stop current integration and restart new orbit.
IF (fixed.or.(out-os-dir*tm).gt.1.d-8)
+go to 246
call OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
Out=out+os
if(iemax.eq.222) return
c: if iemax = 222, stop the current integration and start new orbit
246 continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C Loop 175 is 6 iterations on first sequence and two iterations therafter.
DO 175 M=1,NI
C Loop 174 is for each substep within a sequence.
DO 174 J=2,8
JD=J-1
JDM=J-2
S=H(J)
Q=S
IF(NCL) Q=ONE
C Use Eqs. (2.9) and (2.10) of text to predict positions at each aubstep.
C These collapsed series are broken into two parts because an otherwise
C excellent compiler could not handle the complicated expression.
DO 130 K=1,NV
A=W(3)*B(3,K)+S*(W(4)*B(4,K)+S*(W(5)*B(5,K)+S*(W(6)*B(6,K)+
V S*W(7)*B(7,K))))
Y(K)=X(K)+Q*(T*V(K)+T2*S*(F1(K)*W1+S*(W(1)*B(1,K)+S*(W(2)*B(2,K)
X +S*A))))
IF(NPQ) GO TO 130
C Next are calculated the velocity predictors need for general class II.
A=U(3)*B(3,K)+S*(U(4)*B(4,K)+S*(U(5)*B(5,K)+S*(U(6)*B(6,K)+
T S*U(7)*B(7,K))))
Z(K)=V(K)+S*T*(F1(K)+S*(U(1)*B(1,K)+S*(U(2)*B(2,K)+S*A)))
130 CONTINUE
C Find forces at each substep.
CALL FORCE(Y,Z,TM+S*T,FJ)
NF=NF+1
DO 171 K=1,NV
C Find G-value for the force FJ found at the current substep. This
C section, including the many-branched GOTO, uses Eq. (2.4) of text.
TEMP=G(JD,K)
GK=(FJ(K)-F1(K))/S
GO TO (102,102,103,104,105,106,107,108),J
102 G(1,K)=GK
GO TO 160
103 G(2,K)=(GK-G(1,K))*R(1)
GO TO 160
104 G(3,K)=((GK-G(1,K))*R(2)-G(2,K))*R(3)
GO TO 160
105 G(4,K)=(((GK-G(1,K))*R(4)-G(2,K))*R(5)-G(3,K))*R(6)
GO TO 160
106 G(5,K)=((((GK-G(1,K))*R(7)-G(2,K))*R(8)-G(3,K))*R(9)-
X G(4,K))*R(10)
GO TO 160
107 G(6,K)=(((((GK-G(1,K))*R(11)-G(2,K))*R(12)-G(3,K))*R(13)-
X G(4,K))*R(14)-G(5,K))*R(15)
GO TO 160
108 G(7,K)=((((((GK-G(1,K))*R(16)-G(2,K))*R(17)-G(3,K))*R(18)-
X G(4,K))*R(19)-G(5,K))*R(20)-G(6,K))*R(21)
C Upgrade all B-values.
160 TEMP=G(JD,K)-TEMP
B(JD,K)=B(JD,K)+TEMP
C TEMP is now the improvement on G(JD,K) over its former value.
C Now we upgrade the B-value using this dfference in the one term.
C This section is based on Eq. (2.5).
GO TO (171,171,203,204,205,206,207,208),J
203 B(1,K)=B(1,K)+C(1)*TEMP
GO TO 171
204 B(1,K)=B(1,K)+C(2)*TEMP
B(2,K)=B(2,K)+C(3)*TEMP
GO TO 171
205 B(1,K)=B(1,K)+C(4)*TEMP
B(2,K)=B(2,K)+C(5)*TEMP
B(3,K)=B(3,K)+C(6)*TEMP
GO TO 171
206 B(1,K)=B(1,K)+C(7)*TEMP
B(2,K)=B(2,K)+C(8)*TEMP
B(3,K)=B(3,K)+C(9)*TEMP
B(4,K)=B(4,K)+C(10)*TEMP
GO TO 171
207 B(1,K)=B(1,K)+C(11)*TEMP
B(2,K)=B(2,K)+C(12)*TEMP
B(3,K)=B(3,K)+C(13)*TEMP
B(4,K)=B(4,K)+C(14)*TEMP
B(5,K)=B(5,K)+C(15)*TEMP
GO TO 171
208 B(1,K)=B(1,K)+C(16)*TEMP
B(2,K)=B(2,K)+C(17)*TEMP
B(3,K)=B(3,K)+C(18)*TEMP
B(4,K)=B(4,K)+C(19)*TEMP
B(5,K)=B(5,K)+C(20)*TEMP
B(6,K)=B(6,K)+C(21)*TEMP
171 CONTINUE
174 CONTINUE
IF(NES.OR.M.LT.NI) GO TO 175
C Integration of sequence is over. Next is sequence size control.
HV=ZERO
DO 635 K=1,NV
635 HV=DMAX1(HV,DABS(B(7,K)))
HV=HV*W(7)/TVAL**7
175 CONTINUE
IF (NSF) GO TO 180
IF(.NOT.NES) TP=(SS/HV)**PW*DIR
IF(NES) TP=XL
IF(NES) GO TO 170
IF(TP/T.GT.ONE) GO TO 170
8 FORMAT (A,2X,2I2,2D18.10)
TP=.8D0*TP
NCOUNT=NCOUNT+1
IF(NCOUNT.GT.10) RETURN
C Restart with 0.8x sequence size if new size called for is smaller than
C originally chosen starting sequence size on first sequence.
GO TO 4000
170 NSF=.TRUE.
C Loop 35 finds new X and V values at end of sequence using Eqs. (2.11),(2.12)
180 DO 35 K=1,NV
X(K)=X(K)+V(K)*T+T2*(F1(K)*W1+B(1,K)*W(1)+B(2,K)*W(2)+B(3,K)*W(3)
X +B(4,K)*W(4)+B(5,K)*W(5)+B(6,K)*W(6)+B(7,K)*W(7))
IF(NCL) GO TO 35
V(K)=V(K)+T*(F1(K)+B(1,K)*U(1)+B(2,K)*U(2)+B(3,K)*U(3)
V +B(4,K)*U(4)+B(5,K)*U(5)+B(6,K)*U(6)+B(7,K)*U(7))
35 CONTINUE
TM=TM+T
NS=NS+1
C Return if done.
IF(.NOT.NPER) GO TO 78
CALL OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,NPER,iemax)
RETURN
C Control on size of next sequence and adjust last sequence to exactly
C cover the integration span. NPER=.TRUE. set on last sequence.
78 CALL FORCE (X,V,TM,F1)
NF=NF+1
IF(NES) GO TO 341
TP=DIR*(SS/HV)**PW
IF(TP/T.GT.SR) TP=T*SR
341 IF(NES) TP=XL
IF(DIR*(TM+TP).LT.DIR*TF-1.D-8) GO TO 66
TP=TF-TM
NPER=.TRUE.
66 IF (.not.fixed.or.DIR*(TM+TP).LT.OUT-1.D-8) GO TO 77
TP=DIR*OUT-TM
OUT=OUT+OS
C Now predict B-values for next step. The predicted values from the preceding
C sequence were saved in the E-matrix. Te correction BD between the actual
C B-values found and these predicted values is applied in advance to the
C next sequence. The gain in accuracy is significant. Using Eqs. (2.13):
77 Q=TP/T
DO 39 K=1,NV
IF(NS.EQ.1) GO TO 31
DO 20 J=1,7
20 BD(J,K)=B(J,K)-E(J,K)
31 E(1,K)= Q*(B(1,K)+ 2.D0*B(2,K)+ 3.D0*B(3,K)+
X 4.D0*B(4,K)+ 5.D0*B(5,K)+ 6.D0*B(6,K)+ 7.D0*B(7,K))
E(2,K)= Q**2*(B(2,K)+ 3.D0*B(3,K)+
Y 6.D0*B(4,K)+10.D0*B(5,K)+15.D0*B(6,K)+21.D0*B(7,K))
E(3,K)= Q**3*(B(3,K)+
Z 4.D0*B(4,K)+10.D0*B(5,K)+20.D0*B(6,K)+35.D0*B(7,K))
E(4,K)= Q**4*(B(4,K)+ 5.D0*B(5,K)+15.D0*B(6,K)+35.D0*B(7,K))
E(5,K)= Q**5*(B(5,K)+ 6.D0*B(6,K)+21.D0*B(7,K))
E(6,K)= Q**6*(B(6,K)+ 7.D0*B(7,K))
E(7,K)= Q**7*B(7,K)
DO 39 L=1,7
39 B(L,K)=E(L,K)+BD(L,K)
C Two iterations for every sequence after the first.
NI=2
GO TO 722
END
