solar
deg=phi*180/pi;
fQ = fft(Q')/365;
[qt,qy]=gradient(Q,t,phi*6.37);
fQy = fft(qy')/365;

figure(1);orient tall;preprint;
subplot(2,1,1)
tt=t/(365/12);    
clabel(contour([tt,tt+12],deg,[Q,Q],[0:50:600]),0:100:600);
labels=['21FEB';'21APR';'21JUN';'21AUG';'21OCT';'21DEC'];        
set(gca,'xticklabel',labels,'fontsize',9)
title('Q: Daily Insolation W m^{-2} day','fontsize',12)
ylabel latitude
grid

subplot(2,1,2)
clabel(contour([tt,tt+12],deg,[qy,qy],[-100:20:100]),-100:20:100);
labels=['21FEB';'21APR';'21JUN';'21AUG';'21OCT';'21DEC'];        
set(gca,'xticklabel',labels,'fontsize',9)
title('dQ/dy:  W m^{-2} day (1000km)^{-1}','fontsize',12)
ylabel latitude
grid

figure(2);orient landscape;preprint;

subplot(2,2,1)
plot(deg,abs(fQ(2,:)).^2,deg,abs(fQ(3,:)).^2);
axis([-90 90 0 20000]) 
linewide(2)
grid
xlabel latitude
ylabel ('W^2 m^{-4} day^2')
title( 'Variance of Q: Annual and Semiannual')
legend('A','S');
 
subplot(2,2,3)
plot(deg,real(fQ(2,:)),deg,imag(fQ(2,:)),deg,real(fQ(3,:)),deg,imag(fQ(3,:)))
axis([-90 90 -200 200])
linewide(2)
grid
xlabel latitude
ylabel ('W m^{-2} day')
title( 'Amplitude of Q: Sine and Cosine')
legend('A Cosine','A Sine','S Cosine','S Sine')

subplot(2,2,2)
plot(deg,abs(fQy(2,:)).^2,deg,abs(fQy(3,:)).^2);
axis([-90 90 0 500]) 
linewide(2)
grid
xlabel latitude
ylabel('W^2 m^{-4} day^2 (1000km)^{-2}')
title( 'Variance of dQ/dy')
legend('A','S');
 

subplot(2,2,4)
plot(deg,real(fQy(2,:)),deg,imag(fQy(2,:)),deg,real(fQy(3,:)),deg,imag(fQy(3,:)))
axis([-90 90 -25 25])
linewide(2)
grid
xlabel latitude
ylabel ('W m^{-2} day (1000km)^{-1}')
title( 'Amplitude of dQ/dy')
legend('A Cosine','A Sine','S Cosine','S Sine')

figure(3);orient tall;preprint
subplot(2,1,1)
plot(deg,abs(fQ(2,:)).^2,deg,abs(fQ(3,:)).^2);
axis([-10 10 0 400]) 
linewide(2)
grid
xlabel latitude
ylabel ('W^2 m^{-4} day^2')
title( 'Variance of Insolation Q: Annual and Semiannual')
legend('Annual','Semiannual');
 
subplot(2,1,2)
plot(deg,real(fQ(2,:)),deg,imag(fQ(2,:)),deg,real(fQ(3,:)),deg,imag(fQ(3,:)))
axis([-10 10 -30 30])
linewide(2)
grid
xlabel latitude
ylabel ('W m^{-2} day')
title( 'Amplitude of  Q: Sine and Cosine')
legend('Annual Cosine','Annual Sine','Semiannual Cosine','Semiannual Sine')