tesfatsi AT iastate.edu
NASA Home Page
For additional articles and introductory materials related to NASA, visit the NASA Home Page.
Software Release
Disclaimer:
// EXEC FORTVCLG 00000010
//FORT.SYSIN DD * 00000020
C 00000030
PROGRAM NASA 00000040
C 00000050
C NONLOCAL AUTOMATED SENSITIVITY ANALYSIS OF THE PARAMETERIZED SYSTEM 00000060
C OF EQUATIONS PSI(X,ALPHA) = 0 [COMPUTERS AND MATHEMATICS WITH APPLI- 00000070
C CATIONS, VOL. 20, NO. 2, PP. 53-65, 1990] 00000080
C BY ROBERT KALABA AND LEIGH TESFATSION 00000090
C 00000100
C LAST UPDATED: 5 DECEMBER 1989 00000110
C 00000120
IMPLICIT COMPLEX (A-H,O-Z) 00000130
REAL*4 VSS,TOL,RED 00000140
DIMENSION C(6),FVAL(6),Z(66),T(1110) 00000150
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00000160
COMMON /FLAG/ C,FVAL,IFLGS,ALFA0 00000170
COMMON /RSINIT/ VSS,NSPS,NP,MAXRS,TOL,NR,RED 00000180
C 00000190
C USER INPUTS FOR ALPHA INTEGRATION: 00000200
C 00000210
C N IS THE NUMBER OF PSI COMPONENT FUNCTIONS, ALFA0 IS THE INITIAL 00000220
C VALUE FOR THE PARAMETER ALPHA, H IS THE STEP SIZE FOR ALPHA, NSP 00000230
C IS THE NUMBER OF ALPHA INTEGRATION STEPS TO BE TAKEN, NAP IS THE 00000240
C NUMBER OF ALPHA INTEGRATION STEPS BETWEEN PRINTOUTS, MAXINT IS THE 00000250
C TOTAL NUMBER OF INTEGRATION STEPS ALLOWED IN THE PROGRAM (INCLUDING 00000260
C THE THETA-BETA INITIALIZATION PROCEDURE BELOW), C IS THE INITIAL 00000270
C GUESS FOR X, AND NIT IS THE NUMBER OF THETA-BETA ITERATIONS TAKEN 00000280
C TO IMPROVE ON C PRIOR TO THE ALPHA INTEGRATION. 00000290
C 00000300
N = 2 00000310
ALFA0 = CMPLX(0.5,0.0) 00000320
H = CMPLX(0.001,0.0) 00000330
NSP = 500 00000340
NAP = 50 00000350
MAXINT = 4000 00000360
C(1) = CMPLX(1.2,0.0) 00000370
C(2) = CMPLX(1.1,0.0) 00000380
NIT = 1 00000390
C 00000400
C N IS THE NUMBER OF COMPONENTS OF C AND FVAL WHICH ARE IN ACTUAL 00000410
C USE, NDIM IS THE NUMBER OF COMPONENTS OF Z WHICH ARE IN ACTUAL 00000420
C USE, AND NDIMT IS THE NUMBER OF COMPONENTS OF T WHICH ARE IN 00000430
C ACTUAL USE. (T STORES THE COMPONENTS OF X, THE ADJOINT BY COLUMN, 00000440
C AND THE DETERMINANT, PLUS THEIR DERIVATIVES AND WORKSPACE.) 00000450
NDIM = 1 + (N+1) + (1 + (N+1))*(N+1)/2 00000460
NDIMT = 10*(1 + N + N*N) 00000470
C 00000480
C GETTING INITIAL CONDITIONS VIA A TWO-PHASE COMPLEX CONTINUATION: 00000490
C FIRST SOLVE THE DERIVED SYSTEM OF EQUATIONS 00000500
C THETA*(PSI(X,ALFA0) - PSI(C,ALFA0)) + (1-THETA)*(X - C) = 0 00000510
C FOR THETA VARYING FROM 0 TO 1 TO GET ADJ AND DET OF THE JACOBIAN 00000520
C OF PSI AT (C,ALFA0). NEXT SOLVE THE DERIVED SYSTEM OF EQUATIONS 00000530
C (PSI(X,ALFA0) - PSI(C,ALFA0)) + BETA*(PSI(C,ALFA0)) = 0 00000540
C FOR BETA VARYING FROM 0 TO 1 TO GET APPROX SOLUTION X0 = X(ALFA0) 00000550
C FOR PSI(X,ALFA0) = 0, PLUS THE ADJ AND DET OF THE JACOBIAN OF PSI 00000560
C AT (X0,ALFA0). 00000570
C 00000580
C USER INPUTS FOR THE THETA AND BETA INTEGRATION PATHS THROUGH 00000590
C THE COMPLEX SPIDER-WEB GRID: 00000600
C 00000610
C VSS IS THE STEP SIZE ALONG SPOKES, NSPS IS THE TOTAL NUMBER OF 00000620
C SPOKE-STEPS TO BE TAKEN (WITH VSS*NSPS = 1), NP IS THE NUMBER OF 00000630
C SPOKE-STEPS BETWEEN DETAILED PRINTOUTS, AND MAXRS IS THE MAXIMUM 00000640
C NUMBER OF ONE-DEGREE STEPS TO BE TAKEN ALONG ANY ONE RIM. 00000650
C 00000660
VSS = 0.005 00000670
NSPS = 200 00000680
NP = 20 00000690
MAXRS = 360 00000700
C 00000710
C TOL IS THE MINIMUM TOLERABLE ABSOLUTE VALUE FOR THE JACOBIAN OF 00000720
C THE DETERMINANT OF THE SYSTEM, NR IS THE NUMBER OF ONE-DEGREE RIM 00000730
C STEPS ATTEMPTED BEFORE ANOTHER TOLERANCE TEST AFTER A TOLERANCE 00000740
C TEST FAILURE ALONG A SPOKE, AND RED IS THE PERCENTAGE REDUCTION 00000750
C FACTOR USED TO REDUCE TOL WHEN NR RIM-STEPS OR A SPOKE-STEP AWAY 00000760
C FROM THE CURRENT POINT EACH RESULT IN A TOLERANCE TEST FAILURE, OR 00000770
C THE NUMBER OF STEPS ALONG ANY ONE RIM REACHES MAXRS. 00000780
C 00000790
TOL = .007 00000800
NR = 1 00000810
RED = .9 00000820
C 00000830
C THE REST OF THIS MAIN PROGRAM IS NOT PROBLEM-SPECIFIC. 00000840
C 00000850
C OTHER VARIABLE DEFINITIONS AND INITIALIZATIONS 00000860
NP1 = N + 1 00000870
NDE = 1 + N + N*N 00000880
NVAR = NP1 00000890
NSEC = (NP1 +1)*NP1/2 00000900
INT = 0 00000910
C 00000920
C THE VECTOR T WILL FIRST BE FILLED WITH THE INITIAL GUESS C FOR X, 00000930
C THE ADJOINT I OF THE JACOBIAN OF (X - C) BY COLUMN, AND THE DET 1 00000940
C OF THE JACOBIAN OF (X - C). 00000950
C 00000960
LL = 1 00000970
DO 10 I = 1,N 00000980
T(LL) = C(LL) 00000990
LL = LL + 1 00001000
10 CONTINUE 00001010
DO 1000 IT = 1,NIT 00001020
DO 20 J = 1,N 00001030
DO 30 I = 1,N 00001040
T(LL) = CMPLX(0.0,0.0) 00001050
IF (J.EQ.I) T(LL) = CMPLX(1.0,0.0) 00001060
LL = LL + 1 00001070
30 CONTINUE 00001080
20 CONTINUE 00001090
T(LL) = CMPLX(1.0,0.0) 00001100
C 00001110
C EVALUATION OF PSI(C,ALFA0) 00001120
IFLGS = 2 00001130
ALFCUR = ALFA0 00001140
DO 40 K = 1,N 00001150
CALL FUN(K,C,ALFCUR,Z) 00001160
FVAL(K) = Z(1) 00001170
40 CONTINUE 00001180
WRITE (6,111) 00001190
WRITE (6,112) (C(K),K=1,N) 00001200
111 FORMAT ('1',//1X,'HERE IS THE INITIAL GUESS C FOR X') 00001210
112 FORMAT (1X,3X,4E12.3) 00001220
WRITE (6,113) 00001230
WRITE (6,112) (FVAL(K),K=1,N) 00001240
113 FORMAT (1X,'HERE ARE THE COMPONENT VALUES FOR PSI(C,ALFA0)') 00001250
C 00001260
C THETA AND BETA INTEGRATIONS (IFLGS = 0,1) 00001270
IFLGS = 0 00001280
CALL INIT(T) 00001290
IFLGS = 1 00001300
CALL INIT(T) 00001310
IF (IT.LT.NIT) THEN 00001320
LL = 1 00001330
DO 50 IJ = 1,N 00001340
C(LL) = T(LL) 00001350
LL = LL + 1 00001360
50 CONTINUE 00001370
ENDIF 00001380
1000 CONTINUE 00001390
C 00001400
C ALPHA INTEGRATION (IFLGS = 2) 00001410
IFLGS = 2 00001420
ALFCUR = ALFA0 00001430
BETA = CMPLX(1.0,0.0) 00001440
THETA = CMPLX(1.0,0.0) 00001450
CALL CINT1(T,NDE,ALFCUR,H) 00001460
CALL OUTPUT(T,ALFCUR,BETA,THETA) 00001470
JJ = NSP/NAP 00001480
DO 60 J = 1,JJ 00001490
DO 70 I = 1,NAP 00001500
CALL CINT2(T,NDE,ALFCUR,H) 00001510
INT = INT + 1 00001520
IF (INT.GE.MAXINT) STOP 'ALPHA INTEGRATION' 00001530
70 CONTINUE 00001540
CALL OUTPUT(T,ALFCUR,BETA,THETA) 00001550
60 CONTINUE 00001560
WRITE (6,114) INT 00001570
114 FORMAT (1X,//1X,'THE TOTAL NUMBER OF INTEGRATION STEPS =',I7) 00001580
STOP 00001590
END 00001600
C 00001610
C 00001620
SUBROUTINE PSIFUN(K,X,VCUR,VC,XX,PSIK) 00001630
IMPLICIT COMPLEX (A-H,O-Z) 00001640
DIMENSION X1(66),X2(66),X3(66),X4(66),X5(66),X6(66) 00001650
DIMENSION X(6),C(6),FVAL(6),VC(66),XX(6,66),PSIK(66), 00001660
*ALFCUR(66) 00001670
DIMENSION U1(66),U2(66),U3(66),U4(66),U5(66) 00001680
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00001690
COMMON /FLAG/ C,FVAL,IFLGS,ALFA0 00001700
C 00001710
C THIS SUBROUTINE OBTAINS THE VALUE AND PARTIAL DERIVATIVES FOR THE 00001720
C THE FOLLOWING PROBLEM-SPECIFIC COMPONENTS OF PSI: 00001730
C 00001740
C PSI1(X,ALPHA) = (1/2)(X1**(-1/2))*(X2**(1/3)) - ALPHA 00001750
C PSI2(X,ALPHA) = (1/3)*(X1**(1/2))*(X2**(-2/3)) - (1/3) 00001760
C 00001770
C NOT PROBLEM SPECIFIC: THE INDEPENDENT VARIABLES X AND VCUR 00001780
C ARE VECTORIZED AND ALPHA IS INITIALIZED 00001790
C 00001800
DO 1 I = 1,N 00001810
ZZ = X(I) 00001820
CALL VEC(I,ZZ,U1) 00001830
DO 5 J = 1,NDIM 00001840
XX(I,J) = U1(J) 00001850
5 CONTINUE 00001860
1 CONTINUE 00001870
CALL VEC(NP1,VCUR,VC) 00001880
IF (IFLGS.EQ.0.OR.IFLGS.EQ.1) THEN 00001890
CALL CON(ALFA0,U1) 00001900
DO 10 J = 1,NDIM 00001910
ALFCUR(J) = U1(J) 00001920
10 CONTINUE 00001930
ELSEIF (IFLGS.EQ.2) THEN 00001940
DO 15 J = 1,NDIM 00001950
ALFCUR(J) = VC(J) 00001960
15 CONTINUE 00001970
ENDIF 00001980
C 00001990
C PROBLEM-SPECIFIC STEPS START HERE 00002000
C 00002010
C PUT THE N ROWS OF XX INTO VECTORS X1,...,XN 00002020
DO 20 J = 1,NDIM 00002030
X1(J) = XX(1,J) 00002040
X2(J) = XX(2,J) 00002050
20 CONTINUE 00002060
C 00002070
C SPECIFY ANY PROBLEM-SPECIFIC CONSTANTS USED TO DEFINE PSI 00002080
P1 = CMPLX(0.5,0.0) 00002090
P2 = CMPLX(-0.5,0.0) 00002100
P3 = CMPLX(1.0,0.0)/CMPLX(3.0,0.0) 00002110
P4 = CMPLX(-2.0,0.0)/CMPLX(3.0,0.0) 00002120
C 00002130
IF (K.EQ.1) THEN 00002140
C 00002150
C GIVE THE PROBLEM-SPECIFIC LIST OF CALCULUS SUBROUTINES 00002160
C FOR PSI1(X,ALFCUR) USING X = (X1,...,XN) AND PUT THE 00002170
C FINAL CALCULATION INTO PSIK 00002180
CALL POW(X1,P2,U1) 00002190
CALL CONVE(P1,U1,U2) 00002200
CALL POW(X2,P3,U3) 00002210
CALL MULT(U2,U3,U4) 00002220
CALL SUB(U4,ALFCUR,PSIK) 00002230
RETURN 00002240
C 00002250
ELSEIF (K.EQ.2) THEN 00002260
C 00002270
C GIVE THE PROBLEM-SPECIFIC LIST OF CALCULUS SUBROUTINES 00002280
C FOR PSI2(X,ALFCUR) USING X = (X1,...,XN) AND PUT THE 00002290
C FINAL CALCULATION INTO PSIK 00002300
CALL POW(X1,P1,U1) 00002310
CALL CONVE(P3,U1,U2) 00002320
CALL POW(X2,P4,U3) 00002330
CALL MULT(U2,U3,U4) 00002340
CALL CON(P3,U5) 00002350
CALL SUB(U4,U5,PSIK) 00002360
RETURN 00002370
C 00002380
END IF 00002390
RETURN 00002400
END 00002410
C 00002420
C 00002430
SUBROUTINE OUTPUT(T,ALFCUR,BETA,THETA) 00002440
IMPLICIT COMPLEX (A-H,O-Z) 00002450
DIMENSION T(1110),C(6),FVAL(6),X(6),Z(66) 00002460
DIMENSION PSI(6),PTHETA(6),PBETA(6) 00002470
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00002480
COMMON /FLAG/ C,FVAL,IFLGS,ALFA0 00002490
COMMON /RSINIT/ VSS,NSPS,NP,MAXRS,TOL,NR,RED 00002500
C 00002510
C THIS SUBROUTINE PRINTS OUTPUT FOR THE SYSTEM OF EQUATIONS 00002520
C PSI(X,ALPHA) = 0, AND FOR THE THETA AND BETA CONTINUATIONS OF 00002530
C PSI, WITH COMPONENT FUNCTIONS DEFINED AS IN SUBROUTINE PSIFUN 00002540
C AND SUBROUTINE FUN. 00002550
C 00002560
NDEP1 = NDE + 1 00002570
NDEPN = NDE + N 00002580
NDEMT = 2*NDE 00002590
C 00002600
C HERE ARE THE COMPONENTS OF PSI(X,ALFCUR) 00002610
DO 10 I = 1,N 00002620
X(I) = T(I) 00002630
10 CONTINUE 00002640
IFTEM = IFLGS 00002650
IFLGS = 2 00002660
DO 20 K = 1,N 00002670
CALL FUN(K,X,ALFCUR,Z) 00002680
PSI(K) = Z(1) 00002690
20 CONTINUE 00002700
IFLGS = IFTEM 00002710
C 00002720
C HERE ARE THE COMPONENTS OF PTHETA(X,THETA) = 00002730
C THETA*(PSI(X,ALFA0) - PSI(C,ALFA0) + (1 - THETA)*(X - C) 00002740
C1 = CMPLX(1.0,0.0) 00002750
DO 30 K = 1,N 00002760
PTHETA(K) = THETA*(PSI(K)-FVAL(K))+(C1-THETA)*(T(K)-C(K)) 00002770
30 CONTINUE 00002780
C 00002790
C HERE ARE THE COMPONENTS OF PBETA(X,BETA) = 00002800
C (PSI(X,ALFA0) - PSI(C,ALFA0)) + BETA*PSI(C,ALFA0) 00002810
DO 40 K = 1,N 00002820
PBETA(K) = (PSI(K)-FVAL(K)) + BETA*FVAL(K) 00002830
40 CONTINUE 00002840
C 00002850
IF (IFLGS.EQ.0) THEN 00002860
WRITE (6,100) THETA 00002870
100 FORMAT (1X, //1X,'THETA =',2E12.3) 00002880
WRITE (6,110) 00002890
110 FORMAT (1X ,'HERE ARE THE COMPONENTS OF PTHETA') 00002900
WRITE (6,120) (PTHETA(K),K=1,N) 00002910
120 FORMAT (1X,3X,4E18.6) 00002920
WRITE (6,130) 00002930
130 FORMAT (1X,'HERE ARE THE COMPONENTS OF X AND THE SYSTEM DET') 00002940
WRITE (6,120) (T(I),I=1,N) , T(NDE) 00002950
WRITE (6,140) 00002960
140 FORMAT (1X,'HERE ARE THE DERIVATIVES OF X AND THE DET ', 00002970
* 'WRT THETA') 00002980
WRITE (6,120) (T(J),J=NDEP1,NDEPN), T(NDEMT) 00002990
ELSEIF (IFLGS.EQ.1) THEN 00003000
WRITE (6,200) BETA 00003010
200 FORMAT (1X,//1X,'BETA =',2E12.3) 00003020
WRITE (6,210) 00003030
210 FORMAT (1X,'HERE ARE THE COMPONENTS OF PBETA') 00003040
WRITE (6,120) (PBETA(K),K=1,N) 00003050
WRITE (6,130) 00003060
WRITE (6,120) (T(I),I=1,N), T(NDE) 00003070
WRITE (6,240) 00003080
240 FORMAT (1X,'HERE ARE THE DERIVATIVES OF X AND THE DET ', 00003090
* 'WRT BETA') 00003100
WRITE (6,120) (T(J),J=NDEP1,NDEPN),T(NDEMT) 00003110
ELSEIF (IFLGS.EQ.2) THEN 00003120
WRITE (6,300) ALFCUR 00003130
300 FORMAT (1X, //1X,'ALFA =',2E12.3) 00003140
WRITE (6,310) 00003150
310 FORMAT (1X, 'HERE ARE THE COMPONENTS OF PSI') 00003160
WRITE (6,120) (PSI(K),K=1,N) 00003170
WRITE (6,130) 00003180
WRITE (6,120) (T(I),I=1,N), T(NDE) 00003190
WRITE (6,340) 00003200
340 FORMAT (1X,'HERE ARE THE DERIVATIVES OF X AND THE DET ', 00003210
* 'WRT ALFA') 00003220
WRITE (6,120) (T(J),J=NDEP1,NDEPN), T(NDEMT) 00003230
ENDIF 00003240
RETURN 00003250
END 00003260
C 00003270
SUBROUTINE FUN(K,X,VCUR,Z) 00003280
IMPLICIT COMPLEX (A-H,O-Z) 00003290
DIMENSION X(6),Z(66),C(6),FVAL(6),VC(66),XX(6,66),PSIK(66), 00003300
*XK(66) 00003310
DIMENSION U1(66),U2(66),U3(66),U4(66),U5(66),U6(66),U7(66), 00003320
*U8(66) 00003330
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00003340
COMMON /FLAG/ C,FVAL,IFLGS,ALFA0 00003350
C 00003360
C THIS SUBROUTINE OBTAINS THE VALUE AND PARTIAL DERIVATIVES FOR 00003370
C THE THETA AND BETA CONTINUATIONS PTHETA AND PBETA FOR ANY GIVEN 00003380
C PROBLEM-SPECIFIC FUNCTION PSI GIVEN IN SUBROUTINE PSIFUN, WHERE 00003390
C PTHETAK(X,THETA) = THETA*(PSIK(X,ALFA0)-PSIK(C,ALFA0)) 00003400
C + (1 - THETA)*(XK - CK) 00003410
C PBETAK(X,BETA) = (PSIK(X,ALFA0)-PSIK(C,ALFA0) + BETA*PSIK(C,ALFA0) 00003420
C 00003430
CALL PSIFUN(K,X,VCUR,VC,XX,PSIK) 00003440
IF (IFLGS.EQ.2) THEN 00003450
DO 1111 J = 1,NDIM 00003460
Z(J) = PSIK(J) 00003470
1111 CONTINUE 00003480
RETURN 00003490
ENDIF 00003500
C 00003510
FVK = FVAL(K) 00003520
C 00003530
IF (IFLGS.EQ.0) THEN 00003540
C 00003550
B1 = CMPLX(1.0,0.0) 00003560
CK = C(K) 00003570
DO 2222 I = 1,NDIM 00003580
XK(I) = XX(K,I) 00003590
2222 CONTINUE 00003600
C 00003610
C HERE IS THETA*(PSIK(X,ALFA0)-PSIK(C,ALFA0)) + (1-THETA)*(XK-CK) 00003620
CALL CON(FVK,U1) 00003630
CALL SUB(PSIK,U1,U2) 00003640
CALL MULT(VC,U2,U3) 00003650
CALL CON(B1,U4) 00003660
CALL SUB(U4,VC,U5) 00003670
CALL CON(CK,U6) 00003680
CALL SUB(XK,U6,U7) 00003690
CALL MULT(U5,U7,U8) 00003700
CALL ADD(U3,U8,Z) 00003710
RETURN 00003720
C 00003730
ELSEIF (IFLGS.EQ.1) THEN 00003740
C 00003750
C HERE IS (PSIK(X,ALFA0)-PSIK(C,ALFA0)) + BETA*PSIK(C,ALFA0) 00003760
CALL CON(FVK,U1) 00003770
CALL MULT(VC,U1,U2) 00003780
CALL ADD(PSIK,U2,U3) 00003790
CALL SUB(U3,U1,Z) 00003800
RETURN 00003810
C 00003820
ENDIF 00003830
RETURN 00003840
END 00003850
C 00003860
C 00003870
SUBROUTINE INIT(T) 00003880
IMPLICIT COMPLEX (A-H,O-Z) 00003890
REAL*4 VSS,TOL,RED,RADIUS,ANGLE,DEGREE,PI,SADET,GSADET 00003900
REAL*4 RNUM,MNUM,UANGLE,DANGLE,GANGLE,UADET,DADET 00003910
DIMENSION T(1110),C(6),FVAL(6),STNEW(1110),GSTNEW(1110) 00003920
DIMENSION UTNEW(1110),DTNEW(1110) 00003930
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00003940
COMMON /FLAG/ C,FVAL,IFLGS,ALFA0 00003950
COMMON /RSINIT/ VSS,NSPS,NP,MAXRS,TOL,NR,RED 00003960
C 00003970
C THIS SUBROUTINE CARRIES OUT THE INTEGRATIONS FOR THE THETA AND 00003980
C BETA CONTINUATIONS TO OBTAIN THE INITIAL CONDITIONS NEEDED FOR 00003990
C THE ALPHA INTEGRATION STARTING FROM ALPHA = ALFA0. 00004000
C 00004010
PI = 4.0*ATAN(1.0) 00004020
DEGREE = PI/180.0 00004030
C 00004040
C VARIABLE INITIALIZATIONS 00004050
NRIM = 0 00004060
IFDIR = 0 00004070
NRSTEP = 0 00004080
NSPOKE = 0 00004090
RADIUS = 1.0 00004100
ANGLE = 0.0 00004110
GSADET = 0.0 00004120
GANGLE = ANGLE 00004130
C 00004140
IF (IFLGS.EQ.0) THEN 00004150
THETA = CMPLX(0.0,0.0) 00004160
BETA = CMPLX(0.0,0.0) 00004170
VCUR = THETA 00004180
ELSEIF (IFLGS.EQ.1) THEN 00004190
DO 1 K = 1,N 00004200
T(K) = C(K) 00004210
1 CONTINUE 00004220
THETA = CMPLX(1.0,0.0) 00004230
BETA = CMPLX(0.0,0.0) 00004240
VCUR = BETA 00004250
ENDIF 00004260
HS = CMPLX(VSS,0.0) 00004270
CALL CINT1(T,NDE,VCUR,HS) 00004280
GSHS = HS 00004290
CALL OUTPUT(T,ALFA0,BETA,THETA) 00004300
IF (IFLGS.EQ.0) THEN 00004310
WRITE (6,211) NRIM,NSPOKE,NRSTEP,VCUR 00004320
ELSEIF (IFLGS.EQ.1) THEN 00004330
WRITE (6,221) NRIM,NSPOKE,NRSTEP,VCUR 00004340
ENDIF 00004350
WRITE (6,231) TOL,T(NDE),RADIUS 00004360
C 00004370
C FIRST TEST TO SEE WHETHER A STEP CAN BE MADE INWARD ALONG THE 00004380
C CURRENT SPOKE AT THE CURRENT TOLERANCE LEVEL TOL, I.E., WHETHER 00004390
C THE ABSOLUTE VALUE SADET OF THE DETERMINANT OF THE JACOBIAN 00004400
C IS AT LEAST AS GREAT AS TOL WHEN THIS SPOKE STEP IS TAKEN. 00004410
C 00004420
2000 IF (NRIM.LT.NSPS) THEN 00004430
DO 10 IM = 1,NDIMT 00004440
STNEW(IM) = T(IM) 00004450
10 CONTINUE 00004460
SVNEW = VCUR 00004470
CALL CINT2(STNEW,NDE,SVNEW,HS) 00004480
INT = INT + 1 00004490
IF (INT.GE.MAXINT) STOP 'INITIALIZATION 10' 00004500
DET = STNEW(NDE) 00004510
SADET = CABS(DET) 00004520
C IF SADET IS GREATER OR EQUAL TO TOL, MOVE ALONG THE CURRENT SPOKE 00004530
C TO THE NEXT RIM. 00004540
IF (SADET.GE.TOL) THEN 00004550
DO 20 ID = 1,NDIMT 00004560
T(ID) = STNEW(ID) 00004570
20 CONTINUE 00004580
VCUR = SVNEW 00004590
RADIUS = RADIUS - VSS 00004600
NRIM = NRIM + 1 00004610
IFDIR = 0 00004620
NRSTEP = 0 00004630
GSADET = 0.0 00004640
GSHS = HS 00004650
GANGLE = ANGLE 00004660
ENDIF 00004670
C IF SADET IS LESS THAN TOL, AND THE NUMBER OF RIM-STEPS PREVIOUSLY 00004680
C TAKEN ALONG THE CURRENT RIM IS LESS THAN MAXRS, TEST TO SEE IF NR 00004690
C RIM-STEPS EITHER UPWARD OR DOWNWARD ALONG THE RIM ARE POSSIBLE AT 00004700
C THE CURRENT TOLERANCE LEVEL. IF NOT, MOVE TO THE TESTED SPOKE OR 00004710
C OR RIM POINT WITH THE GREATEST ABSOLUTE VALUE ABS(DET) FOR THE 00004720
C DET OF THE JACOBIAN AND REDUCE TOL TO RED*(GREATEST ABS(DET)). 00004730
IF (SADET.LT.TOL.AND.NRSTEP.LT.MAXRS) THEN 00004740
IF (SADET.GT.GSADET) THEN 00004750
GSADET = SADET 00004760
GSHS = HS 00004770
GANGLE = ANGLE 00004780
DO 30 LT = 1,NDIMT 00004790
GSTNEW(LT) = STNEW(LT) 00004800
30 CONTINUE 00004810
GSVNEW = SVNEW 00004820
ENDIF 00004830
UADET = 0.0 00004840
IF (IFDIR.GE.0) THEN 00004850
DO 40 LK = 1,NDIMT 00004860
UTNEW(LK) = T(LK) 00004870
40 CONTINUE 00004880
UVNEW = VCUR 00004890
UANGLE = ANGLE 00004900
DO 50 JJJ = 1,NR 00004910
UANGLE = UANGLE + DEGREE 00004920
RNUM = 1.0 - RADIUS*COS(UANGLE) 00004930
MNUM = RADIUS*SIN(UANGLE) 00004940
UHS = CMPLX(RNUM,MNUM) - UVNEW 00004950
CALL CINT1(UTNEW,NDE,UVNEW,UHS) 00004960
CALL CINT2(UTNEW,NDE,UVNEW,UHS) 00004970
INT = INT + 1 00004980
IF (INT.GE.MAXINT) STOP 'INITIALIZATION 50' 00004990
50 CONTINUE 00005000
DET = UTNEW(NDE) 00005010
UADET = CABS(DET) 00005020
ENDIF 00005030
DADET = 0.0 00005040
IF (IFDIR.LE.0) THEN 00005050
DO 60 IK = 1,NDIMT 00005060
DTNEW(IK) = T(IK) 00005070
60 CONTINUE 00005080
DVNEW = VCUR 00005090
DANGLE = ANGLE 00005100
DO 70 III = 1,NR 00005110
DANGLE = DANGLE - DEGREE 00005120
RNUM = 1.0 - RADIUS*COS(DANGLE) 00005130
MNUM = RADIUS*SIN(DANGLE) 00005140
DHS = CMPLX(RNUM,MNUM) - DVNEW 00005150
CALL CINT1(DTNEW,NDE,DVNEW,DHS) 00005160
CALL CINT2(DTNEW,NDE,DVNEW,DHS) 00005170
INT = INT + 1 00005180
IF (INT.GE.MAXINT) STOP 'INITIALIZATION 70' 00005190
70 CONTINUE 00005200
DET = DTNEW(NDE) 00005210
DADET = CABS(DET) 00005220
ENDIF 00005230
IF (UADET.GT.SADET.AND.UADET.GE.DADET) THEN 00005240
IF (UADET.LT.TOL) TOL = RED*UADET 00005250
DO 80 LLK = 1,NDIMT 00005260
T(LLK) = UTNEW(LLK) 00005270
80 CONTINUE 00005280
VCUR = UVNEW 00005290
ANGLE = UANGLE 00005300
NRSTEP = NRSTEP + NR 00005310
NSPOKE = NSPOKE + NR 00005320
IFDIR = 1 00005330
RNUM = VSS*COS(ANGLE) 00005340
MNUM = -VSS*SIN(ANGLE) 00005350
HS = CMPLX(RNUM,MNUM) 00005360
CALL CINT1(T,NDE,VCUR,HS) 00005370
ENDIF 00005380
IF (DADET.GT.SADET.AND.DADET.GE.UADET) THEN 00005390
IF (DADET.LT.TOL) TOL = RED*DADET 00005400
DO 90 IIK = 1,NDIMT 00005410
T(IIK) = DTNEW(IIK) 00005420
90 CONTINUE 00005430
VCUR = DVNEW 00005440
ANGLE = DANGLE 00005450
NRSTEP = NRSTEP + NR 00005460
NSPOKE = NSPOKE - NR 00005470
IFDIR = -1 00005480
RNUM = VSS*COS(ANGLE) 00005490
MNUM = -VSS*SIN(ANGLE) 00005500
HS = CMPLX(RNUM,MNUM) 00005510
CALL CINT1(T,NDE,VCUR,HS) 00005520
ENDIF 00005530
IF (SADET.GE.DADET.AND.SADET.GE.UADET) THEN 00005540
TOL = RED*SADET 00005550
DO 100 ILK = 1,NDIMT 00005560
T(ILK) = STNEW(ILK) 00005570
100 CONTINUE 00005580
VCUR = SVNEW 00005590
RADIUS = RADIUS - VSS 00005600
NRIM = NRIM + 1 00005610
IFDIR = 0 00005620
NRSTEP = 0 00005630
GSADET = 0.0 00005640
GSHS = HS 00005650
GANGLE = ANGLE 00005660
ENDIF 00005670
ENDIF 00005680
C IF MAXRS STEPS HAVE BEEN TAKEN AROUND ANY ONE RIM, DECREASE TOL 00005690
C TO RED*(GREATEST ABSOLUTE VALUE OF DET ALONG NEXT INNER RIM) 00005700
C AND MOVE TO THIS NEXT INNER RIM POINT. 00005710
IF (SADET.LT.TOL.AND.NRSTEP.GE.MAXRS) THEN 00005720
TOL = RED*GSADET 00005730
DO 110 KKL = 1,NDIMT 00005740
T(KKL) = GSTNEW(KKL) 00005750
110 CONTINUE 00005760
VCUR = GSVNEW 00005770
RADIUS = RADIUS - VSS 00005780
ANGLE = GANGLE 00005790
NRIM = NRIM + 1 00005800
IFDIR = 0 00005810
NRSTEP = 0 00005820
GSADET = 0.0 00005830
HS = GSHS 00005840
CALL CINT1(T,NDE,VCUR,HS) 00005850
ENDIF 00005860
IF ((NRIM/NP)*NP.EQ.NRIM) THEN 00005870
IF (IFLGS.EQ.0) THEN 00005880
CALL OUTPUT(T,ALFA0,BETA,VCUR) 00005890
WRITE (6,211) NRIM,NSPOKE,NRSTEP,VCUR 00005900
211 FORMAT(1X, /1X,'NRIM =',I6,2X,'NSPOKE =',I6,2X, 00005910
* 'NRSTEP =',I6,2X,'THETA =',2E12.3) 00005920
ELSEIF (IFLGS.EQ.1) THEN 00005930
CALL OUTPUT(T,ALFA0,VCUR,THETA) 00005940
WRITE (6,221) NRIM,NSPOKE,NRSTEP,VCUR 00005950
221 FORMAT(1X, /1X,'NRIM =',I6,2X,'NSPOKE =',I6,2X, 00005960
* 'NRSTEP =',I6,2X,'BETA =',2E12.3) 00005970
ENDIF 00005980
WRITE (6,231) TOL,T(NDE),RADIUS 00005990
231 FORMAT (1X, 'TOL =',E12.3,3X,'DET =',2E12.3,3X,'RADIUS =', 00006000
* F10.3) 00006010
ENDIF 00006020
GO TO 2000 00006030
ENDIF 00006040
WRITE (6,241) INT 00006050
241 FORMAT (1X,//1X, 'NUMBER OF INTEGRATION STEPS SO FAR =', I6) 00006060
RETURN 00006070
END 00006080
C 00006090
C 00006100
SUBROUTINE DAUX(T,VCUR,H) 00006110
C 00006120
C THIS DERIVATIVE AUXILIARY ROUTINE CALCULATES THE DERIVATIVES 00006130
C DX/DVCUR, DADJ/DVCUR, AND DDET/DVCUR, AND STORES THEM IN T. 00006140
C 00006150
IMPLICIT COMPLEX (A-H,O-Z) 00006160
DIMENSION T(1110), A(10,10),PSUBA(10),V(10),DXSUBA(10),B(10,10), 00006170
*D(10,10),E(10,10),F(10,10),ASUBA(10,10),X(6),G(10,10), 00006180
*Z(66),DPSI(10,11),DDPSI(10,11,11),VF(10) 00006190
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00006200
C 00006210
C DICTIONARY OF VALUES T(1) THROUGH T(NDE) 00006220
C 00006230
L = 1 00006240
DO 10 KJ=1,N 00006250
X(KJ)=T(L) 00006260
L=L+1 00006270
10 CONTINUE 00006280
DO 20 J=1,N 00006290
DO 30 I=1,N 00006300
A(I,J)=T(L) 00006310
L=L+1 00006320
30 CONTINUE 00006330
20 CONTINUE 00006340
DET=T(L) 00006350
L = L+1 00006360
C 00006370
C CALCULATE FIRST AND SECOND-ORDER PARTIAL DERIVATIVES OF PSI AND 00006380
C STORE THEM IN DPSI AND DDPSI. 00006390
C 00006400
DO 31 K=1,N 00006410
LL=1 00006420
CALL FUN(K,X,VCUR,Z) 00006430
VF(K)=Z(LL) 00006440
LL=LL+1 00006450
DO 41 I=1,NP1 00006460
DPSI(K,I)=Z(LL) 00006470
LL=LL+1 00006480
41 CONTINUE 00006490
DO 51 I=1,NP1 00006500
DO 61 J=I,NP1 00006510
DDPSI(K,I,J)=Z(LL) 00006520
DDPSI(K,J,I)=Z(LL) 00006530
LL=LL+1 00006540
61 CONTINUE 00006550
51 CONTINUE 00006560
31 CONTINUE 00006570
C 00006580
C CALCULATE DX/DVCUR = -(ADJ/DET)*(DPSI/DVCUR) 00006590
C 00006600
DO 40 K=1,N 00006610
PSUBA(K)=DPSI(K,NP1) 00006620
40 CONTINUE 00006630
CALL MATVEC(A,PSUBA,N,V) 00006640
CC = CMPLX(-1.0,0.0)/DET 00006650
CALL VECON(CC,V,N,DXSUBA) 00006660
C 00006670
C CALCULATE B = D(JACOBIAN)/DVCUR 00006680
C 00006690
DO 50 I=1,N 00006700
DO 60 J=1,N 00006710
SUM=DDPSI(I,J,NP1) 00006720
DO 70 M=1,N 00006730
SUM=SUM+DDPSI(I,J,M)*DXSUBA(M) 00006740
70 CONTINUE 00006750
B(I,J)=SUM 00006760
60 CONTINUE 00006770
50 CONTINUE 00006780
C 00006790
C CALCULATE DADJ/DVCUR = (A*TR(A*B) - A*B*A)/DET 00006800
C 00006810
CALL MULMAT(A,B,N,D) 00006820
CALL TRACE(D,N,TR) 00006830
CALL MULMAT(D,A,N,E) 00006840
CALL CONMAT(TR,A,N,F) 00006850
CALL SUBMAT(F,E,N,G) 00006860
CC = -CC 00006870
CALL CONMAT(CC,G,N,ASUBA) 00006880
C 00006890
C CALCULATE DDET/DVCUR = TR(A*B) 00006900
C 00006910
DDETA=TR 00006920
C 00006930
C MAP THE DERIVATIVES DX/DVCUR, DADJ/DVCUR (BY COLUMN), AND 00006940
C DDET/DVCUR INTO T IN ORDER, STARTING WITH T(NDE + 1). 00006950
C 00006960
DO 80 J=1,N 00006970
T(L)=DXSUBA(J) 00006980
L=L+1 00006990
80 CONTINUE 00007000
DO 81 J=1,N 00007010
DO 82 I=1,N 00007020
T(L)=ASUBA(I,J) 00007030
L=L+1 00007040
82 CONTINUE 00007050
81 CONTINUE 00007060
T(L)=DDETA 00007070
RETURN 00007080
END 00007090
C 00007100
C 00007110
SUBROUTINE CINT1(T,N,X,H) 00007120
C CALCULATES CONSTANTS NEEDED FOR INTEGRATION BY SUB CINT2 00007130
IMPLICIT COMPLEX (A-H,O-Z) 00007140
DIMENSION T(1110) 00007150
COMMON /NSEQ/ N2,N3,N4,N5,N6,N7,N8,N9,NN,KFLAG,IND 00007160
COMMON /HSEQ/ H2,H4,H24,R 00007170
NN=N 00007180
N2=N*2 00007190
N3=N*3 00007200
N4=N*4 00007210
N5=N*5 00007220
N6=N*6 00007230
N7=N*7 00007240
N8=N*8 00007250
N9=N*9 00007260
H2=H*CMPLX(0.5,0.0) 00007270
H4=H2*CMPLX(0.5,0.0) 00007280
H24=H/CMPLX(24.0,0.0) 00007290
R=CMPLX(1.0,0.0)/CMPLX(6.0,0.0) 00007300
C CALCULATE Y PRIME FOR INITIAL CONDITIONS 00007310
CALL DAUX(T,X,H) 00007320
DO 1 I=1,N 00007330
N9I=I+N9 00007340
C TEMPORARY STORAGE FOR Y 00007350
N8I=I+N8 00007360
NNI=I+NN 00007370
T(N9I)=T(I) 00007380
C STORE Y PRIME AT N-3 FOR USE IN ADAMS-MOULTON INTEGRATION 00007390
1 T(N8I)=T(NNI) 00007400
KFLAG=0 00007410
IND=0 00007420
RETURN 00007430
END 00007440
C 00007450
C 00007460
SUBROUTINE CINT2(T,N,X,H) 00007470
C FOURTH-ORDER ADAMS-MOULTON INTEGRATION WITH A RUNGE-KUTTA START. 00007480
IMPLICIT COMPLEX (A-H,O-Z) 00007490
DIMENSION T(1110) 00007500
COMMON /NSEQ/ N2,N3,N4,N5,N6,N7,N8,N9,NN,KFLAG,IND 00007510
COMMON /HSEQ/ H2,H4,H24,R 00007520
C 00007530
C IND=FLAG FOR RUNGE-KUTTA INTEGRATION WHEN EQUAL OR LESS THAN 3. 00007540
C KFLAG=FLAG FOR STORING PAST DERIVATIVES FOR USE IN ADAMS-MOULTON 00007550
C INTEGRATION 00007560
C RUNGE-KUTTA INTEGRATION (6 STEPS AT H/2) 00007570
C NOTE-TWO STEPS OF RUNGE-KUTTA ARE DONE WITH EACH CALL TO THIS 00007580
C SUBROUTINE, SO THAT PRINTOUT POINTS WILL BE AT STEPS OF H. 00007590
C 00007600
IND=IND+1 00007610
KFLAG=KFLAG+1 00007620
IF (IND-3) 20,20,12 00007630
20 DO 9 K=1,2 00007640
C CALCULATE K1 00007650
DO 1 I=1,N 00007660
N2I=I+N2 00007670
NNI=I+NN 00007680
1 T(N2I)=T(NNI)*H2 00007690
C CALCULATE K2 00007700
C STEP UP X 00007710
X=X+H4 00007720
DO 2 I=1,N 00007730
N9I=I+N9 00007740
N2I=I+N2 00007750
2 T(I)=T(N9I)+CMPLX(0.5,0.0)*T(N2I) 00007760
CALL DAUX(T,X,H) 00007770
C STORE K2 00007780
DO 3 I=1,N 00007790
N3I=I+N3 00007800
NNI=I+NN 00007810
3 T(N3I)=T(NNI)*H2 00007820
C CALCULATE K3 00007830
DO 4 I=1,N 00007840
N9I=I+N9 00007850
N3I=I+N3 00007860
4 T(I)=T(N9I)+CMPLX(0.5,0.0)*T(N3I) 00007870
CALL DAUX(T,X,H) 00007880
C STORE K3 00007890
DO 5 I=1,N 00007900
N4I=I+N4 00007910
NNI=I+NN 00007920
5 T(N4I)=T(NNI)*H2 00007930
C CALCULATE K4 00007940
C STEP UP X 00007950
X=X+H4 00007960
DO 6 I=1,N 00007970
N9I=I+N9 00007980
N4I=I+N4 00007990
6 T(I)=T(N9I)+T(N4I) 00008000
CALL DAUX(T,X,H) 00008010
C STORE K4 00008020
DO 7 I=1,N 00008030
N5I=I+N5 00008040
NNI=I+NN 00008050
7 T(N5I)=T(NNI)*H2 00008060
C CALCULATE PREDICTED VALUE OF Y 00008070
DO 8 I=1,N 00008080
N9I=I+N9 00008090
N2I=I+N2 00008100
N3I=I+N3 00008110
N4I=I+N4 00008120
N5I=I+N5 00008130
Y=T(N9I)+R*(T(N2I)+CMPLX(2.0,0.0)*T(N3I) 00008140
*+CMPLX(2.0,0.0)*T(N4I)+T(N5I)) 00008150
C STORE AS CURRENT VALUE OF Y 00008160
T(I)=Y 00008170
C STORE Y IN TEMPORARY STORAGE 00008180
8 T(N9I)=T(I) 00008190
C CALCULATE Y PRIME 00008200
CALL DAUX(T,X,H) 00008210
C STORE DERIVATIVES AS NEEDED FOR ADAMS-MOULTON 00008220
GO TO (10,11,9),KFLAG 00008230
C STORE Y PRIME AT N-2 00008240
10 DO 18 I=1,N 00008250
N7I=I+N7 00008260
NNI=I+NN 00008270
18 T(N7I)=T(NNI) 00008280
GO TO 9 00008290
C STORE Y PRIME AT N-1 00008300
11 DO 19 I=1,N 00008310
N6I=I+N6 00008320
NNI=I+NN 00008330
19 T(N6I)=T(NNI) 00008340
9 CONTINUE 00008350
GO TO 17 00008360
C ADAMS-MOULTON INTEGRATION 00008370
C STORE CURRENT Y AND Y PRIME IN TEMPORARY STORAGE 00008380
12 DO 13 I=1,N 00008390
N9I=I+N9 00008400
N2I=I+N2 00008410
NNI=I+NN 00008420
T(N9I)=T(I) 00008430
13 T(N2I)=T(NNI) 00008440
C CALCULATE PREDICTED VALUE OF Y 00008450
DO 14 I=1,N 00008460
N2I=I+N2 00008470
N6I=I+N6 00008480
N7I=I+N7 00008490
N8I=I+N8 00008500
N9I=I+N9 00008510
YP=T(N9I)+H24*(CMPLX(55.0,0.0)*T(N2I) 00008520
*-CMPLX(59.0,0.0)*T(N6I)+CMPLX(37.0,0.0)* 00008530
*T(N7I)-CMPLX(9.0,0.0)*T(N8I)) 00008540
C STORE AS PREDICTED FUNCTIONAL VALUE OF Y AT X=X+H 00008550
14 T(I)=YP 00008560
C STEP UP X 00008570
X=X+H 00008580
C CALCULATE Y PRIME USING PREDICTED Y 00008590
CALL DAUX(T,X,H) 00008600
C CALCULATE CORRECTED Y 00008610
DO 15 I=1,N 00008620
NNI=I+NN 00008630
N2I=I+N2 00008640
N6I=I+N6 00008650
N7I=I+N7 00008660
N9I=I+N9 00008670
YC=T(N9I)+H24*(CMPLX(9.0,0.0)*T(NNI) 00008680
*+CMPLX(19.0,0.0)*T(N2I)-CMPLX(5.0,0.0)* 00008690
*T(N6I)+T(N7I)) 00008700
C STORE AS NEW CURRENT VALUE OF Y 00008710
15 T(I)=YC 00008720
C CALCULATE Y PRIME TO BE USED IN NEW STEP 00008730
CALL DAUX(T,X,H) 00008740
C REARRANGE STORAGE OF PREVIOUS DERIVATIVES 00008750
DO 16 I=1,N 00008760
N8I=I+N8 00008770
N7I=I+N7 00008780
N6I=I+N6 00008790
N2I=I+N2 00008800
C Y PRIME (N-2) GOES TO (N-3) 00008810
T(N8I)=T(N7I) 00008820
C Y PRIME (N-1) GOES TO (N-2) 00008830
T(N7I)=T(N6I) 00008840
C Y PRIME (N) GOES TO (N-1) 00008850
16 T(N6I)=T(N2I) 00008860
17 RETURN 00008870
END 00008880
C 00008890
C HERE ARE THE CALCULUS SUBROUTINES 00008900
C 00008910
SUBROUTINE VEC(K,V,Z) 00008920
C CALCULUS SUBROUTINE FOR VECTORIZING THE KTH INDEPENDENT VAR V. 00008930
IMPLICIT COMPLEX (A-H,O-Z) 00008940
DIMENSION Z(66) 00008950
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00008960
DO 1 I=1,NDIM 00008970
1 Z(I)=CMPLX(0.0,0.0) 00008980
Z(1)=V 00008990
Z(K+1)=CMPLX(1.0,0.0) 00009000
RETURN 00009010
END 00009020
C 00009030
SUBROUTINE CON(C,Z) 00009040
C CALCULUS SUBROUTINE FOR VECTORIZING A CONSTANT C. 00009050
IMPLICIT COMPLEX (A-H,O-Z) 00009060
DIMENSION Z(66) 00009070
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00009080
DO 1 I=1,NDIM 00009090
1 Z(I)=CMPLX(0.0,0.0) 00009100
Z(1)=C 00009110
RETURN 00009120
END 00009130
C 00009140
SUBROUTINE RECP(X,Z) 00009150
C CALCULUS SUBROUTINE FOR THE RECIPROCAL FUNCTION 1/X. 00009160
IMPLICIT COMPLEX (A-H,O-Z) 00009170
DIMENSION X(66),Z(66),F(5) 00009180
F(1)=CMPLX(1.0,0.0)/X(1) 00009190
F(2)=CMPLX(-1.0,0.0)/(X(1)*X(1)) 00009200
F(3)=CMPLX(2.0,0.0)/(X(1)*X(1)*X(1)) 00009210
F(4)=CMPLX(-6.0,0.0)/(X(1)*X(1)*X(1)*X(1)) 00009220
CALL DER(F,X,Z) 00009230
RETURN 00009240
END 00009250
C 00009260
SUBROUTINE COSS(X,Z) 00009270
C CALCULUS SUBROUTINE FOR THE FUNCTION COS(X). 00009280
IMPLICIT COMPLEX (A-H,O-Z) 00009290
DIMENSION X(66),Z(66),F(5) 00009300
C=X(1) 00009310
F(1)=CCOS(C) 00009320
F(2)=-CSIN(C) 00009330
F(3)=-CCOS(C) 00009340
F(4)=CSIN(C) 00009350
CALL DER(F,X,Z) 00009360
RETURN 00009370
END 00009380
C 00009390
SUBROUTINE SINN(X,Y) 00009400
C CALCULUS SUBROUTINE FOR THE FUNCTION SIN(X). 00009410
IMPLICIT COMPLEX (A-H,O-Z) 00009420
DIMENSION X(66),Y(66),F(5) 00009430
C=X(1) 00009440
F(1)=CSIN(C) 00009450
F(2)=CCOS(C) 00009460
F(3)=-CSIN(C) 00009470
F(4)=-CCOS(C) 00009480
CALL DER(F,X,Y) 00009490
RETURN 00009500
END 00009510
C 00009520
SUBROUTINE POW(X,C,Z) 00009530
C CALCULUS SUBROUTINE FOR X RAISED TO A CONSTANT POWER C. 00009540
IMPLICIT COMPLEX (A-H,O-Z) 00009550
DIMENSION X(66),Z(66),F(5) 00009560
C1=X(1) 00009570
C1=CLOG(C1) 00009580
C1=C1*C 00009590
C2=CEXP(C1) 00009600
F(1)=C2 00009610
F(2)=C*C2/X(1) 00009620
F(3)=(C-CMPLX(1.0,0.0))*F(2)/X(1) 00009630
F(4)=(C-CMPLX(2.0,0.0))*F(3)/X(1) 00009640
CALL DER(F,X,Z) 00009650
RETURN 00009660
END 00009670
C 00009680
SUBROUTINE EXPP(X,Z) 00009690
C CALCULUS SUBROUTINE FOR THE EXPONENTIAL FUNCTION EXP(X). 00009700
IMPLICIT COMPLEX (A-H,O-Z) 00009710
DIMENSION X(66),Z(66),F(5) 00009720
C=X(1) 00009730
F(1)=CEXP(C) 00009740
F(2)=F(1) 00009750
F(3)=F(1) 00009760
F(4)=F(1) 00009770
CALL DER(F,X,Z) 00009780
RETURN 00009790
END 00009800
C 00009810
SUBROUTINE XPY(X,Y,Z) 00009820
C CALCULUS SUBROUTINE FOR X RAISED TO THE POWER Y. 00009830
IMPLICIT COMPLEX (A-H,O-Z) 00009840
DIMENSION X(66),Y(66),Z(66),Z1(66),Z2(66) 00009850
CALL LOGG(X,Z1) 00009860
CALL MULT(Z1,Y,Z2) 00009870
CALL EXPP(Z2,Z) 00009880
RETURN 00009890
END 00009900
C 00009910
SUBROUTINE CPX(C,X,Z) 00009920
C CALCULUS SUBROUTINE FOR A CONSTANT C RAISED TO THE POWER X. 00009930
IMPLICIT COMPLEX (A-H,O-Z) 00009940
DIMENSION X(66),Z(66),F(5) 00009950
C1=X(1) 00009960
C2=CLOG(C) 00009970
C3=C1*C2 00009980
F(1)=CEXP(C3) 00009990
F(2)=F(1)*C2 00010000
F(3)=F(2)*C2 00010010
F(4)=F(3)*C2 00010020
CALL DER(F,X,Z) 00010030
RETURN 00010040
END 00010050
C 00010060
SUBROUTINE LOGG(X,Z) 00010070
C CALCULUS SUBROUTINE FOR THE FUNCTION LOG(X). 00010080
IMPLICIT COMPLEX (A-H,O-Z) 00010090
DIMENSION X(66),Z(66),F(5) 00010100
C=X(1) 00010110
F(1)=CLOG(C) 00010120
F(2)=CMPLX(1.0,0.0)/X(1) 00010130
F(3)=-F(2)*F(2) 00010140
F(4)=CMPLX(2.0,0.0)/(C*C*C) 00010150
CALL DER(F,X,Z) 00010160
RETURN 00010170
END 00010180
C 00010190
SUBROUTINE ADD(X,Y,Z) 00010200
C CALCULUS SUBROUTINE FOR THE FUNCTION Z = X + Y. 00010210
IMPLICIT COMPLEX (A-H,O-Z) 00010220
DIMENSION X(66),Y(66),Z(66) 00010230
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00010240
DO 1 I=1,NDIM 00010250
Z(I)=X(I)+Y(I) 00010260
1 CONTINUE 00010270
RETURN 00010280
END 00010290
C 00010300
SUBROUTINE SUB(X,Y,Z) 00010310
C CALCULUS SUBROUTINE FOR THE FUNCTION Z = X - Y. 00010320
IMPLICIT COMPLEX (A-H,O-Z) 00010330
DIMENSION X(66),Y(66),Z(66) 00010340
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00010350
DO 1 I=1,NDIM 00010360
Z(I)=X(I)-Y(I) 00010370
1 CONTINUE 00010380
RETURN 00010390
END 00010400
C 00010410
SUBROUTINE MULT(X,Y,Z) 00010420
C CALCULUS SUBROUTINE FOR THE FUNCTION Z = X*Y. 00010430
IMPLICIT COMPLEX (A-H,O-Z) 00010440
DIMENSION X(66),Y(66),Z(66) 00010450
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00010460
Z(1)=X(1)*Y(1) 00010470
DO 1 I=1,NVAR 00010480
1 Z(I+1)=X(1)*Y(I+1)+Y(1)*X(I+1) 00010490
KK=NVAR+1 00010500
DO 2 I=1,NVAR 00010510
DO 3 J=I,NVAR 00010520
KK=KK+1 00010530
3 Z(KK)=X(KK)*Y(1)+X(I+1)*Y(J+1)+X(J+1)*Y(I+1)+X(1)*Y(KK) 00010540
2 CONTINUE 00010550
RETURN 00010560
END 00010570
C 00010580
SUBROUTINE DIV(X,Y,Z) 00010590
C CALCULUS SUBROUTINE FOR THE FUNCTION Z = X/Y. 00010600
IMPLICIT COMPLEX (A-H,O-Z) 00010610
DIMENSION X(66),Y(66),Z(66),Z1(66) 00010620
CALL RECP(Y,Z1) 00010630
CALL MULT(X,Z1,Z) 00010640
RETURN 00010650
END 00010660
C 00010670
SUBROUTINE CONVE(C,X,Z) 00010680
C CALCULUS SUBROUTINE FOR A CONSTANT C TIMES A VECTOR X. 00010690
IMPLICIT COMPLEX (A-H,O-Z) 00010700
DIMENSION X(66),Z(66) 00010710
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00010720
DO 1 I=1,NDIM 00010730
Z(I)=C*X(I) 00010740
1 CONTINUE 00010750
RETURN 00010760
END 00010770
C 00010780
SUBROUTINE DER(F,X,Z) 00010790
C CALCULUS SUBROUTINE FOR THE FUNCTION OF A FUNCTION RULE: 00010800
C GIVEN U = F(X), THIS SUBROUTINE OBTAINS DU = F'(X)*DX AND 00010810
C DDU = F''(X)*DX*DX + F'(X)*DDX 00010820
IMPLICIT COMPLEX (A-H,O-Z) 00010830
DIMENSION F(5),X(66),Z(66) 00010840
COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE 00010850
Z(1)=F(1) 00010860
DO 1 I=1,NVAR 00010870
1 Z(I+1)=F(2)*X(I+1) 00010880
KK=NVAR+1 00010890
DO 2 I=1,NVAR 00010900
DO 3 J=I,NVAR 00010910
KK=KK+1 00010920
3 Z(KK)=F(3)*X(I+1)*X(J+1)+F(2)*X(KK) 00010930
2 CONTINUE 00010940
RETURN 00010950
END 00010960
C 00010970
C HERE ARE THE MATRIX SUBROUTINES 00010980
C 00010990
SUBROUTINE VECON(C,V,N,W) 00011000
C CALCULATES THE PRODUCT C*V OF A SCALAR C AND AN N-DIMENSIONAL 00011010
C VECTOR V AND PUTS IT IN W. 00011020
IMPLICIT COMPLEX (A-H,O-Z) 00011030
DIMENSION V(10),W(10) 00011040
DO 10 I=1,N 00011050
W(I)=C*V(I) 00011060
10 CONTINUE 00011070
RETURN 00011080
END 00011090
C 00011100
SUBROUTINE CONMAT(C,B,N,Z) 00011110
C CALCULATES THE PRODUCT C*B OF A CONSTANT C TIMES AN N BY N MATRIX 00011120
C B AND PUTS IT IS Z. 00011130
IMPLICIT COMPLEX (A-H,O-Z) 00011140
DIMENSION B(10,10),Z(10,10) 00011150
DO 10 I=1,N 00011160
DO 20 J=1,N 00011170
Z(I,J)=B(I,J)*C 00011180
20 CONTINUE 00011190
10 CONTINUE 00011200
RETURN 00011210
END 00011220
C 00011230
SUBROUTINE MULMAT(A,B,N,D) 00011240
C MULTIPLIES TWO N BY N MATRICES A AND B AND PUTS THE PRODUCT 00011250
C IN D. 00011260
IMPLICIT COMPLEX (A-H,O-Z) 00011270
DIMENSION A(10,10),B(10,10),D(10,10) 00011280
DO 10 I=1,N 00011290
DO 20 J=1,N 00011300
SUM=CMPLX(0.0,0.0) 00011310
DO 30 K=1,N 00011320
SUM=SUM+A(I,K)*B(K,J) 00011330
30 CONTINUE 00011340
D(I,J)=SUM 00011350
20 CONTINUE 00011360
10 CONTINUE 00011370
RETURN 00011380
END 00011390
C 00011400
SUBROUTINE SUBMAT(A,B,N,D) 00011410
C SUBTRACTS TWO N BY N MATRICES A AND B AND PUTS THEIR DIFFERENCE 00011420
C IN D. 00011430
IMPLICIT COMPLEX (A-H,O-Z) 00011440
DIMENSION B(10,10),D(10,10),A(10,10) 00011450
DO 10 I=1,N 00011460
DO 20 J=1,N 00011470
D(I,J)=A(I,J)-B(I,J) 00011480
20 CONTINUE 00011490
10 CONTINUE 00011500
RETURN 00011510
END 00011520
C 00011530
SUBROUTINE MATVEC(A,W,N,V) 00011540
C CALCULATES THE PRODUCT A*W OF AN N BY N MATRIX A AND AN 00011550
C N-DIMENSIONAL VECTOR W AND PUTS THE PRODUCT IN V. 00011560
IMPLICIT COMPLEX (A-H,O-Z) 00011570
DIMENSION W(10),V(10),A(10,10) 00011580
DO 10 I=1,N 00011590
SUM=CMPLX(0.0,0.0) 00011600
DO 20 J=1,N 00011610
SUM=SUM+A(I,J)*W(J) 00011620
20 CONTINUE 00011630
V(I)=SUM 00011640
10 CONTINUE 00011650
RETURN 00011660
END 00011670
C 00011680
SUBROUTINE TRACE(A,N,TR) 00011690
C CALCULATES THE TRACE OF AN N BY N MATRIX A AND PUTS THE RESULT 00011700
C IN TR 00011710
IMPLICIT COMPLEX (A-H,O-Z) 00011720
DIMENSION A(10,10) 00011730
SUM=CMPLX(0.0,0.0) 00011740
DO 10 I=1,N 00011750
SUM=SUM+A(I,I) 00011760
10 CONTINUE 00011770
TR=SUM 00011780
RETURN 00011790
END 00011800