C
C Subroutine file 'wan.for' to compute quantities along a ray necessary C for computation of the Green function by means of coupling ray-theory C in weakly anisotropic models without interfaces. C C Version: 5.20 C Date: 1998, November 11 C C Coded by: Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: bulant@seis.karlov.mff.cuni.cz C C======================================================================= C SUBROUTINE WAN(LU1,LU2,LU3,GREEN,MGREEN,RPTS,IPTS,MAXPTS,NQ) C C---------------------------------------------------------------------- INTEGER LU1,LU2,LU3 INTEGER MGREEN,MAXPTS,NQ INTEGER NQPTS PARAMETER (NQPTS=21) REAL GREEN(MGREEN),RPTS(NQPTS,MAXPTS/NQPTS) INTEGER IPTS(NQPTS,MAXPTS/NQPTS) C ------------------------------------------------------------------ C Input: C LU1 ... Number of the logical unit connected to the CRT output C file with quantities along rays. C LU2 ... Number of the logical unit connected to the CRT output C file with quantities at the initial points of rays. C LU3 ... Number of the logical unit connected to the CRT output C file with quantities at the storing surface. C MGREEN..Dimension of an output array GREEN. C MAXPTS..Maximum number of records in arrays RPTS and IPTS. C Output: C GREEN...Array containing the Green function for the given ray C and for all the frequencies: C GREEN(1)... Travel time between receiver and source. C GREEN(2)... Imaginary part of the complex-valued travel time C between receiver and source due to attenuation. C GREEN(3:8)... Coordinates of the receiver and coordinates C of the source. C GREEN(9:14)... Derivatives of the travel time with respect C to the coordinates of the receiver and coordinates of the C source. C GREEN(15:(14+NF*18)) ... C for P wave once, for S wave NF times following 18 numbers, C specifying 1 000 000 times enlarged amplitude of the C Green function: contravariant components of the complex- C valued 3*3 matrix Gij in model coordinates, where the C first subscript corresponds to the receiver and the second C subscript corresponds to the source. The components are C ordered as C Re(G11),Im(G11),Re(G21),Im(G21),Re(G31),Im(G31), C Re(G12),Im(G12),Re(G22),Im(G22),Re(G32),Im(G32), C Re(G13),Im(G13),Re(G23),Im(G23),Re(G33),Im(G33). C NQ ... Number of records stored in the array GREEN. C IPTS,RPTS... Auxiliary arrays, redimensioned in each invocation. C C....................................................................... C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C....................................................................... C C Auxiliary storage locations: INTEGER MPTS REAL PIGRA(2,2),PIGIA(2,2),PIGR(2,2),PIGI(2,2) REAL PI PARAMETER (PI=3.1415926) REAL GREENA(32) REAL TTCOR,FREQ,SIGRAY,GAMA,AUX0,AUX1,AUX2,COOR(3) INTEGER NPTS INTEGER I,I1,I2,J INTEGER NFFT,NF REAL DT,FMIN,FMAX,DF,OF C IPTS,RPTS... Quantities in the points on the ray: C IPTS(1,I)... Index of the I-th point, zero for points C added to the ray by interpolation. C RPTS(2,I)... Travel time in I-th point. C RPTS(3-5,I).. Coordinates of the point. C RPTS(6-8,I).. Slowness vector in the point. C RPTS(9-11,I) Polarization vector in the point. C IPTS(12,I)... Index of complex block. C RPTS(13,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qP wave. C RPTS(14,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS2 (faster) wave. C RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS1 (slower) wave. C RPTS(16-19,I) Eigenvectors projected to the plane defined by C polarization vectors. C RPTS(20,I)... Angular difference of eigenvectors of Christoffel C matrix corresponding to qS2 wave in I-th and C in (I-1)-st point. C RPTS(21,I)... Time integral of the difference between qS1 and C qS2 eigenvalues along the ray from the (I-1)-st C to the I-th point. C C----------------------------------------------------------------------- C CALL RSEP3I('NFFT',NFFT,1) CALL RSEP3R('DT ',DT ,1.) CALL RSEP3R('FMIN',FMIN,0.) CALL RSEP3R('FMAX',FMAX,0.5/DT) CALL RSEP3R('DF ',DF ,1./(FLOAT(NFFT)*DT)) CALL RSEP3R('OF ',OF ,DF*NINT(FMIN/DF)) CALL RSEP3I('NF ',NF ,NINT((FMAX-OF)/DF)+1) C C Reading in the quantities stored in individual points on C a ray, computing all the auxiliary quantities in the points, C and, if necessary, adding new points on the ray by interpolation: MPTS=MAXPTS/NQPTS NPTS=0 CALL WAPTS(LU1,LU2,LU3,RPTS,IPTS,MPTS,NPTS) IF (NPTS.EQ.0) THEN C End of rays: NQ=0 RETURN ENDIF C C Computing the values of travel time corrections along the ray: TTCOR=0. DO 30, I1=1,NPTS-1 I2=I1+1 IF (IPTS(12,I2).GT.0.) THEN C P wave: TTCOR=TTCOR+(1./SQRT(RPTS(13,I2))+1./SQRT(RPTS(13,I1)))*0.5* * (RPTS(2,I2)-RPTS(2,I1)) ELSE C S wave: TTCOR=TTCOR+(1./SQRT((RPTS(14,I2)+RPTS(15,I2))*0.5)+ * 1./SQRT((RPTS(14,I1)+RPTS(15,I1))*0.5) )*0.5* * (RPTS(2,I2)-RPTS(2,I1)) ENDIF 30 CONTINUE C IF (IPTS(12,I2).GT.0.) THEN C P wave: !!! OPRAVIT PO ZARAZENI KONVERTOVANYCH VLN !!! Y(1)=TTCOR+YI(1) CALL AP21(GREEN) NQ=32 GOTO 91 ENDIF C C S wave: C Loop over the frequencies: DO 90, I2=0,NF-1 FREQ=OF+I2*DF PIGR(1,1)=1. PIGR(2,1)=0. PIGR(1,2)=0. PIGR(2,2)=1. PIGI(1,1)=0. PIGI(2,1)=0. PIGI(1,2)=0. PIGI(2,2)=0. SIGRAY=0. C Computing the propagator matrix PiGe along the ray: C Loop along points on the ray: DO 40, I1=2,NPTS GAMA=PI*FREQ*RPTS(21,I1) AUX0=SQRT(RPTS(20,I1)**2 + GAMA**2) AUX1=COS(AUX0) IF (AUX0.EQ.0.) THEN AUX2=1. ELSE AUX2=SIN(AUX0)/AUX0 ENDIF C Matrix for this step along the ray: PIGRA(1,1)= AUX1 PIGRA(2,1)=-RPTS(20,I1)*AUX2 PIGRA(1,2)=-PIGRA(2,1) PIGRA(2,2)= AUX1 PIGIA(1,1)=-GAMA*AUX2 PIGIA(2,1)= 0. PIGIA(1,2)= 0. PIGIA(2,2)=-PIGIA(1,1) C Matrix for all the steps along the ray to current point: CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI) 40 CONTINUE C C Computing the Green function: Y(1)=TTCOR+YI(1) G11=RPTS(16,1) G12=RPTS(17,1) G21=RPTS(18,1) G22=RPTS(19,1) PIR11=PIGR(1,1)*G11+PIGR(1,2)*G21 PIR21=PIGR(2,1)*G11+PIGR(2,2)*G21 PIR12=PIGR(1,1)*G12+PIGR(1,2)*G22 PIR22=PIGR(2,1)*G12+PIGR(2,2)*G22 PII11=PIGI(1,1)*G11+PIGI(1,2)*G21 PII21=PIGI(2,1)*G11+PIGI(2,2)*G21 PII12=PIGI(1,1)*G12+PIGI(1,2)*G22 PII22=PIGI(2,1)*G12+PIGI(2,2)*G22 G11=RPTS(16,NPTS) G21=RPTS(17,NPTS) G12=RPTS(18,NPTS) G22=RPTS(19,NPTS) Y(28)=G11*PIR11+G12*PIR21 Y(29)=G11*PII11+G12*PII21 Y(30)=G21*PIR11+G22*PIR21 Y(31)=G21*PII11+G22*PII21 Y(32)=G11*PIR12+G12*PIR22 Y(33)=G11*PII12+G12*PII22 Y(34)=G21*PIR12+G22*PIR22 Y(35)=G21*PII12+G22*PII22 C CALL AP21(GREENA) IF (I2.EQ.0) THEN DO 50 I=1,14 GREEN(I)=GREENA(I) 50 CONTINUE NQ=14 ENDIF J=I2*18 DO 60 I=15,32 GREEN(J+I)=GREENA(I) 60 CONTINUE NQ=NQ+18 C 90 CONTINUE 91 CONTINUE RETURN END C================================================================================ C SUBROUTINE WAPTS(LU1,LU2,LU3,RPTS,IPTS,MPTS,NPTS) C C---------------------------------------------------------------------- C Subroutine to read in the quantities stored in individual points C on the ray, to compute all the auxiliary quantities in the points, C and, if necessary, to add new points on the rays by interpolation. C Reading the files with points on the rays is done by a simple C invocation of subroutine WAREAD. C Computation of elastic parameters is completed by invocation C of subroutine PARM3 of file 'parm.for'. C Then Christoffel matrix is evaluated and its eigenvalues and C eigenvectors are computed by subroutine EIGEN of file 'eigen.for'. C In the next step the angular difference DELTFI is computed for C each subinterval along the ray. If the value of DELTFI is greater than C prescribed limit new points are added using subroutine HIVD2 of C the file 'means.for'. C INTEGER NPTS,NQPTS,MPTS PARAMETER (NQPTS=21) REAL RPTS(NQPTS,MPTS) INTEGER IPTS(NQPTS,MPTS) C Input: C MPTS ... Dimension of arrays RPTS and IPTS. C Output: C RPTS,IPTS ... The arrays are filled with all the quantities for C single two-point ray during one invocation of this C subroutine: C IPTS(1,I) ... Index of the I-th point, zero for points C added to the ray by interpolation. C RPTS(2,I) ... Travel time in I-th point. C RPTS(3-5,I) . Coordinates of the point. C RPTS(6-8,I) . Slowness vector in the point. C RPTS(9-11,I) Polarization vector in the point. C IPTS(12,I) .. Index of complex block. C RPTS(13,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qP wave. C RPTS(14,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS2 (faster) wave. C RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS1 (slower) wave. C RPTS(16-19,I) Eigenvectors projected to the plane defined by C polarization vectors. C RPTS(20,I) .. Angular difference of eigenvectors of Christoffel C matrix corresponding to qS2 wave in I-th and C in (I-1)-st point. C RPTS(21,I) .. Time integral of the difference between qS1 and C qS2 eigenvalues along the ray from the (I-1)-st C to the I-th point. C NPTS ... Number of points stored in RPTS (IPTS). C C External functions required: EXTERNAL WACHAN REAL WACHAN C....................................................................... C C Auxiliary storage locations: INTEGER NPTSE INTEGER MQANT,MNEWP,NNEWP PARAMETER (MNEWP=5) PARAMETER (MQANT=17) REAL ROLDP(MQANT),RNEWP(MQANT,MNEWP) INTEGER KNEWP(MNEWP) INTEGER MAXDIV PARAMETER (MAXDIV=MNEWP-1) REAL AA(10,21),RHO,QQ(21) REAL EE(9),DER(9) REAL DEFI,DELTFI c PARAMETER (DEFI=0.01) PARAMETER (DEFI=0.03) c PARAMETER (DEFI=100.) REAL TTOLD,CROLD(3),DCROLD(3) REAL TTNEW,CRNEW(3),DCRNEW(3) INTEGER I1,I2 SAVE NNEWP DATA NNEWP/0/,KNEWP/MNEWP*0/ C C ROLDP(I),RNEWP(I,J) C I=1 ... Travel time. C 2-4 ... Coordinates. C 5-7 ... Slowness vector. C 8-10 ... Polarization vector E1. C 11 ... qP eigenvalue of Christoffel matrix. C 12 ... qS1 eigenvalue of Christoffel matrix. C 13 ... qS2 eigenvalue of Christoffel matrix. C 14 ... First component of the qS1 eigenvector C projected to the plane perpendicular to the ray, C defined by two basis vectors of ray-centered C coordinate system. C 15 ... Second component of the qS1 eigenvector. C 16 ... First component of the qS2 eigenvector. C 17 ... Second component of the qS2 eigenvector. C KNEWP...Division index of points in RNEWP. C AA.. Values, first and second partial derivatives of real C parts of 21 reduced (divided by the density) elastic C parameters. The order of the value, first and second C partial derivatives of each parameter Aij is: C Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33. C The order of parameters (second array index) is: C A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45, C A55,A16,A26,A36,A46,A56,A66. C RHO... Density at the given point. C QQ ... Imaginary parts of 21 reduced elastic parameters at the C given point, ordered as C Q11,Q12,Q22,Q13,Q23,Q33,Q14,Q24,Q34,Q44,Q15,Q25,Q35,Q45, C Q55,Q16,Q26,Q36,Q46,Q56,Q66. C EE ... Eigenvectors of the Christoffel matrix. C DER ... Derivatives dx/dt. C MAXDIV .. The distance between two points on the rays must not C be smaller than original distance between the points C computed by CRT divided by MAXDIV when adding new C points on the ray. C DEFI .. Maximum allowed angular change for eigenvectors of C Christoffel matrix between neighboring points on the ray. C TTOLD...Travel time in the old point stored for interpolation. C CROLD...Coordinates of the old point stored for interpolation. C DCROLD..Derivatives in the old point stored for interpolation. C TTNEW...Travel time in the new point stored for interpolation. C CRNEW...Coordinates of the new point stored for interpolation. C DCRNEW..Derivatives in the new point stored for interpolation. C NPTSE...Number of points along the ray, where all the quantities C have been checked. C C----------------------------------------------------------------------- C Reading the quantities computed by CRT for one ray: CALL WAREAD(LU1,LU2,LU3,RPTS,IPTS,MPTS,NPTS) IF (NPTS.EQ.0) C End of rays: * RETURN IF (NPTS.LT.2) THEN C WAN-01 CALL ERROR('WAN-01: A ray formed by single point') C This error should not appear. C Each ray should be represented by at least two points. ENDIF C C C Reading the material parameters in the initial point: CALL PARM3(IPTS(12,1),RPTS(3,1),AA,RHO,QQ) C Computing eigenvectors and eigenvalues of C the Christoffel matrix in the initial point: CALL WACHRI(RPTS(6,1),RPTS(7,1),RPTS(8,1), * AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7), * AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14), * AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21), * RPTS(13,1),RPTS(14,1),RPTS(15,1),EE,DER) IF (ABS(RPTS(15,1)-RPTS(14,1)).LT.0.000001) THEN C Isotropic case, no projection of eigenvectors: RPTS(16,1)=0. RPTS(17,1)=0. ELSE C Projecting the qS eigenvectors to the plane C perpendicular to the ray: CALL WAPROJ(RPTS(6,1),RPTS(7,1),RPTS(8,1), * RPTS(9,1),RPTS(10,1),RPTS(11,1), * EE(4),EE(5),EE(6),EE(7),EE(8),EE(9), * RPTS(16,1),RPTS(17,1),RPTS(18,1),RPTS(19,1)) ENDIF RPTS(20,1)=0. RPTS(21,1)=0. C C Quantities for future possible interpolation: TTOLD =RPTS(2,1) CROLD(1) =RPTS(3,1) CROLD(2) =RPTS(4,1) CROLD(3) =RPTS(5,1) IF (IPTS(12,1).GT.0) THEN C P waves: DCROLD(1)=DER(1) DCROLD(2)=DER(2) DCROLD(3)=DER(3) ELSE C S waves: DCROLD(1)=(DER(4)+DER(7))*0.5 DCROLD(2)=(DER(5)+DER(8))*0.5 DCROLD(3)=(DER(6)+DER(9))*0.5 ENDIF C ROLDP( 1)=RPTS( 2,1) ROLDP( 2)=RPTS( 3,1) ROLDP( 3)=RPTS( 4,1) ROLDP( 4)=RPTS( 5,1) ROLDP( 5)=RPTS( 6,1) ROLDP( 6)=RPTS( 7,1) ROLDP( 7)=RPTS( 8,1) ROLDP( 8)=RPTS( 9,1) ROLDP( 9)=RPTS(10,1) ROLDP(10)=RPTS(11,1) ROLDP(11)=RPTS(13,1) ROLDP(12)=RPTS(14,1) ROLDP(13)=RPTS(15,1) ROLDP(14)=RPTS(16,1) ROLDP(15)=RPTS(17,1) ROLDP(16)=RPTS(18,1) ROLDP(17)=RPTS(19,1) NNEWP=0 NPTSE=1 C C C Loop along the ray: 5 CONTINUE IF (NNEWP.LE.0) THEN C Reading the material parameters in a new point on the ray: IF (NNEWP.LT.0) THEN C WAN-02 CALL ERROR('WAN-02: Negative number of points') C This error should not appear. C The number of new points should be zero or positive integer. ENDIF IF (NPTSE+1.GT.NPTS) THEN C WAN-03 CALL ERROR('WAN-03: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C of in the include file ram.inc. ENDIF I1=NPTSE+1 RNEWP( 1,1)=RPTS( 2,I1) RNEWP( 2,1)=RPTS( 3,I1) RNEWP( 3,1)=RPTS( 4,I1) RNEWP( 4,1)=RPTS( 5,I1) RNEWP( 5,1)=RPTS( 6,I1) RNEWP( 6,1)=RPTS( 7,I1) RNEWP( 7,1)=RPTS( 8,I1) RNEWP( 8,1)=RPTS( 9,I1) RNEWP( 9,1)=RPTS(10,I1) RNEWP(10,1)=RPTS(11,I1) NNEWP=1 CALL PARM3(IPTS(12,I1),RNEWP(2,1),AA,RHO,QQ) C Computing the Christoffel matrix in the new point: CALL WACHRI(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1), * AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7), * AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14), * AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20), * AA(1,21),RNEWP(11,1),RNEWP(12,1),RNEWP(13,1),EE,DER) C Projecting the qS eigenvectors to the plane C perpendicular to the ray: CALL WAPROJ(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1),RNEWP(8,1), * RNEWP(9,1),RNEWP(10,1),EE(4),EE(5),EE(6),EE(7),EE(8),EE(9), * RNEWP(14,1),RNEWP(15,1),RNEWP(16,1),RNEWP(17,1)) C C Quantities for future possible interpolation: TTNEW=RNEWP(1,1) CRNEW(1)=RNEWP(2,1) CRNEW(2)=RNEWP(3,1) CRNEW(3)=RNEWP(4,1) IF (IPTS(12,NPTSE+1).GT.0) THEN C P waves: DCRNEW(1)=DER(1) DCRNEW(2)=DER(2) DCRNEW(3)=DER(3) ELSE C S waves: DCRNEW(1)=(DER(4)+DER(7))*0.5 DCRNEW(2)=(DER(5)+DER(8))*0.5 DCRNEW(3)=(DER(6)+DER(9))*0.5 ENDIF C IF (IPTS(12,NPTSE).NE.IPTS(12,NPTSE+1)) THEN C The ray is crossing an interface: C WAN-04 CALL ERROR('WAN-04: A ray is crossing an interface') C This part is not yet coded in this version. ENDIF C IF (ABS(RNEWP(12,1)-RNEWP(13,1)).LT.0.000001) THEN C New point is isotropic, qS eigenvalues are the same, C no change in eigenvectors: RNEWP(14,1)=ROLDP(14) RNEWP(15,1)=ROLDP(15) RNEWP(16,1)=ROLDP(16) RNEWP(17,1)=ROLDP(17) DELTFI=0. GOTO 20 ENDIF IF (ABS(ROLDP(12)-ROLDP(13)).LT.0.000001) THEN C Old point is isotropic (new point is anizotropic): IF ((ROLDP(14).EQ.0.).AND.(ROLDP(15).EQ.0.)) THEN C New point is the first anizotropic point on the ray, C angular change is not to be computed: DELTFI=0. GOTO 20 ENDIF ENDIF ENDIF C C Computing the angular change in eigenvectors, adding new points C on the ray if necessary: DO 10, I1=1,NNEWP DELTFI=WACHAN(ROLDP(14),RNEWP(11,I1)) IF (DELTFI.LE.DEFI) THEN cc IF (DELTFI.LE.( 24/(2.*PI*(OF+NF*DF)*RPTS(2,NPTS)* CC (RNEWP(12,NNEWP)-RNEWP(13,NNEWP)+ROLDP(12)-ROLDP(13)) ) THEN C Angular change is less than prescribed limit, the I1-th point C of array RNEWP will be used as the next point on the ray: NNEWP=I1 GOTO 20 ENDIF 10 CONTINUE C Angular change is greater than prescribed limit for all the points C of array RNEWP. 15 CONTINUE C Loop for adding new points on the ray until the angular change is C small enough: IF (KNEWP(NNEWP).GE.MAXDIV) THEN C WAN-05 CALL ERROR('WAN-05: Maximum number of divisions exceeded') C The angular change of eigenvectors in two consecutive points C is too big. More than MAXDIV divisions of this interval is C needed to keep the change under the prescribed limit. Try to C decrease the parameter STEP, or to C increase MAXDIV or DEFI. ENDIF C Adding new point to the ray: KNEWP(NNEWP)=KNEWP(NNEWP)+1 NNEWP=NNEWP+1 IF (NNEWP.GT.MNEWP) THEN C WAN-06 CALL ERROR('WAN-06: Array KNEWP too small') C This error should not appear, error05 should appear instead. ENDIF C Travel time: RNEWP(1,NNEWP)=(ROLDP(1)+RNEWP(1,NNEWP-1))*0.5 C Coordinates: CALL HIVD2(3,TTOLD,CROLD,DCROLD,TTNEW,CRNEW,DCRNEW, * RNEWP(1,NNEWP),RNEWP(2,NNEWP),DER) C Slowness vector: RNEWP(5,NNEWP)=(ROLDP(5)+RNEWP(5,NNEWP-1))*0.5 RNEWP(6,NNEWP)=(ROLDP(6)+RNEWP(6,NNEWP-1))*0.5 RNEWP(7,NNEWP)=(ROLDP(7)+RNEWP(7,NNEWP-1))*0.5 C Polarization vector: (could be interpolated by HIVD2 with C derivatives obtained from FCT from file 'raycb.for') RNEWP(8,NNEWP) =(ROLDP(8) +RNEWP(8,NNEWP-1)) *0.5 RNEWP(9,NNEWP) =(ROLDP(9) +RNEWP(9,NNEWP-1)) *0.5 RNEWP(10,NNEWP)=(ROLDP(10)+RNEWP(10,NNEWP-1))*0.5 C Material parameters: CALL PARM3(IPTS(12,NPTSE+1),RNEWP(2,NNEWP),AA,RHO,QQ) C Christoffel matrix and eigenvalues: CALL WACHRI(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP), * AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7), * AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14), * AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21), * RNEWP(11,NNEWP),RNEWP(12,NNEWP),RNEWP(13,NNEWP),EE,DER) C Projection of the qS eigenvectors to the plane C perpendicular to the ray: CALL WAPROJ(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP), * RNEWP(8,NNEWP),RNEWP(9,NNEWP),RNEWP(10,NNEWP), * EE(4),EE(5),EE(6),EE(7),EE(8),EE(9), * RNEWP(14,NNEWP),RNEWP(15,NNEWP),RNEWP(16,NNEWP),RNEWP(17,NNEWP)) C Index of the division: KNEWP(NNEWP)=KNEWP(NNEWP-1) C IF (ABS(RNEWP(12,NNEWP)-RNEWP(13,NNEWP)).LT.0.000001) THEN C Isotropic case, qS eigenvalues are the same, C no change in eigenvectors: RNEWP(14,NNEWP)=ROLDP(14) RNEWP(15,NNEWP)=ROLDP(15) RNEWP(16,NNEWP)=ROLDP(16) RNEWP(17,NNEWP)=ROLDP(17) DELTFI=0. ELSE DELTFI=WACHAN(ROLDP(14),RNEWP(11,NNEWP)) ENDIF IF (DELTFI.LE.DEFI) THEN C Angular change is less than prescribed limit, this point C of array RNEWP will be used as the next point on the ray: GOTO 20 ELSE C Angular change is greater than prescribed limit, adding C a new point to the ray: GOTO 15 ENDIF C End of the loop. C 20 CONTINUE C Angular change DELTFI for points ROLDP, RNEWP(i,NNEWP) is less C than prescribed limit. Recording the computed quantities. NPTSE=NPTSE+1 IF (NNEWP.NE.1) THEN C The new point was computed by interpolation. C Shifting the array RPTS: NPTS=NPTS+1 IF (NPTS.GT.MPTS) THEN C WAN-07 CALL ERROR('WAN-07: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C of in the include file ram.inc. ENDIF DO 31, I1=NPTS,NPTSE+1,-1 DO 30, I2=1,NQPTS RPTS(I2,I1)=RPTS(I2,I1-1) 30 CONTINUE 31 CONTINUE C Recording interpolated quantities: RPTS( 1,NPTSE)=0 RPTS( 2,NPTSE)=RNEWP( 1,NNEWP) RPTS( 3,NPTSE)=RNEWP( 2,NNEWP) RPTS( 4,NPTSE)=RNEWP( 3,NNEWP) RPTS( 5,NPTSE)=RNEWP( 4,NNEWP) RPTS( 6,NPTSE)=RNEWP( 5,NNEWP) RPTS( 7,NPTSE)=RNEWP( 6,NNEWP) RPTS( 8,NPTSE)=RNEWP( 7,NNEWP) RPTS( 9,NPTSE)=RNEWP( 8,NNEWP) RPTS(10,NPTSE)=RNEWP( 9,NNEWP) RPTS(11,NPTSE)=RNEWP(10,NNEWP) IPTS(12,NPTSE)=IPTS(12,NPTSE-1) ENDIF C Recording quantities for computation of anisotropic corrections: RPTS(13,NPTSE)=RNEWP(11,NNEWP) RPTS(14,NPTSE)=RNEWP(12,NNEWP) RPTS(15,NPTSE)=RNEWP(13,NNEWP) RPTS(16,NPTSE)=RNEWP(14,NNEWP) RPTS(17,NPTSE)=RNEWP(15,NNEWP) RPTS(18,NPTSE)=RNEWP(16,NNEWP) RPTS(19,NPTSE)=RNEWP(17,NNEWP) RPTS(20,NPTSE)=DELTFI C-LK RPTS(21,NPTSE)= C-LK * (RNEWP(12,NNEWP)-RNEWP(13,NNEWP)+ROLDP(12)-ROLDP(13))*0.25 C-LK * *(RNEWP(1,NNEWP)-ROLDP(1)) RPTS(21,NPTSE)= * (1./SQRT(RNEWP(13,NNEWP))-1./SQRT(RNEWP(12,NNEWP)) * +1./SQRT(ROLDP(13))-1./SQRT(ROLDP(12)))*0.5 * *(RNEWP(1,NNEWP)-ROLDP(1)) C IF (NPTSE.LT.NPTS) THEN C Continuing with the next point on the ray: DO 100, I1=1,MQANT ROLDP(I1)=RNEWP(I1,NNEWP) 100 CONTINUE KNEWP(NNEWP)=0 NNEWP=NNEWP-1 GOTO 5 ENDIF RETURN END C C======================================================================= C SUBROUTINE WAREAD(LU1,LU2,LU3,RPTS,IPTS,MPTS,NPTS) INTEGER LU1,LU2,LU3 C C---------------------------------------------------------------------- C Subroutine reads the unformatted output of program CRT. C Reading the output files is completed by a simple invocation of C subroutine AP00 of file 'ap.for'. Subroutine reads the quantities C stored in individual points on the ray into array RPTS of common C block PTS. C The subroutine expects, that files LU and LUI are already opened C and that at least one point was already read by calling AP00. C INTEGER NPTS,NQPTS,MPTS PARAMETER (NQPTS=21) REAL RPTS(NQPTS,MPTS) INTEGER IPTS(NQPTS,MPTS) C Input: C MPTS ... Dimension of arrays RPTS and IPTS. C Output: C RPTS,IPTS ... The arrays are filled with the quantities stored in C crt output files for single two-point ray during one C invocation of this subroutine: C IPTS(1,I) ... Index of the I-th point, zero for points C added to the ray by interpolation. C RPTS(2,I) ... Travel time in I-th point. C RPTS(3-5,I) . Coordinates of the point. C RPTS(6-8,I) . Slowness vector in the point. C RPTS(9-11,I) Polarization vector in the point. C IPTS(12,I) .. Index of complex block. C NPTS ... Number of points stored in RPTS (IPTS). C C ........................... C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C None of the storage locations of the common block are altered. C....................................................................... C C Auxiliary storage locations: C C----------------------------------------------------------------------- C 5 CONTINUE CALL AP00(LU1,LU2,LU3) IF (IWAVE.LT.1) THEN C End of rays: NPTS=0 RETURN ENDIF IF (IREC.LE.0) C Only two-point rays are to be written into RPTS (IPTS). * GOTO 5 IF ((IPT.NE.0).AND.(IPT.NE.1)) C End of previous ray behind the reference surface: * GOTO 5 C C Initial point of a new two-point ray: NPTS=1 C Index of a point: IPTS(1,NPTS)=NPTS C Travel time: RPTS(2,NPTS)=YI(1) C Coordinates: RPTS(3,NPTS)=YI(3) RPTS(4,NPTS)=YI(4) RPTS(5,NPTS)=YI(5) C Slowness vector: RPTS(6,NPTS)=YI(6) RPTS(7,NPTS)=YI(7) RPTS(8,NPTS)=YI(8) C Polarization vector: RPTS(9,NPTS)= YI(9) RPTS(10,NPTS)=YI(10) RPTS(11,NPTS)=YI(11) C Index of complex block: IPTS(12,NPTS)=ICB1I C C 20 CONTINUE C New point: IF (YF(1).LT.Y(1)) THEN C The point along the ray is before the reference surface, C recording the point along the ray: NPTS=NPTS+1 IF (NPTS.GT.MPTS) THEN C WAN-08 CALL ERROR('WAN-08: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C of in the include file ram.inc. ENDIF C Index of a point: IPTS(1,NPTS)=NPTS C Travel time: RPTS(2,NPTS)=YF(1) C Coordinates: RPTS(3,NPTS)=YF(3) RPTS(4,NPTS)=YF(4) RPTS(5,NPTS)=YF(5) C Slowness vector: RPTS(6,NPTS)=YF(6) RPTS(7,NPTS)=YF(7) RPTS(8,NPTS)=YF(8) C Polarization vector: RPTS(9,NPTS)= YF(9) RPTS(10,NPTS)=YF(10) RPTS(11,NPTS)=YF(11) C Index of complex block: IPTS(12,NPTS)=ICB1F C C Reading the results of the complete ray tracing: CALL AP00(LU1,LU2,LU3) IF ((IWAVE.LT.1).OR.(IPT.LE.1)) THEN C This should not happen, the ray must reach the C reference surface. C WAN-09 CALL ERROR('WAN-09: The ray missed the reference surface') C As only the two-point rays are considered by the subroutine C "WAN", each of the rays should pass the reference surface. C Check, whether you have specified right names C of the input files with points along rays, C points at their initial points and points at the C reference surface in file CRT, C or whether you have correctly specified its name C CRTOUT. ENDIF C GOTO 20 ENDIF C C The point along the ray is at or above the reference surface, C recording the point at the reference surface: NPTS=NPTS+1 IF (NPTS.GT.MPTS) THEN C WAN-10 CALL ERROR('WAN-10: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C of in the include file ram.inc. ENDIF C Index of a point: IPTS(1,NPTS)=NPTS C Travel time: RPTS(2,NPTS)=Y(1) C Coordinates: RPTS(3,NPTS)=Y(3) RPTS(4,NPTS)=Y(4) RPTS(5,NPTS)=Y(5) C Slowness vector: RPTS(6,NPTS)=Y(6) RPTS(7,NPTS)=Y(7) RPTS(8,NPTS)=Y(8) C Polarization vector: RPTS(9,NPTS)= Y(9) RPTS(10,NPTS)=Y(10) RPTS(11,NPTS)=Y(11) C Index of complex block: IPTS(12,NPTS)=ICB1 C RETURN END C C======================================================================= C SUBROUTINE WACHRI(P1,P2,P3,B11,B12,B22,B13,B23,B33, * B14,B24,B34,B44,B15,B25,B35,B45,B55, * B16,B26,B36,B46,B56,B66, * G1,G2,G3,EE,DER) C C---------------------------------------------------------------------- C Subroutine to compute the Christoffel matrix, its eigenvalues C and eigenvectors. REAL P1,P2,P3 REAL B11,B12,B22,B13,B23,B33,B14,B24,B34,B44 REAL B15,B25,B35,B45,B55,B16,B26,B36,B46,B56,B66 REAL G1,G2,G3,EE(9),DER(9) C C Input: C P1,P2,P3... Slowness vector. C Bii ... Values of real parts of 21 reduced C (divided by the density) elastic parameters. C Output: C G1,G2,G3 ... Eigenvalues of the Christoffel matrix. C EE ... Eigenvectors of the Christoffel matrix. C DER ... Derivatives dx/dt=dH/dp=Aijkl Ej Ek pl stored columnwise C for qP, qS1 and qS2 waves. C C----------------------------------------------------------------------- C REAL A11,A12,A13,A14,A21,A22,A23,A24,A31,A32,A33,A34,A44 REAL A15,A25,A35,A45,A55,A16,A26,A36,A46,A56,A66 REAL A1111,A2111,A3111,A1211,A2211,A3211,A1311,A2311,A3311 REAL A1121,A2121,A3121,A1221,A2221,A3221,A1321,A2321,A3321 REAL A1131,A2131,A3131,A1231,A2231,A3231,A1331,A2331,A3331 REAL A1112,A2112,A3112,A1212,A2212,A3212,A1312,A2312,A3312 REAL A1122,A2122,A3122,A1222,A2222,A3222,A1322,A2322,A3322 REAL A1132,A2132,A3132,A1232,A2232,A3232,A1332,A2332,A3332 REAL A1113,A2113,A3113,A1213,A2213,A3213,A1313,A2313,A3313 REAL A1123,A2123,A3123,A1223,A2223,A3223,A1323,A2323,A3323 REAL A1133,A2133,A3133,A1233,A2233,A3233,A1333,A2333,A3333 EQUIVALENCE (A11,A1111) EQUIVALENCE (A22,A2222) EQUIVALENCE (A33,A3333) EQUIVALENCE (A16,A1112,A1121,A1211,A2111) EQUIVALENCE (A26,A2221,A2212,A2122,A1222) EQUIVALENCE (A15,A1113,A1131,A1311,A3111) EQUIVALENCE (A35,A3331,A3313,A3133,A1333) EQUIVALENCE (A24,A2223,A2232,A2322,A3222) EQUIVALENCE (A34,A3332,A3323,A3233,A2333) EQUIVALENCE (A23,A2233,A3322) EQUIVALENCE (A13,A1133,A3311) EQUIVALENCE (A12,A1122,A2211) EQUIVALENCE (A44, A2323,A3232,A2332,A3223) EQUIVALENCE (A55, A1313,A3131,A1331,A3113) EQUIVALENCE (A66, A1212,A2121,A1221,A2112) EQUIVALENCE (A14,A1123,A1132,A2311,A3211) EQUIVALENCE (A25,A2213,A2231,A1322,A3122) EQUIVALENCE (A36,A3312,A3321,A1233,A2133) EQUIVALENCE (A56,A1321,A3112,A2113,A1231,A1213,A2131,A1312,A3121) EQUIVALENCE (A46,A2312,A3221,A1223,A2132,A2123,A1232,A2321,A3212) EQUIVALENCE (A45,A3213,A2331,A1332,A3123,A3132,A1323,A3231,A2313) REAL GAMMA(6),E11,E21,E31,E12,E22,E32,E13,E23,E33 C EQUIVALENCE (GAMMA(1),G1),(GAMMA(3),G2),(GAMMA(6),G3) C EQUIVALENCE (EE(1),E11),(EE(4),E12),(EE(7),E13) C EQUIVALENCE (EE(2),E21),(EE(5),E22),(EE(8),E23) C EQUIVALENCE (EE(3),E31),(EE(6),E32),(EE(9),E33) REAL A111,A112,A121,A122,A113,A123,A131,A132,A133 REAL A211,A212,A221,A222,A213,A223,A231,A232,A233 REAL A311,A312,A322,A313,A321,A323,A331,A332,A333 REAL AUX C C GAMMA,G1,G2,G3...Christoffel matrix, later its eigenvalues. C (E11,E12,E13) C EE=(E21,E22,E23)... Eigenvectors of the christoffel matrix. C (E31,E32,E33) C A111,A211,A311,A112,A212,A312,A122,A222,A322,A113,A213,A313,A123, C A223,A323,A133,A233,A333... A(I,J,K,L)*P(L) summed over L. C A11,A21,A31,A12,A22,A32,A13,A23,A33 ... Aijk*Ek C C....................................................................... C A11=B11 A22=B22 A33=B33 A16=B16 A26=B26 A15=B15 A35=B35 A24=B24 A34=B34 A23=B23 A13=B13 A12=B12 A44=B44 A55=B55 A66=B66 A14=B14 A25=B25 A36=B36 A56=B56 A46=B46 A45=B45 C Christoffel matrix: A111=A1111*P1+A1112*P2+A1113*P3 A112=A1121*P1+A1122*P2+A1123*P3 A121=A1211*P1+A1212*P2+A1213*P3 A122=A1221*P1+A1222*P2+A1223*P3 A113=A1131*P1+A1132*P2+A1133*P3 A123=A1231*P1+A1232*P2+A1233*P3 A131=A1311*P1+A1312*P2+A1313*P3 A132=A1321*P1+A1322*P2+A1323*P3 A133=A1331*P1+A1332*P2+A1333*P3 A211=A2111*P1+A2112*P2+A2113*P3 A212=A2121*P1+A2122*P2+A2123*P3 A221=A2211*P1+A2212*P2+A2213*P3 A222=A2221*P1+A2222*P2+A2223*P3 A213=A2131*P1+A2132*P2+A2133*P3 A223=A2231*P1+A2232*P2+A2233*P3 A231=A2311*P1+A2312*P2+A2313*P3 A232=A2321*P1+A2322*P2+A2323*P3 A233=A2331*P1+A2332*P2+A2333*P3 A311=A3111*P1+A3112*P2+A3113*P3 A312=A3121*P1+A3122*P2+A3123*P3 A322=A3221*P1+A3222*P2+A3223*P3 A313=A3131*P1+A3132*P2+A3133*P3 A321=A3211*P1+A3212*P2+A3213*P3 A323=A3231*P1+A3232*P2+A3233*P3 A331=A3311*P1+A3312*P2+A3313*P3 A332=A3321*P1+A3322*P2+A3323*P3 A333=A3331*P1+A3332*P2+A3333*P3 GAMMA(1)=P1*A111+P2*A211+P3*A311 GAMMA(2)=P1*A112+P2*A212+P3*A312 GAMMA(3)=P1*A122+P2*A222+P3*A322 GAMMA(4)=P1*A113+P2*A213+P3*A313 GAMMA(5)=P1*A123+P2*A223+P3*A323 GAMMA(6)=P1*A133+P2*A233+P3*A333 C C Selecting eigenvalue and eigenvector of the Christoffel matrix: CALL EIGEN(GAMMA,EE,3,0) G1=GAMMA(1) G2=GAMMA(3) G3=GAMMA(6) E11=EE(1) E21=EE(2) E31=EE(3) E12=EE(4) E22=EE(5) E32=EE(6) E13=EE(7) E23=EE(8) E33=EE(9) IF (G3.LT.0.) THEN C WAN-11 CALL ERROR('WAN-11: Negative eigenvalue of Christoffel matrix') C This error should not appear. END IF AUX=E11*E11+E21*E21+E31*E31 IF (ABS(AUX-1.).GT.0.000001) THEN C WAN-12 CALL ERROR('WAN-12: Eigenvector is not normalized') C This error should not appear. ENDIF AUX=E12*E12+E22*E22+E32*E32 IF (ABS(AUX-1.).GT.0.000001) THEN C WAN-13 CALL ERROR('WAN-13: Eigenvector is not normalized') C This error should not appear. ENDIF AUX=E13*E13+E23*E23+E33*E33 IF (ABS(AUX-1.).GT.0.000001) THEN C WAN-14 CALL ERROR('WAN-14: Eigenvector is not normalized') C This error should not appear. ENDIF C C Computation of derivatives dx/dt: A11= A111*E11+A112*E21+A113*E31 A21= A211*E11+A212*E21+A213*E31 A31= A311*E11+A312*E21+A313*E31 A12= A121*E11+A122*E21+A123*E31 A22= A221*E11+A222*E21+A223*E31 A32= A321*E11+A322*E21+A323*E31 A13= A131*E11+A132*E21+A133*E31 A23= A231*E11+A232*E21+A233*E31 A33= A331*E11+A332*E21+A333*E31 DER(1)=A11*E11+ A12*E21+ A13*E31 DER(2)=A21*E11+ A22*E21+ A23*E31 DER(3)=A31*E11+ A32*E21+ A33*E31 A11= A111*E12+A112*E22+A113*E32 A21= A211*E12+A212*E22+A213*E32 A31= A311*E12+A312*E22+A313*E32 A12= A121*E12+A122*E22+A123*E32 A22= A221*E12+A222*E22+A223*E32 A32= A321*E12+A322*E22+A323*E32 A13= A131*E12+A132*E22+A133*E32 A23= A231*E12+A232*E22+A233*E32 A33= A331*E12+A332*E22+A333*E32 DER(4)=A11*E12+ A12*E22+ A13*E32 DER(5)=A21*E12+ A22*E22+ A23*E32 DER(6)=A31*E12+ A32*E22+ A33*E32 A11= A111*E13+A112*E23+A113*E33 A21= A211*E13+A212*E23+A213*E33 A31= A311*E13+A312*E23+A313*E33 A12= A121*E13+A122*E23+A123*E33 A22= A221*E13+A222*E23+A223*E33 A32= A321*E13+A322*E23+A323*E33 A13= A131*E13+A132*E23+A133*E33 A23= A231*E13+A232*E23+A233*E33 A33= A331*E13+A332*E23+A333*E33 DER(7)=A11*E13+ A12*E23+ A13*E33 DER(8)=A21*E13+ A22*E23+ A23*E33 DER(9)=A31*E13+ A32*E23+ A33*E33 C RETURN END C C======================================================================= C SUBROUTINE WAMAT(A,B,C,D) C C---------------------------------------------------------------------- C Subroutine to compute the product of two 2x2 complex matrices. C The second matrix (C+iD) is destroyed in the computation. REAL A(2,2),B(2,2),C(2,2),D(2,2) C Input: C A,B,C,D ... Real and imaginary parts of the two matrices. C Output: C C,D ... Real and imaginary parts of resulting matrix. C C....................................................................... C Auxiliary storage locations: REAL CR11,CR21,CR12,CR22,CI11,CI21,CI12,CI22 C....................................................................... CR11=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1) CR21=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1) CR12=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2) CR22=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2) C CI11=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1) CI21=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1) CI12=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2) CI22=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2) C C(1,1)=CR11 C(2,1)=CR21 C(1,2)=CR12 C(2,2)=CR22 D(1,1)=CI11 D(2,1)=CI21 D(1,2)=CI12 D(2,2)=CI22 C RETURN END C C======================================================================= C REAL FUNCTION WACHAN(G,H) C C---------------------------------------------------------------------- C Subroutine to compute the smallest angle between the two-dimensional C vector G and vectors U,V,-U,-V. C The subroutine reorganizes the vectors U and V in such way, C that the pair U,V is equal to pair G,H rotated with angle WACHAN. C The real numbers R1 (associated with U) and R2 (associated with V) C are reorganized in the same way. C REAL G(4),H(7) C G(1:4)=G1,G2,H1,H2 C H(1:7)=R1,R2,R3,U1,U2,V1,V2 C Input: C G1,G2,H1,H2 ... A pair of two-dimensional orthonormal vectors. C U1,U2,V1,V2 ... A pair of two-dimensional orthonormal vectors. C R1,R2 ... Real numbers associeted to U and V. C Output: C WACHAN ... The smallest one from the angles between vector G C and vectors U,V,-U,-V. C U1,U2,V1,V2 ... Selection from U,V,-U,-V in such way, that C the pair U,V is equal to pair G,H rotated C with angle WACHAN. C R1,R2 ... Real numbers associeted to U and V. C C....................................................................... C C Auxiliary storage locations: REAL SP1,SP2,SP3,SP4,A1,A2,AUX REAL G1,G2,H1,H2,U1,U2,V1,V2,R1,R2,R3 C C----------------------------------------------------------------------- G1=G(1) G2=G(2) H1=G(3) H2=G(4) R3=H(1) R2=H(2) R1=H(3) U1=H(4) U2=H(5) V1=H(6) V2=H(7) C SP1=ABS(1 -( G1*U1+G2*U2)) SP2=ABS(1 -(-G1*U1-G2*U2)) SP3=ABS(1 -( G1*V1+G2*V2)) SP4=ABS(1 -(-G1*V1-G2*V2)) AUX=AMIN1(SP1,SP2,SP3,SP4) IF (AUX.EQ.SP1) THEN C No action. ELSEIF (AUX.EQ.SP2) THEN U1=-U1 U2=-U2 ELSEIF (AUX.EQ.SP3) THEN A1=U1 A2=U2 U1=V1 U2=V2 V1=A1 V2=A2 AUX=R1 R1=R2 R2=AUX ELSEIF (AUX.EQ.SP4) THEN A1=U1 A2=U2 U1=-V1 U2=-V2 V1=A1 V2=A2 AUX=R1 R1=R2 R2=AUX ENDIF SP1=ABS(1 - ( H1*V1+H2*V2)) SP2=ABS(1 - (-H1*V1-H2*V2)) AUX=AMIN1(SP1,SP2) IF (AUX.EQ.SP1) THEN C No action. ELSEIF (AUX.EQ.SP2) THEN V1=-V1 V2=-V2 ENDIF WACHAN=ASIN(0.5*((G1+U1)*(H1-V1)+(G2+U2)*(H2-V2))) H(1)=R3 H(2)=R2 H(3)=R1 H(4)=U1 H(5)=U2 H(6)=V1 H(7)=V2 RETURN END C C C======================================================================= C SUBROUTINE WAPROJ(P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3, * GE1,GE2,HE1,HE2) C C---------------------------------------------------------------------- C Subroutine to project vectors G and H to the plane defined by vector E C and vector PxE. C REAL P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3,GE1,GE2,HE1,HE2 C Input: C P1,P2,P3,E1,E2,E3 ... Vectors defining the plane. C G1,G2,G3,H1,H2,H3 ... Vectors to be projected. C Output: C GE1,GE2,HE1,HE2 ... Projected vectors. C C....................................................................... C C Auxiliary storage locations: REAL F1,F2,F3,AUX C C----------------------------------------------------------------------- C Second vector defining the plane: F1=P2*E3-E2*P3 F2=P3*E1-E3*P1 F3=P1*E2-E1*P2 AUX=SQRT(F1*F1+F2*F2+F3*F3) F1=F1/AUX F2=F2/AUX F3=F3/AUX C Projecting vectors G and H to the plane defined by E and F: GE1=E1*G1+E2*G2+E3*G3 GE2=F1*G1+F2*G2+F3*G3 HE1=E1*H1+E2*H2+E3*H3 HE2=F1*H1+F2*H2+F3*H3 C RETURN END C======================================================================= C INCLUDE 'eigen.for' C eigen.for INCLUDE 'means.for' C means.for C C======================================================================= C