C
C Program CRTLEW to calculate the directions of the projections of the C ray segments into the given 2-D section and to collect them according C to the given angular intervals. C C Version: 5.70 C Date: 2003, May 18 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 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 the input and output files: C CRTOUT='string'... File with the names of the output files of C program CRT. A single set of CRT output files is read C from CRTOUT. Only the files 'CRT-R', and 'CRT-I' are C read by CRTLEW. C For general description of file CRTOUT refer to C file writ.for. C Default: CRTOUT=' ' which means 'CRT-R'='r01.out', C 'CRT-I'='r01i.out' C MODLEW='string' ... Name of the output formatted file with the C distribution of the ray segments directions according C to the given angular intervals. C Description of file MODLEW C Default: MODLEW='modlew.out' C Switch between all rays and two-point rays only: C KALL=integer: C KALL.LE.0: only two-point rays are considered, C KALL.GE.1: all rays are considered. C Default: KALL=0 C Specification of the 2-D section for the projection of the rays: C KOOR1=integer... Index of the first coordinate determining the C 2-D section for the calculation of the ray directions. C KOOR1 must be 1, 2 or 3. C Default: KOOR1=1 C KOOR2=integer... Index of the second coordinate determining the C 2-D section for the calculation of the ray directions. C KOOR2 must be 1, 2 or 3 and must differ from KOOR1. C Default: KOOR2=2 C Data to specify the angular intervals, according which the ray C directions should be sorted: C NA=integer... Number of angular directions. C Default: NA=90 C DA... Angular step. C The default value is recommended. C Default: DA=3.141592/NA C OA... The first angle. The angles are defined in radians, C -pi/2 for positive half-axis KOOR1, 0 for positive C half-axis KOOR2, pi/2 for negative half-axis KOOR1. C The default value is usually sufficient. C Default: OA=-1.570796 C C C Output formatted file MODLEW: C The program calculates lenghts of the ray segments projected C to the plane defined by COOR1 and COOR2. (By the ray segment C we understand the line connecting two successive points in the C file CRT-R.) For each projected ray segment, the slowness vector C at its endpoints is averaged, and the average is considered to be C the direction of the segment. For each angular interval given by C OA, DA and NA the lengths of the ray segments of the direction C corresponding to the interval are accumulated. Finally, the C accumulated lengths are divided by the length af all the ray C segments, multiplied by NA, and the resulting number is written to C the file MODLEW. The file has form C LINes: C (1) String identifying the file, terminated by / (a slash). C (2) The header 2.1 and for each angular interval (NA times) C the values 2.2: C (2.1) 'W' / ... Reference text, terminated by / (a slash). C (2.2) A W / ... Mean angle A of the interval, and the relative length C of the ray segments corresponding to the interval. The data C 2.2 repeat NA times. C (3) / Checksum CHESUM ... CHESUM is the sum of all the lenghts W. C (4) / Average angular change AVDFI ... AVDFI is the average difference C of the directions of the ray along ray segments. 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: C AP00... File 'ap.for'. C C----------------------------------------------------------------------- C CHARACTER*80 FILSEP,FCRT,CH INTEGER LU0 PARAMETER (LU0=1) C C Auxiliary storage locations: INTEGER LU1,LU2,LU3,MA PARAMETER (LU1=1,LU2=2,LU3=3,MA=360) CHARACTER*80 FILE1,FILE2,FILE3 INTEGER KALL,KOOR1,KOOR2,NA,NFI,I1 REAL RAYLEN(MA),FI,CFI,DFI,SUMLEN,SUMDFI,SUMSEG,SEGLEN,AVDFI, * FIOLD,YOLD1,YOLD2,OA,DA,A,CHESUM DATA RAYLEN/MA*0./ C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+CRTLEW: Enter input filename: ' FILSEP=' ' READ(*,*) FILSEP WRITE(*,'(A)') '+CRTLEW: Working ... ' C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LU0,FILSEP) ELSE C CRTLEW-01 CALL ERROR('CRTLEW-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 C Reading input parameters from the SEP file: CALL RSEP3T('CRTOUT',FCRT,' ') FILE1='r01.out' FILE2='r01i.out' IF (FCRT.NE.' ') THEN OPEN (LU1,FILE=FCRT,FORM='FORMATTED',STATUS='OLD') READ (LU1,*) FILE1,CH,FILE2,CH CLOSE (LU1) ENDIF CALL RSEP3T('MODLEW',FILE3,'modlew.out') C C Switch between all rays and only two-point rays: CALL RSEP3I('KALL',KALL,0) C Coordinates of the section: CALL RSEP3I('KOOR1',KOOR1,1) CALL RSEP3I('KOOR2',KOOR2,2) IF (KOOR1.EQ.KOOR2) THEN C CRTLEW-02 CALL ERROR('CRTLEW-02: KOOR1 equal to KOOR2') C KOOR2 must differ from KOOR1. ENDIF C Angular intervals: CALL RSEP3I('NA',NA,90) IF(NA.GT.MA) THEN C CRTLEW-03 CALL ERROR('CRTLEW-03: Too many directions') C The dimension MA of the array RAYLEN should be increased. END IF CALL RSEP3R('DA',DA,3.141592/FLOAT(NA)) CALL RSEP3R('OA',OA,-1.570796+0.5*DA) IF(OA.LT.-1.570796.OR.OA+DA*FLOAT(NA-1).GT.1.570796) THEN C CRTLEW-04 CALL ERROR('CRTLEW-04: Wrong direction') END IF C OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=FILE3) WRITE(LU3,'(9(A,F8.5))') '''Directions and relative weights'' /' C C....................................................................... C C Loop for the points of rays SUMDFI=0. SUMSEG=0. 10 CONTINUE CALL AP00(0,LU1,LU2) IF (IWAVE.LT.1) THEN C End of rays: GOTO 80 ELSEIF (KALL.GT.0.OR.IREC.GT.0) THEN C This ray is to be processed: C Calculating the direction C at the point of the ray: IF (Y(5+KOOR2).EQ.0.) THEN IF (Y(5+KOOR1).LT.0.) THEN FI=1.570796 ELSE FI=-1.570796 ENDIF ELSE FI=ATAN(-Y(5+KOOR1)/Y(5+KOOR2)) ENDIF IF (IPT.GT.1) THEN C Length of the ray segment: SEGLEN=SQRT((YOLD1-Y(2+KOOR1))**2 + (YOLD2-Y(2+KOOR2))**2) C Angular change along the ray segment and the direction C in the center of the segment: DFI=ABS(FI-FIOLD) CFI=(FI+FIOLD)/2. IF (DFI.GT.1.570796) THEN DFI=3.141592-DFI IF (CFI.LT.0.) THEN CFI=1.570796+CFI ELSE CFI=-1.570796+CFI ENDIF ENDIF SUMDFI=SUMDFI+DFI SUMSEG=SUMSEG+1. C Number of the interval for the center of the ray segment: NFI=NINT((CFI-OA)/DA)+1 IF (NFI.EQ.NA+1) NFI=NA IF (NFI.EQ.0) NFI=1 RAYLEN(NFI)=RAYLEN(NFI)+SEGLEN SUMLEN=SUMLEN+SEGLEN ENDIF C Storing the quantities at the point of the ray: FIOLD=FI YOLD1=Y(2+KOOR1) YOLD2=Y(2+KOOR2) ENDIF GOTO 10 C 80 CONTINUE IF ((SUMSEG.EQ.0.).OR.(SUMLEN.EQ.0.)) GOTO 100 AVDFI=SUMDFI/SUMSEG WRITE(LU3,'(9(A,F8.5))') '''W'' /' CHESUM=0. DO 90, I1=1,NA RAYLEN(I1)=RAYLEN(I1)/SUMLEN*FLOAT(NA) A=OA+DA*FLOAT(I1-1) WRITE(LU3,'(9(A,F8.5))') ' ',A,' ',RAYLEN(I1),' /' CHESUM=CHESUM+RAYLEN(I1) 90 CONTINUE C 100 CONTINUE WRITE(LU3,'(9(A,F8.5))') '/ Checksum ',CHESUM WRITE(LU3,'(9(A,F8.5))') '/ Average angular change ',AVDFI,' rad,' WRITE(LU3,'(9(A,F8.3))') '/ i.e. ',180*AVDFI/3.141592,' deg.' CLOSE(LU3) CLOSE(LU2) CLOSE(LU1) WRITE(*,'(A)') '+CRTLEW: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'length.for' C length.for INCLUDE 'ap.for' C ap.for C C======================================================================= C