 function FitData3(DataFolder,location,p)
%Fits meridional wind data at location of 4 day periods with the function a0+a1*sin(2*pi*(x+c1))+a2*sin(2*pi*(x*2+c2))+a3*sin(2*pi*(x/2+c3))+a4*sin(2*pi*(x*3+c4))and weights the fitting with the number of meteors for each data point.
%FitData + exclude data if two following 48 hour periods have simultaniousdata gaps bigger than 12 hours
%
%INPUT:
%locatoin = radar
%p=1 => with plots
%p=0 => without plots  
%
%
%OUTPUT:
%Data.Fit.DaysWithFit:      Dates with a fit  [datenum]
%Data.Fit.DaysWithoutFit:   no fit could be produced for this Dates (too few datapoints) [datenum]
%
%Data.Fit.a0:               a0 (Mean) for all fits [size(size(DaysWithFit)]
%Data.Fit.a1:               a1 (Amplitude 24h tide) for all fits [size(size(DaysWithFit)]
%Data.Fit.a2:               a2 (Amplitude 12h tide) for all fits [size(size(DaysWithFit)]
%Data.Fit.a3:               a3 (Amplitude 2 day wave ) for all fits [size(size(DaysWithFit)]
%Data.Fit.c1:               c1 (Phase 24h tide) for all fits [size(size(DaysWithFit)]
%Data.Fit.c2:               c2 (Phase 12h tide) for all fits [size(size(DaysWithFit)]
%Data.Fit.c3:               c3 (Phase  2 day wave) for all fits [size(size(DaysWithFit)]
%
%Data.Fit.a0_95             95% confidence level for a0
%Data.Fit.a1_95             95% confidence level for a1
%Data.Fit.a2_95             95% confidence level for a2
%Data.Fit.a3_95             95% confidence level for a3
%Data.Fit.c1_95             95% confidence level for c1
%Data.Fit.c2_95             95% confidence level for c2
%Data.Fit.c3_95             95% confidence level for c3

%Data.Fit.gof.sse:
%Data.Fit.gof.rsquare:      Rsquare for all fits [size(size(DaysWithFit)]
%Data.Fit.gof.dfe:
%Data.Fit.gof.adjrsquare:
%Data.Fit.gof.rmse:
%Data.Fit.errorflag:        # errorflags for all fits
%Data.Fit.iter:             # Iterations for all fits [size(size(DaysWithFit)]


%%
if p==1
    figure2 = figure('NumberTitle','off','Name',strcat(location,'_FitParameter'),'units','centimeters','Position',[0 0 10 8],...                    % setter str på figuren du vil se på skjermen. VIKIT Gat units kommer FØR Position
        'PaperUnits','centimeters','PaperPosition', [0 0 10 8],...% setter str på figuren du vil lagre (uavhengig av format) VIKITG at units kommer FØR Position
         'DefaultAxesFontSize',12,'DefaultAxesFontWeight','b');               %setter str på font i figure 
    r2=0.8;
else
    r2=1;
end     

%Check teh location exists
if strcmp(location,'GooseBay')||strcmp(location,'Kapuskasing')||strcmp(location,'Pykkvibaer')||strcmp(location,'PrinceGeorge')||strcmp(location,'RankinInlet')||strcmp(location,'Saskatoon')||strcmp(location,'Stokkseyri')||strcmp(location,'Halley')|| strcmp(location,'Hankasalmi')||strcmp(location,'Salmon')||strcmp(location,'Kodiak')||strcmp(location,'Sanae')
    
    load(fullfile(DataFolder,location,strcat(location,'Clean')))
%%
Vx=zeros(4*24,1);
Weight=zeros(size(Vx));
Time=zeros(size(Vx));

start=round(Data.DateNum(1))+2; 
stop=round(Data.DateNum(end))-2;
j=0;
n=0;
Data.Fit.DaysWithFit=zeros(size(stop-start+1,1));
Data.Fit.DaysWithoutFit=zeros(size(stop-start+1,1));
Data.Fit.a0=zeros(size(stop-start+1,1));
Data.Fit.a1=zeros(size(stop-start+1,1));
Data.Fit.a2=zeros(size(stop-start+1,1));
Data.Fit.a3=zeros(size(stop-start+1,1));
Data.Fit.a4=zeros(size(stop-start+1,1));

Data.Fit.c1=zeros(size(stop-start+1,1));
Data.Fit.c2=zeros(size(stop-start+1,1));
Data.Fit.c3=zeros(size(stop-start+1,1));
Data.Fit.c4=zeros(size(stop-start+1,1));

Data.Fit.a0_95=zeros(size(stop-start+1,1));
Data.Fit.a1_95=zeros(size(stop-start+1,1));
Data.Fit.a2_95=zeros(size(stop-start+1,1));
Data.Fit.a3_95=zeros(size(stop-start+1,1));
Data.Fit.a4_95=zeros(size(stop-start+1,1));

Data.Fit.c1_95=zeros(size(stop-start+1,1));
Data.Fit.c2_95=zeros(size(stop-start+1,1));
Data.Fit.c3_95=zeros(size(stop-start+1,1));
Data.Fit.a4_95=zeros(size(stop-start+1,1));

Data.Fit.gof.sse=zeros(size(stop-start+1,1));
Data.Fit.gof.rsquare=zeros(size(stop-start+1,1));
Data.Fit.gof.dfe=zeros(size(stop-start+1,1));
Data.Fit.gof.adjrsquare=zeros(size(stop-start+1,1));
Data.Fit.gof.rmse=zeros(size(stop-start+1,1));
Data.Fit.errorflag=zeros(size(stop-start+1,1));
Data.Fit.iter=zeros(size(stop-start+1,1));
%%

%%
for i=start:stop 
    test=1;
    MaxGabSize=0;
    MaxGabSize2=0;
    dh=zeros(24,1);
    d2d=zeros(48,1);
    Time1=datenum(datevec(i-2:1/24:i+1+23/24))';
    
    [Date,iTime,iDateNum]=intersect(Time1,Data.DateNum);
    if ~isempty(iDateNum)
    s=iDateNum(1);
    e=iDateNum(end);
    
    %datevec(i)
    %datevec(Data.DateNum(s))
    %datevec(Data.DateNum(e))
    
    %Count number of different hours in 4 day data periode
    
    for h=0:23
        check=find(Data.Hour(s:e)==h,1);
        if ~isempty(check)
            dh(h+1)=1;
            MaxGabSize=0;
        else
          MaxGabSize=MaxGabSize+1;
            if MaxGabSize==4% more than 4 houers in the row a missing in dh
                test=0;
                %disp('MaxGabSize>=4')
            end
        end
    end
    Time2=datenum(datevec(i-2:1/24:i-1+23/24))';
    Time3=datenum(datevec(i:1/24:i+1+23/24))';
    [Date2,iTime2]=intersect(Time2,Data.DateNum);
    [Date3,iTime3]=intersect(Time3,Data.DateNum);
    
    d2d(iTime2)=1;
    d2d(iTime3)=1;
    for k=1:48
        if d2d(k)==0
            MaxGabSize2=MaxGabSize2+1;
            if MaxGabSize2==13
                test=0;
                %disp('MaxGabSize2>=12')
            end
        else
            MaxGabSize2=0;
        end
    end
   
 %%   
 
    if  (sum(dh)<16) % Check that the used 4 day data periode consists of at least 16 differen hours (0-23)
        test=0;
        %disp('sum(dh)<16')
    end
    if  (e-s)<48 %Checks that the used 4 day data periode consist of at least 48 different datapoints
       test=0;
       %disp('(e-s)<48')
    end
    
    if test==1
        j=j+1;
       Data.Fit.DaysWithFit(j)=i;
       
       Vx(1:e-s)=Data.Vx(s:e-1);
       
       Weight(1:e-s)=Data.NMeteor(s:e-1);% Weights by the number of Meteors for each dataset.
            [Year,Month,Day,Hour] = AllocateTime(Data,s,e-1); % Creates a date vector of the right size/ dates
       Time(1:e-s)=datenum(Year,Month,Day,Hour,0,0); % transforms the date to datenumbers
       
       
       [cf_,gof,output]=createFitVx3(Time(1:e-s),Vx(1:e-s),Weight(1:e-s),p); %opens the fitting routin
       Data.Fit.errorflag(j)=output.exitflag;
       Data.Fit.iter(j)=output.iterations;
       ci=confint(cf_);
       
       
       Data.Fit.a0(j)=cf_.a0;
       Data.Fit.a0_95(j)=cf_.a0-ci(1,1);
       Data.Fit.a1(j)=cf_.a1;
       Data.Fit.a1_95(j)=cf_.a1-ci(1,2);
       Data.Fit.a2(j)=cf_.a2;
       Data.Fit.a2_95(j)=cf_.a2-ci(1,3);
       Data.Fit.a3(j)=cf_.a3;
       Data.Fit.a3_95(j)=cf_.a3-ci(1,4);
       Data.Fit.a4(j)=cf_.a4;
       Data.Fit.a4_95(j)=cf_.a4-ci(1,5);
       Data.Fit.c1(j)=mod(cf_.c1,1);
       Data.Fit.c1_95(j)=cf_.c1-ci(1,6);
       Data.Fit.c2(j)=mod(cf_.c2,1);
       Data.Fit.c2_95(j)=cf_.c2-ci(1,7);
       Data.Fit.c3(j)=mod(cf_.c3,1);
       Data.Fit.c3_95(j)=cf_.c3-ci(1,8);
       Data.Fit.c4(j)=mod(cf_.c4,1);
       Data.Fit.c4_95(j)=cf_.c4-ci(1,9);
       
       
       Data.Fit.gof.sse(j)=gof.sse;
       Data.Fit.gof.rsquare(j)=gof.rsquare;
       Data.Fit.gof.dfe(j)=gof.dfe;
       Data.Fit.gof.adjrsquare(j)=gof.adjrsquare;
       Data.Fit.gof.rmse(j)=gof.rmse;
       %%
       
       if gof.rsquare>r2
       hold on
       errorbar(Time(1:e-s),Vx(1:e-s),1000./Data.NMeteor(1:e-s))
       datetick('x','yymmdd','keeplimits');
       ylabel('meridional wind [m/s]')
       legend('data','fit')
       title(datestr(Time(1),'dd-mmm-yy'))
       set(gca,'xticklabel',[1:4])
       hold off
       pause(1)
       
      %figure2;
       % x=(Time(1):1/24:Time(e-s));
       %plot(x,Data.Fit.a1(j)*sin(2*pi*x-Data.Fit.c1(j)))
       %hold on
      % plot(x,Data.Fit.a2(j)*sin(2*pi*2*x-Data.Fit.c2(j)),'r')
      % plot(x,Data.Fit.a3(j)*sin(2*pi*3*x-Data.Fit.c3(j)),'g')
      % datetick('x','yymmdd','keepticks');
      % hold off
       else
           close
       end
    else
        n=n+1;
        Data.Fit.DaysWithoutFit(n)=i;
    end
       %%
    if p==1
    if j==5
        pause
    end
    end
end
end
%%
Data.Fit.DaysWithFit=Data.Fit.DaysWithFit(Data.Fit.DaysWithFit~=0);
Data.Fit.DaysWithoutFit=Data.Fit.DaysWithoutFit(1:n);

Data.Fit.a0=Data.Fit.a0(1:j);
Data.Fit.a1=Data.Fit.a1(1:j);
Data.Fit.a2=Data.Fit.a2(1:j);
Data.Fit.a3=Data.Fit.a3(1:j);
Data.Fit.a4=Data.Fit.a4(1:j);
Data.Fit.c1=Data.Fit.c1(1:j);
Data.Fit.c2=Data.Fit.c2(1:j);
Data.Fit.c3=Data.Fit.c3(1:j);
Data.Fit.c4=Data.Fit.c4(1:j);
Data.Fit.gof.sse=Data.Fit.gof.sse(1:j);
Data.Fit.gof.rsquare=Data.Fit.gof.rsquare(1:j);
Data.Fit.gof.dfe=Data.Fit.gof.dfe(1:j);
Data.Fit.gof.adjrsquare=Data.Fit.gof.adjrsquare(1:j);
Data.Fit.gof.rmse=Data.Fit.gof.rmse(1:j);
Data.Fit.errorflag=Data.Fit.errorflag(1:j);
Data.Fit.iter=Data.Fit.iter(1:j);

fprintf(['# Days with fit: ',int2str(j),'\n'])
fprintf(['# Days without fit: ',int2str(n),'\n'])

%%
save(fullfile(DataFolder,location,strcat(location,'_Fit3')),'Data')
filename=[location,'_Fit3.fig'];

%%
%Plots fitting parameters a0,a1,a2,a3,c1,c2,c3, #of iterations, errorflag
%and r^2 for all fits and the data coverage of teh fitting
PlotFitParameter2(DataFolder,location,filename,Data.Fit.DaysWithFit,...
                 Data.Fit.errorflag,...
                 Data.Fit.iter,...
                 Data.Fit.gof.rsquare,...
                 Data.Fit.a0,...
                 Data.Fit.a1,...
                 Data.Fit.c1,...
                 Data.Fit.a2,...
                 Data.Fit.c2,...
                 Data.Fit.a3,...
                 Data.Fit.c3,...
                 Data.Fit.a4,...
                 Data.Fit.c4,...
                 Data.Fit.DaysWithFit,j,...
                 Data.Fit.DaysWithoutFit,n)

 else
     fprintf([location,' does not exsists \n']);
 end
 end