Trailing-Edge
-
PDP-10 Archives
-
decuslib10-01
-
43,50144/ampbx.f4
There are no other files named ampbx.f4 in the archive.
C THE PROGRAM AMPBD.F4 IS A SAMPLE PROGRAM, WHEN USED IN
C CONJUNCTION AMPBX(LOADED AND SAVED AS AMPBD.SAV) WILL
C PRODUCE A RUNNING EXAMPLE OF THE USE OF AMPBX.
SUBROUTINE AMPB2(IND, TEMP,X,DX,Y,F,N,ICOUNT,NITE
+ R,MTST)
DIMENSION TEMP(500),Y(500),F(500),AK(2)
C INTEGRATE TO NEXT POINT
401 M=N+1
IF (ICOUNT-1) 5001, 406, 411
406 N1=NITER+1
AK(1)=.92962963
AK(2)=-.070370370
C RESTORE PREVIOUS VALUES
411 X=TEMP(2*M+1)
DO 421 I=2,M
IP2M=I+2*M
Y(I-1)=TEMP(IP2M)
IP3M=IP2M+M
421 F(I-1)=TEMP(IP3M)
IF (IND) 431, 5001, 501
431 IF (ICOUNT-4) 601, 701, 701
C REDUCE DX
501 DX=DX/2.0
ICOUNT=1
C INTEGRATION BY RUNGE KUTTA
601 CALL RKPB2( TEMP,X,DX,Y,F,N)
GO TO 1001
C INTEGRATION BY ADAMS-MOULTON
C INDEPENDENT VARIABLE
701 X=X+DX
TEMP(M+1)=X
C PREDICTOR VALUES
801 DO 821 I=2,M
IPM=I+M
IP2M=IPM+M
IP3M=IP2M+M
IP4M=IP3M+M
IP5M=IP4M+M
IP6M=IP5M+M
Y(I-1)=TEMP(IP2M)+DX*(55.*TEMP(IP3M)-59.*TEMP(IP4M)+37
+.*TEMP(IP5M)-9.*TEMP(IP6M))/24.
TEMP(IPM)=Y(I-1)
IF (MTST) 811, 821,811
811 IF(ICOUNT-4) 5001,901,821
821 Y(I-1)=Y(I-1)+AK(1)/AK(2)*TEMP(1)
C CORRECTOR VALUES
901 DO 991 J1=1,N1
CALL DERIV
DO 921 I=2,M
IPM=I+M
IP2M=IPM+M
IP3M=IP2M+M
IP4M=IP3M+M
IP5M=IP4M+M
N=I-1
Y(N)=TEMP(IP2M)+DX*(9.*F(I-1)+19.*TEMP(IP3M)-5.*TEMP(I
+ P4M)
+ +TEMP(IP5M))/24.
TEMP(I)=AK(2)*(Y(N)-TEMP(IPM))
IF (MTST) 911, 921,911
911 Y(N)=Y(N)+TEMP(I)
921 CONTINUE
991 CONTINUE
C RESTORE NORMAL MODE
1001 IND=-1
5001 RETURN
END
C AMPB1************ADAMS-MOULTON INTEGRATION-1**********
+ ***
SUBROUTINE AMPB1(IND, TEMP,X,DX,Y,F,N,ICOUNT,NITE
+ R,MTST)
DIMENSION TEMP(500),Y(500),F(500)
1 IF (IND) 101, 11, 101
11 ICOUNT=0
GO TO 301
101 M=N+1
DO 191 I=1,M
DO 121 J1=1,ICOUNT
I1=I+(ICOUNT+3-J1)*M
I2=I1+M
121 TEMP(I2)=TEMP(I1)
191 CONTINUE
201 IF (IND) 301,301,302
302 IF (ICOUNT-6) 301, 211, 5001
211 DX=2.*DX
DO 261 J1=1,3
DO 241 I=1,M
I1=I+(J1+3)*M
IF (I-1) 5001,221,222
222 I2=I+(2*J1+3)*M
TEMP(I1)=TEMP(I2)
GO TO 241
221 TEMP(I1)=TEMP(I1)*2.
241 CONTINUE
261 CONTINUE
ICOUNT=3
301 ICOUNT= MIN0(ICOUNT+1,6)
CALL RKPB1( TEMP,X,DX,Y,F,N)
1001 IND=-1
5001 RETURN
END