clear all
close all
%% read f107
f107=load('f10.7_monthly_1955-2005.txt');
[m,n]=size(f107);
f107=reshape(f107',1,m*n);

%% read ap
ap=load('ap_monthly_1955-2005.txt');
[m,n]=size(ap);
ap=reshape(ap',1,m*n);


%% temperature data
Tn=load('temperature_280km_monthly_1955-2005.txt');
[m,n]=size(Tn);
Tn=reshape(Tn',1,m*n);



n=1:length(Tn);
yy=n/12;

figure; 
set(gcf,'position',[440    67   589   731])

subplot(3,1,1)
plot(yy,f107)
set(gca,'ylim',[2000 10000],'fontsize',12)
ylabel('solar flux')

subplot(3,1,2)
plot(yy,Tn)


%% linearfit

aa0=ones(size(yy));
aa1=f107; aa2=ap; aa3=f107.*f107;
AA=[aa0' aa1' aa2' aa3'];
coef=AA\Tn';
Tn_fit=(AA*coef)';

hold on;
plot(yy,Tn_fit,'r')

legend('Tn','f107&ap-fitting')
legend boxoff

set(gca,'ylim',[600 1200],'fontsize',12)
ylabel('Tn (K)')

subplot(313)

plot(yy,Tn-Tn_fit)

%% fitting the trend
P=polyfit(yy(yy>15),Tn(yy>15)-Tn_fit(yy>15),1);
trend=polyval(P,yy(yy>15));
hold on;
plot(yy(yy>15),trend,'r')

P2=fix(P(1)*100)/10; %% 10 year trend
text(45,-50,['trend=' num2str(P2) 'K/10 yr'],'color','r')

set(gca,'ylim',[-100 100],'fontsize',12)
ylabel('\deltaTn (K)')

%% save fig
set(gcf,'PaperPositionMode','auto');
Fname=['trend'];
saveas(gcf,Fname,'epsc2')

