function[result,ConfidenceInterval] = FitDatasTrd(startyear,startmonth,startday,endyear,endmonth,endday)

warning('off','curvefit:prepareFittingData:removingNaNAndInf');

begindate=datenum(startyear,startmonth,startday);
enddate=datenum(endyear,endmonth,endday);

%Making matrix that will contain results (global)YearlyConfidenceInterval(:,i,j)
%(a0,a1,a2,a3,a4,c1,c2,c3,c4) per amount of days
%the number of data day will be the row.
% a0:mean wind amplitude
% a1:48h tide amplitude
% a2:24h tide amplitude
% a3:12h tide amplitude
% a4:8h  tide amplitude
% c1:48h tide phase
% c2:24h tide phase
% c3:12h tide phase
% c4:8h  tide phase

global result;
result = zeros((enddate-begindate-2),9);
ConfidenceInterval = zeros((enddate-begindate-2),9);

%Curve fitting function
function [fitresult, gof] = Fit(datax, datay, day)
[xData, yData] = prepareCurveData(datax,datay);

% Option setting
ft = fittype('a0+a1*sin((2*pi/48)*(x+c1))+a2*sin((2*pi/24)*(x+c2))+a3*sin((2*pi/12)*(x+c3))+a4*sin((2*pi/8)*(x+c4))','independent','x','dependent','y');
opts = fitoptions('Method','NonlinearLeastSquares');
opts.Display = 'Off';

if day==1 || sum(isnan(result(day-1,:)))==9
   opts.StartPoint = [nanmean(datay) nanstd(datay) nanstd(datay) nanstd(datay) nanstd(datay) 0 0 0 0];
else
   opts.StartPoint = result(day-1,:);
end %if

% fitting the model to the data
[fitresult, gof] = fit(xData, yData, ft, opts);

end %Fit function


for day=1:(enddate-begindate-2) %loop through all files to enddate-3
    
    % Making 4 days segments of V (Zonal wind)
    % time series of zonal wind at 82 km altitude (height bin 1)
    
    StartOfSeg = datevec(begindate+day-1);
    EndOfSeg = datevec(begindate+day+2);
    MRdata=read_hwd_data_trd(StartOfSeg(1),StartOfSeg(2),StartOfSeg(3),EndOfSeg(1),EndOfSeg(2),EndOfSeg(3));
    
    x1 = MRdata.time(1,:,1);
    x2 = MRdata.time(1,:,2)+ 24*ones(size(x1));
    x3 = MRdata.time(1,:,3)+ 48*ones(size(x2));
    x4 = MRdata.time(1,:,4)+ 72*ones(size(x3));
    
    y1 = MRdata.v(1,:,1);
    y2 = MRdata.v(1,:,2);
    y3 = MRdata.v(1,:,3);
    y4 = MRdata.v(1,:,4);
    
    dataX = [x1 x2 x3 x4];
    dataY = [y1 y2 y3 y4];
    
    %Quality CK function checks whether the 4 days segment data are within
    %the quality that is ok to wave fit. Test=1 is good quality and Test=0
    %is not good enough quality. When Test is 1 we will proceed wave
    %fitting and when Test is 0, all the amplitude and phase will be NaN.
    
    if QualityCK(y1,y2,y3,y4) == 1
        DataFit = Fit(dataX, dataY, day);
        result(day,:) = clean(DataFit.a0,DataFit.a1,DataFit.a2,DataFit.a3,DataFit.a4,DataFit.c1,DataFit.c2,DataFit.c3,DataFit.c4);
        %clean data function corrects the negative amplitudes and
        %phases and too big phases
        ci = confint(DataFit);
        ConfidenceInterval(day,:) = [abs(ci(2,1)-ci(1,1)) abs(ci(2,2)-ci(1,2)) abs(ci(2,3)-ci(1,3)) abs(ci(2,4)-ci(1,4)) abs(ci(2,5)-ci(1,5)) abs(ci(2,6)-ci(1,6)) abs(ci(2,7)-ci(1,7)) abs(ci(2,8)-ci(1,8)) abs(ci(2,9)-ci(1,9))];
       
    else
        result(day,:) = [NaN NaN NaN NaN NaN NaN NaN NaN NaN];
        ConfidenceInterval(day,:) = [NaN NaN NaN NaN NaN NaN NaN NaN NaN];
        disp(datevec(begindate+day-1));
    end 
end %for
% 
% plotting the results
% time = (1:enddate-begindate-2)';
%  
%  figure(1);
%  plot(time,result(:,2));
%  xlabel('days');ylabel('amp 48h');
%  
%  figure(2);
%  plot(time,result(:,3));
%  xlabel('days');ylabel('amp 24h');
%  
%  figure(3);
%  plot(time,result(:,4));
%  xlabel('days');ylabel('amp 12h');
%  
%  figure(4);
%  plot(time,result(:,5));
%  xlabel('days');ylabel('amp 8h');
%  
%  figure(5);
%  plot(time,result(:,1));
%  xlabel('days');ylabel('mean');

warning('on','curvefit:prepareFittingData:removingNaNAndInf');

end
