%% tempjet3.m is my third implementation of J.E.Titheridge temperature 
%% model as set out in Feb 1998 JGR. tempjet3.m is different from the 
%% previous tempjets in that it only uses the Titheridge model and 
%% not the IRI at low altitudes to determine the electron temperatures.
%%

function [tmpi,tmpe,tmpn,tmpic,tmpec,tmpnc]=tempjet3(lonmat,latmat,dismat,timmat,year,month,day)

% lonmat = geographic longitude of points of interest (rads).
% latmat = geographic latitude of points of interest (rads).
% dismat = distance from center of Earth of points of interest (km).
% timmat = univerisal time (UT) of points of interest. 
% year = year matrix.
% month = month matrix.
% day = day matrix.
% tmpi = ion temperatures at point of interest (K).
% tmpe = electron temperatures at point of interest (K).
% tmpn = neutral temperatures at point of interest (K).
% tmpic = ion temperatures at conjugate point of interest (K).
% tmpec = electron temperatures at conjugate point of interest (K).
% tmpnc = neutral temperatures at conjugate point of interest (K).

sz=size(lonmat);                                       %size of input

%Convert all input to column vectors
lonmat=lonmat(:);
latmat=latmat(:);
dismat=dismat(:);
timmat=timmat(:);
year=year(:);
month=month(:);
day=day(:);

%Calculate required solar parameters
dno=daynum(day,month,year);                            %day number
[kp,apday,f107,f107a]=solarval(timmat,dno,year);

re=6371;                                               %Earth radi (km)
altmat=dismat-re;                                      %altitude of points of interest
[lonmatm,latmatm]=geo2mag3(lonmat,latmat);             %magnetic co-ordinates
[lonmatc,latmatc]=mag2geo3(lonmatm,-latmatm);          %conjugate geographic co-ordinates

lsh=lshpar(-ones(size(dismat)),latmatm,dismat);        %L-shell of point of interest

h03=300;                                               %base height from which to use JET model
h04=400;                                               %2nd base height        
  
% Estimate the required G0 and T0 values using the values in Table 2 and equation 19
t0di=[1230 2200 -3290 -0.26 -0.68];                    %T0 day co-efficents
g0di=[4.65 -8.55 4.14 -2.16 1.45];                     %G0 day co-efficents
t0ni=[985 -963 1125 -0.6 0.1];                         %T0 night co-efficents
g0ni=[0.756 -0.88 0.29 -2.63 1.84];                    %G0 night co-efficents
 
sl=1-((re+h03)./(re+altmat)).*cos(latmatm).^2;         %sl parameter
 
t0d=(t0di(1)+t0di(2).*sl+t0di(3).*sl.^2)./(1+t0di(4).*sl+t0di(5).*sl.^2);  %required T0 day values (eqn 19)
g0d=(g0di(1)+g0di(2).*sl+g0di(3).*sl.^2)./(1+g0di(4).*sl+g0di(5).*sl.^2);  %required G0 day values (eqn 19)
t0n=(t0ni(1)+t0ni(2).*sl+t0ni(3).*sl.^2)./(1+t0ni(4).*sl+t0ni(5).*sl.^2);  %required T0 night values (eqn 19)
g0n=(g0ni(1)+g0ni(2).*sl+g0ni(3).*sl.^2)./(1+g0ni(4).*sl+g0ni(5).*sl.^2);  %required G0 night values (eqn 19)

% Modify the values of G0 and T0 to reflect changes in solar activity
t0dc=3.4*(f107-120);                                   %T0 day correction (eqn 16a)
t0nc=3.4*(f107-120);                                   %T0 night correction (eqn 16a)
g0d=g0d.*t0d./(t0d+t0dc);                              %modifed G0 day (eqn 16b)
g0n=g0n.*t0n./(t0n+t0nc);                              %modifed G0 night (eqn 16b)
t0d=t0d+t0dc;                                          %modifed T0 day
t0n=t0n+t0nc;                                          %modifed T0 night
  
