%% temppaw2.m is just temppaw.m, but requires the input of the neutral temperature
%% at the point of interest and its conjugate point, as well as the solar parameters
%% for 3-hourly Kp and F10.7 flux for the previous day. 
%%
%% Author: Phillip A. Webb
%% Last updated: June 9 2004                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
%%

function [tmpi,tmpe,tmpic,tmpec,tdiur,tdiurc]=temppaw2(lonmat,latmat,dismat,timmat,year,month,day,tmpn,tmpnc,kp,f107)

% 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 (fractional hours). 
% year = year matrix.
% month = month matrix.
% day = day matrix.
% tmpn = neutral temperatures at point of interest (K).
% tmpnc = neutral temperatures at conjugate point of interest (K).
% kp = 3-hourly Kp
% f107 = F10.7 flux
%
% tmpi = ion temperatures at point of interest (K).
% tmpe = electron temperatures at point of interest (K).
% tmpic = ion temperatures at conjugate point of interest (K).
% tmpec = electron temperatures at conjugate point of interest (K).
% tdiur = diurinal parameter
% tdiurc = conjugate diurinal parameter

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(:);
tmpn=tmpn(:);
tmpnc=tmpnc(:);
kp=kp(:);
f107=f107(:);

nop=prod(sz);                                          %number of points of interest
[lonmatm,latmatm]=geo2mag3(lonmat,latmat);             %magnetic co-ordinates
[lonmatc,latmatc]=mag2geo3(lonmatm,-latmatm);          %conjugate geographic co-ordinates
lsh=lshpar(-ones(nop,1),latmatm,dismat);               %L-shell of point of interest
h03=300;                                               %base height from which to use JET model
h04=400;                                               %2nd base height        
re=6371;                                               %Earth radi (km)
altmat=dismat-re;                                      %altitude of points of interest

dno=daynum(day,month,year);                            %day number

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

[long3,latg3]=mag2geo3(lonmatm,latm3);                 %geographic co-ordinates of field lines at 300km
[long3c,latg3c]=mag2geo3(lonmatm,-latm3);              %conjugate geographic co-ordinates of field lines at 300km
  
sl=1-((re+h03)./dismat).*cos(latmatm).^2;              %sl parameter
  
%%Calculate required G0 and T0 values.
[t0d,t0n,g0d,g0n]=temppawt0g0(sl);                     

%% Modify the values of G0 and T0 to reflect changes in magnetic activity
l3i=abs(rad2deg(latm3))>46;                            %index of latm3 values above 46 degrees
if any(l3i)
 l3=abs(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)
 t0dif=dk.*(kp(l3i)-1.5);                              %T0 correction
 g0ddif=(t0d(l3i)./(t0d(l3i)+t0dif)).^(2.5);           %day G0 correction term
 g0ndif=(t0n(l3i)./(t0n(l3i)+t0dif)).^(2.5);           %night G0 correction term 
 g0d(l3i)=g0d(l3i).*g0ddif;                            %modifed G0 day (eqn 18b)
 g0n(l3i)=g0n(l3i).*g0ndif;                            %modifed G0 night (eqn 18b)
 t0d(l3i)=t0d(l3i)+t0dif;                              %modifed T0 day (eqn 18a)
 t0n(l3i)=t0n(l3i)+t0dif;                              %modifed T0 night (eqn 18a)
end 

%% Modify the values of G0 and T0 to reflect changes in solar activity.
[fct0d,fct0n,fcg0d,fcg0n]=temppawf107(sl,f107);        %calculate fractional change due to F10.7
t0d=t0d.*fct0d;                                        %corrected day T0 
t0n=t0n.*fct0n;                                        %corrected night T0 
g0d=g0d.*fcg0d;                                        %corrected day G0 
g0n=g0n.*fcg0n;                                        %corrected night G0 

% Modify the values of G0 and T0 to reflect changes due to season. Because
% the seasonal change is different in each hemisphere, need to separate the
% day and night into northern and southern hemisphere values.
t0dc=t0d;                          %conjugate copy of t0d
t0nc=t0n;                          %conjugate copy of t0n
g0dc=g0d;                          %conjugate copy of g0d
g0nc=g0n;                          %conjugate copy of g0n

%Calculate seasonal variation from mean.
[dt0d,dt0n,dg0d,dg0n,dt0dc,dt0nc,dg0dc,dg0nc]=temppawsea(sl,dno);

%Modify the values of T0 and G0 due to the season
t0d=t0d+dt0d;
t0dc=t0dc+dt0dc;
t0n=t0n+dt0n;
t0nc=t0nc+dt0nc;
g0d=g0d+dg0d;
g0dc=g0dc+dg0dc;
g0n=g0n+dg0n;
g0nc=g0nc+dg0nc;

