C
C Subroutine file 'means.for' containing some utility programs helpful C when dealing with the model. C C Date: 2006, June 8 C Coded by Ludek Klimes C C....................................................................... 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 the algorithm. C CDE C CROSS...Subroutine designed to find the point of intersection of a C curve with a surface, see C C.R.T.5.8.4b. C CROSS 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 C.R.T.5.8.4a. C HIVD2 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 SMVPRD C C======================================================================= 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 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 C the boundaries X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX of the C computational volume. C KBOUND..Indices of the quantities corresponding to the C isosurfaces. The computational volume is limited by the C 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... Copy of X,Y,D. 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 Date: 1999, May 17 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C auxmod.inc C C....................................................................... C C Other auxiliary storage locations: INTEGER NSRF1,ISRF(2) 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 XB=X DO 39 I=1,NDIM YB(I)=Y(I) DB(I)=D(I) 39 CONTINUE C C (d) Check for crossing the boundary of the complex block: C Note: IY(4)=ISB1, IY(5)=ICB1, IY(6)=ISRF. IF(NOUGHT.EQ.0) THEN NSRF1=1 ELSE NSRF1=2 ISRF(2)=NOUGHT END IF IF(X.GT.X0) THEN ISRF(1)=0 CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2) 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 5.30: CALL SRFC2(IABS(ISRF2),Y(KDIM1),FAUX) 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 ISRF(1)=0 CALL BLOCKS(Y01(KDIM1),NSRF1,ISRF,ISB1,ISRF2,K1,K2) 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) ISRF(1)=0 CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2) IF(ISRF2.EQ.IY(6).AND.X.NE.X1) THEN C 587 CALL ERROR('587 in CDE: Boundary point out of block') C This error should not appear. Contact the authors. END IF C 5.30: IF(ISRF2.NE.0) THEN CALL SRFC2(IABS(ISRF2),Y(KDIM1),FAUX) END IF IF(ISRF2.NE.0.AND.ISRF2.NE.IY(6).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) ISRF(1)=0 CALL BLOCKS(YB(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2) IF(ISB2.EQ.IY(4)) THEN CALL SRFC2(IABS(IY(6)),YB(KDIM1),FAUX) IF(FAUX(1).NE.0.) THEN C 582 CALL ERROR('582 in CDE: Excluded program branch') C This error should not appear. Contact the authors. END IF C Point XB is situated exactly at the surface IY(6) ISRF(1)=IY(6) CALL BLOCKS(YB(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2) END IF 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. ISRF(1)=IY(6) CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,K1,K2) 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: ISRF(1)=IY(6) CALL BLOCKS(YB(KDIM1),NSRF1,ISRF,IY(4),ISRF2,K1,K2) IF(ISRF2.NE.0) THEN C Surface ISRF2 separates the point XB from the simple C block IY(4): ISRF(1)=ISRF2 CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),K0,K1,K2) 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 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 45 I=1,NDIM Y(I)=YAUX(I) D(I)=DAUX(I) 45 CONTINUE IF(ISB2.EQ.ISBAUX) THEN C Boundary of the complex block has not been crossed C during the last step of numerical integration XB=X DO 46 I=1,NDIM YB(I)=Y(I) DB(I)=D(I) 46 CONTINUE IY(6)=ISRAUX GO TO 49 END IF ISRF(1)=0 CALL BLOCKS * (YAUX(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2) IF(ISRF2.EQ.0) THEN C ISRF2 can be zero only if ISBAUX.EQ.0: IF(ISBAUX.EQ.0) THEN C 5.30: CALL BLOCKS * (Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2) IF(ISRF2.EQ.0) THEN GO TO 49 END IF ELSE WRITE(*,'(20(A,I3))') ' ISRFC1=',ISRAUX, * ' ISB1=',ISBAUX, * ' ICB1=',IY(5), * ' ISRFC2=',ISRF2, * ' ISB2=',IY(4), * ' ICB2=',ICB2 C 583 CALL ERROR('583 in CDE: Excluded program branch') C This error should not appear. Contact the authors. END IF END IF C 5.30: CALL SRFC2(IABS(ISRF2),YAUX(KDIM1),FAUX) 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 C 581 CALL ERROR('581 in CDE: 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. 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 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 C.R.T.5.8.4b. C The curve is parametrized by an independent variable X and evaluated C by the Hermite interpolation from the two given points. The surface C is specified in an implicit way by subroutine SRFC2 which is described C elsewhere, or may coincide with an isosurface of a computed quantity C (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, accidentally, the function has the same C sign at the points X1 and XA, the value at X1 is assumed C to 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 F(2) to F(10) need not be defined. 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 Date: 1996, July 10 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,ITER REAL FAOLD,X,FX,DFX,FB,DFB,FA,DFA,XC,XCB *old REAL ,XCA *!!! REAL XD,XD2,XE,AUX0,AUX1,AUX2,AUX3 *out real xam(50),fam(50),dfam(50),xbm(50),fbm(50),dfbm(50) *out integer icross,jcross,mcross *out data icross,jcross,mcross/0,0,0/ 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 FAOLD=FA XB=X1 X =X1 DO 2 I=KDIM1,KDIM2 YB(I)=Y1(I) DB(I)=D1(I) 2 CONTINUE C In the beginning, both points XB and X are identical with point X1 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 C 585 CALL ERROR('585 in CROSS: 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. END IF C C Iterations: DO 9 ITER=1,50 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 FA=FX ELSE IF(FA*FX.GE.0.) THEN IF(ITER.EQ.1) THEN C Input points X1 and XA are situated at the same side IF(FA.EQ.0.) THEN X=XA FX=FA ELSE XA=X FA=FX END IF ELSE C Here FA, FX and FB should be non-zero due to previous checks XA=XB FA=FB DFA=DFB END IF END IF XB=X FB=FX DFB=DFX C *out xam(iter)=xa *out fam(iter)=fa *out dfam(iter)=dfa *out xbm(iter)=xb *out fbm(iter)=fb *out dfbm(iter)=dfb 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(FB*FAOLD.LT.0.) THEN C Point XA is situated at the other side of the surface than C the point X1 - changing XA and XB X =XA XA=XB XB=X END IF IF((XB-XA)*(XB-X1).LT.0.) THEN C Point XB is closer to X1 than point XA IF(FA*FB.LT.0.) THEN C Points XA and XB cannot be changed C 586 CALL ERROR('586 in CROSS: Reverse order of points') C A pair of close points XA and XB situated at different C sides of the surface has been found, but point XA situated C at the same side as point X1 is far from X1 than point XB. C This error should not appear. Contact the authors. ELSE C Points XA and XB are situated at the same side of the C surface (and, hopefully, FA and FB should be zero) C - changing XA and XB X =XA XA=XB XB=X END IF 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) *out if(icross.eq.0) then *out open(57,file='cross.out') *out end if *out icross=icross+1 *out jcross=jcross+iter-1 *out mcross=max0(mcross,iter-1) *out if(mod(icross,100).eq.0) then *out write(57,*) icross,jcross,mcross,float(jcross)/float(icross) *out end if *out if(iter.gt.20) then *out open(58,file='error.out') *out do 7 I=1,iter *out write(58,*) xam(i),xbm(i),fam(i),fbm(i),dfam(i),dfbm(i) *out7 continue *out endfile (58) *out backspace(58) *out pause 'Error new in CROSS: More than 20 iterations' *out end if RETURN END IF C cccc (a) Odstranit *old cccc (b) Odstranit *old, odstranit *??? cccc (c) Odstranit *old, odstranit *???, odstranit *new (nedulezite) cccc (d) Vratit *old, vratit *???, vratit *new cccc (e) Odstranit *!!! cccc (f) Odstranit *!!!, odstranit *old cccc C New approximation: IF(MOD(ITER,2).EQ.1) THEN C Regula falsi: X=(FA*XB-FB*XA)/(FA-FB) C!!! CHECKING THE ACCURACY OF THE REGULA FALSI METHOD: *!!! IF(ITER.EQ.1) THEN C!!! COEFFICIENTS OF THE CUBIC TAYLOR EXPANSION AT XC=(XA+XB)/2.: *!!! XD=XB-XA *!!! XD2=XD*XD/2. *!!! AUX3=(DFB+DFA)/2. *!!! AUX2=(DFB-DFA)/XD *!!! AUX1=(FB-FA)/XD *!!! AUX0=(FB+FA)/2.-AUX2*XD2 *!!! AUX1=1.5*AUX1-0.5*AUX3 *!!! AUX3=(AUX3-AUX1)/XD2 C!!! COEFFICIENTS OF THE QUADRATIC TAYLOR EXPANSION AT X: *!!! XD=X-XC *!!! XD=X-XC *!!! XD2=XD*XD/2. *!!! XE=SIGN(1.,AUX1) *!!! AUX0=AUX0+AUX1*XD+AUX2*XD2+AUX3*XD*XD2/3. *!!! AUX1=AUX1+AUX2*XD+AUX3*XD2 *!!! IF(ABS(AUX0).GT.0.5*ERR*ABS(AUX1)) THEN C!!! THE ACCURACY SHOULD BE IMPROVED: *!!! AUX2=AUX2+AUX3*XD *!!! AUX3=AUX1**2-2.*AUX0*AUX2 C!!! HERE 2**(24/2)=4096 ASSUMES 24 BIT FLOATING-POINT ACCURACY *!!! IF(4096.*ABS(AUX0*AUX2).GT.AUX1**2.AND.AUX3.GE.0.) THEN *!!! XE=X+(XE*SQRT(AUX3)-AUX1)/AUX2 *!!! ELSE *!!! XE=X-AUX0/AUX1 *!!! END IF *!!! IF((XE-XA)*(XE-XB).LE.0.) THEN *!!! X=XE *!!! END IF *!!! END IF *!!! END IF ELSE C Modified Newton-Raphson: XC=(XA+XB)/2. *old XCA=XA-FA/DFA+SIGN(ERR/50.,XA-XB) XCB=XB-FB/DFB+SIGN(ERR/50.,XA-XB) *old IF((XCA-XC)*(XCA-XB).LT.0.) THEN *old IF((XCB-XC)*(XCB-XB).LT.0.) THEN *old IF(ABS(XCA-XB).LT.ABS(XCB-XB)) THEN *old X=XCA *old ELSE *old X=XCB *old END IF *old ELSE *old X=XCA *old END IF *old 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) *new IF((X-XC)*(X-XB).GE.0.) THEN *new X=XC *new pause 'New warning in CROSS' *new END IF *??? END IF END IF *old 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 *out open(58,file='error.out') *out do 8 iter=1,50 *out write(58,*) xam(iter),xbm(iter),fam(iter),fbm(iter), *out * dfam(iter),dfbm(iter) *out8 continue C 584 CALL ERROR('584 in CROSS: Too many iterations') C More than 50 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 data C DCRT - C less than rounding errors. C Usually, this error should not appear. Contact the authors. END C C======================================================================= 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 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 designed to evaluate symmetric matrix by vector C product. It may be, e.g., called after the invocation of the METRIC C subroutine to transform the covariant vectors to the contravariant C ones and vice versa. 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