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))&(ygamma_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