You need to use element-wise operations in several places in ‘integrand’ and ‘grsp_fun’. With those corrected, your code is:
% input value
v=1520; % Phase speed
n=1; % Mode number
f=30;
% Pekeris waveguide parameters
H=100; % Source & Receiver Depth (m)
cb=1800; % Half-space speed (m/s)
cw=1500; % Water Column isospeed (m/s)
m=2.2; % Density ratio
v1=1560;v2=1570;
integrand=@(v) (dept_fun(v,n,f,2.2,1800,1500).*grsp_fun(v,n,dept_fun(v,n,f,2.2,1800,1500),2.2,1800,1500))./( v-grsp_fun(v,n,dept_fun(v,n,f,2.2,1800,1500),2.2,1800,1500).*sqrt(v2.^2-v.^2));
ycs1=2*cot(0.05)*integral(@(v)integrand(v),v1,v2)
function dept=dept_fun(v,n,f,m,cb,cw)
dum1=sqrt(cw.^(-2)-v.^(-2));
dum2=sqrt((cw.^(-2)-v.^(-2))/(v.^(-2)-cb.^(-2)));
dept=1./(2*pi*f*dum1)*(pi*n-atan(m*dum2));
end
function grsp=grsp_fun(v,n,H,m,cb,cw)
dum2_1=sqrt((cw.^(-2)-v.^(-2))/(v.^(-2)-cb.^(-2)));
dum1=1./(v.^(2)*cw.^(-2)-1);
dum_2_up=m.*(cw.^(-2)-cb.^(-2)).*(v.^2.*cw.^(-2)-1).^(-.5).*(-(v.^2).*cb.^-2+1).^(-.5);
dum_2_low=(pi*n-atan(m*dum2_1))*(m.^2*cw.^-2-cb.^-2-(m.^2-1)*v.^-2);
dum_all=dum1+dum_2_up./dum_2_low;
grsp=v.*(1+(dum_all).^-1).^-1;
end
However it throws this Warning:
Warning: Minimum step size reached near x = 1570. There may be a singularity, or the tolerances may be
too tight for this problem.
eventually returning:
ycs1 =
-2.920669390162302e+02
