#ifdef __OPENCL_VERSION__ # define USE_OPENCL #elif defined(_OPENMP) # define USE_OPENMP #endif // If opencl is not available, then we are compiling a C function // Note: if using a C++ compiler, then define kernel as extern "C" #ifdef USE_OPENCL typedef int int32_t; # if defined(USE_SINCOS) # define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar) # else # define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) # endif // Intel CPU on Mac gives strange values for erf(); on the verified // platforms (intel, nvidia, amd), the cephes erf() is significantly // faster than that available in the native OpenCL. #define NEED_ERF // OpenCL only has type generic math #define expf exp #ifndef NEED_ERF # define erff erf # define erfcf erfc #endif #else // !USE_OPENCL // Use SAS_DOUBLE to force the use of double even for float kernels # define SAS_DOUBLE dou ## ble # ifdef __cplusplus #include #include using namespace std; #if defined(_MSC_VER) #include #include #define kernel extern "C" __declspec( dllexport ) inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } inline double fmin(double x, double y) { return x>y ? y : x; } inline double fmax(double x, double y) { return x::quiet_NaN()) // non-signalling NaN #define INFINITY (std::numeric_limits::infinity()) #define NEED_ERF #define NEED_EXPM1 #define NEED_TGAMMA #else #define kernel extern "C" #include #endif inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); } # else // !__cplusplus #include // C99 guarantees that int32_t types is here #include #if defined(__TINYC__) typedef int int32_t; #include // TODO: check isnan is correct inline double _isnan(double x) { return x != x; } // hope this doesn't optimize away! #undef isnan #define isnan(x) _isnan(x) // Defeat the double->float conversion since we don't have tgmath inline SAS_DOUBLE trunc(SAS_DOUBLE x) { return x>=0?floor(x):-floor(-x); } inline SAS_DOUBLE fmin(SAS_DOUBLE x, SAS_DOUBLE y) { return x>y ? y : x; } inline SAS_DOUBLE fmax(SAS_DOUBLE x, SAS_DOUBLE y) { return x // C99 type-generic math, so sin(float) => sinf #endif // MSVC doesn't support C99, so no need for dllexport on C99 branch #define kernel #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) # endif // !__cplusplus # define global # define local # define constant const // OpenCL powr(a,b) = C99 pow(a,b), b >= 0 // OpenCL pown(a,b) = C99 pow(a,b), b integer # define powr(a,b) pow(a,b) # define pown(a,b) pow(a,b) #endif // !USE_OPENCL #if defined(NEED_EXPM1) // TODO: precision is a half digit lower than numpy on mac in [1e-7, 0.5] // Run "explore/precision.py sas_expm1" to see this (may have to fiddle // the xrange for log to see the complete range). static SAS_DOUBLE expm1(SAS_DOUBLE x_in) { double x = (double)x_in; // go back to float for single precision kernels // Adapted from the cephes math library. // Copyright 1984 - 1992 by Stephen L. Moshier if (x != x || x == 0.0) { return x; // NaN and +/- 0 } else if (x < -0.5 || x > 0.5) { return exp(x) - 1.0; } else { const double xsq = x*x; const double p = ((( +1.2617719307481059087798E-4)*xsq +3.0299440770744196129956E-2)*xsq +9.9999999999999999991025E-1); const double q = (((( +3.0019850513866445504159E-6)*xsq +2.5244834034968410419224E-3)*xsq +2.2726554820815502876593E-1)*xsq +2.0000000000000000000897E0); double r = x * p; r = r / (q - r); return r+r; } } #endif // Standard mathematical constants: // M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4, // M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2) // OpenCL defines M_constant_F for float constants, and nothing if double // is not enabled on the card, which is why these constants may be missing #ifndef M_PI # define M_PI 3.141592653589793 #endif #ifndef M_PI_2 # define M_PI_2 1.570796326794897 #endif #ifndef M_PI_4 # define M_PI_4 0.7853981633974483 #endif #ifndef M_E # define M_E 2.718281828459045091 #endif #ifndef M_SQRT1_2 # define M_SQRT1_2 0.70710678118654746 #endif // Non-standard function library // pi/180, used for converting between degrees and radians // 4/3 pi for computing sphere volumes // square and cube for computing squares and cubes #ifndef M_PI_180 # define M_PI_180 0.017453292519943295 #endif #ifndef M_4PI_3 # define M_4PI_3 4.18879020478639 #endif inline double square(double x) { return x*x; } inline double cube(double x) { return x*x*x; } inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } // CRUFT: support old style models with orientation received qx, qy and angles // To rotate from the canonical position to theta, phi, psi, first rotate by // psi about the major axis, oriented along z, which is a rotation in the // detector plane xy. Next rotate by theta about the y axis, aligning the major // axis in the xz plane. Finally, rotate by phi in the detector plane xy. // To compute the scattering, undo these rotations in reverse order: // rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi // The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit // vector in the q direction. // To change between counterclockwise and clockwise rotation, change the // sign of phi and psi. #if 1 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ SINCOS(phi*M_PI_180, sn, cn); \ q = sqrt(qx*qx + qy*qy); \ cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180)); \ sn = sqrt(1 - cn*cn); \ } while (0) #else // SasView 3.x definition of orientation #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ SINCOS(theta*M_PI_180, sn, cn); \ q = sqrt(qx*qx + qy*qy);\ cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ sn = sqrt(1 - cn*cn); \ } while (0) #endif #if 1 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ q = sqrt(qx*qx + qy*qy); \ const double qxhat = qx/q; \ const double qyhat = qy/q; \ double sin_theta, cos_theta; \ double sin_phi, cos_phi; \ double sin_psi, cos_psi; \ SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ zhat = qxhat*(-sin_theta*cos_phi) \ + qyhat*(-sin_theta*sin_phi); \ } while (0) #else // SasView 3.x definition of orientation #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ q = sqrt(qx*qx + qy*qy); \ const double qxhat = qx/q; \ const double qyhat = qy/q; \ double sin_theta, cos_theta; \ double sin_phi, cos_phi; \ double sin_psi, cos_psi; \ SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ } while (0) #endif