hurricane='isabel';
rep='meanless\temp\filtered\2000\gridded_data\corre\';
if strcmp(hurricane,'frances')==1
    petit=2;
    grand=3;
    el=3;
else
    petit=1;
    grand=2;
    el=1.5
end
u=[];
v=[];
dt=[];
r=[];
times=[];
h=[];
for boxe=1:3
    load(['C:\Documents and Settings\sylvie lorsolo\My Documents\research\',hurricane,'\data\modified\',rep,'corre',num2str(boxe),'.mat']);
    if boxe==1 & strcmp(hurricane,'frances')==1
        x2=4500;  %Isabel       -4000
        x1=2000;   % isabel      -6000
        y2= 2500;   %isabel      -1500
        y1=1000;   % isabel        -3500
        
    elseif boxe==2 & strcmp(hurricane,'frances')==1
        x2=6000;   %isabel   -2000
        x1=4000;  %isabel   -4000
        y2= 2500;  %isabel -5000
        y1=1000;   %isabel  -6500
    elseif boxe==3 & strcmp(hurricane,'frances')==1
        x2=4000;    %isabel   -3000 
        x1=2000;    %isabel    -5000
        y2=4500;    %isabel     -4000
        y1=3000;   %isabel   -5500
        
    elseif boxe==1 & strcmp(hurricane,'isabel')==1
        x2= -4000;
        x1=-6000;
        y2= -1500;
        y1= -3500;
        
    elseif boxe==2 & strcmp(hurricane,'isabel')==1
        x2=-2000;
        x1=-4000;
        y2= -5000;
        y1=-6500;
    elseif boxe==3 & strcmp(hurricane,'isabel')==1
        x2=-3000;
        x1=-5000;
        y2=-4000;
        y1=-5500;
    end
    x_mil2=((x2-x1)/2)+x1;
    y_mil2=((y2-y1)/2)+y1;
    h2=sin(deg2rad(el))*sqrt(x_mil2^2 + y_mil2^2)
    
    for i=1:length(corre)
        for j=i+petit:i+grand % for isabel i+1:i+2, for frances i+2:i+3
            u=[u;corre(i).v_x(j)];
            v=[v;corre(i).v_y(j)];   
            dt=[dt;corre(i).delta_t(j)];
            r=[r;corre(i).r(j)];
            times=[times;corre(i).time(j,:)];
            x_mil1=x_mil2+corre(i).v_x(j);
            y_mil1=y_mil2+corre(i).v_y(j);
            h1=sin(deg2rad(el))*sqrt(x_mil1^2 + y_mil1^2);
            h_mil=(abs(h2-h1)/2)+min(h2,h1);
            h=[h;h_mil];
        end
    end
end
t=datenum(times);
[temp1,temp2]=sort(t);
u=u(temp2);
v=v(temp2);
dt=dt(temp2);
r=r(temp2);
h=h(temp2);
t=temp1;
ws=sqrt((u.^2)+(v.^2))./dt;

%--------------------------------------- Wind direction computation---------%
[wd]=wind_direction_deg(u,v);

% ---------------------fin wind direction computation------------------%



%----------------Quality control---------------%
u_corre=find(r>0.5);
u=u(u_corre);
v=v(u_corre);
ws=ws(u_corre);
wd=wd(u_corre);
t=t(u_corre);
dt=dt(u_corre);
h=h(u_corre);

mu_1=mean(u);
mu_2=mean(v);
sigma1=std(u);
sigma2=std(v);
outliers = abs(u - mu_1) > 1.5*sigma1;
u_out=find(outliers==1);
u(u_out)=NaN;
v(u_out)=NaN;
ws(u_out)=NaN;
t(u_out)=NaN;
dt(u_out)=NaN;
wd(u_out)=NaN;
h(u_out)=NaN;
outliers = abs(v - mu_2) > 1.5*sigma2;
u_out=find(outliers==1);
u(u_out)=NaN;
v(u_out)=NaN;
ws(u_out)=NaN;
t(u_out)=NaN;
dt(u_out)=NaN;
wd(u_out)=NaN;
h(u_out)=NaN;
u_NaN=find(isnan(t));
u(u_NaN)=[];
v(u_NaN)=[];
ws(u_NaN)=[];
t(u_NaN)=[];
dt(u_NaN)=[];
wd(u_NaN)=[];
h(u_NaN)=[];
mu=mean(ws);
sigma=std(ws);
outliers = abs(ws - mu) > 1.5*sigma;
u_out=find(outliers==1);
u(u_out)=[];
v(u_out)=[];
ws(u_out)=[];
t(u_out)=[];
wd(u_out)=[];
dt(u_out)=[];
h(u_out)=[];
%------------------------Fin quality control------------%
[h_mean,ws_mean,wd_mean,u_mean,v_mean]=mean_profile(hurricane);
diff_ang=[];
for ii=1:length(h)
    diff1=h(ii)-h_mean;
    u_up=find(diff1<=0);
    h_max=h_mean(u_up(1));
    u_max=u_mean(u_up(1));
    v_max=v_mean(u_up(1));
    u_low=find(diff1>=0);
    h_min=h_mean(u_low(length(u_low)));
    u_min=u_mean(u_low(length(u_low)));
    v_min=v_mean(u_low(length(u_low)));
    input_ord_u=[u_min u_max];
    input_ord_v=[v_min v_max];
    input_abs=[h_min h_max];
    hi=[h_min;h(ii);h_max];
    
    ui = interp1(input_abs,input_ord_u,hi);
    vi=  interp1(input_abs,input_ord_v,hi);
    [wdi]=wind_direction_deg(ui,vi);
    %temp(ii)=vi(2);
    diff_ang=[diff_ang;wdi(2)-wd(ii)];
end
    u_diff_ang=find(diff_ang > 0 & diff_ang >180);
    diff_ang(u_diff_ang)=diff_ang(u_diff_ang)-360;
    u_diff_ang=find(diff_ang < 0 & diff_ang >180);
    diff_ang(u_diff_ang)=diff_ang(u_diff_ang)+360;

