function mroots = dispersion_free_surface(alpha,N,h) % mroots = dispersion_free_surface(alpha,N,h) % calculates the positive imaginary and N first positive real solutions of -alpha = k*tan(k h) % for complex alpha. It uses three methods - homotopy for % starting with alpha =1, guess with linear expansion % and a linear expansion. % The first roots is positive imaginary and the next are the first N positive real % ordered from smallest. % % If the value for h is not given the default value is h = 1; % it would be easy to write a much faster program for only real alpha if nargin == 2 h = 1; else alpha = h*alpha; % scaling for the solution end %first of all we calculate the root for alpha = 1; mroots = zeros(1,N+1); if N ==0; % the N = 0 case does not involve any of the special methods and is treated separately. count = 0; mroots(count+1) = homotopy(alpha,count); else count = 0; mroots(count+1) = homotopy(alpha,count); count = count + 1; while 0 <= 1 %first of all we calculate the root for alpha = 1; mroots(count+1) = homotopy(alpha,count); if abs(mroots(count + 1) - (i*count*pi + alpha/(i*count*pi))) < 0.01 %sprintf('Now we can use the close guess') while 0 <=1 mroots(count + 1) = oneroot(alpha,i*count*pi + alpha/(i*count*pi)); %abs(mroots(count + 1) - (i*count*pi + alpha/(i*count*pi))) if abs(mroots(count + 1) - (i*count*pi + alpha/(i*count*pi))) < 1e-8 % sprintf('Now we can work the rest out easily') mroots(count+1:N+1) = i*(count:N)*pi + alpha./(i*(count:N)*pi); count = N; break end if count ==N break end count = count + 1; end end if count == N break end count = count +1; end end mroots = -i/h*mroots;mroots(1) = -mroots(1); function mroot = homotopy(alpha,N) %calculates the Nth root using the homotopy method if N == 0; mroot = oneroot(1,1); else mroot = oneroot(1,i*N*pi); end step =0.043; if abs(alpha) < 1 alphastep = ([1:-step:abs(alpha),abs(alpha)]); else alphastep = ([1:step:abs(alpha),abs(alpha)]); end for k=2:length(alphastep) mroot = oneroot(alphastep(k),mroot); end if angle(alpha) > 0 alphastep = abs(alpha)*exp(i*[0:pi/30:angle(alpha),angle(alpha)]); else alphastep = abs(alpha)*exp(i*[0:-pi/30:angle(alpha),angle(alpha)]); end for k=2:length(alphastep) mroot = oneroot(alphastep(k),mroot); end function out = oneroot(alpha,guess) %calculates the root nearest the root guess. ans1 = guess+1; out = guess; while abs(ans1 - out) > 1e-9 ans1 = out; out = ans1 - f(ans1,alpha)/difff(ans1); end function out = f(z,alpha) out = z*tanh(z) - alpha; function out = difff(z) out = tanh(z) + z*sech(z).^2;