%% pwtpfit11.m is designed to fit Titheridges eqn 13a and 13b temperatures
%% obtained from both electron temperature models found within the 
%% IRI - the Brace and Theis model used in my PhD work, and the more
%% recent Truhlik et al. empirical model.
%% pwtpfit11.m is based on pwtpfit9.m, which uses the F10.7 uncorrected 
%% temperature values from the Truhlik et al. model and then corrects them 
%% to F10.7=120 using the same approach as used with the Brace and Theis 
%% model. However, unlike pwtpfit9.m the 400km Aloutette data is included. 
%% As before the EXOS-D are still ignored. Also, the ISIS-1 data is 
%% corrected from the original by assume only a 10% reduction in value, 
%% instead of using the Kohnlein correction. The Brace and Theis 
%% temperatures are calculated directly from the IRI.
%% For ease of programming, the AE-C data is assumed to occur at 401km, to 
%% seperate it from the 400km Aloutte data.
%% pwtpfit11.m is designed to address concerns raised by of the referees
%% of my temperature paper.
%%
%% Author: Phillip A. Webb
%% Last updated: April 1 2003 
%%

clear all
close all

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

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

%Load in Titheridge satellite data
load DayTemps.asc -ascii
load NightTemps.asc -ascii
latd=DayTemps(:,1);
dayt1t=DayTemps(:,[2 3]);
nitt1t=NightTemps(:,[2 3]);

%Load in the Brace and Theis data
load pwtpfit4
dayt1bt=daytemps;
nitt1bt=nighttemps;

%reset ISIS-1 data 
nop=length(latd);
dayt1bt(5:nop,4)=0.9*dayt1bt(5:nop,4);
nitt1bt(5:nop,4)=0.9*nitt1bt(5:nop,4);

%Combine the Titheridge and Brace and Theis data by altitude
dayt1=[dayt1t(:,1) dayt1bt(:,1:2) dayt1t(:,2) dayt1bt(:,3:4)];
nitt1=[nitt1t(:,1) nitt1bt(:,1:2) nitt1t(:,2) nitt1bt(:,3:4)];
altv1=[400 401 600 1000 1400 3000];

%Replace zero terms with NaN
dayt1(dayt1==0)=nan;
nitt1(nitt1==0)=nan;

%Load in the required Truhlik et al. results.
load pwtpfit9b
altv2=[550 900 1500 2500];         %altitudes of model values
noa=length(altv2);                 %number of altitude bins
dayt2=daytemps;
nitt2=nighttemps;

%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.
f107ik=[mean([ik19f107 ik24f107 ik25f107]) mean([ik19f107 ik24f107 ik25f107]) mean([ik24f107 ik25f107]) mean([ik24f107 ik25f107])];

%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, as well as the Titheridge, 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

 dayt1i=dayt1(i,:);                     %i'th day temperature values
 nitt1i=nitt1(i,:);                     %i'th night temperature values

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

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

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

 dayt1(i,:)=dayt1i.*(tpedi(1:6)./tpedi(7:12));    %corrected i'th day temperature values
 nitt1(i,:)=nitt1i.*(tpeni(1:6)./tpeni(7:12));    %corrected i'th night temperature values
end

%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

 dayt2i=dayt2(i,:);                     %i'th day temperature values
 nitt2i=nitt2(i,:);                     %i'th night temperature values

 %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,[altv2 altv2],lsh(i)*sv8,h0*sv8);         %required day temperatures
 tpeni=tehfn(t0ni,g0ni,[altv2 altv2],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

 dayt2(i,:)=dayt2i.*(tpedi(1:4)./tpedi(5:8));     %corrected i'th day temperature values
 nitt2(i,:)=nitt2i.*(tpeni(1:4)./tpeni(5:8));     %corrected i'th night temperature values
end

%Combine the two data sets into one and order by altitude
altv=[altv1 altv2];                     %combine altitude vector
dayt=[dayt1 dayt2];                     %combine day temperatures
nitt=[nitt1 nitt2];                     %combine night temperatures
[altv,ogi]=sort(altv);                  %sort vectors by altitude
dayt=dayt(:,ogi);                       %re-order day temperatures
nitt=nitt(:,ogi);                       %re-order night temperatures

%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,1);
g0md=t0md;
t0mn=t0md;
g0mn=t0md;

