For additional articles and introductory materials related to NASA, visit the NASA Home Page.
// 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