function sg2 = sg2_calculation(j_ut, lat, lon, h, varargin)
% sg2 = SG2_calculation(j_ut, lat, lon, [h, 'param_1','val_1',...])
% j_ut : vector of julian dates UT            (decimal day)
% lat  : vector of latitude' geopoints        (degree)
% lon  : vector of longitude' geopoints       (degree)
% [h]  : vector of elevations' geopoint       (m)              [0 m]
% ['param_k', 'val_k'] : optional parameters
% - 'P'       : Pressure                      (mbar)           [1010 mbar]
% - 'T'       : Temperature                   (°C)             [10 °C]
% - 'Ellps'   : Reference ellips
%               {'WGS84','NTF','AA','SPA','RGF83','SPHERE'}    ['SPA']
% - 'Delta_t' : TT-UT1                        (s)              [Polynomial approximation +/- 0.5 s]
% - 'alpha'   : surface azimuth               (degree)         [0°]
% - 'beta'    : surface slope                 (degree)         [0°]

%% Constants

deg2rad = pi/180;
ua = 149597870691;

%% Managment of input optional parameters
param_opt.TZ = 0;
param_opt.P = 1010;
param_opt.T = 10;
param_opt.ellps = 'SPA';
param_opt.delta_t = [];
param_opt.alpha = 0;
param_opt.beta = 0;
param_opt.atmos_refract = 0.5667;
param_opt.sun_radius = 0.2667;
param_opt.gamma_S_sunrise = -0.8333;
param_opt.corr_refract = 1;
param_opt.corr_refract_method = 'sae';

if (nargin == 0)
    sg2 = param_opt;
    return
end

ni = length(varargin);

if (ni == 1)
   param_opt = varargin{1};
   fn = fieldnames(param_opt);
   for k = 1:length(fn),
       varargin{2*k-1} = fn{k};
       varargin{2*k} = param_opt.(fn{k});
   end
   ni = 2*length(fn);
end

if (ni >= 2) & (mod(ni,2)==0)
    for ki = 1:2:ni,
        name_param = lower(varargin{ki});
        val_param = varargin{ki+1};
        switch(name_param)
            case 'p'
                param_opt.P = val_param;
            case 't'
                param_opt.T = val_param;
            case 'ellps'
                param_opt.ellps = val_param;
            case 'delta_t'
                param_opt.delta_t = val_param;
            case 'alpha'
                param_opt.alpha = val_param;
            case 'beta'
                param_opt.beta = val_param;
            case 'atmos_refract'
                param_opt.atmos_refract = val_param;
            case 'sun_radius'
                param_opt.sun_radius = val_param;
            case 'gamma_S_sunrise'
                param_opt.gamma_S_sunrise = val_param;
            case 'corr_refract'
                param_opt.corr_refract = (val_param==1);
            case 'corr_refract_method'
                param_opt.corr_refract_method = val_param;
        end
    end
end
sg2.info0 = '** Optional parameters';
sg2.param_opt = param_opt;

P = param_opt.P;
T = param_opt.T;
ellps = param_opt.ellps;
delta_t = param_opt.delta_t;
alpha = param_opt.alpha*deg2rad;
beta = param_opt.beta*deg2rad;
atmos_refract = param_opt.atmos_refract*deg2rad;
sun_radius = param_opt.sun_radius*deg2rad;
gamma_S_sunrise = param_opt.gamma_S_sunrise*deg2rad;
corr_refract = param_opt.corr_refract;
corr_refract_method = param_opt.corr_refract_method;

%% Initialisation

approx_R_1.j0 = 2444239.5; % ymd_to_jd([1980 1 1 0]);
approx_R_1.N = 1;
approx_R_1.f0 = 1/365.254902;
approx_R_1.a0 = 2*pi*approx_R_1.f0;
approx_R_1.ro = 0.016704;
approx_R_1.phi = -3.091159;
approx_R_1.a = 0.;
approx_R_1.b = 1.000140;

approx_L_10.j0 = 2444239.5; % ymd_to_jd([1980 1 1 0]);
approx_L_10.N = 10;
approx_L_10.f0 = 1./[365.261278 182.632412 29.530634 399.529850 291.956812 ...
    583.598201 4652.629372 1450.236684 199.459709 365.355291];