% Modify the values of G0 and T0 to reflect changes in magnetic activity
latc=(re+h03)./(lsh*re);                               %latitude check term
latc(latc>1)=1;                                        %reset latc terms greater than 1 to 1      
latm3=acos(sqrt(latc));                                %magnetic latitude at 300km altitude on L shell
latm3(latmatm<-1e-4)=-latm3(latmatm<-1e-4);            %correction for terms in southern hemisphere

l3i=rad2deg(latm3)>46;                                 %index of latm3 values above 46 degrees
if any(l3i)
 l3=latm3(l3i);                                        %latm3 values above 46 degrees
 l3=rad2deg(l3);                                       %convert l3 to degrees
 dk=37+1.33*(l3-46)-37*cos(0.135*(l3-46));             %Dk parameter (eqn 17)
 t0c=dk.*(kp(l3i)-1.5);                                %T0 correction
 g0c=(t0d(l3i)./(t0d(l3i)+t0c)).^(2.5);                %G0 correction term
 g0d(l3i)=g0d(l3i).*g0c;                               %modifed G0 day (eqn 18b)
 g0n(l3i)=g0n(l3i).*g0c;                               %modifed G0 night (eqn 18b)
 t0d(l3i)=t0d(l3i)+t0c;                                %modifed T0 day (eqn 18a)
 t0n(l3i)=t0n(l3i)+t0c;                                %modifed T0 night (eqn 18a)
end 
 
%Calculate the average temperature with the JET model using eqn 13a, 13b.                                                        
r0=1+h04/re;
rh=1+altmat/re;
hlog=log(altmat./h04);
bh=(0.05./(2*lsh-r0)).*(88+hlog.*(10.5-hlog));
dhr=((lsh-1)*re-h04)./r0^2-((lsh-1)*re-altmat)./rh.^2;
 
dhr(altmat<h04)=dhr(altmat<h04)./(1-0.01*dhr(altmat<h04)*r0^2);
TDay=t0d.*(1+bh.*(g0d./t0d).*dhr).^(2/7);                        %average electron day temperature
TNit=t0n.*(1+bh.*(g0n./t0n).*dhr).^(2/7);                        %average electron night temperature


%Determine time of local ground sunrise for point below field line at 300km altitude
[long3,latg3]=mag2geo3(lonmatm,latm3);                           %geographic co-ordinates of field lines at 300km
amon=month+(day-1)/31;                                           %day of year in months        
sundec=0.41015*sin((amon-3.7)*0.5236);                           %sun declination
cdlat=cos(latg3-sundec);
slong=1-cdlat./(cos(latg3).*cos(sundec));
srise=zeros(size(long3));                                             %sunrise matrix (0=no night)
srise(slong>=1)=12;                                                   %no day        
srise(slong>-1 & slong<1)=12-acos(slong(slong>-1 & slong<1))*3.82;    %time of ground sunrise
 
%Conjugate point values
[long3c,latg3c]=mag2geo3(lonmatm,-latm3);                             %geographic co-ordinates of field lines at 300km
cdlat=cos(latg3c-sundec);
slongc=1-cdlat./(cos(latg3c).*cos(sundec));
srisec=zeros(size(lonmat));                                                %sunrise matrix (0=no night)
srisec(slongc>=1)=12;                                                      %no day        
srisec(slongc>-1 & slongc<1)=12-acos(slongc(slongc>-1 & slongc<1))*3.82;   %time of ground sunrise
 
%Determine the required Titheridge correction parameters
polz=sqrt(cos(latm3));                       %polz correction parameter
dayz=0.7./(max(0.7,abs(slong)));             %dayz correction parameter
slz=(0.8*sl-1.4).*dayz;                      %slz correction parameter
aslong=abs(slong);                           %absolute values of slong
dayzc=0.7./(max(0.7,abs(slongc)));           %dayz conjugate correction parameter
slzc=(0.8*sl-1.4).*dayzc;                    %slz conjugate correction parameter
aslongc=abs(slongc);                         %conjugate absolute values of slong
 
