function [minor, major, major_tilt,rotation]=hodograph_sat(t,u,v)

% Sarah Gille's code to least-square fit wind data to an ellipse
% (modifieded)

date=datevec(t(:,1));


edges=datenum([0 0 0 0 0 0]):abs(datenum([0 0 0 0 0 0])-datenum([0 0 0 1 0 0])):datenum([0 0 0 23 59 59]);
edges=[edges 0.999999]; % bin edges every hour
[axxx,bxxx]=histc(t(:,1)-floor(t(:,1)),edges);

for m=1:6;%12; % 1-12 months
    ind=find(date(:,2)==m);
    um=u(ind);
    vm=v(ind);
    tm=t(ind);
    bx=bxxx(ind);
    u_data(:,m)=accumarray([bx],um,[],@nanmean); % daily variance
    v_data(:,m)=accumarray([bx],vm,[],@nanmean);
    u_var(:,m)=accumarray([bx],um,[],@nanstd); % daily variance
    v_var(:,m)=accumarray([bx],vm,[],@nanstd);
    time(:,m)=accumarray([bx],tm-floor(tm),[],@nanmean);
    clear ind um vm bx
    clear ind um vm bx
end


w=2*pi/24; % diurnal frequeny

u_data=u_data(u_data(:,1)~=0,:);
v_data=v_data(v_data(:,1)~=0,:);
u_var=u_var(u_var(:,1)~=0,:);
v_var=v_var(v_var(:,1)~=0,:);
time=time(time(:,1)~=0,:);
%%%defining matrix A and determining least squares parameters
%xx=find(u_var(1,:)>0); % days where there is at least 2 measurements


