%clear all

varrho=1.3;
S=1.;
alpha=.9;
beta=.3;
alpha0=3;
NPdz=200;
V=30.;
x=linspace(-30,30,NPdz);
pb=linspace(1,2,NPdz);
ener=linspace(1,2,NPdz);
for i=1:NPdz
    alphaa=-atan(x(i)/V);
    qc=2*x(i)*funccz(alpha,alpha0,alphaa)-(V^2+x(i)^2)*(funcdcz(alpha,alpha0,alphaa)+funccx(beta,alphaa));
    qs=2*x(i)*funccx(beta,alphaa)+(V^2+x(i)^2)*(funccz(alpha,alpha0,alphaa)-funcdcx(beta,alphaa));
    pb(i)=(V*varrho*S/2)*(qc*cos(alphaa)+qs*sin(alphaa));
    ener(i)=((varrho*S*V^2)/2)*(sin(alphaa)*funccx(beta,alphaa)+cos(alphaa)*funccz(alpha,alpha0,alphaa));
    if(i == 0)
         ener(i)=1;
    else
        ener(i)=200.*ener(i)/x(i);
    end
end
figure
plot(x,pb,x,ener)
hold on
plot([-25 -25],[-40000 40000],'-')
hold on
plot([25 25],[-40000 40000],'-')
xlabel('vitesse en z, (domaine quasi-statique entre les barres)')
ylabel('criteres pb(bleu) et ener(rouge)')
title('Zones: D+,D-,E+ et E-')