PROGRAM WKBZ IMPLICIT REAL*8(A-H, O-Z) PARAMETER ( NIU = 1, NPU = 6, NOU = 2, & MAXLYR = 100, MAXLYR1 = 101, MAXSSP = 10, & MAXMODE = 2000, MAXNRC = 1000, & MAXSOU = 2, MAXREC = 4 ) REAL*4 RANGE00(MAXSSP), DEPTH00(MAXSSP), & Z00(MAXSSP, MAXLYR1), C00(MAXSSP, MAXLYR1) INTEGER NLAYER00(MAXSSP) REAL*8 ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR), & DEPSOU(MAXSOU), DEPREC(MAXREC), RANREC(MAXNRC), & EIGVL1(MAXMODE), EIGVL2(MAXMODE), UNRMT(MAXREC), & UNRM1(MAXMODE,MAXSOU), UNRM2(MAXMODE,MAXREC) COMPLEX*16 YFACTOR, YY(MAXSOU,MAXREC,MAXNRC) REAL*4 RTL, YYTL, TLS(MAXSOU, MAXREC) CHARACTER*80 TITLE COMMON NLAYER COMMON /EIGEN/ PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS, & WFREQ, WFREQ2, WKZ0, WKSUR2 PI = DATAN( 1.D0 ) * 4. PREC = 1.D-8 CCCCCC PART 1 CCCCCC OPEN(UNIT=NIU, FILE='WKBZ.IN', FORM='FORMATTED', STATUS='OLD') OPEN(UNIT=NOU, FILE='WKBZ.OUT', FORM='FORMATTED', STATUS='NEW') READ(NIU, *) TITLE WRITE(NPU, *) TITLE READ(NIU, *) FREQ, RMKM, DRKM WFREQ = FREQ * PI * 2. WFREQ2 = WFREQ * WFREQ NRCALC = RMKM / DRKM + 1./2. IF ( NRCALC .GT. MAXNRC ) THEN WRITE(NPU, 1002) MAXNRC, NRCALC STOP ENDIF RMKM = DRKM * NRCALC DR = DRKM * 1000. DO 105 I = 1, NRCALC RANREC(I) = DR * I 105 CONTINUE WRITE(NPU,1000) FREQ WRITE(NPU,1001) RMKM READ(NIU, *) NSOU IF ( NSOU .GT. MAXSOU ) THEN WRITE(NPU, 1003) MAXSOU, NSOU STOP ENDIF READ(NIU, *) (DEPSOU(I), I=1,NSOU) READ(NIU, *) NREC IF ( NREC .GT. MAXREC ) THEN WRITE(NPU, 1004) MAXREC, NREC STOP ENDIF READ(NIU, *) (DEPREC(I), I=1,NREC) 1000 FORMAT( 4X, 'Frequency = ', G11.4, 'Hz', / ) 1001 FORMAT( 4X, 'MAX. RANGE = ', G11.4, 'km', // ) 1002 FORMAT( 4X, 'MAXNRC =', I5, 4X, 'NRCALC = ', I5, / ) 1003 FORMAT( 4X, 'MAXSOU =', I5, 4X, 'NSOU = ', I5, / ) 1004 FORMAT( 4X, 'MAXREC =', I5, 4X, 'NREC = ', I5, / ) CCCCCC PART 2 CCCCCC NSSP = 0 110 CONTINUE READ(NIU, *, END=120) RANGE1, DPH11 WRITE(NPU, 1010) RANGE1, DPH11 1010 FORMAT(8X, 'SSP RANGE =',G11.4, 'km, ', 4X, 'DEPTH =',G11.4, 'm',/) NSSP = NSSP + 1 RANGE00(NSSP) = RANGE1 * 1000. DEPTH00(NSSP) = DPH11 DO 115 I = 1, MAXLYR1 READ(NIU, *) ZTMP, CTMP IF ( DABS( ZTMP - DPH11 ) .LT. PREC ) THEN Z00(NSSP,I) = DPH11 C00(NSSP,I) = CTMP NLAYER00(NSSP) = I - 1 GOTO 116 ELSE Z00(NSSP,I) = ZTMP C00(NSSP,I) = CTMP ENDIF 115 CONTINUE 116 CONTINUE IF( RANGE1 .LT. RMKM ) GOTO 110 120 CONTINUE CLOSE(NIU) CCCCCC PART 3 CCCCCC 1020 FORMAT( /4X, 'ISSP = ', I4, 4X, 'RANGE = ', G11.4, 'km', // ) NMODE = MAXMODE ISSP = 1 RANGE1 = 0. RANGE2 = 0. WRITE(NPU,1020) ISSP, RANGE2/1000. DPH11 = DEPTH00(ISSP) NLAYER = NLAYER00(ISSP) NPLUS1 = NLAYER + 1 DO 125 I = 1, NPLUS1 ZI( I ) = Z00(ISSP,I) CI0(I) = C00(ISSP,I) 125 CONTINUE LZ0AXIS = 1 DO 129 I = 1, NLAYER I1 = I + 1 IF( CI0(I1) .EQ. CI0(I) ) CI0(I1) = CI0(I1) - 0.01 AI(I) = ( CI0(I1) - CI0(I) )/( ZI(I1) - ZI(I) ) IF ( CI0(I1) .LT. CI0(LZ0AXIS) ) LZ0AXIS = I1 129 CONTINUE Z0AXIS = ZI(LZ0AXIS) !!!! CHANNEL AXIS CZ0AXIS = CI0(LZ0AXIS) WKZ0 = WFREQ / CZ0AXIS WKSUR2 = WFREQ2 / ( CI0(1) * CI0(1) ) IF( LZ0AXIS .GT. 1 ) THEN ULIN = WFREQ/CI0( LZ0AXIS-1 ) ELSE ULIN = WFREQ/CI0( CI0(1) + CI0(2) ) * 2. ENDIF LTHVAL = 0 130 CONTINUE LTHVAL1 = LTHVAL + 1 CALL EIGENVAL( LTHVAL, ULIN, ULOUT, XSP, PHADJ, & ZREA, LREA, ZREB, LREB, MMSUR, MMBOT ) IF ( MMBOT .NE. 0 ) THEN NMODE = LTHVAL GOTO 140 ENDIF EIGVL1(LTHVAL1) = 0. EIGVL2(LTHVAL1) = ULOUT write(*,*) lthval1, ulout, wfreq/ulout DO 135 I = 1, NSOU ZDPH = DEPSOU(I) CALL EIGENFUN( LTHVAL, ULOUT, XSP, PHADJ, & ZREA, LREA, ZREB, LREB, ZDPH, YL ) UNRM1( LTHVAL1, I) = YL 135 CONTINUE DO 139 I = 1, NREC ZDPH = DEPREC(I) CALL EIGENFUN( LTHVAL, ULOUT, XSP, PHADJ, & ZREA, LREA, ZREB, LREB, ZDPH, YL ) UNRM2( LTHVAL1, I) = YL 139 CONTINUE LTHVAL = LTHVAL1 IF( LTHVAL .LT. NMODE ) THEN ULIN = ULOUT - PI * 2. /XSP GOTO 130 ENDIF CCCCCC PART 4 CCCCCC 140 CONTINUE NEEDBG = 1 NEEDLF = NRCALC - NEEDBG + 1 150 CONTINUE IF ( ISSP .EQ. NSSP ) GOTO 200 CCCCCC PART 5 CCCCCC ISSP = ISSP + 1 RANGE1 = RANGE2 RANGE2 = RANGE00(ISSP) WRITE(NPU,1020) ISSP, RANGE2/1000. DPH11 = DEPTH00(ISSP) NLAYER = NLAYER00(ISSP) NPLUS1 = NLAYER + 1 INR = NEEDBG NEEDCA = 0 161 CONTINUE IF ( RANREC( INR ) .GT. RANGE2 ) GOTO 162 INR = INR + 1 NEEDCA = NEEDCA + 1 IF ( INR .LE. NRCALC ) THEN GOTO 161 ELSE NEEDLF = 0 GOTO 163 ENDIF 162 NEEDLF = NRCALC - INR + 1 163 CONTINUE DO 165 I = 1, NPLUS1 ZI( I ) = Z00(ISSP,I) CI0(I) = C00(ISSP,I) 165 CONTINUE LZ0AXIS = 1 DO 169 I = 1, NLAYER I1 = I + 1 IF( CI0(I1) .EQ. CI0(I) ) CI0(I1) = CI0(I1) - 0.01 AI(I) = ( CI0(I1) - CI0(I) )/( ZI(I1) - ZI(I) ) IF ( CI0(I1) .LT. CI0(LZ0AXIS) ) LZ0AXIS = I1 169 CONTINUE Z0AXIS = ZI(LZ0AXIS) !!!! CHANNEL AXIS CZ0AXIS = CI0(LZ0AXIS) WKZ0 = WFREQ / CZ0AXIS IF( LZ0AXIS .GT. 1 ) THEN ULIN = WFREQ/CI0( LZ0AXIS-1 ) ELSE ULIN = WFREQ/CI0( CI0(1) + CI0(2) ) * 2. ENDIF LTHVAL = 0 170 CONTINUE LTHVAL1 = LTHVAL + 1 CALL EIGENVAL( LTHVAL, ULIN, ULOUT, XSP, PHADJ, & ZREA, LREA, ZREB, LREB, MMSUR, MMBOT ) IF ( MMBOT .NE. 0 ) THEN NMODE = LTHVAL GOTO 190 ENDIF write(*,*) lthval1, ulout, wfreq/ulout DO 172 I = 1, NREC ZDPH = DEPREC(I) CALL EIGENFUN( LTHVAL, ULOUT, XSP, PHADJ, & ZREA, LREA, ZREB, LREB, ZDPH, YL ) UNRMT(I) = YL 172 CONTINUE IF ( NEEDCA .EQ. 0 ) GOTO 180 DO 179 INR = 1, NEEDCA INRR = INR + NEEDBG - 1 RANRR = RANREC(INRR) - RANGE1 RANRR1 = RANRR / ( RANGE2 - RANGE1 ) ULOUTR = EIGVL2(LTHVAL1) + ( ULOUT - EIGVL2(LTHVAL1) ) * RANRR1 WNPHA = EIGVL1(LTHVAL1) + (EIGVL2(LTHVAL1)+ULOUTR) * RANRR / 2. YFACTOR = CDEXP( DCMPLX(0.D0, WNPHA) ) / DSQRT( WNPHA ) DO 175 IR = 1, NREC YLR = UNRM2(LTHVAL1, IR) + & ( UNRMT(IR)-UNRM2(LTHVAL1, IR) )*RANRR1 DO 174 IS = 1, NSOU YRINRR = UNRM1(LTHVAL1, IS) * YLR YY(IS, IR, INRR) = YY(IS, IR, INRR) + YFACTOR * YRINRR 174 CONTINUE 175 CONTINUE 179 CONTINUE 180 CONTINUE EIGVL1(LTHVAL1) = EIGVL1(LTHVAL1) + & ( EIGVL2(LTHVAL1) + ULOUT ) * (RANGE2-RANGE1) / 2. EIGVL2(LTHVAL1) = ULOUT DO 182 I = 1, NREC UNRM2(LTHVAL1, I) = UNRMT(I) 182 CONTINUE LTHVAL = LTHVAL1 IF( LTHVAL .LT. NMODE ) THEN ULIN = ULOUT - PI * 2. /XSP GOTO 170 ENDIF 190 CONTINUE IF ( NEEDLF .EQ. 0 ) THEN GOTO 300 ELSE NEEDBG = NRCALC - NEEDLF + 1 GOTO 150 ENDIF CCCCCC PART 6 CCCCCC 200 CONTINUE RANGE1 = RANGE2 DO 210 INRR = NEEDBG, NRCALC RANRR = RANREC( INRR ) - RANGE1 DO 209 LTHVAL = 1, NMODE WNPHA = EIGVL1(LTHVAL) + EIGVL2(LTHVAL) * RANRR YFACTOR = CDEXP( DCMPLX(0.D0, WNPHA) ) / DSQRT( WNPHA ) DO 205 IR = 1, NREC YLR = UNRM2(LTHVAL, IR) DO 204 IS = 1, NSOU YRINRR = UNRM1(LTHVAL, IS) * YLR YY(IS, IR, INRR) = YY(IS, IR, INRR) + YFACTOR * YRINRR 204 CONTINUE 205 CONTINUE 209 CONTINUE 210 CONTINUE CCCCCC PART 7 CCCCCC 300 CONTINUE WRITE(NOU,*) TITLE WRITE(NOU, 1000) FREQ WRITE(NOU, 2000) NSOU WRITE(NOU, 2001) (DEPSOU(IS), IS=1, NSOU) WRITE(NOU, 2002) NREC WRITE(NOU, 2003) (DEPREC(IR), IR=1, NREC) 2000 FORMAT(/4X, 'SOURCE', 6X, I5, 6x, 'DEPTH IN METERS') 2001 FORMAT(/4X, 4G11.4) 2002 FORMAT(/4X, 'RECEIVER', 4X, I5, 6x, 'DEPTH IN METERS') 2003 FORMAT(/4X, 4G11.4) WRITE(NOU, 2004) 2004 FORMAT(//X, 'RANGE (km) & TRANSMISSION LOSS (dB)', //) FACTOR = 20. / LOG(10.) YFACTOR = CDEXP( DCMPLX(0.D0, PI/4.) ) * DSQRT( PI * 2. ) IF ( (NSOU + NREC) .GT. 2 ) GOTO 320 DO 310 INR = 1, NRCALC RTL = RANREC(INR) / 1000. YY( 1, 1, INR) = YFACTOR * YY( 1, 1, INR ) TLS( 1, 1 ) = FACTOR * LOG( CDABS( YY( 1, 1, INR ) ) ) YYTL = CDABS( YY( 1, 1, INR ) ) IF ( YYTL .LE. 1.D-10 ) THEN TLS( 1, 1 ) = -200. ELSE TLS( 1, 1 ) = FACTOR * LOG( YYTL ) ENDIF WRITE(NOU, *) RTL, TLS(1,1) 310 CONTINUE GOTO 350 320 CONTINUE DO 329 INR = 1, NRCALC RTL = RANREC(INR) / 1000. DO 329 IS = 1, NSOU DO 328 IR = 1, NREC YY( IS, IR, INR) = YFACTOR * YY( IS, IR, INR ) YYTL = CDABS( YY( IS, IR, INR ) ) IF ( YYTL .LE. 1.D-10 ) THEN TLS( IS, IR ) = -200. ELSE TLS( IS, IR ) = FACTOR * LOG( YYTL ) ENDIF 328 CONTINUE WRITE(NOU,*) RTL, ( TLS(IS,IR), IR=1, NREC ) 329 CONTINUE 350 CLOSE(NOU) STOP END SUBROUTINE EIGENVAL( LTHVAL, ULIN, ULOUT, XSP, PHADJSUR, & ZREA, LREA, ZREB, LREB, MMSUR, MMBOT) IMPLICIT REAL*8(A-H, O-Z) PARAMETER ( MAXLYR = 100, MAXLYR1 = 101, PREC = 1.D-8 ) REAL*8 ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR) COMMON NLAYER COMMON /EIGEN/ PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS, & WFREQ, WFREQ2, WKZ0, WKSUR2 NITER = 0 UL0 = ULIN 110 NITER = NITER + 1 CALL SPANPHASE( UL0,XSP,XPH,ZREA,LREA,ZREB,LREB,MMSUR,MMBOT ) TVARAI0 = WKSUR2 - UL0 * UL0 IF ( MMSUR .NE. 0 ) THEN OMEGA0 = DABS( TVARAI0 ) OMEGA0 = -( CI0(1)*( OMEGA0**(3./2.) ) )/(AI(1) * WKSUR2 * 3.) ELSE OMEGA0 = 0. ENDIF DGRAD23 = ( DABS( -AI(1) * WKSUR2 * 2. / CI0(1) ) ) ** (2./3.) TVARAI0 = - TVARAI0 / DGRAD23 CALL AIRYFUNC( TVARAI0, THETA ) PHADJSUR = PI / 2. + DATAN( DTAN( ( THETA - OMEGA0 ) * 2. ) ) UL1 = UL0 - ( ( PI*LTHVAL*2. + PI/2. + PHADJSUR ) - XPH ) / XSP IF ( DABS( UL1 - UL0 )/UL1 .GT. PREC ) THEN UL0 = UL1 GOTO 110 ENDIF ULOUT = UL1 RETURN END SUBROUTINE EIGENFUN( LTHVAL, UL, XSP, PHADJSUR, & ZREA, LREA, ZREB, LREB, ZD, YZ) IMPLICIT REAL*8(A-H, O-Z) PARAMETER ( MAXLYR = 100, MAXLYR1 = 101, PREC = 1.D-8 ) REAL*8 ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR) COMMON NLAYER COMMON /EIGEN/ PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS, & WFREQ, WFREQ2, WKZ0, WKSUR2 PARAMETER ( BPARA = 2.152D0, DPARA = 1.619D0 ) C LZD = 0 110 LZD = LZD + 1 IF ( ZD .GT. ZI( LZD + 1 ) ) GOTO 110 CZD = CI0( LZD ) + AI( LZD ) * ( ZD - ZI( LZD ) ) C UKZ = WFREQ / CZD SQUKZ = UKZ * UKZ EKSQZ = DABS( -SQUKZ * AI( LZD ) * 2. / CZD ) EKS23 = EKSQZ ** ( 2./3. ) EKS43 = EKSQZ ** ( 4./3. ) SQK = ( UKZ + UL ) * ( UKZ - UL ) SQK2 = SQK * SQK C IF ( ZD .GE. ZREB ) GOTO 400 IF ( ZD .GE. ZREA ) GOTO 300 CCCCCC [ 0, ZREA ) CCCCCC CONST = WFREQ/UL XCOSA = CONST / CZD XTANA = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA XANGL = DATAN( XTANA ) XXP0 = XTANA - XANGL XPHASE = 0. IF ( LZD .GE. LREA ) GOTO 210 DO II = LZD+1, LREA XCOSA = CONST / CI0( II ) XTANA = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA XANGL = DATAN( XTANA ) XXPP = XTANA - XANGL XPHASE = XPHASE + ( XXPP - XXP0 ) / AI( II-1 ) XXP0 = XXPP END DO 210 XPHASE = XPHASE - XXP0 / AI ( LREA ) XPHASE = WFREQ * XPHASE C FRACT = XPHASE FRACT = DEXP( -FRACT ) * DCOS( PHADJSUR / 2. ) DENOM = BPARA * EKS43 - DPARA * EKS23 * SQK + SQK2 * 16. GOTO 500 CCCCCC [ ZREA, ZREB ) CCCCCC 300 CONTINUE CONST = UL / WFREQ IF (ZREA .EQ. 0. ) THEN XSINA = CI0(1) * CONST XSINA = DSQRT( DABS( 1.D0 - XSINA*XSINA ) ) ELSE XSINA = 0. END IF XXP0 = DLOG ( ( 1.D0 + XSINA )/( 1.D0 - XSINA ) ) - XSINA * 2. XPHASE = 0 IF ( LREA .GE. LZD ) GOTO 310 DO II = LREA+1, LZD XSINA = CI0( II ) * CONST XSINA = DSQRT( DABS( 1.D0 - XSINA*XSINA ) ) XXPP = DLOG( ( 1.D0 + XSINA )/( 1.D0 - XSINA ) ) - XSINA * 2. XPHASE = XPHASE + ( XXP0-XXPP ) / AI( II-1 ) XXP0 = XXPP END DO 310 XSINA = CZD * CONST XSINA = DSQRT( DABS( 1.D0 - XSINA * XSINA ) ) XXPP = DLOG( ( 1.D0 + XSINA )/( 1.D0 - XSINA ) ) - XSINA * 2. XPHASE = XPHASE + ( XXP0-XXPP ) / AI(LZD) XPHASE = WFREQ * XPHASE / 2. C FRACT = XPHASE + ( PI - PHADJSUR ) / 2. FRACT = DSIN( FRACT ) DENOM = BPARA * EKS43 - DPARA * EKS23 * SQK + SQK2 GOTO 500 CCCCCC [ZREB, DPH11 ] CCCCCC 400 CONTINUE CONST = WFREQ/UL IF ( ZREB .EQ. DPH11 ) THEN XCOSA = CONST / CI0( NLAYER + 1 ) XTANA = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA XANGL = DATAN( XTANA ) ELSE XTANA=0. XANGL=0. ENDIF XXP0 = XTANA - XANGL XPHASE = 0 IF ( LREB .GE. LZD ) GOTO 410 DO II = LREB+1, LZD XCOSA = CONST / CI0( II ) XTANA = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA XANGL = DATAN( XTANA ) XXPP = XTANA - XANGL XPHASE = XPHASE + DABS( ( XXPP - XXP0 ) / AI( II-1 ) ) XXP0 = XXPP END DO 410 XCOSA = CONST / CZD XTANA = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA XANGL = DATAN( XTANA ) XXPP = XTANA - XANGL XPHASE = XPHASE + DABS( ( XXPP - XXP0 ) / AI( LZD ) ) XPHASE = WFREQ * XPHASE C FRACT = DFLOAT((-1)**LTHVAL) * DEXP( -XPHASE ) / DSQRT( 2.D0 ) DENOM = BPARA * EKS43 - DPARA * EKS23 * SQK + SQK2 * 16. CCCCCC 500 FRACTOR = DSQRT( UL * 4. / XSP ) YZ = FRACTOR * FRACT / DENOM ** ( 1./8. ) RETURN END SUBROUTINE SPANPHASE( UL, XSPAN, XPHASE, ZREA, LREA, & ZREB, LREB, MMSUR, MMBOT ) IMPLICIT REAL*8(A-H, O-Z) PARAMETER ( MAXLYR = 100, MAXLYR1 = 101 ) REAL*8 ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR) REAL*8 XSINA(MAXLYR1) COMMON NLAYER COMMON /EIGEN/ PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS, & WFREQ, WFREQ2, WKZ0, WKSUR2 C MMSUR=0 MMBOT=0 C CALCULATE THE SOUND SPEED AT THE TURNING-POINT DEPTH CTURN = WFREQ/UL C ZREA -- THE RAY TURNING-DEPTH ABOVE THE CHANNEL AXIS LREA = LZ0AXIS - 1 110 CONTINUE IF ( LREA .EQ. 0 ) THEN MMSUR = 1 ZREA = 0. GOTO 120 ENDIF IF ( CTURN .LT. CI0( LREA ) ) THEN GOTO 111 ELSE LREA = LREA - 1 GOTO 110 ENDIF 111 ZREA = ZI( LREA ) + ( CTURN - CI0( LREA ) ) / AI( LREA ) 120 CONTINUE C ZREB -- THE RAY TURNING-DEPTH ABOVE THE CHANNEL AXIS LREB = LZ0AXIS 130 CONTINUE IF ( LREB .GT. NLAYER ) THEN MMBOT = 1 ZREB = DPH GOTO 140 ENDIF IF ( CTURN .LT. CI0( LREB+1 ) ) THEN GOTO 131 ELSE LREB = LREB + 1 GOTO 130 ENDIF 131 ZREB = ZI( LREB ) + ( CTURN - CI0( LREB ) ) / AI( LREB ) 140 CONTINUE C SOUND RAY SPAN CONST = UL / WFREQ XSPAN = 0. XPHASE = 0. IF ( LREA .EQ. 0) THEN LREA = 1 XXSINA = CI0(1) * CONST XSINA(1) = DSQRT( DABS( 1.D0 - XXSINA * XXSINA ) ) ELSE XSINA(LREA) = 0. ENDIF IF ( LREB .GT. NLAYER ) THEN LREB = NLAYER XXSINA = CI0( NLAYER+1 ) * CONST XSINA(NLAYER+1) = DSQRT( DABS( 1.D0 - XXSINA * XXSINA ) ) ELSE XSINA( LREB+1 ) = 0. END IF XXP0= DLOG( ( 1.D0 + XSINA(LREA) ) / ( 1.D0 - XSINA(LREA) ) ) 1 - XSINA( LREA ) * 2. DO II = LREA+1, LREB XXSINA = CI0( II ) * CONST XSINA(II) = DSQRT( DABS( 1.D0 - XXSINA * XXSINA ) ) XXPP = DLOG( ( 1.D0 + XSINA(II) ) / ( 1.D0 - XSINA(II) ) ) & - XSINA( II ) * 2. XSPAN = XSPAN + ( XSINA( II-1 ) - XSINA( II ) ) / AI( II-1 ) XPHASE = XPHASE + ( XXP0 - XXPP ) / AI( II-1 ) XXP0 = XXPP END DO XSPAN = XSPAN + ( XSINA( LREB ) - XSINA( LREB+1 ) ) / AI( LREB ) XSPAN = XSPAN * 2. / CONST XXPP = DLOG( ( 1.D0 + XSINA(LREB+1) )/( 1.D0 - XSINA(LREB+1) ) ) 1 - XSINA(LREB+1) * 2. XPHASE = XPHASE + ( XXP0 - XXPP ) / AI( LREB ) XPHASE = WFREQ * XPHASE RETURN END SUBROUTINE AIRYFUNC( XTIN, THETA ) IMPLICIT REAL*8(A-H, O-Z) PARAMETER ( PREC = 1.D-6 ) IF ( XTIN .EQ. 0.D0 ) THEN UAIX=1.0899290710 !!XTIN = 0. VAIX=0.6292708425 GOTO 2500 ENDIF C C IF ( ( XTIN .LT. 5. ) .AND. ( XTIN .GT. -8. ) ) GOTO 2460 IF ( XTIN .LT. 0. ) GOTO 2450 C C XW = ( XTIN**(3./2.) )*2./3. !! XTIN > 0. EXW = DEXP( -XW * 2. ) XW72 = 1.D0 / (XW*72.) NAIITER=1 XXAI=DFLOAT(5)*XW72 UAIX=DFLOAT(1)+XXAI VAIX=DFLOAT(1)-XXAI 2401 NAIITER = NAIITER+1 NAIITER6 = 6*NAIITER XXAI = XXAI*XW72*(NAIITER6-1)*(NAIITER6-5)/NAIITER UAIX = UAIX+XXAI VAIX = VAIX+XXAI*((-1)**NAIITER) IF ( DABS(XXAI/UAIX) .GE. PREC ) GOTO 2401 THETA = DATAN( VAIX/UAIX * EXW * 0.5 ) GOTO 2500 C C 2450 XTIN = DABS(XTIN) !! XTIN < 0. XW = (XTIN**(3./2.))*2./3. XWSS = 1.D0 / ( XTIN**(1./4.) ) PI4 = DATAN( 1.D0 ) XWCOS = DCOS( XW+PI4 ) XWSIN = DSIN( XW+PI4 ) XW72 = 1.D0 / ( XW*72. ) NAIITER = 1 XXUAI = 1.D0 XXVAI = XW72 * 5. UAIX = 1.D0 VAIX = XXVAI 2451 NAIITER = NAIITER+1 NAIITER6 = 6*NAIITER XXUAI = -XXVAI*XW72*(NAIITER6-1)*(NAIITER6-5)/NAIITER NAIITER = NAIITER+1 NAIITER6 = 6*NAIITER XXVAI = XXUAI*XW72*(NAIITER6-1)*(NAIITER6-5)/NAIITER UAIX = UAIX + XXUAI VAIX = VAIX + XXVAI IF((DABS(XXUAI/UAIX).GE. PREC ).OR.(DABS(XXVAI/VAIX).GE. PREC )) 1 GOTO 2451 XXUAI = XWCOS*UAIX + XWSIN*VAIX VAIX = XWSIN*UAIX - XWCOS*VAIX UAIX = XXUAI THETA = DATAN( VAIX/UAIX ) GOTO 2500 C 2460 UAIX0 = 1.0899290710 !!!! VAIX0 = 0.6292708425 UAIX10 = 0.7945704238 VAIX10 = -0.4587454481 XTIN3 = XTIN ** 3. XXUAI = 1.D0 XXVAI = 1.D0 UAIX = 1.D0 VAIX = 1.D0 NAIITER = 0 2461 NAIITER = NAIITER + 1 NAI3 = NAIITER * 3 XXUAI = XXUAI * XTIN3 / ( ( NAI3-1 ) * NAI3 ) XXVAI = XXVAI * XTIN3 / ( NAI3 * ( NAI3+1 ) ) UAIX = UAIX + XXUAI VAIX = VAIX + XXVAI IF ( DABS( XXVAI/VAIX ).GT. PREC ) GOTO 2461 XXUAI = UAIX0 * UAIX + UAIX10 * VAIX * XTIN VAIX = VAIX0 * UAIX + VAIX10 * VAIX * XTIN UAIX = XXUAI THETA = DATAN( VAIX/UAIX ) 2500 CONTINUE RETURN END