%% To include conjugate heating effects, need to determine the temperature at the point
%% of interest, as well as at the conjugate point along the field line. Then take a 
%% weighted average of the two based on the height of the point of interest. 
dtr=1.2+0.5*sl;                              %delta tr (eqn 20b)
dts=dtr+0.9;                                 %delta ts (eqn 20b)
 
tdiur=zeros(length(srise(:)),2);             %tdiur results matrix
TNv=tdiur;                                   %night time results matrix
TDv=tdiur;                                   %day time results matrix
 
for k=1:2  
 if k==1
  sr=srise;                                  %sunrise of point of interest
  timl=timut2l(timmat,long3);                %local time of points of interest
  dayzk=dayz;
  slzk=slz;
  slongk=slong;
  aslongk=aslong;
 else 
  sr=srisec;
  timl=timut2l(timmat,long3c);               %local time of points of interest
  dayzk=dayzc;
  slzk=slzc;
  slongk=slongc;
  aslongk=aslongc;
 end 
 ts=min(1,0.15*sr);                          %ts parameter (eqn 20a)
  
 %Determine the local time at 300km altitude along the field line
 dt=timl;                                    %time difference from midnight
 dt(dt>12)=dt(dt>12)-24;                     %correct dt value for hours 'before midnight'
 
 %Use equations 20 to 22 to determine the temperature at the required local time
 dval=zeros(size(timl));                                                                       %required D values
 dval(timl<=12)=(sr(timl<=12)-0.5*ts(timl<=12)-timl(timl<=12))./dtr(timl<=12);                 %D values for morning (eqn 21b)
 dval(timl>12)=(timl(timl>12)-ts(timl>12)+sr(timl>12)-24)./(dtr(timl>12)+0.9*polz(timl>12));   %D values for afternoon (eqn 21c)
 
 tdiur(:,k)=1./(1+exp(3.2*dval.*polz));      %diurinal parameter (eqn 21a)
 if any(aslongk>0.7)
  tdiur((aslongk>0.7),k)=tdiur((aslongk>0.7),k).*dayzk(aslongk>0.7);  %=0 if all night
 end 
 if any(slongk<-0.7)
  tdiur((slongk<-0.7),k)=tdiur((slongk<-0.7),k)+1-dayzk(slongk<-0.7); %=1 if all day
 end 
  
 %% I have decided to add in a small correction factor to force the temperatures
 %% at local midnight and midday to be continious. Use the correction between 21
 %% and 3, 9 and 15 hours. For midnight, force TDiurn equal to 0 and for midday
 %% equal to 1.
 tc1=timl>21 | timl<3;                                      %check for any times between 21 and 3 hours
 if any(tc1)
  tdiur(tc1,k)=tdiur(tc1,k).*sin((pi/6)*timl(tc1)).^2;      %add in correction
 end
 tc2=timl>9 & timl<15;                                      %check for any times between 9 and 15 hours
 if any(tc2)
  tdiur(tc2,k)=tdiur(tc2,k).*sin((pi/6)*(timl(tc2)-12)).^2+cos((pi/6)*(timl(tc2)-12)).^2;      %add in correction
 end
  
 TDa=(timl-12.5)./(max((12.5-sr),3));                       %day correction factor
 TDv(:,k)=0.97+0.22*TDa.^2.*dayzk.*polz;                    %required 'day' temperature (eqn 22a)
 
 dtN=exp(dt./max(5,sr-3));                                  %night correction factor
 TNa=0.83+2.1*polz.*dtN./(14+dtN);                          %required 'night' temperature (eqn 22b)
 TNv(:,k)=min(TNa,2);                                       %eqn 22b
end 
 
