C
C Program TCGREEN to compute propagator matrices U of equation (39) of C the paper by Klimes (1999). The propagator matrices are written C to the output formatted files in the form of file C GREEN. C C References: C Klimes L.(1999): 'Analytical one-way plane-wave solution in the C 1-D anisotropic "twisted crystal" model'. Report 8, pp.103-118, C Charles University, Prague. C C Bulant, P., & Klimes L.(1999): 'Comparison of ray methods C with the exact solution in the 1-D anisotropic C "twisted crystal" model'. Report 8, pp.119-126, C Charles University, Prague. C C Version: 5.80 C Date: 2004, June 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 http://sw3d.cz/staff/bulant.htm C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of a SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Data describing the twisted-crystal model: C SIN2TH=real ... Second power of sinus theta C (equation 15 of the paper by Klimes). C Used to specify the default values of EPSILON and V0. C Has no meaning if both EPSILON and V0 are specified. C Default: SIN2TH=0. C GAMMA=real ... Gamma of equation 15 of the paper by Klimes. C Used to specify the default values of EPSILON and V0. C Has no meaning if both EPSILON and V0 are specified. C Default: GAMMA=0. C TCK=real ... Uppercase K of equation 9 of the paper by Klimes. C Default: TCK=0. C EPSILON=real ... Epsilon of equation 9 of the paper by Klimes. C Default: EPSILON=-(GAMMA*SIN2TH)/(1.+GAMMA*SIN2TH) C A44=real ... a44 of equation 3 of the paper by Bulant & Klimes. C Used to specify the default value of V0, and to specify C the reference stress for output file TCSTRESS. C Default: A44=0. C V0=real ... v0 of equation 8 of the paper by Klimes. C Default: V0=SQRT(A44*(1.+GAMMA*SIN2TH)) C VREF2=real ... second power of reference velocity C (see equation 111 of the paper by Klimes). C Default: VREF2=V0*V0 C Data describing source and receiver: C REC='string'... If non-blank, the name of the file with the names C and coordinates of the receiver points. Only the first C point in the file is taken into account. Its name is C used in output files. If X3 is not defined, x3 coordinate C of the point is used instead. C Description of file C REC C Default: REC=' ' C SRC='string'... If non-blank, the name of the file with the name C and coordinates of the source point. The name is used in C the output files, the coordinates are not considered. C Description of file C SRC C Default: SRC=' ' C X3=real ... Parameter specifying the source-receiver C distance X3=X3receiver-X3source in the direction of the C X3 axis. If not specified, coordinates of the source and C receiver are taken from files SRC and REC (default). C Default: X3=X3receiver-X3source C Data describing the frequency domain: C DF=real ... Frequency step. C Default: DF=0. C OF=real ... The propagator matrix is calculated for NF C frequencies OF,OF+DF,OF+2*DF,...,OF+(NF-1)*DF. C Default: OF=0. C NF=integer ... Number of frequencies. C Default: NF=1 C Names of the output files: C The output files are not generated if the corresponding C name equals ' '. C TCGREENE='string'... Name of the output formatted file with the C exact propagator matrix U (equation 39) written in the C form of the file GREEN. C Default: TCGREENE='tcgreene.out' C TCGREENW='string'... Name of the output formatted file with the C coupling ray theory propagator matrix (equation 39 with C coefficients F0 to F3 from equation 106) written in the C form of the file GREEN. C Default: TCGREENW='tcgreenw.out' C TCGREENQ='string'... Name of the output formatted file with the C quasi-isotropic propagator matrix (equation 39 with C coefficients F0 to F3 from equation 111) written in the C form of the file GREEN. C Default: TCGREENQ='tcgreenq.out' C TCGREENA='string'... Name of the output formatted file with the C anisotropic propagator matrix (equation 39 with C coefficients F0 to F3 from equation 116) written in the C form of the file GREEN. C Default: TCGREENA='tcgreena.out' C TCGREENI='string'... Name of the output formatted file with the C isotropic propagator matrix (equation 39 with C coefficients F0 to F3 from equation 121) written in the C form of the file GREEN. C Default: TCGREENI='tcgreeni.out' C TCSTRESS='string'... Name of the output formatted file with the C stress components corresponding to the exact propagator C matrix at the initial point. The stress is divided by C "reference stress"=2*pi*F*SQRT(A44), C where F is frequency. C Default: TCSTRESS='tcstress.out' C Rotation of the model: C TCE1=real, TCE2=real, TCE3=real ... Three components of C the rotation vector describing the rotation of the model. C The model described in the papers referenced above may be C further rotated by specifying the rotation vector. C The model is rotated along the axis given by the rotation C vector, the angle of the rotation in radians is given by C the Cartesian length of the vector. C Only the exact solution may be calculated in the rotated C model, all the parameters TCGREENW(Q,A,I) must equal ' ' C if the rotation vector is specified. C Default: TCE1=TCE2=TCE3=0. (no rotation). C C Description of the output files TCGREENE(W,Q,A,I): C Formatted files containing the computed propagator matrix. C The files have the same format as the file GREEN generated by C program GREEN. Strings 'Src' and 'Rec' are written in place C of names of the source and the receiver, coordinates of the source C are [0.,0.,0.], coordinates of the receiver are [0.,0.,X3], C travel time is X3/V0, zeros are written in place of derivatives C of the travel time with respect to the coordinates of the receiver C and coordinates of the source. Components of the 2x2 propagator C matrix are written as corresponding amplitudes, zeros are written C instead of the remaining amplitudes. C Description of file C GREEN C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,RSEP1,RSEP3I,RSEP3T,RSEP3R,FORM2,LENGTH, * EQ39,TCGMAT,TCGTRA INTEGER LENGTH C ERROR ... File C error.for. C RSEP1,RSEP3I,RSEP3T,RSEP3R ... C File sep.for. C FORM2 ... File C forms.for. C LENGTH ... File C length.for. C EQ39,TCGMAT,TCGTRA ... This file. C C----------------------------------------------------------------------- C C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc REAL GREEN(MRAM) EQUIVALENCE (GREEN,RAM) C C Input and output data files: CHARACTER*80 FSEP,FOUTE,FOUTW,FOUTQ,FOUTA,FOUTI,FOUTS, * FILREC,FILSRC,RECNAM,SRCNAM CHARACTER*260 FORMAT INTEGER LU PARAMETER (LU=1) REAL PI PARAMETER (PI=3.1415926) C C Auxiliary storage locations: INTEGER I1,I2,I,NF,NGREEN REAL EPS,S2TH,GAMMA,OF,DF,F,VK,K0,A44,V0,X3,K02,VK2, * E1,E2,E3,AE,AE2,COSAE,SINAE,AUX1,AUX2, * K02VK2,EPS2,JEPS2,SJEPS2,SQJMEP,SQJPEP, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,RECIT, * REJM,IMJM,AUX,REAUX,IMAUX,REPZ,IMPZ,REDZ,IMDZ, * REF12,IMF12,REF02,IMF02,AF02,FIF02,AF0,FIF0, * REU11,IMU11,REU21,IMU21,REU31,IMU31, * REU12,IMU12,REU22,IMU22,REU32,IMU32, * REU13,IMU13,REU23,IMU23,REU33,IMU33,REU(3,3),IMU(3,3), * REE11,IME11,REE21,IME21,REE31,IME31, * REE12,IME12,REE22,IME22,REE32,IME32, * REE13,IME13,REE23,IME23,REE33,IME33,REE(3,3),IME(3,3), * CFIF02,SFIF02,TT,VR,VR2,R1,R2,R3, * REA(3,3),IMA(3,3),REET(3,3),IMET(3,3) LOGICAL LROT EQUIVALENCE (REU11,REU(1,1)),(IMU11,IMU(1,1)),(REU21,REU(2,1)), * (IMU21,IMU(2,1)),(REU31,REU(3,1)),(IMU31,IMU(3,1)), * (REU12,REU(1,2)),(IMU12,IMU(1,2)),(REU22,REU(2,2)), * (IMU22,IMU(2,2)),(REU32,REU(3,2)),(IMU32,IMU(3,2)), * (REU13,REU(1,3)),(IMU13,IMU(1,3)),(REU23,REU(2,3)), * (IMU23,IMU(2,3)),(REU33,REU(3,3)),(IMU33,IMU(3,3)), * (REE11,REE(1,1)),(IME11,IME(1,1)),(REE21,REE(2,1)), * (IME21,IME(2,1)),(REE31,REE(3,1)),(IME31,IME(3,1)), * (REE12,REE(1,2)),(IME12,IME(1,2)),(REE22,REE(2,2)), * (IME22,IME(2,2)),(REE32,REE(3,2)),(IME32,IME(3,2)), * (REE13,REE(1,3)),(IME13,IME(1,3)),(REE23,REE(2,3)), * (IME23,IME(2,3)),(REE33,REE(3,3)),(IME33,IME(3,3)) DATA REE,IME/18*0./ C C----------------------------------------------------------------------- C C Main input data: FSEP=' ' WRITE(*,'(A)') '+TCGREEN: Enter input filename: ' READ(*,*) FSEP IF (FSEP.EQ.' ') THEN C TCGREEN-01 CALL ERROR('TCGREEN-01: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C WRITE(*,'(A)') '+TCGREEN: Working... ' C C Reading all the data from file FSEP to the memory C (SEP parameter file form): CALL RSEP1(LU,FSEP) C C Reading the names of the output files: CALL RSEP3T('TCGREENE',FOUTE,'tcgreene.out') CALL RSEP3T('TCGREENW',FOUTW,'tcgreenw.out') CALL RSEP3T('TCGREENQ',FOUTQ,'tcgreenq.out') CALL RSEP3T('TCGREENA',FOUTA,'tcgreena.out') CALL RSEP3T('TCGREENI',FOUTI,'tcgreeni.out') CALL RSEP3T('TCSTRESS',FOUTS,'tcstress.out') C C Recalling the data: C (arguments: Name of value in input data, Variable, Default): CALL RSEP3R('EPSILON',EPS,999999.) CALL RSEP3R('SIN2TH',S2TH,0.) CALL RSEP3R('GAMMA',GAMMA,0.) CALL RSEP3R('OF',OF,0.) CALL RSEP3R('DF',DF,1.) CALL RSEP3I('NF',NF,1) CALL RSEP3R('TCK',VK,0.) CALL RSEP3R('V0',V0,999999.) CALL RSEP3R('A44',A44,0.) CALL RSEP3R('TCE1',E1,0.) CALL RSEP3R('TCE2',E2,0.) CALL RSEP3R('TCE3',E3,0.) AE2=E1*E1+E2*E2+E3*E3 IF (AE2.EQ.0.) THEN LROT=.FALSE. ELSE LROT=.TRUE. AE=SQRT(AE2) COSAE=COS(AE) AUX1=(1-COSAE)/AE2 SINAE=SIN(AE) AUX2=SINAE/AE REE11=COSAE+E1*E1*AUX1 REE21= E2*E1*AUX1 -E3*AUX2 REE31= E3*E1*AUX1 +E2*AUX2 REE12= E1*E2*AUX1 +E3*AUX2 REE22=COSAE+E2*E2*AUX1 REE32= E3*E2*AUX1-E1*AUX2 REE13= E1*E3*AUX1 -E2*AUX2 REE23= E2*E3*AUX1+E1*AUX2 REE33=COSAE+E3*E3*AUX1 CALL TCGTRA(REE,IME,REET,IMET) IF ((FOUTW.NE.' ').OR.(FOUTQ.NE.' ').OR.(FOUTA.NE.' ').OR. * (FOUTI.NE.' ')) THEN C TCGREEN-06 CALL WARN('TCGREEN-06: Only file TCGREENE may be calculated') C If the model is to be rotated, only the exact solution may be C calculated. All other output files will not be generated. FOUTW=' ' FOUTQ=' ' FOUTA=' ' FOUTI=' ' ENDIF ENDIF C Name of the source: CALL RSEP3T('SRC',FILSRC,' ') IF (FILSRC.NE.' ') THEN OPEN(LU,FILE=FILSRC,STATUS='OLD') READ(LU,*) (SRCNAM,I=1,20) R3=0. READ(LU,*) SRCNAM,R1,R2,R3 CLOSE(LU) ENDIF C Name of the receiver: CALL RSEP3T('REC',FILREC,' ') IF (FILREC.NE.' ') THEN OPEN(LU,FILE=FILREC,STATUS='OLD') READ(LU,*) (RECNAM,I=1,20) X3=0. READ(LU,*) RECNAM,R1,R2,X3 CLOSE(LU) ENDIF R3=X3-R3 CALL RSEP3R('X3',X3,R3) IF (EPS.EQ.999999.) THEN AUX=GAMMA*S2TH EPS=-AUX/(1.+AUX) ENDIF EPS2=EPS*EPS JEPS2=1.-EPS2 SJEPS2=SQRT(JEPS2) SQJMEP=SQRT(1.-EPS) SQJPEP=SQRT(1.+EPS) IF (V0.EQ.999999.) THEN AUX=GAMMA*S2TH V0=SQRT(A44*(1+AUX)) ENDIF IF (V0.EQ.0.) THEN C TCGREEN-03 CALL ERROR('TCGREEN-03: Wrong value of V0') C V0 must not equal zero. ENDIF IF (VK.EQ.0.) THEN C TCGREEN-04 CALL ERROR('TCGREEN-04: Wrong value of TCK') C TCK must not equal zero. ENDIF IF (FOUTS.NE.' ') THEN IF (A44.LE.0.) THEN C TCGREEN-05 CALL ERROR('TCGREEN-05: Wrong value of A44') C A44 must be positive if TCSTRESS is not blank. ENDIF ENDIF DO 10, I2=1,14 GREEN(I2)=0. 10 CONTINUE TT=X3/V0 GREEN(1)=TT GREEN(5)=X3 NGREEN=14+18*NF IF (NGREEN.GT.MRAM) THEN C TCGREEN-02 CALL ERROR('TCGREEN-02: Small array RAM.') C Parameter MRAM of the file ram.inc should be greater than C 18 times number of frequencies plus 14. C File ram.inc. ENDIF C C C Exact solution: IF (FOUTE.NE.' '.OR.FOUTS.NE.' ') THEN IF (FOUTS.NE.' ') THEN OPEN(LU,FILE=FOUTS) END IF I2=15 DO 100, I1=0,NF-1 F=OF+I1*DF K0=2.*PI*F/V0 C Eq 63: VK2=VK*VK K02=K0*K0 K02VK2=K02-VK2 AUX=1.-EPS2*((VK2/K02VK2)**2) REJM=K02VK2*SQRT(JEPS2) IF (AUX.GE.0.) THEN IMJM=0. REJM=REJM*(SJEPS2+SQRT(AUX)) ELSE IMJM=REJM*SQRT(-AUX) REJM=REJM*SJEPS2 ENDIF RECIT=EPS*VK*K02 IF (IMJM.EQ.0.) THEN REF1=RECIT/REJM IMF1=0. ELSE AUX=REJM*REJM-IMJM*IMJM REF1=RECIT*REJM/AUX IMF1=-RECIT*IMJM/AUX ENDIF C Eq 52: REF12=REF1*REF1-IMF1*IMF1 IMF12=2.*REF1*IMF1 REPZ=K02/JEPS2+REF12 IMPZ=IMF12 REAUX=1.+REF12*JEPS2/VK2 IMAUX= IMF12*JEPS2/VK2 AUX=REAUX*REAUX+IMAUX*IMAUX REDZ= REAUX/AUX IMDZ=-IMAUX/AUX REF02=REPZ*REDZ-IMPZ*IMDZ IMF02=REPZ*IMDZ+IMPZ*REDZ AF02=SQRT(REF02*REF02+IMF02*IMF02) IF (AF02.EQ.0.) THEN REF0=0. IMF0=0. ELSE CFIF02=REF02/AF02 SFIF02=IMF02/AF02 IF (CFIF02.GE.0.) THEN FIF02=ASIN(SFIF02) ELSEIF(SFIF02.GE.0.) THEN FIF02=ACOS(CFIF02) ELSE FIF02=-ACOS(CFIF02) ENDIF AF0=SQRT(AF02) FIF0=FIF02/2. REF0=AF0*COS(FIF0) IMF0=AF0*SIN(FIF0) IF (REF0.LT.0.) THEN REF0=-REF0 IMF0=-IMF0 ENDIF IF (IMF0.LT.0.) THEN IMF0=-IMF0 IMF1=-IMF1 ENDIF ENDIF C Eq 49: REAUX=-EPS+REF1*JEPS2/VK IMAUX= IMF1*JEPS2/VK REF3=REF0*REAUX-IMF0*IMAUX IMF3=IMF0*REAUX+REF0*IMAUX C Eq 47: REF2=VK+EPS*REF1 IMF2= EPS*IMF1 C C Eq 32,33,39: CALL EQ39(VK,X3,TT,F, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3, * REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22) C Remaining components of the 3x3 matrix: REU13=0. IMU13=0. REU23=0. IMU23=0. REU31=0. IMU31=0. REU32=0. IMU32=0. REU33=0. IMU33=0. C C Rotation (multiplication E x U x ET): IF (LROT) THEN CALL TCGMAT(REE,IME,REU,IMU,REA,IMA) CALL TCGMAT(REA,IMA,REET,IMET,REU,IMU) ENDIF C GREEN(I2 )=REU11*1000000. GREEN(I2+ 1)=IMU11*1000000. GREEN(I2+ 2)=REU21*1000000. GREEN(I2+ 3)=IMU21*1000000. GREEN(I2+ 4)=REU31*1000000. GREEN(I2+ 5)=IMU31*1000000. GREEN(I2+ 6)=REU12*1000000. GREEN(I2+ 7)=IMU12*1000000. GREEN(I2+ 8)=REU22*1000000. GREEN(I2+ 9)=IMU22*1000000. GREEN(I2+10)=REU32*1000000. GREEN(I2+11)=IMU32*1000000. GREEN(I2+12)=REU13*1000000. GREEN(I2+13)=IMU13*1000000. GREEN(I2+14)=REU23*1000000. GREEN(I2+15)=IMU23*1000000. GREEN(I2+16)=REU33*1000000. GREEN(I2+17)=IMU33*1000000. I2=I2+18 C IF (FOUTS.NE.' ') THEN C Writing the density normalized initial stress, divided by C "reference stress" = 2.*PI*F*SQRT(A44) = K0*V0*SQRT(A44): IF(F.EQ.0.) THEN S11R= 0. S11I= SQRT(1.-EPS**2)*V0/SQRT(A44) S12R= 0. S12I= 0. S22R= 0. S22I= SQRT(1.-EPS**2)*V0/SQRT(A44) ELSE S11R=-(1.+EPS)*(IMF0+IMF3)*V0/K0/SQRT(A44) S11I= (1.+EPS)*(REF0+REF3)*V0/K0/SQRT(A44) S12R=-(1.-EPS**2)*REF1 *V0/K0/SQRT(A44) S12I=-(1.-EPS**2)*IMF1 *V0/K0/SQRT(A44) S22R=-(1.-EPS)*(IMF0-IMF3)*V0/K0/SQRT(A44) S22I= (1.-EPS)*(REF0-REF3)*V0/K0/SQRT(A44) END IF WRITE(LU,'(7(F9.6,1X))') F,S11R,S11I,S12R,S12I,S22R,S22I ENDIF 100 CONTINUE C IF (FOUTS.NE.' ') THEN CLOSE(LU) ENDIF C IF (FOUTE.NE.' ') THEN C Writing: OPEN(LU,FILE=FOUTE) WRITE(LU,'(A)') ' /' FORMAT(1:4)='(6A,' CALL FORM2(14,GREEN,GREEN,FORMAT(5:260)) WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)), * ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''', * (' ',GREEN(I),I=1,14) DO 110 I2=15,NGREEN-18,18 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17) 110 CONTINUE I2=NGREEN-17 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /' WRITE(LU,'(A)') ' /' CLOSE(LU) ENDIF ENDIF C C C Coupling ray theory: IF (FOUTW.NE.' ') THEN OPEN(LU,FILE=FOUTW) WRITE(LU,'(A)') ' /' I2=15 DO 200, I1=0,NF-1 F=OF+I1*DF K0=2.*PI*F/V0 C Eq 106: REF0=K0*(SQJPEP+SQJMEP)/2./SJEPS2 IMF0=0. REF1=0. IMF1=0. REF2=VK IMF2=0. REF3=-K0*(SQJPEP-SQJMEP)/2./SJEPS2 IMF3=0. C Eq 32,33,39: CALL EQ39(VK,X3,TT,F, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3, * REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22) C GREEN(I2 )=REU11*1000000. GREEN(I2+ 1)=IMU11*1000000. GREEN(I2+ 2)=REU21*1000000. GREEN(I2+ 3)=IMU21*1000000. GREEN(I2+ 4)=0. GREEN(I2+ 5)=0. GREEN(I2+ 6)=REU12*1000000. GREEN(I2+ 7)=IMU12*1000000. GREEN(I2+ 8)=REU22*1000000. GREEN(I2+ 9)=IMU22*1000000. GREEN(I2+10)=0. GREEN(I2+11)=0. GREEN(I2+12)=0. GREEN(I2+13)=0. GREEN(I2+14)=0. GREEN(I2+15)=0. GREEN(I2+16)=0. GREEN(I2+17)=0. I2=I2+18 200 CONTINUE C Writing: FORMAT(1:4)='(6A,' CALL FORM2(14,GREEN,GREEN,FORMAT(5:260)) WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)), * ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''', * (' ',GREEN(I),I=1,14) DO 210 I2=15,NGREEN-18,18 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17) 210 CONTINUE I2=NGREEN-17 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /' WRITE(LU,'(A)') ' /' CLOSE(LU) ENDIF C C C Quasi-isotropic approach: IF (FOUTQ.NE.' ') THEN OPEN(LU,FILE=FOUTQ) WRITE(LU,'(A)') ' /' I2=15 CALL RSEP3R('VREF2',VR2,V0*V0) VR=SQRT(VR2) DO 300, I1=0,NF-1 F=OF+I1*DF K0=2.*PI*F/V0 C Eq 111: REF0=K0*V0/VR*(3./2.-1./2.*V0*V0/VR/VR) IMF0=0. REF1=0. IMF1=0. REF2=VK IMF2=0. REF3=-K0/2.*EPS*V0*V0*V0/VR/VR/VR IMF3=0. C Eq 32,33,39: CALL EQ39(VK,X3,TT,F, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3, * REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22) C GREEN(I2 )=REU11*1000000. GREEN(I2+ 1)=IMU11*1000000. GREEN(I2+ 2)=REU21*1000000. GREEN(I2+ 3)=IMU21*1000000. GREEN(I2+ 4)=0. GREEN(I2+ 5)=0. GREEN(I2+ 6)=REU12*1000000. GREEN(I2+ 7)=IMU12*1000000. GREEN(I2+ 8)=REU22*1000000. GREEN(I2+ 9)=IMU22*1000000. GREEN(I2+10)=0. GREEN(I2+11)=0. GREEN(I2+12)=0. GREEN(I2+13)=0. GREEN(I2+14)=0. GREEN(I2+15)=0. GREEN(I2+16)=0. GREEN(I2+17)=0. I2=I2+18 300 CONTINUE C Writing: FORMAT(1:4)='(6A,' CALL FORM2(14,GREEN,GREEN,FORMAT(5:260)) WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)), * ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''', * (' ',GREEN(I),I=1,14) DO 310 I2=15,NGREEN-18,18 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17) 310 CONTINUE I2=NGREEN-17 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /' WRITE(LU,'(A)') ' /' CLOSE(LU) ENDIF C C C Anisotropic ray theory: IF (FOUTA.NE.' ') THEN OPEN(LU,FILE=FOUTA) WRITE(LU,'(A)') ' /' I2=15 DO 400, I1=0,NF-1 F=OF+I1*DF K0=2.*PI*F/V0 C Eq 116: REF0=K0*(SQJPEP+SQJMEP)/2./SJEPS2 IMF0=0. REF1=0. IMF1=0. REF2=0. IMF2=0. REF3=-K0*(SQJPEP-SQJMEP)/2./SJEPS2 IMF3=0. C Eq 32,33,39: CALL EQ39(VK,X3,TT,F, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3, * REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22) C GREEN(I2 )=REU11*1000000. GREEN(I2+ 1)=IMU11*1000000. GREEN(I2+ 2)=REU21*1000000. GREEN(I2+ 3)=IMU21*1000000. GREEN(I2+ 4)=0. GREEN(I2+ 5)=0. GREEN(I2+ 6)=REU12*1000000. GREEN(I2+ 7)=IMU12*1000000. GREEN(I2+ 8)=REU22*1000000. GREEN(I2+ 9)=IMU22*1000000. GREEN(I2+10)=0. GREEN(I2+11)=0. GREEN(I2+12)=0. GREEN(I2+13)=0. GREEN(I2+14)=0. GREEN(I2+15)=0. GREEN(I2+16)=0. GREEN(I2+17)=0. I2=I2+18 400 CONTINUE C Writing: FORMAT(1:4)='(6A,' CALL FORM2(14,GREEN,GREEN,FORMAT(5:260)) WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)), * ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''', * (' ',GREEN(I),I=1,14) DO 410 I2=15,NGREEN-18,18 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17) 410 CONTINUE I2=NGREEN-17 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /' WRITE(LU,'(A)') ' /' CLOSE(LU) ENDIF C C C Isotropic ray theory: IF (FOUTI.NE.' ') THEN OPEN(LU,FILE=FOUTI) WRITE(LU,'(A)') ' /' I2=15 DO 500, I1=0,NF-1 F=OF+I1*DF K0=2.*PI*F/V0 C Eq 121: REF0=K0*(SQJPEP+SQJMEP)/2./SJEPS2 IMF0=0. REF1=0. IMF1=0. REF2=VK IMF2=0. REF3=0. IMF3=0. C Eq 32,33,39: CALL EQ39(VK,X3,TT,F, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3, * REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22) C GREEN(I2 )=REU11*1000000. GREEN(I2+ 1)=IMU11*1000000. GREEN(I2+ 2)=REU21*1000000. GREEN(I2+ 3)=IMU21*1000000. GREEN(I2+ 4)=0. GREEN(I2+ 5)=0. GREEN(I2+ 6)=REU12*1000000. GREEN(I2+ 7)=IMU12*1000000. GREEN(I2+ 8)=REU22*1000000. GREEN(I2+ 9)=IMU22*1000000. GREEN(I2+10)=0. GREEN(I2+11)=0. GREEN(I2+12)=0. GREEN(I2+13)=0. GREEN(I2+14)=0. GREEN(I2+15)=0. GREEN(I2+16)=0. GREEN(I2+17)=0. I2=I2+18 500 CONTINUE C Writing: FORMAT(1:4)='(6A,' CALL FORM2(14,GREEN,GREEN,FORMAT(5:260)) WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)), * ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''', * (' ',GREEN(I),I=1,14) DO 510 I2=15,NGREEN-18,18 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17) 510 CONTINUE I2=NGREEN-17 FORMAT(1:4)='(1A,' CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260)) WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /' WRITE(LU,'(A)') ' /' CLOSE(LU) ENDIF C WRITE(*,'(A)') '+TCGREEN: Done. ' STOP END C C======================================================================= C SUBROUTINE EQ39(VK,X3,TT,F, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3, * REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22) C REAL VK,X3,TT,F, * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3, * REU11,IMU11,REU12,IMU12,REU21,IMU21,REU22,IMU22 C REAL PI PARAMETER (PI=3.1415926) REAL AUX,REAUX,IMAUX,REF12,IMF12,REF22,IMF22,REF32,IMF32, * REFI,IMFI,AFI,FIFI, * REFI2,IMFI2,AFI2,FIFI2, * RVF11,IVF11,RVF12,IVF12,RVF21,IVF21,RVF22,IVF22, * REVF11,IMVF11,REVF12,IMVF12,REVF21,IMVF21,REVF22,IMVF22, * RA11,IA11,RA12,IA12,RA21,IA21,RA22,IA22, * RB11,IB11,RB12,IB12,RB21,IB21,RB22,IB22, * RC11,IC11,RC12,IC12,RC21,IC21,RC22,IC22, * RD11,ID11,RD12,ID12,RD21,ID21,RD22,ID22, * CFIFI2,SFIFI2 C C Eq 32: REF12=REF1*REF1-IMF1*IMF1 IMF12=2.*REF1*IMF1 REF22=REF2*REF2-IMF2*IMF2 IMF22=2.*REF2*IMF2 REF32=REF3*REF3-IMF3*IMF3 IMF32=2.*REF3*IMF3 REFI2=REF32+REF22-REF12 IMFI2=IMF32+IMF22-IMF12 AFI2=SQRT(REFI2*REFI2+IMFI2*IMFI2) IF (AFI2.EQ.0.) THEN REFI=0. IMFI=0. ELSE CFIFI2=REFI2/AFI2 SFIFI2=IMFI2/AFI2 IF (CFIFI2.GE.0.) THEN FIFI2=ASIN(SFIFI2) ELSEIF(SFIFI2.GE.0.) THEN FIFI2=ACOS(CFIFI2) ELSE FIFI2=-ACOS(CFIFI2) ENDIF AFI=SQRT(AFI2) FIFI=FIFI2/2. REFI=AFI*COS(FIFI) IMFI=AFI*SIN(FIFI) ENDIF C Eq 33: AUX=REFI*REFI+IMFI*IMFI IF (AUX.NE.0.) THEN REAUX= REFI/AUX IMAUX=-IMFI/AUX RVF11=REF3 IVF11=IMF3 RVF12=-IMF1+IMF2 IVF12= REF1-REF2 RVF21=-IMF1-IMF2 IVF21= REF1+REF2 RVF22=-REF3 IVF22=-IMF3 REVF11=RVF11*REAUX-IVF11*IMAUX IMVF11=RVF11*IMAUX+IVF11*REAUX REVF12=RVF12*REAUX-IVF12*IMAUX IMVF12=RVF12*IMAUX+IVF12*REAUX REVF21=RVF21*REAUX-IVF21*IMAUX IMVF21=RVF21*IMAUX+IVF21*REAUX REVF22=RVF22*REAUX-IVF22*IMAUX IMVF22=RVF22*IMAUX+IVF22*REAUX ELSE REVF11=0. IMVF11=0. REVF12=0. IMVF12=0. REVF21=0. IMVF21=0. REVF22=0. IMVF22=0. ENDIF C Eq 39: RA11=COS(VK*X3) RA12=-SIN(VK*X3) RA21=-RA12 RA22=RA11 AUX=SIN(REFI*X3) RB11=COS(REFI*X3)-IMVF11*AUX IB11= REVF11*AUX RB12= -IMVF12*AUX IB12= REVF12*AUX RB21= -IMVF21*AUX IB21= REVF21*AUX RB22=COS(REFI*X3)-IMVF22*AUX IB22= REVF22*AUX RC11=RA11*RB11+RA12*RB21 IC11=RA11*IB11+RA12*IB21 RC12=RA11*RB12+RA12*RB22 IC12=RA11*IB12+RA12*IB22 RC21=RA21*RB11+RA22*RB21 IC21=RA21*IB11+RA22*IB21 RC22=RA21*RB12+RA22*RB22 IC22=RA21*IB12+RA22*IB22 AUX=REF0*X3 - 2.*PI*F*TT REAUX=COS(AUX) IMAUX=SIN(AUX) RD11=REAUX*RC11-IMAUX*IC11 ID11=REAUX*IC11+IMAUX*RC11 RD12=REAUX*RC12-IMAUX*IC12 ID12=REAUX*IC12+IMAUX*RC12 RD21=REAUX*RC21-IMAUX*IC21 ID21=REAUX*IC21+IMAUX*RC21 RD22=REAUX*RC22-IMAUX*IC22 ID22=REAUX*IC22+IMAUX*RC22 AUX=SINH(IMFI*X3) RA11=COSH(IMFI*X3)-REVF11*AUX IA11= -IMVF11*AUX RA12= -REVF12*AUX IA12= -IMVF12*AUX RA21= -REVF21*AUX IA21= -IMVF21*AUX RA22=COSH(IMFI*X3)-REVF22*AUX IA22= -IMVF22*AUX IF (IMF0.NE.0.) THEN AUX=-IMF0*X3 AUX=EXP(-IMF0*X3) RB11=AUX*RA11 IB11=AUX*IA11 RB12=AUX*RA12 IB12=AUX*IA12 RB21=AUX*RA21 IB21=AUX*IA21 RB22=AUX*RA22 IB22=AUX*IA22 ELSE RB11=RA11 IB11=IA11 RB12=RA12 IB12=IA12 RB21=RA21 IB21=IA21 RB22=RA22 IB22=IA22 ENDIF REU11=RD11*RB11-ID11*IB11 + RD12*RB21-ID12*IB21 REU12=RD11*RB12-ID11*IB12 + RD12*RB22-ID12*IB22 REU21=RD21*RB11-ID21*IB11 + RD22*RB21-ID22*IB21 REU22=RD21*RB12-ID21*IB12 + RD22*RB22-ID22*IB22 IMU11=RD11*IB11+ID11*RB11 + RD12*IB21+ID12*RB21 IMU12=RD11*IB12+ID11*RB12 + RD12*IB22+ID12*RB22 IMU21=RD21*IB11+ID21*RB11 + RD22*IB21+ID22*RB21 IMU22=RD21*IB12+ID21*RB12 + RD22*IB22+ID22*RB22 RETURN END C C======================================================================= C SUBROUTINE TCGMAT(A,B,C,D,E,F) C C---------------------------------------------------------------------- C Subroutine to compute the product of two 3x3 complex matrices C (E+iF)=(A+iB)*(C+iD) REAL A(3,3),B(3,3),C(3,3),D(3,3),E(3,3),F(3,3) C Input: C A,B,C,D ... Real and imaginary parts of the two matrices. C Output: C E,F ... Real and imaginary parts of the resulting matrix. C C....................................................................... E(1,1)=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1)+ * A(1,3)*C(3,1)-B(1,3)*D(3,1) E(2,1)=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1)+ * A(2,3)*C(3,1)-B(2,3)*D(3,1) E(3,1)=A(3,1)*C(1,1)-B(3,1)*D(1,1)+A(3,2)*C(2,1)-B(3,2)*D(2,1)+ * A(3,3)*C(3,1)-B(3,3)*D(3,1) E(1,2)=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2)+ * A(1,3)*C(3,2)-B(1,3)*D(3,2) E(2,2)=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2)+ * A(2,3)*C(3,2)-B(2,3)*D(3,2) E(3,2)=A(3,1)*C(1,2)-B(3,1)*D(1,2)+A(3,2)*C(2,2)-B(3,2)*D(2,2)+ * A(3,3)*C(3,2)-B(3,3)*D(3,2) E(1,3)=A(1,1)*C(1,3)-B(1,1)*D(1,3)+A(1,2)*C(2,3)-B(1,2)*D(2,3)+ * A(1,3)*C(3,3)-B(1,3)*D(3,3) E(2,3)=A(2,1)*C(1,3)-B(2,1)*D(1,3)+A(2,2)*C(2,3)-B(2,2)*D(2,3)+ * A(2,3)*C(3,3)-B(2,3)*D(3,3) E(3,3)=A(3,1)*C(1,3)-B(3,1)*D(1,3)+A(3,2)*C(2,3)-B(3,2)*D(2,3)+ * A(3,3)*C(3,3)-B(3,3)*D(3,3) C F(1,1)=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1)+ * A(1,3)*D(3,1)+B(1,3)*C(3,1) F(2,1)=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1)+ * A(2,3)*D(3,1)+B(2,3)*C(3,1) F(3,1)=A(3,1)*D(1,1)+B(3,1)*C(1,1)+A(3,2)*D(2,1)+B(3,2)*C(2,1)+ * A(3,3)*D(3,1)+B(3,3)*C(3,1) F(1,2)=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2)+ * A(1,3)*D(3,2)+B(1,3)*C(3,2) F(2,2)=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2)+ * A(2,3)*D(3,2)+B(2,3)*C(3,2) F(3,2)=A(3,1)*D(1,2)+B(3,1)*C(1,2)+A(3,2)*D(2,2)+B(3,2)*C(2,2)+ * A(3,3)*D(3,2)+B(3,3)*C(3,2) F(1,3)=A(1,1)*D(1,3)+B(1,1)*C(1,3)+A(1,2)*D(2,3)+B(1,2)*C(2,3)+ * A(1,3)*D(3,3)+B(1,3)*C(3,3) F(2,3)=A(2,1)*D(1,3)+B(2,1)*C(1,3)+A(2,2)*D(2,3)+B(2,2)*C(2,3)+ * A(2,3)*D(3,3)+B(2,3)*C(3,3) F(3,3)=A(3,1)*D(1,3)+B(3,1)*C(1,3)+A(3,2)*D(2,3)+B(3,2)*C(2,3)+ * A(3,3)*D(3,3)+B(3,3)*C(3,3) C RETURN END C C======================================================================= C SUBROUTINE TCGTRA(A,B,C,D) C C---------------------------------------------------------------------- C Subroutine to compute the transpose matrix to the 3x3 complex matrix C (C+iD)=(A+iB)^T REAL A(3,3),B(3,3),C(3,3),D(3,3) C Input: C A,B ... Real and imaginary parts of the input matrix. C Output: C C,D ... Real and imaginary parts of the transposed matrix. C C....................................................................... C(1,1)=A(1,1) C(2,1)=A(1,2) C(3,1)=A(1,3) C(1,2)=A(2,1) C(2,2)=A(2,2) C(3,2)=A(2,3) C(1,3)=A(3,1) C(2,3)=A(3,2) C(3,3)=A(3,3) C D(1,1)=B(1,1) D(2,1)=B(1,2) D(3,1)=B(1,3) D(1,2)=B(2,1) D(2,2)=B(2,2) D(3,2)=B(2,3) D(1,3)=B(3,1) D(2,3)=B(3,2) D(3,3)=B(3,3) C RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C======================================================================= C