% Emacs mode header: -*-Octave-*- source. %

function [xres,zres,ih,rh,ah] = tw_lin_M(omega);
% clear all
global N; N=10000; %10000
global L; L=10;
global alpha; alpha=1;
global omega0; omega0=1000;
global d; d=2;
global omega; %omega=11.1; %20/ 286.5, 82, 23.51/ 11
global gamma; gamma=-0.5; %-0.5, 0, 
rho=1;
global mu; mu=rho*omega*omega;

function zdot = func(x0,z) %Matlab
     %global mu;
     %global d;
     %global gamma;
     A=lin_res_A(x0)+i*gamma;
     Lr = -real(mu/A);
     Li = -imag(mu/A);
     zzd(1)=z(3);
     zzd(2)=z(4);
     zzd(3)=+Lr*z(1)-Li*z(2);
     zzd(4)=+Lr*z(2)+Li*z(1);
      zdot = [zzd(1); zzd(2); zzd(3); zzd(4);];
end %--- end function

function A = lin_res_A(x0)
     %global alpha;
     %global omega0;
     %global L;
     %global omega;
     %global d;
     A=alpha*(omega0*exp(-(L-x0)/d)-omega);
end %--- end function     

k0=omega*sqrt(rho/(lin_res_A(0)));
q0r=real(k0) %+
q0i=imag(k0) %+

p0abs=1e-5; %1e-5;
pdabs=5e-6; %new
%z0=[q0r*pdabs,0*pdabs,q0r*p0abs,q0i*p0abs];
z0=[0,0,q0r*p0abs,q0i*p0abs];
x=linspace(0,L,N);
%zres=lsode("func",z0,x);
[xres, zres] = ode45(@func,x,z0); %,L/N);

p=zres(:,1)+i*zres(:,2);
q=zres(:,3)+i*zres(:,4);
h=p./(lin_res_A(x)'+i*gamma);

rp=zres(:,1);
ip=zres(:,2);
rq=zres(:,3);
iq=zres(:,4);

norm=abs(h(N));
rh=real(h);
ih=imag(h);
ah=abs(h);
data=[x',ah];

% gset logscale y;
% gplot data w line

p2=abs(p(N)*p(N))
h2=abs(h(N)*h(N))

end %--- end function main