%Loop through each latitude
for j=1:nol
 altvdj=altv;                      %j'th day altitude vector
 altvnj=altv;                      %j'th night altitude vector
 lshdj=lsh(j)*ones(1,noa);         %j'th day L-shell vector
 lshnj=lshdj;                      %j'th night L-shell vector
 tmpdj=dayt(j,:);                  %j'th day time temperatures
 tmpnj=nitt(j,:);                  %j'th night time temperatures

 altvdj(isnan(tmpdj))=[];          %remove day altitude terms when no temperature information exists
 altvnj(isnan(tmpnj))=[];          %remove night altitude terms when no temperature information exists
 lshdj(isnan(tmpdj))=[];           %remove day L-shell terms when no temperature information exists
 lshnj(isnan(tmpnj))=[];           %remove night L-shell terms when no temperature information exists
 tmpdj(isnan(tmpdj))=[];           %remove day temperature terms when no temperature information exists
 tmpnj(isnan(tmpnj))=[];           %remove night temperature terms when no temperature information exists

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

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

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

end 

%Fit polynominals to fitted coefficents
lli=2;
t0dp=polyfit(sl(lli:nol),t0md(lli:nol),4);
t0np=polyfit(sl(lli:nol),t0mn(lli:nol),4);
g0dp=polyfit(sl(lli:nol),g0md(lli:nol),4);
g0np=polyfit(sl(lli:nol),g0mn(lli:nol),4);

%Calculate fitted curves
slf=linspace(0.05,1);
latdf=rad2deg(asin(sqrt(slf)));
t0df=polyval(t0dp,slf);
t0nf=polyval(t0np,slf);
g0df=polyval(g0dp,slf);
g0nf=polyval(g0np,slf);

disp(' ')
fprintf('t0dp=[%6.3e %6.3e %6.3e %6.3e %6.3e]\n',t0dp)
fprintf('t0np=[%6.3e %6.3e %6.3e %6.3e %6.3e]\n',t0np)
fprintf('g0dp=[%6.3e %6.3e %6.3e %6.3e %6.3e]\n',g0dp)
fprintf('g0np=[%6.3e %6.3e %6.3e %6.3e %6.3e]\n',g0np)
disp(' ')

%Save results to file
save pwtpfit11
save pwt0g0 t0dp t0np g0dp g0np