approx_L_10.a0 = 2*pi*approx_L_10.f0;
approx_L_10.ro =  [3.401508e-002 3.486440e-004 3.136227e-005 3.578979e-005 ...
    2.676185e-005 2.333925e-005 1.221214e-005 1.217941e-005 1.343914e-005 8.499475e-004];
approx_L_10.phi = [1.600780 1.662976 -1.195905 -1.042052 2.012613 -2.867714 ...
    1.225038 -0.828601 -3.108253 -2.353709];
approx_L_10.a = 1/58.130101;
approx_L_10.b = 1.742145;

approx_Dpsi_1.j0 = 2444239.5; % ymd_to_jd([1980 1 1 0]);
approx_Dpsi_1.N = 1;
approx_Dpsi_1.f0 = 1/6791.164405;
approx_Dpsi_1.a0 = 2*pi*approx_Dpsi_1.f0;
approx_Dpsi_1.ro = 8.329092e-5;
approx_Dpsi_1.phi = -2.052757;
approx_Dpsi_1.a = 0;
approx_Dpsi_1.b = 0;

approx_Dtau_cst = -20.4898/3600*deg2rad;

approx_Epsilon_1.j0 = 2444239.5; % ymd_to_jd([1980 1 1 0]);
approx_Epsilon_1.N = 1;
approx_Epsilon_1.f0 = 1/6791.164405;
approx_Epsilon_1.a0 = 2*pi*approx_Epsilon_1.f0;
approx_Epsilon_1.ro = 4.456183e-5;
approx_Epsilon_1.phi = 2.660352;
approx_Epsilon_1.a = -6.216374e-9;
approx_Epsilon_1.b = 4.091383e-1;

approx_nu0_0.j0 = 2444239.5; % ymd_to_jd([1980 1 1 0]);
approx_nu0_0.N = 0;
approx_nu0_0.f0 = [];
approx_nu0_0.a0 = [];
approx_nu0_0.ro = [];
approx_nu0_0.phi = [];
approx_nu0_0.a = 6.300388099;
approx_nu0_0.b = 1.742079140;

approx_M_0.j0 = 2444239.5; % ymd_to_jd([1980 1 1 0]);
approx_M_0.N = 0;
approx_M_0.f0 = [];
approx_M_0.a0 = [];
approx_M_0.ro = [];
approx_M_0.phi = [];
approx_M_0.a = 1/58.130099904;
approx_M_0.b = -1.399410798;

approx_DeltaT_MSc.nP = 3;
approx_DeltaT_MSc.y = [1961 1986 ; 1986 2005 ; 2005 2050];
approx_DeltaT_MSc.P(1).y = 1975;
approx_DeltaT_MSc.P(1).a = fliplr([45.45 1.067 -1/260 -1/718]);
approx_DeltaT_MSc.P(2).y = 2000;
approx_DeltaT_MSc.P(2).a = fliplr([63.86 0.3345 -0.060374 ...
    0.0017275 0.000651814 0.00002373599]);
approx_DeltaT_MSc.P(3).y = 2000;
approx_DeltaT_MSc.P(3).a = fliplr([62.8938127 0.32100612 0.005576068]);

%% Date manipulations

[nl,nc] = size(j_ut);
j_ut = j_ut(:);
nj = length(j_ut);

ymd_ut = sg2_jd_to_ymd(j_ut);
year_mm_dec = ymd_ut(:,1) + (ymd_ut(:,2)-0.5)/12;

sg2.info1 = '** Dates and hours';
sg2.nj = nj;
sg2.j_ut = j_ut;
sg2.ymd_ut = ymd_ut;
sg2.year_mm_dec = year_mm_dec;

%% Computation of delta_t = TT - UT1 => j_tt

if (isempty(delta_t))
    delta_t = sg2_calc_delta_t(year_mm_dec, approx_DeltaT_MSc);
end
delta_t_j = delta_t/86400;
j_tt = j_ut + delta_t_j;

sg2.delta_t = delta_t;
sg2.delta_t_j = delta_t_j;
sg2.j_tt = j_tt;

