linked=ValueQ[linked`$callingNotebook];Clear[linked`$callingNotebook]

NUMLAB2
Francis Auclair, Pierre-Antoine Dumont

Acoustic - Gravity Waves
Kernel: dispersion relations

Model

Only 2 kinds of AGWAVES (acoustic and surface) are conserved in the present non-stratified ocean. This means that \hat{\rho}_h=\rho_0 or  D-> \infinity or (N^2) incompressible=0.

Remove[RelP0,RelFS,RelRL,α12,ω,w,F,FF,FF1,λζ]

System of equations (simplest model , 4 unknowns and 4 equations)

Assumptions:
- no stratification: ρ=(Subscript[ρ, 0]),
- several possible surface BC: free-surface,RL, DT, P0...
- no penetration bottom BC,
- sediment layer,
- continuity of pressure and orthogonal velocity at the interface ocean / sediment (when sediment layer),
- 1st order vertical advection “only”,
- NBQ / NHINT CROCO kernels.

Inner relations

Eq1[x_,z_,t_]=(∂_t)u[x,z,t]+(∂_x)p[x,z,t]+λζ g(∂_x)ζ[x,t];
Eq2[x_,z_,t_]=(∂_t)w[x,z,t]+(∂_z)p[x,z,t]+ g ρ[x,z,t] ;
Eq3[x_,z_,t_]=(∂_t)ρ[x,z,t]+((∂_x)u[x,z,t]+(∂_z)w[x,z,t]);
Eq4[x_,z_,t_]=(∂_t)p[x,z,t]+ λζ g(∂_t)ζ[x,t]+(cs^2)((∂_x)u[x,z,t]+(∂_z)w[x,z,t])-g w[x,z,t];

Boundary relations (surface)

