%% temppaw.m is my latest implementation of J.E.Titheridge temperature 
%% model as set out in Feb 1998 JGR. Based on tempjet3.m, it is different
%% in that it uses my fitted values for base temperature and gradient that
%% does not include EXOS-D, and 400km Alouette data. The base temperature 
%% and gradient term results are fitted with polynomial functions, not the
%% rationalized polynominals that Titheridge used (eqn 19). The polynominal
%% fit is obtained from tempfit7.m.
%% Corrections are included for seasons night and for the effects of 
%% averaging around midday and midnight of the orignal data sets is included.
%% The neutral temperatures are calculated by the MSIS-90 directly. Note that
%% for points above 500km, that the 500km value is returned.
%% Due to disagreement between the calculated sunrise/sunset times of
%% John T.'s model (eqn 20a) and the IRI95, a modified version of 20a is 
%% used.
%% The model is valid only above 350km, and values at high latitudes (75+ 
%% degrees) should be considered approximate only. Between 350km and 400km,
%% Titheridge correction (see below eqn 13b) is not used. Reset to zero any
%% term in the ^(2/7) bracket that is negative instead. 
%%
%% NOTE - (4/7/99) I have decided to completely redo J.T. correction values
%%        at high latitudes. They seem overly complicated and are really
%%        designed for the case where magnetic and geographic co-ords are the
%%        same, which is not what I wanted. (6/7/99) Have managed to replace
%%        them with some equally complicated functions, but at least I
%%        understand exactly what each of them is doing!
%%
%% NOTE - (13/3/00) I have realised that J.T. approach to the effect of F107
%%        is not correct. A variety of observations and the FLIP clearly show
%%        that the electron temperature drops as the electron density increases
%%        from solar minimum to maximum. While the energy input has increased
%%        (F10.7 up by a factor of ~3 and heating proportional to density) the
%%        cooling rate has increased more, since it goes as density squared.
%%        Therefore I have run the FLIP to estimate this variation as a function 
%%        of F107 (see fliptpf.m to fliptpf6.m).
%%
%% NOTE - (13/3/00) Add in seasonal variation (both day and night) in place of 
%%        of correction based on electron density at 350km. See tempfit8.m for
%%        determination of seasonal corrections polynominal co-efficents.
%%
%% NOTE - (13/3/00) Decided to not include the averaging corrections for day and
%%        night values. 
%%
%% NOTE - (16/5/00) Use more rigid method to determine function describing variation
%%        of tr and ts with local sunrise. Used to replace eqn 20a (see tempss.m).
%%
%% NOTE - (9/6/00) Simplified the polar region weighting functions, which no longer
%%        use the slong parameter, which was causing a few problems with very large
%%        values near the pole.
%%
%% NOTE - (29/10/2001) Modify to use the 7 value Ap vector for MSIS.
%%
%% NOTE - (5/3/2002) Complete overall of the code! Now use T0 and G0 obtained from 
%%        fits to both electron temperature models in the IRI - the original Brace
%%        and Theis model used previoulsy and the newer Truhlik et al. empirical 
%%        model. The later is now used to predict the variation with F10.7, which
%%        seperately calculates changes in both the gradient and base temperature.
%%        Seasonal corrections are now also seperately calculated for gradient and
%%        and base temperature.
%%
%% NOTE - (22/3/2002) Ok, have decided to go back to FLIP derived corrections for 
%%        variations in F10.7. This was because the Truhlik et al. ones result in
%%        an overall decrease in the heat flux with increasing F10.7 - not good!  
%%        
%% NOTE - (3/4/2003) Based on some comments from John Titheridge on my temperature
%%        paper I have decided to use his correction for Te at low altitudes. I now
%%        realise that this gives better results at low altitudes were (espically 
%%        during the day) Te will not equal Tn until ~ 100 - 200 km.
%%
%% April 15 2003 - change the Go correction so that it is fractional change
%% rather than a absolute change.
%%
%% Author: Phillip A. Webb
%% Last updated: October 21 2003 
%%
    
function [tmpi,tmpe,tmpn,tmpic,tmpec,tmpnc,tdiur,tdiurc,amat,amatc,t0out,t0outc,g0out,g0outc]=temppaw(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).
% tdiur = diurinal parameter
% tdiurc = conjugate diurinal parameter
% amat = Titheridge 'a' parameter
% amatc = conjugate Titheridge 'a' parameter
% t0out = estimate for final T0 value (K)
% t0outc = estimate for final conjugate T0 value (K)
% g0out = estimate for final G0 value (K)
% g0outc = estimate for final conjugate G0 value (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(:);

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

[kp,kpmin,apday,f107,f107p,f107a,apmat]=solarval(timmat,dno,year); %calculate required solar parameters

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;

%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;
dhr(poid)=dhr(poid)./(1-0.01*dhr(poid)*r0^2);

TDay=t0d.*(1+bh.*(g0d./t0d).*dhr).^(2/7);              %average electron day temperature (eqn 13)
TDayc=t0dc.*(1+bh.*(g0dc./t0dc).*dhr).^(2/7);          %average conjugate electron day temperature (eqn 13)
TNit=t0n.*(1+bh.*(g0n./t0n).*dhr).^(2/7);              %average electron night temperature (eqn 13)
TNitc=t0nc.*(1+bh.*(g0nc./t0nc).*dhr).^(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) 

%Determine estimates of final T0 and G0 values for output
t0out=t0d+(t0d-t0n).*TDiurn;
t0outc=t0dc+(t0dc-t0nc).*TDiurnc;
g0out=g0d+(g0d-g0n).*TDiurn;
g0outc=g0dc+(g0dc-g0nc).*TDiurnc;

%% To determine the ion temperature, use the method suggested by Titheridge (1998). 
%% Calculate the neutral temperature using the MSIS model. To determine the 
%% neutral temperature, calculate the temperature using the geographic locations
%% of the two bases of the field line. For points above 750km, calculate the 
%% temperature at 750km (since constant by this point). Then use the weightings used to
%% mix the conjugate electron temperatures to mix the neutral temperatures.
%% 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,hent,n2nt,tmpnt]=modmsis1([long3 long3c],[latg3 latg3c],dism,timmat*svt,day*svt,month*svt,year*svt,f107a*svt,f107p*svt,[apmat apmat],1);
tmpn=xh1.*tmpnt(:,1)+xh2.*tmpnt(:,2);                            %corrected neutral temperature
tmpnc=xh2.*tmpnt(:,1)+xh1.*tmpnt(:,2);                           %corrected conjugate neutral temperature

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
amatc=0.08+0.18*TDiurnc;                                                   %a parameter
tmpic=tmpnc+(tmpec-tmpnc)./(1+amatc.*(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.*tmpioc;            %corrected ion temperature
tmpic=xh2.*tmpio+xh1.*tmpioc;           %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);
tmpn=reshape(tmpn,sz);
tmpnc=reshape(tmpnc,sz);
tdiur=reshape(tdiur,sz);
tdiurc=reshape(tdiurc,sz);