%% Computation of heliocentric coordinates

R = sg2_calc_approxS(j_tt, approx_R_1);
L = sg2_calc_approxS(j_tt, approx_L_10);

sg2.info2 = '** Heliocentric coordinates';
sg2.R = sg2_setrange(R);
sg2.L = sg2_setrange(L);

%% Computation of geocentric coordinates

Delta_psi = sg2_setrange(sg2_calc_approxS(j_tt, approx_Dpsi_1),0);
Delta_tau = approx_Dtau_cst;
Theta = L + pi;
Theta_a = Theta + Delta_psi + Delta_tau;
epsilon = sg2_calc_approxS(j_tt, approx_Epsilon_1);

cos_epsilon = cos(epsilon);
sin_Theta_a = sin(Theta_a);
delta_g = asin(sin_Theta_a.*sin(epsilon));
r_alpha_g = atan2(sin_Theta_a.*cos_epsilon, cos(Theta_a));

sg2.info3 = '** Geocentric coordinates';
sg2.Delta_tau = Delta_tau;
sg2.Delta_psi = sg2_setrange(Delta_psi,0);
sg2.Theta = sg2_setrange(Theta);
sg2.Theta_a = sg2_setrange(Theta_a);
sg2.epsilon = sg2_setrange(epsilon);
sg2.delta_g = sg2_setrange(delta_g, 0);  % [-pi, pi]
sg2.r_alpha_g = sg2_setrange(r_alpha_g, 1); % [0, 2pi]

%% Temps sidéraux et EOT
nu0 = sg2_calc_approxS(j_ut, approx_nu0_0);
Delta_psi_cos_epsilon = Delta_psi.*cos_epsilon;
nu = nu0 + Delta_psi_cos_epsilon;

M = sg2_calc_approxS(j_tt, approx_M_0);
EOT = M - 0.0001 - r_alpha_g + Delta_psi_cos_epsilon;

sg2.info4 = '** Sideral times and EOT';
sg2.nu0 = sg2_setrange(nu0);
sg2.nu = sg2_setrange(nu);
sg2.M = sg2_setrange(M);
sg2.EOT = sg2_setrange(EOT);

%% Local observers

if (isempty(lat) | isempty(lon))
    return;
end
if (isempty(h))
    h = 0;
end

phi = lat(:)'*deg2rad;
lambda = lon(:)'*deg2rad;
def_ellps = sg2_ref_ellps(ellps);

np = length(phi);
if ((length(h) == 1) && (np > 1))
    h = repmat(h,1,np);
else
    h = h(:)';
end

cos_phi = cos(phi);
sin_phi = sin(phi);
tan_phi = sin_phi./cos_phi;
h_a = h/def_ellps.a;
app = 1-def_ellps.f;

u = atan(app*tan_phi);
x = cos(u) + h_a.*cos_phi;
y = app*sin(u) + h_a.*sin_phi;

sg2.info5 = '** Geographic and geocentric coordinates of the local observers';
sg2.phi = phi;
sg2.lambda = lambda;
sg2.h = h;
sg2.a = def_ellps.a;
sg2.f = def_ellps.f;
sg2.u = u;
sg2.x = x;
sg2.y = y;

%% Angles topocentriques
xi = (def_ellps.a/ua);

omega_g = zeros(nj,np);
Delta_r_alpha = zeros(nj,np);
r_alpha = zeros(nj,np);
delta = zeros(nj,np);
omega = zeros(nj,np);
gamma_S0 = zeros(nj,np);
alpha_S = zeros(nj,np);

cos_delta_g = cos(delta_g);
omega_g0 = nu - r_alpha_g;

