function [Data] = ProcessedMonthlyWeightedMeanLST(DataFolder, stations, location,f,p)
%Calculates Monthly Mean of FitParameters weigthed with error for all years between yearStart and yearEnd
% For more info on FitParameters see FitData. For more info on structure Data see readData   
%
% Data.Processed.Monthly.Mean
%.a0 / .a1 / .a2 / .a3 / .c1 / .c2 /.c3 
% Monthly Weighted Mean for each FitPrameter

% Data.Processed.Monthly.STD
%.a0 / .a1 / .a2 / .a3 / .c1 / .c2 /.c3 
%error of weighttes mean
%
% Saves as Data.Name=[location,'_Weighted_Monthly_',int2str(yearStart),'-',int2str(yearEnd)];
% at C:\Users\kleinkne\Documents\PHD\DATA
%
% Calls PlotMonthly for plotting off all monthly means
if f==1 || f==2 || f==3

if strmatch(location,stations)>0%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')
   
   
    if f==1
    load(fullfile(DataFolder,location,[location,'_Fit']));
    elseif f==2
    load(fullfile(DataFolder,location,[location,'_Fit2']));
    elseif f==3
    load(fullfile(DataFolder,location,[location,'_Fit3']));
    end
    
DV=datevec(Data.Fit.DaysWithFit);
yearStart=DV(1,1);
yearEnd=DV(end,1)

%if yearStart<DV(1,1)|| yearEnd>DV(end,1)
%    fprintf(['Datarange ',location,' is ', int2str(DV(1,1)),' - ', int2str(DV(end,1)),'\n'])
   
%else    

    if f==1
        
            Data.Name=[location,'_Weighted_Monthly_LST',int2str(yearStart),'-', int2str(yearEnd)];
   
    elseif f==2
       
            Data.Name=[location,'_Weighted_Monthly_LST',int2str(yearStart),'-', int2str(yearEnd),'_2'];
    else
            Data.Name=[location,'_Weighted_Monthly_LST',int2str(yearStart),'-', int2str(yearEnd),'_3'];
    end
Data.Processed.Monthly.Mean.a0=zeros(1,12);
Data.Processed.Monthly.Mean.a1=zeros(1,12);


Data.Processed.Monthly.Mean.a2=zeros(1,12);
Data.Processed.Monthly.Mean.a3=zeros(1,12);
Data.Processed.Monthly.Mean.c1=zeros(1,12);
Data.Processed.Monthly.Mean.c1_test=zeros(1,12);
Data.Processed.Monthly.Mean.c2=zeros(1,12);
Data.Processed.Monthly.Mean.c3=zeros(1,12);
if f==3
Data.Processed.Monthly.Mean.c4=zeros(1,12);
end
Data.Processed.Monthly.STD.a0=zeros(1,12);
Data.Processed.Monthly.STD.a1=zeros(1,12);
Data.Processed.Monthly.STD.a2=zeros(1,12);
Data.Processed.Monthly.STD.a3=zeros(1,12);
Data.Processed.Monthly.STD.c1=zeros(1,12);
Data.Processed.Monthly.STD.c2=zeros(1,12);
Data.Processed.Monthly.STD.c3=zeros(1,12);
if f==3
   Data.Processed.Monthly.STD.c4=zeros(1,12);
end

% figure;
% plot(Data.Fit.DaysWithFit,Data.Fit.a0_95,'.')
% title([location,' a0 95'])
% datetick('x','ddmmyy','keepticks');
% figure;
% plot(Data.Fit.DaysWithFit,Data.Fit.a1_95,'.')
% title([location,' a1 95'])
% datetick('x','ddmmyy','keepticks');
% figure;
% plot(Data.Fit.DaysWithFit,Data.Fit.a2_95,'.')
% title([location,' a2 95'])
% datetick('x','ddmmyy','keepticks');
% figure;
% plot(Data.Fit.DaysWithFit,Data.Fit.a3_95,'.')
% title([location,' a3 95'])
% datetick('x','ddmmyy','keepticks');
% figure;
% plot(Data.Fit.DaysWithFit,Data.Fit.c1_95,'.')
% title([location,' c2 95'])
% datetick('x','ddmmyy','keepticks');
% figure;
% plot(Data.Fit.DaysWithFit,Data.Fit.a0_95,'.')
% title([location,' c3 95'])
% datetick('x','ddmmyy','keepticks');


indStart=find(DV(:,1)==yearStart);
indStart=indStart(1);
indEnd=find(DV(:,1)==yearEnd);
if strmatch(location,stations{6})>0
indEnd=find(Data.Fit.DaysWithFit<=datenum(2009,3,3));
end
indEnd=indEnd(end);

