C
C Program CRTCART converting the rays calculated by program CRT from C curvilinear to Cartesian coordinates. C C Version: 5.80 C Date: 2003, November 23 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: klimes@seis.karlov.mff.cuni.cz C C....................................................................... C C C Description of data files: C C Main input data read from external interactive device (*): C The data consist of character strings, read by list directed (free C format) input. The strings have thus to be enclosed in C apostrophes. The interactive * external unit may be redirected to C the file containing the data. C (1) 'MODEL','RAYALL','INIALL','RAY2P','INI2P',/ C 'MODEL'... Formatted file with the input data for the model. C Only integer KOORS specifying the type of the coordinate C system is read from data file MODEL. C 'RAYALL'... Input file 'CRT-R' with the quantities stored along C rays, see C C.R.T.5.5.1, C or input file 'CRT-S' with the quantities stored C at the specified surfaces, see C C.R.T.5.5.2. C 'INIALL'... Input file 'CRT-I' with the quantities at the initial C points of rays, corresponding to file RAYALL, see C C.R.T.6.1. C 'RAY2P'... Output file 'CRT-R' or 'CRT-S', with the quantities C transformed into Cartesian coordinates. C 'INI2P'... Output file 'CRT-I' with the quantities at the initial C points of rays, corresponding to file RAY2P, see C C.R.T.6.1. C Default: 'MODEL'='model.dat', 'RAYALL'='r01.out', 'INIALL'='r01i.out', C 'RAY2P'='ray2p.out', 'INI2P'='ini2p.out'. C C Formatted file MODEL: C See the description within source code file 'model.for'. C Description of data MODEL C C Unformatted files RAYALL and RAY2P: C See the description within source code file 'writ.for'. C Description of files 'CRT-R' C Description of files 'CRT-S' C C Unformatted files INIALL and INI2P: C See the description within source code file 'writ.for'. C Description of files 'CRT-I' C C----------------------------------------------------------------------- C C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C pointc.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL AP00 C AP00... File 'ap.for'. C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER LU1,LU2,LU3,LU4,I PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4) CHARACTER*80 FILE0,FILE1,FILE2,FILE3,FILE4 CHARACTER*1 TEXTM C C....................................................................... C C Opening input and output files: FILE0='model.dat' FILE1='r01.out' FILE2='r01i.out' FILE3='ray2p.out' FILE4='ini2p.out' WRITE(*,'(2A)') * ' Enter 5 filenames', * ' (model.dat, r01.out, r01i.out, ray2p.out, ini2p.out): ' READ(*,*) FILE0,FILE1,FILE2,FILE3,FILE4 OPEN(LU1,FILE=FILE0,STATUS='OLD') READ(LU1,*) TEXTM I=0 READ(LU1,*) I CALL METR1(I) CLOSE(LU1) OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED') OPEN(LU4,FILE=FILE4,FORM='UNFORMATTED') C C Loop for the points of rays 10 CONTINUE C Reading the results of the complete ray tracing CALL AP00(0,LU1,LU2) IF(IWAVE.LT.1)THEN C End of rays GO TO 80 END IF C Conversion to Cartesian coordinates CALL TOCART(Y,YL) WRITE(LU3) ISRC,IWAVE,IRAY,NY,ICB1,ISRF,KMAH, * X,UEBRAY,YL,(Y(I),I=1,NY) IF(IPT.LE.1)THEN C New ray - recording the initial point CALL TOCART(YI,YLI) WRITE(LU4) ISRC,-IWAVE,IRAY,ICB1I,IEND,ISHEET,IREC,YLI,YI END IF GO TO 10 C 80 CONTINUE CLOSE(LU1) CLOSE(LU2) CLOSE(LU3) CLOSE(LU4) STOP END C C======================================================================= C C C SUBROUTINE TOCART(Y,YL) REAL Y(11),YL(6) C C This subroutine transforms points, slowness vectors, R.C.C.S. basis C vector and velocity gradient from the model coordinates to Cartesian C coordinates. C C Input and output: C Y... Array containing basic quantities computed along the ray, C C.R.T.5.2.1. C Description of Y C YL... Array containing local quantities at the point of the ray C C.R.T.5.5.4. C Description of YL C C Subroutines and external functions required: EXTERNAL CARTES C C Date: 1998, February 28 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL CART(3),CART1,CART2,CART3 REAL PDER(9),PDER1,PDER2,PDER3,PDER4,PDER5,PDER6,PDER7,PDER8,PDER9 EQUIVALENCE (CART(1),CART1),(CART(2),CART2),(CART(3),CART3) EQUIVALENCE (PDER(1),PDER1),(PDER(2),PDER2),(PDER(3),PDER3) EQUIVALENCE (PDER(4),PDER4),(PDER(5),PDER5),(PDER(6),PDER6) EQUIVALENCE (PDER(7),PDER7),(PDER(8),PDER8),(PDER(9),PDER9) C C....................................................................... C C Conversion of model coordinates to Cartesian coordinates CART: CALL CARTES(Y(3),.TRUE. ,CART,PDER) C Transposed transformation matrix PDER of covariant vectors: CALL CARTES(Y(3),.FALSE.,CART,PDER) C C Transformations: Y( 3)=CART1 Y( 4)=CART2 Y( 5)=CART3 CART1=PDER1*Y( 6)+PDER2*Y( 7)+PDER3*Y( 8) CART2=PDER4*Y( 6)+PDER5*Y( 7)+PDER6*Y( 8) CART3=PDER7*Y( 6)+PDER8*Y( 7)+PDER9*Y( 8) Y( 6)=CART1 Y( 7)=CART2 Y( 8)=CART3 CART1=PDER1*Y( 9)+PDER2*Y(10)+PDER3*Y(11) CART2=PDER4*Y( 9)+PDER5*Y(10)+PDER6*Y(11) CART3=PDER7*Y( 9)+PDER8*Y(10)+PDER9*Y(11) Y( 9)=CART1 Y(10)=CART2 Y(11)=CART3 CART1=PDER1*YL(4)+PDER2*YL(5)+PDER3*YL(6) CART2=PDER4*YL(4)+PDER5*YL(5)+PDER6*YL(6) CART3=PDER7*YL(4)+PDER8*YL(5)+PDER9*YL(6) YL(4)=CART1 YL(5)=CART2 YL(6)=CART3 END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'metric.for' C metric.for INCLUDE 'ap.for' C ap.for C C======================================================================= C