figure(1)
subplot(2,1,1)
h1=plot(sl,t0md,'ro',slf,t0df,'-r',sl,t0mn,'^b',slf,t0nf,'b-');
set(h1([2 3 4]),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
grid
set(gca,'fontsize',14,'fontweight','bold')
xlabel('{\it s_L}','fontsize',14,'fontweight','bold')
ylabel('{\it T_o} (K)','fontsize',14,'fontweight','bold')
%legend(h1([1 3]),'Day','Night',2)

subplot(2,1,2)
h1=plot(sl,g0md,'ro',slf,g0df,'-r',sl,g0mn,'^b',slf,g0nf,'b-');
set(h1([2 3 4]),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
grid
set(gca,'fontsize',14,'fontweight','bold')
xlabel('{\it s_L}','fontsize',14,'fontweight','bold')
ylabel('{\it G_o} (K/km)','fontsize',14,'fontweight','bold')
legend(h1([1 3]),'Day','Night',1)

figure(2)
subplot(2,1,1)
h1=plot(latd,t0md,'ro',latdf,t0df,'-r',latd,t0mn,'^b',latdf,t0nf,'b-');
set(h1([2 3 4]),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
grid
set(gca,'fontsize',10,'fontweight','bold')
xlabel('Latitude at 300km','fontsize',10,'fontweight','bold')
ylabel('T0 (K)','fontsize',10,'fontweight','bold')
legend(h1([1 3]),'Day','Night',2)

subplot(2,1,2)
h1=plot(latd,g0md,'ro',latdf,g0df,'-r',latd,g0mn,'^b',latdf,g0nf,'b-');
set(h1([2 3 4]),'linewidth',1.5)
set(h1(1),'markerfacecolor','r')
grid
set(gca,'fontsize',10,'fontweight','bold')
xlabel('Latitude at 300km','fontsize',10,'fontweight','bold')
ylabel('G0 (K/km)','fontsize',10,'fontweight','bold')
legend(h1([1 3]),'Day','Night',2)

%Required latitude values to plot
latdv=20:10:90;
lati=0;

%Plot the results for each 10^o
for pi=1:2
 for k=1:4

  if k==1
   figure(pi+2)
  end

  %Determine required j'th values
  lati=lati+1;
  j=find(latd==latdv(lati));
  
  altvdj=altv;                      %j'th day altitude vector
  altvnj=altv;                      %j'th night altitude vector
  lshdj=lsh(j)*ones(1,noa);         %j'th day L-shell vector
  lshnj=lshdj;                      %j'th night L-shell vector
  tmpdj=dayt(j,:);                  %j'th day time temperatures
  tmpnj=nitt(j,:);                  %j'th night time temperatures

  altvdj(isnan(tmpdj))=[];          %remove day altitude terms when no temperature information exists
  altvnj(isnan(tmpnj))=[];          %remove night altitude terms when no temperature information exists
  lshdj(isnan(tmpdj))=[];           %remove day L-shell terms when no temperature information exists
  lshnj(isnan(tmpnj))=[];           %remove night L-shell terms when no temperature information exists
  tmpdj(isnan(tmpdj))=[];           %remove day temperature terms when no temperature information exists
  tmpnj(isnan(tmpnj))=[];           %remove night temperature terms when no temperature information exists

  %Determine the fit profiles
  [t0d,t0n,g0d,g0n]=temppawt0g0(sl(j));
  altf=linspace(350,5000,100);
  svf=ones(size(altf));
  tmped=tehfn(t0d*svf,g0d*svf,altf,lshdj(1)*svf,h0*svf);
  tmpen=tehfn(t0n*svf,g0n*svf,altf,lshnj(1)*svf,h0*svf);

  altv1a=[altv1(2:3) altv1(5:6)];
  altv1b=[altv1(1) altv1(4)];
  bti=[];
  tit=[];
  ti=[];

  for i=1:4
   poi=find(altvdj==altv1a(i));
   if ~isempty(poi)
    bti=[bti poi];
   end 
  end
  
  for i=1:4
   poi=find(altvdj==altv2(i));
   if ~isempty(poi)
    ti=[ti poi];
   end 
  end

  for i=1:2
   poi=find(altvdj==altv1b(i));
   if ~isempty(poi)
    tit=[tit poi];
   end 
  end
  
  altv1a2=[altv1(2:3) altv1(5:6)];
  altv1b2=[altv1(1) altv1(4)];
  bti2=[];
  tit2=[];
  ti2=[];
  
  for i=1:4
   poi=find(altvnj==altv1a(i));
   if ~isempty(poi)
    bti2=[bti2 poi];
   end 
  end
  
  for i=1:4
   poi=find(altvnj==altv2(i));
   if ~isempty(poi)
    ti2=[ti2 poi];
   end 
  end

  for i=1:2
   poi=find(altvnj==altv1b(i));
   if ~isempty(poi)
    tit2=[tit2 poi];
   end 
  end

  subplot(2,2,k)
  h1=semilogy(tmped,altf,'g-',tmpdj(bti),altvdj(bti),'r^',tmpdj(ti),altvdj(ti),'bo',tmpdj(tit),altvdj(tit),'md');
  hold on
  h2=semilogy(tmpen,altf,'g-',tmpnj(bti2),altvnj(bti2),'r^',tmpnj(ti2),altvnj(ti2),'bo',tmpnj(tit2),altvnj(tit2),'md');
  hold off
  axis([500 5000 300 5000])
  grid
  set(gca,'ytick',[500 1000 1500 2000:1000:5000])
  set(h1(1),'linewidth',2)
  set(h1(2),'markerfacecolor','r')
  set(h1(3),'markerfacecolor','b')
  if length(h1)==4
   set(h1(4),'markerfacecolor','m')
  end 
  set(h2,'linewidth',2)
  set(gca,'fontweight','bold','fontsize',14)
  xlabel('{\it T_e} (K)','fontweight','bold','fontsize',14)
  ylabel('Altitude (km)','fontweight','bold','fontsize',14)
  text(3100,500,['Lat_{m300}=',int2str(latdv(lati)),'^o'],'fontweight','bold','fontsize',14)
  if k==1 & pi==1
   if length(h1)==3
    %legend([h1(2:3);h2(2:3)],'Brace and Theis (day)','Truhlik et al. (day)','Brace and Theis (night)','Truhlik et al. (night)',1)
    legend([h1(2:3)],'BTM','ISTM',1)
   else
    %legend([h1(2:4);h2(2:4)],'Brace and Theis (day)','Truhlik et al. (day)','Titheridge (day)','Brace and Theis (night)','Truhlik et al. (night)','Titheridge (night)',1)
    legend([h1(2:4)],'BTM','ISTM','Alouette',1)
   end
  end

 end
end