e=1;
% for i=yearStart:yearEnd 
% s=e;
% e=e+size(find(DV(:,1)==i),1); 
% DVJ(s:e-1)=find(DV(:,1)==i);
% end
for i=1:12
%M=find(DV(DVJ(1:e-1),2)==i);
M=find(DV(indStart:indEnd,2)==i)+indStart-1;
%%
Data.Processed.Monthly.Mean.a0(i)=sum(Data.Fit.a0(M)./(Data.Fit.a0_95(M)./2).^2)/sum(1./(Data.Fit.a0_95(M)./2).^2);
Data.Processed.Monthly.STD.a0(i)=(1/sum(1./(Data.Fit.a0_95(M)./2).^2))^(1/2);
% Data.Processed.Monthly.Mean.a1(i)=sum(Data.Fit.a1(M)./(Data.Fit.a1_95(M)./2).^2)/sum(1./(Data.Fit.a1_95(M)./2).^2);
% Data.Processed.Monthly.STD.a1(i)=(1/sum(1./(Data.Fit.a1_95(M)./2).^2))^(1/2);
% Data.Processed.Monthly.Mean.c1(i)=sum(Data.Fit.c1(M)./(Data.Fit.c1_95(M)./2).^2)/sum(1./(Data.Fit.c1_95(M)./2).^2);
% Data.Processed.Monthly.STD.c1(i)=(1/sum(1./(Data.Fit.c1_95(M)./2).^2))^(1/2);
% Data.Processed.Monthly.Mean.a2(i)=sum(Data.Fit.a2(M)./(Data.Fit.a2_95(M)./2).^2)/sum(1./(Data.Fit.a2_95(M)./2).^2);
% Data.Processed.Monthly.STD.a2(i)=(1/sum(1./(Data.Fit.a2_95(M)./2).^2))^(1/2);
% Data.Processed.Monthly.Mean.c2(i)=sum(Data.Fit.c2(M)./(Data.Fit.c2_95(M)./2).^2)/sum(1./(Data.Fit.c2_95(M)./2).^2);
%Data.Processed.Monthly.STD.c2(i)=(1/sum(1./(Data.Fit.c2_95(M)./2).^2))^(1/2);
Data.Processed.Monthly.Mean.a3(i)=sum(Data.Fit.a3(M)./(Data.Fit.a3_95(M)./2).^2)/sum(1./(Data.Fit.a3_95(M)./2).^2);
%Data.Processed.Monthly.STD.a3(i)=(1/sum(1./(Data.Fit.a3_95(M)./2).^2))^(1/2);

a1=Data.Fit.a1(M); %fitted amplitudes for all days in the month
c1=Data.Fit.c1(M); %fitted phase for all days in the month
da1=Data.Fit.a1_95(M)./2; %error of fitted amplitudes for all days in the month
dc1=Data.Fit.c1_95(M)./2;  %error of fitted phase for all days in the month
A1=Data.Fit.a1(M).*cos(2*pi()*Data.Fit.c1(M)); %transfer Amplitude an phase to a and b
B1=Data.Fit.a1(M).*sin(2*pi()*Data.Fit.c1(M)); %transfer Amplitude an phase to a and b
dA1=((da1.*cos(2*pi.*c1)).^2+(a1.*sin(2*pi.*c1)*2*pi().*dc1).^2); %calc error of a
dB1=((da1.*sin(2*pi.*c1)).^2+(a1.*cos(2*pi.*c1)*2*pi().*dc1).^2); %calc error of b
MA1=sum(A1./dA1)/sum(1./dA1); %weighted monthly mean of a
dMA1=1/sum(1./dA1); %weighted error of the mean for a
MB1=sum(B1./dB1)/sum(1./dB1); %weighted monthly mean of b
dMB1=1/sum(1./dB1);  %weighted error of the mean for b

Data.Processed.Monthly.Mean.a1(i)=(MA1^2+MB1^2)^(1/2); % transfer a and b back to amplitude
Data.Processed.Monthly.STD.a1(i)=((MA1^2*dMA1^2+MB1^2+dMB1^2)/(MA1^2+MB1^2))^(1/2); % transfer error of a and b back to error of amplitude 
Data.Processed.Monthly.Mean.c1(i)=mod(((-atan2(MB1,MA1)+pi()/2)/(2*pi())),1)*24+Data.Lon(1)/15;  % transfer a and b back to phase % converts to maximum position in UT
Data.Processed.Monthly.STD.c1(i)=((MA1^2*dMA1^2+MB1^2+dMB1^2)/(MA1^2+MB1^2)^2)^(1/2)*24/(2*pi);  % transfererror of  a and b back to error of phase %converts from -pi->pi into 0 to 1!!                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     

