%% temppawsea.m calculates the requred seasonal variation in the base 
%% temperature and gradients required by the Titheridge model. These 
%% are variations away from the mean values produced from temppawt0g0.m.
%% The polynominal values given here are derived by pwtpfit6.m and were 
%% produced using a base height of 400km.
%% Note - points with sl<0.05 are reset to sL=0.05 before the variation
%%        is calculated. However, the original sl value is used when 
%%        calculating the required correction near the equator.
%%
%%
%% April 15 2003 - change the Go correction so that it is fractional change
%% rather than a absolute change.
%%
%% Author: Phillip A. Webb
%% Last updated: June 9 2004                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
%%

function [dt0d,dt0n,dg0d,dg0n,dt0dc,dt0nc,dg0dc,dg0nc]=temppawsea(sl,dno)

% sl = sL values of the points of interest
% dno = day number of points of interest
% dt0d = northern hemisphere seasonal variation in day base temperature (K)
% dt0n = northern hemisphere seasonal variation night base temperature (K)
% dg0d = northern hemisphere seasonal variation day base gradient (K/km)
% dg0n = northern hemisphere seasonal variation night base gradient (K/km)
% dt0dc = southern hemisphere seasonal variation in day base temperature (K)
% dt0nc = southern hemisphere seasonal variation night base temperature (K)
% dg0dc = southern hemisphere seasonal variation day base gradient (K/km)
% dg0nc = southern hemisphere seasonal variation night base gradient (K/km)

%Polynominal coefficents from pwtpfit6.m
load pwsea 

%Calculate required values
tde=polyval(dt0dp1,sl);                      %day equinox temperature correction
tds=polyval(dt0dp2,sl);                      %day summer temperature correction
tdw=polyval(dt0dp3,sl);                      %day winter temperature correction
tne=polyval(dt0np1,sl);                      %night equinox temperature correction
tns=polyval(dt0np2,sl);                      %night summer temperature correction
tnw=polyval(dt0np3,sl);                      %night winter temperature correction

gde=polyval(dg0dp1,sl);                      %day equinox gradient correction
gds=polyval(dg0dp2,sl);                      %day summer gradient correction
gdw=polyval(dg0dp3,sl);                      %day winter gradient correction
gne=polyval(dg0np1,sl);                      %night equinox gradient correction
gns=polyval(dg0np2,sl);                      %night summer gradient correction
gnw=polyval(dg0np3,sl);                      %night winter gradient correction

yev=1.2*sin((pi/182.5)*(dno+11.25)).^2;      %variation in equiox value
yev2=1.2-yev;                                

dt0d=tde.*yev;                     %set season day temperature correction to day equinox correction
dt0dc=tde.*yev;                    %set season conjugate temperature day correction to day equinox correction
dt0n=tne.*yev;                     %set season night temperature correction to night equinox correction
dt0nc=tne.*yev;                    %set season conjugate night temperature correction to night equinox correction

dg0d=gde.*yev;                     %set season day gradient correction to day equinox correction
dg0dc=gde.*yev;                    %set season conjugate gradient day correction to day equinox correction
dg0n=gne.*yev;                     %set season night gradient correction to night equinox correction
dg0nc=gne.*yev;                    %set season conjugate gradient night correction to night equinox correction

%Calculate the required day and night corrections for both hemispheres
pois=dno>80 & dno<=264;

dt0d(pois)=dt0d(pois)+(tds(pois).*yev2(pois));
dt0d(~pois)=dt0d(~pois)+(tdw(~pois).*yev2(~pois));
dt0dc(pois)=dt0dc(pois)+(tdw(pois).*yev2(pois));
dt0dc(~pois)=dt0dc(~pois)+(tds(~pois).*yev2(~pois));

