%% pwtpfit10.m calculates the average temperatures from the Truhlik et al. 
%% and Brace and Theis electron temperature model as a function of magnetic 
%% latitude for day (0900MLT-1600MLT) and night (2100MLT-0400MLT) at the
%% heights of their respective satellite data bases. The calculations are 
%% undertaken in each season. 
%%
%% NOTE - pwtpfit10.m is based on the earlier pwtpfit6.m. It is designed to 
%%        address concerns raised by of the referees of my temperature 
%%        paper.
%%
%% April 1 2003 - changed the base F10.7 to 120 to agree with original
%% Titheridge temperature model.
%%
%% April 15 2003 - change the Go correction so that it is fractional change
%% rather than a absolute change.
%%
%% Author: Phillip A. Webb
%% Last updated: April 15 2003 
%%
%%

clear all
close all

%Base corrected F10.7 to be used.
f107b=120;

%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
day=21;
month=3;
year=2001;
dno=daynum(day,month,year);
tmpetr1=temptruh(latmatm,timmlt,dno*svm,f107b*svm,1,0);
tmpebt1=tempbrac(latmatm,timmlt,1*svm);

%Calculate the summer profiles
month=6;
dno=daynum(day,month,year);
tmpetr2=temptruh(latmatm,timmlt,dno*svm,f107b*svm,1,0);
tmpebt2=tempbrac(latmatm,timmlt,2*svm);

%Calculate the winter profiles
month=12;
dno=daynum(day,month,year);
tmpetr3=temptruh(latmatm,timmlt,dno*svm,f107b*svm,1,0);
tmpebt3=tempbrac(latmatm,timmlt,4*svm);

%Determine the average day and night profiles for each altitude
altvtr=[550 900 1500 2500];
for i=1:4
 alts=int2str(altvtr(i));
 for j=1:3
  eval(['tmpe',alts,'_',int2str(j),'=tmpetr',int2str(j),'(:,i);'])
  eval(['tmpe',alts,'_',int2str(j),'=reshape(tmpe',alts,'_',int2str(j),',size(svm));'])
  eval(['tmpe',alts,'_',int2str(j),'ni=mean(tmpe',alts,'_',int2str(j),'(:,[1:5 14:16]),2);'])
  eval(['tmpe',alts,'_',int2str(j),'di=mean(tmpe',alts,'_',int2str(j),'(:,6:13),2);'])
 end
end

altvbt=[400 600 1400 3000];
for i=1:4
 alts=int2str(altvbt(i));
 for j=1:3
  eval(['tmpe',alts,'_',int2str(j),'=tmpebt',int2str(j),'(:,i);'])
  eval(['tmpe',alts,'_',int2str(j),'=reshape(tmpe',alts,'_',int2str(j),',size(svm));'])
  eval(['tmpe',alts,'_',int2str(j),'ni=mean(tmpe',alts,'_',int2str(j),'(:,[1:5 14:16]),2);'])
  eval(['tmpe',alts,'_',int2str(j),'di=mean(tmpe',alts,'_',int2str(j),'(:,6:13),2);'])
 end
end

%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(altvtr(i));
 eval(['lat',alts,'=lshpar(lsh,-sv,(re+altvtr(i))*sv);'])
 eval(['imck=imag(lat',alts,');'])
 eval(['lat',alts,'(abs(imck)>0)=-1;'])
end

