function [bh,bw,ch,cw] = vmf3o_b_c(doy,lat,lon)

% *************************************************************************
%
% vmf3o_b_c.m
% 
% Computing b and c coefficients of VMF3o
%
% (c) Department of Geodesy and Geoinformation, Vienna University of
%     Technology, 2018
%     created by Janina Boisits
%
% This subroutine determines the coefficients bh, bw, ch, and cw of VMF3o.
% These coefficients are of empirical nature containing a geographical and
% temporal dependence, represented in spherical harmonics. The spherical
% harmonics coefficients are developed to degree and order 12 and are based 
% on a 5°x5° grid containing ray-tracing data from 2001-2010.
%
% -------------------------------------------------------------------------
%
% History:
%
% 2018-12-07: m-file created by Janina Boisits
%
% -------------------------------------------------------------------------
%
% INPUT: 
%     - doy   ...day of year: epoch of observation in days (i.e. January 2,
%                6am -> doy = 1.25)
%                format: [n x 1] vector with n being the number of
%                observations
%     - lat   ...latitude of station in degree
%                format: scalar
%     - lon   ...longitude of station in degree
%                format: scalar
%
% OUTPUT:
%     - bh    ...b coefficient of hydrostatic MF
%                format: [n x 1] vector with n being the number of
%                observations
%     - bw    ...b coefficient of wet MF
%                format: [n x 1] vector with n being the number of
%                observations
%     - ch    ...c coefficient of hydrostatic MF
%                format: [n x 1] vector with n being the number of
%                observations
%     - cw    ...c coefficient of wet MF
%                format: [n x 1] vector with n being the number of
%                observations
%
% *************************************************************************

%% unit vector of station
lon = lon * pi / 180; % rad
lat = lat * pi / 180; % rad
pd  = pi/2 - lat;     % polar distance in rad

x = sin(pd)*cos(lon);
y = sin(pd)*sin(lon);
z = cos(pd);

%% read in coefficients of spherical harmonics
% bh
fid = fopen('anm_bh.txt');
anm_bh = textscan(fid,'%f%f%f%f%f');
fclose(fid);
fid = fopen('bnm_bh.txt');
bnm_bh = textscan(fid,'%f%f%f%f%f');
fclose(fid);
% bw
fid = fopen('anm_bw.txt');
anm_bw = textscan(fid,'%f%f%f%f%f');
fclose(fid);
fid = fopen('bnm_bw.txt');
bnm_bw = textscan(fid,'%f%f%f%f%f');
fclose(fid);
% ch
fid = fopen('anm_ch.txt');
anm_ch = textscan(fid,'%f%f%f%f%f');
fclose(fid);
fid = fopen('bnm_ch.txt');
bnm_ch = textscan(fid,'%f%f%f%f%f');
fclose(fid);
% cw
fid = fopen('anm_cw.txt');
anm_cw = textscan(fid,'%f%f%f%f%f');
fclose(fid);
fid = fopen('bnm_cw.txt');
bnm_cw = textscan(fid,'%f%f%f%f%f');
fclose(fid);

% bh
anm_bh_A0 = anm_bh{1}; anm_bh_A1 = anm_bh{2}; anm_bh_B1 = anm_bh{3}; anm_bh_A2 = anm_bh{4}; anm_bh_B2 = anm_bh{5};
bnm_bh_A0 = bnm_bh{1}; bnm_bh_A1 = bnm_bh{2}; bnm_bh_B1 = bnm_bh{3}; bnm_bh_A2 = bnm_bh{4}; bnm_bh_B2 = bnm_bh{5};
% bw
anm_bw_A0 = anm_bw{1}; anm_bw_A1 = anm_bw{2}; anm_bw_B1 = anm_bw{3}; anm_bw_A2 = anm_bw{4}; anm_bw_B2 = anm_bw{5};
bnm_bw_A0 = bnm_bw{1}; bnm_bw_A1 = bnm_bw{2}; bnm_bw_B1 = bnm_bw{3}; bnm_bw_A2 = bnm_bw{4}; bnm_bw_B2 = bnm_bw{5};
% ch
anm_ch_A0 = anm_ch{1}; anm_ch_A1 = anm_ch{2}; anm_ch_B1 = anm_ch{3}; anm_ch_A2 = anm_ch{4}; anm_ch_B2 = anm_ch{5};
bnm_ch_A0 = bnm_ch{1}; bnm_ch_A1 = bnm_ch{2}; bnm_ch_B1 = bnm_ch{3}; bnm_ch_A2 = bnm_ch{4}; bnm_ch_B2 = bnm_ch{5};
% cw
anm_cw_A0 = anm_cw{1}; anm_cw_A1 = anm_cw{2}; anm_cw_B1 = anm_cw{3}; anm_cw_A2 = anm_cw{4}; anm_cw_B2 = anm_cw{5};
bnm_cw_A0 = bnm_cw{1}; bnm_cw_A1 = bnm_cw{2}; bnm_cw_B1 = bnm_cw{3}; bnm_cw_A2 = bnm_cw{4}; bnm_cw_B2 = bnm_cw{5};

%% compute Legendre polynomials
% degree/order of SH
nmax = 12;

% Legendre polynomials
V(1,1) = 1;
W(1,1) = 0;
V(2,1) = z * V(1,1);
W(2,1) = 0;

for n = 2:nmax
    V(n+1,1) = ((2*n-1) * z * V(n,1) - (n-1) * V(n-1,1)) / n;
    W(n+1,1) = 0;
end

