CTITLESANRAP -- NEWTON - RAPHSHON APPROXIMATION OF THETA-1 00010001 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** C 00020000 CA AUTHOR ARCO 00030000 CA DESIGNER ARCO 00040000 CA LANGUAGE FORTRAN H 00050000 CA WRITTEN 09-01-84 00060000 C REVISED MO-DA-YR BY PROGRAMMER 00070000 C REVISED 04-02-87 JMP - PUT IN PRODUCTION SPARC. 00080001 C 00090000 CA CALLING SEQUENCE: 00100000 CA CALL SANRAP ( OFFSET, KL, D, DEPTH, ALPHA, ANGLE) 00110001 CA 00120000 CA IN/OUT ARGUMENT TYPE DESCRIPTION 00130000 CA 00140000 CA IN OFFSET R4 OFFSET DISTANCE FROM SOURCE TO RECEIVER 00150000 CA CURRENTLY BEING CALCULATED. 00160000 CA IN KL I4 POINTER TO LAYER CURRENTLY BEING WORKED ON 00170000 CA IN/OUT D R4 DEPTH TO ITH LAYER 00180000 CA 00230000 CA PURPOSE OF PROGRAM: 00240000 CA THIS SUBROUTINE IS USED TO CALCULATE THE ANGLE "THETA" AT WHICH 00250000 CA A SEISMIC RAY WILL REFLECT FROM THE INTERFACE IN A GIVEN LAYER 00260000 CA FOR A MODEL WHICH IS DESCRIBED BY A P-WAVE SONIC, S-WAVE SONIC, 00270000 CA AND BULK DENSITY FOR EACH LAYER. THIS APPROXIMATION IS AN 00280000 CA ITERATIVE PROCESS FOR APPROXIMATING THE ANGLE. MOST CALCULATIONS 00290000 CA ARE DONE USING DOUBLE PRECISION VARIABLES TO AID IN CONVERGING ON 00300000 CA A SOLUTION. 00310000 CA 00320000 CA ROUTINES CALLED: DASIN - DOUBLE PRECISION ARCSINE FUNCTION 00330001 CA DATAN - DOUBLE PRECISION ARCTAN FUNCTION 00340001 CA SADIST - FUNCTION FOR COMPUTING DISTANCE 00350001 CA SADRIV - FUNCTION FOR COMPUTING DERIVATIVE 00360001 CA 00370000 C--------------------------------------------------------------------- 00380000 SUBROUTINE SANRAP ( OFFSET, KL, D, DEPTH, ALPHA, ANGLE) 00390001 C 00400000 REAL DEPTH (1) 00650001 C 00690000 REAL*8 ALPHA (1) 00850001 3 , ANGLE (1) 00860001 4 , PI 00870000 C 00890000 REAL *8 DEL 01150000 1 , D 01160000 2 , THETA 01170000 3 , XR 01180000 4 , DXR 01190000 5 , ALPHMX 01200000 6 , THICK 01210000 REAL *4 OFFSET 01230000 C 01240000 INTEGER KL 01340001 C 01350000 C____________________________________________________________ 01360000 C 01370000 C PROGRAM STARTS HERE 01380000 C 01390000 C____________________________________________________________ 01400000 C 01410000 PI = 4.0 * ATAN (1.0) 01411001 RADIAN = 180.0 / PI 01412001 DEL = 0.0 01420000 D = D + DEPTH(KL) 01430000 THICK = DEPTH(KL) 01440000 IF (KL .EQ. 1) GO TO 300 01450000 C 01460000 C**************************************************** 01470000 C 01480000 C FOR MORE THAN ONE LAYER, MAKE A FIRST GUESS AT THE VALUE THETA 01490000 C TO GET THE NEWTON-RAPHSHON ROLLING. 01500000 C 01510000 C**************************************************** 01520000 C 01530000 IF (ALPHA(KL).GE.ALPHA(KL-1)) GO TO 300 01540001 C 01550000 C**************************************************** 01560000 C 01570000 C ONE APPROXIMATION IS IF THE KL-TH VELOCITY IS LESS THAN THE 01580000 C (KL-1)-TH VELOCITY. 01590000 C 01600000 C**************************************************** 01610000 C 01620000 THETA = DASIN (DSIN (0.5 * DATAN (0.5 * OFFSET / D)) /ALPHA(KL)) 01630001 GO TO 400 01640000 C 01650000 C**************************************************** 01660000 C 01670000 C AND THE OTHER APPROXIMATION IS FOR AN INCREASING VELOCITY 01680000 C PROFILE. 01690000 C 01700000 C**************************************************** 01710000 C 01720000 300 THETA = DASIN (DSIN (0.5 * (DATAN (0.5 * OFFSET / D) + 01730001 1 DATAN( 0.5 * OFFSET / THICK))) / ALPHA(KL)) 01740000 C 01750000 C**************************************************** 01760000 C 01770000 C NOW THE ITERATION REALLY STARTS. INITIALIZE THETA TO THE FIRST 01780000 C GUESS WE MADE ABOVE, AND ALTER IT BY THE DEL CALCULATED BELOW. 01790000 C 01800000 C**************************************************** 01810000 C 01820000 400 THETA = THETA + DEL 01830000 C 01840000 C**************************************************** 01850000 C 01860000 C THE SHOT-RECEIVER DISTANCE FUNCTION, XR, AND ITS DERIVATIVE, 01870000 C DXR, MUST BE STARTED AT ZERO FOR EVERY RAY WHOSE THETA WE 01880000 C CALCULATE. 01890000 C 01900000 C**************************************************** 01910000 C 01920000 XR = 0.0 01930000 DXR = 0.0 01940000 C 01950000 C**************************************************** 01960000 C 01970000 C DO THE SUMMATION IN THE XR AND DXR EQUATIONS USING THE SADIST 01980001 C AND SADRIV FUNCTIONS. THUS XR AND DXR ARE CALCULATED FOR I 01990001 C LAYERS. WHEN THETA IS THE ROOT OF XR, THEN XR = OFFSET. THE 02000000 C FUNCTION XR IS THE SUMMATION TO I OF THE TWICE THE THICKNESS 02010000 C OF LAYER J TIMES THE TANGENT OF THE RAY'S ANGLE IN LAYER J. 02020000 C BUT THE ANGLE IN LAYER J IS IN TERMS OF THETA-1, THE RAY'S 02030000 C ANGLE IN LAYER 1. DXR IS THE DERIVATIVE OF THIS FUNCTION. 02040000 C 02050000 C**************************************************** 02060000 C 02070000 DO 500 J = 1, KL 02080000 XR = XR + SADIST ( THETA, ALPHA(J), DEPTH(J)) 02090001 DXR = DXR + SADRIV ( THETA, ALPHA(J), DEPTH(J)) 02100001 500 CONTINUE 02110000 DEL = - (XR - OFFSET) / DXR 02120000 C 02130000 C**************************************************** 02140000 C 02150000 C FIND THE MAXIMUM ELEMENT OF THE ALPHA ARRAY UP TO THE I-TH 02160000 C ELEMENT. 02170000 C 02180000 C**************************************************** 02190000 C 02200000 ALPHMX = ALPHA(1) 02210000 IF (KL .EQ. 1) GO TO 1000 02220000 DO 700 L = 1, KL 02230000 ALPHMX = DMAX1 (ALPHMX, ALPHA(L)) 02240000 700 CONTINUE 02250000 C 02260000 C**************************************************** 02270000 C 02280000 C NOW WE MAKE SURE THAT SIN(THETA+DEL) IS NOT GREATER THAN THE 02290000 C INVERSE OF THE BIGGEST ALPHA THAT WE ARE USING AT THIS MOMENT. 02300000 C THIS IS BECAUSE WE ARE TAKING THE SQUARE ROOT OF (1-(ALPHA(I)* 02310000 C SIN(THETA))**2). IF, INDEED, SIN(THETA) > 1/ALPHA(I) THEN WE 02320000 C WILL BE TRYING TO TAKE THE SQUARE ROOT OF A NEGATIVE NUMBER. 02330000 C BUT, IF WE ARE SUMMING XR AND DXR OVER MORE THAN ONE LAYER, 02340000 C WE WANT TO BE SURE THAT SIN(THETA) IS NOT GREATER THAN THE 02350000 C LARGEST ALPHA(I), NAMELY ALPHMX. 02360000 C 02370000 C**************************************************** 02380000 C 02390000 IF ( DSIN(THETA+DEL) .LT. (1 / ALPHMX)) GO TO 1000 02400001 DEL = DASIN(0.5 * (1 / ALPHMX + DSIN(THETA))) - THETA 02410001 C 02420000 C**************************************************** 02430000 C 02440000 C PLACE THE LIMITATION ON DEL THAT WILL EITHER CONTINUE OR STOP 02450000 C THE NEWTON_RAPHSHON ITERATION. 02460000 C 02470000 C**************************************************** 02480000 C 02490000 1000 IF (DABS(DEL) .GT. 0.000001) GO TO 400 02500000 C 02510000 C**************************************************** 02520000 C 02530000 C MAKE SURE THETA IS NO GREATER THAN PI. 02540000 C 02550000 C**************************************************** 02560000 C 02570000 IF (THETA .GE. PI) THETA = PI * 02580001 1 ((THETA / PI) - (DINT(THETA / PI))) 02590000 C 02600000 C**************************************************** 02610000 C 02620000 C AND, MAKE SURE THETA IS NO GREATER THAN PI/2. 02630000 C 02640000 C**************************************************** 02650000 C 02660000 IF (THETA .GE. (PI / 2.0)) THETA = PI - THETA 02670000 C 02680000 C**************************************************** 02690000 C 02700000 C CONVERT THETA FROM RADIANS TO DEGREES. 02710000 C 02720000 C AND STORE THIS THETA AS THE INCIDENT ANGLE FROM THE SHOT 02730000 C FOR THE RAY REFLECTING OFF THE BOTTOM OF LAYER I. 02740000 C 02750000 C**************************************************** 02760001 C 02770001 ANGLE(KL) = THETA 02790000 C 02800000 C**************************************************** 03480000 C 03490000 C RETURN TO CALLING PROGRAM 03500000 C 03510000 C**************************************************** 03520000 C 03530000 RETURN 03540001 END 03600000