C
C Subroutine file 'rp2d.for' to control 1-parametric shooting. C C Date: 1998, October 21 C Coded by Ludek Klimes C C....................................................................... C C This file consists of the following subroutines: C RP2D... Auxiliary subroutine to RPAR2, designed to control C one-parametric boundary-value ray tracing (as well as the C initial-value ray tracing). C RP2D C RP2DM...Auxiliary subroutine to RPAR4 designed to store the C quantities dependent on the ray history, required for C one-parametric boundary-value ray tracing. C RP2DM C RP2DG...Auxiliary subroutine to RPAR4 designed to determine the C area GG and inverse specific moment G11,G12,G22 of the C element of the ray-parameter surface, corresponding C to the ray at one-parametric boundary-value ray tracing. C RP2DG C C....................................................................... C C Storage in the memory: C Common block /RP2DC/ to store the quantities needed for C one-parametric boundary-value ray tracing is defined in the C include file 'rp2d.inc'. C rp2d.inc C C======================================================================= C C C SUBROUTINE RP2D(IRAY0,JRAY,G1NEW) INTEGER IRAY0,JRAY REAL G1NEW C C This subroutine controls one-parametric boundary-value ray tracing. C Auxiliary subroutine to RPAR2. C C Input: C IRAY0...Difference IRAY0=IRAY-JRAY between index of a ray within C the elementary wave, and the index within one shooting C domain. C JRAY... Number of the already calculated rays within the shooting C domain. JRAY=0 when starting to cover a new shooting C domain (e.g. when beginning the computation of a new C elementary wave. Otherwise, the output from the previous C invocation of this subroutine. C C Output: C JRAY... JRAY=0 if all rays are computed and the computation of the C elementary wave will be terminated. C Otherwise, input value increased by 1. C G1NEW...Shooting parameter of the ray to be calculated. C The domain of the shooting parameter is specified by the C two first executive statements of this subroutine. C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RP2DC/: INCLUDE 'rp2d.inc' C rp2d.inc C All storage locations of the common block may be altered. C C Subroutines and external functions required: EXTERNAL LUWARN,WRITA INTEGER LUWARN C WRITA... File 'writ.for'. C LUWARN.. File 'crt.for'. C C Date: 1997, November 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,IAUX REAL G1MIN,G1MAX,G1STEP,XAUX CHARACTER*80 TEXT C C G1MIN,G1MAX... Boundaries of the shooting domain (domain of the C ray parameter selected to control the shooting). C G1STEP... Basic step in the shooting parameter. C C Specification of shooting domain and basic step: G1MIN=0. G1MAX=ANUM G1STEP=1. C C....................................................................... C IF(JRAY.EQ.0) THEN NEW=1 A (NEW)=G1MIN-G1STEP IND(NEW)=0 INDOLD =-1 END IF C C Control centre of the one-parametric boundary-value ray tracing 10 CONTINUE IF(A(1).GT.G1MIN+G1STEP/2.) THEN C At least 2 rays of the elementary wave are computed - tests may C be performed: IF(IND(NEW).NE.INDOLD) THEN C Rays with different histories - looking for the boundary ray: IF(NEW.LT.NEWMAX.AND.ABS(A(NEW)-AOLD).GT.AERR) THEN A(NEW+1)=(A(NEW)+AOLD)/2. IRE(NEW+1)=0 ITER=0 GO TO 80 END IF ELSE IF(INDOLD.GT.0) THEN IF(ISRFX(1).NE.0) THEN C Two successful rays - looking for the profiles in the C interval: XAUX=X(NEW) XLEFT=XOLD IF(IREOLD.NE.0) THEN IF((XREC(1,IREOLD)-X(NEW))*(XREC(1,IREOLD)-XOLD).LE.0.) * THEN XLEFT=XREC(1,IREOLD) END IF END IF IAUX=0 C Loop for the given profiles: DO 20 I=1,NREC IF((XREC(1,I)-X(NEW))*(XREC(1,I)-XLEFT).LE.0.) THEN C There is a profile in the interval: boundary-value ray C tracing C5.10 IF(ABS(XREC(1,I)-XOLD).LE.XERR) THEN IF(IREOLD.EQ.I) THEN C Boundary-value ray successfully found - endpoint XOLD ITER=0 ELSE IF(NEW.GE.NEWMAX) THEN C 661 WRITE(TEXT,'(A,I4,A)') * 'Warning 661 in RPAR2: Boundary-value ray to', * I,'th receiver probably inaccurate' CALL WARN(TEXT) C The dimension NEWMAX of the arrays in the common block C /RP2DC/ is too small to allow for sufficiently fine C division of the basic step in the ray parameters. C Boundary-value ray is thus found with greater error C than XERR. ITER=0 ELSE IF((XREC(1,I)-XAUX)*(XREC(1,I)-XLEFT).LE.0.) THEN C There is a boundary-value ray to be found XAUX=XREC(1,I) IAUX=I END IF END IF END IF 20 CONTINUE C5.00 IF(XAUX.NE.X(NEW)) THEN IF(IAUX.NE.0) THEN C There is a boundary-value ray to be found C5.00 IF(ABS(XAUX-X(NEW)).LE.XERR) THEN IF(IRE(NEW).EQ.IAUX) THEN C Boundary-value ray successfully found ITER=0 ELSE C Shooting a new ray within the interval IRE(NEW+1)=0 IF(MOD(ITER,2).EQ.0) THEN C Regula falsi: A(NEW+1)=((XAUX-XOLD)*A(NEW)+(X(NEW)-XAUX)*AOLD) * /(X(NEW)-XOLD) ELSE C Newton-Raphson (paraxial approximation): IF(ABS(XAUX-XOLD).LE.ABS(XAUX-X(NEW))) THEN A(NEW+1)=AOLD+(XAUX-XOLD)/XAOLD ELSE A(NEW+1)=A(NEW)+(XAUX-X(NEW))/XA(NEW) END IF END IF IF(A(NEW+1).LE.AOLD.OR.A(NEW).LE.A(NEW+1)) THEN A(NEW+1)=(A(NEW)+AOLD)/2. END IF IF(A(NEW+1).LE.AOLD.OR.A(NEW).LE.A(NEW+1)) THEN C 662 WRITE(LUWARN(0),'(A,I4,A)') * 'Warning 662 in RPAR2: Boundary-value ray to', * IAUX,'th receiver inaccurate' C Rounding error does not allow for sufficiently fine C division of the basic step in the ray parameters. C Numerically, there is no ray between the two rays C surrounding the profile. C It is recommended to decrease the allowed error UEB of C the computation of the ray (see line (3) of the input C data set dcrt for the complete ray tracing described C in the file 'ray.for'), or to increase the error XERR C (line (2) of the input data for this file) of the C boundary-value ray. IF(ABS(XAUX-XOLD).LE.ABS(XAUX-X(NEW))) THEN A(NEW+1)=AOLD ELSE A(NEW+1)=A(NEW) END IF IRE(NEW+1)=IAUX C-5.10 ITER=0 C-5.10 GO TO 70 END IF ITER=ITER+1 GO TO 80 END IF END IF END IF END IF END IF 70 CONTINUE C Tests are performed - changing to the next interval: C Note: now, DAOLD is an irremovable discrepancy in the element of C the take-off parameter A, corresponding to the old ray. We are C going to forget it. IF(IND(NEW).EQ.INDOLD) THEN C Leaving a homogeneous interval: CALL WRITA(IRAY0+IRAOLD,IRAY0+IRA(NEW),0) END IF INDOLD=IND(NEW) IREOLD=IRE(NEW) IRAOLD=IRA(NEW) AOLD=A(NEW) DAOLD=DA(NEW) XOLD=X(NEW) XAOLD=XA(NEW) NEW=NEW-1 IF(NEW.EQ.0) THEN C New basic interval: A(1)=AOLD+G1STEP IRE(1)=0 ITER=0 IF(A(1).GT.G1MAX+0.001*G1STEP) THEN C End of one-parametric shooting JRAY=0 RETURN END IF GO TO 80 END IF GO TO 10 C C New ray 80 CONTINUE JRAY=JRAY+1 NEW=NEW+1 IRA(NEW)=JRAY G1NEW=A(NEW) RETURN END C C======================================================================= C C C SUBROUTINE RP2DM(IRAY,ISHEET,X1,X1A,IREC) INTEGER IRAY,ISHEET,IREC REAL X1,X1A C C This subroutine stores the quantities dependent on the ray history and C required for one-parametric boundary-value ray tracing. C Designed to be called by subroutine RPAR4. C C Input: C IRAY... Index of the ray within a single shooting domain. C ISHEET..Ray-history index: positive - successful ray, C otherwise - unsuccessful ray. C X1... Value of the first X-function at the point of intersection C of the ray with the reference surface ISRFR. C Meaningful only for successful rays, i.e. ray crossing the C reference surface ISRFR at least IABS(IPOINT) times. C X1A... Derivative of the first X-function with respect to the C shooting parameter A. C Meaningful only for successful rays. C C Output: C IREC... Index of the receiver for a two-point ray, C otherwise zero. C C Common block /RPARD/ (defined in file 'rpar.for'): INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RP2DC/: INCLUDE 'rp2d.inc' C rp2d.inc C Storage locations IND,X,XA of the common block may be altered. C C No subroutines and external functions required. C C Date: 1997, November 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I REAL XAUX C C....................................................................... C IF(IRE(NEW).NE.0) THEN IREC=IRE(NEW) ELSE C Determination of two-point rays: IREC=0 IF(ISRFX(1).NE.0) THEN IF(ISHEET.GT.0) THEN XAUX=XERR C Loop for the given profiles: DO 20 I=1,NREC IF((X1-XOLD)*(XREC(1,I)-XOLD).GE.0.) THEN IF(ABS(XREC(1,I)-X1).LE.XAUX) THEN C Ray endpoint is close to the profile IF(IRAY.GT.1) THEN IF(INDOLD.EQ.ISHEET) THEN IF(IREOLD.EQ.I) THEN GO TO 10 END IF END IF END IF IF(NEW.GT.1) THEN IF(IND(NEW-1).EQ.ISHEET) THEN IF(IRE(NEW-1).EQ.I) THEN GO TO 10 END IF END IF IF((X1-X(NEW-1))*(XREC(1,I)-X(NEW-1)).LT.0.) THEN GO TO 10 END IF END IF IREC=I XAUX=ABS(XREC(1,I)-X1) END IF END IF 10 CONTINUE 20 CONTINUE END IF END IF END IF C C Storing: IND(NEW)=ISHEET IRE(NEW)=IREC X(NEW)=X1 XA(NEW)=X1A C RETURN END C C======================================================================= C C C SUBROUTINE RP2DG(GG,GG11,GG12,GG22) REAL GG,GG11,GG12,GG22 C C This subroutine determines the area GG and inverse specific moment C G11,G12,G22 of the element of the ray-parameter surface, corresponding C to the ray at one-parametric boundary-value ray tracing. C Auxiliary subroutine to RPAR4. C C No input. C C Output: C GG... Area of the element of the ray-parameter surface, C corresponding to the ray. C GG11,GG12,GG22... Components of the symmetric matrix inverse to C the specific moment of the element of the ray-parameter C surface corresponding to the ray, see C.R.T.6.1. C C Common block /RPARD/ (defined in file 'rpar.for'): INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RP2DC/: INCLUDE 'rp2d.inc' C rp2d.inc C Storage locations DAOLD,DA of the common block may be altered. C C No subroutines and external functions required. C C Date: 1994, January 3 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL A1,A2,B1,B2,AB C C A1,A2...Step along the ray-parameter surface in the direction A. C B1,B2...Step along the ray-parameter surface in the direction B. C AB... Area of the basic element of the ray-parameter surface for C a two-parametric system of rays, or the square of its C width for a one-parametric system of rays. The basic C element of the ray-parameter surface is given by the basic C steps in the directions A and B. C C....................................................................... C C Determination of the element of the parameter a corresponding to C the ray: IF(NEW.EQ.1) THEN DA(NEW)=1. ELSE DA(NEW)=(A(NEW-1)-AOLD)/2. DA(NEW-1)=DA(NEW-1)-(A(NEW)-AOLD)/2. DAOLD=DAOLD-(A(NEW-1)-A(NEW))/2. IF(IND(NEW).EQ.IND(NEW-1)) THEN DA(NEW)=DA(NEW)+DA(NEW-1) DA(NEW-1)=0. END IF END IF IF(IND(NEW).EQ.INDOLD) THEN DA(NEW)=DA(NEW)+DAOLD DAOLD=0. END IF C C Determination of the area GG and inverse specific moment G11, G12, C G22 of the element of the ray-parameter surface, corresponding to C the ray: IF(ANUM.NE.0.) THEN A1=(PAR1A-PAR1L)/ANUM A2=(PAR2A-PAR2L)/ANUM IF(BNUM.NE.0.) THEN C Two-parametric system of rays B1=(PAR1B-PAR1L)/BNUM B2=(PAR2B-PAR2L)/BNUM AB=A1*B2-A2*B1 GG =DA(NEW)*ABS(AB) GG11= 12.*(A2*A2+B2*B2)/(AB*AB) GG12=-12.*(A1*A2+B1*B2)/(AB*AB) GG22= 12.*(A1*A1+B1*B1)/(AB*AB) ELSE C One-parametric system of rays AB=A1*A1+A2*A2 GG =DA(NEW)*SQRT(AB) GG11=12.*(A1*A1)/(AB*AB) GG12=12.*(A1*A2)/(AB*AB) GG22=12.*(A2*A2)/(AB*AB) END IF ELSE IF(BNUM.NE.0.) THEN C One-parametric system of rays B1=(PAR1B-PAR1L)/BNUM B2=(PAR2B-PAR2L)/BNUM AB=B1*B1+B2*B2 GG =SQRT(AB) GG11=12.*(B1*B1)/(AB*AB) GG12=12.*(B1*B2)/(AB*AB) GG22=12.*(B2*B2)/(AB*AB) ELSE C Single ray GG =1. GG11=0. GG12=0. GG22=0. END IF END IF DA(NEW)=0. RETURN END C C======================================================================= C