for m = 1:nmax
    V(m+1,m+1) = (2*m-1) * (x*V(m,m) - y*W(m,m));
    W(m+1,m+1) = (2*m-1) * (x*W(m,m) + y*V(m,m));
    if (m < nmax) 
        V(m+2,m+1) = (2*m+1) * z * V(m+1,m+1);
        W(m+2,m+1) = (2*m+1) * z * W(m+1,m+1);
    end
    for n = m+2:nmax
        V(n+1,m+1) = ((2*n-1)*z*V(n,m+1) - (n+m-1)*V(n-1,m+1)) / (n-m);
        W(n+1,m+1) = ((2*n-1)*z*W(n,m+1) - (n+m-1)*W(n-1,m+1)) / (n-m);
    end
end

%% evaluate SH
bh_A0 = 0; bh_A1 = 0; bh_B1 = 0; bh_A2 = 0; bh_B2 = 0;
bw_A0 = 0; bw_A1 = 0; bw_B1 = 0; bw_A2 = 0; bw_B2 = 0;
ch_A0 = 0; ch_A1 = 0; ch_B1 = 0; ch_A2 = 0; ch_B2 = 0;
cw_A0 = 0; cw_A1 = 0; cw_B1 = 0; cw_A2 = 0; cw_B2 = 0;
i = 0;

for n = 0:nmax
    for m = 0:n
        
        i = i+1;
        
        bh_A0 = bh_A0 + (anm_bh_A0(i)*V(n+1,m+1) + bnm_bh_A0(i)*W(n+1,m+1));
        bh_A1 = bh_A1 + (anm_bh_A1(i)*V(n+1,m+1) + bnm_bh_A1(i)*W(n+1,m+1));
        bh_B1 = bh_B1 + (anm_bh_B1(i)*V(n+1,m+1) + bnm_bh_B1(i)*W(n+1,m+1));
        bh_A2 = bh_A2 + (anm_bh_A2(i)*V(n+1,m+1) + bnm_bh_A2(i)*W(n+1,m+1));
        bh_B2 = bh_B2 + (anm_bh_B2(i)*V(n+1,m+1) + bnm_bh_B2(i)*W(n+1,m+1));
        
        bw_A0 = bw_A0 + (anm_bw_A0(i)*V(n+1,m+1) + bnm_bw_A0(i)*W(n+1,m+1));
        bw_A1 = bw_A1 + (anm_bw_A1(i)*V(n+1,m+1) + bnm_bw_A1(i)*W(n+1,m+1));
        bw_B1 = bw_B1 + (anm_bw_B1(i)*V(n+1,m+1) + bnm_bw_B1(i)*W(n+1,m+1));
        bw_A2 = bw_A2 + (anm_bw_A2(i)*V(n+1,m+1) + bnm_bw_A2(i)*W(n+1,m+1));
        bw_B2 = bw_B2 + (anm_bw_B2(i)*V(n+1,m+1) + bnm_bw_B2(i)*W(n+1,m+1));
        
        ch_A0 = ch_A0 + (anm_ch_A0(i)*V(n+1,m+1) + bnm_ch_A0(i)*W(n+1,m+1));
        ch_A1 = ch_A1 + (anm_ch_A1(i)*V(n+1,m+1) + bnm_ch_A1(i)*W(n+1,m+1));
        ch_B1 = ch_B1 + (anm_ch_B1(i)*V(n+1,m+1) + bnm_ch_B1(i)*W(n+1,m+1));
        ch_A2 = ch_A2 + (anm_ch_A2(i)*V(n+1,m+1) + bnm_ch_A2(i)*W(n+1,m+1));
        ch_B2 = ch_B2 + (anm_ch_B2(i)*V(n+1,m+1) + bnm_ch_B2(i)*W(n+1,m+1));
        
        cw_A0 = cw_A0 + (anm_cw_A0(i)*V(n+1,m+1) + bnm_cw_A0(i)*W(n+1,m+1));
        cw_A1 = cw_A1 + (anm_cw_A1(i)*V(n+1,m+1) + bnm_cw_A1(i)*W(n+1,m+1));
        cw_B1 = cw_B1 + (anm_cw_B1(i)*V(n+1,m+1) + bnm_cw_B1(i)*W(n+1,m+1));
        cw_A2 = cw_A2 + (anm_cw_A2(i)*V(n+1,m+1) + bnm_cw_A2(i)*W(n+1,m+1));
        cw_B2 = cw_B2 + (anm_cw_B2(i)*V(n+1,m+1) + bnm_cw_B2(i)*W(n+1,m+1));
        
    end
end

bh = bh_A0 + bh_A1*cos(doy./365.25*2*pi) + bh_B1*sin(doy./365.25*2*pi) + bh_A2*cos(doy./365.25*4*pi) + bh_B2*sin(doy./365.25*4*pi);
bw = bw_A0 + bw_A1*cos(doy./365.25*2*pi) + bw_B1*sin(doy./365.25*2*pi) + bw_A2*cos(doy./365.25*4*pi) + bw_B2*sin(doy./365.25*4*pi);
ch = ch_A0 + ch_A1*cos(doy./365.25*2*pi) + ch_B1*sin(doy./365.25*2*pi) + ch_A2*cos(doy./365.25*4*pi) + ch_B2*sin(doy./365.25*4*pi);
cw = cw_A0 + cw_A1*cos(doy./365.25*2*pi) + cw_B1*sin(doy./365.25*2*pi) + cw_A2*cos(doy./365.25*4*pi) + cw_B2*sin(doy./365.25*4*pi);

end