a2=Data.Fit.a2(M);
c2=Data.Fit.c2(M);
da2=Data.Fit.a2_95(M)./2;
dc2=Data.Fit.c2_95(M)./2;
A2=Data.Fit.a2(M).*cos(2*pi()*Data.Fit.c2(M));
B2=Data.Fit.a2(M).*sin(2*pi()*Data.Fit.c2(M));
dA2=((da2.*cos(2*pi.*c2)).^2+(a2.*sin(2*pi.*c2)*2*pi().*dc2).^2);
dB2=((da2.*sin(2*pi.*c2)).^2+(a2.*cos(2*pi.*c2)*2*pi().*dc2).^2);
MA2=sum(A2./dA2)/sum(1./dA2);
dMA2=1/sum(1./dA2);
MB2=sum(B2./dB2)/sum(1./dB2);
dMB2=1/sum(1./dB2);

Data.Processed.Monthly.Mean.a2(i)=(MA2^2+MB2^2)^(1/2);
Data.Processed.Monthly.STD.a2(i)=((MA2^2*dMA2^2+MB2^2+dMB2^2)/(MA2^2+MB2^2))^(1/2);
Data.Processed.Monthly.Mean.c2(i)=mod(((-atan2(MB2,MA2)+pi()/2)/(2*pi())),1)*12+Data.Lon(1)/15; % converts to maximum position in UT
Data.Processed.Monthly.STD.c2(i)=((MA2^2*dMA2^2+MB2^2+dMB2^2)/(MA2^2+MB2^2)^2)^(1/2)*12/(2*pi); %convert from -pi->pi into 0 to 1!!  


%ind=find(mod(Data.Fit.DaysWithFit,2)==0);
%Data.Fit.c3(ind)=Data.Fit.c3(ind)+0.5;
a3=Data.Fit.a3(M);
c3=Data.Fit.c3(M);
da3=Data.Fit.a3_95(M)./2;
dc3=Data.Fit.c3_95(M)./2;
A3=Data.Fit.a3(M).*cos(2*pi()*Data.Fit.c3(M));
B3=Data.Fit.a3(M).*sin(2*pi()*Data.Fit.c3(M));
dA3=((da3.*cos(2*pi.*c3)).^2+(a3.*sin(2*pi.*c3)*2*pi().*dc3).^2);
dB3=((da3.*sin(2*pi.*c3)).^2+(a3.*cos(2*pi.*c3)*2*pi().*dc3).^2);
MA3=sum(A3./dA3)/sum(1./dA3);
dMA3=1/sum(1./dA3);
MB3=sum(B3./dB3)/sum(1./dB3);
dMB3=1/sum(1./dB3);
%Data.Processed.Monthly.Mean.a3(i)=(MA3^2+MB3^2)^(1/2);
Data.Processed.Monthly.STD.a3(i)=((MA3^2*dMA3^2+MB3^2+dMB3^2)/(MA3^2+MB3^2))^(1/2);
Data.Processed.Monthly.Mean.c3(i)=mod(((-atan2(MB3,MA3)+pi()/2)/(2*pi())),1)*48+Data.Lon(1)/15; % converts to maximum position in UT   c3 always from 1.day of 4 day intervall!! must be changed to get a right calculation
Data.Processed.Monthly.STD.c3(i)=((MA3^2*dMA2^2+MB3^2+dMB2^2)/(MA3^2+MB3^2)^2)^(1/2)*48/(2*pi); %convert from -pi->pi into 0 to 1!!  


if f==3
a4=Data.Fit.a4(M);
c4=Data.Fit.c4(M);
da4=Data.Fit.a4_95(M)./2;
dc4=Data.Fit.c4_95(M)./2;
A4=Data.Fit.a4(M).*cos(2*pi()*Data.Fit.c4(M));
B4=Data.Fit.a4(M).*sin(2*pi()*Data.Fit.c4(M));
dA4=((da4.*cos(2*pi.*c4)).^2+(a4.*sin(2*pi.*c4)*2*pi().*dc4).^2);
dB4=((da4.*sin(2*pi.*c4)).^2+(a4.*cos(2*pi.*c4)*2*pi().*dc4).^2);
MA4=sum(A4./dA4)/sum(1./dA4);
dMA4=1/sum(1./dA4);
MB4=sum(B4./dB4)/sum(1./dB4);
dMB4=1/sum(1./dB4);

