clear
close all
load tempPoland
%load tempEquator
czas=TIME(:,1)+day2jd(TIME(:,1),TIME(:,2),TIME(:,3))/365;
plot(czas,TEMP)
title('Przebieg srednich dobowych temperatur powietrza')
pause

k=0;
for i=1948:2006
   k=k+1; 
   ind=find(TIME(:,1)==i);
   Tyear(k)=mean(TEMP(ind));
   year(k)=i;
end    
plot(year,Tyear)
title('Srednia roczna temperatura powietrza')
A=polyfit(year,Tyear,1)
VAL=polyval(A,year);
hold on
plot(year,VAL,'r')
pause


k=0
for i=1948:2006
    k=k+1;
    for j=1:12
     ind=find((TIME(:,1)==i)&(TIME(:,2)==j));
     Tmonave(k,j)=mean(TEMP(ind));
     Amp(k,j)=max(TEMP(ind))-min(TEMP(ind));
    end 
end
figure
Tmon=mean(Tmonave);
Tstd=std(Tmonave);
mon=[1:12];
errorbar(mon,Tmon,Tstd)
title('Srednie miesieczne temperatury powietrza')
legend('TEMP','1 STD')

pause

figure
hist(TEMP,50)
title('Histogram temperatury powietrza')
pause

figure

k=0;
for i=1948:2006
   k=k+1; 
   ind=find((TIME(:,1)==i)&(TIME(:,2)==7));
   Tave(k)=mean(TEMP(ind));
end    
ind=find(TIME(:,2)==1);
[a,b]=hist(TEMP(ind),50);
bar(b,a);
hold on
ind=find(TIME(:,2)==6);
[a,b]=hist(TEMP(ind),50);
h=bar(b,a,'r');
set(h,'EdgeAlpha',0.0);
title('Histogram temperatury powietrza')
legend('Styczen','Czerwiec')
pause
%hist(Tave)
disp('-------------------------')
ind1=find(TIME(:,1)<1980);
ind2=find(TIME(:,1)>1980);
index=find(TEMP(ind1)<-10); 100*length(index)/length(ind1)
index=find(TEMP(ind2)<-10); 100*length(index)/length(ind2)
index=find(TEMP(ind1)>22); 100*length(index)/length(ind1)
index=find(TEMP(ind2)>22); 100*length(index)/length(ind2)
[prob1,t1]=hist(TEMP(ind1),75);
dT1=diff(t1);
[prob2,t2]=hist(TEMP(ind2),75);
dT2=diff(t2);
prob1=prob1/(length(ind1)*dT1(1));
prob2=prob2/(length(ind2)*dT2(1));

figure
bar(t1,prob1)
hold on
h=bar(t2,prob2,'r')
%set(h,'FaceAlpha',0.0);
set(h,'EdgeAlpha',0.0);
legend('1948-1979','1980-2006')
title('Rozklad pra-stwa wystepowania srednich dobowych temp')
xlabel('temperatura')
ylabel('gestosc prawdopodobienstwa')
[prob,t]=hist(TEMP,75);
dt=diff(t);
prob=prob/(dt(1)*length(TEMP));
val0=[mean(TEMP),2];
val=fminsearch('fun',val0,[],t1,prob1)
PROB=gauss(t,val(1),val(2));
hold on
plot(t,PROB,'g')
legend('GAUSS Fit','1948-1979','1980-2006')

return
figure
TT=[-40:0.1:40];
ave=0
STD=5
Tgaus=gauss(TT,ave,STD);
plot(TT,Tgaus)
ind=find((TT>ave-STD)&(TT<ave+STD));
trapz(TT(ind),Tgaus(ind))