% Emacs mode header: -*-Octave-*- source. %
1;
clear
N=10000;
global L=10;
global alpha=1;
global omega_0=1000;
global d=2;
global omega=20;
global gamma=-0.5;
rho=1;
global mu=rho*omega*omega;
function zdot = func(z,x0)
     global mu;
     global d;
     global gamma;
     A=lin_res_A(x0)+i*gamma;
     Lr = -real(mu/A);
     Li = -imag(mu/A);
     zdot(1)=z(3);
     zdot(2)=z(4);
     zdot(3)=Lr*z(1)-Li*z(2);
     zdot(4)=+Lr*z(2)+Li*z(1);
endfunction

function A = lin_res_A(x0)
     global alpha;
     global omega_0;
     global L;
     global omega;
     global d;
     A=alpha*(omega_0*exp(-(L-x0)/d)-omega);
endfunction     

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

p0abs=1e-5;
z0=[0,0,q0r*p0abs,q0i*p0abs];
x=linspace(0,L,N);
zres=lsode("func",z0,x);

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);
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));