for kp = 1:np,

    omega_g(:,kp) = omega_g0 + lambda(kp);

    sin_phi_kp = sin_phi(kp);
    cos_phi_kp = cos_phi(kp);

    Delta_r_alpha(:,kp) = (-x(kp).*sin(omega_g(:,kp))./cos_delta_g*xi);
    r_alpha(:,kp) = r_alpha_g + Delta_r_alpha(:,kp);

    delta(:,kp) = delta_g + (x(kp)*cos(omega_g(:,kp)).*sin(delta_g)-y(kp)*cos_delta_g)*xi;

    omega(:,kp) = omega_g(:,kp)-Delta_r_alpha(:,kp);

    cos_omega_kp = cos(omega(:,kp));
    cos_delta_kp = cos(delta(:,kp));
    sin_delta_kp = sin(delta(:,kp));
    tan_delta_kp = sin_delta_kp./cos_delta_kp;

    gamma_S0(:,kp) = asin(sin_phi_kp.*sin_delta_kp + ...
        cos_phi_kp.*cos_delta_kp.*cos_omega_kp);

    alpha_S(:,kp) = atan2(sin(omega(:,kp)),cos_omega_kp.*sin_phi_kp-tan_delta_kp.*cos_phi_kp)+pi;

end

if (corr_refract)
    Dgamma_S_refract = sg2_refract_correction(gamma_S0, P, T, corr_refract_method);
else
    Dgamma_S_refract = zeros(size(gamma_S0));
end

gamma_S = gamma_S0 + Dgamma_S_refract;
theta_S = pi/2-gamma_S;
theta_S0 = pi/2-gamma_S0;

sg2.info6 = '** Topocentric solar angles';

sg2.omega_g = sg2_setrange(omega_g, 1); % [0, 2pi]
sg2.xi = xi;
sg2.Delta_r_alpha = sg2_setrange(Delta_r_alpha, 0);
sg2.r_alpha = sg2_setrange(r_alpha);
sg2.delta = sg2_setrange(delta,0);
sg2.omega = sg2_setrange(omega);
sg2.gamma_S0 = sg2_setrange(gamma_S0,0);
sg2.Dgamma_S_refract = sg2_setrange(Dgamma_S_refract,1);
sg2.gamma_S = sg2_setrange(gamma_S,0);
sg2.theta_S = sg2_setrange(theta_S,0);
sg2.theta_S0 = sg2_setrange(theta_S0,1);

idx_lat_pos = find(lat >= 0);
idx_lat_neg = find(lat < 0);

sg2.alpha_S(:,idx_lat_pos) = sg2_setrange(alpha_S(:,idx_lat_pos),1);
sg2.alpha_S(:,idx_lat_neg) = sg2_setrange(alpha_S(:,idx_lat_neg),0);

sg2.info7 = '** Incidence sur plan incline';
sg2.alpha = alpha;
sg2.beta = beta;
cos_theta_I = cos(theta_S).*cos(beta) + sin(beta).*sin(theta_S).*cos(alpha_S-alpha);
sg2.cos_theta_I = cos_theta_I;
sg2.theta_I = acos(cos_theta_I);

%% Julian date to YMD conversion
function ymd = sg2_jd_to_ymd(jd)

jd = jd(:);

H = (jd + 0.5 - floor(jd + 0.5)) * 24;
L = floor(jd + 0.5) + 68569;
N = floor(4 * L / 146097);
L = L - floor((146097 * N + 3) / 4);
I = floor(4000 * (L + 1) / 1461001);
L = L - floor(1461 * I / 4) + 31;

J = floor(80 * L / 2447);
K = L - floor(2447 * J / 80);
L = floor(J / 11);
J = J + 2 - 12 * L;
I = 100 * (N - 49) + I + L;

ymd = [I J K H];

%% Set range of angle matrix
function y = sg2_setrange(x, mode, K)

if (nargin < 2)
    mode = 1;
end
if (nargin < 3)
    K = 2*pi;
end

if (mode == 0)
    y = x - round(x/K)*K;
else
    y = x - floor(x/K)*K;
end

%% Delta_t computation with polynomial approximation
function delta_t = sg2_calc_delta_t(y, approx_DeltaT_MSc)

n = length(y);
delta_t = nan(n,1);

for kP = 1:approx_DeltaT_MSc.nP,
    idx_kP = find((y>=approx_DeltaT_MSc.y(kP,1))&(y<approx_DeltaT_MSc.y(kP,2)));
    if (~isempty(idx_kP))
        delta_t(idx_kP) = polyval(approx_DeltaT_MSc.P(kP).a, ...
            y(idx_kP) - approx_DeltaT_MSc.P(kP).y);
    end
end