%%Calculate the average temperature with the JET model using eqn 13a, 13b.                                                        
r0=1+h04/re;
rh=1+altmat/re;
hlog=log10(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;

dhrd=1+bh.*(g0d./t0d).*dhr;
dhrdc=1+bh.*(g0dc./t0dc).*dhr;
dhrn=1+bh.*(g0n./t0n).*dhr;
dhrnc=1+bh.*(g0nc./t0nc).*dhr;

%Rather than use the Titheridge correction for points below h0,
%reset the points inside the 2/7 term to 0, which will return
%a temperature of 0. This will be later reset to the neutral
%temperature. This required to allow temperatures below h0 to 
%be obtained down to 350km.
dhrd(dhrd<0)=0;
dhrdc(dhrdc<0)=0;
dhrn(dhrn<0)=0;
dhrnc(dhrnc<0)=0;

TDay=t0d.*dhrd.^(2/7);                  %average electron day temperature (eqn 13)
TDayc=t0dc.*dhrdc.^(2/7);               %average conjugate electron day temperature (eqn 13)
TNit=t0n.*dhrn.^(2/7);                  %average electron night temperature (eqn 13)
TNitc=t0nc.*dhrnc.^(2/7);               %average conjugate electron night temperature (eqn13)

%Determine time of local ground sunrise for point below field line at 300km altitude
[sr,ss,slong,sund]=sunrs(latg3,day,month,year);        %local ground sunrise          
[src,ssc,slongc,sundc]=sunrs(latg3c,day,month,year);   %conjugate local ground sunrise          

slg=sin(latg3).^2;                           %geographic equivalent of sl
slgc=sin(latg3c).^2;                         %conjugate geographic equivalent of sl
 
%% 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*slg;                             %delta tr (eqn 20b)
dtrc=1.2+0.5*slgc;                           %conjugate delta tr (eqn 20b)
   
% New variation in tr, determined by tempss.m. 
tr=-0.61*sr+3.12;                            %modified tr parameter (eqn 20a)
trc=-0.61*src+3.12;                          %modified conjugate tr parameter (eqn 20a)
ts=-2*tr;                                    %modified ts parameter (eqn 20a)
tsc=-2*trc;                                  %modified conjuagte ts parameter (eqn 20a)
 
%Determine the local time at 300km altitude along the field line
timl=timut2l(timmat,long3);                  %local time of points of interest
timlc=timut2l(timmat,long3c);                %local time of conjugate points of interest
 
%Use equations 20 to 22 to determine the temperature at the required local time
dval=zeros(nop,1);                                                    %required D values
dvalc=dval;                                                           %conjugate D values
poi1=timl<=12;                                                        %points before midday
poi2=timl>12;                                                         %points after midday
dval(poi1)=(sr(poi1)+tr(poi1)-timl(poi1))./dtr(poi1);                 %D values for morning (eqn 21b)
dval(poi2)=(timl(poi2)-ts(poi2)+sr(poi2)-24)./(dtr(poi2)+0.9);        %D values for afternoon (eqn 21c)
poi1=timlc<=12;                                                       %conjugate points before midday
poi2=timlc>12;                                                        %conjugate points after midday
dvalc(poi1)=(src(poi1)+trc(poi1)-timlc(poi1))./dtrc(poi1);            %D values for morning (eqn 21b)
dvalc(poi2)=(timlc(poi2)-tsc(poi2)+src(poi2)-24)./(dtrc(poi2)+0.9);   %D values for afternoon (eqn 21c)

tdiur=1./(1+exp(3.2*dval));                  %diurinal parameter (eqn 21a)
tdiurc=1./(1+exp(3.2*dvalc));                %conjugate diurinal parameter (eqn 21a)

%% To allow the electron temperature profiles to smoothly join at north and south poles
%% need to introduce several functions to modify the profiles. One is a polar correction
%% based on geographihc latitude. Finally, to correctly know which combination of day and night
%% temperatures to use at the poles, use a function based on the solar declination.
latbd=deg2rad(65);                                     %maximum geographic latitude before polar corrections starts coming into effect
latg3p=abs(latg3);                                     %polar copy of absolute of latg3
wtp=(latg3p-latbd)/(pi/2-latbd);                       %linear polar weighting
wtp(wtp<0)=0;                                          %set minimum polar weighting to 0
pol=cos(0.5*pi*wtp);                                   %polar correction term
latg3pc=abs(latg3c);                                   %conjugate polar copy of absolute of latg3
wtpc=(latg3pc-latbd)/(pi/2-latbd);                     %conjugate linear polar weighting
wtpc(wtpc<0)=0;                                        %set minimum conjugate polar weighting to 0
polc=cos(0.5*pi*wtpc);                                 %conjugate polar correction term

TDv=0.97+0.22*((timl-12.5)./min(-4,(sr-12.5))).^2;     %'day' temperature scaling factor (eqn 22a)
TDvc=0.97+0.22*((timlc-12.5)./min(-4,src-12.5)).^2;    %conjugate 'day' temperature  scaling factor (eqn 22a)

dt=timl;                                               %time difference from midnight
dt(dt>12)=dt(dt>12)-24;                                %correct dt value for hours 'before midnight'
dtc=timlc;                                             %time difference from midnight
dtc(dtc>12)=dtc(dtc>12)-24;                            %correct dt value for hours 'before midnight'

TNv=0.83+0.15*exp(-0.2*dt.*(1.4-0.8*slg));             %'night' temperature  scaling factor (eqn 22b)
TNvc=0.83+0.15*exp(-0.2*dtc.*(1.4-0.8*slgc));          %conjugate 'night' temperature  scaling factor (eqn 22b)

%Set upper limits on the allowed scaling factors of 1.2
TDv(TDv>1.2)=1.2;
TDvc(TDvc>1.2)=1.2;
TNv(TNv>1.2)=1.2;
TNvc(TNvc>1.2)=1.2;

%Modify the scaling factors so they are 1 at the poles 
TDv=pol.*TDv+(1-pol);                   %smoothed day factor
TNv=pol.*TNv+(1-pol);                   %smoothed night factor
TDvc=polc.*TDvc+(1-polc);               %smoothed conjugate day factor
TNvc=polc.*TNvc+(1-polc);               %smoothed conjugate night factor

%Modify tdiur so that at the poles it has a value that is an average
%of day and night values. This is based on the current solar declination and
%the pol correction factor.
dmx=0.41015;                            %maximum positive solar declination
decwt=0.5*(sund/dmx+1);                 %declination weighting term
decwt(latg3<0)=1-decwt(latg3<0);        %correction for southern hemisphere terms
tdiur=decwt.*(1-pol)+pol.*tdiur;        %rescale tdiur values
decwtc=0.5*(sundc/dmx+1);               %conjugate declination weighting term
decwtc(latg3c<0)=1-decwtc(latg3c<0);    %correction for southern hemisphere terms
tdiurc=decwtc.*(1-polc)+polc.*tdiurc;   %rescale conjugate tdiur values

%% Combine point of interest with conjugate point to determine
%% required electron temperatures.
xh2=0.5*(dismat-re)./((lsh-1)*re);      %weight for conjugate term
xh1=1-xh2;                              %weight for local term 
 
TDiurn=tdiur.*xh1+tdiurc.*xh2;               %=0, night to 1, day 
TeN=(TNit.*TNv.*xh1)+(TNitc.*TNvc.*xh2);     %eqn 22a (modified)
TeD=(TDay.*TDv.*xh1)+(TDayc.*TDvc.*xh2);     %eqn 22b (modified) 
tmpe=TeN+(TeD-TeN).*TDiurn;                  %estimated electron temperature from JET model (eqn 21a) 
 
TDiurnc=tdiur.*xh2+tdiurc.*xh1;              %=0, night to 1, day 
TeNc=(TNit.*TNv.*xh2)+(TNitc.*TNvc.*xh1);    %eqn 22a (modified)
TeDc=(TDay.*TDv.*xh2)+(TDayc.*TDvc.*xh1);    %eqn 22b (modified) 
tmpec=TeNc+(TeDc-TeNc).*TDiurnc;             %estimated electron temperature from JET model (eqn 21a) 

%Need to use electron temperature correction method to smoothly join the neutral temperatures
tmpno=tmpn;                             %copy of neutral temperature                                       
tmpnoc=tmpnc;                           %copy of conjugate neutral temperature
tmpn=xh1.*tmpno+xh2.*tmpnc;             %corrected neutral temperature
tmpnc=xh2.*tmpno+xh1.*tmpnc;            %corrected conjugate neutral temperature

%% To determine the ion temperature, use the method suggested by Titheridge (1998). 
ht=0.75*tmpn-75;                                                 %ht parameter
amat=0.08+0.18*TDiurn;                                           %a parameter
tmpi=tmpn+(tmpe-tmpn)./(1+amat.*(ht./altmat)+(ht./altmat).^5);   %calculate ion temperature

ht=0.75*tmpnc-75;                                                     %ht parameter
amat=0.08+0.18*TDiurnc;                                               %a parameter
tmpic=tmpnc+(tmpec-tmpnc)./(1+amat.*(ht./altmat)+(ht./altmat).^5);    %calculate conjugate ion temperature

%Need to use electron temperature correction method to smoothly join the ion temperatures
tmpio=tmpi;                             %copy of ion temperature                                       
tmpioc=tmpic;                           %copy of conjugate ion temperature
tmpi=xh1.*tmpio+xh2.*tmpic;             %corrected ion temperature
tmpic=xh2.*tmpio+xh1.*tmpic;            %corrected conjugate ion temperature
tmpi(tmpi>tmpe)=tmpe(tmpi>tmpe);
tmpic(tmpic>tmpec)=tmpec(tmpic>tmpec);

%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);
tdiur=reshape(tdiur,sz);
tdiurc=reshape(tdiurc,sz);