dt0n(pois)=dt0n(pois)+(tns(pois).*yev2(pois));
dt0n(~pois)=dt0n(~pois)+(tnw(~pois).*yev2(~pois));
dt0nc(pois)=dt0nc(pois)+(tnw(pois).*yev2(pois));
dt0nc(~pois)=dt0nc(~pois)+(tns(~pois).*yev2(~pois));

dg0d(pois)=dg0d(pois)+(gds(pois).*yev2(pois));
dg0d(~pois)=dg0d(~pois)+(gdw(~pois).*yev2(~pois));
dg0dc(pois)=dg0dc(pois)+(gdw(pois).*yev2(pois));
dg0dc(~pois)=dg0dc(~pois)+(gds(~pois).*yev2(~pois));

dg0n(pois)=dg0n(pois)+(gns(pois).*yev2(pois));
dg0n(~pois)=dg0n(~pois)+(gnw(~pois).*yev2(~pois));
dg0nc(pois)=dg0nc(pois)+(gnw(pois).*yev2(pois));
dg0nc(~pois)=dg0nc(~pois)+(gns(~pois).*yev2(~pois));

%Add in a correction so the that day and night values agree at sL=1 in both
%hemispheres (this is not currently done by the fit from pwtpfit6.m). Use a 
%sine weighting for points above sL>0.9.
poic=sl>0.9;
if any(any(poic))

 dt0d0=dt0d;
 dt0n0=dt0n;
 dg0d0=dg0d;
 dg0n0=dg0n;

 dt0dc0=dt0dc;
 dt0nc0=dt0nc;
 dg0dc0=dg0dc;
 dg0nc0=dg0nc;

 wt=sin(((sl(poic)-0.9)/0.2).*(pi/2)).^2;

 dt0d(poic)=(1-wt).*dt0d0(poic)+wt.*dt0n0(poic);
 dt0n(poic)=(1-wt).*dt0n0(poic)+wt.*dt0d0(poic);
 dg0d(poic)=(1-wt).*dg0d0(poic)+wt.*dg0n0(poic);
 dg0n(poic)=(1-wt).*dg0n0(poic)+wt.*dg0d0(poic);

 dt0dc(poic)=(1-wt).*dt0dc0(poic)+wt.*dt0nc0(poic);
 dt0nc(poic)=(1-wt).*dt0nc0(poic)+wt.*dt0dc0(poic);
 dg0dc(poic)=(1-wt).*dg0dc0(poic)+wt.*dg0nc0(poic);
 dg0nc(poic)=(1-wt).*dg0nc0(poic)+wt.*dg0dc0(poic);

end

%Add in a correction so the that northern and southern hemisphere values
%agree at the equator (sL=0) (this is not currently done by the fit from 
%pwtpfit6.m). Use a sine weighting for points below sL<0.2.
poic=sl<0.2;
if any(any(poic))

 dt0d0=dt0d;
 dt0n0=dt0n;
 dg0d0=dg0d;
 dg0n0=dg0n;

 dt0dc0=dt0dc;
 dt0nc0=dt0nc;
 dg0dc0=dg0dc;
 dg0nc0=dg0nc;

 wt=sin(((1-sl(poic))/0.4).*(pi/2)).^2;

 dt0d(poic)=(1-wt).*dt0d0(poic)+wt.*dt0dc0(poic);
 dt0n(poic)=(1-wt).*dt0n0(poic)+wt.*dt0nc0(poic);
 dt0dc(poic)=(1-wt).*dt0dc0(poic)+wt.*dt0d0(poic);
 dt0nc(poic)=(1-wt).*dt0nc0(poic)+wt.*dt0n0(poic);

 dg0d(poic)=(1-wt).*dg0d0(poic)+wt.*dg0dc0(poic);
 dg0n(poic)=(1-wt).*dg0n0(poic)+wt.*dg0nc0(poic);
 dg0dc(poic)=(1-wt).*dg0dc0(poic)+wt.*dg0d0(poic);
 dg0nc(poic)=(1-wt).*dg0nc0(poic)+wt.*dg0n0(poic);

end