% flux_ADCP2.m
% Calcul du flux suivant la methode 1 : à travers la surface orthogonale à
% la direction principale du courant.
% Possibilité de faire le calcul suivant la methode 2 : à travers la section bateau.
% Pour calculer la direction principale de la veine de courant il y a 4
% facons (choixs):
%          1/ en sélectionnant les courants à gauche de la section
%          2/ en sélectionnant les courants ouest (u<0)
%          3/ en sélectionnant les courants Nord (v>0)
%          4/ aucun critère : tous les courants de la section
% fait appel a :
%   - umod.m
%   - direction_princ_ADCP2.m
%   - anglesection.m
%   - calculflux_ADCP2.m
% angtr angle en deg
% FF  flux en Sv (1Sv=10^6 m3/s)
% SS  surface en km2
% FU flux surfacique en Sv.km^-2
% II incertitude liée au clacul de flux

clear all
ici=pwd;

% recherche des courants superieurs a une certaine valeur pour calculer le flux
%val=input('A partir de quelle vitesse selectionne-t-on les courants? val :');
val=0.05

% chargement de u et v 
%cd d:\golts\Survey\Survey291102\ADCPTethys\tronca\
cd d:\golts\Survey\Survey_0602_anne
%cd e:\Gatti2\golts\Survey\Survey110602\ADCPTethys\tronca
%cd d:\golts\Valor\Valor070405\ADCPTethys\transect
load g0602tr11t2a.mat
%load v0405tr2t9a.mat
u=uv(1:2:end-1,:);
v=uv(2:2:end,:);
zc=flipud(zc);

% calcul de (u,v) en (r,tt) avec tt /Nord clockwise (deg)
cd (ici)
[r,tt]=umod(u,v);

% choix d'un critère pour sélectionner la veine de courant
choixs=menu('Voulez vous sélectionner la veine de courant avec le critère?','courant à gauche de la surface','u<0','v>0','aucun','courant à gauche pour dir');

% recherche des profondeurs (<0) correspondant au CN
profmin=input('profmin (<0) :');
profmax=input('profmax (<0) :');
indzc=find(zc<=profmin & zc>=profmax);  % indices des profs 

choixl=menu('Voulez vous la largeur du courant en :','latitude','longitude','nEnsemble');
if choixl==1
    % recherche des latitudes correspondant au CN
    latmin=input('latmin :');
    latmax=input('latmax :');
    indt=find(xyt(:,2)<=latmax & xyt(:,2)>=latmin);
elseif choixl==2
    % recherche des longitudes correspondant au CN
    lonmin=input('lonmin :');
    lonmax=input('lonmax :');
    indt=find(xyt(:,1)<=lonmax & xyt(:,1)>=lonmin);
elseif choixl==3
    % recherche des numéros d'ensemble correspondant au CN
    indmin=input('indmin :');
    indmax=input('indmax :');
    indt=(indmin:indmax)';
end

choix=input('Quel calcul de flux?\n 1/a travers la surface orthogonale au courant principal\n ou 2/a travers la section bateau\n choix :');
if choix==1
    % CALCUL DE LA DIRECTION PRINCIPALE DU COURANT = DIRECTION DU FLUX
    %pour ce calcul on choisit des courants dont la vitesse est superieure a Vmin
    Vmin=0.1;
    direction_princ_ADCP2
    
    %pour que direction comprise dans [0,180] et que le flux soit positif ou negatif suivant
    %la convention
    if direction>180 & direction<=360
        direction=direction-180;
    end
    
    %angle de la surface normale au courant par rapport au Nord (en deg)
    ang=direction-90;
    
elseif choix==2
    dtest=1; % pour pouvoir faire choixs=1 : courants à gauche 
    % CALCUL DE LA DIRECTION DU FLUX
    % calcul de ang = angle de la surface normale à la direction du flux / Nord (en deg)
    % ici : ang = angle de la section totale /Nord (en deg)
    cd d:\cluster\matlab
    [ang,Dlon,Dlat]=anglesection(xyt(indt(1),1),xyt(indt(1),2),xyt(indt(end-1),1),xyt(indt(end-1),2));
    cd (ici)
    
    % direction du flux
    direction=ang+90;
end

% CALCUL DU FLUX
unit=111190;   %car 1 degre de latitude = 111.19km
dz=4;   %metres pour GOLTS
FF=0;
SS=0;
FU=0;
II=0;
nf=0;
for i=1:length(indt)-1
    clear ind2 indz inter ind1
    ind2=indt(i);
    %angle vitesse/direction du flux (en rad)
    alpha=pi/180*(tt(:,ind2)-direction);
    if choixs==1 
       indz=find(cos(alphasect(:,ind2))<0 & r(:,ind2)>val);
    elseif choixs==2
       indz=find(u(:,ind2)<0 & r(:,ind2)>val);
    elseif choixs==3
       indz=find(v(:,ind2)>0 & r(:,ind2)>val);
    elseif choixs==4
       indz=find(r(:,ind2)>val);   
    end
    %indz=find(u(:,ind2)>0 & r(:,ind2)>val); pour flux du Rhone à SOFI
    inter=find(indz>=indzc(1) & indz<=indzc(end));
    indz=indz(inter);
    for p=1:length(indz)
        ind1=indz(p);
        if isnan(tt(ind1,ind2))==0 & isnan(r(ind1,ind2))==0
            %angle du transect local (en deg) /Nord
            cd d:\cluster\matlab
            [angtr,Dloni,Dlati]=anglesection(xyt(ind2,1),xyt(ind2,2),xyt(ind2+1,1),xyt(ind2+1,2));
            cd (ici)
            calculflux_ADCP2
            nf=nf+1;
        end
    end %boucle en profondeur sur p
end %boucle en espace sur i


% ENREGISTREMENT DES DONNEES
if choixl==1      %latitude
    save ../../flux/flux_g0602tr11t2ast.dat choix choixs val direction profmin profmax latmin latmax FF SS FU II -ascii
elseif choixl==2  %longitude
    save ../../flux/flux_g0602tr11t2ast.dat choix choixs val direction profmin profmax lonmin lonmax FF SS FU II -ascii    
elseif choixl==3  %nEnsemble
    save ../../flux/flux_g0602tr11t2ast.dat choix choixs val direction profmin profmax indmin indmax FF SS FU II -ascii
end
