%% pwtpfit11c.m is based on pwtpfit9d.m, but includes the 400km Alouette data
%% and thus allows the user to see the effects of including this data set to 
%% the fits obtained.
%% For ease of programming, the AE-C data is assumed to occur at 401km, to 
%% seperate it from the 400km Aloutte data.
%%
%% Author: Phillip A. Webb
%% Last updated: April 4 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

%Determine the required latitude and F10.7
latd
latdj=input('Choose one latitude? (has to be 35^o or greater) ');
j=latd==latdj;

%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);

%Loop through each latitude
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
 poi=altvdj<3000; 
 cdvd=fminsearch('tmin',c0,options,altvdj(poi),tmpdj(poi),lshdj(poi),h0);
end

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

%Record results
t0d=cdvd(1);
g0d=cdvd(2);
t0n=cdvn(1);
g0n=cdvn(2);

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);


%Plot the required results
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

%Required ISIS-1 temps
isisdj=(10/9)*tmpdj(10);
isisnj=(10/9)*tmpnj(10);

figure(1)
subplot(1,2,1)
h1=semilogy(tmped,altf,'g-',tmpdj(1:9),altvdj(1:9),'ro',isisdj,3000,'b^',0.9*isisdj,3000,'bd',0.8*isisdj-400,3000,'bs');
axis([0 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:5]),'markerfacecolor','b')
set(gca,'fontweight','bold','fontsize',14)
xlabel('{\it T_e} (K)','fontweight','bold','fontsize',14)
ylabel('Altitude (km)','fontweight','bold','fontsize',14)
title(['Day Lat=',int2str(latdj),'^o'],'fontweight','bold','fontsize',14)

subplot(1,2,2)
h1=semilogy(tmpen,altf,'g-',tmpnj(1:9),altvnj(1:9),'ro',isisnj,3000,'b^',0.9*isisnj,3000,'bd',0.8*isisnj-400,3000,'bs');
axis([0 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:5]),'markerfacecolor','b')
set(gca,'fontweight','bold','fontsize',14)
xlabel('{\it T_e} (K)','fontweight','bold','fontsize',14)
ylabel('Altitude (km)','fontweight','bold','fontsize',14)
title(['Night Lat=',int2str(latdj),'^o'],'fontweight','bold','fontsize',14)
legend(h1,'Fit','Other','ISIS','ISIS 10%','ISIS Koh.',4)

figure(2)
subplot(1,2,1)
h1=semilogy(tmped,altf,'g-',tmpdj(1:9),altvdj(1:9),'ro',isisdj,3000,'b^',0.9*isisdj,3000,'bd',0.8*isisdj-400,3000,'bs');
axis([0 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:5]),'markerfacecolor','b')
set(gca,'fontweight','bold','fontsize',14)
xlabel('{\it T_e} (K)','fontweight','bold','fontsize',14)
ylabel('Altitude (km)','fontweight','bold','fontsize',14)
title(['Day Lat=',int2str(latdj),'^o'],'fontweight','bold','fontsize',14)
legend(h1,'Fit','Other','ISIS','ISIS 10%','ISIS Koh.',2)
text(200+tmpdj(1),altvdj(1),'AE-C','fontweight','bold','fontsize',12)
text(tmpdj(2)-1200,altvdj(2)+25,'Alouette','fontweight','bold','fontsize',12)
text(200+tmpdj(3),altvdj(3),'AEROS','fontweight','bold','fontsize',12)
text(tmpdj(4)-350,altvdj(4)+50,'IK','fontweight','bold','fontsize',12)
text(200+tmpdj(5),altvdj(5),'IK','fontweight','bold','fontsize',12)
text(200+tmpdj(6),altvdj(6),'Alouette','fontweight','bold','fontsize',12)
text(200+tmpdj(7),altvdj(7)-100,'ISIS-2','fontweight','bold','fontsize',12)
text(-300+tmpdj(8),altvdj(8)+100,'IK','fontweight','bold','fontsize',12)
text(-500+tmpdj(9),altvdj(9),'IK','fontweight','bold','fontsize',12)
text(tmpdj(10)-1000,altvdj(10)+300,'ISIS-2','fontweight','bold','fontsize',12)

subplot(1,2,2)
h1=semilogy(tmpen,altf,'g-',tmpnj(1:9),altvnj(1:9),'ro',isisnj,3000,'b^',0.9*isisnj,3000,'bd',0.8*isisnj-400,3000,'bs');
axis([0 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:5]),'markerfacecolor','b')
set(gca,'fontweight','bold','fontsize',14)
xlabel('{\it T_e} (K)','fontweight','bold','fontsize',14)
ylabel('Altitude (km)','fontweight','bold','fontsize',14)
title(['Night Lat=',int2str(latdj),'^o'],'fontweight','bold','fontsize',14)
text(200+tmpnj(1),altvdj(1),'AE-C','fontweight','bold','fontsize',12)
text(tmpnj(2)-1100,altvdj(2)+25,'Alouette','fontweight','bold','fontsize',12)
text(200+tmpnj(3),altvdj(3),'AEROS','fontweight','bold','fontsize',12)
text(tmpnj(4)-350,altvdj(4)+50,'IK','fontweight','bold','fontsize',12)
text(200+tmpnj(5),altvdj(5),'IK','fontweight','bold','fontsize',12)
text(200+tmpnj(6),altvdj(6),'Alouette','fontweight','bold','fontsize',12)
text(-800+tmpnj(7),altvdj(7)-100,'ISIS-2','fontweight','bold','fontsize',12)
text(200+tmpnj(8),altvdj(8)+100,'IK','fontweight','bold','fontsize',12)
text(200+tmpnj(9),altvdj(9),'IK','fontweight','bold','fontsize',12)
text(tmpnj(10),altvdj(10)+300,'ISIS-2','fontweight','bold','fontsize',12)