%% Sinusoidal approximation
function s = sg2_calc_approxS(j, approx)

j = j(:);
np = length(j);

jc = j - approx.j0;

s = approx.a*jc + approx.b;

for kf = 1:approx.N,
    s = s + approx.ro(kf)*cos(approx.a0(kf)*jc-approx.phi(kf));
end

%% Reference ellipse definition
function ellps = sg2_ref_ellps(nom_ellps)

switch(lower(nom_ellps))
    case {'ntf' 'clarke1880'}
        ellps.a = 6378249.2;
        ellps.b = 6356515.0;
        ellps.e = sqrt((ellps.a^2-ellps.b^2)/ellps.a^2);
        ellps.f = 1-ellps.b/ellps.a;
    case 'wgs84'
        ellps.a = 6378137.0;
        ellps.f = 1/298.257223563;
        ellps.b = ellps.a*(1-ellps.f);
        ellps.e = sqrt((ellps.a^2-ellps.b^2)/ellps.a^2);
    case 'rgf83'
        ellps.a = 6378137.0;
        ellps.f = 1/298.257222101;
        ellps.b = ellps.a*(1-ellps.f);
        ellps.e = sqrt((ellps.a^2-ellps.b^2)/ellps.a^2);
    case 'aa'
        ellps.a = 6378136.6;
        ellps.f = 1/298.25642;
        ellps.b = ellps.a*(1-ellps.f);
        ellps.e = sqrt((ellps.a^2-ellps.b^2)/ellps.a^2);
    case 'spa'
        ellps.a = 6378140.0;
        ellps.f = 1-0.99664719;
        ellps.b = ellps.a*(1-ellps.f);
        ellps.e = sqrt((ellps.a^2-ellps.b^2)/ellps.a^2);
    case 'sphere'
        ellps.a = 637103.0;
        ellps.f = 0.0;
        ellps.b = ellps.a*(1-ellps.f);
        ellps.e = sqrt((ellps.a^2-ellps.b^2)/ellps.a^2);
    otherwise
        error('Type ELLPS non prévu');
end

%% Atmospheric corrections of solar elevation
function Dgamma = sg2_refract_correction(gamma_S0, P, T, method)

deg2rad = pi/180;
method = strtrim(method);

gamma_S0 = gamma_S0 - round(gamma_S0/(2*pi))*2*pi;

switch(lower(method(1:3)))
    case 'sae',

        gamma_S0_seuil = -0.010035643198967;
        K0 = (P/1010)*(283/(273+T))*2.96706e-4;

        Dgamma = zeros(size(gamma_S0));
        idx = find(gamma_S0>gamma_S0_seuil);
        Dgamma(idx) = K0./(tan(gamma_S0(idx)+0.0031376./(gamma_S0(idx)+0.089186)));
        idx = find(gamma_S0<=gamma_S0_seuil);
        K = K0/(tan(gamma_S0_seuil+0.0031376/(gamma_S0_seuil+0.089186)));
        Dgamma(idx) = K*tan(gamma_S0_seuil)./tan(gamma_S0(idx));

    case 'zim',

        K = (P/1013.0)*(283/(273+T))*pi/180/3600;
        tan_gamma_S0 = tan(gamma_S0);
        Dgamma = zeros(size(gamma_S0));

        idx = find((gamma_S0>=0.087266462599716)&(gamma_S0<=1.483529864195180)); %% [5°, 85[
        Dgamma(idx) = K*(58.1./tan_gamma_S0(idx) ...
            - 0.07./tan_gamma_S0(idx).^3 + 0.000086./tan_gamma_S0(idx).^5);

        idx = find((gamma_S0>=-0.010036)&(gamma_S0<=0.087266)); %% [5°, 85[
        Dgamma(idx) = (1735-2.969067e4*gamma_S0(idx)+...
            3.394422e5*gamma_S0(idx).^2 + ...
            -2.405683e6*gamma_S0(idx).^3 + ...
            7.66231727e6*gamma_S0(idx).^4)*K;
        idx = find(gamma_S0<=-0.010036);
        Dgamma(idx) = (-20.774./tan_gamma_S0(idx))*K;

    otherwise
        error('Unknow correction method');
end

