Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib20-01 - decus/20-0025/ampbx.for
There is 1 other file named ampbx.for in the archive. Click here to see a list.
```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,NITER,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  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(IP4M)
+    +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
SUBROUTINE AMPB1(IND,TEMP,X,DX,Y,F,N,ICOUNT,NITER,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
```