for i=1:4
 alts=int2str(altvbt(i));
 eval(['lat',alts,'=lshpar(lsh,-sv,(re+altvbt(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(altvtr(i));
 for j=1:3
  eval(['tmpe',alts,'_',int2str(j),'n=interp1(latmv,tmpe',alts,'_',int2str(j),'ni,lat',alts,');'])
  eval(['tmpe',alts,'_',int2str(j),'d=interp1(latmv,tmpe',alts,'_',int2str(j),'di,lat',alts,');'])
 end
end

for i=1:4
 alts=int2str(altvbt(i));
 for j=1:3
  eval(['tmpe',alts,'_',int2str(j),'n=interp1(latmv,tmpe',alts,'_',int2str(j),'ni,lat',alts,');'])
  eval(['tmpe',alts,'_',int2str(j),'d=interp1(latmv,tmpe',alts,'_',int2str(j),'di,lat',alts,');'])
 end
end

%Combine results for each season by altitude
altv=sort([altvbt altvtr]);
for j=1:3
 flsd=[];
 flsn=[];
 for k=1:length(altv)
  eval(['flsdi=tmpe',int2str(altv(k)),'_',int2str(j),'d;'])
  eval(['flsni=tmpe',int2str(altv(k)),'_',int2str(j),'n;'])
  flsd=[flsd;flsdi];
  flsn=[flsn;flsni];
 end
 eval(['tmped',int2str(j),'=flsd;']) 
 eval(['tmpen',int2str(j),'=flsn;']) 
end

%Delete non required arrays
clear tmpe*_1 tmpe*_2 tmpe*_3 tmpe*_1*i tmpe*_2*i tmpe*_3*i tmpe*_*d tmpe*_*n

%Base altitude at which to conduct Titheridge fit
h0=400;

%Reset ISIS-1 data 
nop=length(latd);
for j=1:3
 eval(['tmped',int2str(j),'(length(altv),:)=0.9*tmped',int2str(j),'(length(altv),:);'])
 eval(['tmpen',int2str(j),'(length(altv),:)=0.9*tmpen',int2str(j),'(length(altv),:);'])
end

%Calculate the L-shell at 300km latitude
re=6371;                      %Earth radi (km)
sv=ones(length(latd),1);
lsh=lshpar(-sv,deg2rad(latd'),(re+300)*sv);

sl=sin(deg2rad(latd')).^2;              %sin(lat)^2 values

%Loop through each sL value and correct required Brace and Theis values.
for i=1:nop
 sli=sl(i);                             %ith sL value
 sli(sli<0.1)=0.1;                      %reset sL values less than 0.1 to 0.1

 %Determine the required base temperatures and gradients
 sv8=ones(1,8);
 f107p=[f107b*ones(1,4) 90 85 115 150];   %required F107 values
 [t0dc,t0nc,g0dc,g0nc,t0di,t0ni,g0di,g0ni]=temppawf107(sli*sv8,f107p);

 tpedi=tehfn(t0di,g0di,[altvbt altvbt],lsh(i)*sv8,h0*sv8);      %required day temperatures
 tpeni=tehfn(t0ni,g0ni,[altvbt altvbt],lsh(i)*sv8,h0*sv8);      %required night temperatures

 tpedi(tpedi==0)=nan;                             %replace day zeros with NaN
 tpeni(tpeni==0)=nan;                             %replace night zeros with NaN

 tpedi=tpedi';
 tpeni=tpeni';

 poi=[1 3 5 8]';                        %Brace and Theis temperatures
 for j=1:3
  eval(['tmped',int2str(j),'(poi,i)=tmped',int2str(j),'(poi,i).*(tpedi(1:4)./tpedi(5:8));'])   %corrected i'th day temperature values
  eval(['tmpen',int2str(j),'(poi,i)=tmpen',int2str(j),'(poi,i).*(tpeni(1:4)./tpeni(5:8));'])   %corrected i'th night temperature values
 end
end

%Estimate the mean F10.7 for each altitude band from the Truhlik et al. results.
%These values are taken to be equal to the mean of the F10.7 values calculated 
%by pwtpfit9b.m, with IK19 values only used at the lower two altitude bands.
load pwtpfit9b
f107ik=[mean([ik19f107 ik24f107 ik25f107]) mean([ik19f107 ik24f107 ik25f107]) mean([ik24f107 ik25f107]) mean([ik24f107 ik25f107])];

%Loop through each sL value and correct required Truhlik et al. values.
for i=1:nop
 sli=sl(i);                             %ith sL value
 sli(sli<0.1)=0.1;                      %reset sL values less than 0.1 to 0.1

 %Determine the required base temperatures and gradients
 sv8=ones(1,8);
 f107p=[f107b*ones(1,4) f107ik];          %required F107 values
 [t0dc,t0nc,g0dc,g0nc,t0di,t0ni,g0di,g0ni]=temppawf107(sli*sv8,f107p);

 tpedi=tehfn(t0di,g0di,[altvtr altvtr],lsh(i)*sv8,h0*sv8);      %required day temperatures
 tpeni=tehfn(t0ni,g0ni,[altvtr altvtr],lsh(i)*sv8,h0*sv8);      %required night temperatures

 tpedi(tpedi==0)=nan;                   %replace day zeros with NaN
 tpeni(tpeni==0)=nan;                   %replace night zeros with NaN

 tpedi=tpedi';
 tpeni=tpeni';

 poi=[2 4 6 7]';                        %Truhlik et al. temperatures
 for j=1:3
  eval(['tmped',int2str(j),'(poi,i)=tmped',int2str(j),'(poi,i).*(tpedi(1:4)./tpedi(5:8));'])   %corrected i'th day temperature values
  eval(['tmpen',int2str(j),'(poi,i)=tmpen',int2str(j),'(poi,i).*(tpeni(1:4)./tpeni(5:8));'])   %corrected i'th night temperature values
 end
end


%Loop through each value latitude, conducting a non-linear fit using fmins.
c0=[1200 1];                            %starting parameters
options=optimset('display','off');

%Result matrices for G0 and T0 for day and night
nol=length(latd);
noa=length(altv);
t0md=zeros(nol,3);
g0md=t0md;
t0mn=t0md;
g0mn=t0md;

%Loop through each latitude and season
for i=1:nol
 for j=1:3
  altvdi=altv;                                    %i'th day altitude vector
  altvni=altv;                                    %i'th night altitude vector
  lshdi=lsh(i)*ones(1,noa);                       %i'th day L-shell vector
  lshni=lshdi;                                    %i'th night L-shell vector
  eval(['tmpdi=tmped',int2str(j),'(:,i);'])       %i'th day time temperatures
  eval(['tmpni=tmpen',int2str(j),'(:,i);'])       %i'th day time temperatures

  tmpdi=tmpdi';
  tmpni=tmpni';

  altvdi(isnan(tmpdi))=[];          %remove day altitude terms when no temperature information exists
  altvni(isnan(tmpni))=[];          %remove night altitude terms when no temperature information exists
  lshdi(isnan(tmpdi))=[];           %remove day L-shell terms when no temperature information exists
  lshni(isnan(tmpni))=[];           %remove night L-shell terms when no temperature information exists
  tmpdi(isnan(tmpdi))=[];           %remove day temperature terms when no temperature information exists
  tmpni(isnan(tmpni))=[];           %remove night temperature terms when no temperature information exists

  %Conduct fits to day data
  if isempty(altvdi) | length(altvdi)<2
   cdvd=[nan nan];
  else 
   cdvd=fminsearch('tmin',c0,options,altvdi,tmpdi,lshdi,h0);
  end

  %Conduct fits to night data
  if isempty(altvni) | length(altvni)<2
   cdvn=[nan nan];
  else 
   cdvn=fminsearch('tmin',c0,options,altvni,tmpni,lshni,h0);
  end

  %Record results
  t0md(i,j)=cdvd(1);
  g0md(i,j)=cdvd(2);
  t0mn(i,j)=cdvn(1);
  g0mn(i,j)=cdvn(2);
 end
end 

%Fit polynominals to fitted coefficents
lli=2;
for j=1:3
 eval(['t0dp',int2str(j),'=polyfit(sl(lli:nol),t0md(lli:nol,j),5);'])
 eval(['t0np',int2str(j),'=polyfit(sl(lli:nol),t0mn(lli:nol,j),5);'])
 eval(['g0dp',int2str(j),'=polyfit(sl(lli:nol),g0md(lli:nol,j),5);'])
 eval(['g0np',int2str(j),'=polyfit(sl(lli:nol),g0mn(lli:nol,j),5);'])
end

%Calculate fitted curves
slf=linspace(0,1);
for j=1:3
 eval(['t0df',int2str(j),'=polyval(t0dp',int2str(j),',slf);'])
 eval(['t0nf',int2str(j),'=polyval(t0np',int2str(j),',slf);'])
 eval(['g0df',int2str(j),'=polyval(g0dp',int2str(j),',slf);'])
 eval(['g0nf',int2str(j),'=polyval(g0np',int2str(j),',slf);'])
end

%Calculate the difference between each season and the seasonal mean.
t0d=0.25*(2*t0md(:,1)+t0md(:,2)+t0md(:,3));
t0n=0.25*(2*t0mn(:,1)+t0mn(:,2)+t0mn(:,3));
g0d=0.25*(2*g0md(:,1)+g0md(:,2)+g0md(:,3));
g0n=0.25*(2*g0mn(:,1)+g0mn(:,2)+g0mn(:,3));
dt0md=t0md-t0d*ones(1,3);
dt0mn=t0mn-t0n*ones(1,3);
dg0md=g0md./(g0d*ones(1,3));
dg0mn=g0mn./(g0n*ones(1,3));

for j=1:3
 eval(['dt0dp',int2str(j),'=polyfit(sl(lli:nol),dt0md(lli:nol,j),5);'])
 eval(['dt0np',int2str(j),'=polyfit(sl(lli:nol),dt0mn(lli:nol,j),5);'])
 eval(['dg0dp',int2str(j),'=polyfit(sl(lli:nol),dg0md(lli:nol,j),5);'])
 eval(['dg0np',int2str(j),'=polyfit(sl(lli:nol),dg0mn(lli:nol,j),5);'])
end

disp(' ')
fprintf('dt0dp1=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dt0dp1)
fprintf('dt0np1=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dt0np1)
fprintf('dg0dp1=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dg0dp1)
fprintf('dg0np1=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dg0np1)
fprintf('dt0dp2=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dt0dp2)
fprintf('dt0np2=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dt0np2)
fprintf('dg0dp2=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dg0dp2)
fprintf('dg0np2=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dg0np2)
fprintf('dt0dp3=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dt0dp3)
fprintf('dt0np3=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dt0np3)
fprintf('dg0dp3=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dg0dp3)
fprintf('dg0np3=[%6.3e %6.3e %6.3e %6.3e %6.3e %6.3e]\n',dg0np3)
disp(' ')

%Calculate fitted curves
for j=1:3
 eval(['dt0df',int2str(j),'=polyval(dt0dp',int2str(j),',slf);'])
 eval(['dt0nf',int2str(j),'=polyval(dt0np',int2str(j),',slf);'])
 eval(['dg0df',int2str(j),'=polyval(dg0dp',int2str(j),',slf);'])
 eval(['dg0nf',int2str(j),'=polyval(dg0np',int2str(j),',slf);'])
end

%Save results to file
save pwtpfit10
save pwsea d*0*p* 

figure(1)
subplot(2,2,1)
h1=plot(sl,t0md(:,1),'ro',sl,t0md(:,2),'g^',sl,t0md(:,3),'bs',slf,t0df1,'-r',slf,t0df2,'-g',slf,t0df3,'-b');
set(h1(4:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
set(h1(2),'markerfacecolor','g')
set(h1(3),'markerfacecolor','b')
grid
set(gca,'fontsize',10,'fontweight','bold')
xlabel('sL parameter','fontsize',10,'fontweight','bold')
ylabel('Day - T0 (K)','fontsize',10,'fontweight','bold')
legend(h1(1:3),'Equinox','Summer','Winter',2)

subplot(2,2,2)
h1=plot(sl,t0mn(:,1),'ro',sl,t0mn(:,2),'g^',sl,t0mn(:,3),'bs',slf,t0nf1,'-r',slf,t0nf2,'-g',slf,t0nf3,'-b');
set(h1(4:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
set(h1(2),'markerfacecolor','g')
set(h1(3),'markerfacecolor','b')
grid
set(gca,'fontsize',10,'fontweight','bold')
xlabel('sL parameter','fontsize',10,'fontweight','bold')
ylabel('Night - T0 (K)','fontsize',10,'fontweight','bold')

subplot(2,2,3)
h1=plot(sl,g0md(:,1),'ro',sl,g0md(:,2),'g^',sl,g0md(:,3),'bs',slf,g0df1,'-r',slf,g0df2,'-g',slf,g0df3,'-b');
set(h1(4:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
set(h1(2),'markerfacecolor','g')
set(h1(3),'markerfacecolor','b')
grid
set(gca,'fontsize',10,'fontweight','bold')
xlabel('sL parameter','fontsize',10,'fontweight','bold')
ylabel('Day - G0 (K/km)','fontsize',10,'fontweight','bold')

subplot(2,2,4)
h1=plot(sl,g0mn(:,1),'ro',sl,g0mn(:,2),'g^',sl,g0mn(:,3),'bs',slf,g0nf1,'-r',slf,g0nf2,'-g',slf,g0nf3,'-b');
set(h1(4:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
set(h1(2),'markerfacecolor','g')
set(h1(3),'markerfacecolor','b')
grid
set(gca,'fontsize',10,'fontweight','bold')
xlabel('sL parameter','fontsize',10,'fontweight','bold')
ylabel('Night - G0 (K/km)','fontsize',10,'fontweight','bold')

figure(2)
subplot(2,2,1)
h1=plot(sl,dt0md(:,1),'ro',sl,dt0md(:,2),'g^',sl,dt0md(:,3),'bs',slf,dt0df1,'-r',slf,dt0df2,'-g',slf,dt0df3,'-b');
set(h1(2:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
axis([0 1 -800 800])
grid
set(gca,'xtick',0:0.2:1)
set(gca,'fontsize',12,'fontweight','bold')
xlabel('{\it s_L}','fontsize',14,'fontweight','bold')
ylabel('Day  \Delta{\itT_o} (K)','fontsize',14,'fontweight','bold')

subplot(2,2,2)
h1=plot(sl,dt0mn(:,1),'ro',sl,dt0mn(:,2),'g^',sl,dt0mn(:,3),'bs',slf,dt0nf1,'-r',slf,dt0nf2,'-g',slf,dt0nf3,'-b');
set(h1(2:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
axis([0 1 -800 800])
grid
set(gca,'xtick',0:0.2:1)
set(gca,'fontsize',12,'fontweight','bold')
xlabel('{\it s_L}','fontsize',14,'fontweight','bold')
ylabel('Night  \Delta{\itT_o} (K)','fontsize',14,'fontweight','bold')
legend(h1(1:3),'Equinox','Summer','Winter',3)

subplot(2,2,3)
h1=plot(sl,dg0md(:,1),'ro',sl,dg0md(:,2),'g^',sl,dg0md(:,3),'bs',slf,dg0df1,'-r',slf,dg0df2,'-g',slf,dg0df3,'-b');
set(h1(2:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
axis([0 1 0 2])
grid
set(gca,'xtick',0:0.2:1)
set(gca,'fontsize',12,'fontweight','bold')
xlabel('{\it s_L}','fontsize',14,'fontweight','bold')
ylabel('Day  \Delta{\itG_o} (frac. change)','fontsize',14,'fontweight','bold')

subplot(2,2,4)
h1=plot(sl,dg0mn(:,1),'ro',sl,dg0mn(:,2),'g^',sl,dg0mn(:,3),'bs',slf,dg0nf1,'-r',slf,dg0nf2,'-g',slf,dg0nf3,'-b');
set(h1(2:6),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
axis([0 1 0 2])
grid
set(gca,'xtick',0:0.2:1)
set(gca,'fontsize',12,'fontweight','bold')
xlabel('{\it s_L}','fontsize',14,'fontweight','bold')
ylabel('Night  \Delta{\itG_o} (frac. change)','fontsize',14,'fontweight','bold')


