C SUBROUTINE FILE 'MEANS.FOR' CONTAINING SOME UTILITY PROGRAMS HELPFUL C WHEN DEALING WITH THE MODEL. C C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK C C THIS FILE CONSISTS OF THE FOLLOWING SUBROUTINES: C CDE... SUBROUTINE DESIGNED TO SEARCH FOR THE POINT OF C INTERSECTION OF THE GIVEN CURVE WITH THE BOUNDARIES OF THE C COMPLEX BLOCK. THE CURVE IS INTERPOLATED FROM THE TWO C GIVEN POINTS. THE SUBROUTINE PERFORMS THE STEPS 5.8.3(C), C (D) AND (E) OF OF THE ALGORITHM. C CROSS...SUBROUTINE DESIGNED TO FIND THE POINT OF INTERSECTION OF A C CURVE WITH A SURFACE (SEE C.R.T.5.8.4B). C HIVD2...SUBROUTINE PERFORMING THE HERMITE INTERPOLATION OF A C VECTOR AND ITS DERIVATIVES USING FUNCTIONAL VALUES AND C DERIVATIVES AT 2 GIVEN POINTS (SEE C.R.T.5.8.4A). C SMVPRD..SUBROUTINE DESIGNED TO EVALUATE SYMMETRIC MATRIX BY VECTOR C PRODUCT. IT MAY BE, E.G., CALLED AFTER THE INVOCATION OF C THE METRIC SUBROUTINE TO TRANSFORM THE COVARIANT VECTORS C TO THE CONTRAVARIANT ONES AND VICE VERSA. C C DATE: 1992, NOVEMBER 2 C C======================================================================= C SUBROUTINE CDE(NOUGHT,NEND,KEND,NBOUND,KBOUND,BOUND, * KDIM1,KDIM2,NDIM,IY,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB) INTEGER NOUGHT,NEND,KEND(*),NBOUND,KBOUND(*) INTEGER IY(8),KDIM1,KDIM2,NDIM,MY PARAMETER (MY=35) REAL BOUND(*),ERR,X0,X1,Y1(NDIM),D1(NDIM),X2,Y2(NDIM),D2(NDIM) REAL X,Y(NDIM),D(NDIM),XB,YB(NDIM),DB(NDIM) C C THIS SUBROUTINE DETERMINES THE POINT OF INTERSECTION OF THE GIVEN C CURVE ELEMENT WITH THE BOUNDARY OF THE COMPLEX BLOCK. THE CURVE IS C INTERPOLATED FROM THE TWO GIVEN POINTS. THE SUBROUTINE PERFORMS THE C STEPS 5.8.3(C), (D) AND (E) OF OF THE ALGORITHM. C C GENERAL MEANING OF SOME ARGUMENTS USED: C X1,X2,X,XB... INDEPENDENT VARIABLE ALONG THE CURVE. C Y1,Y2,Y,YB,D1,D2,D,DB... ARRAYS OF THE DIMENSION NDIM. C Y1(KDIM1:KDIM2),Y2(KDIM1:KDIM2),Y(KDIM1:KDIM2),YB(KDIM1:KDIM2)... C COORDINATES. C D1(1:NDIM),D2(1:NDIM),D(1:NDIM),DB(1:NDIM)... DERIVATIVES OF C ARRAYS Y1,Y2,Y,YB WITH RESPECT TO THE INDEPENDENT C VARIABLE. C C INPUT: C NOUGHT..ZERO IF THE CURVE IS NOT SITUATED ALONG THE STRUCTURAL C INTERFACE, C OTHERWISE THE INDEX OF THE SURFACE ON WHICH THE CURVE IS C SITUATED. C NEND... NUMBER OF END SURFACES LIMITING THE COMPUTATIONAL VOLUME. C USUALLY NEND=0. C KEND... CONTAINS THE INDICES OF END SURFACES. C ARRAY OF DIMENSION NEND, NOT USED IF NEND=0. C NBOUND..NUMBER OF ISOSURFACES OF COMPUTED QUANTITIES, LIMITING THE C COMPUTATIONAL VOLUME. THESE ISOSURFACES WILL LIKELY BE THE C BOUNDARIES X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX OF THE C COMPUTATIONAL VOLUME. C KBOUND..INDICES OF THE QUANTITIES CORRESPONDING TO THE ISOSURFACES. C THE COMPUTATIONAL VOLUME IS LIMITED BY THE INEQUALITIES C Y(IABS(KBOUND(I))).LE.BOUND(I) FOR KBOUND(I).LT.0, C Y(IABS(KBOUND(I))).GE.BOUND(I) FOR KBOUND(I).GT.0. C ARRAY OF DIMENSION NBOUND, NOT USED IF NBOUND=0. C BOUND...VALUES OF THE ISOSURFACES, LIKELY COINCIDING WITH C BOUNDARIES X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX OF THE C COMPUTATIONAL VOLUME. IN SUCH A CASE, KBOUND WOULD TAKE C THE VALUES KDIM1,-KDIM1,KDIM1+1,-KDIM1-1,KDIM2,-KDIM2. C ARRAY OF DIMENSION NBOUND, NOT USED IF NBOUND=0. C KDIM1,KDIM2... KDIM1-TH TO KDIM2-TH ELEMENTS OF ARRAYS Y1,Y2,Y,YB C CONTAIN COORDINATES. KDIM2 IS ASSUMED TO BE KDIM1+2. C NDIM... DIMENSION OF THE ARRAYS Y1,Y2,Y,YB,D1,D2,D,DB. IT MUST C NOT EXCEED THE PARAMETER MY. C IY... INTEGER ARRAY OF THE DIMENSION AT LEAST 8. C IY(4)=ISB1... INDEX OF THE SIMPLE BLOCK CONTAINING THE POINT X1. C IY(5)=ICB1... INDEX OF THE COMPLEX BLOCK CONTAINING THE POINT X1. C IT MAY (BUT NEED NOT) BE SUPPLEMENTED BY A SIGN '+' FOR P C WAVE AND SIGN '-' FOR S WAVE. C ERR... MAXIMUM ERROR IN INDEPENDENT VARIABLE FOR THE C DETERMINATION OF THE POINT OF INTERSECTION. C X0... FOR X.LE.X0 JUST THE INTERSECTION WITH THE BOUNDARY OF C THE COMPUTATIONAL VOLUME IS CHECKED, NOT WITH THE C STRUCTURAL INTERFACES. C FOR X1.LE.X0, THE INITIAL POINT X1 IS CHECKED FOR LOCATION C OUTSIDE THE COMPLEX BLOCK, BUT VERY CLOSE TO ITS BOUNDARY. C IN SUCH A CASE, X0 SHOULD BE CLOSE TO X1 (WITHIN AN C INTERVAL OF ERR) AND MAY BE LOCATED INSIDE THE COMPLEX C BLOCK. C X1,Y1,D1... QUANTITIES AT THE INITIAL POINT OF THE CURVE ELEMENT. C X2,Y2,D2... QUANTITIES AT ANOTHER POINT OF THE CURVE. THE CURVE C IS INTERPOLATED USING THE VALUES AND THEIR DERIVATIVES AT C THE POINTS X1 AND X2. C X,Y,D...QUANTITIES AT THE ENDPOINT OF THE CURVE ELEMENT. C XB,YB,DB... VALUES IGNORED. C C OUTPUT: C NOUGHT,KDIM1,KDIM2,NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2... UNCHANGED C INPUT. C IY(1),IY(2),IY(3),IY(5)... UNCHANGED INPUT. C IY(4)=ISB1... INDEX OF THE SIMPLE BLOCK CONTAINING THE POINT X. C OUTPUT IF THE ENDPOINT OF THE CURVE ELEMENT IS SITUATED IN THE COMPLEX C BLOCK IY(5): C IY(6),IY(7),IY(8)... UNCHANGED INPUT. C X,Y,D...UNCHANGED INPUT. C XB,YB,DB... UNDEFINED. C OUTPUT IF THE ENDPOINT OF THE CURVE ELEMENT IS SITUATED IN ANOTHER C COMPLEX BLOCK OR OUTSIDE THE COMPUTATIONAL VOLUME: C IY(6)=ISRF... INDEX OF THE SURFACE AT WHICH THE POINT OF C INTERSECTION WITH THE BOUNDARY OF THE COMPLEX BLOCK IS C SITUATED, SUPPLEMENTED BY A SIGN '+' OR '-' FOR THE POINT C SITUATED AT THE POSITIVE OR NEGATIVE SIDE OF THE SURFACE, C RESPECTIVELY. C IY(7)=ISB2... INDEX OF THE SIMPLE BLOCK TOUCHING THE COMPLEX C BLOCK ICB1 FROM THE OTHER SIDE OF THE SURFACE ISRF AT C THE POINT OF INTERSECTION. C ISB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF ISRF. C IY(8)=ICB2... INDEX OF THE COMPLEX BLOCK TOUCHING THE COMPLEX C BLOCK ICB1 FROM THE OTHER SIDE OF THE SURFACE ISRF AT C THE POINT OF INTERSECTION. C ICB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF ISRF. C X,Y,D...VALUES CORRESPONDING TO THE POINT OF INTERSECTION OF THE C CURVE WITH THE BOUNDARY OF THE COMPLEX BLOCK OR THE C COMPUTATIONAL VOLUME. IF POSSIBLE, THIS POINT IS SITUATED C INSIDE THE GIVEN COMPLEX BLOCK CLOSE TO ITS BOUNDARY, OR C DIRECTLY ON ITS BOUNDARY. C XB,YB,DB... VALUES CORRESPONDING TO ANOTHER APPROXIMATION OF THE C POINT OF INTERSECTION. THIS POINT IS SITUATED OUTSIDE THE C THE GIVEN COMPLEX BLOCK, CLOSE TO ITS BOUNDARY OR DIRECTLY C ON ITS BOUNDARY. C NOTE THAT XB-X SHOULD BE LESS THAN ERR. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL NSRFC,BLOCK,SRFC2,CROSS INTEGER NSRFC C NSRFC,BLOCK... FILE 'MODEL.FOR'. C SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR'). C CROSS,HIVD2... THIS FILE. C C ERROR MESSAGES: C 581... TOO MANY FICTITIOUS INTERFACES: C MORE THAN 100 FICTITIOUS INTERFACES CROSSED DURING ONE C STEP OF THE NUMERICAL INTEGRATION. C THIS ERROR SHOULD NOT APPEAR. CONTACT THE AUTHORS. C 582,583... EXCLUDED PROGRAM BRANCH: C THIS ERROR SHOULD NOT APPEAR. CONTACT THE AUTHORS. C C DATE: 1992, DECEMBER 31 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS: REAL G(12),GAMMA(18),GSQRD,FAUX(10) REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL C THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A C COMMON BLOCK. THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY C COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY. COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL C C G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'. C FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS C DERIVATIVES. C UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'. C UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE C 'MODEL.FOR'. C C....................................................................... C C OTHER AUXILIARY STORAGE LOCATIONS: INTEGER ISRAUX,ISBAUX,ISB1,ISRF2,ISB2,ICB2,K0,K1,K2,K,I REAL XAUX,YAUX(MY),DAUX(MY),X01,Y01(MY),D01(MY) C C ISB1,X01,Y01,D01,VD... USED ONLY TO CHECK FOR THE LOCATION OF THE C POINT X1. C C....................................................................... C C (C) CHECK FOR CROSSING THE COORDINATE BOUNDARIES OF THE C COMPUTATIONAL VOLUME: DO 31 I=1,NBOUND K=IABS(KBOUND(I)) IF((KBOUND(I).GT.0.AND.Y(K).LT.BOUND(I)).OR. * (KBOUND(I).LT.0.AND.Y(K).GT.BOUND(I))) THEN IY(6)=100+I FAUX(1)=BOUND(I) CALL CROSS(SRFC2,IY(6),K,K,NDIM,ERR,X1,Y1,D1,X2,Y2,D2, * X,Y,D,XB,YB,DB,FAUX) END IF 31 CONTINUE IF(IY(6).GT.100) THEN IY(7)=IY(4) IY(8)=IABS(IY(5)) END IF ISB1=IY(4) C C (D) CHECK FOR CROSSING THE BOUNDARY OF THE COMPLEX BLOCK: C NOTE: IY(4)=ISB1, IY(5)=ICB1, IY(6)=ISRF. IF(X.GT.X0) THEN CALL BLOCK(Y(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2,FAUX) 40 CONTINUE IF(ISRF2.NE.0) THEN C BOUNDARY OF THE SIMPLE BLOCK IS CROSSED C NOTE: IN THIS ROUTINE, UNLIKE IN THE PAPER ON C.R.T., THE C POINT OF INTERSECTION WITH THE BOUNDARY OF THE SIMPLE BLOCK IS C FOUND EVEN IF THE BOUNDARY OF THE COMPLEX BLOCK IS NOT C CROSSED. C (D1) ISRAUX=IY(6) ISBAUX=ISB2 XAUX=X DO 41 I=1,NDIM YAUX(I)=Y(I) DAUX(I)=D(I) 41 CONTINUE C FOLLOWING LOOP IS INCLUDED TO AVOID INFINITE REPEATING OF THE C STEPS (D2) AND (D3) OF THE ALGORITHM DO 47 K=1,100 C (D2) IY(6)=ISRF2 C (D3) C CHECK FOR THE LOCATION OF THE POINT X1 BEFORE CALLING CROSS IF(X1.LE.X0) THEN C X1 MAY BE LOCATED OUTSIDE THE COMPLEX BLOCK CALL SRFC2(IABS(IY(6)),Y1(KDIM1),VD) IF(VD(1)*FAUX(1).GT.0.) THEN C POINTS X1 AND X ARE LOCATED OUTSIDE THE COMPLEX BLOCK, C LOOKING FOR THE POINT X01 BETWEEN X0 AND X, SITUATED C INSIDE THE COMPLEX BLOCK. X01=X0 42 CONTINUE CALL HIVD2(KDIM2-KDIM1+1,X1,Y1(KDIM1),D1(KDIM1),X2, * Y2(KDIM1),D2(KDIM1),X01,Y01(KDIM1),D01(KDIM1)) CALL SRFC2(IABS(IY(6)),Y01(KDIM1),VD) IF(VD(1)*FAUX(1).LE.0.) THEN C POINT X01 IS LIKELY LOCATED INSIDE THE COMPLEX BLOCK CALL BLOCK(Y01(KDIM1),NOUGHT,ISB1,ISRF2,K1,K2,VD) IF(ISRF2.EQ.0) THEN C POINT X01 IS LOCATED INSIDE THE SIMPLE BLOCK ISB1, C POINT X1 MAY BE REPLACED BY X01. CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,X01,Y01,D01) CALL CROSS(SRFC2,IABS(IY(6)),KDIM1,KDIM2,NDIM,ERR, * X01,Y01,D01,X2,Y2,D2,X,Y,D,XB,YB,DB,FAUX) GO TO 43 END IF ELSE C TRYING A NEW POINT X01 IF(X01.EQ.X0) THEN X01=X01+ERR ELSE X01=X01+(X01-X0) END IF IF(X01.LT.X) THEN GO TO 42 END IF END IF END IF END IF CALL CROSS(SRFC2,IABS(IY(6)),KDIM1,KDIM2,NDIM,ERR, * X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB,FAUX) 43 CONTINUE C X AND XB ARE THE APPROXIMATIONS OF THE POINT OF INTERSECTION C WITH THE SURFACE IY(6) CALL BLOCK(Y(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2,FAUX) IF(ISRF2.NE.0.AND.(ISRF2.NE.IY(6).OR.X.NE.X1).AND. * FAUX(1)*(D(KDIM1)*FAUX(2)+D(KDIM1+1)*FAUX(3) * +D(KDIM2)*FAUX(4)).GT.0.) THEN C (D3-I) C POINT X IS NOT SITUATED AT THE BOUNDARY OF THE SIMPLE C BLOCK. IT IS SEPARATED FROM THE SIMPLE BLOCK IY(4) BY C THE SURFACE ISRF2 SITUATED BEFORE (NOT AFTER) THE POINT X. C GO TO (D2) ELSE C X IS SITUATED AT THE BOUNDARY OF THE SIMPLE BLOCK IY(4) CALL BLOCK(YB(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2,FAUX) IF(ISB2.EQ.IY(4)) THEN CALL SRFC2(IABS(IY(6)),YB(KDIM1),FAUX) IF(FAUX(1).NE.0.) THEN PAUSE 'ERROR 582 IN CDE: EXCLUDED PROGRAM BRANCH' STOP END IF C POINT XB IS SITUATED EXACTLY AT THE SURFACE IY(6) CALL BLOCK(YB(KDIM1),IY(6),IY(4),ISRF2,ISB2,ICB2,FAUX) IF(NOUGHT.NE.0.AND.NOUGHT.EQ.ISRF2) THEN K1=ISB2 CALL BLOCK(YB(KDIM1),IY(6),K1,ISRF2,ISB2,ICB2,FAUX) END IF END IF IF(NOUGHT.EQ.0) THEN C NEAR THE EDGE OF A SIMPLE BLOCK, TWO DIFFERENT SURFACES C BOUNDING SIMPLE BLOCK IY(4) MAY BE SITUATED BETWEEN C POINTS X AND XB. IN SUCH A CASE, ONLY ONE OF THEM C SEPARATES SIMPLE BLOCK IY(4) FROM ISB2 AND IY(6) MAY C BE THE SECOND. CALL BLOCK(Y(KDIM1),IY(6),IY(4),ISRF2,K1,K2,FAUX) IF(ISRF2.EQ.0.AND.ISB2.NE.K1) THEN C SURFACE IY(6) DOES NOT FORM THE INTERFACE BETWEEN C SIMPLE BLOCKS IY(4) AND ISB2: CALL BLOCK(YB(KDIM1),IY(6),IY(4),ISRF2,K1,K2,FAUX) IF(ISRF2.NE.0) THEN C SURFACE ISRF2 SEPARATES THE POINT XB FROM THE SIMPLE C BLOCK IY(4): CALL BLOCK(Y(KDIM1),ISRF2,IY(4),K0,K1,K2,FAUX) IF(K0.EQ.0.AND.ISB2.EQ.K1) THEN C SURFACE ISRF2 FORMS THE INTERFACE BETWEEN SIMPLE C BLOCKS IY(4) AND ISB2: IY(6)=ISRF2 END IF END IF END IF END IF IF(ICB2.EQ.IABS(IY(5))) THEN C (D3-II) C X,Y IS SITUATED AT THE BOUNDARY OF THE SIMPLE BLOCK C BUT NOT SITUATED AT THE BOUNDARY OF THE COMPLEX BLOCK IY(4)=ISB2 X=XAUX DO 46 I=1,NDIM Y(I)=YAUX(I) D(I)=DAUX(I) 46 CONTINUE IF(ISB2.EQ.ISBAUX) THEN C BOUNDARY OF THE COMPLEX BLOCK HAS NOT BEEN CROSSED C DURING THE LAST STEP OF NUMERICAL INTEGRATION IY(6)=ISRAUX GO TO 49 END IF CALL BLOCK(YAUX(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2 * ,FAUX) IF(ISRF2.EQ.0) THEN C ISRF2 CAN BE ZERO ONLY IF ISB2.EQ.ISBAUX: PAUSE 'ERROR 583 IN CDE: EXCLUDED PROGRAM BRANCH' STOP END IF C GO TO (D2) ELSE C (D3-III) C X IS SITUATED AT THE BOUNDARY OF THE SIMPLE BLOCK AND C X IS SITUATED AT THE BOUNDARY OF THE COMPLEX BLOCK GO TO 48 END IF END IF 47 CONTINUE PAUSE 'ERROR 581 IN CDE: TOO MANY FICTITIOUS INTERFACES' STOP 48 CONTINUE IY(7)=ISB2 IY(8)=ICB2 END IF 49 CONTINUE END IF C C (E) CHECK FOR CROSSING THE END SURFACES DO 51 I=1,NEND IF(IABS(KEND(I)).GT.NSRFC()) THEN CALL SRFC2(IABS(KEND(I)),Y(KDIM1),FAUX) IF(FAUX(1)*FLOAT(KEND(I)).LE.0.) THEN IY(6)=IABS(KEND(I)) CALL CROSS(SRFC2,IY(6),KDIM1,KDIM2,NDIM,ERR, * X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB,FAUX) END IF END IF 51 CONTINUE C RETURN END C C======================================================================= C SUBROUTINE CROSS(SRFC2,ISRF,KDIM1,KDIM2,NDIM, * ERR,X1,Y1,D1,X2,Y2,D2,XA,YA,DA,XB,YB,DB,F) EXTERNAL SRFC2 INTEGER ISRF,KDIM1,KDIM2,NDIM REAL ERR,X1,Y1(NDIM),D1(NDIM),X2,Y2(NDIM),D2(NDIM) REAL XA,YA(NDIM),DA(NDIM),XB,YB(NDIM),DB(NDIM),F(10) C C THIS SUBROUTINE FINDS THE POINT OF INTERSECTION OF A CURVE WITH C A SURFACE (SEE C.R.T.5.8.4B). THE CURVE IS PARAMETRIZED BY AN C INDEPENDENT VARIABLE X AND EVALUATED BY THE HERMITE INTERPOLATION FROM C THE TWO GIVEN POINTS. THE SURFACE IS SPECIFIED IN AN IMPLICIT WAY BY C SUBROUTINE SRFC2 WHICH IS DESCRIBED ELSWHERE, OR MAY COINCIDE WITH AN C ISOSURFACE OF A COMPUTED QUANTITY (E.G. WITH A COORDINATE PLANE). C C INPUT: C SRFC2...NAME OF THE EXTERNAL PROCEDURE EVALUATING THE FUNCTION C DESCRIBING THE SURFACE ISRF, FOR ISRF.LE.100. C ISRF... INDEX OF THE SURFACE. C FOR ISRF.LE.100: C THE SURFACE COINCIDES WITH THE ZERO ISOSURFACE OF THE C FUNCTION NO. ISRF EVALUATED BY THE SUBROUTINE SRFC2. C FOR ISRF.GT.100: C THE SURFACE COINCIDES WITH AN ISOSURFACE Y(KDIM1)=F(1). C KDIM1,KDIM2... INDICES OF QUANTITIES ON WHICH THE FUNCTION C DESCRIBING THE SURFACE IS DEPENDENT: C FOR ISRF.LE.100: C KDIM1-TH TO KDIM2-TH ELEMENTS OF ARRAYS Y1, Y2, YA, AND C YB CONTAIN COORDINATES. KDIM2 IS ASSUMED TO BE KDIM1+1 C (2 COORDINATES) OR KDIM1+2 (3 COORDINATES). C THE FUNCTION DESCRIBING THE SURFACE ISRF IS DETERMINED C BY THE SUBROUTINE SRFC2. C FOR ISRF.GT.100: C THE SURFACE COINCIDES WITH AN ISOSURFACE Y(KDIM1)=F(1). C KDIM2 IS ASSUMED TO EQUAL KDIM1. C WHEN SEARCHING FOR THE POINT OF INTERSECTION, ONLY C QUANTITIES Y(KDIM1:KDIM2) AND THEIR DERIVATIVES ARE C INTERPOLATED ALONG THE CURVE. C NDIM... DIMENSION OF ARRAYS Y1,D1,Y2,D2,YA,DA,YB,DB. C AT THE POINT OF INTERSECTION, QUANTITIES Y(1:NDIM) AND C THEIR DERIVATIVES ARE INTERPOLATED. C ERR... MAXIMUM ERROR IN INDEPENDENT VARIABLE FOR THE C DETERMINATION OF THE POINT OF INTERSECTION. C X1... INDEPENDENT VARIABLE CORRESPONDING TO THE FIRST POINT C GIVEN FOR THE INTERPOLATION OF THE CURVE. C Y1... ARRAY CONTAINING DEPENDENT VARIABLES AT POINT X1. C Y1(KDIM1) TO Y1(KDIM2) MUST CONTAIN THE COORDINATES OF C POINT X1. C D1... ARRAY CONTAINING THE DERIVATIVES OF THE DEPENDENT C VARIABLES AT POINT X1. C X2... INDEPENDENT VARIABLE CORRESPONDING TO THE SECOND POINT C GIVEN FOR THE INTERPOLATION OF THE CURVE. C Y2... ARRAY CONTAINING DEPENDENT VARIABLES AT POINT X2. C Y2(KDIM1) TO Y2(KDIM2) MUST CONTAIN THE COORDINATES OF C POINT X2. C D2... ARRAY CONTAINING THE DERIVATIVES OF THE DEPENDENT C VARIABLES AT POINT X2. C XA... INDEPENDENT VARIABLE CORRESPONDING TO THE POINT OF THE C CURVE AT WHICH THE FUNCTION SPECIFYING THE SURFACE HAS THE C OPPOSITE SIGN THAN AT X1. THEN THE POINT OF INTERSECTION C IS BEING FOUND BETWEEN THE POINTS X1 AND XA. THE FOUND C APPROXIMATIONS XA AND XB OF THE POINT OF INTERSECTION ARE C SITUATED CLOSE TO THE SURFACE, XA AT THE SAME SIDE AS THE C GIVEN POINT X1, XB AT THE OPOSITE SIDE THAN THE GIVEN C POINT X1. IF, ACCIDENTALY, THE FUNCTION HAS THE SAME SIGN C AT THE POINTS X1 AND XA, THE VALUE AT X1 IS ASSUMED TO C EQUAL ZERO. THEN XA=X1 AND XB=X1 ON OUTPUT. C YA... ARRAY OF THE DIMENSION AT LEAST NDIM. YA(KDIM1) TO C YA(KDIM2) MUST CONTAIN THE COORDINATES OF THE POINT. C OTHER STORAGE LOCATIONS MAY BE UNDEFINED. C DA... ARRAY OF THE DIMENSION AT LEAST NDIM, CONTAINING THE C DERIVATIVES OF THE DEPENDENT VARIABLES AT POINT X. C DA(KDIM1) TO DA(KDIM2) MUST CONTAIN THE DERIVATIVES OF C COORDINATES WITH RESPECT TO X, AT THE POINT XA. OTHER C STORAGE LOCATIONS MAY BE UNDEFINED. C YB... ARRAY OF THE DIMENSION AT LEAST NDIM. C DB... ARRAY OF THE DIMENSION AT LEAST NDIM. C F... FOR ISRF.LE.100: C ARRAY CONTAINING THE VALUE, AND AT LEAST FIRST C DERIVATIVES OF THE FUNCTION ISRF SPECIFYING THE SURFACE, C AT POINT XA. C FOR ISRF.GT.100: C THE SURFACE COINCIDES WITH AN ISOSURFACE Y(KDIM1)=F(1). C ISRF, KDIM1, KDIM2, NDIM, ERR, X1, Y1, D1, X2, Y2, D2 ARE UNALTERED. C C OUTPUT: C XA... INDEPENDENT VARIABLE CORRESPONDING TO THE APPROXIMATION OF C THE POINT OF INTERSECTION, SITUATED CLOSE TO THE SURFACE C AT THE SAME SIDE AS THE GIVEN POINT X1. XA=X1 IF THE C FUNCTION HAD AT POINT X1 THE SAME SIGN AS F(1) ON INPUT. C YA... ARRAY CONTAINING DEPENDENT VARIABLES AT THE POINT XA. C DA... ARRAY CONTAINING THE DERIVATIVES OF THE DEPENDENT C VARIABLES AT THE POINT XA. C XB... INDEPENDENT VARIABLE CORRESPONDING TO THE APPROXIMATION OF C THE POINT OF INTERSECTION, SITUATED CLOSE TO THE SURFACE C AT THE OPOSITE SIDE THAN THE GIVEN POINT X1. XB=X1 IF THE C FUNCTION HAD AT POINT X1 THE SAME SIGN AS F(1) ON INPUT. C YB... ARRAY CONTAINING DEPENDENT VARIABLES AT THE POINT XB. C DB... ARRAY CONTAINING THE DERIVATIVES OF THE DEPENDENT C VARIABLES AT THE POINT XB. C F... UNDEFINED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL HIVD2 C SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') - DUMMY ARGUMENT. C HIVD2... THIS FILE. C C ERROR MESSAGES: C 584... TOO MANY ITERATIONS: C MORE THAN 20 ITERATIONS WHEN DETERMINING THE POINT OF C INTERSECTION OF THE RAY WITH A SURFACE. THIS MAY BE C CAUSED BY TOO SMALL UPPER ERROR BOUND UEB IN THE INPUT C DATA 'DCRT.DAT' - LESS THAN ROUNDING ERRORS. USUALLY, C THIS ERROR SHOULD NOT APPEAR. CONTACT THE AUTHORS. C 585... ZERO-LENGTH INTERVAL: C TWO POINTS X1 AND X2 USED FOR INTERPOLATION OF THE GIVEN C LINE (E.G. THE VELOCITY ISOLINE, OR THE RAY), WHEN C SEARCHING FOR ITS INTERSECTION WITH THE INTERFACE, ARE C IDENTICAL AND NO INTERPOLATION IS POSSIBLE. LIKELY A BUG C IN THE PROCEDURE CALLING THIS SUBROUTINE. C C DATE: 1994, JANUARY 9 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER I,ITER REAL X,FX,DFX,FB,DFB,FA,DFA,XC,XCB,XCA C C....................................................................... C C INITIAL VALUES: IF(ISRF.LE.100) THEN FA=F(1) DFA=F(2)*DA(KDIM1) DO 1 I=1,KDIM2-KDIM1 DFA=DFA+F(2+I)*DA(KDIM1+I) 1 CONTINUE ELSE FA=YA(KDIM1)-F(1) DFA=DA(KDIM1) END IF XB=X1 X =X1 DO 2 I=KDIM1,KDIM2 YB(I)=Y1(I) DB(I)=D1(I) 2 CONTINUE C C CHECK FOR ZERO INTERVALS: IF(XA.EQ.X1) THEN DO 3 I=1,NDIM YA(I)=Y1(I) DA(I)=D1(I) YB(I)=Y1(I) DB(I)=D1(I) 3 CONTINUE RETURN END IF IF(X2.EQ.X1) THEN PAUSE 'ERROR 585 IN CROSS: ZERO-LENGTH INTERVAL' STOP END IF C C ITERATIONS: DO 9 ITER=1,20 C C FUNCTIONAL VALUE AND DERIVATIVE AT X: IF(ISRF.LE.100) THEN CALL SRFC2(ISRF,YB(KDIM1),F) FX=F(1) DFX=F(2)*DB(KDIM1) DO 5 I=1,KDIM2-KDIM1 DFX=DFX+F(2+I)*DB(KDIM1+I) 5 CONTINUE ELSE FX=YB(KDIM1)-F(1) DFX=DB(KDIM1) END IF C C SELECTION OF POINTS: IF(FX.EQ.0.) THEN XA=X ELSE IF(FA*FX.GE.0.) THEN IF(ITER.EQ.1) THEN IF(FA.EQ.0.) THEN X=XA ELSE FX=0. END IF ELSE XA=XB FA=FB DFA=DFB END IF END IF XB=X FB=FX DFB=DFX C C NEW POINT OR END OF ITERATIONS: IF(ABS(XB-XA).LE.ERR) THEN C POINT OF INTERSECTION IS FOUND WITHIN THE SPECIFIED ERROR ERR C *** END OF ITERATIONS *** IF(ABS(XB-X1).LT.ABS(XA-X1)) THEN C POINT XA IS SITUATED AT THE OTHER SIDE OF THE SURFACE THAN C THE POINT X1 XB=XA XA=X END IF CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,XA,YA,DA) CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,XB,YB,DB) RETURN END IF C C NEW APPROXIMATION: IF(MOD(ITER,2).EQ.1) THEN C REGULA FALSI: X=(FA*XB-FB*XA)/(FA-FB) ELSE C MODIFIED NEWTON-RAPHSON: XC=(XA+XB)/2. XCA=XA-FA/DFA+SIGN(ERR/50.,XA-XB) XCB=XB-FB/DFB+SIGN(ERR/50.,XA-XB) IF((XCA-XC)*(XCA-XB).LT.0.) THEN IF((XCB-XC)*(XCB-XB).LT.0.) THEN IF(ABS(XCA-XB).LT.ABS(XCB-XB)) THEN X=XCA ELSE X=XCB END IF ELSE X=XCA END IF ELSE IF((XCB-XC)*(XCB-XB).LT.0.) THEN X=XCB ELSE X=XC IF(ABS(XCB-XB).LT.SQRT(ABS(XC-XB)*ERR)) THEN C ATTEMPT TO HALVE THE NUMBER OF ITERATIONS X=XB+SIGN(SQRT(ABS(XC-XB)*ERR),XA-XB) END IF END IF END IF END IF C C INTERPOLATION OF THE RAY: CALL HIVD2(KDIM2-KDIM1+1,X1,Y1(KDIM1),D1(KDIM1), * X2,Y2(KDIM1),D2(KDIM1),X,YB(KDIM1),DB(KDIM1)) C 9 CONTINUE C END OF LOOP FOR ITERATIONS PAUSE 'ERROR 584 IN CROSS: TOO MANY ITERATIONS' STOP END C C======================================================================= C SUBROUTINE HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,X,Y,D) INTEGER NDIM REAL X1,Y1(NDIM),D1(NDIM),X2,Y2(NDIM),D2(NDIM) REAL X,Y(NDIM),D(NDIM) C C THIS SUBROUTINE PERFORMS HERMITE INTERPOLATION OF A VECTOR AND ITS C DERIVATIVES USING FUNCTIONAL VALUES AND DERIVATIVES AT 2 GIVEN POINTS. C C INPUT: C NDIM... DIMENSION OF ARRAYS Y1,D1,Y2,D2,Y,D (THE NUMBER OF C INDEPENDENT VARIABLES). C X1... INDEPENDENT VARIABLE CORRESPONDING TO THE FIRST GIVEN C POINT. C Y1... ARRAY CONTAINING FUNCTIONAL VALUES AT THE FIRST GIVEN C POINT. C D1... ARRAY CONTAINING THE DERIVATIVES AT THE FIRST GIVEN POINT. C X2... INDEPENDENT VARIABLE CORRESPONDING TO THE SECOND GIVEN C POINT. C Y2... ARRAY CONTAINING FUNCTIONAL VALUES AT THE SECOND GIVEN C POINT. C D2... ARRAY CONTAINING THE DERIVATIVES AT THE SECOND GIVEN POINT C X... INDEPENDENT VARIABLE OF THE POINT AT WHICH THE C INTERPOLATED VECTOR IS TO BE EVALUATED. C NONE OF THE INPUT PARAMETERS ARE ALTERED. C C OUTPUT: C Y... ARRAY CONTAINING INTERPOLATED FUNCTIONAL VALUES AT X. C D... ARRAY CONTAINING THE DERIVATIVES OF THE INTERPOLATED C FUNCTIONAL VALUES AT X. C C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED. C C DATE: 1989, OCTOBER 20 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER I REAL A,B,A1,A2,B1,B2,DA1,DB1,DB2 C C SUBSTITUTIONS: A=(X-X2)/(X1-X2) B=(A-1.)*A C C BASIC FUNCTIONS: A1=(A-B-B)*A A2=1.-A1 B1=B*(X-X2) B2=B*(X-X1) C DERIVATIVES OF BASIC FUNCTIONS: DA1=6.*B/(X2-X1) DB1=3.*B+A DB2=3.*B+1.-A C C INTERPOLATION: DO 1 I=1,NDIM Y(I)=A1*Y1(I)+A2*Y2(I)+B1*D1(I)+B2*D2(I) D(I)=DA1*(Y1(I)-Y2(I))+DB1*D1(I)+DB2*D2(I) 1 CONTINUE C RETURN END C C======================================================================= C SUBROUTINE SMVPRD(G,A1,A2,A3,B1,B2,B3) REAL G(6),A1,A2,A3,B1,B2,B3 C C THIS SUBROUTINE IS AN AUXILIARY ROUTINE TO FCT. IT EVALUATES SYMMETRIC C MATRIX BY VECTOR PRODUCT. C C INPUT: C G... ARRAY CONTAINING COMPONENTS G11, G12, G22, G13, G23, G33 C OF THE 3*3 SYMMETRIC MATRIX. C A1,A2,A3... COMPONENTS OF THE 3-VECTOR. C C OUTPUT: C B1,B2,B3... COMPONENTS OF VECTOR B=G*A. C C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED. C C DATE: 1989, SEPTEMBER 5 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C B1=G(1)*A1+G(2)*A2+G(4)*A3 B2=G(2)*A1+G(3)*A2+G(5)*A3 B3=G(4)*A1+G(5)*A2+G(6)*A3 RETURN END C C======================================================================= C