%% pwtpfit9b.m is designed to estimate the mean F10.7 of the Intercosmos
%% satellites used by the Truhlik et al. empirical model ie IK 19, 24 and 
%% 25. Furthermore, pwtpfit9b.m calculates the average temperatures from 
%% the Truhlik et al. electron temperature model as a function of magnetic 
%% latitude at the altitudes 550km, 900km, 1500km and 2500km for day 
%% (0900MLT-1600MLT) and night (2100MLT-0400MLT). The calculations are 
%% undertaken using the uncorrected F10.7 option, unlike pwtpfit1.m.
%%

clear all
close all

%Dates of IK observations (based on Fig. 1 from Advances Space Res. Truhlik 
%et al. 2000) - format is day, month and year.
ik19s=[1 1 1979];             %IK19 start date
ik19f=[31 12 1980];           %IK19 finsih date
ik24s=[1 10 1989];            %IK24 start date
ik24f=[30 9 1991];            %IK24 finsih date
ik25s=[1 1 1992];             %IK25 start date
ik25f=[30 6 1993];            %IK25 finsih date

ikstr=['19';'24';'25'];       %IK ID string

%Combined array
ikmat=[ik19s ik19f;ik24s ik24f;ik25s ik25f];

%Loop through for each satellite
nos=3;                        %number of satellites
for i=1:nos
 
 dnos=daynum(ikmat(i,1),ikmat(i,2),ikmat(i,3));        %initial day of i'th satellite data
 dnof=daynum(ikmat(i,4),ikmat(i,5),ikmat(i,6));        %final day of i'th satellite data
 years=ikmat(i,3);                                     %initial year of i'th satellite data
 yearf=ikmat(i,6);                                     %final year of i'th satellite data
 
 dc=daynum4(dnos,years,dnof,yearf);                    %number of days in i'th satellite data
 
 %Loop through each day to create required day number and year vector
 dnov=zeros(1,dc);            %day number vector
 yearv=zeros(1,dc);           %year number vector
 for j=1:dc
  if j==1
   dnov(j)=dnos;              %1st day is initial day number
   yearv(j)=years;            %1st year is initial year
  else
   [dnov(j),yearv(j)]=daynum3(dnov(j-1),yearv(j-1),1);      %calculate the next day number passed previous day number
  end
 end
 
 %Using the daynumber and year vector calculate the daily F10.7 at 0UT for each day 
 f107=f107day(zeros(1,dc),dnov,yearv);

 %Calculate the mean F10.7 for the period of observations from the i'th satellite
 eval(['ik',ikstr(i,:),'f107=mean(f107);'])
  
end


%% Calculate the base temperature profile from Truhlik et al. without F10.7 or seasonal effects.
%% Note that there are still variations in the model depending on day number (even, apparently,
%% when the seasonal correction option is turned off!). So, decided to average the profiles 
%% across each season with seasonal correction deselected.

%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);
tmpe1=temptruh(latmatm,timmlt,dno*svm,0*svm,0,0);

%Calculate the summer profiles
month=6;
dno=daynum(day,month,year);
tmpe2=temptruh(latmatm,timmlt,dno*svm,0*svm,0,0);

%Calculate the winter profiles
month=12;
dno=daynum(day,month,year);
tmpe3=temptruh(latmatm,timmlt,dno*svm,0*svm,0,0);

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

%Determine the average day and night profiles for each altitude
altv=[550 900 1500 2500];
for i=1: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

%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=[tmpe550d' tmpe900d' tmpe1500d' tmpe2500d']; 
nighttemps=[tmpe550n' tmpe900n' tmpe1500n' tmpe2500n'];
 
%Save the required results
save pwtpfit9b ik19f107 ik24f107 ik25f107 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)')

