ssol_ranst_sun_dir.c (11570B)
1 /* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016, 2018 CNRS 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #include "ssol.h" 18 #include "ssol_ranst_sun_dir.h" 19 20 #include <star/ssp.h> 21 22 #include <rsys/double33.h> 23 #include <rsys/math.h> 24 #include <rsys/mem_allocator.h> 25 #include <rsys/rsys.h> 26 #include <rsys/ref_count.h> 27 28 /******************************************************************************* 29 * Distributions types 30 ******************************************************************************/ 31 struct ran_buie_state { 32 double thetaSD; 33 double deltaThetaCSSD; 34 double gamma; 35 double etokTimes1000toGamma; 36 double alpha; 37 double hRect1; 38 double hRect2; 39 double proba_rect1; 40 double basis[9]; 41 }; 42 43 struct ran_pillbox_state { 44 double sin2_theta_max; 45 double basis[9]; 46 }; 47 48 struct ran_gaussian_state { 49 double std_dev; 50 double basis[9]; 51 }; 52 53 struct ran_dirac_state { 54 double dir[3]; 55 }; 56 57 /* One single type for all distributions. Only the state type depends on the 58 * distribution type */ 59 struct ranst_sun_dir { 60 double*(*get) 61 (const struct ranst_sun_dir* ran, struct ssp_rng* rng, double dir[3]); 62 union { 63 struct ran_buie_state buie; 64 struct ran_pillbox_state pillbox; 65 struct ran_gaussian_state gaussian; 66 struct ran_dirac_state dirac; 67 } state; 68 69 ref_T ref; 70 struct mem_allocator* allocator; 71 }; 72 73 /******************************************************************************* 74 * Helper functions 75 ******************************************************************************/ 76 static void 77 distrib_sun_release(ref_T* ref) 78 { 79 struct ranst_sun_dir* ran; 80 ASSERT(ref); 81 ran = CONTAINER_OF(ref, struct ranst_sun_dir, ref); 82 MEM_RM(ran->allocator, ran); 83 } 84 85 /******************************************************************************* 86 * Buie distribution 87 ******************************************************************************/ 88 static FINLINE double 89 chi_value(const double csr) 90 { 91 if (csr > 0.145) { 92 return -0.04419909985804843 + 93 csr * (1.401323894233574 + 94 csr * (-0.3639746714505299 + 95 csr * (-0.9579768560161194 + 1.1550475450828657 * csr))); 96 } 97 if (csr > 0.035) { 98 return 0.022652077593662934 + 99 csr * (0.5252380349996234 + 100 csr * (2.5484334534423887 - 0.8763755326550412 * csr)); 101 } 102 return 0.004733749294807862 + 103 csr * (4.716738065192151 + 104 csr * (-463.506669149804 + 105 csr * (24745.88727411664 + 106 csr * (-606122.7511711778 + 5521693.445014727 * csr)))); 107 } 108 109 static FINLINE double 110 phi_solar_disk(const double theta) 111 { 112 /* The parameter theta is the zenith angle in radians */ 113 return cos(326 * theta) / cos(308 * theta); 114 } 115 116 static FINLINE double 117 phi_circum_solor_region(const double theta, const struct ran_buie_state* state) 118 { 119 /* The parameter theta is the zenith angle in radians */ 120 return state->etokTimes1000toGamma * pow(theta, state->gamma); 121 } 122 123 static FINLINE double 124 phi(const double theta, const struct ran_buie_state* state) 125 { 126 /* The parameter theta is the zenith angle in radians */ 127 if (theta < state->thetaSD) return phi_solar_disk(theta); 128 else return phi_circum_solor_region(theta, state); 129 } 130 131 static FINLINE double 132 pdf_theta(const double theta, const struct ran_buie_state* state) 133 { 134 /* The parameter theta is the zenith angle in radians */ 135 return state->alpha * phi(theta, state) * sin(theta); 136 } 137 138 static FINLINE double 139 gamma_value(const double chi) 140 { 141 /* gamma is the gradient of the 2nd part of the curve in log/log space*/ 142 return 2.2 * log(0.52 * chi) * pow(chi, 0.43) - 0.1; 143 } 144 145 static FINLINE double 146 k_value(const double chi) 147 { 148 /* k is the intercept of the 2nd part of the curve at an angular displacement 149 * of zero */ 150 return 0.9 * log(13.5 * chi) * pow(chi, -0.3); 151 } 152 153 static FINLINE double 154 integral_B 155 (const double k, 156 const double gamma, 157 const double thetaCS, 158 const double thetaSD) 159 { 160 const double g2 = gamma + 2.0; 161 return exp(k) * pow(1000, gamma) / g2 * (pow(thetaCS, g2) - pow(thetaSD, g2)); 162 } 163 164 static FINLINE double 165 proba_rect1 166 (const double widthR1, 167 const double heightR1, 168 const double widthR2, 169 const double heightR2) 170 { 171 const double areaR1 = widthR1 * heightR1; 172 const double areaR2 = widthR2 * heightR2; 173 return areaR1 / (areaR1 + areaR2); 174 } 175 176 static FINLINE double 177 zenith_angle(struct ssp_rng* rng, const struct ran_buie_state* state) 178 { 179 double theta; 180 double value; 181 ASSERT(state); 182 183 do { 184 if (ssp_rng_canonical(rng) < state->proba_rect1) { 185 theta = ssp_rng_canonical(rng) * state->thetaSD; 186 value = ssp_rng_canonical(rng) * state->hRect1; 187 } 188 else { 189 theta = state->thetaSD + ssp_rng_canonical(rng) * state->deltaThetaCSSD; 190 value = ssp_rng_canonical(rng) * state->hRect2; 191 } 192 } while (value > pdf_theta(theta, state)); 193 194 return theta; 195 } 196 197 static void 198 fill_buie_state(struct ran_buie_state* s, const double p) 199 { 200 double thetaCS, integralA, chi, k, integralB; 201 ASSERT(s); 202 ASSERT(0 <= p && p < 1); 203 204 s->thetaSD = 0.00465; 205 thetaCS = 0.0436; 206 s->deltaThetaCSSD = thetaCS - s->thetaSD; 207 integralA = 9.224724736098827E-6; 208 chi = chi_value(p); 209 k = k_value(chi); 210 s->gamma = gamma_value(chi); 211 s->etokTimes1000toGamma = exp(k) * pow(1000, s->gamma); 212 integralB = integral_B(k, s->gamma, thetaCS, s->thetaSD); 213 s->alpha = 1 / (integralA + integralB); 214 /* The value 0.0038915695846209047 radians is the value for which the 215 * probability density function is a maximum. The value of the maximum varies 216 * with the circumsolar ratio, but the value of theta where the maximum is 217 * obtained remains fixed. */ 218 s->hRect1 = 1.001 * pdf_theta(0.0038915695846209047, s); 219 s->hRect2 = pdf_theta(s->thetaSD, s); 220 s->proba_rect1 = proba_rect1 221 (s->thetaSD, s->hRect1, s->deltaThetaCSSD, s->hRect2); 222 } 223 224 static double* 225 ran_buie_get 226 (const struct ranst_sun_dir* ran, struct ssp_rng* rng, double dir[3]) 227 { 228 double phi, theta, sinTheta, cosTheta, cosPhi, sinPhi; 229 ASSERT(ran->state.buie.thetaSD > 0); 230 phi = ssp_rng_uniform_double(rng, 0, 2 * PI); 231 theta = zenith_angle(rng, &ran->state.buie); 232 sinTheta = sin(theta); 233 cosTheta = cos(theta); 234 cosPhi = cos(phi); 235 sinPhi = sin(phi); 236 dir[0] = sinTheta * sinPhi; 237 dir[1] = sinTheta * cosPhi; 238 dir[2] = cosTheta; 239 d33_muld3(dir, ran->state.buie.basis, dir); 240 return dir; 241 } 242 243 /******************************************************************************* 244 * Pillbox random variate 245 ******************************************************************************/ 246 static double* 247 ran_pillbox_get 248 (const struct ranst_sun_dir* ran, 249 struct ssp_rng* rng, 250 double dir[3]) 251 { 252 double pt[3]; 253 double phi, sin2_theta, sin_theta, cos_theta; 254 255 ASSERT(ran && 0 <= ran->state.pillbox.sin2_theta_max 256 && ran->state.pillbox.sin2_theta_max <= 1); 257 sin2_theta = ssp_rng_uniform_double(rng, 0, ran->state.pillbox.sin2_theta_max); 258 sin_theta = sqrt(sin2_theta); 259 cos_theta = sqrt(1 - sin2_theta); 260 phi = ssp_rng_uniform_double(rng, 0, 2 * PI); 261 pt[0] = cos(phi) * sin_theta; 262 pt[1] = sin(phi) * sin_theta; 263 pt[2] = cos_theta; 264 d33_muld3(dir, ran->state.pillbox.basis, pt); 265 d3_normalize(dir, dir); 266 return dir; 267 } 268 269 /******************************************************************************* 270 * Gaussian random variate 271 ******************************************************************************/ 272 static double* 273 ran_gaussian_get 274 (const struct ranst_sun_dir* ran, 275 struct ssp_rng* rng, 276 double dir[3]) 277 { 278 double pt[3]; 279 double phi, cos_theta, sin_theta; 280 281 /* The following is not truly a gaussian sunshape, 282 * but an accurate enough approximation for small angles */ 283 ASSERT(ran && 0 <= ran->state.gaussian.std_dev); 284 sin_theta 285 = ran->state.gaussian.std_dev * sqrt(-2 * log(ssp_rng_canonical(rng))); 286 sin_theta = MMIN(1, sin_theta); /* macro: don't merge with previous line! */ 287 cos_theta = sin2cos(sin_theta); 288 phi = ssp_rng_uniform_double(rng, 0, 2 * PI); 289 pt[0] = cos(phi) * sin_theta; 290 pt[1] = sin(phi) * sin_theta; 291 pt[2] = cos_theta; 292 d33_muld3(dir, ran->state.gaussian.basis, pt); 293 d3_normalize(dir, dir); 294 return dir; 295 } 296 297 /******************************************************************************* 298 * Dirac distribution 299 ******************************************************************************/ 300 static double* 301 ran_dirac_get 302 (const struct ranst_sun_dir* ran, 303 struct ssp_rng* rng, 304 double dir[3]) 305 { 306 (void) rng; 307 ASSERT(d3_is_normalized(ran->state.dirac.dir)); 308 d3_set(dir, ran->state.dirac.dir); 309 return dir; 310 } 311 312 /******************************************************************************* 313 * Local functions 314 ******************************************************************************/ 315 res_T 316 ranst_sun_dir_create 317 (struct mem_allocator* allocator, 318 struct ranst_sun_dir** out_ran) 319 { 320 struct ranst_sun_dir* ran = NULL; 321 322 if (!out_ran) return RES_BAD_ARG; 323 324 allocator = allocator ? allocator : &mem_default_allocator; 325 326 ran = MEM_CALLOC(allocator, 1, sizeof(struct ranst_sun_dir)); 327 if (!ran) return RES_MEM_ERR; 328 329 ref_init(&ran->ref); 330 ran->allocator = allocator; 331 *out_ran = ran; 332 333 return RES_OK; 334 } 335 336 res_T 337 ranst_sun_dir_ref_get(struct ranst_sun_dir* ran) 338 { 339 if (!ran) return RES_BAD_ARG; 340 ref_get(&ran->ref); 341 return RES_OK; 342 } 343 344 res_T 345 ranst_sun_dir_ref_put(struct ranst_sun_dir* ran) 346 { 347 if (!ran) return RES_BAD_ARG; 348 ref_put(&ran->ref, distrib_sun_release); 349 return RES_OK; 350 } 351 352 double* 353 ranst_sun_dir_get 354 (const struct ranst_sun_dir* ran, 355 struct ssp_rng* rng, 356 double dir[3]) 357 { 358 ASSERT(ran); 359 return ran->get(ran, rng, dir); 360 } 361 362 res_T 363 ranst_sun_dir_buie_setup 364 (struct ranst_sun_dir* ran, 365 const double param, 366 const double dir[3]) 367 { 368 const double minCRSValue = 1E-6; 369 const double maxCRSValue = 0.849; 370 if (!ran || !dir || param < minCRSValue || param > maxCRSValue) 371 return RES_BAD_ARG; 372 373 ran->get = ran_buie_get; 374 d33_basis(ran->state.buie.basis, dir); 375 fill_buie_state(&ran->state.buie, param); 376 return RES_OK; 377 } 378 379 res_T 380 ranst_sun_dir_pillbox_setup 381 (struct ranst_sun_dir* ran, 382 const double theta_max, 383 const double dir[3]) 384 { 385 double s; 386 if (!ran || !dir || theta_max <= 0 || theta_max > PI ) 387 return RES_BAD_ARG; 388 ran->get = ran_pillbox_get; 389 s = sin(theta_max); 390 ran->state.pillbox.sin2_theta_max = s * s; 391 d33_basis(ran->state.pillbox.basis, dir); 392 return RES_OK; 393 } 394 395 res_T 396 ranst_sun_dir_gaussian_setup 397 (struct ranst_sun_dir* ran, 398 const double std_dev, 399 const double dir[3]) 400 { 401 if(!ran || !dir || std_dev <= 0) return RES_BAD_ARG; 402 ran->get = ran_gaussian_get; 403 ran->state.gaussian.std_dev = std_dev; 404 d33_basis(ran->state.pillbox.basis, dir); 405 return RES_OK; 406 } 407 408 res_T 409 ranst_sun_dir_dirac_setup 410 (struct ranst_sun_dir* ran, 411 const double dir[3]) 412 { 413 if (!ran || !dir) return RES_BAD_ARG; 414 if (0 == d3_normalize(ran->state.dirac.dir, dir)) { 415 return RES_BAD_ARG; /* zero vector */ 416 } 417 ran->get = ran_dirac_get; 418 return RES_OK; 419 } 420