Upload New File

DV4ORK34_shock_deep.m

+function [u4,err] = DV4ORK34_shock_deep(hx,u,omega,lambda0,lambda1,lambda2,nu0,nu1,nu2,alpha0,alpha1,alpha2,alpha40,alpha41,alpha42,beta0,beta1,beta2,gamma0,gamma1,gamma2,muintexp0,muintexp1,muintexp2,dythe_hilb0,dythe_hilb1,dythe_hilb2,shoalflag,dis,kh0,kh1,kh2,omega_0)
+% adaptive embedded RK4(3) method for solving the variable depth HONLS in the Interaction Picture
+% INPUT
+% hz: current step size
+% u: complex field
+% omega: frequency bins
+% lambda0, 1, 2 are the values of lambda at x, x+hx/2, x + hx
+% mn0, 1, 2 are the values of nu at x, x+hx/2, x + hx
+% alpha0, 1, 2 are the values of alpha (TOD) at x, x+hx/2, x + hx
+% alpha40, 1, 2 are the values of alpha (FOD) at x, x+hx/2, x + hx 
+% beta0, 1, 2 are the values of beta (|A|^2 A_t) at x, x+hx/2, x + hx
+% gamma0, 1, 2 are the values of beta (A^2 A*_t) at x, x+hx/2, x + hx
+%
+% We limit ourselves to the Djordjevic Redekopp shoaling
+% muintexp0, 1, 2 are the values of exponential of the integral of shoaling at x, x+hx/2, x + hx
+% shoalflag allows to switch shoaling off (it could be useful to compare to fibres)
+%
+% BEWARE the IP integrating factor now involves an integral of lambda in
+% the expontential and the multiplication by a factor on account of
+% shoaling.
+%
+% OUTPUT: function value at z + hz, correcting hz0 accordingly
+%
+% Andrea Armaroli 26.02.19 
+% Coefficients come from Sedletsky work. No nonlocality is present in the
+% finite depth case
+%Added hilbert term
+% 02.04.19 Included 4OD to suppress the resonant-radiation
+%%
+persistent k5
+n = size(omega,1);
+mask = ones(n,1);
+% mask(round(6*n/14:8*n/12-1)) = 0;
+
+if shoalflag == 1
+    % integrating factor x -> x + hx/2
+    phaselin1 = muintexp0./muintexp1.*exp(-1i.*(omega.^2.*(lambda0+lambda1)./2 + omega.^3.*(alpha0 + alpha1)/2 + omega.^4.*(alpha40 + alpha41)/2 ).*hx/2);
+    
+    % integrating factor x+ hx/2 -> x + hx
+    phaselin2 = muintexp1./muintexp2.*exp(-1i.*(omega.^2.*(lambda2+lambda1)./2 + omega.^3.*(alpha1+alpha2)/2 + omega.^4.*(alpha41 + alpha42)/2 ).*hx/2);
+
+else 
+     % integrating factor x -> x + hx/2
+    phaselin1 = exp(-1i.*(omega.^2.*(lambda0+lambda1)./2 + omega.^3.*(alpha0 + alpha1)/2 + omega.^4.*(alpha40 + alpha41)/2 ).*hx/2);
+    
+    % integrating factor x+ hx/2 -> x + hx
+    phaselin2 = exp(-1i.*(omega.^2.*(lambda2+lambda1)./2 + omega.^3.*(alpha1+alpha2)/2 + omega.^4.*(alpha41 + alpha42)/2 ).*hx/2);
+
+end
+    
+uIP = fft((phaselin1.*(ifft(u).*mask)));
+if isempty(k5)
+    Ncurr = nonlinear(u,omega, nu0, beta0, gamma0,dis,dythe_hilb0,kh0,omega_0);
+    k1 = fft(phaselin1.*(ifft(Ncurr)));
+else
+    k1 = fft(phaselin1.*(ifft(k5)));
+end
+
+	
+k2 = nonlinear(uIP + hx/2.*k1 ,omega, nu1, beta1, gamma1,dis,dythe_hilb1,kh1,omega_0);
+k3 = nonlinear(uIP + hx/2.*k2 ,omega,  nu1, beta1, gamma1,dis,dythe_hilb1,kh1,omega_0);
+u2 = fft(phaselin2.*(ifft(uIP + hx.*k3)));
+k4 = nonlinear(u2, omega, nu2, beta2, gamma2,dis,dythe_hilb2,kh2,omega_0);
+beta = fft(phaselin2.*ifft(uIP + hx/6.*(k1 + 2*k2 + 2*k3)));
+u4 = beta + hx/6*k4;
+k5 = nonlinear(u4, omega, nu2, beta2, gamma2,dis,dythe_hilb2,kh2,omega_0);
+u3 = beta + hx/30.*(2*k4+3*k5);
+err = norm(u3-u4);
+
+
+end 
+
+function [uNL] = nonlinear(u, omega, nu, beta, gamma,d,dythe_hilb,kh,omega_0)
+
+sigma = tanh(kh);
+k_hil = omega_0.^2./sigma;
+h_hil = kh.*sigma./omega_0.^2;
+% group and phase speed
+cp = omega_0./k_hil;
+cg = (sigma + kh.*(1-sigma.^2))/2./omega_0;
+D_hilb_num_1=2+(1-sigma.^2).*cg./cp;
+D_hilb_den_1=kh-sigma.*cg.^2./(cp.^2);
+D_hilb_general_1=kh.*(D_hilb_num_1./D_hilb_den_1).*omega_0./(4.*sigma);
+	
+% if nu<0
+%     NLflag = 0;
+% else 
+    NLflag = 1;
+	gen_hilb_flag = 1;
+% end
+ee=1
+uNL3O = nu*abs(u).^2.*u;
+uNL4O = beta.*abs(u).^2.*fft(-1i.*omega.*ifft(u)) + gamma.*u.^2.*fft(-1i.*omega.*ifft(conj(u)))+dythe_hilb.*i.*u.*(fft(abs(omega).*ifft(abs(u).^2)))   ;
+uNL = -NLflag*(1i.*uNL3O - uNL4O)- d.*u;
+end