C
C Program WFSRF to generate wavefronts composed of polygons for display C purposes C C Version: 5.60 C Preliminary version tested only within history file len-wf.h. C len-wf.h C Date: 2001, October 24 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 C The post processing program to the CRT program performs linear C interpolation of wavefronts. Wavefronts are approximated by triangles C composed of points located at the rays of individual ray tubes. If C such a triangle would intersect a structural interface, the wavefront C is locally approximated by a triangle and a tetragon. Similarly in C more complicated situations. C C This version of the program generates a single wavefront specified by C travel time WFTIME during one invocation of WFSRF. Only the wavefront C corresponding to the first elementary wave stored in the CRT output C file CRT-T is generated. C C The interpolation is not performed in demarcation belts between C different ray histories, similarly as in program C mtt.for. C The maximum width of the belts is controlled C by input parameter AERR of the CRT C program. The typical width of the belts, measured in take-off angles, C is 0.0001 radians in order of magnitude. The belts may become wide C in areas of large geometrical spreading. The division of rays into C different histories is, by default, very detailed and is controlled C by input parameter PRM0(3) of the C CRT program. Such a behaviour is useful for two-point ray tracing but C has some awkward consequences on the interpolation. For example, the C demarcation belt between the rays incident at different sides of the C computational volume for ray tracing extends across the whole model, C causing an artificial gap within the continuous wavefronts. C If the ray history does not account for the termination of a ray at C different sides of the computational volume, the gaps corresponding C to the edges of the computational volume are removed, but the corners C along the edges are not filled with the ray cells any more, and C the wavefront may disappear at the corners. Then the C computational volume for ray tracing should sufficiently exceed the C extent of target volume. C Input parameter PRM0(3) of the CRT C program may thus considerably influence the results of the C interpolation. However, gaps due to the demarcation belts C corresponding to structural interfaces cannot be removed within C the present interpolation algorithm. 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 the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Names of input files related to ray tracing (output of CRT): C CRTOUT='string'... File with the names of the output files of C program CRT. For general description of file CRTOUT refer C to file writ.for. C A single set of CRT output files is read from file CRTOUT. C Only files CRT-R, C CRT-I and C CRT-T are read by WFSRF. C File CRT-R must contain all rays traced by CRT. C Default: CRTOUT=' ' which means 'CRT-R'='r01.out', C 'CRT-I'='r01i.out', 'CRT-T'='t01.out' C Wavefront travel time: C WFTIME=real ... Travel time corresponding to the wavefront to C be covered by polygons. C Default: WFTIME=1 C Names of the output files: C VRTX='string'... Name of the output file with vertices of the C polygons. C Form of file VRTX C Default: VRTX='vrtx.out' C PLGN='string'... Name of the output file describing the polygons. C If blank, the file is not generated. C Form of file PLGN C Default: PLGN='plgn.out' C PLGNS='string'... Name of the file describing the polygons in C terms of the names of the vertices. C Form of file PLGNS C Default: PLGNS=' ' (the file is not generated) C Parameters specifying the quantities to be written into the file VRTX. C TEXTP='string' ... First part of names of vertices. The second C part of the name of a vertex is formed by number giving C its position in the file VRTX. C Default: TEXTP='V' C COLUMN01 to COLUMN69, POWER01 to POWER69, IVALUE01 to IVALUE69: C IVALUEii=integer ... An integer value required for some special C values of COLUMNii. Reserved for future extensions. C POWERii=real ... Power of the quantity to be written in column ii. C COLUMNii='string' ... String which specifies the quantity to be C written to the column ii of the file VRTX. First six C columns usually contain coordinates of the vertices and C the normals. Column zero is reserved for names of the C vertices. Following strings are allowed: C ' ' (a space) ... Nothing is to be written to the column C ii and to all the following columns. C 'X1' ... First coordinate of the vertex. C 'X2' ... Second coordinate of the vertex. C 'X3' ... Third coordinate of the vertex. C 'NORM1'. First component of the normal to the C wavefront at the vertex. Slowness vectors are used C as the normals in the current version of WFSRF. C 'NORM2'. Second component of the normal. C 'NORM3'. Third component of the normal. C 'ICB' .. Index of complex block in which the vertex C is located. C 'MTT' .. Travel time of the wavefront. C 'MHIST'. Ray history in the vertex. C 'MPQ11'. Component 1,1 of the 4x4 ray propagator matrix. C 'MPQ21'. Component 2,1 of the 4x4 ray propagator matrix. C 'MPQ31'. Component 3,1 of the 4x4 ray propagator matrix. C 'MPQ41'. Component 4,1 of the 4x4 ray propagator matrix. C 'MPQ12'. Component 1,2 of the 4x4 ray propagator matrix. C 'MPQ22'. Component 2,2 of the 4x4 ray propagator matrix. C 'MPQ32'. Component 3,2 of the 4x4 ray propagator matrix. C 'MPQ42'. Component 4,2 of the 4x4 ray propagator matrix. C 'MPQ13'. Component 1,3 of the 4x4 ray propagator matrix. C 'MPQ23'. Component 2,3 of the 4x4 ray propagator matrix. C 'MPQ33'. Component 3,3 of the 4x4 ray propagator matrix. C 'MPQ43'. Component 4,3 of the 4x4 ray propagator matrix. C 'MPQ14'. Component 1,4 of the 4x4 ray propagator matrix. C 'MPQ24'. Component 2,4 of the 4x4 ray propagator matrix. C 'MPQ34'. Component 3,4 of the 4x4 ray propagator matrix. C 'MPQ44'. Component 4,4 of the 4x4 ray propagator matrix. C 'GBW11'. Sum of squares of Gaussian beam widths C corresponding to GBY11. The sum of squares of Gaussian C beam widths is defined here as the trace of the C inverse imaginary part of the matrix of second C derivatives of the travel time of the Gaussian beam. C The sum of squares of Gaussian beam widths depends on C parameters GBEij, GBR11, GBR12, GBR22, GBY11 and GBY22 C described below. It may be expressed as sum W11+W22, C where W11 is dependent on GBY11 but independent C on GBY22, and W22 is dependent on GBY22 but C independent on GBY11. Values of W11 are written to C the column given by COLUMNii='GBW11'. C Note that positive GBY11 must be specified for C this option. C 'GBW1'.. Values of SQRT(W11). See the description of C COLUMNii='GBW11'. Positive GBY11 must be specified. C 'GBW22'.. Values of W22. See the description of C COLUMNii='GBW11'. Positive GBY22 must be specified. C 'GBW2'.. Values of SQRT(W22). See the description of C COLUMNii='GBW11'. Positive GBY22 must be specified. C 'AMP' .. Amplitude of the Green function in the vertex. C All strings may be entered either in uppercase or in C lowercase letters. C Defaults: COLUMN01='X1', COLUMN02='X2', COLUMN03='X3', C COLUMN04='NORM1', COLUMN05='NORM2', COLUMN06='NORM3', C COLUMN07 to COLUMN69=' ', C POWER01 to POWER69=1, IVALUE01 to IVALUE69=1. C Numerical parameters more precisely specifying the above quantities: C Refer to the list in program mtt.for. C Numerical parameters specifying optional ray-parameter grid: C Refer to the list in program mtt.for. C C----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL WFEROR,WFREAD,WFERAS,WFTUBE,WFINTP, * WFQD,WFQRI,WFQRA,WFSIDE,WFLINE, *MTTQ, *ERROR,WARN,RSEP1,RSEP3T,RSEP3R,RSEP3I,INDEXX,AP00, *FORM1,LOWER,LENGTH INTEGER LENGTH C WFEROR,WFREAD,WFERAS,WFTUBE,WFINTP, C WFQD,WFQRI,WFQRA,WFSIDE,WFLINE ... This file. C MTTQ ... File mttq.for. C ERROR,WARN ... File error.for. C RSEP1,RSEP3I,RSEP3T,RSEP3R ... C File sep.for. C INDEXX ... File indexx.for. C AP00 ... File ap.for. C FORM1,LOWER ... forms.for. C LENGTH ... length.for. C C Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C wfsrf.inc. C....................................................................... C Auxiliary storage locations: INTEGER KRAMP PARAMETER (KRAMP=2) INTEGER IRAY1,IRAY2,IRAY3 INTEGER IT1,IIT1 INTEGER I INTEGER LU1,LU2,LU3,LU4,LU5 PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=5) CHARACTER*80 FILSEP,FILCRT,FILE1,FILE2,FILE3,CH INTEGER I1,I2,I3,I4,J1,J2,J3,NPOINT,LEN1,LEN2,LENG INTEGER MLJPOL PARAMETER (MLJPOL=10) REAL AUX,Z(MOUT),OUTMIN,OUTMAX CHARACTER*20 FORMAT,FORMA1,FORMA2,TEXTS(MLJPOL) C C----------------------------------------------------------------------- C C Reading a name of the file with the input data: FILSEP=' ' WRITE(*,'(A)') '+WFSRF: Enter input filename: ' READ(*,*) FILSEP IF (FILSEP.EQ.' ') THEN C WFSRF-01 CALL ERROR('WFSRF-01: No input file specified') 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 WRITE(*,'(A)') '+WFSRF: Reading input data. ' C C Reading all the data from the SEP file into the memory: CALL RSEP1(LU1,FILSEP) C C Reading filenames of the files with computed rays C and triangles: CALL RSEP3T('CRTOUT',FILCRT,' ') FILE1='r01.out' FILE2='r01i.out' FILE3='t01.out' IF (FILCRT.NE.' ') THEN OPEN (LU2,FILE=FILCRT,FORM='FORMATTED',STATUS='OLD') READ (LU2,*) FILE1,CH,FILE2,FILE3 CLOSE (LU2) ENDIF C C Opening the CRT output files: OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=FILE3,FORM='FORMATTED',STATUS='OLD') C C Reading travel time(s) of the wavefront(s) to be triangulated, C filename(s) of the output file(s), C recording which quantities are to be written into the file(s): CALL WFQD C C Reading file with computed triangles, C sorting the rays in triangles: WRITE(*,'(A)') '+WFSRF: Reading the file with triangles.' NT=0 NRAMP=0 IRAYMI=999999 IRAYMA=0 READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3 10 CONTINUE READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3 IF (IRAY1.EQ.0) GOTO 20 IF (IRAY1.GT.IRAYMA) IRAYMA=IRAY1 IF (IRAY2.GT.IRAYMA) IRAYMA=IRAY2 IF (IRAY3.GT.IRAYMA) IRAYMA=IRAY3 IF (IRAY1.LT.IRAYMI) IRAYMI=IRAY1 IF (IRAY2.LT.IRAYMI) IRAYMI=IRAY2 IF ((IRAY3.NE.0).AND.(IRAY3.LT.IRAYMI)) IRAYMI=IRAY3 IF (IRAY1.LT.IRAY2) THEN I=IRAY1 IRAY1=IRAY2 IRAY2=I ENDIF IF (IRAY2.LT.IRAY3) THEN I=IRAY2 IRAY2=IRAY3 IRAY3=I ENDIF IF (IRAY1.LT.IRAY2) THEN I=IRAY1 IRAY1=IRAY2 IRAY2=I ENDIF IF ((IRAY1.LE.0).OR.(IRAY2.LE.0).OR.(IRAY3.LT.0)) THEN C WFSRF-02 CALL ERROR('WFSRF-02: Disorder in file CRT-T.') C The 'triangles' in the input file CRT-T should be written C as three nonnegative integers (indices of rays), two of them C being positive (the third may equal zero in 2-D computations). ENDIF IF (NRAMP+3.GT.MAXRAM) CALL WFEROR(1) NRAMP=NRAMP+1 IRAM(NRAMP)=IRAY1 NRAMP=NRAMP+1 IRAM(NRAMP)=IRAY2 NRAMP=NRAMP+1 IRAM(NRAMP)=IRAY3 NT=NT+1 GOTO 10 20 CONTINUE CLOSE(LU3) NR=IRAYMA-IRAYMI+1 C C Check for the 'triangles' in case of 2D computation: IF (IRAM(3).EQ.0) THEN DO 25, I1=6,3*NT,3 IF (IRAM(I1).NE.0) THEN C WFSRF-03 CALL ERROR('WFSRF-03: Disorder in triangles.') C The first 'triangle' is formed by two rays with zero C instead of the third ray. This identifies 2-D calculation C and all the 'triangles' should be formed by two rays and C zero instead of the third ray. This is not the case of the C current triangle. Either input file CRT-T C is wrong, or there is a bug in the code. In the latter case, C please, contact the author. ENDIF 25 CONTINUE L3D=.FALSE. ELSE L3D=.TRUE. ENDIF C C Sorting the triangles: WRITE(*,'(A)') '+WFSRF: Sorting the triangles. ' IF (NRAMP+2*NT.GT.MAXRAM) CALL WFEROR(1) DO 30, I1=NRAMP+NT+1,NRAMP+NT+NT RAM(I1)=FLOAT(IRAM((I1-NRAMP-NT-1)*3+1)) 30 CONTINUE CALL INDEXX(NT,RAM(NRAMP+NT+1),IRAM(NRAMP+1)) NRAMP=NRAMP+NT C C Setting the value of MRAMP: MRAMP=MAXRAM-KRAMP*(NWFT*NT*5 + NWFT*NR*NQ) IF (MRAMP.LE.NRAMP) CALL WFEROR(1) NPOI=MRAMP NPLGN=MAXRAM+1 C C C Forming an auxiliary array with information about when can be rays C erased from the memory ("deleting array"): IF (NRAMP+NR.GT.MRAMP) CALL WFEROR(1) DO 50, I1=NRAMP+1,NRAMP+NR IRAM(I1)=0 50 CONTINUE NRAMP=NRAMP+NR ORAYE=IRAYMI-1-4*NT DO 60, I2=3*NT+1,4*NT I1=(IRAM(I2)-1)*3+1 IRAM(IRAM(I1 )-ORAYE)=IRAM(I1) IRAM(IRAM(I1+1)-ORAYE)=IRAM(I1) IF (IRAM(I1+2).NE.0) IRAM(IRAM(I1+2)-ORAYE)=IRAM(I1) 60 CONTINUE C C C Forming an auxiliary array with information about addresses C of the ends of records for rays in array RAM ("addressing array"): C "Ray" IRAYMI-1: NRAMP=NRAMP+1 IF (NRAMP.GT.MRAMP) CALL WFEROR(1) IRAM(NRAMP)=NRAMP+NR+NWFT*NR C All other rays: IF (NRAMP+NR.GT.MRAMP) CALL WFEROR(1) DO 70, I1=NRAMP+1,NRAMP+NR IRAM(I1)=0 70 CONTINUE NRAMP=NRAMP+NR ORAYA=IRAYMI-1-4*NT-NR-1 C C C Forming an auxiliary array(s) with information about addresses C of the ends of records for points on wavefronts according to rays C in array RAM ("points at wavefronts addressing array(s)"): IF (NRAMP+NWFT*NR.GT.MRAMP) CALL WFEROR(1) DO 80, I1=NRAMP+1,NRAMP+NWFT*NR IRAM(I1)=0 80 CONTINUE NRAMP=NRAMP+NWFT*NR OWFT(1)=IRAYMI-1-4*NT-NR-NR-1 C OWFT(2)=IRAYMI-1-4*NT-NR-NR-1-NR C OWFT(3)=IRAYMI-1-4*NT-NR-NR-1-NR-NR C C C Loop for all the triangles: IIT1=1 I1=-1 100 CONTINUE IT1=(IRAM(3*NT+IIT1)-1)*3+1 I2=NINT((100.*IIT1)/(NT)) IF (I2.NE.I1) THEN WRITE(*,'(''+'',79('' ''))') WRITE(*,'(A,1I4,A)') * '+WFSRF: Interpolating. ',I2,'% of ray tubes completed.' I1=I2 ENDIF C C If necessary, reading new rays: IF * (((IRAM(IRAM(IT1)-ORAYA+1).EQ.0).AND.(IRAM(IT1).NE.IRAYMA)) .OR. * ((IRAM(IRAM(IT1)-ORAYA ).EQ.0).AND.(IRAM(IT1).EQ.IRAYMA))) * CALL WFREAD(LU1,LU2,IT1) C C Empting the array RAM to enable writing of new points C and polygons at wavefronts: IF ((MRAMP-NRAMP).LT.(NINT(MAXRAM/10.))) CALL WFERAS C C Creating polygons at wavefronts along the rays of the triangle: CALL WFTUBE(IT1) C IIT1=IIT1+1 IF (IIT1.LE.NT) GOTO 100 C End of the loop for all the triangles. CLOSE(LU1) CLOSE(LU2) C C C Storing points and polygons on wavefronts: WRITE(*,'(A)') *'+WFSRF: Writing. ' C C Unifying character strings in TEXTC and CHOUT: DO 110, I1=1,NQ CALL LOWER(CHOUT(I1)) IF (CHOUT(I1).EQ.'mp1') CHOUT(I1)='norm1' IF (CHOUT(I1).EQ.'mp2') CHOUT(I1)='norm2' IF (CHOUT(I1).EQ.'mp3') CHOUT(I1)='norm3' 110 CONTINUE C C Storing points: OPEN(LU1,FILE=VRTX,FORM='FORMATTED') WRITE(LU1,'(A)') '/' NPOINT=(NPOI-MRAMP)/NQ C Length of the names of points: LEN1=LENGTH(TEXTP) LEN2=0 202 CONTINUE IF (NPOINT.GE.10**LEN2) THEN LEN2=LEN2+1 GOTO 202 ENDIF LENG=LEN1+LEN2 C Loop over points: DO 218, I1=1,NPOINT C Address in RAM just before the current point: J1=MRAMP+(I1-1)*NQ C Name of the point: DO 204, I2=0,LEN2-1 TEXTP(LENG-I2:LENG-I2)= * CHAR(ICHAR('0')+MOD(I1,10**(I2+1))/10**I2) 204 CONTINUE C The quantities at the point: J2=0 C Loop over the quantities to be stored: 206 CONTINUE IF ((TEXTC(J2+1).NE.' ').AND.(J2+1.LE.MOUT)) THEN J2=J2+1 DO 210, I2=1,NQ IF (TEXTC(J2).EQ.CHOUT(I2)) THEN J3=I2 GOTO 211 ENDIF 210 CONTINUE C WFSRF-23 CALL ERROR ('WFSRF-23: No CHOUT found for TEXTC.') C This error should not appear. C Please contact the author. 211 CONTINUE IF ((TEXTC(J2).EQ.'icb').OR.(TEXTC(J2).EQ.'mhist')) THEN AUX=FLOAT(IRAM(J1+J3)) ELSE AUX=RAM(J1+J3) ENDIF IF (POWER(J2).EQ.1.) THEN Z(J2)=AUX ELSE Z(J2)=AUX**POWER(J2) ENDIF GOTO 206 ENDIF C End of the loop for quantities to be stored. C Writing the quantities at the point: C Setting output format for the array: FORMAT='(3A,00(F00.0,1X),A)' FORMAT(6:6)=CHAR(ICHAR('0')+MOD(J2,10)) FORMAT(5:5)=CHAR(ICHAR('0')+J2/10) OUTMIN=0. OUTMAX=0. DO 214, I2=1,J2 IF (OUTMIN.GT.Z(I2)) OUTMIN=Z(I2) IF (OUTMAX.LT.Z(I2)) OUTMAX=Z(I2) 214 CONTINUE CALL FORM1(OUTMIN,OUTMAX,FORMAT(8:15)) FORMAT(14:17)= '1X),' C Output format is set. WRITE(LU1,FORMAT) '''',TEXTP(1:LENG),''' ',(Z(I2),I2=1,J2),'/' 218 CONTINUE C End of the loop over all points. WRITE(LU1,'(A)') '/' CLOSE(LU1) C C Storing polygons: IF(PLGN.NE.' ') OPEN(LU1,FILE=PLGN,FORM='FORMATTED') IF(PLGNS.NE.' ')OPEN(LU2,FILE=PLGNS,FORM='FORMATTED') FORMA1='(00(I8,1X),A)' FORMA2='(00(3A,1X),A)' I1=MAXRAM-1 C Loop over polygons: 220 CONTINUE IF (I1.GT.NPLGN) THEN J1=IRAM(I1) I3=MOD(J1,10) FORMA1(3:3)=CHAR(ICHAR('0')+I3) FORMA2(3:3)=CHAR(ICHAR('0')+I3) I3=J1/10 FORMA1(2:2)=CHAR(ICHAR('0')+I3) FORMA2(2:2)=CHAR(ICHAR('0')+I3) J2=IRAM(I1-1) IF (PLGN.NE.' ') THEN WRITE(LU1,FORMA1) * ((IRAM(I2)-MRAMP)/NQ,I2=I1-1,I1-J1,-1),'/' ENDIF IF (PLGNS.NE.' ') THEN IF (J1.GT.MLJPOL) THEN C WFSRF-04 CALL ERROR ('WFSRF-04: Small array for names of polygons.') C This error may be caused by small dimension MLJPOL of C array TEXTS defined in this program. ENDIF DO 224, I2=I1-1,I1-J1,-1 C Index of the point (from 1 to NPOINT): J2=(IRAM(I2)-MRAMP)/NQ C Name of the point: I4=I1-I2 DO 222, I3=0,LEN2-1 TEXTS(I4)(1:LEN1)=TEXTP(1:LEN1) TEXTS(I4)(LENG-I3:LENG-I3)= * CHAR(ICHAR('0')+MOD(J2,10**(I3+1))/10**I3) 222 CONTINUE 224 CONTINUE WRITE(LU2,FORMA2) * ('''',TEXTS(I2)(1:LENG),'''',I2=1,J1),'/' ENDIF I1=I1-(J1+2) GOTO 220 ENDIF C End of the loop over polygons. IF (PLGN.NE.' ') THEN WRITE(LU1,'(A)') '/' CLOSE(LU1) ENDIF IF (PLGNS.NE.' ') THEN WRITE(LU2,'(A)') '/' CLOSE(LU2) ENDIF C WRITE(*,'(A)') *'+WFSRF: Done. ' STOP END C C C======================================================================= C SUBROUTINE WFREAD(LU1,LU2,IT1) C C---------------------------------------------------------------------- C Subroutine to read the unformatted output of program CRT and C to write it into array (I)RAM used in program WFSRF. C Reading the output files is completed by a simple invocation of C subroutine AP00 of file 'ap.for'. C INTEGER LU1,LU2,IT1 C Input: C LU1 ... Number of logical unit corresponding to the file with C the quantities stored along rays (see C.R.T.5.5.1). C LU2 ... Number of logical unit corresponding to the file with C the quantities at the initial points of rays, C corresponding to file 'RAYPOINTS' (see C.R.T.6.1). C IT1 ... Position of the first ray of the actually processed C triangle in the array IRAM. C No output. 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 Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C wfsrf.inc. C....................................................................... C C Auxiliary storage locations: INTEGER IRAY0,I1 C IRAY0.. Index of the last ray read in into the array RAM. C C----------------------------------------------------------------------- C Loop for the points of rays: 10 CONTINUE IF ((NRAMP+2*NQ).GT.MRAMP) THEN C Empting the memory: CALL WFERAS IF ((NRAMP+2*NQ).GT.MRAMP) CALL WFEROR(1) ENDIF C Reading the results of the complete ray tracing: IRAY0=IRAY IWAVE0=IWAVE IF ((IRAY0.EQ.0).AND.(IWAVE0.EQ.0)) THEN C Reading the first point on a ray of the first wave: CALL AP00(0,LU1,LU2) IF (IWAVE.LT.1) GOTO 20 ELSEIF (IWAVE.EQ.IWAVE0) THEN C Reading the next point on a ray of the actual wave: CALL AP00(0,LU1,LU2) IF (IWAVE.NE.IWAVE0) GOTO 20 ENDIF IF (IRAY.LT.IRAYMI) GOTO 10 IF (IRAY.GT.IRAYMA) GOTO 10 IF ((IRAY-IRAY0).GE.2) THEN C Some rays skipped by AP00: DO 15, I1=IRAY0+1,IRAY-1 IF (I1.GE.IRAYMI) THEN IRAM(I1-ORAYA)=IRAM(IRAY0-ORAYA) ENDIF 15 CONTINUE ENDIF IF (IRAM(IRAY-ORAYE).NE.0) THEN C Writing the results of the complete ray tracing - recording C new point on a ray: IF (IPT.LE.1) THEN IF ((ISHEET.EQ.0).AND.(IRAY.EQ.IRAYMI)) THEN C WFSRF-05 CALL WARN * ('WFSRF-05: a ray with history equal to 0 was observed') C All the rays are probably of the same history 0. Only the C initial-value ray tracing was performed. Width and shape C of ray tubes were not controlled. The interpolation C in such a case makes sense only in smooth models. C Check the value of the parameter IPOINT in CRT input data C RPAR (4). ENDIF C New ray - recording the initial point: CALL WFQRI ENDIF C Recording the new point: CALL WFQRA ENDIF IRAM(IRAY-ORAYA)=NRAMP IF ((IPT.LE.1).AND.(IRAY.GT.IRAM(IT1)).AND.(IRAY.NE.IRAYMA)) * RETURN GOTO 10 C 20 CONTINUE C End of rays: IF (IRAY0.LT.IRAYMA) THEN C Some rays at the end of the CRT output file missing: DO 22, I1=IRAY0+1,IRAYMA IF (I1.GE.IRAYMI) THEN IRAM(I1-ORAYA)=NRAMP ENDIF 22 CONTINUE ENDIF RETURN END C C C======================================================================= C SUBROUTINE WFERAS C C---------------------------------------------------------------------- C Subroutine for empting the array (I)RAM. All the parameters C of all the rays, which will no more be used, are erased. C C No input. C No output. C C ........................... C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C IRAY .. Index of the ray being actually read in by CIREAD. C This procedure supposes, that any ray with higher C index than IRAY was not read in. C None of the storage locations of the common block are altered. C ........................... C Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C wfsrf.inc. C....................................................................... C Auxiliary storage locations: INTEGER I1,I2,J1 INTEGER IADDRP C I1 ... Controls main loop over rays. C I2 ... Controls the loop over parameters of ray I1. C J1 ... address of the last used record of array RAM. C C----------------------------------------------------------------------- J1=IRAM(IRAYMI-1-ORAYA) IADDRP=J1 C Loop for the rays: DO 20, I1=IRAYMI,IRAY IF (IRAM(I1-ORAYE).GE.(IRAY-1)) THEN C This ray is not to be erased: DO 10, I2=IADDRP+1,IRAM(I1-ORAYA) J1=J1+1 IRAM(J1)=IRAM(I2) 10 CONTINUE ELSE C This ray is to be erased: IRAM(I1-ORAYE)=0 ENDIF IADDRP=IRAM(I1-ORAYA) IRAM(I1-ORAYA)=J1 20 CONTINUE NRAMP=J1 RETURN END C======================================================================= C SUBROUTINE WFTUBE(IT1) C C---------------------------------------------------------------------- C Subroutine for interpolation within ray tube formed by the rays C IRAM(IT1), IRAM(IT1+1), IRAM(IT1+2). C INTEGER IT1 C Input: C IT1 ... The address of the index of the first ray of the ray tube, C in which the interpolation is to be performed. C No output. C C ........................... C Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C wfsrf.inc. C....................................................................... INTEGER IRAY1,IRAY2,IRAY3 INTEGER I1MI,I2MI,I3MI,I1MA,I2MA,I3MA,I1,I2,J1,J2 INTEGER KRAY1,KRAY2,KRAY3,ICBR1,ICBR2,ICBR3,NEXT INTEGER NPOL,MPOL PARAMETER (MPOL=10) INTEGER KPOL(MPOL) LOGICAL LSHIFT C J1 ... C J2 ... C C----------------------------------------------------------------------- IRAY1=IRAM(IT1 ) IRAY2=IRAM(IT1+1) IRAY3=IRAM(IT1+2) IF ( (IRAM(IRAY1-ORAYA).EQ.0).OR. * (IRAM(IRAY2-ORAYA).EQ.0).OR. * (L3D.AND.(IRAM(IRAY3-ORAYA).EQ.0))) THEN C WFSRF-06 CALL ERROR * ('WFSRF-06: Parameters of a ray not found in WFTUBE.') C This error may be caused by C K2P C not equal to zero, then only two-point rays are stored in C output files of CRT. We recommend to run CRT with file C writall.dat. ENDIF C I1MI= IRAM(IRAY1-1-ORAYA) I2MI= IRAM(IRAY2-1-ORAYA) IF (L3D) I3MI=IRAM(IRAY3-1-ORAYA) I1MA= IRAM(IRAY1-ORAYA)-NQ I2MA= IRAM(IRAY2-ORAYA)-NQ IF (L3D) I3MA=IRAM(IRAY3-ORAYA)-NQ IF ((I1MI.GT.I1MA).OR.(I2MI.GT.I2MA).OR.(L3D.AND.(I3MI.GT.I3MA))) * RETURN C C Loop for wavefronts, searching for a polygon C which is an intersection of the tube with the wavefront: DO 100, I1=1,NWFT C KRAY1=0 KRAY2=0 KRAY3=0 ICBR1=0 ICBR2=0 ICBR3=0 C C Computing point on the first ray: IF (IRAM(IRAY1-OWFT(I1)).EQ.0) THEN C The point is not computed, computing the point: IF ((RAM(I1MI+6).LE.WFTIME(I1)).AND. * (RAM(I1MA+6).GE.WFTIME(I1))) THEN C The point exists, computing the point: DO 10, I2=I1MI+NQ,I1MA,NQ IF (RAM(I2+6).GE.WFTIME(I1)) THEN CALL WFINTP(I2-NQ,I2,WFTIME(I1)) KRAY1=NPOI ICBR1=IRAM(KRAY1-NQ+5) IRAM(IRAY1-OWFT(I1))=NPOI GOTO 11 ENDIF 10 CONTINUE 11 CONTINUE ELSE C The point does not exist: IRAM(IRAY1-OWFT(I1))=-1 ENDIF ELSEIF (IRAM(IRAY1-OWFT(I1)).EQ.-1) THEN C The point does not exist: CONTINUE ELSE C The point is already computed: KRAY1=IRAM(IRAY1-OWFT(I1)) ICBR1=IRAM(KRAY1-NQ+5) ENDIF C C Computing point on the second ray: IF (IRAM(IRAY2-OWFT(I1)).EQ.0) THEN C The point is not computed, computing the point: IF ((RAM(I2MI+6).LE.WFTIME(I1)).AND. * (RAM(I2MA+6).GE.WFTIME(I1))) THEN C The point exists, computing the point: DO 20, I2=I2MI+NQ,I2MA,NQ IF (RAM(I2+6).GE.WFTIME(I1)) THEN CALL WFINTP(I2-NQ,I2,WFTIME(I1)) KRAY2=NPOI ICBR2=IRAM(KRAY2-NQ+5) IRAM(IRAY2-OWFT(I1))=NPOI GOTO 21 ENDIF 20 CONTINUE 21 CONTINUE ELSE C The point does not exist: IRAM(IRAY2-OWFT(I1))=-1 ENDIF ELSEIF (IRAM(IRAY2-OWFT(I1)).EQ.-1) THEN C The point does not exist: CONTINUE ELSE C The point is already computed: KRAY2=IRAM(IRAY2-OWFT(I1)) ICBR2=IRAM(KRAY2-NQ+5) ENDIF C C Computing point on the third ray: IF (L3D) THEN IF (IRAM(IRAY3-OWFT(I1)).EQ.0) THEN C The point is not computed, computing the point: IF ((RAM(I3MI+6).LE.WFTIME(I1)).AND. * (RAM(I3MA+6).GE.WFTIME(I1))) THEN C The point exists, computing the point: DO 30, I2=I3MI+NQ,I3MA,NQ IF (RAM(I2+6).GE.WFTIME(I1)) THEN CALL WFINTP(I2-NQ,I2,WFTIME(I1)) KRAY3=NPOI ICBR3=IRAM(KRAY3-NQ+5) IRAM(IRAY3-OWFT(I1))=NPOI GOTO 31 ENDIF 30 CONTINUE 31 CONTINUE ELSE C The point does not exist: IRAM(IRAY3-OWFT(I1))=-1 ENDIF ELSEIF (IRAM(IRAY3-OWFT(I1)).EQ.-1) THEN C The point does not exist: CONTINUE ELSE C The point is already computed: KRAY3=IRAM(IRAY3-OWFT(I1)) ICBR3=IRAM(KRAY3-NQ+5) ENDIF ENDIF C C Recording point on the first ray: IF (KRAY1.NE.0) THEN KPOL(1)=KRAY1 NPOL=1 ELSE NPOL=0 ENDIF C C Computing and recording points C between the first and the second ray: IF (ICBR1.NE.ICBR2) THEN NEXT=0 40 CONTINUE CALL WFSIDE(IRAY1,KRAY1,ICBR1,IRAY2,KRAY2,ICBR2, * I1MI,I1MA,I2MI,I2MA,WFTIME(I1),NEXT) NPOL=NPOL+1 IF (NPOL.GT.MPOL) CALL WFEROR(4) KPOL(NPOL)=NPOI IF (NEXT.NE.0) GOTO 40 ENDIF C C Recording point on the second ray: IF (KRAY2.NE.0) THEN NPOL=NPOL+1 IF (NPOL.GT.MPOL) CALL WFEROR(4) KPOL(NPOL)=KRAY2 ENDIF C IF (L3D) THEN C Computing and recording points C between the second and the third ray: IF (ICBR2.NE.ICBR3) THEN NEXT=0 50 CONTINUE CALL WFSIDE(IRAY2,KRAY2,ICBR2,IRAY3,KRAY3,ICBR3, * I2MI,I2MA,I3MI,I3MA,WFTIME(I1),NEXT) NPOL=NPOL+1 IF (NPOL.GT.MPOL) CALL WFEROR(4) KPOL(NPOL)=NPOI IF (NEXT.NE.0) GOTO 50 ENDIF C C Recording point on the third ray: IF (KRAY3.NE.0) THEN NPOL=NPOL+1 IF (NPOL.GT.MPOL) CALL WFEROR(4) KPOL(NPOL)=KRAY3 ENDIF C C Computing and recording points C between the third and the first ray: IF (ICBR3.NE.ICBR1) THEN NEXT=0 60 CONTINUE CALL WFSIDE(IRAY3,KRAY3,ICBR3,IRAY1,KRAY1,ICBR1, * I3MI,I3MA,I1MI,I1MA,WFTIME(I1),NEXT) NPOL=NPOL+1 IF (NPOL.GT.MPOL) CALL WFEROR(4) KPOL(NPOL)=NPOI IF (NEXT.NE.0) GOTO 60 ENDIF ENDIF C C The polygon is found, shifting it in case that it contains C points in more than one complex block: LSHIFT=.FALSE. IF (L3D) THEN DO 70, I2=2,NPOL IF (IRAM(KPOL(1)-NQ+5).NE.IRAM(KPOL(I2)-NQ+5)) * LSHIFT=.TRUE. 70 CONTINUE ENDIF IF (LSHIFT) THEN 72 CONTINUE IF (IRAM(KPOL(1)-NQ+5).EQ.IRAM(KPOL(NPOL)-NQ+5)) THEN IF (NPOL+1.GT.MPOL) CALL WFEROR(4) KPOL(NPOL+1)=KPOL(1) DO 74, I2=1,NPOL KPOL(I2)=KPOL(I2+1) 74 CONTINUE GOTO 72 ENDIF ENDIF C IF (NPOL.GE.2) THEN C Recording the polygon to RAM: J2=1 C Loop over polygons: 80 CONTINUE NPLGN=NPLGN-1 IF (NPLGN.LE.NPOI) CALL WFEROR(3) RAM(NPLGN)=WFTIME(I1) NPLGN=NPLGN-1 IF (NPLGN.LE.NPOI) CALL WFEROR(3) J1=NPLGN IRAM(J1)=1 NPLGN=NPLGN-1 IF (NPLGN.LE.NPOI) CALL WFEROR(3) IRAM(NPLGN)=KPOL(J2) C Loop over points of the polygon: 82 CONTINUE J2=J2+1 IF (J2.GT.NPOL) GOTO 89 IF (IRAM(KPOL(J2)-NQ+5).NE.IRAM(KPOL(J2-1)-NQ+5)) * GOTO 80 IRAM(J1)=IRAM(J1)+1 NPLGN=NPLGN-1 IF (NPLGN.LE.NPOI) CALL WFEROR(3) IRAM(NPLGN)=KPOL(J2) GOTO 82 89 CONTINUE ENDIF C 100 CONTINUE C End of the loop over wavefronts. C RETURN END C C======================================================================= C SUBROUTINE WFQD C C....................................................................... C C Subroutine designed to read the input data describing the names C of the output files with the interpolated quantities, to read the C quantities computed by CRT, to interpolate and to write the output C files with the interpolated quantities. C C....................................................................... C Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C For the description of the parameters stored in the common block C WFSRFC refer to the file wfsrf.inc. 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 Auxiliary storage locations for all the entries: INTEGER I1,I2 CHARACTER*20 FORMA1,FORMA2,FORMA3 CHARACTER*5 QNAMEL(MOUT),TEXT,TEXTL C....................................................................... C C Reading the wavefront travel time: CALL RSEP3R('WFTIME',WFTIME(1),1.) NWFT=1 C C Recalling the output filenames: CALL RSEP3T('VRTX',VRTX,'vrtx.out') CALL RSEP3T('PLGN',PLGN,'plgn.out') CALL RSEP3T('PLGNS',PLGNS,' ') C C Recalling the first part of names of points in output file VRTX: CALL RSEP3T('TEXTP',TEXTP,'V') C C Reading the list of the quantities which may be interpolated: CALL MTTQ DO 3, I1=1,MOUT QNAMEL(I1)=QNAMES(I1) IF (QNAMEL(I1).NE.' ') THEN CALL LOWER(QNAMEL(I1)) ENDIF 3 CONTINUE C C Recalling the data specifying the quantities to be written into C the individual columns of the output file VRTX with points C wavefronts: FORMA1='COLUMN01' FORMA2='POWER01' FORMA3='IVALUE01' I1=1 5 CONTINUE IF (I1.EQ.1) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'X1') ELSEIF (I1.EQ.2) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'X2') ELSEIF (I1.EQ.3) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'X3') ELSEIF (I1.EQ.4) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'NORM1') ELSEIF (I1.EQ.5) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'NORM2') ELSEIF (I1.EQ.6) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'NORM3') ELSE CALL RSEP3T(FORMA1,TEXTC(I1),' ') ENDIF CALL RSEP3R(FORMA2,POWER(I1),1.) CALL RSEP3I(FORMA3,IVALUE(I1),1) IF (TEXTC(I1).NE.' ') THEN CALL LOWER(TEXTC(I1)) IF ((TEXTC(I1).NE.'x1').AND.(TEXTC(I1).NE.'x2').AND. * (TEXTC(I1).NE.'x3').AND.(TEXTC(I1).NE.'norm1').AND. * (TEXTC(I1).NE.'norm2').AND.(TEXTC(I1).NE.'norm3').AND. * (TEXTC(I1).NE.'icb').AND. * (TEXTC(I1).NE.'mtt').AND. * (TEXTC(I1).NE.'mhist')) THEN DO 7, I2=1,MOUT IF (TEXTC(I1).EQ.QNAMEL(I2)) GOTO 8 7 CONTINUE C WFSRF-07 CALL ERROR('WFSRF-07: Wrong value of COLUMNii.') C See allowed values of COLUMNii in the C description of file SEP. ENDIF 8 CONTINUE I1=I1+1 FORMA1(8:8)=CHAR(ICHAR('0')+MOD(I1,10)) FORMA2(7:7)=FORMA1(8:8) FORMA3(8:8)=FORMA1(8:8) FORMA1(7:7)=CHAR(ICHAR('0')+I1/10) FORMA2(6:6)=FORMA1(7:7) FORMA3(7:7)=FORMA1(7:7) IF (I1.GT.MOUT) THEN CALL RSEP3T(FORMA1,TEXT,' ') IF (TEXT.NE.' ') THEN C WFSRF-08 CALL ERROR('WFSRF-08: More than MOUT COLUMNs.') C Currently up to MOUT values of COLUMNii may be specified. C See allowed values of COLUMNii in the C description of file SEP. ENDIF ELSE GOTO 5 ENDIF ENDIF C End of the loop over future columns of the output file. C C C Recording which quantities are to be computed and stored C in the points on the rays and in the points on the wavefronts: CHOUT(1)='X1' CHOUT(2)='X2' CHOUT(3)='X3' CHOUT(4)='ISRF' CHOUT(5)='ICB' CHOUT(6)='MTT' CHOUT(7)='MP1' CHOUT(8)='MP2' CHOUT(9)='MP3' NQ=9 DO 10, I1=1,MOUT IF (TEXTC(I1).EQ.'mhist') THEN NQ=NQ+1 IF (NQ.GT.MOUT) CALL WFEROR(5) CHOUT(NQ)='MHIST' ENDIF 10 CONTINUE DO 20, I2=1,MOUT IF (QNAMES(I2).EQ.' ') GOTO 21 TEXT=QNAMES(I2) TEXTL=QNAMEL(I2) DO 15, I1=1,MOUT IF (TEXTC(I1).EQ.TEXTL) THEN NQ=NQ+1 IF (NQ.GT.MOUT) CALL WFEROR(5) CHOUT(NQ)=TEXT CALL MTTQD(TEXT) ENDIF 15 CONTINUE 20 CONTINUE 21 CONTINUE C RETURN C C======================================================================= C ENTRY WFQRI C C....................................................................... C C Entry designed to read the output of the CRT at an initial point of C a ray. C C Coordinates: RAM(NRAMP+1)=YI(3) RAM(NRAMP+2)=YI(4) RAM(NRAMP+3)=YI(5) C Index of surface: IRAM(NRAMP+4)=0 C Index of complex block: IRAM(NRAMP+5)=ICB1I C Travel time: RAM(NRAMP+6)=YI(1) C Three components of the slowness vector: RAM(NRAMP+7)=YI(6) RAM(NRAMP+8)=YI(7) RAM(NRAMP+9)=YI(8) C Other quantities to be interpolated: DO 30, I1=10,NQ IF (CHOUT(I1).EQ.'MHIST') THEN C Ray history: IRAM(NRAMP+I1)=ISHEET ELSE C User-defined quantity: RAM(NRAMP+I1)=0. CALL MTTQRI(CHOUT(I1),RAM(NRAMP+I1)) ENDIF 30 CONTINUE NRAMP=NRAMP+NQ RETURN C======================================================================= C ENTRY WFQRA C C....................................................................... C C Entry designed to read the output of the CRT at other than initial C point of a ray. C C Coordinates: RAM(NRAMP+1)=Y(3) RAM(NRAMP+2)=Y(4) RAM(NRAMP+3)=Y(5) C Index of surface: IRAM(NRAMP+4)=ISRF C Index of complex block: IRAM(NRAMP+5)=ICB1 C Travel time: RAM(NRAMP+6)=Y(1) C Three components of the slowness vector: RAM(NRAMP+7)=Y(6) RAM(NRAMP+8)=Y(7) RAM(NRAMP+9)=Y(8) C Other quantities to be interpolated: DO 40, I1=10,NQ IF (CHOUT(I1).EQ.'MHIST') THEN C Ray history: IRAM(NRAMP+I1)=ISHEET ELSE C User-defined quantity: RAM(NRAMP+I1)=0. CALL MTTQRA(CHOUT(I1),RAM(NRAMP+I1)) ENDIF 40 CONTINUE NRAMP=NRAMP+NQ RETURN END C C======================================================================= C SUBROUTINE WFINTP(IPA,IPB,TI) C C---------------------------------------------------------------------- INTEGER IPA,IPB REAL TI C C Subroutine for linear interpolation between two points A and B. C C---------------------------------------------------------------------- C C Input: C IPA,IPB ... Indices of positions of points A and B in array RAM. C The points must be in the same complex block. C TI ... Time for interpolation. Time TI is supposed to be between C the travel times TA and TB in points A and B. C C No output. C The quantities in the interpolated point are written into the C array (I)RAM. C C....................................................................... C C Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C For the description of the parameters stored in the common block C WFSRFC refer to the file wfsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I1 REAL WA,WB,AUX REAL TA,TB C----------------------------------------------------------------------- C IF (.NOT.(IRAM(IPA+5).EQ.IRAM(IPB+5))) THEN C WFSRF-09 CALL ERROR('WFSRF-09: Different indices of complex blocks.') C This error should not appear. C Please contact the author. ENDIF C IF (NPOI+NQ.GT.NPLGN) CALL WFEROR(2) C TA =RAM(IPA+6) TB =RAM(IPB+6) C AUX=TB-TA IF ((TB.EQ.TI).OR.(AUX.EQ.0.)) THEN WA=0. WB=1. ELSEIF (TA.EQ.TI) THEN WA=1. WB=0. ELSE WB=(TI-TA)/AUX WA=1.-WB ENDIF C DO 10, I1=1,NQ IF (CHOUT(I1).EQ.'ISRF') THEN IF (IRAM(IPA+I1).EQ.IRAM(IPB+I1)) THEN IRAM(NPOI+I1)=IRAM(IPA+I1) ELSE IRAM(NPOI+I1)=0 ENDIF ELSEIF (CHOUT(I1).EQ.'ICB') THEN IRAM(NPOI+I1)=IRAM(IPA+I1) ELSEIF (CHOUT(I1).EQ.'MHIST') THEN IRAM(NPOI+I1)=IRAM(IPA+I1) ELSE RAM(NPOI+I1)=WA*RAM(IPA+I1)+WB*RAM(IPB+I1) ENDIF 10 CONTINUE C NPOI=NPOI+NQ RETURN END C C======================================================================= C SUBROUTINE WFSIDE(IRAY1,KRAY1,ICBR1,IRAY2,KRAY2,ICBR2, * I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT) C C....................................................................... C C Subroutine for search for points between rays at a side of a ray tube. C C---------------------------------------------------------------------- C C Input: C IRAY1,IRAY2 ... Indices of the two rays. C KRAY1,KRAY2 ... Addresses of the interpolated points on the two C rays, or zeros if the points do not exist. C ICBR1,ICBR2 ... Indices of complex blocks of the above two points. C I1MI,I1MA,I2MI,I2MA ... Minimum and maximum addresses of the points C of the two rays IRAY1,IRAY2. C WFTIM ... Wawefront travel time. C NEXT ... Indicator of the direction of the search for the C interpolated points at the side of the tube. C NEXT=0 ... The search begins. C NEXT=1 ... Upwards search. C NEXT=-1 ... Downwards search. C Output: C NEXT ... Indicator of the direction of the search for the C interpolated points at the side of the tube. C NEXT=0 ... The search is terminated. C NEXT=1 ... Upwards search. C NEXT=-1 ... Downwards search. C C---------------------------------------------------------------------- C C Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C For the description of the parameters stored in the common block C WFSRFC refer to the file wfsrf.inc. C ........................... INTEGER IRAY1,KRAY1,ICBR1,IRAY2,KRAY2,ICBR2, * I1MI,I1MA,I2MI,I2MA,NEXT REAL WFTIM REAL TMIN,TMAX INTEGER IPA,IPB SAVE IPA,IPB C....................................................................... C IF (NEXT.EQ.0) THEN C Search for first point: IF (KRAY1.EQ.0) THEN C The point must be at the top of the tube: IPA=IRAM(IRAY1-ORAYA)-NQ IPB=IRAM(IRAY2-ORAYA)-NQ TMIN=AMIN1(RAM(IPA+6),RAM(IPB+6)) TMAX=AMAX1(RAM(IPA+6),RAM(IPB+6)) IF (.NOT.((TMIN.LE.WFTIM).AND.(WFTIM.LE.TMAX))) THEN C WFSRF-10 CALL ERROR('WFSRF-10: Wrong invocation of WFSIDE.') C This error should not appear. C Please contact the author. ENDIF CALL WFINTP(IPA,IPB,WFTIM) IF (IRAM(IPA+5).NE.ICBR2) THEN C The new point is in another block than the point on ray 2, C there must be some other points: NEXT=-1 ENDIF ELSE C The point must be within (or at the top of) the tube. C Starting from the source: IPA=I1MI IPB=I2MI NEXT=1 CALL WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT) IF (IRAM(IPA+5).NE.ICBR1) THEN C Starting from the top of the tube: IPA=I1MA IPB=I2MA NEXT=-1 CALL WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT) IF (IRAM(IPA+5).NE.ICBR1) THEN C WFSRF-11 CALL ERROR('WFSRF-11: Point in propper block not found.') C The point in the same complex block were found neither C in the upwards search, nor in the downwards search. C This error should not appear. C Please contact the author. ENDIF ENDIF CALL WFINTP(IPA,IPB,WFTIM) IF ((IRAM(IPA+5).EQ.ICBR2).OR. * ((ICBR2.EQ.0).AND.(IRAM(IPA+5).EQ.IRAM(I2MA+5)))) THEN C The new point is in the same block as the point on ray B: NEXT=0 ENDIF ENDIF ELSE C Search upwards or downwards the tube: IPA=IPA+NEXT*NQ IPB=IPB+NEXT*NQ CALL WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT) CALL WFINTP(IPA,IPB,WFTIM) IF ((IRAM(IPA+5).EQ.ICBR2).OR. * ((ICBR2.EQ.0).AND.(IRAM(IPA+5).EQ.IRAM(I2MA+5)))) THEN C The new point is in the same block as the point on ray B: NEXT=0 ENDIF ENDIF RETURN END C C======================================================================= C SUBROUTINE WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT) C C....................................................................... C C Subroutine for a search for a line at a side of a ray tube. C C Common block /WFSRFC/ to store triangles and rays: INCLUDE 'wfsrf.inc' C For the description of the parameters stored in the common block C WFSRFC refer to the file wfsrf.inc. C............................ INTEGER IPA,IPB,I1MI,I1MA,I2MI,I2MA,NEXT REAL WFTIM C INTEGER NQL REAL TMIN,TMAX C IF ((NEXT.NE.1).AND.(NEXT.NE.-1)) THEN C WFSRF-12 CALL ERROR('WFSRF-12: Wrong invocation of WFLINE.') C This error should not appear. C Please contact the author. ENDIF NQL=NEXT*NQ C 10 CONTINUE TMIN=AMIN1(RAM(IPA+6),RAM(IPB+6)) TMAX=AMAX1(RAM(IPA+6),RAM(IPB+6)) IF ((TMIN.LE.WFTIM).AND.(WFTIM.LE.TMAX)) THEN C The point is between IPA and IPB: RETURN ENDIF 20 CONTINUE IPA=IPA+NQL IF ((IPA.LT.I1MI).OR.(IPA.GT.I1MA)) THEN C WFSRF-13 CALL ERROR('WFSRF-13: Point A not reached.') C The point on first ray not reached. C This error should not appear. C Please contact the author. ENDIF IF (.NOT.((IRAM(IPA+4).NE.0).OR. * (IPA.EQ.I1MI).OR.(IPA.EQ.I1MA))) GOTO 20 30 CONTINUE IPB=IPB+NQL IF ((IPB.LT.I2MI).OR.(IPB.GT.I2MA)) THEN C WFSRF-14 CALL ERROR('WFSRF-14: Point B not reached.') C The point on second ray not reached. C This error should not appear. C Please contact the author. ENDIF IF (.NOT.((IRAM(IPB+4).NE.0).OR. * (IPB.EQ.I2MI).OR.(IPB.EQ.I2MA))) GOTO 30 GOTO 10 END C C======================================================================= C SUBROUTINE WFEROR(IERR) C C....................................................................... C C Subroutine to write error messages. C INTEGER IERR IF (IERR.EQ.1) THEN C WFSRF-15 CALL ERROR('WFSRF-15: Small array RAM.') C Dimensions of array RAM exceeded when writing points C on rays. C Try to enlarge the dimension MRAM in the file C ram.inc. ELSEIF (IERR.EQ.2) THEN C WFSRF-16 CALL ERROR('WFSRF-16: Small array RAM.') C Dimensions of array RAM exceeded when writing points C at wavefronts. C Try to enlarge the dimension MRAM in the file C ram.inc. ELSEIF (IERR.EQ.3) THEN C WFSRF-17 CALL ERROR('WFSRF-17: Small array RAM.') C Dimensions of array RAM exceeded when writing polygons C at wavefronts. C Try to enlarge the dimension MRAM in the file C ram.inc. ELSEIF (IERR.EQ.4) THEN C WFSRF-18 CALL ERROR('WFSRF-18: Small array KPOL.') C This error may be caused by small dimension MPOL of array C KPOL defined in subroutine WFTUBE. Try to enlarge MPOL. ELSEIF (IERR.EQ.5) THEN C WFSRF-20 CALL ERROR('WFSRF-20: Small array CHOUT.') C This error may be caused by small dimension MOUT of array CHOUT. C Try to enlarge MOUT in the include file C wfsrf.inc. ELSE C WFSRF-19 CALL ERROR('WFSRF-19: Wrong IERR.') C This subroutine was invocated with wrong value of IERR. C This error should not appear. C Please contact the author. ENDIF END C C C======================================================================= C INCLUDE 'mttq.for' C mttq.for INCLUDE 'error.for' C error.for INCLUDE 'ap.for' C ap.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'indexx.for' C indexx.for C C======================================================================= C