1 /* Copyright (c) 2019, The University of Oxford. See LICENSE file. */ 2 3 /* See "Spherical Harmonic Lighting: The Gritty Details", 4 * Robin Green, GDC, 2003. */ 5 6 #define OSKAR_SPHERICAL_HARMONIC_SUM_REAL(NAME, FP)\ 7 KERNEL(NAME) (\ 8 const int num_points,\ 9 GLOBAL_IN(FP, theta),\ 10 GLOBAL_IN(FP, phi),\ 11 const int l_max,\ 12 GLOBAL_IN(FP, coeff),\ 13 const int stride,\ 14 const int offset_out,\ 15 GLOBAL_OUT(FP, surface))\ 16 {\ 17 KERNEL_LOOP_PAR_X(int, i, 0, num_points)\ 18 FP sum = (FP)0;\ 19 const FP phi_ = phi[i], theta_ = theta[i];\ 20 if (phi_ != phi_) sum = phi_; /* Propagate NAN. */\ 21 else {\ 22 FP sin_t, cos_t;\ 23 SINCOS(theta_, sin_t, cos_t);\ 24 for (int l = 0; l <= l_max; ++l) {\ 25 FP val;\ 26 const int ind0 = l * l + l;\ 27 const FP f_ = (2 * l + 1) / (4 * (FP)M_PI);\ 28 const FP c0 = coeff[ind0];\ 29 if (c0 != (FP)0) {\ 30 OSKAR_LEGENDRE1(FP, l, 0, cos_t, sin_t, val)\ 31 sum += c0 * val * sqrt(f_);\ 32 }\ 33 for (int abs_m = 1; abs_m <= l; ++abs_m) {\ 34 FP sin_p, cos_p, d_fact = (FP)1, s_fact = (FP)1;\ 35 const FP cmm = coeff[ind0 - abs_m], cpm = coeff[ind0 + abs_m];\ 36 if ((cpm == (FP)0) && (cmm == (FP)0)) continue;\ 37 OSKAR_LEGENDRE1(FP, l, abs_m, cos_t, sin_t, val)\ 38 const int d_ = l - abs_m, s_ = l + abs_m;\ 39 for (int i_ = 2; i_ <= d_; ++i_) d_fact *= i_;\ 40 for (int i_ = 2; i_ <= s_; ++i_) s_fact *= i_;\ 41 val *= sqrt(2 * f_ * d_fact / s_fact);\ 42 const FP p = abs_m * phi_;\ 43 SINCOS(p, sin_p, cos_p);\ 44 sum += (val * (cpm * cos_p + cmm * sin_p));\ 45 }\ 46 }\ 47 }\ 48 surface[i * stride + offset_out] = sum;\ 49 KERNEL_LOOP_END\ 50 }\ 51 OSKAR_REGISTER_KERNEL(NAME) 52