%% temppaw35.m is designed to be called by temppaw34.m. See it for input
%% and output description.
%%
%% Author: Phillip A. Webb
%% Last updated: April 3 2003 
%%

function [tmpe,tmpi,tmpec,tmpic,tdiur,tdiurc]=temppaw35(t0d,t0n,g0d,g0n,t0dc,t0nc,g0dc,g0nc,bh,dhr,lsh,long3,latg3,long3c,latg3c,dismat,latmatm,tmpn,tmpnc,timmat,dno,year)

[day,month]=daynum2(dno,year);                         %day and month
re=6371;                                               %Earth radi (km)
h04=400;                                               %2nd base height        
altmat=dismat-re;                                      %altitude of points of interest
svy=ones(size(dismat,1),1);
svx=ones(1,size(dismat,2));

latm3=latmatm(:,1);                                    %magnetic latitude at 350km
latm3d=rad2deg(abs(latm3));                            %350km altitude point in degrees

%Use the Titheridge correction at low altitudes. For points with
%dhr<0 reset equal to dhr/(1-0.01*dhr*r0^2).
poid=dhr<0;
r0=1+h04/re;
dhr(poid)=dhr(poid)./(1-0.01*dhr(poid)*r0^2);
TDay=(t0d*svx).*(1+bh.*((g0d./t0d)*svx).*dhr).^(2/7);                 %average electron day temperature (eqn 13)
TDayc=(t0dc*svx).*(1+bh.*((g0dc./t0dc)*svx).*dhr).^(2/7);             %average conjugate electron day temperature (eqn 13)
TNit=(t0n*svx).*(1+bh.*((g0n./t0n)*svx).*dhr).^(2/7);                 %average electron night temperature (eqn 13)
TNitc=(t0nc*svx).*(1+bh.*((g0nc./t0nc)*svx).*dhr).^(2/7);             %average conjugate electron night temperature (eqn13)

%Determine time of local ground sunrise for point below field line at 350km altitude
[sr,ss,slong,sund]=sunrs(latg3,day*svy,month*svy,year*svy);           %local ground sunrise          
[src,ssc,slongc,sundc]=sunrs(latg3c,day*svy,month*svy,year*svy);      %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(length(timl),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*altmat./(re*(lsh-1)*svx);       %weight for conjugate term
xh1=1-xh2;                              %weight for local term 

TDiurn=(tdiur*svx).*xh1+(tdiurc*svx).*xh2;                  %=0, night to 1, day 
TeN=(TNit.*(TNv*svx).*xh1)+(TNitc.*(TNvc*svx).*xh2);        %eqn 22a (modified)
TeD=(TDay.*(TDv*svx).*xh1)+(TDayc.*(TDvc*svx).*xh2);        %eqn 22b (modified) 
tmpe=TeN+(TeD-TeN).*TDiurn;                                 %estimated electron temperature from JET model (eqn 21a) 
 
TDiurnc=(tdiur*svx).*xh2+(tdiurc*svx).*xh1;                 %=0, night to 1, day 
TeNc=(TNit.*(TNv*svx).*xh2)+(TNitc.*(TNvc*svx).*xh1);       %eqn 22a (modified)
TeDc=(TDay.*(TDv*svx).*xh2)+(TDayc.*(TDvc*svx).*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);