surfFS[x_,t_]=(∂_t)p[x,0,t]-g w[x,0,t];
surfDT[x_,t_]=(∂_t) ζ[x,t]+ (∂_x)(∫_(-Hi)^0)u[x,z,t]\[DifferentialD]z  ;
(* surfP0[x_,t_]: W'[0]=0 & ζ=0    *)
(* surfRL[x_,t_]: W[0] = 0         *)

Notes:
λζ=0: as in AGWaves 2021 in a non-stratified ocean. p is total pressure pressure (Subscript[p, t]).
λζ=1: the (Subscript[ρ, 0]) surface-induced term is treated separately both in p (Eq1 and Eq4) but ρ is still (Subscript[ρ, 0])+δρ. In this case (Subscript[p, t])= (Subscript[ρ, 0])gζ+ p where p is hydrostatic due to (Subscript[ρ, 0]) (no surface elevation anomaly) and to compressibility δp. Caution this is different from (Subscript[p, t])=(Subscript[p, 0]) +δp with (Subscript[p, 0])=(Subscript[ρ, 0])g(ζ+z) since (Subscript[p, 0])is here separated in two parts. When implemented with surfDT this relation requires to add a surface-induced component to w to obtain the total vertical velocity  => (Subscript[w, t])= w +(Subscript[w, h]) where (Subscript[w, h]) is \partial ζ / \partial t at the surface and 0 at the bottom (see bellow).

Cancellations / associations :
-  weight (Eq2) with vertical advection (1st order Eq4) and W'[z] in W propagation equation.
-  surface pressure force (Eq1) with tendency of surface anomaly (Eq4).

Decoupling NHINT/HCOMP/H3D:
- w(tot)=w(comp)+w(surface hydro)... 
- w=w(comp) prognosticated with (Eq2)
- this approximate formulation decouples the inner dispersion relation and the surface boundary dispersion relation (no ζ in the inner, no compressibility in the surface boundary).
- since dispersion relations are independent: (Ω, kx) can be different for surface and inner-field waves.

Substitution of Fourier mode into 3-equation models

(*Print["Caution!"]*)
Eq1[x,z,t]/.{u[x,z,t]->U[z]Exp[I(kx*x-Ω*t)],w[x,z,t]->W[z]Exp[I(kx*x-Ω*t)],p0[x,z,t]->P0[z]Exp[I(kx*x-Ω*t)],p[x,z,t]->P[z]Exp[I(kx*x-Ω*t)],ρ[x,z,t]->R[z]Exp[I(kx*x-Ω*t)],wh[x,z,t]->WH[z]Exp[I(kx*x-Ω*t)]};
%/.Derivative[a1_,a2_,a3_][u][x,z,t]:>D[U[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][w][x,z,t]:>D[W[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][wh][x,z,t]:>D[WH[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a3_][ζ][x,t]:>D[Z Exp[I(kx*x-Ω*t)],{x,a1},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p][x,z,t]:>D[P[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
E1[z_]=(%(-I (E^(-I (kx x-t Ω))))//Expand)(*-kx P[0]*)

(*Print["Caution!"]*)
Eq2[x,z,t]/.{u[x,z,t]->U[z]Exp[I(kx*x-Ω*t)],w[x,z,t]->W[z]Exp[I(kx*x-Ω*t)],p0[x,z,t]->P0[z]Exp[I(kx*x-Ω*t)],p[x,z,t]->P[z]Exp[I(kx*x-Ω*t)],ρ[x,z,t]->R[z]Exp[I(kx*x-Ω*t)],wh[x,z,t]->WH[z]Exp[I(kx*x-Ω*t)]};
%/.Derivative[a1_,a2_,a3_][u][x,z,t]:>D[U[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][w][x,z,t]:>D[W[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][w][x,0,t]:>D[W[0]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][wh][x,z,t]:>D[WH[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p0][x,z,t]:>D[P0[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p][x,z,t]:>D[P[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
(*%/.Derivative[a1_,a2_,a3_][p][x,0,t]:>D[P[Z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];*)
E2[z_]=(%(-I (E^(-I (kx x-t Ω))))//Expand) (*+I P'[0]*)

Eq3[x,z,t]//.{u[x,z,t]->U[z]Exp[I(kx*x-Ω*t)],w[x,z,t]->W[z]Exp[I(kx*x-Ω*t)],p0[x,z,t]->P0[z]Exp[I(kx*x-Ω*t)],p[x,z,t]->P[z]Exp[I(kx*x-Ω*t)],ρ[x,z,t]->R[z]Exp[I(kx*x-Ω*t)],wh[x,z,t]->WH[z]Exp[I(kx*x-Ω*t)]};
%/.Derivative[a1_,a2_,a3_][u][x,z,t]:>D[U[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][w][x,z,t]:>D[W[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][wh][x,z,t]:>D[WH[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p0][x,z,t]:>D[P0[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p][x,z,t]:>D[P[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][ρ][x,z,t]:>D[R[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
E3[z_]=%(-I (E^(-I (kx x-t Ω))))//Expand

Eq4[x,z,t]//.{u[x,z,t]->U[z]Exp[I(kx*x-Ω*t)],w[x,z,t]->W[z]Exp[I(kx*x-Ω*t)],p0[x,z,t]->P0[z]Exp[I(kx*x-Ω*t)],p[x,z,t]->P[z]Exp[I(kx*x-Ω*t)],ρ[x,z,t]->R[z]Exp[I(kx*x-Ω*t)],wh[x,z,t]->WH[z]Exp[I(kx*x-Ω*t)]};
%/.Derivative[a1_,a2_,a3_][u][x,z,t]:>D[U[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][w][x,z,t]:>D[W[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][wh][x,z,t]:>D[WH[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p0][x,z,t]:>D[P0[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p][x,z,t]:>D[P[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a3_][ζ][x,t]:>D[Z Exp[I(kx*x-Ω*t)],{x,a1},{t,a3}];
E4[z_]=%(-I (E^(-I (kx x-t Ω))))//Expand

surfFS[x,t];
surfFS[x,t]//.{u[x,z,t]->U[z]Exp[I(kx*x-Ω*t)],w[x,z_,t]->W[z]Exp[I(kx*x-Ω*t)],p0[x,z,t]->P0[z]Exp[I(kx*x-Ω*t)],p[x,z,t]->P[z]Exp[I(kx*x-Ω*t)],wh[x,z,t]->WH[z]Exp[I(kx*x-Ω*t)]};
%/.Derivative[a1_,a2_,a3_][u][x,z,t]:>D[U[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][w][x,0,t]:>D[W[0]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][wh][x,z,t]:>D[WH[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p0][x,z,t]:>D[P0[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p][x,0,t]:>D[P[0]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
EsurfFS=%(-I (E^(-I (kx x-t Ω))))//Expand

Tricky way to deal with integraL... I did not find better treatment so far.

surfDT[x,t]
surfDT[x,t]/. (∫_(-Hi)^0)Derivative[1,0,0][u][x,z,t]\[DifferentialD]z->I kx Exp[I(kx*x-Ω*t)] (∫_(-Hi)^0)U[z]\[DifferentialD]z;
%/.(∫_(-Hi)^0)Derivative[0,0,1][ρ][x,z,t]\[DifferentialD]z-> -I Ω Exp[I(kx*x-Ω*t)] (∫_(-Hi)^0)R[z]\[DifferentialD]z;
%/.{u[x,z,t]->U[z]Exp[I(kx*x-Ω*t)],w[x,z,t]->W[z]Exp[I(kx*x-Ω*t)],p0[x,z,t]->P0[z]Exp[I(kx*x-Ω*t)],p[x,z,t]->P[z]Exp[I(kx*x-Ω*t)],wh[x,z,t]->WH[z]Exp[I(kx*x-Ω*t)]};
%/.Derivative[a1_,a2_,a3_][u][x,z,t]:>D[U[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][w][x,z,t]:>D[W[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a2_,a3_][wh][x,z,t]:>D[WH[z]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
%/.Derivative[a1_,a3_][ζ][x,t]:>D[Z Exp[I(kx*x-Ω*t)],{x,a1},{t,a3}];
%/.Derivative[a1_,a2_,a3_][p][x,0,t]:>D[P[0]Exp[I(kx*x-Ω*t)],{x,a1},{z,a2},{t,a3}];
EsurfDT=%(-I (E^(-I (kx x-t Ω))))//Expand

Surface induced component of vertical velocity

WHp[z_]=(-I Ω) λζ Z(((Hi+z))/Hi)

U, DP as functions of W: 1st stage

U1[z_]=U[z]/.First@Solve[(E1[z]==0),U[z]];
US[z_]=U[z]/.First@Solve[(E3[z]==0),U[z]];
R1[z_] = R[z] /. First@Solve[E3[z] == 0 /. U[z] -> U1[z]/.WH'[z]-> WHp'[z], R[z]];

Collect[((E4[z]  )/. U[z] -> U1[z]/.WH[z]-> WHp[z]/.WH'[z]-> WHp'[z])//Simplify,P[z]]== 0

P1[z_] = Collect[Collect[P[z] /. First[Solve[E4[z] == 0 /. U[z] -> U1[z]/.WH[z]-> WHp[z]/.WH'[z]-> WHp'[z], P[z]]] ,{W[z], Derivative[1][W][z], Derivative[2][W][z]}]//ExpandAll//Simplify,{Z,λζ,g}]
P1P[z_]=P1'[z];

U1[z]
US[z]
R1[z]
P1[z]
P1P[z]

2nd-stage relations

To derive the relations above (and in particular the relation for P1[z]), one had to supposed the relation above was not satisfied. This means that this particular case will have to be examined below...

U2[z_]=Collect[U1[z]/.P[z_]->P1[z],{W[z],W'[z],W''[z]}]//Simplify
U2b[z_]=Collect[U1[z]/.P[z_]->0,{W[z],W'[z],W''[z]}]//Simplify

Remove[R2]
R2[z_]=Collect[R1[z]/.P[z]->P1[z],{W[z],W'[z],W''[z]}]//Simplify

P2[z_]=P1[z]

P2P[z_]=P2'[z]

Summary

Model equations

E1[z]
E2[z]
E3[z]

Boundary conditions:

EsurfFS
EsurfDT

Propagation equation for W

Recall functions...

E1[z]
E2[z]
E3[z]
E4[z]

P1[z]
P2[z]
P2P[z]

R1[z]/.P[z]-> P1[z]
R2[z]

Define new W1

W1[z_]=((W[z]/.First@Solve[E2[z]==0,W[z]])/.R[z]-> 0)+Wh[z]

W1[z] and Wh[z]

W1'[z]

(*Whp[z_]=(-I Ω) Z(((Hi+z))/Hi)*)
WHp[0]
WHp[-Hi]

Evolution equation for R[z]

E3[z]
{U1[z],W1[z]}

E3R[z_]=E3[z]/.{U[z]-> U1[z],W'[z]-> W1'[z]}//Simplify

E3R[z]/.{P[z]-> (cs^2)R[z],P''[z]-> (cs^2)R''[z]}//Expand
%/.Wh'[z]-> Whp'[z] ;
EqR2=Collect[(((Ω %)/(cs^2))//Simplify),{R[z],R''[z],Z}]

Due to the separation of variables, the horizontal derivative of the Laplacian of ρ is associated to R[z] !

Coefficient[EqR2,Z]//Together

Conclusion: no pb for long surface gravity waves, this term vanishes !

Evolution equation for P[z]

E4[z]
(cs^2)E3[z]//Expand

Coherence is necessary! Caution with the weight!

Propagation relation for W[z] ==> InnerRel

Model is not isochore and is non-hydrostatic, start from E2[z]

E2[z]

Inject

P2P[z]

and that’s it:

E2[z]
P2[z]
P2P[z]
R2[z]

E5[z_]=Collect[Simplify[(-(cs^2)(kx^2)+(Ω^2))I Ω E2[z]/.{WH[z_]-> WHp[z]}/.{P'[z_]->P2P[z],R[z]->R2[z],P[z]->P2[z]}/.WH'[z]-> WHp'[z]],{W[z],W'[z],W''[z]}];
E6a[z_]=Collect[Simplify[((E5[z])/(  (Ω^2)  (cs^2)  ))],{W[z],W'[z],W''[z]}];
E6[z_]=Collect[I E6a[z],{W[z],W'[z],W''[z],Z,Hi+z}];
E7[z_]=E6[z]//Simplify

(*((g^2) (kx^2)+(cs^2) (kx^2) (Ω^2)-(Ω^4))//Simplify//Expand
titi[z_]=Ω2/.First@Solve[(%/.Ω-> (Sqrt[Ω2]))==0,Ω2]//Simplify
titi[0]*)

Conclusion: this clearly confirms that the surface component of the pressure force does not influence the inner propagation in W !!! This component is treated with Wh[z].

------ to be re-derived ------
Same exact solution based on an equation for P and P’[0 and -H]=0

(*((W[z]/.First@DSolve[{E7[z]==0 && W[-H]==0},W[z],z])Sin[kz z]//Simplify)/Sin[kz z]//Simplify*)

Change of Variable: propagation relation for F[z] ==> InnerRel-F

- same kx and Ω (modes)
- same w and dP at the interface

duma=Coefficient[E7[z],W'[z]];
D1[z_]=Quiet[1/duma];
WP[z_]=F[z](E^((∫_z^0)(1/(2 D1[z1]))\[DifferentialD]z1));
Collect[Simplify@Collect[(E7[z]/.(E^((∫_z^0)(1/(2 D1[z1]))\[DifferentialD]z1))->1),{F[z],F'[z],F''[z]}],{F[z],F'[z],F''[z]}];
InnerRel[z_]=%;
dum1=InnerRel[z]/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]}//Simplify;
duma=Coefficient[%,F''[z]];
dum2=dum1/duma//Simplify
kz2p=Coefficient[%,F[z]]
InnerRelF[z_]=Collect[dum2,F[z]]/.Coefficient[dum2,F[z]]-> (kz^2)//Simplify

kz2p

Boundary conditions for W[z] and F[z] ==> BCFS-F,  BCDT-F

EsurfFS
EsurfDT

P2[0]

BCFS=EsurfFS((((cs^2) (kx^2)-(Ω^2)))/(I (cs^2) ))/.{P[0]-> P2[0]}//Simplify
BCFSF=Collect[(BCFS/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]})//Simplify,F[0]]

R2[z]
dum1r=R2[z]
dum2r=Coefficient[dum1r,W'[z]]
dum3r=Coefficient[dum1r,W[z]]
dum4r=Coefficient[dum1r,Z]

U2b[z]
dum1=U2b[z](*/.W'[z_]-> F'[z];*)/.cs-> cs1;
dum2=Coefficient[dum1,W'[z]]
dum3=Coefficient[dum1,W[z]]
dum4=Coefficient[dum1,Z]

dum0=EsurfDT//Simplify

BCDT=dum0((cs1^2) (kx^2)-(Ω^2))/I/.{(∫_(-Hi)^0)U[z]\[DifferentialD]z->(*dum3(∫_(-Hi)^0)W[z]\[DifferentialD]z+dum2(W[0]-W[-Hi])*)+dum4 Z (*λζ*) Hi(* ,(∫_(-Hi)^0)R[z]\[DifferentialD]z->(dum3r(∫_(-Hi)^0)W[z]\[DifferentialD]z+dum2r(W[0]-W[-Hi])+dum4r Z (*λζ*) Hi)*) }/.cs-> cs1//Simplify
BCDT/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]};
BCDTF=(%)//Simplify

Summary and update

Propagation relation (inner column)

InnerRel[z]
InnerRelF[z]

Boundary conditions

BCFS//Simplify
BCFSF

BCDT//Simplify
BCDTF

Inner solutions in ocean and sediment layers

δz: two ratios

(Subscript[ϵ, z]) ?

Work on ϵz small parameter (normalized form... or not)

(*=(kz2/kz1)((Hb-Hi)/Hi)*)
ϵzp=ϵh Sqrt[((kz2p/.cs-> cs2)/(kz2p/.cs-> cs1))]
%/.cs2-> ϵ12 cs1/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.δz2->δz1 ϵz/.kzc-> kzc1/.kzc1-> δzc1 / Hi//Simplify//Expand;
ϵzpn=FullSimplify[%,Assumptions{ϵ12∈Reals,ϵ12>0}]//PowerExpand

Particular cases:

ϵzpn/.{ϵ12-> 1}//Simplify

ϵzpn/.{ϵa1-> 0}//Simplify;
FullSimplify[%,Assumptions{ϵ12∈Reals,ϵ12>0}]//PowerExpand

ϵzc ratio ?

ϵzcp= Sqrt[(((Ω^2)/(cs1^2)-(kx^2))/((Ω^2)/(cs2^2)-(kx^2)))]//Simplify
ϵzcpn= ϵzcp/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb-> Hi /ϵhi//Simplify

Particular case:

ϵzcpn/.ϵ12-> 1

ϵzcpn/.{ϵa1-> 0}//Simplify

Inner relations in ocean (1) and sediment (2) layers (EI1 & EI2)

Pressure as a function of F

P2F[z_]=P2[z]/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]}

Inner relations for ocean (1) and sediment (2) layers

InnerRelF[z]

EI1=InnerRelF[z]/.{cs->cs1,kz->kz1,F[z]-> F1[z],F''[z]-> F1''[z],F'[0]->F1'[0],F''[0]->F1''[0]}
EI2=InnerRelF[z]/.{cs->cs2,kz->kz2,F[z]-> F2[z],F''[z]-> F2''[z],F''[0]-> 0,F'[0]-> 0,Z-> 0}

P1[z]
P2[z]

PP1[z_]=P2F[z]/.{cs->cs1,F[z]-> F1[z],F'[z]->F1'[z]}//Simplify
PP2[z_]=(P2F[z]/.{cs->cs2,F[z]-> F2[z],F'[z]->F2'[z]})/.g-> 0/.Z->0//Simplify

Surface dispersion relation (if necessary)

Surface dispersion: no influence of P[z] (no acoustic, shallow waters)

E1[z]
E1[z]/.P[z]-> PP1[z]0 (* get rid of pressure gradient here *)
solU=U[z]/.First@Solve[%==0,U[z]]
dum1U=Coefficient[solU,Z] 
dum2U=Coefficient[solU,F1'[z]] 0

solR=R2[z]/.{cs->cs1,W[z]-> F1[z],W'[z]->F1'[z]}
dum1R=Coefficient[solR,Z] 0
dum2R=Coefficient[solR,F1'[z]]  0

dum0=EsurfDT/kx//Simplify
%/.{(∫_(-Hi)^0)U[z]\[DifferentialD]z-> dum1U Z Hi+dum2U(F1[0]-F1[-Hi])};
%/.{(∫_(-Hi)^0)R[z]\[DifferentialD]z-> dum1R Z Hi+dum2R(F1[0]-F1[-Hi])};
dum2=%
dum3=kx Ω dum2/Z/.λζ-> 1//Simplify//Together
dum3ss=%;

Ω2/.First@Solve[(dum3/.Ω-> (Sqrt[Ω2]))==0,Ω2]

Surface dispersion: influence of P[z] (acoustic, NBQ-like shallow water): useless...

E4[z]

E1[z]
EsurfDT

PP1[z]

E1[z]/.P[z]-> PP1[z]
solU=U[z]/.First@Solve[%==0,U[z]]
dum1U=Coefficient[solU,Z] 
dum2U=Coefficient[solU,F1'[z]]

(*solR=R2[z]/.{cs->cs1,W[z]-> F1[z],W'[z]->F1'[z]}
dum1R=Coefficient[solR,Z] 
dum2R=Coefficient[solR,F1'[z]]  *)

(*R2[z]
solR=R2[z]/.{cs->cs1,W[z]-> F1[z],W'[z]->F1'[z]}
dum1R=Coefficient[solR,Z] 
dum2R=Coefficient[solR,F1'[z]]  *)

dum0=EsurfDT/kx//Simplify
%/.{(∫_(-Hi)^0)U[z]\[DifferentialD]z->  Z (∫_(-Hi)^0)dum1U\[DifferentialD]z+dum2U(F1[0]-F1[-Hi])}
%/.{(∫_(-Hi)^0)R[z]\[DifferentialD]z-> dum1R Z Hi+dum2R(F1[0]-F1[-Hi])}
dum2=%/.{F1[-Hi]-> 0,F1[0]->0}
dum3=((kx^2) Ω dum2/Z)/.λζ-> 1//Simplify//Together
dum3wth=%;

Ω2/.Solve[(dum3/.Ω-> (Sqrt[Ω2]))==0,Ω2]

Comparison

dum3ss
dum3wth/(Ω^2)//Simplify

Only long surface gravity waves can propagate in both versions...

Inner solutions in both layers (FF1 and FF1ζ, FF2)

Deal with particular solution (if any):

EI1/.{F1[z]->0,F1'[z]->0,F1''[z]->0}
solP1=-%/Coefficient[EI1,F1[z]]//Expand

Solution in the ocean layer (cas λζ=0)

FF1[z_]=B1 Cos[kz1 z]+A1 Sin[kz1 z]
(*FF1[z_]=Collect[FF1[z]+(solP1/.{F1'[0]-> FF1'[0],F1''[0]-> FF1''[0]}),B1]*)

{FF1[0],FF1'[0],FF1''[0]}
PP1[0]
(*PP1[0]/.{F1'[0]-> FF1'[0],F1[0]-> FF1[0]}//Simplify//Expand
Simplify[%/.kz1-> Sqrt[kz2p]/.cs-> cs1,Assumptions->-(kx^2)+((Ω^2)/(cs1^2))>0]*)

Solution in the ocean layer (cas λζ=1)

FF1ζ[z_]=B1 Cos[kz1 z]+A1 Sin[kz1 z];
FF1ζ[z_]=Collect[FF1ζ[z]+(solP1/.{F1'[0]-> FF1'[0],F1''[0]-> FF1''[0]}),B1]

{FF1ζ[0],FF1ζ'[0],FF1ζ''[0]}
PP1[0]
(*PP1[0]/.{F1'[0]-> FF1'[0],F1[0]-> FF1[0]}//Simplify//Expand
Simplify[%/.kz1-> Sqrt[kz2p]/.cs-> cs1,Assumptions->-(kx^2)+((Ω^2)/(cs1^2))>0]*)

Solution in the sediment layer:

FF2[z_]=B2 Cos[kz2 z]+A2 Sin[kz2 z]
solBCB2=B2/.First@Solve[FF2[-Hb]==0,B2]

FF2bbc[z_]=((FF2[z]/.B2-> solBCB2)//TrigReduce)/. Sec[Hb kz2]-> 1
(*FF2bbc[z_]=A2 Sin[kz2 (z+Hb)]*)

Interfacial condition

- The interfacial condition is based on the continuity of pressure and on the continuity of the normal component of velocity.
-  In the present case (flat bottom), this is just the vertical velocity.
 - The surface induced component of pressure is not matched at the interface.

Work on pressure: total pressure is supposed to be continuous. Caution: continuity of pressure at the interface is not equivalent to continuity of W’, even when the rho0 surface-induced term is removed. 
- The component -i g Ω F1[z] is non zero and this term is associated to the advection of pressure (see Eq 2 in AGWaves 2021, 1st term in (N^2) parenthesis). It is canceled out the W’ coefficient of the W propagation equation by the weight.
- The rho0 surface-induced term does not appear in the W’ coefficient of the W propagation equation. Note that the corresponding component of pressure (((-(cs^2) g (kx^2) Z λζ+g Z λζ (ω^2)-i g ω W[z]])/((cs^2) (kx^2)-(ω^2)))) is z-dependant through W[z] ! (see Eq1 and Eq4 for the origin of the first two terms which are included in P (=Pt) in AGaves 2021). Here the impact of surface pressure can be studied separately.

- The additional term in the first-order advection term in Eq 11 of AGWaves 2021is associated to the stratification (ie \hat{\rho}_{h}). With these notations, stratification is proportional to 1/D, which means that no stratification <-> D-> \infinity. This term leads to the non-vanishing coefficient for W’ in the W propagation equation as a consequence of the advection of \rho this time in the EOS (Eq 11 of AGWaves). It also leads to the first component of pressure in Eq 18.c in this study ( (cs^2)/D W term).

Continuity of pressure & w at the interface

Recall first pressure and vertical component of velocity.

P1[z]//Together
W1[z]
WHp'[z]

(*PP1[z]-Coefficient[PP1[z],F1'[z]] WHp'[z]//Simplify*)
DPP1[z_]=PP1[z]/.Z-> 0

Interfacial relations corresponding to pressure and vertical velocity continuity

DPP1[z]-PP2[z]

(((DPP1[z]-PP2[z])/.z-> -Hi)/(-I Ω)//Simplify)/.(cs2^2) (kx^2)-(Ω^2)->- (cs2^2)(kzc2^2)/.-(cs1^2) (kx^2) Ω+(Ω^3)-> Ω (cs1^2)(kzc1^2)/.(cs1^2) (kx^2)-(Ω^2)->- (cs1^2)(kzc1^2)/.kzc2->kzc1/ϵzc //Simplify
INTP=%/.Z-> 0
INTW=(F1[z]-F2[z]/ϵρ)/.z-> -Hi// Simplify

FF1[z]
FF2[z]

[λζ=0] Pressure & w continuity +  BC(-Hb)==> F2  ==> IntRel[A1,B1]

P & W continuity:

INTPAB=INTP/.{F1[z_]-> FF1[z],F2[z_]-> FF2bbc[z],F1'[z_]-> FF1'[z],F2'[z_]-> FF2bbc'[z]}//Simplify
INTWAB=INTW/.{F1[z_]-> FF1[z],F2[z_]-> FF2bbc[z],F1'[z_]-> FF1'[z],F2'[z_]-> FF2bbc'[z]}

Compute A2 from W continuity:

solA2=A2/.First@Solve[INTWAB==0,A2]

Compute interface relation:

INTPAB2=INTPAB/.A2-> solA2/.kzc2-> kzc1 /ϵzc

INTPAB2;
IntRel0=Collect[(%Tan[Hb kz2-Hi kz2]/Cos[Hi kz1](kzc1^2)(cs1^2))//Simplify,{Tan[Hi kz1],Tan[Hb kz1-Hi kz1]}]

Final interfacial relation for F (normalized w) once bottom BC in sediment and interfacial relations have been substituted : IntRel depends on A1 and B1. The choice made here is to work with (A1,B1) couple of constant associated to the ocean layer.

IntRel0/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb-> Hi /ϵhi//Simplify;
IntRel=Simplify[%((Hi ϵh ((ϵ12)^2) )/(((ϵ12)^2)(cs1^2) )),Assumptions->{cs1>0,ϵa1>0,Hi>0}]

Exploration...

dum10=B1 δz1 ϵz ((ϵzc)^2)+(-A1 δz1 ϵh+B1 ((ϵa1)^2) ϵh) Tan[δz1 ϵz]-Tan[δz1] (A1 δz1 ϵz ((ϵzc)^2)+(B1 δz1+A1 ((ϵa1)^2)) ϵh Tan[δz1 ϵz]);

Coefficient[dum10,A1]//Expand
Coefficient[IntRel,A1]//Expand

Coefficient[dum10,B1]//Expand
Coefficient[IntRel,B1]//Expand

[λζ=1] Pressure & w continuity +  BC(-Hb)==> F2  ==> IntRel[A1,B1]

P & W continuity:

INTPABζ=INTP/.{F1[z_]-> FF1ζ[z],F2[z_]-> FF2bbc[z],F1'[z_]-> FF1ζ'[z],F2'[z_]-> FF2bbc'[z]}//Simplify
INTWABζ=INTW/.{F1[z_]-> FF1ζ[z],F2[z_]-> FF2bbc[z],F1'[z_]-> FF1ζ'[z],F2'[z_]-> FF2bbc'[z]}

Compute A2 from W continuity:

solA2=A2/.First@Solve[INTWABζ==0,A2]

Compute interface relation:

INTPAB2ζ=INTPABζ/.A2-> solA2/.kzc2-> kzc1 /ϵzc

INTPAB2ζ;
IntRel0ζ=Collect[(%Tan[Hb kz2-Hi kz2]/Cos[Hi kz1](kzc1^2)(cs1^2))//Simplify,{Tan[Hi kz1],Tan[Hb kz1-Hi kz1]}]

Final interfacial relation for F (normalized w) once bottom BC in sediment and interfacial relations have been substituted : IntRel depends on A1 and B1. The choice made here is to work with (A1,B1) couple of constant associated to the ocean layer.

IntRel0ζ/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb-> Hi /ϵhi//Simplify;
IntRelζ=Simplify[%((Hi ϵh ((ϵ12)^2) )/(((ϵ12)^2)(cs1^2) )),Assumptions->{cs1>0,ϵa1>0,Hi>0}]

Exploration...

dum10=B1 δz1 ϵz ((ϵzc)^2)+(-A1 δz1 ϵh+B1 ((ϵa1)^2) ϵh) Tan[δz1 ϵz]-Tan[δz1] (A1 δz1 ϵz ((ϵzc)^2)+(B1 δz1+A1 ((ϵa1)^2)) ϵh Tan[δz1 ϵz]);

Coefficient[dum10,A1]//Expand
Coefficient[IntRelζ,A1]//Expand

Coefficient[dum10,B1]//Expand
Coefficient[IntRelζ,B1]//Expand

Coefficient[IntRelζ,Z]//Expand

Surface boundary dispersion relations

Work directly on A1 and B1 coefficients to discuss surface boundary condition

[λζ=0] Rigid Lid (B1=0)                                                       ==> RelRLn[ω, δx, δz1]

Recall functions:

FF1[z]

Simplify interfacial relation:

dum0=(IntRel/.B1-> 0)/A1//Simplify

Simplified expression:

RelRLn=dum0//Expand//Simplify

Note that Ω does not appear explicitly.

Full expression:

RelRLn/.ϵzc-> ϵzcpn//Simplify
%(((-(kx^2)+((Ω^2)/(cs2^2))))^2)//Expand//Simplify;
%  (cs2^4)/.cs2-> ϵ12 cs1//Simplify
solω2RL=ω2/.First@Solve[(%/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify

Same acoustic velocity

solω2RLcs=solω2RL/.{ϵ12->1}

No sediment:

solω2RL/.{ϵh-> 0}

These are Lamb Waves

Jansen as a limit for infinite sediment layer depth:

dum0=RelRLn/. ϵz-> ϵh  (kz2/kz1);
dum1=dum0//Simplify;
dum2=dum1((kz2 δz2i)/(H2 kz1));
(*dum3=Normal@Series[dum2,{δz2i,0,2}]*)
dum4=((dum2/( δz2i))/.ϵzc->(( kz1)/( kz2))/.ϵh-> (H2/Hi)/.δz2->  kz2 H2/.δz1-> kz1 Hi/.ϵa1-> 0)//Simplify
dum5=dum4/.Tan[H2 kz2]-> I;
dum6=dum5

dum5/.{kz1-> δz1/Hi,kz2-> δz2/H2/I }/.H2-> ϵh Hi;
RelRLninf=% Hi ϵh//Simplify

[λζ=0] P=0 (A1=0) ==> B1sol                                         ==> RelP0n[ω, δx, δz1]

Caution... We only mean here that the acoustic component of pressure is vanishing. If we carefully look below the formulation of P1[0], the first term is associated to the free-surface, the second to the vertical “advection” of the pressure field and the third is the acoustic component.
In the configuration λζ=0 (the free-surface pressure component is included in p), we cancel the W’[0] term anyway (free surface anomaly is supposed to remain small)...

Recall functions

FF1[z]
FF1'[z]
P1[0]
P2F[0]/.{F[z_]-> FF1[z],F'[z_]-> FF1'[z]}
P2F[0]/.{F[z_]-> FF1[z],F'[z_]-> FF1'[z]}/.λζ-> 0/.A1-> 0

Simplify interfacial relation:

(IntRel/.A1-> 0)/B1//Simplify

Simplified expression:

RelP0n=% //Simplify

Note that Ω does not appear explicitly.

Full expression:

ϵzcpn

RelP0n/.ϵzc-> ϵzcpn
%(((-(kx^2)+((Ω^2)/(cs2^2))))^2)//Expand//Simplify;
%  (cs2^4)/.cs2-> ϵ12 cs1//Simplify
solω2P0=ω2/.First@Solve[(%/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify

Same acoustic velocity

solω2P0cs=solω2P0/.{ϵ12->1}

No sediment:

solω2P0/.{ϵh-> 0}

These are supposed to be Lamb Waves...

Jensen as a limit for infinite sediment layer depth:

dum0=RelP0n/. ϵz-> ϵh  (kz2/kz1);
dum1=dum0//Simplify;
dum2=dum1
(*dum3=Normal@Series[dum2,{δz2i,0,2}]*)
dum4=((dum2/( δz2i))((kz2 δz2i)/(H2 kz1))/.ϵzc->(( kz1)/( kz2))/.ϵh-> (H2/Hi)/.δz2->  kz2 H2/.δz1-> kz1 Hi/.ϵa1-> 0)//Simplify
dum5=dum4/.Tan[H2 kz2]-> I;
dum6=dum5

dum5/.{kz1-> δz1/Hi,kz2-> δz2/H2/I }/.H2-> ϵh Hi;
RelP0ninf=% Hi ϵh//Simplify

[λζ=0] Free surface (no sediment) (A1...B1...=0) ==> B1sol ==> RelFS1n[ω, δx, δz1]

FF1NS[z_]=A1 Sin[δz1(1+(z/Hi))]
(*B1sol=B1/.First@Solve[FF1[-Hi]==0,B1]/.Z-> 0//Expand//Simplify;
FF1NS[z_]=FF1[z]/.B1-> B1sol/.Z-> 0*)

BCFS;
BCFSF/.{F[z_]-> FF1NS[z],F'[z_]-> FF1NS'[z]}/.cs->cs1/.λζ-> 0

dum1=BCFSF/.{F[z_]-> FF1NS[z],F'[z_]-> FF1NS'[z]}//Simplify;
dum2=(% (Hi^3) /(A1 Cos[δz1](cs1^2) ((ϵa1)^2))) /.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.δz2->δz1 ϵz/.λζ-> 0//Simplify

RelFS1n=Collect[dum2,{Tan[δz1 ϵz],Tan[δz1]}]//Simplify

solω2FS1=ω2/.First@Solve[(RelFS1n/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify

kzcp2=((Ω^2)/(cs1^2))-(kx^2)/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}

Simplifications (1)

(*Limit[solω2FS1,δz1-> 0]*)

Simplifications (2)

[λζ=0] Free surface                                 (A1...B1...=0) ==> B1sol ==> RelFSn[ω, δx, δz1]

BCFSF

BCFSF/.λζ-> 0
%/.{F[z_]-> FF1[z],F'[z_]-> FF1'[z]}/.cs->cs1/.λζ-> 0
%/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb->Hi/ ϵhi
B1surf=B1/.First@Solve[%==0,B1]

IntRel

dum1=(IntRel/.B1-> B1surf);
dum2=(% ((δx)^2) ((ϵhi)^2)/A1) //Expand//Simplify

RelFSn=Collect[dum2//Simplify,{Tan[δz1 ϵz],Tan[δz1]}]//Simplify

Approximation... as if ϵzc did not depend on Ω (cs1=cs2)

solω2FSapp=ω2/.First@Solve[(RelFSn/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify

which means when same cs :

solω2FSapp/.{ϵzc->1,ϵ12->1}

under this approximation and without sediment

solω2FSapp/.{ϵhi->1,ϵh-> 0}

Jansen as a limit for infinite sediment layer depth:

dum0=RelFSn/. ϵz-> ϵh  (kz2/kz1);
dum1=dum0//Simplify;
dum2=dum1
(*dum3=Normal@Series[dum2,{δz2i,0,2}]*)
dum4=(((-dum2)/( δz2i))((kz2 δz2i)/(((ϵhi)^2)H2 kz1))/.ϵzc->(( kz1)/( kz2))/.ϵh-> (H2/Hi)/.δz2->  kz2 H2/.δz1-> kz1 Hi/.ϵa1-> 0)//Simplify
dum5=dum4/.Tan[H2 kz2]-> I;
dum6=dum5

dum5/.{kz1-> δz1/Hi,kz2-> δz2/H2/I }/.H2-> ϵh Hi;
RelFSninf=% Hi ϵh//Simplify

RelRLn
RelP0n
RelFSn

Recover P0 and RL when H2 is infinit...

RelRLninf
RelP0ninf
RelFSninf

RelRLninf
(-RelFSninf/.ω-> 0)/((δx)^2)//Simplify

RelP0ninf
-RelFSninf/(ω^2)//Simplify;
(%/.ω-> 1/ωi)//Simplify;
(%/.ωi-> 0)/δz1//Simplify

Full expression 1

RelFSn/.ϵzc-> ϵzcpn;
% (((((δx)^2)-((ϵa1)^2) ((ϵ12)^2) (ω^2)))/((cs1^4) ((ϵhi)^2) ((ϵ12)^4)))/.cs2-> ϵ12 cs1//Simplify;
RelFSnω=Collect[% (cs1^(4 ))((ϵ12)^4),ω]//Simplify

Collect[Normal@Series[RelFSnω,{ϵa1,0,0}]//Expand,ω]

RelFSnω;
%/. ϵhi ->( 1/(1+ϵh));
% (((1+ϵh))^2)(-1+Tan[δz1] Tan[ϵz])/.Tan[δz1 ϵz]-> ((Tan[δz1 ]+Tan[ϵz])/(1-Tan[δz1 ]Tan[ ϵz]))//Expand//Simplify
Collect[%,{Tan[δz1] ,Tan[ϵz]}]

Mathematica...

Full expression 2

dum0=RelFSn(((((δx)^2)-((ϵa1)^2) ((ϵ12)^2) (ω^2)))/((ϵhi)^2))/.ϵzc-> ϵzcpn//Simplify

((δz1)^2) ϵz (ω^2) (-((δx)^2)+((ϵa1)^2) (ω^2))+((δx)^2) δz1 ϵz (((δx)^2)-((ϵa1)^2) (ω^2)) Tan[δz1]+ϵh (((δx)^2)-((ϵa1)^2) ((ϵ12)^2) (ω^2)) (δz1 (((δx)^2)-((ϵa1)^2) (ω^2))+(((δx)^2) ((ϵa1)^2)+((δz1)^2) (ω^2)) Tan[δz1]) Tan[δz1 ϵz]

dum1=Normal@Series[dum0,{ϵ12,1,0}]//Simplify
ω2/.Solve[(dum1/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify

[λζ=1] Divergence of transport (no sediment)   ==> B1sol ==> RelDT1n[ω, δx, δz1]

Z here is not associated to the continuity of pressure at the interface but rather to the decomposition of pressure in the ocean layer and so to the formulation of surface boundary condition itself.

Particular case... No sediment

First cancel B1 (the surface-induced component of the pressure is treated independently meaning that the present w component (w=(Subscript[w, t])-(Subscript[w, h])) is associated to a RL surface boundary condition.

FF1ζ[z]/.B1-> 0
FF1ζ[-Hi]/.B1-> 0
(*B1sol=B1/.First@Solve[(FF1ζ[-Hi]/.A1->0)==0,B1]/.λζ->1//Expand*)
FF1NS[z_]=FF1ζ[z]/.B1-> 0//Simplify

BBC ==> A1=0, W is just Wh

BCDTF
BCDTF/.{F[z_]-> FF1NS[z],F'[z_]-> FF1NS'[z]}/.λζ->1//Simplify
%/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.δz2->δz1 ϵz/.λζ->1;
%  δz1 (Hi^3)/((cs1^2) ((ϵa1)^2))//Simplify
RelDT1n=Simplify[% (( ((ϵa1)^2)ω)/(-I cs1 Z δz1 ϵa1   )) ,Assumptions->{cs1>0,ϵa1>0,Hi>0}]//Simplify

(*solω2DT1=ω2/.First@Solve[(RelDT1n/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify*)

Simplification 0

(*(Numerator[solω2DT1]/Cos[δz1]//TrigReduce)/(Denominator[solω2DT1]/Cos[δz1]//TrigReduce)
Normal@Series[%,{ϵa1,0,2}]*)

Simplification 1 ( Hi -> 0 )

(*Limit[solω2DT1,{δz1-> 0,ϵa1->0}]*)

Simplification 2

(*Limit[solω2DT1,ϵa1-> 0]*)

Simplification 3

(*Limit[solω2DT1,δz1-> 0]
%/.ϵa1-> 0/.ϵhi-> 1*)

To be compared to:

(*RelFSn
RelRLn
RelP0n*)

[λζ=1] Divergence of transport                                   ==> B1sol ==> RelDTn[ω, δx, δz1]

First cancel B1 (the surface-induced component of the pressure is treated independently meaning that the present w component (w=(Subscript[w, t])-(Subscript[w, h])) is associated to a RL surface boundary condition.

FF1ζ[z]/.B1-> 0

BCDTF;
((BCDTF/.{F[z_]-> FF1ζ[z],F'[z_]-> FF1ζ'[z]})/(cs1^2))//Simplify
%/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb->  Hi/ϵhi//Simplify;
dum0=%(( (Hi^2)  ω)/( I  Z  (Sqrt[(((cs1^2) ((ϵa1)^2))/(Hi^2))])))
B1sol=0;

Work directly on A1 (and B1=0) coefficients to discuss surface boundary condition

(*B1sol=(B1/.First@Solve[dum0==0,B1])//Simplify*)

IntRel/.B1-> B1sol
dum1= %//Simplify;
RelDTn=%/A1//Simplify

Similar to rigid lid

(*(ω2/.First@Solve[(RelDTn/.ω-> (Sqrt[ω2]))==0,ω2])//Simplify;
solω2DT=Simplify[%//Expand,Assumptions->{ϵz>0,δx>0,δz1>0,δz1-((ϵa1)^2) Tan[δz1]<0,ϵhi>0}] (* TBJ *)*)

(*solω2DT/.{ϵh->0,ϵhi-> 1}//Expand
Simplify[%//Expand,Assumptions->{ϵz>0,δx>0,δz1>0,ϵhi>0,δz1-((ϵa1)^2) Tan[δz1]>0}] *)

Simplification 1:
Without sediment (no dev...)

(*dum10=Limit[solω2DT,{ϵz-> 0,ϵh->0,ϵhi-> 1} ]
dum11=(Numerator[dum10] Cos[δz1]//Expand//Simplify)/(Denominator[dum10] Cos[δz1]//Expand//Simplify)*)

To be compared with

(*solω2DT1*)

Long waves (Hi -> 0 )

(*Limit[dum11,{δz1-> 0,ϵa1->0}]*)

To be compared to:

RelFSn
RelRLn
RelP0n

Recover RL and P0 as asymptotic limits of Free-Surface model

Collect[RelRLn//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[(RelFSn/((ϵhi)^2)/((δx)^2)/.ω->0)//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[(RelFSnϵ/((ϵhi)^2)/((δx)^2)/.ϵ->0)//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]

Collect[RelP0n,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[Limit[(RelFSn/δz1/(ω^2)/((ϵhi)^2))//Expand,ω-> Infinity]//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[Limit[(RelFSnϵ/δz1/(ϵ^2)/((ϵhi)^2))//Expand,ϵ-> Infinity]//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]

DT (no sediment): limit w(0)

FF1NS[0]

To be verified:

EsurfFS

or

BCFSF
BCFSF/.{F[z_]-> FF1NS[z],F'[z_]-> FF1NS'[z]}/.cs->cs1;
(% (Hi^3)/A1/(cs1^2)/((ϵa1)^2))/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi};
%/ω/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.δz2->δz1 ϵz/.kzc-> kzc1/.kzc1-> δzc1 / Hi//Simplify

Compute the BC dispersion:

P2F[0]
(-ω P[0]+I g W[0])/.P[z_]-> P2[0]
(-ω P[0]+I g W[0])/.P[z_]-> P2[0]/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]};
%/.{F[z_]->FF1NS[z],F'[z_]->FF1NS'[z]}
dum20=(%/.kzc->Sqrt[ (ω^2)/(cs^2)-(kx^2)])//Together//Simplify

dum20/.cs->cs1//Simplify
%/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.δz2->δz1 ϵz/.kzc-> kzc1/.kzc1-> δzc1 / Hi;
% ((Hi (((δx)^2)-((ϵa1)^2) (ω^2)))/(I A1 (cs1^2) ((ϵa1)^2)))
dum21=Simplify[%,Assumptions->{cs1>0,ϵa1>0,Hi>0}]

Recover the BQ dispersion relation:

dum0=dum21//Simplify;
dum1=dum0/.g-> ((ϵa1)^2)(cs1^2)/Hi//Simplify;
dum2=-dum1 //Simplify;
dum3=ω2/.First@Solve[(dum2/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify;
dum4=(Numerator[dum3]/Cos[δz1]//Simplify)/(Denominator[dum3]/Cos[δz1]//Simplify)

EsurfFS

P2[0]

EsurfFS((((cs^4) (kx^2)-(cs^2)(Ω^2)))/(I (cs^2) ))/.{P[0]-> P2[0]}//Simplify
Collect[(%/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]})//Simplify,F[0]]/.kzc->Sqrt[ (Ω^2)/(cs^2)-(kx^2)] //Simplify

Remember indeed that:

solω2FS1

However from the divergence of transport:

(*solω2DT1*)

So compute the error made on EsurfFS: the relation to be satisfied is:

(*dum10=(solω2DT1-dum4)//Together//Simplify*)

The numerator must vanish:

(*dum11=Numerator[dum10]*)

Long wave limit:

(*Normal@Series[dum11,{δz1,0,2}]*)

Acoustic limit :

(*Normal@Series[dum11,{ϵa1,0,2}]*)

EsurfDT

Eq3[x,z,t]
EsurfFS

-Ω Integrate[R[z],{z,-Hi,0}]+kx Integrate[ U[z],{z,-Hi,0}]-I( W[-Hi]-W[0])//Expand

(*solω2DT1*)

[λζ=0] Summary: recover RelRLn & RelP0n=0  from Free Surface RelFSn

Recover RL and P0 as asymptotic limits of Free-Surface model

Collect[RelRLn//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[(RelFSn/.ω->0)/(-g (kx^2) ((ϵ12)^2))//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[(RelFSnϵ/.ϵ->0)/(-g (kx^2) ((ϵ12)^2))//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]

Collect[RelP0n,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[Limit[(RelFSn/(ω^2)/kz1/((ϵ12)^2))//Expand,ω-> Infinity]//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]
Collect[Limit[(RelFSnϵ/(ϵ^2)/kz1/((ϵ12)^2))//Expand,ϵ-> Infinity]//Simplify,{Tan[δz1 ϵz],Tan[δz1 ]}]

DT (no sediment): limit w(0)

FF1NS[0]

To be verified:

EsurfFS

or

BCFSF
BCFSF/.{F[z_]-> FF1NS[z],F'[z_]-> FF1NS'[z]}/.cs->cs1/.λζ-> 0;
(% (Hi^3)/A1/(cs1^2)/((ϵa1)^2))/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi};
%/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.δz2->δz1 ϵz/.kzc-> kzc1/.kzc1-> δzc1 / Hi//Simplify

Compute the BC dispersion:

P2F[0]/.λζ-> 0
(-Ω P[0]+I g W[0])/.P[z_]-> P2[0]/.λζ-> 0
(-Ω P[0]+I g W[0])/.P[z_]-> P2[0]/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]}/.λζ-> 0;
%/.{F[z_]->FF1NS[z],F'[z_]->FF1NS'[z]}/.λζ-> 0
dum20=(%/.kzc->Sqrt[ (Ω^2)/(cs^2)-(kx^2)])//Together//Simplify

dum20/.cs->cs1//Simplify
%/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.δz2->δz1 ϵz/.kzc-> kzc1/.kzc1-> δzc1 / Hi;
% ((Hi (((δx)^2)-((ϵa1)^2) (ω^2)))/(I A1 (cs1^2) ((ϵa1)^2)))
dum21=Simplify[%,Assumptions->{cs1>0,ϵa1>0,Hi>0}]

Recover the BQ dispersion relation:

dum0=dum21//Simplify;
dum1=dum0/.g-> ((ϵa1)^2)(cs1^2)/Hi//Simplify;
dum2=-dum1 //Simplify;
dum3=ω2/.First@Solve[(dum2/.ω-> (Sqrt[ω2]))==0,ω2]//Simplify;
dum4=(Numerator[dum3]/Cos[δz1]//Simplify)/(Denominator[dum3]/Cos[δz1]//Simplify)

EsurfFS

P2[0]/.λζ-> 0

EsurfFS((((cs^2) (kx^2)-(Ω^2)))/(I (cs^2) ))/.{P[0]-> P2[0]}/.λζ-> 0//Simplify
Collect[(%/.{W''[z_]-> WP''[z],W'[z_]-> WP'[z],W[z_]-> WP[z]})//Simplify,F[0]]/.kzc->Sqrt[ (ω^2)/(cs^2)-(kx^2)] //Simplify

Remember indeed that:

solω2FS1

However from the divergence of transport:

(*solω2DT1*)

So compute the error made on EsurfFS: the relation to be satisfied is:

(*dum10=(solω2DT1-dum4)//Together//Simplify*)

The numerator must vanish:

dum11=Numerator[dum10]

Long wave limit:

Normal@Series[dum11,{δz1,0,2}]

Acoustic limit :

Normal@Series[dum11,{ϵa1,0,2}]

EsurfDT

Eq3[x,z,t]
EsurfFS

-Ω Integrate[R[z],{z,-Hi,0}]+kx Integrate[ U[z],{z,-Hi,0}]-I( W[-Hi]-W[0])//Expand

solω2DT1

System solving

Define root function (Findδz)

Remove[Findδz];
Findδz[func_Symbol,δz1min_,δz1max_,Nδz1_]:=
Module[{Deltaz1},
listedz1 = {}; 
Δδz1=N@ ((δz1max-δz1min)/(Nδz1));
ϵδz1=(10^(-6));
maxiter=(10^12);
ϵnum=(10^(-12));
For [i=0,i≤Nδz1,i++,
Quiet[booleen = Check[FindRoot[{func[δz1]==0},{{δz1,δz1min+i Δδz1}},MaxIterations->maxiter],$Failed]];
Quiet[If[Head @ booleen == List && Abs@func[δz1 /. booleen]<ϵnum
,(*Print[δz1 /. booleen,"  ",Abs@func[δz1 /. booleen]];*)
listedz1 = Insert[listedz1, δz1 /. booleen,-1];
(*listedz1 = Insert[listedz1, δz1 /. booleen,-1];*)
,0]];
];
listedz1bis = DeleteDuplicates[IntegerPart[listedz1 (1/(ϵδz1))]ϵδz1]; 
Deltaz1 = N@Sort[DeleteCases[DeleteCases[listedz1bis,x_/;x<δz1min],x_/;x>δz1max]];
Return[Deltaz1];
]

Test function

Remove[functest]
functest[δz1_]=(δz1-2)(δz1-10)
Findδz[functest,1,7,10]

Initialize (or recall) δz functions

(kz^2) from inner dispersion in ocean (1) and sediment (2) layers

kz2p

ϵhp=((Hb-Hi)/Hi);
ϵhip=(Hi/Hb);
δz2p1=(Hi^2)kz2p/.cs-> cs1/.Ω->ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb->Hi/ ϵhi//Simplify//Expand
δz2p2=(((Hb-Hi))^2)kz2p/.cs-> cs2/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb->Hi/ ϵhi//Simplify//Expand

Num1=((δz1)^2)-δz2p1
Num2=((δz2)^2)-δz2p2

(kzc^2) from pure acoustic wave dispersion in ocean (1) and sediment (2) layers

kz2cp=((Ω^2)/(cs^2))-(kx^2);
ϵz2cp=(((kz2cp/.cs-> cs1))/((kz2cp/.cs-> cs2)))/.Ω-> ω Sqrt[(g/Hi)]/. {kz1-> δz1/Hi,kz2-> δz2/(Hb-Hi),kx->δx/Hi,g-> ((ϵa1)^2)(cs1^2)/Hi}/.kzc2->kzc1/ϵzc/.cs2-> cs1 /ϵ12/.(Hb-Hi)-> ϵh Hi/.δz2->δz1 ϵz/.Hb->Hi/ ϵhi//Simplify

Numerical evaluation

Small parameter for plots

ϵn=0.1;

Discretization in δz1

Nδz1=10;

First-stage numerical evaluations:

Remove[fn]
gn=9.81;
Hin=100;
Hbn=2000;
cs1n=1500;
cs2n=1800;
ϵρn=1800/1000;

Second-stage numerical evaluations:

ωn=( Sqrt[Hi/g] 2π fn)/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n}  ;
ϵhn=ϵhp/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n} ;
ϵhin=ϵhip/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n} ;
ϵ12n=(cs1/cs2)/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n} ;
ϵa1n=((Sqrt[g Hi])/cs1)/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n} ;

Rigid Lid

Initialize interface relation

Recall δz in both layers

δz2p1
δz2p2

Relation to be tested (modifier ϵr ici !!!!)

RelRLn/. ((ϵzc)^2)-> ϵz2cp/.ϵρ-> ϵρn;
(*(ω^2)-((((δx)^2)Tan[δz1])/(δz1))*)
%/.ϵz-> ((δz2)/(δz1));
%/.{((δz1)^2)-> δz2p1,((δz2)^2)-> ((ϵr)^2)δz2p2}/.{δz1-> Sqrt@δz2p1,δz2-> ϵr Sqrt@δz2p2};
NumInt=%;

Recall system (3 dispersion relations for 3 unknowns)

Num1
Num2
NumInt

Method δz1

Remove[Numn]
Solδx=Sqrt@δx2/.First@Solve[(Num1/.δx->Sqrt@δx2)==0,δx2];
NumInt/.δx-> Solδx;
%/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n};
Numn[δz1_]=Simplify[%,δz1>0]//Simplify

Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn[δz1]/.fn-> (10^(-1))/.ϵr->-I},{δz1,0,π+ϵ}]

ipmin=-5;
ipmax=5;
dip=1;

RootRL1={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,ϵn,2π+ϵn,Nδz1];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootRL1=Insert[RootRL1,{N@(10^ip),sol[[jj]]},-1];];
]
(*RootRL2={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> I/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,π/2,π,100];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootRL2=Insert[RootRL2,{N@(10^ip),sol[[jj]]},-1];];
]*)

RootRL1;
(*RootRL2*)
Export["/home/dumpa/NUMLAB/svg_txt/RL.txt",RootRL1];
floc=RootRL1[[All,1]];dz=RootRL1[[All,2]];
deltaxFS = Sqrt[-dz^2+Hin^2*(2*π*floc/cs1n)^2];
cphiFs = 2*π*floc/(deltaxFS/Hin);
deltaz2FS = Sqrt[-deltaxFS^2*Hbn^2/Hin^2+Hbn^2*(2*π*floc/cs2n)^2];
Export["/home/dumpa/NUMLAB/svg_txt/RL_deltax.txt",deltaxFS];
Export["/home/dumpa/NUMLAB/svg_txt/RL_cphi.txt",cphiFs];
Export["/home/dumpa/NUMLAB/svg_txt/RL_deltaz2.txt",deltaz2FS];

Show[
ListLogLinearPlot[RootRL1]
(*,ListLogLinearPlot[RootRL2,PlotStyle->Red]*)
]

Search for solutions...

ϵ=0.1;
Plot[{Numn[δz1]/.fn-> (10^(-8))/.ϵr-> 1,Numn2/.fn-> (10^(-8))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-3))/.ϵr-> 1,Numn2/.fn-> (10^(-3))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-2))/.ϵr-> 1,Numn2/.fn-> (10^(-2))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Plot[{Numn2/.fn-> 1/.ϵr-> 1,Numn2/.fn-> 1/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10/.ϵr-> 1,Numn2/.fn-> 10/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 100/.ϵr-> 1,Numn2/.fn-> 100/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 1000/.ϵr-> 1,Numn2/.fn-> 1000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10000/.ϵr-> 1,Numn2/.fn-> 10000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Use 2D contourplot to search for solutions

ContourPlot[Numn2/.ϵr-> 1,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]
ContourPlot[Numn2/.ϵr-> -I,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]

Test particular singularity (not a solution)

Plot[{Numn2/.fn-> (10^(-1))/.ϵr-> 1},{δz1,π/2-ϵ,π/2+ϵ}]
Plot[{Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π/2+ϵ}]

Study solutions (and limits) for δx

toto=Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}
Plot[toto,{fn,0,100}]

Particular solution for δz1

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 2/.ϵr->1)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 2/.ϵr->1
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 2/.ϵr->1*)

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 1/.ϵr->-I)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 1/.ϵr->-I
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 1/.ϵr->-I*)

Rigid Lid (inf)

Initialize interface relation

Recall δz in both layers

δz2p1
δz2p2

Relation to be tested (modifier ϵr ici !!!!)

RelRLninf/. ((ϵzc)^2)-> ϵz2cp/.ϵρ-> ϵρn;
(*(ω^2)-((((δx)^2)Tan[δz1])/(δz1))*)
%/.ϵz-> ((δz2)/(δz1));
%/.{((δz1)^2)-> δz2p1,((δz2)^2)-> ((ϵr)^2)δz2p2}/.{δz1-> Sqrt@δz2p1,δz2-> ϵr Sqrt@δz2p2};
NumInt=%;

Recall system (3 dispersion relations for 3 unknowns)

Num1
Num2
NumInt

Method δz1

Remove[Numn]
Solδx=Sqrt@δx2/.First@Solve[(Num1/.δx->Sqrt@δx2)==0,δx2];
NumInt/.δx-> Solδx;
%/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n};
Numn[δz1_]=Simplify[%,δz1>0]//Simplify

Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn[δz1]/.fn-> (10^(-1))/.ϵr->-I},{δz1,0,π+ϵ}]

ipmin=-5;
ipmax=4;
dip=1(*75*);

RootRL1={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,ϵn,2π+ϵn,Nδz1];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootRL1=Insert[RootRL1,{N@(10^ip),sol[[jj]]},-1];];
]
(*RootRL2={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> I/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,π/2,π,100];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootRL2=Insert[RootRL2,{N@(10^ip),sol[[jj]]},-1];];
]*)

RootRL1
(*RootRL2*)
Export["/home/dumpa/NUMLAB/svg_txt/RL_INF.txt",RootRL1];
floc=RootRL1[[All,1]];dz=RootRL1[[All,2]];
deltaxFS = Sqrt[-dz^2+Hin^2*(2*π*floc/cs1n)^2];
Export["/home/dumpa/NUMLAB/svg_txt/RL_INF_deltax.txt",deltaxFS];

Show[
ListLogLinearPlot[RootRL1]
(*,ListLogLinearPlot[RootRL2,PlotStyle->Red]*)
]

Search for solutions...

ϵ=0.1;
Plot[{Numn[δz1]/.fn-> (10^(-8))/.ϵr-> 1,Numn2/.fn-> (10^(-8))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-3))/.ϵr-> 1,Numn2/.fn-> (10^(-3))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-2))/.ϵr-> 1,Numn2/.fn-> (10^(-2))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Plot[{Numn2/.fn-> 1/.ϵr-> 1,Numn2/.fn-> 1/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10/.ϵr-> 1,Numn2/.fn-> 10/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 100/.ϵr-> 1,Numn2/.fn-> 100/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 1000/.ϵr-> 1,Numn2/.fn-> 1000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10000/.ϵr-> 1,Numn2/.fn-> 10000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Use 2D contourplot to search for solutions

ContourPlot[Numn2/.ϵr-> 1,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]
ContourPlot[Numn2/.ϵr-> -I,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]

Test particular singularity (not a solution)

Plot[{Numn2/.fn-> (10^(-1))/.ϵr-> 1},{δz1,π/2-ϵ,π/2+ϵ}]
Plot[{Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π/2+ϵ}]

Study solutions (and limits) for δx

toto=Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}
Plot[toto,{fn,0,100}]

Particular solution for δz1

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 2/.ϵr->1)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 2/.ϵr->1
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 2/.ϵr->1*)

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 1/.ϵr->-I)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 1/.ϵr->-I
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 1/.ϵr->-I*)

Pressure Vanish (P0)

Initialize interface relation

Relation to be tested

δz2p2/.δx-> Solδx/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n}/.fn-> 1000

RelP0n/. ((ϵzc)^2)-> ϵz2cp /.ϵρ-> ϵρn;
(*(ω^2)-((((δx)^2)Tan[δz1])/(δz1))*)
%/.ϵz-> ((δz2)/(δz1));
%/.{((δz1)^2)-> δz2p1,((δz2)^2)-> ((ϵr)^2)δz2p2}/.{δz1-> Sqrt@δz2p1,δz2-> ϵr Sqrt@δz2p2};
NumInt=%;

Recall system (3 dispersion relations for 3 unknowns)

Num1
Num2
NumInt

Method δz1

Remove[Numn]
Solδx=Sqrt@δx2/.First@Solve[(Num1/.δx->Sqrt@δx2)==0,δx2];
NumInt/.δx-> Solδx;
%/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n};
Numn[δz1_]=Simplify[%,δz1>0]//Simplify

Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn[δz1]/.fn-> (10^(-1))/.ϵr->-I},{δz1,0,π+ϵ}]

ipmin=-5;
ipmax=4;
dip=1;

Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->20//Simplify;
sol=Findδz[Numnn,0.,4π+ϵn,Nδz1];
deltax = Sqrt[-sol^2+Hin^2*(2*π*fn)^2/cs1n^2] /. {fn->20};
kz = sol/Hin
kx = deltax/Hin
kz2=Sqrt[-kx^2*(1+gn^2/(cs2n^2*(2π fn)^2))+(2π fn)^2/cs2n^2] /. {fn->20}
cg = {};
For [i=1.,i<Length[kx]+1,i++,
Δ=cs1n^4 (kkz^2+kkx^2)^2-4 gn^2 /. {kkz->kz[[i]]};
Ω = Sqrt[cs1n^2/2*(kkz^2+kkx^2)+.5*Sqrt[Δ]] /. {kkz->kz[[i]]};
cg = Insert[cg,D[Ω,kkx] /. {kkx->kx[[i]]},-1];];
cg

(*RootP01={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,ϵn,2π+ϵn,Nδz1];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootP01=Insert[RootP01,{N@(10^ip),sol[[jj]]},-1];];
]*)
(*RootP02={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> I/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,π/2,π,100];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootP02=Insert[RootP02,{N@(10^ip),sol[[jj]]},-1];];
]*)

RootP01
(*RootP02*)
Export["/home/dumpa/NUMLAB/svg_txt/P0.txt",RootP01];
floc=RootP01[[All,1]];dz=RootP01[[All,2]];
deltaxFS = Sqrt[-dz^2+Hin^2*(2*π*floc)^2/cs1n^2];
cphiFs = 2*π*floc/(deltaxFS/Hin);
deltaz2FS = Sqrt[-deltaxFS^2*Hbn^2/Hin^2+Hbn^2*(2*π*floc/cs2n)^2];
Export["/home/dumpa/NUMLAB/svg_txt/P0_deltax.txt",deltaxFS];
Export["/home/dumpa/NUMLAB/svg_txt/P0_cphi.txt",cphiFs];
Export["/home/dumpa/NUMLAB/svg_txt/P0_deltaz2.txt",deltaz2FS];

Show[
ListLogLinearPlot[RootP01]
(*,ListLogLinearPlot[RootP02,PlotStyle->Red]*)
]

Search for solutions...

ϵ=0.1;
Plot[{Numn[δz1]/.fn-> (10^(-8))/.ϵr-> 1,Numn2/.fn-> (10^(-8))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-3))/.ϵr-> 1,Numn2/.fn-> (10^(-3))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-2))/.ϵr-> 1,Numn2/.fn-> (10^(-2))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Plot[{Numn2/.fn-> 1/.ϵr-> 1,Numn2/.fn-> 1/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10/.ϵr-> 1,Numn2/.fn-> 10/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 100/.ϵr-> 1,Numn2/.fn-> 100/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 1000/.ϵr-> 1,Numn2/.fn-> 1000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10000/.ϵr-> 1,Numn2/.fn-> 10000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Use 2D contourplot to search for solutions

ContourPlot[Numn2/.ϵr-> 1,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]
ContourPlot[Numn2/.ϵr-> -I,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]

Test particular singularity (not a solution)

Plot[{Numn2/.fn-> (10^(-1))/.ϵr-> 1},{δz1,π/2-ϵ,π/2+ϵ}]
Plot[{Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π/2+ϵ}]

Study solutions (and limits) for δx

toto=Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}
Plot[toto,{fn,0,100}]

Particular solution for δz1

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 2/.ϵr->1)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 2/.ϵr->1
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 2/.ϵr->1*)

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 1/.ϵr->-I)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 1/.ϵr->-I
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 1/.ϵr->-I*)

Pressure Vanish (P0 inf)

Initialize interface relation

Relation to be tested

δz2p2/.δx-> Solδx/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n}/.fn-> 1000

RelP0n
RelP0ninf

RelP0ninf(*/.δz2-> I δz2*)/. ((ϵzc)^2)-> ϵz2cp /.ϵρ-> ϵρn;
(*(ω^2)-((((δx)^2)Tan[δz1])/(δz1))*)
%/.ϵz-> ((δz2)/(δz1));
%/.{((δz1)^2)-> δz2p1,((δz2)^2)-> ((ϵr)^2)δz2p2}/.{δz1-> Sqrt@δz2p1,δz2-> ϵr Sqrt@δz2p2};
NumInt=%;

Recall system (3 dispersion relations for 3 unknowns)

Num1
Num2
NumInt

Method δz1

Remove[Numn]
Solδx=Sqrt@δx2/.First@Solve[(Num1/.δx->Sqrt@δx2)==0,δx2];
NumInt/.δx-> Solδx;
%/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n};
Numn[δz1_]=Simplify[%,δz1>0]//Simplify

Plot[{Numn[δz1]/.fn->20/.ϵr-> 1,Numn[δz1]/.fn-> 20/.ϵr->-I},{δz1,2.,3.}]

ipmin=0;
ipmax=2;
dip=20(*75;*)

Numn[δz1] /. ϵr->-I /. fn-> 20

(*Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->20//Simplify;
sol=Findδz[Numnn,2.1,2.5,5];
deltax = Sqrt[-sol^2+Hin^2*(2*π*fn)^2/cs1n^2] /. {fn->20};
kz = sol/Hin
kx = deltax/Hin
kz2=Sqrt[-kx^2*(1+gn^2/(cs2n^2*(2π fn)^2))+(2π fn)^2/cs2n^2] /. {fn->20}
cg = {};
For [i=1.,i<Length[kx]+1,i++,
Δ=cs1n^4 (kkz^2+kkx^2)^2-4 gn^2 /. {kkz->kz[[i]]};
Ω = Sqrt[cs1n^2/2*(kkz^2+kkx^2)+.5*Sqrt[Δ]] /. {kkz->kz[[i]]};
cg = Insert[cg,D[Ω,kkx] /. {kkx->kx[[i]]},-1];];
cg*)

RootP01={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,ϵn,2π+ϵn,Nδz1];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootP01=Insert[RootP01,{N@(10^ip),sol[[jj]]},-1];];
]
(*RootP02={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> I/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,π/2,π,100];
(*If[Length[sol]>0,
RootP0=Insert[RootP0,{N@(10^ip),Last@sol},-1];];*)
For[jj=1,jj≤ Length[sol],jj++,
RootP02=Insert[RootP02,{N@(10^ip),sol[[jj]]},-1];];
]*)

(*RootP01
(*RootP02*)
Export["/home/dumpa/NUMLAB/svg_txt/P0_INF.txt",RootP01];
floc=RootP01[[All,1]];dz=RootP01[[All,2]];
deltaxFS = Sqrt[-dz^2+Hin^2*(2*π*floc)^2/cs1n^2];
Export["/home/dumpa/NUMLAB/svg_txt/P0_INF_deltax.txt",deltaxFS];*)

Show[
ListLogLinearPlot[RootP01]
(*,ListLogLinearPlot[RootP02,PlotStyle->Red]*)
]

Search for solutions...

ϵ=0.1;
Plot[{Numn[δz1]/.fn-> (10^(-8))/.ϵr-> 1,Numn2/.fn-> (10^(-8))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-3))/.ϵr-> 1,Numn2/.fn-> (10^(-3))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-2))/.ϵr-> 1,Numn2/.fn-> (10^(-2))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Plot[{Numn2/.fn-> 1/.ϵr-> 1,Numn2/.fn-> 1/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10/.ϵr-> 1,Numn2/.fn-> 10/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 100/.ϵr-> 1,Numn2/.fn-> 100/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 1000/.ϵr-> 1,Numn2/.fn-> 1000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10000/.ϵr-> 1,Numn2/.fn-> 10000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Use 2D contourplot to search for solutions

ContourPlot[Numn2/.ϵr-> 1,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]
ContourPlot[Numn2/.ϵr-> -I,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]

Test particular singularity (not a solution)

Plot[{Numn2/.fn-> (10^(-1))/.ϵr-> 1},{δz1,π/2-ϵ,π/2+ϵ}]
Plot[{Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π/2+ϵ}]

Study solutions (and limits) for δx

toto=Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}
Plot[toto,{fn,0,100}]

Particular solution for δz1

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 2/.ϵr->1)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 2/.ϵr->1
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 2/.ϵr->1*)

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 1/.ϵr->-I)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 1/.ϵr->-I
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 1/.ϵr->-I*)

Free Surface

Initialize interface relation (RelIntTest)

Relation to be tested

δz2p2

RelFSn/. ((ϵzc)^2)-> ϵz2cp /.ϵρ-> ϵρn
%/.ϵz-> ((δz2)/(δz1));
%/.{((δz1)^2)-> δz2p1,((δz2)^2)-> ((ϵr)^2) δz2p2}/.{δz1-> Sqrt@δz2p1,δz2-> ϵr Sqrt@δz2p2};
NumInt=%;

Recall system (3 dispersion relations for 3 unknowns)

Num1
Num2
NumInt

Method δz1

Remove[Numn]
Solδx=Sqrt@δx2/.First@Solve[(Num1/.δx->Sqrt@δx2)==0,δx2]
NumInt/.δx-> Solδx;
%/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n};
Numn[δz1_]=Simplify[%,δz1>0]//Simplify

ipmin=-5;
ipmax=4;
dip=200;

Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->20//Simplify;
sol=Findδz[Numnn,ϵn,4π+ϵn,Nδz1];
deltax = Sqrt[-sol^2+Hin^2*(2*π*fn)^2/cs1n^2] /. {fn->20};
kz = sol/Hin
kx = deltax/Hin
kz2=Sqrt[-kx^2*(1+gn^2/(cs2n^2*(2π fn)^2))+(2π fn)^2/cs2n^2] /. {fn->20}
cg = {};
For [i=1.,i<Length[kx]+1,i++,
Δ=cs1n^4 (kkz^2+kkx^2)^2-4 gn^2 /. {kkz->kz[[i]]};
Ω = Sqrt[cs1n^2/2*(kkz^2+kkx^2)+.5*Sqrt[Δ]] /. {kkz->kz[[i]]};
cg = Insert[cg,D[Ω,kkx] /. {kkx->kx[[i]]},-1];];
cg

(*RootFS1={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,ϵn,2π+ϵn,Nδz1];
For[jj=1,jj≤ Length[sol],jj++,
RootFS1=Insert[RootFS1,{N@(10^ip),sol[[jj]]},-1];];
]*)

(*RootFS2={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]= Numn[δz1]/.ϵr->I/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,ϵn,2π+ϵn,Nδz1];
For[jj=1,jj≤ Length[sol],jj++,
RootFS2=Insert[RootFS2,{N@(10^ip),sol[[jj]]},-1];];
]*)

RootFS1
(*RootFS2*)
Export["/home/dumpa/NUMLAB/svg_txt/FS.txt",RootFS1];
floc=RootFS1[[All,1]];dz=RootFS1[[All,2]];
deltaxFS = Sqrt[-dz^2+Hin^2*(2*π*floc)^2/cs1n^2];
cphiFs = 2*π*floc/(deltaxFS/Hin);
deltaz2FS = Sqrt[-deltaxFS^2*Hbn^2/Hin^2+Hbn^2*(2*π*floc/cs2n)^2];
Export["/home/dumpa/NUMLAB/svg_txt/FS_deltax.txt",deltaxFS];
Export["/home/dumpa/NUMLAB/svg_txt/FS_cphi.txt",cphiFs];
Export["/home/dumpa/NUMLAB/svg_txt/FS_deltaz2.txt",deltaz2FS];

Abs@Numn[3.36889`]/.fn-> 0.00001`/.ϵr-> 1

Show[
ListLogLinearPlot[RootFS1]
(*,ListLogLinearPlot[RootFS2,PlotStyle->Red]*)
]

Search for solutions...

ϵ=0.1;
Plot[{Numn[δz1]/.fn-> (10^(-8))/.ϵr-> 1,Numn2/.fn-> (10^(-8))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-3))/.ϵr-> 1,Numn2/.fn-> (10^(-3))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-2))/.ϵr-> 1,Numn2/.fn-> (10^(-2))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Plot[{Numn2/.fn-> 1/.ϵr-> 1,Numn2/.fn-> 1/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10/.ϵr-> 1,Numn2/.fn-> 10/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 100/.ϵr-> 1,Numn2/.fn-> 100/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 1000/.ϵr-> 1,Numn2/.fn-> 1000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10000/.ϵr-> 1,Numn2/.fn-> 10000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Use 2D contourplot to search for solutions

ContourPlot[Numn2/.ϵr-> 1,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]
ContourPlot[Numn2/.ϵr-> -I,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]

Test particular singularity (not a solution)

Plot[{Numn2/.fn-> (10^(-1))/.ϵr-> 1},{δz1,π/2-ϵ,π/2+ϵ}]
Plot[{Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π/2+ϵ}]

Study solutions (and limits) for δx

toto=Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}
Plot[toto,{fn,0,100}]

Particular solution for δz1

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 2/.ϵr->1)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 2/.ϵr->1
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 2/.ϵr->1*)

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 1/.ϵr->-I)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 1/.ϵr->-I
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 1/.ϵr->-I*)

Free Surface (H2 inf)

Initialize interface relation (RelIntTest)

Relation to be tested

δz2p2

RelFSninf/. ((ϵzc)^2)-> ϵz2cp /.ϵρ-> ϵρn
%/.ϵz-> ((δz2)/(δz1));
%/.{((δz1)^2)-> δz2p1,((δz2)^2)-> ((ϵr)^2) δz2p2}/.{δz1-> Sqrt@δz2p1,δz2-> ϵr Sqrt@δz2p2};
NumInt=%;

Recall system (3 dispersion relations for 3 unknowns)

Num1
Num2
NumInt

Method δz1

Remove[Numn]
Solδx=Sqrt@δx2/.First@Solve[(Num1/.δx->Sqrt@δx2)==0,δx2]
NumInt/.δx-> Solδx;
%/.{g-> gn,Hi-> Hin,Hb->Hbn,cs1-> cs1n,cs2-> cs2n,ω-> ωn,ϵh-> ϵhn,ϵhi-> ϵhin,ϵ12-> ϵ12n,ϵa1-> ϵa1n};
Numn[δz1_]=Simplify[%,δz1>0]//Simplify

ipmin=-5;
ipmax=4;
dip=0;

RootFS1={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]=Numn[δz1]/.ϵr-> 1/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn, ϵn,(2π+ϵn),Nδz1];
For[jj=1,jj≤ Length[sol],jj++,
RootFS1=Insert[RootFS1,{N@(10^ip),sol[[jj]]},-1];];
]

(*RootFS2={};
For [ii=0,ii≤dip,ii++,
ip=ipmin+(((ipmax-ipmin))/dip)ii;
Numnn[δz1_]= Numn[δz1]/.ϵr->I/.fn->(10^ip)//Simplify;
sol=Findδz[Numnn,ϵn,2π+ϵn,Nδz1];
For[jj=1,jj≤ Length[sol],jj++,
RootFS2=Insert[RootFS2,{N@(10^ip),sol[[jj]]},-1];];
]*)

(*RootFS1
(*RootFS2*)
(*Export["/home/dumpa/NUMLAB/svg_txt/FS_INF.txt",RootFS1];*)
floc=RootFS1[[All,1]]
dz=RootFS1[[All,2]];
deltaxFS = Sqrt[-dz^2+Hin^2*(2*π*floc)^2/cs1n^2]
(*Export["/home/dumpa/NUMLAB/svg_txt/FS_INF_deltax.txt",deltaxFS];*)*)

Abs@Numn[3.36889`]/.fn-> 0.00001`/.ϵr-> 1

Show[
ListLogLinearPlot[RootFS1]
(*,ListLogLinearPlot[RootFS2,PlotStyle->Red]*)
]

Search for solutions...

ϵ=0.1;
Plot[{Numn[δz1]/.fn-> (10^(-8))/.ϵr-> 1,Numn2/.fn-> (10^(-8))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-3))/.ϵr-> 1,Numn2/.fn-> (10^(-3))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn-> (10^(-2))/.ϵr-> 1,Numn2/.fn-> (10^(-2))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn[δz1]/.fn->(10^(-1))/.ϵr-> 1,Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Plot[{Numn2/.fn-> 1/.ϵr-> 1,Numn2/.fn-> 1/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10/.ϵr-> 1,Numn2/.fn-> 10/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 100/.ϵr-> 1,Numn2/.fn-> 100/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 1000/.ϵr-> 1,Numn2/.fn-> 1000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]
Plot[{Numn2/.fn-> 10000/.ϵr-> 1,Numn2/.fn-> 10000/.ϵr->-I},{δz1,π/2-ϵ,π+ϵ}]

Use 2D contourplot to search for solutions

ContourPlot[Numn2/.ϵr-> 1,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]
ContourPlot[Numn2/.ϵr-> -I,{δz1,π/2-ϵ,π+ϵ},{fn,0,50},PlotRange->{-5 (10^(-1)),5 (10^(-1))}]

Test particular singularity (not a solution)

Plot[{Numn2/.fn-> (10^(-1))/.ϵr-> 1},{δz1,π/2-ϵ,π/2+ϵ}]
Plot[{Numn2/.fn-> (10^(-1))/.ϵr->-I},{δz1,π/2-ϵ,π/2+ϵ}]

Study solutions (and limits) for δx

toto=Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}
Plot[toto,{fn,0,100}]

Particular solution for δz1

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 2/.ϵr->1)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 2/.ϵr->1
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 2/.ϵr->1*)

(*δz1n=δz1/.FindRoot[(Numn/.fn-> 1/.ϵr->-I)==0 ,{δz1,π-ϵ,π/2,π}](*,MaxIterations->1000000*)
Numn/.δz1-> δz1n/.fn-> 1/.ϵr->-I
Solδx/.{δz1-> δz1n,ω-> ωn,ϵa1-> ϵa1n}/.fn-> 1/.ϵr->-I*)

Gather

Select only part of P0 and RL dispersions

RootRL1LF=DeleteCases[RootRL1,{x_,y_}/;x>1];
RootP01HF=DeleteCases[RootP01,{x_,y_}/;x<1];

Show[
ListLogLinearPlot[RootFS1,PlotStyle->Red,Ticks->{Automatic,{(π/2),π,((6π)/4),((8π)/4)}},PlotRange->{{(10^ipmin),(10^(ipmax-1))}, {ϵn,2π+ϵn}},GridLines->{Automatic,{π/4,π/2,3π/4,π,((5π)/4),((6π)/4),((7π)/4),((8π)/4)}}],
(*ListLogLinearPlot[RootFS2,PlotStyle->Green],*)
ListLogLinearPlot[RootRL1LF,PlotStyle->Black,PlotMarkers->"-",PlotRange->{{(10^ipmin),1},Automatic}],
ListLogLinearPlot[RootP01HF,PlotStyle->Blue,PlotMarkers->"-",PlotRange->{{1,(10^ipmax)},Automatic}],
AxesLabel->{f,(Subscript[δz, 1])}
]

(*Show[
ListLogLinearPlot[RootFS2,PlotStyle->Green],
ListLogLinearPlot[RootRL2,PlotStyle->Blue,PlotMarkers->"-"],
ListLogLinearPlot[RootP02,PlotStyle->Black,PlotMarkers->"-"],
PlotRange-> {π/2,π}
]*)