clear all
rhoe=1000;%mass volumique de l'eau
rhoa=1.2;%masse volumique de l'air
sa=900;%surface de référence des oiles fois le Cz des voile
se=.1;% surface de référence de la structure immergée
V=20; %vitesse du vent absolu
cx=.1;% trainée des foils
%cx=.3; % trainée des flotteurss
%se=1;%surface des flotteurss
%=========
npvit=600; %nombre de points pour la vitesse du bateau
vitmin=.1*V; %vitesse mini
vitmax=3.*V;%vitesse maxi
dx=(vitmax-vitmin)/npvit; % pas de vitesse
%========
nptheta=3; % nombre de points dans la discretisation en allure
thetamin=.3; %allure mini en radian
thetamax=1.4; % allure maxi en radian

dtheta=(thetamax-thetamin)/nptheta;% pas en allure
%===========
npalpha=600;% nombre de points en bordage
alphamin=.1;% bordage mini
alphamax=thetamax;% bordage maxi avant d'être à contre

dalpha=(alphamax-alphamin)/npalpha;% pas en bordage

%=======
f=zeros(npalpha,npvit);%fonction Fx-Tx
ord=linspace(alphamin,alphamax,npalpha);
ab=linspace(vitmin,vitmax,npvit);
figure
%============
for i=1:nptheta
    figure
    theta=thetamin+dtheta*(i-1);
    
    for k=1:npvit
        x=(vitmin+dx*(k-1))/V;
        ab(k)=x*V;
        for j=1:npalpha
        alpha=min(theta,alphamin+dalpha*(j-1));
        ord(j)=alpha;
        cosmua=max(0.,(sin(theta-alpha)-x*sin(alpha))/(sqrt(1+2*x*cos(theta)+x^2)));
        cnumua=1.2*cosmua^2;% coefficient de propulsion vélique
        xi=rhoa*sa*cnumua*sin(alpha)/(rhoe*se*cx);%poussée vélique
        f(j,k)=1*min(.00001,max(-.00001,(xi-1)*x^2+2*x*cos(theta)*xi));
        end
    end
    
    surf(ab,ord,f)
    shading interp
    view([0,90])
    hold on
    title(['Vitesse du vent: V=',num2str(V),'  Angle amure theta=',num2str(theta)])
    xlabel('vitesse u du bateau');
    ylabel('alpha<theta qui est le réglage du plan de voilure');
end