%% pwtpfit4.m calculates the average temperatures from the Brace and Theis 
%% electron temperature model as a function of magnetic latitude at the
%% altitudes 400km, and 1400km and 3000km for day (0900MLT-1600MLT) and
%% night (2100MLT-0400MLT). The results are saved in the file pwtpfit4.mat.
%% Also, the day and night temperatures at 600km are calculated as per
%% the IRI model, which is based on AEROS satellite data. Note that this
%% profile has no seasonal variation.
%%

clear all
close all

%Set up the required latitude and time arrays
latd=10:5:90;
latmv=deg2rad(latd);
timv=[0:4 9:16 21:23];
[timmlt,latmatm]=meshgrid(timv,latmv);

svm=ones(size(timmlt));

%Calculate the spring and fall profiles
tmpe1=tempbrac(latmatm,timmlt,1*svm);

%Calculate the summer profiles
tmpe2=tempbrac(latmatm,timmlt,2*svm);

%Calculate the winter profiles
tmpe3=tempbrac(latmatm,timmlt,4*svm);

%Average the seasonal profiles
tmpe=0.25*(2*tmpe1+tmpe2+tmpe3);

%Select out required results for 400km, 1400km and 3000km.
tmpe1=tmpe1(:,2:4);
tmpe2=tmpe2(:,2:4);
tmpe3=tmpe3(:,2:4);

%Determine the average day and night profiles for each altitude
altv=[400 600 1400 3000];
for i=[1 3 4]
 alts=int2str(altv(i));
 eval(['tmpe',alts,'=tmpe(:,i);'])
 eval(['tmpe',alts,'=reshape(tmpe',alts,',size(svm));'])
 eval(['tmpe',alts,'ni=mean(tmpe',alts,'(:,[1:5 14:16]),2);'])
 eval(['tmpe',alts,'di=mean(tmpe',alts,'(:,6:13),2);'])
end

%Calculate the 600km day and night profiles
[tmpe600di,tmpe600ni]=tempaeros(latmv');

%Determine the latitudes of each required altitude when given as a function of latitude at 300km.
%Note that for some points -1 will be returned at low latitude to represent non-physical points.
re=6371;
sv=ones(size(latmv));
lsh=lshpar(-sv,latmv,(re+300)*sv);
for i=1:4
 alts=int2str(altv(i));
 eval(['lat',alts,'=lshpar(lsh,-sv,(re+altv(i))*sv);'])
 eval(['imck=imag(lat',alts,');'])
 eval(['lat',alts,'(abs(imck)>0)=-1;'])
end

%Use interpolation to determine the required temperatures at each altitude
for i=1:4
 alts=int2str(altv(i));
 eval(['tmpe',alts,'n=interp1(latmv,tmpe',alts,'ni,lat',alts,');'])
 eval(['tmpe',alts,'d=interp1(latmv,tmpe',alts,'di,lat',alts,');'])
end

%Record the required results
daytemps=[tmpe400d' tmpe600d' tmpe1400d' tmpe3000d']; 
nighttemps=[tmpe400n' tmpe600n' tmpe1400n' tmpe3000n'];
 
%Save the required results
save pwtpfit4 latd daytemps nighttemps

figure(1)
subplot(2,1,1)
plot(latd',daytemps(:,1),'o-',latd',daytemps(:,2),'^-',latd',daytemps(:,3),'s-',latd',daytemps(:,4),'d-')
grid
xlabel('Latitude at 300km')
ylabel('Day Temperature (K)')
subplot(2,1,2)
plot(latd',nighttemps(:,1),'o-',latd',nighttemps(:,2),'^-',latd',nighttemps(:,3),'s-',latd',nighttemps(:,4),'d-')
grid
xlabel('Latitude at 300km')
ylabel('Night Temperature (K)')