xh2=0.5*(dismat-re)./((lsh-1)*re);                          %weight for conjugate term
xh1=1-xh2;                                                  %weight for local term 
 
TDiurn=tdiur(:,1).*xh1+tdiur(:,2).*xh2;                     %=0, night to 1, day 
TeN=TNit.*(TNv(:,1).*xh1+TNv(:,2).*xh2);                    %eqn 22a
TeD=TDay.*(TDv(:,1).*xh1+TDv(:,2).*xh2);                    %eqn 22b  
tmpe=TeN+(TeD-TeN).*TDiurn;                                 %estimated electron temperature from JET model (eqn 21a) 
 
TDiurnc=tdiur(:,1).*xh2+tdiur(:,2).*xh1;                    %=0, night to 1, day 
TeNc=TNit.*(TNv(:,1).*xh2+TNv(:,2).*xh1);                   %eqn 22a
TeDc=TDay.*(TDv(:,1).*xh2+TDv(:,2).*xh1);                   %eqn 22b  
tmpec=TeNc+(TeDc-TeNc).*TDiurnc;                            %estimated electron temperature from JET model (eqn 21a) 

%% To determine the ion temperature, use the method suggested by Titheridge (1998). 
%% Calculate the neutral temperature using the MSIS-90 model. For points above 750km
%% it is assumed that the neutral temperature at the point of interest is equal to 
%% the value at 750km altitude directly below the point of interest. Assume
%% that the parameter A varies between the night time value of 0.08 and the day time
%% value of 0.26 with the same variation given by eqn 21a.
svt=[1 1];
dism=dismat*svt;
dism(dism>(re+750))=re+750;
[oxnt,hynt,tmpnt]=modmsis9([lonmat lonmatc],[latmat latmatc],dism,timmat*svt,day*svt,month*svt,year*svt,f107a*svt,f107*svt,apday*svt);
tmpn=tmpnt(:,1);
tmpnc=tmpnt(:,2);

ht=0.75*tmpn-75;                                                      %ht parameter
timl=timut2l(timmat,lonmat);                                          %local time
amat=0.08*ones(size(tmpn));                                           %a parameter
amat(timl>=10 & timl<=17)=0.26;
amat(timl>4 & timl<10)=0.08+0.18*(timl(timl>4 & timl<10)-4)/6;
amat(timl>17 & timl<22)=0.26-0.18*(timl(timl>17 & timl<22)-17)/6;
tmpi=tmpn+(tmpe-tmpn)./(1+amat.*(ht./altmat)+(ht./altmat).^5);        %calculate ion temperature

ht=0.75*tmpnc-75;                                                     %ht parameter
timl=timut2l(timmat,lonmatc);                                         %local time
amat=0.08*ones(size(tmpn));                                           %a parameter
amat(timl>=10 & timl<=17)=0.26;
amat(timl>4 & timl<10)=0.08+0.18*(timl(timl>4 & timl<10)-4)/6;
amat(timl>17 & timl<22)=0.26-0.18*(timl(timl>17 & timl<22)-17)/6;
tmpic=tmpnc+(tmpec-tmpnc)./(1+amat.*(ht./altmat)+(ht./altmat).^5);    %calculate conjugate ion temperature

%Force any points whose electron or ion temperature is less than the neutral 
%temperature to be equal to the neutral temperature at that point.
tmpi(tmpi<tmpn)=tmpn(tmpi<tmpn);
tmpe(tmpe<tmpn)=tmpn(tmpe<tmpn);
tmpic(tmpic<tmpnc)=tmpnc(tmpic<tmpnc);
tmpec(tmpec<tmpnc)=tmpnc(tmpec<tmpnc);

%Reshape to original size
tmpe=reshape(tmpe,sz);
tmpec=reshape(tmpec,sz);
tmpi=reshape(tmpi,sz);
tmpic=reshape(tmpic,sz);
tmpn=reshape(tmpn,sz);
tmpnc=reshape(tmpnc,sz);



