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

global N; N=1000;

omega=82; %11.1, 20, 23.51, 82, 286, gamma=10, 0.1
%[x,zres,ih,rh,ah] = tw_lin_M(omega);
[x,P,Q,ih,rh,eh,data_ph] = tw_nonlin_4M4(23.51,1e-9,0.1);
%(82,5.4556e-26,1);

ah=eh;
%tw_nonlin_4
yy=abs(rh);
F_int_convex_M
F_lambda_M
recc=1000./(3.5*omega.*lambda(:));  %(1:8) -> (:)
latcum=cumtrapz(-x(ilist(:)),recc);
NL_1 = length(ilist)-1;
NL_2 = length(ilist)-2;
hold on
grid on
%plot(-x,sqrt(rh.*rh+ih.*ih)/40,'g')
%ccc=1000*500; gamma=10;
ccc=500;
plot(-x,ah*ccc,'b') %20: /40, 23.. / 800
plot(-x,-ah*ccc,'b') %20: /40, 23.. / 800
plot(-x,rh*ccc,'g') %20: /40, 23.. / 800
%figure
plot(-x(ilist(1:NL_2)),3.5*omega.*lambda(1:NL_2)/50,'r') %(1:8) -> (:)
plot(-x(ilist(1:NL_2)),3.5*omega.*lambda(1:NL_2)/50,'r+')
%plot(-x(ilist(:)),recc,'m')   %(1:8) -> (:), etc
%plot(-x(ilist(:)),recc,'mo')
plot(-x(ilist),zeros(NL,1),'ko') %y-axis was "list"
plot(-x(ilist(:)),latcum,'k')
plot(-x(ilist),latcum,'kx')
%title(' Numerical calculation of latency and velocity')
axis([-10 0 -25 +100])
%axis([-10 -5 -25 25])