\documentstyle[12pt]{article} \batchmode \textwidth = 210 mm \textheight = 150 mm \topmargin = 0 mm \oddsidemargin = 0 mm \evensidemargin = 0 mm \parskip = 5 mm \pagestyle{empty} \newcommand{\heading}[1]{ \vspace{3 mm} \noindent {\bf {#1}} \vspace{5 mm} } \newcommand{\be}{\[} \newcommand{\ee}{\]} \newcommand{\bea}{\begin{eqnarray*}} \newcommand{\eea}{\end{eqnarray*}} \newcommand{\gee}{\mbox{g}} \newcommand{\bfa}{{\bf a}} \newcommand{\bfad}{{\bf \dot{a}}} \newcommand{\bfadd}{{\bf \ddot{a}}} \newcommand{\bfb}{{\bf b}} \newcommand{\bfc}{{\bf c}} \newcommand{\bfd}{{\bf d}} \newcommand{\bfe}{{\bf e}} \newcommand{\bff}{{\bf f}} \newcommand{\bfg}{{\bf g}} \newcommand{\bfh}{{\bf h}} \newcommand{\bfk}{{\bf k}} \newcommand{\bfkp}{{\bf k}^{\prime}} \newcommand{\bfn}{{\bf n}} \newcommand{\bfp}{{\bf p}} \newcommand{\bfq}{{\bf q}} \newcommand{\bfr}{{\bf r}} \newcommand{\bfs}{{\bf s}} \newcommand{\bft}{{\bf t}} \newcommand{\bfu}{{\bf u}} \newcommand{\bfv}{{\bf v}} \newcommand{\bfvhat}{\widehat{\bf v}} \newcommand{\bfx}{{\bf x}} \newcommand{\bfA}{{\bf A}} \newcommand{\bfB}{{\bf B}} \newcommand{\bfC}{{\bf C}} \newcommand{\bfD}{{\bf D}} \newcommand{\bfE}{{\bf E}} \newcommand{\bfF}{{\bf F}} \newcommand{\bfG}{{\bf G}} \newcommand{\bfH}{{\bf H}} \newcommand{\bfI}{{\bf I}} \newcommand{\bfK}{{\bf K}} \newcommand{\bKs}{\bf K_{s}} \newcommand{\bKn}{\bf K_{n}} \newcommand{\bfL}{{\bf L}} \newcommand{\bfM}{{\bf M}} \newcommand{\bfN}{{\bf N}} \newcommand{\bfP}{{\bf P}} \newcommand{\bfR}{{\bf R}} \newcommand{\bfS}{{\bf S}} \newcommand{\bfT}{{\bf T}} \newcommand{\bfU}{{\bf U}} \newcommand{\bfV}{{\bf V}} \newcommand{\bfVhat}{\widehat{\bf V}} \newcommand{\bfrp}{{\bf r}^{\prime}} \newcommand{\bfrho}{\rho} \newcommand{\bfrpp}{{\bf r}^{\prime\prime}} \newcommand{\pmb}[1]{\setbox0=\hbox{#1}% \kern-.25em\copy0\kern-\wd0 \kern.05em\copy0\kern-\wd0 \kern-.25em\raise.0433em\box0 } \newcommand{\gamp}{\gamma^{\prime}} %\newcommand{\gamp}{\boldmath \mbox{$\gamma$}^{\prime}} \newcommand{\Rev}{\mbox{Re}} \newcommand{\Ai}{{\rm Ai}\,} \newcommand{\Bi}{{\rm Bi}\,} \newcommand{\pa}{\partial} \newcommand{\tri}{\Delta} \newcommand{\pr}{\prime} \newcommand{\del}{\Delta} \newcommand{\zw}{\rm\,} \newcommand{\vek}[1]{ {\bf #1} } \newcommand{\MK}[1]{ {\bf #1} } \begin{document} \Large \vspace{20 mm} %---------------------------------------------------------------------- \vspace{-10 mm} \title{Ocean Acoustic Modeling in MATLAB} \author{ Michael B. Porter \\ Science Applications International Corp. \\ 1299 Prospect St., Suite 305 \\ La Jolla, CA 92037 } \date{} \maketitle \begin{itemize} \item Rays \item Modes \item Wavenumber integration (FFP) \item Parabolic Equation \end{itemize} \newpage %---------------------------------------------------------------------- Governing equation (the wave equation): \[ \nabla ^2 p + { 1 \over c^2(r,z) } p_{tt} = {-\delta (r-r_s) \delta (z-z_s) \over r} \] Helmholtz Equation: \[ \nabla ^2 p + { \omega^2 \over c^2(\vek{x}) }\, p = {-\delta ( \vek{x} - \vek{x}_s )} \:, \] \newpage %---------------------------------------------------------------------- \heading{Ray Theory} Seek a solution of the Helmholtz equation in the following form: \[ p( \vek{x} ) = e^{i \omega\,\tau ( \vek{x} )}\, \sum _{j=0}^\infty {A_j( \vek{x} ) \over (i \omega)^j} \] where, $k = \omega /c_0$. \begin{itemize} \item called the {\em Ray Series}. \item generally divergent. \item provides an {\em asymptotic} approximation \end{itemize} \newpage %---------------------------------------------------------------------- \heading{Differentiate the ray series ...} \vspace{5 mm} \[ p_x = e^{ i \omega \tau } \left[ i \omega \tau _x \sum _{j=0}^\infty { A_j \over ( i \omega )^j} + \sum _{j=0}^\infty { A_{j,x} \over ( i \omega )^j} \right] \,, \] and \[ p_{xx} = e^{ i \omega \tau } \left\{ \left[ -\omega^2 (\tau _x)^2 + i \omega \tau _{xx} \right] \sum _{j=0}^\infty { A_j \over ( i \omega )^j} + 2 i \omega \tau_x \sum _{j=0}^\infty { A_{j,x} \over ( i \omega )^j} + \sum _{j=0}^\infty { A_{j,xx} \over ( i \omega )^j} \right\} . \] \newpage %---------------------------------------------------------------------- \heading{Substitute back in the Helmholtz equation ...} \vspace{5 mm} \be \nabla ^2 p = e^{ i \omega \tau } \left\{ \left[ -\omega^2 | \nabla \tau |^2 + i \omega \nabla ^2 \tau \right] \sum _{j=0}^\infty { A_j \over ( i \omega )^j} + 2 i \omega \nabla \tau \cdot \sum _{j=0}^\infty { \nabla A_j \over ( i \omega )^j} + \sum _{j=0}^\infty { \nabla ^2 A_j \over ( i \omega )^j} \right\} . \ee Equate terms of like order in $\omega$ \be \begin{array}{lrcl} O( \omega^2)\,\ \ : & \left| \nabla \tau \right|^2 & = & c^{-2}(\vek{x}) \hspace*{10mm} \mbox{Eikonal} \\ O( \omega )\:\ \ \ : & 2 \nabla \tau \cdot \nabla A_0 + (\nabla^2 \tau )A_0 & = & 0 \hspace*{10mm} \mbox{Transport}\\ O( \omega^{1-j} ): & 2 \nabla \tau \cdot \nabla A_j + (\nabla^2 \tau )A_j & = & -\nabla^2 A_{j-1}\,, \qquad j = 1, 2, \ldots \end{array} \ee \newpage %---------------------------------------------------------------------- \heading{Solving the Eikonal Equation (Method of Characteristics)} Define {\em rays} as curves perpendicular to the {\em wavefronts} of $\tau (\vek{x})$: \be { d \vek{x} \over ds } = c \,\nabla \tau \ee But, phase is still unknown. Lots of work \ldots \be { d \over ds } \left( { 1 \over c }\, { d \vek{x} \over ds } \right) = -{ 1 \over c^2 }\, \nabla c \ee Rays are now defined just in terms of $c( \vek{x} )$! \newpage %---------------------------------------------------------------------- \heading{Ray equations in cylindrical coordinates} \bea {dr \over ds} & = & c\, \xi ( s ) \,, \qquad {d \xi \over ds } = -{ 1 \over c^2 }\, { dc \over dr }\:, \\ {dz \over ds} & = & c\, \zeta ( s ) \,, \qquad {d \zeta \over ds } = -{ 1 \over c^2 }\, { dc \over dz } \eea $[ r(s), z(s) ]$ is the trajectory of the ray. Initial Conditions \bea r & = & r_s\,,\qquad \xi = { \cos \theta \over c( 0 ) } \:,\\ z & = & z_s\,,\qquad \zeta = { \sin \theta \over c( 0 ) } \eea \newpage %---------------------------------------------------------------------- \heading{Solving the eikonal equation} Eikonal equation: \be \nabla \tau \cdot \nabla \tau = { 1 \over c^2 } \:. \ee Therefore \be \nabla \tau \cdot { 1 \over c }\,{ d \vek{x} \over ds } = { 1 \over c^2 } \:, \ee or \be { d \tau \over ds } = { 1 \over c } \:. \ee (eikonal equation in ray coordinate $s$). Original nonlinear PDE is now a linear ODE Solution: \be \tau ( s ) = \tau ( 0 ) + \int _0 ^s { 1 \over c(s') }\, ds'\:. \ee Phase of the wave is delayed by its travel time. \newpage %---------------------------------------------------------------------- \heading{NORMAL MODES} Helmholtz equation (2-D): \be {1 \over r}\, {\partial \over \partial r} \left( r\, {\partial p \over \partial r} \right) + \, {\partial^2 p \over \partial z^2} + { \omega ^2 \over c^2(z) }\, p =- { \delta (r)\,\delta (z-z_s) \over 2 \pi r} \ee Solve by {\it separation of variables}. Seek $p(r,z) = \Phi(r) \Psi(z)$. \newpage %---------------------------------------------------------------------- \heading{Depth-separated equation} $\Psi( z )$ must satisfy: \bea { d^2 \Psi_m( z ) \over dz^2 } + \left[{ \omega ^2 \over c^2( z ) } - k_{ rm }^2 \right] \Psi_m(z) = 0 \\ \Psi(0) = 0 \\ \left.{d\Psi \over dz}\right|_{z=D} = 0 \eea This must be solved {\em numerically} for $k_m$ and $\Psi_m( z )$ (the eigenvalues and eigenfunctions of this Sturm-Liouville problem). \begin{itemize} \item There are an infinite number of modes! \item Fortunately, we can make do with a finite number \item They can be scaled arbitrarily \end{itemize} \newpage %---------------------------------------------------------------------- \heading{Depth-separated equation continued ...} Orthogonality property: \be \int_0^D \Psi_m(z)\, \Psi_n(z) \, dz = 0 \,, \qquad \mbox{for} \qquad m \neq n \ee Normalization: \be \int_0^D \Psi_m^2(z) \, dz = 1 \ee \newpage %---------------------------------------------------------------------- \heading{Range-separated equation} $\Phi( r )$ must satisfy: \be {1 \over r}\, {d \over dr} \left[ r\, {d \Phi_n(r) \over dr} \right] + k_{n}^2\, \Phi_n(r) = 0 \:. \ee The range functions are easy: \[ \Phi_n( r ) = H_0^{(1)}(k_{n} r). \] We can evaluate the Hankel functions exactly but if we're more than $\lambda$ away, we use: \be H_{0}^{(1)}( k r ) \rightarrow \left( \frac{2}{\pi k r} \right)^{1/2} e^{i(kr-\frac{\pi}{4})} \ee \newpage %---------------------------------------------------------------------- \heading{A summing up ...} \be p(r,z) = {i \over 4 } \sum _{m=1}^\infty \Psi_m(z_s)\, \Psi_m(z)\, H_0^{(1)}(k_{m} r) \ee or, using the asymptotic approximation to the Hankel function, \be p(r,z) \simeq {i \over \sqrt{8 \pi r} }\, e^{-i \pi / 4} \sum _{m=1}^\infty \Psi_m(z_s)\,\Psi_m(z)\, { e ^{i k_{m} r} \over \sqrt { k_{m} } } \ee Transmission Loss \be \mbox{TL}(r,z) \simeq -20 \log \left| \sqrt{{2 \pi \over r} } \sum _{m=1}^\infty \Psi_m(z_s)\, \Psi_m(z)\, { e ^{i k_{m} r} \over \sqrt { k_{m} } } \right| \:. \ee \newpage %---------------------------------------------------------------------- \heading{Example: Isovelocity profile (c(z) = 1500 m/s) } \begin{itemize} \item Solve for the modes \bea { d^2 \Psi_m( z ) \over dz^2 } + \left[{ \omega ^2 \over c^2 } - k_{ m }^2 \right] \Psi_m(z) = 0 \\ \Psi(0) = 0 \\ \left.{d\Psi \over dz}\right|_{z=D} = 0 \eea General solution: \be \Psi_m(z) = A \sin (k_z z) + B \cos (k_z z) \:, \ee where the vertical wavenumber $k_z$ is: \be k_z = \sqrt{\left( {\omega \over c }\right)^2 - k^2 } \:. \ee %---------------------------------------------------------------------- \item Top BC implies $B = 0$. \item Bottom BC implies \be Ak_z\, \cos (k_z D) = 0 \:, \ee requiring \be k_z D = \left( m-{1 \over 2} \right) \pi\,, \qquad m=1, 2, \ldots \,, \ee Thus, \be k_{m} = \sqrt{ \left({\omega \over c} \right)^2 -\left[ \left( m-{1 \over 2} \right) { \pi \over D} \right]^2}\,, \qquad m=1, 2, \ldots \ee \item Corresponding eigenfunctions are given by \be \Psi_m(z) = \sqrt{{2 \over D}}\, \sin (k_{zm} z) \ee \item Sum up the modes to get the pressure field \be p(r,z) = {i \over 2 D} \sum _{m=1}^\infty \sin(k_{zm} z_s) \sin(k_{zm} z)\, H_0^{(1)}(k_{m} r) \ee \end{itemize} \newpage %---------------------------------------------------------------------- \heading{Numerics} Recall the modal equation: \bea \Psi''(z) + \left[ { \omega ^2 \over c^2(z) } - k^2 \right] \Psi(z) & = & 0 \,,\\ \Psi(0) & = & 0 \,, \\ \Psi(D) & = & 0 \,. \eea Finite-difference approximation ($\psi_j = \Psi( z_j )$): \be \psi''_{j} = { \psi_{j-1} - 2\psi_{j} + \psi_{j+1} \over h^2 } + O(h^2) \:. \ee \bea { \psi_{j-1} \over h^2 } + \left\{ {-2 \over h^2 } + {\omega^2 \over c^2(z_j)} - k^2 \right\} \psi_{j} + { \psi_{j+1} \over h^2 } & = & 0 \,, \qquad j = 1, \ldots N-1\,, \nonumber\\ \psi_{0} & = & 0 \nonumber \\ \psi_{N+1} & = & 0 \eea \newpage %---------------------------------------------------------------------- \heading{Matrix form} \be \left[ \MK{A} -k^2 \MK{I} \right] \,\vek{\psi} = 0 \:. \ee This equation has $N$ eigenvalues $k^{(m)}$ and corresponding eigenvectors ${\bf \psi^{(m)} }$. $\vek{\psi^m}$ is the vector with components $\Psi_1^{m}, \Psi_2^{m}, \ldots \Psi_N^{m}$. Thus, $\psi_j^{(m)} \simeq \Psi_m( z_j )$. \be \MK{A} = \left[ \begin{array}{cccccccc} d_1 & e_2 & & & & & \\ e_2 & d_3 & e_3 & & & & \\ & e_3 & d_2 & e_4 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & e_{N-2} & d_{N-2} & e_{N-1} & \\ & & & & e_{N-1} & d_{N-1} & e_N\\ & & & & & e_N & d_N\\ \end{array} \right] \ee where, \be d_j = {-2 \over h^2 } + { \omega^2 \over c^2(z_j) }, \hspace*{10mm} e_j= {1 \over h^2 } \ee \newpage %---------------------------------------------------------------------- \heading{Mode Normalization} The integral term can be evaluated by the trapezoidal rule. That is, \be \int _0^D \Psi_m^2 (z) \, dz \simeq h \left( {1 \over 2} \phi_0 + \phi_1 + \phi_2 + \cdots + \phi_{N-1} + {1 \over 2}\phi_N \right) \:, \ee where \be \phi_j = \left( \psi_j^{(m)} \right)^2 \:. \ee \newpage %---------------------------------------------------------------------- \heading{Say it all in matrices:} \begin{enumerate} \item Solve the algebraic eigenvalue problem ${ \bf A \psi} = \lambda {\bf \psi}$ (${\bf A}$ defined above). Result: ${\bf \psi^{(m)}}, \lambda^{(m)}, m = 1, \ldots , N$. Note: $k_m \leftarrow \sqrt{ \lambda }$. \item Sort $k_m$ and keep the modes with the largest real part. \item Normalize the eigenfunctions \[ \psi^{(m)} \leftarrow { \psi^{(m)} \over \| \psi^{(m)} \| } \] \item Assemble the eigenvectors into a matrix \[ {\bf \psi} = \left[ \begin{array}{cccc} {\bf \psi^{(1)} } & {\bf \psi^{(2)} } & \ldots & {\bf \psi^{(M)} } \end{array} \right] \] ($\psi_{im}$ is the $i$th element of the $m$th mode.) \item Calculate mode excitation coefficients ${\bf C}$: \[ {\bf C} \leftarrow \left[ \begin{array}{c} \psi_{isd,1} \\ \psi_{isd,2} \\ \vdots \\ \psi_{isd,M} \\ \end{array} \right] \] \item Calculate ${\bf \tilde{\psi}}$, a mode matrix scaled by the mode excitation: \[ {\bf \tilde{\psi} } \leftarrow {\bf \psi} \left[ \begin{array}{ccc} 1/C_1 & & \\ & \ddots & \\ & & 1/C_M \\ \end{array} \right] \] \item Calculate the $\bf{\Phi}$ is a phase matrix: \[ {\bf \Phi} \leftarrow \left[ \begin{array}{ccc} 1/\sqrt{k_1} & & \\ & \ddots & \\ & & 1/\sqrt{k_M}\\ \end{array} \right] \left[ \begin{array}{ccccc} e^{i k_1 r_1} & e^{i k_1 r_2} & \cdots & \cdots & e^{i k_1 r_{NR}} \\ e^{i k_2 r_1} & e^{i k_2 r_2 }& & & e^{i k_2 r_{NR}} \\ \vdots & \vdots & & & \vdots \\ e^{i k_M r_1} & e^{i k_M r_2 }& \cdots & \cdots & e^{i k_M r_{NR}} \\ \end{array} \right] \left[ \begin{array}{ccc} 1/\sqrt{r_1} & & \\ & \ddots & \\ & & 1/\sqrt{r_{NR}}\\ \end{array} \right] \] \item Sum up the modes: \[ p = \tilde{ {\bf \psi} } {\bf \Phi} \] Recall, \be p(r,z) \simeq {i \over \sqrt{8 \pi r} }\, e^{-i \pi / 4} \sum _{m=1}^\infty \Psi_m(z_s)\,\Psi_m(z)\, { e ^{i k_{m} r} \over \sqrt { k_{m} } } \ee \end{enumerate} where, \newpage %---------------------------------------------------------------------- \heading{Adiabatic Modes} As the modes propagate in a range-dependent environment, they continually change shape. As they change shape, energy is continously re-distributed amongst the local modes. However, for sufficiently slowly-varying environments the energy mostly stays in the same mode. Then, the acoustic field can easily be computed using the {\em adiabatic mode formula}: \be p(r,z) \simeq { i \over \sqrt{8 \pi r} }\, e^{-i \pi / 4} \sum _{m=1}^\infty \Psi_m(z_s)\,\Psi_m(r,z)\, { e^{ i\int _0^r k_{m} (r') \, dr'} \over \sqrt{k_{m}(r)} } \:. \ee \newpage %---------------------------------------------------------------------- \heading{FAST FIELD PROGRAM} Apply a Fourier-Bessel transform to the Helmholtz equation: \[ g(k,z) = \int _0 ^\infty p(r,z) J_0( kr ) r dr \] \begin{eqnarray*} \Rightarrow {d^2 g \over dz^2} + \left( { \omega ^2 \over c^2(z) } - k^2 \right) g & = & \delta ( z - z_s ), \\ g(0) = 0, \hspace{10 mm} {dg \over dz} (D) & = & 0 \end{eqnarray*} The pressure field is reconstructed using the inverse transform: \[ p(r,z) = \int _0 ^\infty g(k,z) J_0(kr) k dk \] which can be efficiently evaluated using an FFT (Marsh, DiNapoli(1967) ): \[ p(r,z) \approx {e^{i\pi /4} \over \sqrt{2 \pi r} } \int _0 ^{K_{max}} g(k,z) \sqrt{k} e ^{ikr} dk . \] \newpage %---------------------------------------------------------------------- \heading{Matrix form} \be \left[ \MK{A} -k^2 \MK{I} \right] \,\vek{g} = \delta \:. \ee \be \MK{A} = \left[ \begin{array}{cccccccc} d_0 & e_1 & & & & & \\ e_1 & d_1 & e_2 & & & & \\ & e_2 & d_2 & e_3 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & e_{N-2} & d_{N-2} & e_{N-1} & \\ & & & & e_{N-1} & d_{N-1} & e_N\\ & & & & & e_N & d_N\\ \end{array} \right] \ee where, \bea d_j & = & {-2 \over h^2 } + { \omega^2 \over c^2(z_j) } \\ e_j & = & {1 \over h^2 } \eea \newpage %---------------------------------------------------------------------- \begin{itemize} \item This is a linear system of $N$ equations. \item $\delta$ is a function that is $1$ at the source depth and $0$ elsewhere. \item For each choice of the horizontal wavenumber $k$ we get a response $g(k,z)$. \end{itemize} \newpage %---------------------------------------------------------------------- \heading{Overview of the numerical procedure} \begin{enumerate} \item Set up the matrix ${\bf A}$. \item Define a sequence of wavenumbers \bea \{ k_j & = & k_{min} + j \Delta k + i \alpha, j = 1, NK \} \\ \Delta k & = & { k_{max} - k_{min} \over NK - 1 } \\ \alpha & = & \Delta k \eea \item Solve $ [ A - k_{j} I ] g = \Delta$ for each wavenumber $k_j$. \\ (Result is a matrix ${\bf g}$ with $g_{ij} \simeq g( z_i ; k_j )$ \item Multiply by $\sqrt{k}$: \[ {\bf g^T} \leftarrow {\bf g} \left[ \begin{array}{ccc} 1/\sqrt{k_1} & & \\ & \ddots & \\ & & 1/\sqrt{k_{NK}}\\ \end{array} \right] \] \item Do an FFT: $g( r, z ) = FFT( g( k, z) )$ \item Put in the cylindrical spreading and damping factor. \[ p = {\bf g} \left[ \begin{array}{ccc} e^{\alpha r_1} / \sqrt{r_1} & & \\ & \ddots & \\ & & e^{\alpha r_{NR}} / \sqrt{r_{NR}}\\ \end{array} \right] \] \end{enumerate} \newpage %---------------------------------------------------------------------- \heading{PARABOLIC EQUATION MODELING} (following Tappert and Hardin (1973)) Recall, Helmholtz equation: \be {\pa ^2 p\over{\pa r^2}} + {1\over r}{\pa p\over{\pa r}} + {\pa ^2 p\over{\pa z^2}} + k_0^2 n^2 p = 0, \ee Seek, \be p(r,z)=\psi(r,z) H_0^{(1)} (k_0 r), \ee where $H_0^{(1)} (k_0 r)$ is the Hankel function and satisfies: \be {\pa^2 H_0^{(1)}\over{\pa r^2}} + {1\over r}{\pa H_0^{(1)}\over{\pa r}} + k_0^2H_0^{(1)}=0, \ee We find, \be {\pa^2\psi \over{\pa r^2}} + \left({2\over{H_0^{(1)}}}{\pa H_0^{(1)} \over{\pa r}} + {1\over r}\right){\pa\psi \over\pa r} + {\pa^2\psi \over{\pa z^2}} + k_0^2(n^2-1)\psi=0. \ee \newpage %---------------------------------------------------------------------- \heading{PE derivation continued ...} Far-field approximation to Hankel function: \be {\pa^2\psi \over{\pa r^2}} + {{2ik_0}{\pa\psi \over\pa r}} + {\pa^2\psi \over{\pa z^2}} + k_0^2( n^2-1 )\psi=0. \ee (still elliptic and not marchable) Small-angle approximation (this is not well-motivated yet): \be {\pa^2\psi\over{\pa r^2}} \ll 2ik_0{\pa\psi\over{\pa r}}. \ee Gives the {\it standard parabolic equation} \be 2ik_0{\pa\psi\over{\pa r}} + {\pa^2\psi\over{\pa z^2}} + k_0^2(n^2-1)\psi=0, \ee \newpage %---------------------------------------------------------------------- \heading{Starting Fields} We can't start the PE with a delta function because it radiates uniformly in angle and the PE abuses the high-angle paths. Choices: \begin{itemize} \item Gaussian starter (just a smoothed delta function) \be \psi(0,z)=\sqrt{k_0}\,e^{-{k_0^2\over 2}(z-z_s)^2}, \ee \item Modal starter \item Greene's source \be \psi(0,z)=\sqrt{k_0}\left[1.4467-0.4201k_0^2(z-z_s)^2\right]e^{-{k_0^2\left( z-z_s\right)^2\over3.0512}}, \ee \end{itemize} \newpage %---------------------------------------------------------------------- \heading{Numerics} Let $\psi_j^i = \psi( r_i, z_j)$, i.e. ${\bf \psi}^i$ is a vector that approximates the field at range $r_i$. \[ {\pa^2\psi\over{\pa z^2}} + k_0^2(n^2-1)\psi \] is approximately ${\bf A \psi}$ where \[ \MK{A} = \left[ \begin{array}{cccccccc} d_1 & e_2 & & & & & \\ e_2 & d_3 & e_3 & & & & \\ & e_3 & d_3 & e_4 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & e_{N-2} & d_{N-2} & e_{N-1} & \\ & & & & e_{N-1} & d_{N-1} & e_N\\ & & & & & e_N & d_N\\ \end{array} \right] \] and, \bea d_j & = & {-2 \over h^2 } + k_0^2 ( n^2(z_j) - 1 ) \\ e_j & = & {1 \over h^2 } \eea \newpage %---------------------------------------------------------------------- \heading{Numerics (continued)} Recall, \be 2ik_0{\pa\psi\over{\pa r}} + {\pa^2\psi\over{\pa z^2}} + k_0^2(n^2-1)\psi=0, \ee Thus, \[ 2 i k_0 { \psi^{i+1} - \psi^i \over \Delta r } + A { \psi^{i+1} + \psi^i \over 2 } = 0 \] Rearrange, \[ \left[ { 2 i k_0 \over \Delta r } + { {\bf A} \over 2 } \right] {\bf \psi}^{i+1} = \left[ { 2 i k_0 \over \Delta r } - { {\bf A} \over 2 } \right] {\bf \psi}^{i } \] Or, \[ {\bf C \psi}^{i+1} = {\bf B \psi}^{i } \] \end{document}