% ws=10 + (45-10)*rand(100,1);
% wd=0 + (359 -0)*rand(100,1);
% wdr=deg2rad(wd);
% vt=-0.5 + (7-(-0.5))*rand(100,1);;
az=0:11.25:359;
azr=deg2rad(az);
inc=deg2rad(30);
for i=1:32
    vr(:,i)=(-ws.*cos(azr(i)-wdr)*sin(deg2rad(30))+vt*cos(deg2rad(30)))+ randn(100,1);

end

ftype=fittype('-a*cos(x-b)*sin(deg2rad(30))+ c*cos(deg2rad(30))');
%ftype=fittype('0.5*a0+a1*cos(x)+b1*sin(x)+a2*cos(2*x)+b2*sin(2*x)');
opts=fitoptions(ftype);
opts.MaxFunEvals=100000;
opts.MaxIter=100000;
opts.Algorithm='Levenberg-Marquardt';
opts.Robust='on';
for i=1:100
    i
    fresult=fit(azr',vr(i,:)',ftype,opts);
    a(i)=fresult.a;
    b(i)=fresult.b;
    c(i)=fresult.c;
end
b2=(b);

u=find(a.*b2 >=0 & a > 0);
n=fix(b2(u)/(2*pi));
b2(u)=b2(u)-n*(2*pi);

u=find(a.*b2 <=0 & a < 0);
n=fix(b2(u)/(pi));
b2(u)=b2(u)-n*(pi);
v=find(n==0);
b2(u(v))=b2(u(v))+pi;

u=find(a.*b2 >=0 & a < 0);
n=fix(b2(u)/(pi));
b2(u)=b2(u)-n*(pi);
v=find(n==0);
b2(u(v))=b2(u(v))+pi;

u=find(a.*b2 <=0 & a > 0);
n=fix(b2(u)/(2*pi));
b2(u)=b2(u)-n*(2*pi);

u=find(b2<0);
b2(u)=b2(u)+ (2*pi);
bd=(rad2deg(b2))';