function [] = LonFit5_long(DataFolder,RadarsWorking,p,l)
%Creates a longitudal fit for anomlay wind (weighted smoothed climatology removed):
%S0+S1*sin(2*Pi*lon/360+c1)+S2*sin(2*Pi*lon/180+c2)over RadarsWorking radars in 2 steps
location0='TimeLonVec';
location='LonFit';
k=0;
%% Max number of Radars in fit
if RadarsWorking>10
    disp('Max number off radars =10')
    disp('Number of radars included in fit is set to 10')
    RadarsWorking=10;
else
    disp(['Number of radars included in fit is ', int2str(RadarsWorking)]) 
end

%% Load TimeLonVecWind and initalies Fit structure

load(fullfile(DataFolder,location0,'TimeLonVecWind_long3')) %includ 8hour and 2 day


TimeLonVec.Fit.S0=TimeLonVec.Time*NaN;
TimeLonVec.Fit.S1=TimeLonVec.Fit.S0;
TimeLonVec.Fit.S2=TimeLonVec.Fit.S0;
TimeLonVec.Fit.c1=TimeLonVec.Fit.S0;
TimeLonVec.Fit.c2=TimeLonVec.Fit.S0;

TimeLonVec.Fit.S0_95=TimeLonVec.Fit.S0;
TimeLonVec.Fit.S1_95=TimeLonVec.Fit.S0;
TimeLonVec.Fit.S2_95=TimeLonVec.Fit.S0;
TimeLonVec.Fit.c1_95=TimeLonVec.Fit.S0;
TimeLonVec.Fit.c2_95=TimeLonVec.Fit.S0;

TimeLonVec.Fit.gof.sse=TimeLonVec.Fit.S0;
TimeLonVec.Fit.gof.rsquare=TimeLonVec.Fit.S0;
TimeLonVec.Fit.gof.dfe=TimeLonVec.Fit.S0;
TimeLonVec.Fit.gof.adjrsquare=TimeLonVec.Fit.S0;
TimeLonVec.Fit.gof.rmse=TimeLonVec.Fit.S0;
TimeLonVec.Fit.errorflag=TimeLonVec.Fit.S0;
TimeLonVec.Fit.iter=TimeLonVec.Fit.S0;

TimeLonVec.Fit.MeanTest=TimeLonVec.Fit.S0;
yyy=TimeLonVec.Fit.S0;
yyyp=TimeLonVec.Fit.S0;
%% Weight Fit with % error of a0. TODO!
%TimeLonVec.Weight=1./(TimeLonVec.MeanWind_95P).^2;
TimeLonVec.Weight=1./TimeLonVec.MeanWindWeighted_err.^2;
TimeLonVec.MeanWind=TimeLonVec.MeanWindWeighted;
%% Find times where choosen number of radars are operating
INT=zeros(size(TimeLonVec.Time,2),1);
for i=1:size(TimeLonVec.Time,2)
    LonS=TimeLonVec.Lon(~isnan(TimeLonVec.MeanWind(:,i)));
    %LonS=LonS(abs(LonS)>0)
  
    if (sum(isnan(TimeLonVec.MeanWind(:,i)))==10-RadarsWorking) && (isnan(TimeLonVec.MeanWind(9,i))+isnan(TimeLonVec.MeanWind(10,i))<=1) && (isnan(TimeLonVec.MeanWind(1,i))==0) && (max(abs(LonS(2:length(LonS))-LonS(1:length(LonS)-1)))<80)
    INT(i)=1;
    if RadarsWorking==6 
        iii=0;
        for ii=1:10
             if isnan(TimeLonVec.MeanWind(ii,i))==1
                 iii=iii+1;
             else iii=0;
             end
             if iii==2 %2 radars next to each other do not function=> no fit
                datevec(TimeLonVec.Time(i))
                INT(i)=0;
             end
        end
    end
    end
end
ind=find(INT==1);
indn=find(INT==0);
disp(['Number of fits: ',int2str(size(ind,1))])



for i=1:size(ind,1)
    
