prerequisite conditions
put velocity potential function as follow \[ \begin{equation} \phi=U x+\varphi-\frac{1}{2} U^2 tz \end{equation} \]
Preliminary: \[ \begin{array}{lll} \partial_t \phi=\partial_t \varphi-\frac{1}{2} U^2 ; & & \partial_t^2 \phi=\partial_t^2 \varphi \\ \partial_x \phi=U+\partial_x \varphi ; & & \partial_x^2 \phi=\partial_x^2 \varphi \\ \partial_y \phi=\partial_y \varphi ; & & \partial_y^2 \phi=\partial_y^2 \varphi \\ \partial_z \phi=\partial_z \varphi ; & & \partial_z^2 \phi=\partial_z^2 \varphi \end{array} \]
Dynamic boundary condition (on \(z= \eta\)) \[ \begin{equation} g\eta+\partial_t \phi+\frac{1}{2}\left[\left(\partial_x \phi\right)^2+\left(\partial_y \phi\right)^2+\left(\partial_z \phi\right)^2\right]-\frac{T}{\rho}\left(\partial_x^2 \eta+\partial_y^2 \eta\right)+p=0 \end{equation} \] where \[ \begin{equation} p=\bar p H(b-r)e^{i\omega_e t} \end{equation} \] \(H\) is Heaviside function, we posit \(\varepsilon\) as a small perturbation parameter, and define \(\phi,\eta,p\) as follows, \[ \begin{equation} \begin{aligned} & \phi=U x-\frac{1}{2} U^2 t+\varepsilon \varphi_1+\varepsilon^2 \varphi_2+\cdots \\ & \eta=\varepsilon \eta_1+\varepsilon^2 \eta_2+\cdots\\ &p=\varepsilon p_1+\varepsilon^2 p_2+\cdots \end{aligned} \end{equation} \] subsequently \(p_1=\bar p H(b-r)e^{i\omega_e t}\), and \(\varphi\) satisfies the Laplace equation \[ \begin{equation} \nabla^2 \phi=0 \end{equation} \] We consider only first-order terms, neglecting terms of higher order in '\(\varepsilon\)', thus, the water surface is described by \(z=0\). The kinematic boundary conditions are specified: \[ \begin{equation} \partial_z \phi=\partial_t \eta+\partial_x \phi \partial_x \eta+\partial_y \phi \partial_y \eta \quad \text{on } z=0 \end{equation} \]
Core Methoodology
Substituting equations (5) into equations (3), (4) and (6), and considering solely the first-order terms in \(\varepsilon\) yields
\(\varepsilon\): \[ \begin{equation} \begin{array}{l} \partial_x^2 \varphi_1+\partial_y^2 \varphi_1+\partial_z^2 \varphi_1=0\\ \varepsilon g \eta_1+\varepsilon \partial_t \varphi_1-\frac{1}{2} U^2+\frac{1}{2}\left(U+\varepsilon \partial_x \varphi_1\right)^2-\frac{T}{\rho}{\varepsilon}\left(\partial_x^2 \eta_1+\partial_y^2 \eta_1\right)+\varepsilon p_1=0\\ \varepsilon \partial_t \eta_1+\left(U+\varepsilon \partial_x \varphi_1\right)\left(\varepsilon \partial_x \eta_1\right)+\varepsilon \partial_y \varphi_1 \varepsilon \partial_y \eta_1=\varepsilon \partial_z \varphi_1 \end{array} \end{equation} \] Let \((\eta,\varphi)\) denote \((\eta_1,\varphi_1)\) and thereafter simplify equation \((8)\) leads to \[ \begin{equation} \begin{aligned} &\partial_x^2 \varphi+\partial_y^2 \varphi+\partial_z^2 \varphi=0 \\ &\left.\begin{array}{l} \partial_t \varphi+U \partial_x \varphi+g \eta-\frac{T}{\rho}\left(\partial_x^2 \eta+\partial_y^2 \eta\right)+p_1=0 \\ \partial_t \eta+U \partial_x \eta-\partial_z \varphi=0 \end{array}\right\} \text { on } z= 0 \end{aligned} \end{equation} \] we employ Fourier transforms in spatial time domain \[ \begin{equation} \begin{aligned} & \hat{\varphi}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \varphi e^{-i\left(k_x x+k_y y\right)} e^{-i \omega t} d x d y d t \\ & \hat{\eta}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \eta e^{-i\left(k_x x+k_y y\right)} e^{-i \omega t} d x d y d t \\ & \hat{p}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} p_1 e^{-i\left(k_x x+k_y y\right)} e^{-i \omega t} d x d y d t \end{aligned} \end{equation} \]
Calculations:
performing Fourier transform to \(\partial_x,\left(\partial_x^2+\partial_y^2\right), \partial_t\) \[ \begin{aligned} & F\left[\partial_x f\right]=i k_x \hat{f} \\ & F\left[\left(\partial_x^2+\partial_y^2\right) f\right]=-\left(k_x^2+k_y^2\right) \hat{f} \\ & F\left[\partial_t f\right]=i \omega f \end{aligned}\tag{c1} \]
Upon application of the Fourier transform, equation \((10)\) becomes \[ \begin{equation} \begin{aligned} &-\left(k_x^2+k_y^2\right) \hat{\varphi}+\partial_z^2 \hat{\varphi}=0 \\ &\left.\begin{array}{l} i \omega \hat{\varphi}+U i k_x \hat{\varphi}+g \hat{\eta}+\frac{T}{\rho}\left(k_x^2+k_y^2\right) \hat{\eta}+\hat{p}=0 \\ i \omega \hat{\eta}+U i k \hat{\eta}-\partial_z \hat{\varphi}=0 \end{array}\right\} \quad \text { on } z=0 \end{aligned} \end{equation} \]
Calculation:
we put \(k=\sqrt{k_x^2+k_y^2}\) then \[ -k^2 \hat{\varphi}+\partial_z^2 \hat{\varphi}=0 \tag{c2} \] According to the method of characteristic functions, the solution to this differential equation is readily obtained: \[ \hat{\varphi}=B(k, \omega) e^{-k z}\tag{c3} \] It follows that on \(z=0\) we obtain \[ \partial_z \hat{\varphi}=-k B; \quad \hat{\varphi}=B\tag{c4} \] Also, substituting it into the kinematic boundary condition yields \[ \hat{\eta}= \frac{\partial_z \hat{\varphi}}{i \omega+i U k_x} \quad (\text {on } z=0)\tag{c5} \] Substituting this into dynamic boundary condition on \(z=0\) results in \[ \left(i \omega+i U k_x\right) \hat{\varphi}+\left(g+\frac{T}{\rho} k^2\right) \hat{\eta}+\hat{p}=0\tag{c6} \] we put \(\omega_0=gk+\frac{T}{\rho}k^2\), the preceding expression can be simplified to \[ (i\omega +iUk_x)\hat\varphi+\frac{\omega_0}{k}\frac{\partial_z \hat\varphi}{i\omega+iUk_x}+\hat{p}=0\tag{c7} \] then \(\Rightarrow\) \[ \begin{align*} & \left(i \omega+i U k_x\right)^2 \hat{\varphi}+\frac{\omega_0^2}{k} \partial_z \hat{\varphi}+\left(i \omega+i U k_x\right) \hat{p}=0 \\ \Rightarrow & \\ & \frac{\omega_0^2}{k} \partial_z \hat{\varphi}-\left(\omega+U k_x\right)^2 \hat{\varphi}+\left(i \omega+i U k_x\right) \hat{p}=0 \end{align*}\tag{c8} \] substituting equation \((c3)\) into \((c8)\) yields \[ \begin{align*} & \frac{\omega_0^2}{k}(-k B)-\left(\omega+U k_x\right)^2 B+\left(i \omega+i U k_x\right) \hat{p}=0 \\ & {\left[\omega_0^2-\left(\omega+U k_x\right)^2\right] B+\left(i \omega+i U k_x\right) \hat{p}=0} \\ & B=\frac{\left(i \omega+i U k_x\right) \hat{p}}{\left(\omega+U k_x\right)^2-\omega_0^2} \end{align*} \] It follows that \[ \begin{align*} &\begin{aligned} & \hat{\varphi}=\frac{i\left(\omega+U k_x\right) \hat{p}}{\left(\omega+U k_x\right)^2-\omega_0^2} e^{-k z} \\ & \hat{\eta}=\frac{1}{i \omega+i U k_x}(-k) \frac{i\left(\omega+U k_x\right) \hat{p}}{\left(\omega+U k_x\right)^2-\omega_0^2} \end{aligned}\\ &=-\frac{k \hat{p}}{\left(\omega+U k_x\right)^2-\omega_0^2} \end{align*} \]
It's readily obtained that: \[ \begin{equation} \begin{aligned} & \hat{\varphi}=\frac{i\left(\omega+U k_x\right) \hat{p}}{\left(\omega+U k_x\right)^2-\omega_0^2} e^{-k z} \\ & \hat{\eta}=-\frac{k \hat{p}}{\left(\omega+U k_x\right)^2-\omega_0^2} \end{aligned} \end{equation} \] we now consider the form of \(\hat{p}\), since \(p\) has been previously assumed in equation \((3)\). We perform a spatiotemporal Fourier transform on \(p\) (spatial coordinate \(r\) and temporal coordinate \(t\)), resulting in: \[ \begin{equation} \int_{-\infty}^{\infty} e^{i \omega_e t} e^{-i \omega t} d t=2 \pi \delta(\omega-\omega_e) \end{equation} \] and \[ \begin{equation} \iint_{R^2} \bar{p} H(r-b) e^{-i\left(k_x x+k_y y\right)} d x d y=\frac{2 \pi b}{k} \bar{p} J_1(k b) \end{equation} \]
Heaviside function and Fourier transform
The initial hypotheses of pressure formulated as follows \[ \begin{equation} p=\bar{p} H(r-b) \end{equation} \]
Fourier transform for \(H(r-b)\)
We perform Fourier transform on function \(\bar{p}H(r-b)\) \[ \begin{equation} \iint_{R^2} \bar{p} H(r-b) e^{-i\left(k_x x+k_y y\right)} d x d y \end{equation} \] We transform the integral into polar coordinates, the coordinate transformation relations take the following form: \[ \begin{equation} \left\{\begin{array} { l } { k _ { x } = - k \operatorname { c o s } \beta } \\ { k _ { y } = k \operatorname { s i n } \beta } \end{array} \left\{\begin{array}{l} x=r \cos \theta \\ y=r \sin \theta \end{array}\right.\right. \end{equation} \]
Calculation \[ \begin{aligned} k_x x+k_y y & =-k r \cos \theta \cos \beta+k r \sin \theta \sin \beta \\ & =k r \cos (\beta+\theta) \end{aligned} \] and \[ d x d y=r d r d \theta \]
therefore equation (15) becomes \[ \begin{equation} \int_R \int_0^{2 \pi} \bar{p} H(r-b) e^{-i k r \cos (\beta+\theta)} r d r d \theta \end{equation} \] Over a complete integration period \(2\pi\), it follows immediately that \(\cos(\beta+\theta)=\cos \theta\), then equation (17) becomes \[ \begin{equation} \begin{aligned} & \int_R \bar{p} H(r-b) \int_0^{2 \pi} e^{-i k r \cos \theta} d \theta d r \\ \Rightarrow \Rightarrow & \int_0^{2 \pi} e^{-i k r \cos \theta} d \theta=2 \pi J_0(-k r)=2 \pi J_0(k r) \end{aligned} \end{equation} \]
other part
the water surface is (old) \[ \eta=i \rho \frac{1}{T} \int_{\frac{\pi}{2}-\theta}^{\frac{\pi}{2}} \frac{b J_1\left(k_g b\right)}{k_g-k_T} e^{-i r k_g \cos (\beta+\theta)} d \beta \]
\[ \begin{aligned} & \eta=\rho \frac{1}{T} \frac{b J_1\left(k_g b\right)}{k_g-k_T}\left\{\sqrt{\frac{2 \pi}{r\left|\psi^{\prime \prime}\left(\beta_0\right)\right|}} \sin \left[\psi\left(\beta_0\right) r+\frac{\pi}{4} \operatorname{sgn} \psi^{\prime \prime}\left(\beta_0\right)\right]\right\} \\ & k_T-k_g=\frac{1}{T} \sqrt{\rho^2 U^4 \cos ^4 \beta_0-4 g \rho T}\\ & \rho^2 U^4 \cos ^4 \beta_0-4 g \rho T>0 \\ & k_g=\frac{\rho U^2 \cos ^2 \beta_0-\sqrt{\rho^2 U^4 \cos ^4 \beta_0-4 g T \rho}}{2 T} \\ & \psi(\beta)=k_g \cos (\beta+\theta) \\ & \psi^{\prime}\left(\beta_0\right)=0 \\ & \psi^{\prime \prime}\left(\beta_0\right) \neq 0 \\ & \frac{\pi}{2}-\theta<\beta_0<\frac{\pi}{2} \\ & \frac{\pi}{2}-\theta<\beta<\frac{\pi}{2} \\ & 0<\theta<\frac{\pi}{2}\\ & r>0, r \text{ is large} \end{aligned} \] these are constant values \[ \left\{\begin{array}{l} T: 0.07 \mathrm{~N} / \mathrm{m} \\ \rho: 10^3 \mathrm{~kg} / \mathrm{m}^3 \\ b: 10 \mathrm{~m} \end{array}\right. \] now I want to plot \(\eta\) with respect to \(r\) and \(\theta\) can you give the c++ code? you can using some extra package or library
\[ \Phi_t+U \Phi_{x_1}+g \eta-\frac{T}{\rho} \eta_{x_1 x_1}=-f(\mathbf{x}) \]