Data.Processed.Monthly.Mean.a4(i)=(MA4^2+MB4^2)^(1/2);
Data.Processed.Monthly.STD.a4(i)=((MA4^2*dMA4^2+MB4^2+dMB4^2)/(MA4^2+MB4^2))^(1/2);
Data.Processed.Monthly.Mean.c4(i)=mod(((-atan2(MB4,MA4)+pi()/2)/(2*pi())),1)*8+Data.Lon(1)/15; % converts to maximum position in UT
Data.Processed.Monthly.STD.c4(i)=((MA4^2*dMA4^2+MB4^2+dMB4^2)/(MA4^2+MB4^2)^2)^(1/2)*8/(2*pi); %convert from -pi->pi into 0 to 1!!  
end


%%
end
%figure
%C=[0 0 0;0 0.5 0; 0 1 0; 1 0 0; 0 0 1; 1 0 1; 0 1 1;1 1 0; 1 0.5 0;0.5 0 0.5;0 0.5 1; 0.5 1 0];
%x=0:1/24:1;
%for i=1:12
%W=sin(2*pi()*(x+Data.Processed.Monthly.Mean.c1(i)));
%plot(0:24,W,'Color',C(i,:))
%hold on
%end
%legend('Jan','Feb','Mar','Apr','Mai','Jun','Jul','Aug','Sep','Oct','Nov','Dec')
%title([location,'Phase UT 24 h'])

%figure
%C=[0 0 0;0 0.5 0; 0 1 0; 1 0 0; 0 0 1; 1 0 1; 0 1 1;1 1 0; 1 0.5 0;0.5 0 0.5;0 0.5 1; 0.5 1 0];
%x=0:1/24:1;
%for i=1:12
%W=sin(2*pi()*(x+Data.Processed.Monthly.Mean.c2(i)));
%plot(0:24,W,'Color',C(i,:))
%hold on
%end
%legend('Jan','Feb','Mar','Apr','Mai','Jun','Jul','Aug','Sep','Oct','Nov','Dec')
%title([location,'Phase UT 12 h'])

save(fullfile(DataFolder,location,Data.Name),'Data')



if p>0
    if f<3
        PlotMonthlyLST(DataFolder,Data.Name,location,1:12,...
                        Data.Processed.Monthly.Mean.a0,Data.Processed.Monthly.STD.a0,... 
                        Data.Processed.Monthly.Mean.a1,Data.Processed.Monthly.STD.a1,...
                        Data.Processed.Monthly.Mean.c1,Data.Processed.Monthly.STD.c1,...
                        Data.Processed.Monthly.Mean.a2,Data.Processed.Monthly.STD.a2,...
                        Data.Processed.Monthly.Mean.c2,Data.Processed.Monthly.STD.c2,...
                        Data.Processed.Monthly.Mean.a3,Data.Processed.Monthly.STD.a3,...
                        Data.Processed.Monthly.Mean.c3,Data.Processed.Monthly.STD.c3)
    else
        PlotMonthlyLST2(DataFolder,Data.Name,location,1:12,...
                        Data.Processed.Monthly.Mean.a0,Data.Processed.Monthly.STD.a0,... 
                        Data.Processed.Monthly.Mean.a1,Data.Processed.Monthly.STD.a1,...
                        Data.Processed.Monthly.Mean.c1,Data.Processed.Monthly.STD.c1,...
                        Data.Processed.Monthly.Mean.a2,Data.Processed.Monthly.STD.a2,...
                        Data.Processed.Monthly.Mean.c2,Data.Processed.Monthly.STD.c2,...
                        Data.Processed.Monthly.Mean.a3,Data.Processed.Monthly.STD.a3,...
                        Data.Processed.Monthly.Mean.c3,Data.Processed.Monthly.STD.c3,...
                        Data.Processed.Monthly.Mean.a4,Data.Processed.Monthly.STD.a4,...
                        Data.Processed.Monthly.Mean.c4,Data.Processed.Monthly.STD.c4)
    end
end
if p>1
figure
errorbar(1:12,Data.Processed.Monthly.Mean.a0,Data.Processed.Monthly.STD.a0)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.a1,Data.Processed.Monthly.STD.a1)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.a2,Data.Processed.Monthly.STD.a2)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.a3,Data.Processed.Monthly.STD.a3)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.a4,Data.Processed.Monthly.STD.a4)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.c1,Data.Processed.Monthly.STD.c1)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.c2,Data.Processed.Monthly.STD.c2)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.c3,Data.Processed.Monthly.STD.c3)
figure
errorbar(1:12,Data.Processed.Monthly.Mean.c4,Data.Processed.Monthly.STD.c4)
end
%end
else fprintf('Location does not exist \n');
end
else
disp('f=1 || f=2')
end
end

