function [tec,dtde]=tecmodel(time,posp,enm,model,pos0)
% TECMODEL : TEC model
%
% Compute vertical TEC (Total Electron Content) by the global or local tec model.
% The tec models comply with as follows.
%
% local  : sum(sum(enm*((lat-lat0)/s)^n*((lon-lon0)/s)^m)) (n,m=0,1,...,nmax)
% global : sum(sum(enm*Ynm(lat,UT+lon-pi))) (n=0,1,...,nmax,m=0,1,...,n)
%
%   (lat,lon=latitude/longitude, lat0,lon0=origin of development, s=scale,
%    Ynm=(n,m) spheric harmonic function)
%
%          Copyright (C) 2006 by T.TAKASU, All rights reserved.
%
% argin  : time  = time relative to 0:00 (sec)
%          posp  = latitude/longitude [lat,lon] (rad)
%          enm   = serialized tec model coefficients
%              local  : [e00,e01,e02,...,e10,e11,e12,...,e20,e21,e22,...]
%              global : [e00,e10,e11c,e11s,e20,e21c,e21s,e22c,e22s,...]
%                       (length(enm)=(nmax+1)*(nmax+1))
%          model = tec model (0:local,1:global)
%          pos0  = origin of local tec model [lat0,lon0] (rad)
%
% argout : tec   = vertical TEC (TECU)
%          dtde  = partial drivatives d(tec)/d(enm)
%
% version: $Revision: 3 $ $Date: 06/02/25 21:35 $
% history: 2006/02/23 1.1 new
%          2006/02/25 1.2 fix bug of spheric harmonic function Yn0

nmax=round(sqrt(length(enm)))-1;
if ~model, [tec,dtde]=teclocal(posp,enm,nmax,pos0);
else       [tec,dtde]=tecglobal(time,posp,enm,nmax); end

% local tec model --------------------------------------------------------------
function [tec,dtde]=teclocal(posp,enm,nmax,pos0)
dtde=zeros(1,(nmax+1)*(nmax+1)); i=1; s=10*pi/180; % scale
for n=0:nmax
    for m=0:nmax
        dtde(i)=((posp(1)-pos0(1))/s)^n*((posp(2)-pos0(2))/s)^m; i=i+1;
    end
end
tec=enm*dtde';

% global tec model -------------------------------------------------------------
function [tec,dtde]=tecglobal(time,posp,enm,nmax)
dtde=sphfun(posp(1),pi*time/43200+posp(2)-pi,nmax);
tec=enm*dtde';

% spheric harmonic function : y=[Y00,Y10,Y11c,Y11s,Y20,Y21c,Y21s,...] ----------
function y=sphfun(lat,lon,nmax)
y=zeros(1,(nmax+1)*(nmax+1));
p=nlegendre(sin(lat),nmax); i=1;
for n=0:nmax
    y(i)=p(n+1,1); i=i+1;
    for m=1:n, y(i:i+1)=p(n+1,m+1)*[cos(m*lon),sin(m*lon)]; i=i+2; end
end

% normalized legendre function : p(n+1,m+1)=Pnm --------------------------------
function p=nlegendre(x,nmax)
p=zeros(nmax+1,nmax+1); c=ones(size(p)); p(1:2,1:2)=[1,0;x,sqrt(1-x*x)]; 
for n=2:nmax
    p(n+1,n+1)=(2*n-1)*sqrt(1-x*x)*p(n,n);
    for m=0:n-1, p(n+1,m+1)=((2*n-1)*x*p(n,m+1)-(n+m-1)*p(n-1,m+1))/(n-m); end
end
for n=1:nmax
    c(n+1,1)=sqrt(2*n+1);
    for m=1:n, c(n+1,m+1)=sqrt(2*(2*n+1)*factorial(n-m)/factorial(n+m)); end
end
p=c.*p;