ind2=find(isnan(TimeLonVec.MeanWind(:,ind(i)))==0);

    %% Fit Data
    [cf_,gof,output,MeanTest]=createFitLonS1S2_0(TimeLonVec.Lon(ind2),TimeLonVec.MeanWind(ind2,ind(i)),TimeLonVec.Weight(ind2,ind(i)),0);
    %[cf_,gof,output,MeanTest]=createFitLonS2(TimeLonVec.Lon(ind2),TimeLonVec.MeanWind(ind2,ind(i)),TimeLonVec.Weight(ind2,ind(i)),p);
    %% Sort fit into structure
    
  
    TimeLonVec.Fit.MeanTest(ind(i))=MeanTest;
    TimeLonVec.Fit.errorflag(ind(i))=output.exitflag;
    TimeLonVec.Fit.iter(ind(i))=output.iterations;
    TimeLonVec.Fit.gof.rsquare(ind(i))=gof.rsquare;
    ci=confint(cf_);
    
    TimeLonVec.Fit.S2(ind(i))=cf_.S2;
    TimeLonVec.Fit.S2_95(ind(i))=cf_.S2-ci(1,1);
    TimeLonVec.Fit.c2(ind(i))=mod(cf_.c2,1);
    TimeLonVec.Fit.c2_95(ind(i))=cf_.c2-ci(1,2);

    yyy(ind2,ind(i))=TimeLonVec.MeanWind(ind2,ind(i))-cf_.S2*sin(2*pi*(TimeLonVec.Lon(ind2)/180+cf_.c2));
    %[cf_,gof,output]=createFitLonS1_0(TimeLonVec.Lon(ind2),yyy(ind2,ind(i)),TimeLonVec.Weight(ind2,ind(i)),0);
    %% Sort fit into structure
    
    
    
    
     TimeLonVec.Fit.S1(ind(i))=cf_.S1;
    TimeLonVec.Fit.S1_95(ind(i))=cf_.S1-ci(1,3);
     TimeLonVec.Fit.c1(ind(i))=mod(cf_.c1,1);
    TimeLonVec.Fit.c1_95(ind(i))=cf_.c1-ci(1,4);
    
    TimeLonVec.Fit.S0(ind(i))=mean(TimeLonVec.MeanWind(ind2,ind(i))-...
            TimeLonVec.Fit.S2(ind(i))*sin(2*pi*(TimeLonVec.Lon(ind2)/180+TimeLonVec.Fit.c2(ind(i))))-...
            TimeLonVec.Fit.S1(ind(i))*sin(2*pi*(TimeLonVec.Lon(ind2)/360+TimeLonVec.Fit.c1(ind(i)))));
    
   TimeLonVec.Fit.S0_95(ind(i))=std(TimeLonVec.MeanWind(ind2,ind(i))-...
            TimeLonVec.Fit.S2(ind(i))*sin(2*pi*(TimeLonVec.Lon(ind2)/180+TimeLonVec.Fit.c2(ind(i))))-...
            TimeLonVec.Fit.S1(ind(i))*sin(2*pi*(TimeLonVec.Lon(ind2)/360+TimeLonVec.Fit.c1(ind(i))))); 
   
        
    if p==1 
    get_figure([20,20],14,'n')
    lon=[-180:1:180];
    errorbar(TimeLonVec.Lon(ind2),TimeLonVec.MeanWind(ind2,ind(i)),TimeLonVec.MeanWind_95P(ind2,ind(i))) %data
    hold on 
    Fitplot=TimeLonVec.Fit.S0(ind(i))+...
            TimeLonVec.Fit.S1(ind(i))*sin(2*pi*(lon/360+TimeLonVec.Fit.c1(ind(i))))+...
            TimeLonVec.Fit.S2(ind(i))*sin(2*pi*(lon/180+TimeLonVec.Fit.c2(ind(i))));
        plot(lon,Fitplot,'r') %fit
       hold on
    
    
    % plot(TimeLonVec.Lon,yyy(:,ind(i)),'k')%data-S2
    yyyp2=cf_.S1*sin(2*pi*(lon/360+cf_.c1)); 
    plot(lon,yyyp2,'g-.','linewidth',1) %S1
    yyyp=TimeLonVec.Fit.S2(ind(i))*sin(2*pi*(lon/180+TimeLonVec.Fit.c2(ind(i))));
    plot(lon,yyyp,'m-.','linewidth',1)%S2
    
    plot([-180,180],[TimeLonVec.Fit.S0(ind(i)),TimeLonVec.Fit.S0(ind(i))],'k-.')
   
    legend('data','fit','S1+mean','S2+mean','mean')
    xlim([-180,180])
    A=datevec(TimeLonVec.Time(ind(i)));
    title(int2str(A(1:3)))
   
    if l==2
    saveas(gca,fullfile(DataFolder,location,[int2str(A(1:3)),'Fit5_S2S1']),'fig')
    elseif l==3
        saveas(gca,fullfile(DataFolder,location,[int2str(A(1:3)),'Fit5_S2S1_3']),'fig') %uses TimeLonVec were also 8 hour tide is removed
    end
    
     if l==2
    makePDF(gca,[20,20],[int2str(A(1:3)),'Fit5_S2S1'],fullfile(DataFolder,location))
    elseif l==3
    makePDF(gca,[20,20],[int2str(A(1:3)),'Fit5_S2S1_3'],fullfile(DataFolder,location))  
    end
    end
    
    k=k+1;
    if k==10
        p=0;
    end
end

if l==2
save(fullfile(DataFolder,location,['TimeLonVecWind_Fit5',int2str(RadarsWorking),'_long']),'TimeLonVec')
filename=[location,'_Fit5',int2str(RadarsWorking),'_long','.fig']
elseif l==3
    save(fullfile(DataFolder,location,['TimeLonVecWind_Fit5',int2str(RadarsWorking),'_long3']),'TimeLonVec')
    filename=[location,'_Fit5',int2str(RadarsWorking),'_long3','.fig']
end


%% Plots fitting parameters S0,S1,S2,c1,c2, #of iterations, errorflag
%and r^2 for all fits and the data coverage of teh fitting
PlotFitParameterLon(DataFolder,location,filename,TimeLonVec.Time,...
                 TimeLonVec.Fit.errorflag,...
                 TimeLonVec.Fit.iter,...
                 TimeLonVec.Fit.gof.rsquare,...
                 TimeLonVec.Fit.S0,...
                 TimeLonVec.Fit.S1,...
                 TimeLonVec.Fit.c1,...
                 TimeLonVec.Fit.S2,...
                 TimeLonVec.Fit.c2,...
                 TimeLonVec.Time(ind),size(ind,1),...
                 TimeLonVec.Time(indn),size(indn,1))
close all
end

