% Emacs mode header: -*-Octave-*- source. %
%function [p_out,xres,zres,rh,eh] = tw_nonlin_4M2(p_inp);
  function [x,P,Q,rh,eh,data_ph] = tw_nonlin_4M2(o_inp,p_inp,g_inp);
tic
1;
%%%clear
%x0=10;
N=1000; % N=10000;
npoint=length(p_inp);

P=zeros(N,npoint);
Q=zeros(N,npoint);

global L; L=10;
global alpha; alpha=1;
global omega0; omega0=1000;
global d; d=2;
global omega; omega=o_inp; %orig: 20, 286.5, 82, 23.51
omega
global gamma; gamma=g_inp; %orig 0;/ -10;
global B; B=1e-7;  %orig B=1e-7;/ 2e-7; 0..error
rho=1;
global mu; mu=rho*omega*omega;

function zdot = func(x0,z) %Matlab - reversed order of arguments...
     %global mu;
     %global d;
     %global gamma;
     %global B;
     A=lin_res_A(x0);
     p2=z(1)*z(1) + z(2)*z(2);
     u=u_nlin(x0,p2);
     Lr = -real(mu/(A+i*B*u+i*gamma));
     Li = -imag(mu/(A+i*B*u+i*gamma));
     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%--------endfunction

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%--------endfunction 

function u = u_nlin(x0,p2)
     %global B;
     %global gamma;
     if (B==0),
      u=0;
     else
        a0=-p2/B^2;
        a1=(lin_res_A(x0)^2+gamma^2)/(B^2);
        a2=2*gamma/B;
     	s=a1/3-a2*a2/9;
     	r=-a0/2+a1*a2/6-a2^3/27;
        u = (r+sqrt(s^3+r^2))^(1/3)+...
        sign(r-sqrt(s^3+r^2))*(abs(r-sqrt(s^3+r^2)))^(1/3);
     end%--if
end%--------endfunction

%npoint=9; %orig 9; 3, 9
p0=[1e-9,5e-9,1e-8,5e-8,1e-7,5e-7,1e-6,5e-6,1e-5];

%k0=omega*sqrt(-rho/(lin_res_A(0)))

for l=1:npoint %orig 1:npoint

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

%pabs0=p0(l); %---
pabs0=p_inp(l);
z0=[0,0,pabs0*q0r,pabs0*q0i]';
%pause %z0=[0,0,0,0]';
%z0
x=linspace(0,L,N)';
%zres=lsode("func",z0,x);
%[xres, zres]=ode23('func',0,L,z0,L/N);
[xres, zres] = ode45(@func,x,z0); %,L/N);

%^error('test')

p=zres(:,1)+i*zres(:,2);
q=zres(:,3)+i*zres(:,4);

for j=1:N
 p2=abs(p(j)*p(j));
 x0=x(j);
 h(j)=p(j)/(lin_res_A(x0)+i*B*u_nlin(x0,p2)+i*gamma);
 u(j)=u_nlin(x0,p2);
  P(j,l)=p(j); %--- P
  Q(j,l)=q(j); %--- Q
end%---for

%sizeh=size(h)

xr=L-d*log(omega0/omega);
n=floor((N-1)*(xr/L)+1); %----------%
p20(l)=abs(p(N)*p(N));
p2x(l)=abs(p(n)*p(n));
h2x(l)=abs(h(n)*h(n));

rh(l,:)=real(h(:)); %--indexing by l
ih(l,:)=imag(h(:));
eh(l,:)=abs(rh(l,:)+i*ih(l,:));

% gset nologscale;
% data = [x',rh];
% gplot data w l

%filename=["data",num2str(l),".dat"];
%myfile=fopen(filename,"w");
%fprintf(myfile,"%10.5e %10.5e \n",rh,ih);
%fclose(myfile);
%%% figure;
%%% plot(-x,rh/10,'b',-x,eh/10,'g', -x,zres,'r');
%%% close
end%---for %l
%rhh=1;
toc
data_ph=[p20;h2x];
%gset logscale xy
%gplot data w line
%%% p20
%%% h2x

end%--- end function main
