function [mode, offset]=skew_fit_no_plot3(data); 
%fit a series of observations to a skew normal distribution
%mode = mode of the distribution
%offset = the difference between the mode and the mean; due to the skewness

%calculate the first three moments of the data
meaN=mean(data);
va=var(data);
ske=mean((data-meaN).^3)/(va)^(3/2) ;


if ske > 1 
   warning(['Skewness has value ', num2str(ske),' and will be adjusted to 1 to calculate mode and offset']); 
    ske=.9;
end

if ske <= -1
    warning(['Skewness has value ', num2str(ske),' and will be adjusted to -1 to calculate mode and offset']); 
    ske=-.9;
end
    

%reverse engineer the parameters of the distribution
gamma=(ske/abs(ske))*((pi/2)^.5)*((abs(ske))^(1/3)/((abs(ske)^(2/3)+((4-pi)/2)^(2/3))^(1/2)));

%alpha2=tan(asin(gamma));
alpha=gamma/((1-gamma^2)^.5);
gamma2=alpha/(1+alpha^2)^.5;
sig=(va/(1-2*(gamma^2)/pi))^.5;
mea=meaN-sig*gamma*(2/pi)^.5;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NP=500; %number of points;
x=min(data)-.3*(max(data)-min(data)):(max(data)-min(data))/NP:max(data)+.3*(max(data)-min(data));

G=(1/((2*pi*sig^2)^.5))*exp((-(x-mea).^2)/(2*sig^2));
F=.5*(1+erf((alpha/sig)*((x-mea)/2^.5)));
Gs=2*G.*F;
% NH=20; %number of divisions in histogram
% [num,xout]=hist(data,NH);
% figure(11); clf;
% plot(xout,num, 'b'); hold on;
% plot(x,Gs*trapz(xout,num), 'r')
% 
% 
% %calculate the first three orders
% nMEA=mean(x.*Gs)/mean(Gs);



%now search for the mode, numerically.

xmode1=x(find(Gs==max(Gs)));

%now make the search zoom in;
interva=(max(data)-min(data))/NP;
Nzoom=100;
x2=xmode1-interva:interva/Nzoom:xmode1+interva;
Gzoom=(1/((2*pi*sig^2)^.5))*exp((-(x2-mea).^2)/(2*sig^2));
Fzoom=.5*(1+erf((alpha/sig)*((x2-mea)/2^.5)));
Gszoom=2*Gzoom.*Fzoom;


mode=x2(find(Gszoom == max(Gszoom)));
% figure(69);clf; plot(x2, Gszoom,'b', mode, max(Gszoom), 'kp')
% title('zoomed in max')

% figure(11); hold on;
% plot([mode mode],[0 (1.02*max(Gs))*trapz(xout,num)], 'b--', [meaN meaN], [0 1.02*(max(Gs))*trapz(xout,num)], 'k--');
% legend('data', 'skew-normal', 'MODE', 'MEAN')

offset=meaN-mode;
%gamma=alpha/((1+alpha^2)^.5);
%aMEA=mea+sig*gamma*((2/pi)^.5)
%nVAR=mean(Gs.*(x-nMEA).^2)/mean(Gs)
%aVAR=(sig^2)*(1-(2*gamma^2)/pi)

%%skewness
%aSKE=((4-pi)/2)*((gamma*(2/pi)^.5)^3)/(1-2*(gamma^2)/pi)^(3/2)
%nSKE=(mean(Gs.*(x-nMEA).^3)/mean(Gs))./(nVAR^(3/2))