%%%defining matrix A and determining least squares parameters
xx=find(u_var(1,:)>0 & u_var(2,:)>0 & u_var(3,:)>0 & u_var(4,:)>0 & u_var(5,:)>0 & u_var(6,:)>0);
n=6;
for i=1:length(xx);
Amat=[ones(n,1) cos(2*pi/24.*time(:,xx(i))) sin(2*pi/24.*time(:,xx(i)))]; 
Amat=Amat./(u_var(:,xx(i))*ones(1,3));
Ainv=inv(Amat'*Amat)*Amat';
u_fit=Ainv*(u_data(:,xx(i))./u_var(:,xx(i)));
u_err=sqrt(diag(inv(Amat'*Amat)));
Amat=[ones(n,1) cos(2*pi/24.*time(:,xx(i))) sin(2*pi/24.*time(:,xx(i)))]; 
Amat=Amat./(v_var(:,xx(i))*ones(1,3));
Ainv=inv(Amat'*Amat)*Amat';
v_fit=Ainv*(v_data(:,xx(i))./v_var(:,xx(i)));
v_err=sqrt(diag(inv(Amat'*Amat)));

for j=1:1000
Amat=[ones(n,1) cos(2*pi/24.*time(:,xx(i))) sin(2*pi/24.*time(:,xx(i)))]; 
Amat=Amat./(u_var(:,xx(i))*ones(1,3));
Ainv=inv(Amat'*Amat)*Amat';
ud=u_data(:,xx(i))+randn(n,1).*u_var(:,xx(i));
vd=v_data(:,xx(i))+randn(n,1).*v_var(:,xx(i));
u_fe(:,j)=Ainv*(ud./u_var(:,xx(i)));
u_ee(:,j)=sqrt(diag(inv(Amat'*Amat)));
Amat=[ones(n,1) cos(2*pi/24.*time(:,xx(i))) sin(2*pi/24.*time(:,xx(i)))]; 
Amat=Amat./(v_var(:,xx(i))*ones(1,3));
Ainv=inv(Amat'*Amat)*Amat';
v_fe(:,j)=Ainv*(vd./v_var(:,xx(i)));
v_ee(:,j)=sqrt(diag(inv(Amat'*Amat)));
end

%%%determining ellipse value for ++ and --
 u_bar=u_fit(1);
 v_bar=v_fit(1);
 a=u_fit(2);
 b=u_fit(3);
 c=v_fit(2);
 d=v_fit(3);
 gamma= b.^2 + d.^2 - a.^2 - c.^2;
 beta= -2.*(a.*b+c.*d);
 u_bare=u_fe(1,:);
 v_bare=v_fe(1,:);
 ae=u_fe(2,:);
 be=u_fe(3,:);
 ce=v_fe(2,:);
 de=v_fe(3,:);
 gammae= be.^2 + de.^2 - ae.^2 - ce.^2;
 betae= -2.*(ae.*be+ce.*de);

% --- phase information (compute t_2 standard deviation relative to mean)
 t_1 = atan2(beta, gamma)/(2*w);
 t_1 = mod(t_1+12,12);
 t_2 = atan2(-beta,-gamma)/(2*w);
 t_2 = mod(t_2+12,12);

 t_1e = atan2(betae, gammae)/(2*w);
 t_1e = mod(t_1e+12,12);
 t_2e = atan2(-betae,-gammae)/(2*w);
 t_2e = mod(t_2e+12,12);

 time_err(xx(i))=std(rem(t_2e-t_2+6,12)-6);
 time_out(xx(i))=t_2;


% ---- major and minor axis and errors
 minor(xx(i)) = sqrt((a.*cos(w*t_1) + b.*sin(w*t_1)).^2 + (c.*cos(w*t_1) + d.*sin(w*t_1)).^2);
 major(xx(i)) = sqrt((a.*cos(w*t_2) + b.*sin(w*t_2)).^2 + (c.*cos(w*t_2) + d.*sin(w*t_2)).^2);

 minore = sqrt((ae.*cos(w*t_1e) + be.*sin(w*t_1e)).^2 + (ce.*cos(w*t_1e) + de.*sin(w*t_1e)).^2);
 majore = sqrt((ae.*cos(w*t_2e) + be.*sin(w*t_2e)).^2 + (ce.*cos(w*t_2e) + de.*sin(w*t_2e)).^2);
 minor_error2(xx(i))=std(minore);
 major_error2(xx(i))=std(majore);

da=u_err(2);
db=u_err(3);
dc=v_err(2);
dd=v_err(3);
minor_error_a=(2*a.*(cos(w*t_1)).^2+2*b.*cos(w*t_1).*sin(w*t_1)).*da;
minor_error_b=(2*b.*(sin(w*t_1)).^2+2*a.*cos(w*t_1).*sin(w*t_1)).*db;
minor_error_c=(2*c.*(cos(w*t_1)).^2+2*d.*cos(w*t_1).*sin(w*t_1)).*dc;
minor_error_d=(2*d.*(sin(w*t_1)).^2+2*c.*cos(w*t_1).*sin(w*t_1)).*dd;

major_error_a=(2*a.*(cos(w*t_2)).^2+2*b.*cos(w*t_2).*sin(w*t_2)).*da;
major_error_b=(2*b.*(sin(w*t_2)).^2+2*a.*cos(w*t_2).*sin(w*t_2)).*db;
major_error_c=(2*c.*(cos(w*t_2)).^2+2*d.*cos(w*t_2).*sin(w*t_2)).*dc;
major_error_d=(2*d.*(sin(w*t_2)).^2+2*c.*cos(w*t_2).*sin(w*t_2)).*dd;

minor_error(xx(i))=(sqrt(minor_error_a.^2 + minor_error_b.^2 + minor_error_c.^2 + minor_error_d.^2))./(2*minor(xx(i)));
major_error(xx(i))=(sqrt(major_error_a.^2 + major_error_b.^2 + major_error_c.^2 + major_error_d.^2))./(2*major(xx(i)));

% ---- tilt
 minor_tilt(xx(i)) = atan2(c.*cos(w*t_1)+d.*sin(w*t_1), a.*cos(w*t_1)+b.*sin(w*t_1));
 major_tilt(xx(i)) = atan2(c.*cos(w*t_2)+d.*sin(w*t_2), a.*cos(w*t_2)+b.*sin(w*t_2));

 minor_tilte = atan2(ce.*cos(w*t_1e)+de.*sin(w*t_1e), ae.*cos(w*t_1e)+be.*sin(w*t_1e));
 major_tilte = atan2(ce.*cos(w*t_2e)+de.*sin(w*t_2e), ae.*cos(w*t_2e)+be.*sin(w*t_2e));
 major_tilt_err(xx(i))=std(rem(major_tilte-major_tilt(xx(i))+3*pi,2*pi)-pi);
 minor_tilt_err(xx(i))=std(rem(minor_tilte-minor_tilt(xx(i))+3*pi,2*pi)-pi);

% ----- reversal
 reversal_minor(xx(i))= (u_bar*2.*(a.*cos(w*t_1)+b.*sin(w*t_1))+v_bar*2.*(c.*cos(w*t_1)+d.*sin(w*t_1)));
 reversal_major(xx(i))= (u_bar*2.*(a.*cos(w*t_2)+b.*sin(w*t_2))+v_bar*2.*(c.*cos(w*t_2)+d.*sin(w*t_2)));

% ----- rotation
theta_zero=atan2(c,a);
theta_six=atan2(d,b);
rindex=rem(theta_six-theta_zero+3*pi,2*pi)-pi;
rotation(xx(i))=1;    %%%Counter clockwise
if(rindex==0)rotation(xx(i))=0; end
if(rindex<0)rotation(xx(i))=-1; end   %%%Clockwise

theta_zero=atan2(ce,ae);
theta_six=atan2(de,be);
rindex=rem(theta_six-theta_zero+3*pi,2*pi)-pi;
rotatione=ones(100,1);  %%%Counter clockwise
z=find(rindex==0); if(length(z)>=1) rotatione(z)=0; end
z=find(rindex<0); if(length(z)>=1) rotatione(z)=-1; end %%% Clockwise

rotation_e(xx(i))=sum(rotatione)/100;

end



%%% 
% % if fitting an elipse to each day, use daily mean variance in vector
% % format:
% xx=find(u_var(1,:)>0); % days where there is at least 2 measurements
% u_var=ones(axxx(xx(i)),1)*u_var(xx(i));
% v_var=ones(axxx(xx(i)),1)*v_var(xx(i)); and
% Amat=[ones(axxx(xx(i)),1) cos(w.*time(bxxx==xx(i))) sin(w.*time(bxxx==xx(i)))];

n=24; %number of data points per day; or use axxx(xx(i)) if not constant
% % for i=1:length(xx);
% % % defining matrix for least squares
% % Amat=[ones(n,1) cos(w.*time(:,xx(i))) sin(w.*time(:,xx(i)))];
% % 
% % Amat=Amat./(u_var(:,xx(i))*ones(1,3));
% % 
% %         Ainv=inv(Amat'*Amat)*Amat'; % SINGULAR!! only diagonal values are diff from zero
% % %         [U,S,V]=svd(Amat'*Amat); % Use SVD to invert the matrix Amat
% % %         InvS=inv(S); 
% % %         InvS(InvS>10)=0; % threshold in S is 0.01, replace w zero where sing value is smaller than 0.01
% % %         InvAmat=(U*InvS*V');
% % %         Ainv=InvAmat*Amat';
% % 
% % u_fit=Ainv*(u_data(:,xx(i))./(u_var(:,xx(i))));
% % u_err=sqrt(diag(inv(Amat'*Amat)));
% % Amat=[ones(n,1) cos(w.*time(:,xx(i))) sin(w.*time(:,xx(i)))]; 
% % Amat=Amat./(v_var(:,xx(i))*ones(1,3));
% % 
% % Ainv=inv(Amat'*Amat)*Amat'; % SINGULAR!! only diagonal values are diff from zero
% % %         [U,S,V]=svd(Amat'*Amat); % Use SVD to invert the matrix Amat
% % %         InvS=inv(S); 
% % %         InvS(InvS>10)=0; % threshold in S is 0.01, replace w zero where sing value is smaller than 0.01
% % %         InvAmat=(U*InvS*V');
% % %         Ainv=InvAmat*Amat';
% % 
% % v_fit=Ainv*(v_data(:,xx(i))./(v_var(:,xx(i))));
% % v_err=sqrt(diag(inv(Amat'*Amat)));
% % 
% % for j=1:100
% % Amat=[ones(n,1) cos(w.*time(:,xx(i))) sin(w.*time(:,xx(i)))]; 
% % Amat=Amat./(u_var(:,xx(i))*ones(1,3));
% % 
% % Ainv=inv(Amat'*Amat)*Amat'; % SINGULAR!! only diagonal values are diff from zero
% % %         [U,S,V]=svd(Amat'*Amat); % Use SVD to invert the matrix Amat
% % %         InvS=inv(S); 
% % %         InvS(InvS>10)=0; % threshold in S is 0.01, replace w zero where sing value is smaller than 0.01
% % %         InvAmat=(U*InvS*V');
% % %         Ainv=InvAmat*Amat';
% % 
% % ud=u_data(:,xx(i))+randn(n,1).*(u_var(:,xx(i)));
% % vd=v_data(:,xx(i))+randn(n,1).*(v_var(:,xx(i)));
% % u_fe(:,j)=Ainv*(ud./(u_var(:,xx(i))));
% % u_ee(:,j)=sqrt(diag(inv(Amat'*Amat)));
% % Amat=[ones(n,1) cos(w.*time(:,xx(i))) sin(w.*time(:,xx(i)))]; 
% % Amat=Amat./(v_var(:,xx(i))*ones(1,3));
% % 
% % Ainv=inv(Amat'*Amat)*Amat'; % SINGULAR!! only diagonal values are diff from zero
% % %         [U,S,V]=svd(Amat'*Amat); % Use SVD to invert the matrix Amat
% % %         InvS=inv(S); 
% % %         InvS(InvS>10)=0; % threshold in S is 0.01, replace w zero where sing value is smaller than 0.01
% % %         InvAmat=(U*InvS*V');    
% % %         Ainv=InvAmat*Amat';
% % 
% % v_fe(:,j)=Ainv*(vd./(v_var(:,xx(i))));
% % v_ee(:,j)=sqrt(diag(inv(Amat'*Amat)));
% % end
% % 
% % %%%determining ellipse value for ++ and --
% %  u_bar=u_fit(1);
% %  v_bar=v_fit(1);
% %  a=u_fit(2);
% %  b=u_fit(3);
% %  c=v_fit(2);
% %  d=v_fit(3);
% %  gamma= b.^2 + d.^2 - a.^2 - c.^2;
% %  beta= -2.*(a.*b+c.*d);
% %  u_bare=u_fe(1,:);
% %  v_bare=v_fe(1,:);
% %  ae=u_fe(2,:);
% %  be=u_fe(3,:);
% %  ce=v_fe(2,:);
% %  de=v_fe(3,:);
% %  gammae= be.^2 + de.^2 - ae.^2 - ce.^2;
% %  betae= -2.*(ae.*be+ce.*de);
% % 
% % % --- phase information (compute t_2 standard deviation relative to mean)
% %  t_1 = atan2(beta, gamma)/(2*w);
% %  t_1 = mod(t_1+12,12);
% %  t_2 = atan2(-beta,-gamma)/(2*w);
% %  t_2 = mod(t_2+12,12);
% % 
% %  t_1e = atan2(betae, gammae)/(2*w);
% %  t_1e = mod(t_1e+12,12);
% %  t_2e = atan2(-betae,-gammae)/(2*w);
% %  t_2e = mod(t_2e+12,12);
% % 
% %  time_err(xx(i))=std(rem(t_2e-t_2+6,12)-6);
% %  time_out(xx(i))=t_2;
% % 
% % 
% % % ---- major and minor axis and errors
% %  minor(xx(i)) = sqrt((a.*cos(w*t_1) + b.*sin(w*t_1)).^2 + (c.*cos(w*t_1) + d.*sin(w*t_1)).^2);
% %  major(xx(i)) = sqrt((a.*cos(w*t_2) + b.*sin(w*t_2)).^2 + (c.*cos(w*t_2) + d.*sin(w*t_2)).^2);
% % 
% %  minore = sqrt((ae.*cos(w*t_1e) + be.*sin(w*t_1e)).^2 + (ce.*cos(w*t_1e) + de.*sin(w*t_1e)).^2);
% %  majore = sqrt((ae.*cos(w*t_2e) + be.*sin(w*t_2e)).^2 + (ce.*cos(w*t_2e) + de.*sin(w*t_2e)).^2);
% %  minor_error2(xx(i))=std(minore);
% %  major_error2(xx(i))=std(majore);
% % 
% % da=u_err(2);
% % db=u_err(3);
% % dc=v_err(2);
% % dd=v_err(3);
% % minor_error_a=(2*a.*(cos(w*t_1)).^2+2*b.*cos(w*t_1).*sin(w*t_1)).*da;
% % minor_error_b=(2*b.*(sin(w*t_1)).^2+2*a.*cos(w*t_1).*sin(w*t_1)).*db;
% % minor_error_c=(2*c.*(cos(w*t_1)).^2+2*d.*cos(w*t_1).*sin(w*t_1)).*dc;
% % minor_error_d=(2*d.*(sin(w*t_1)).^2+2*c.*cos(w*t_1).*sin(w*t_1)).*dd;
% % 
% % major_error_a=(2*a.*(cos(w*t_2)).^2+2*b.*cos(w*t_2).*sin(w*t_2)).*da;
% % major_error_b=(2*b.*(sin(w*t_2)).^2+2*a.*cos(w*t_2).*sin(w*t_2)).*db;
% % major_error_c=(2*c.*(cos(w*t_2)).^2+2*d.*cos(w*t_2).*sin(w*t_2)).*dc;
% % major_error_d=(2*d.*(sin(w*t_2)).^2+2*c.*cos(w*t_2).*sin(w*t_2)).*dd;
% % 
% % minor_error(xx(i))=(sqrt(minor_error_a.^2 + minor_error_b.^2 + minor_error_c.^2 + minor_error_d.^2))./(2*minor(xx(i)));
% % major_error(xx(i))=(sqrt(major_error_a.^2 + major_error_b.^2 + major_error_c.^2 + major_error_d.^2))./(2*major(xx(i)));
% % 
% % % ---- tilt
% %  minor_tilt(xx(i)) = atan2(c.*cos(w*t_1)+d.*sin(w*t_1), a.*cos(w*t_1)+b.*sin(w*t_1));
% %  major_tilt(xx(i)) = atan2(c.*cos(w*t_2)+d.*sin(w*t_2), a.*cos(w*t_2)+b.*sin(w*t_2));
% % 
% %  minor_tilte = atan2(ce.*cos(w*t_1e)+de.*sin(w*t_1e), ae.*cos(w*t_1e)+be.*sin(w*t_1e));
% %  major_tilte = atan2(ce.*cos(w*t_2e)+de.*sin(w*t_2e), ae.*cos(w*t_2e)+be.*sin(w*t_2e));
% %  major_tilt_err(xx(i))=std(rem(major_tilte-major_tilt(xx(i))+3*pi,2*pi)-pi);
% %  minor_tilt_err(xx(i))=std(rem(minor_tilte-minor_tilt(xx(i))+3*pi,2*pi)-pi);
% % 
% % % ----- reversal
% %  reversal_minor(xx(i))= (u_bar*2.*(a.*cos(w*t_1)+b.*sin(w*t_1))+v_bar*2.*(c.*cos(w*t_1)+d.*sin(w*t_1)));
% %  reversal_major(xx(i))= (u_bar*2.*(a.*cos(w*t_2)+b.*sin(w*t_2))+v_bar*2.*(c.*cos(w*t_2)+d.*sin(w*t_2)));
% % 
% % % ----- rotation
% % theta_zero=atan2(c,a);
% % theta_six=atan2(d,b);
% % rindex=rem(theta_six-theta_zero+3*pi,2*pi)-pi;
% % rotation(xx(i))=1;    %%%Counter clockwise
% % if(rindex==0)rotation(xx(i))=0; end
% % if(rindex<0)rotation(xx(i))=-1; end   %%%Clockwise
% % 
% % theta_zero=atan2(ce,ae);
% % theta_six=atan2(de,be);
% % rindex=rem(theta_six-theta_zero+3*pi,2*pi)-pi;
% % rotatione=ones(100,1);  %%%Counter clockwise
% % z=find(rindex==0); if(length(z)>=1) rotatione(z)=0; end
% % z=find(rindex<0); if(length(z)>=1) rotatione(z)=-1; end %%% Clockwise
% % 
% % rotation_e(xx(i))=sum(rotatione)/100;
% % end


