source: sasmodels/sasmodels/kernel_iq.c @ 7c35fda

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 7c35fda was 7c35fda, checked in by Adam Washington <adam.washington@…>, 6 years ago

Document reason for change in I(q)

  • Property mode set to 100644
File size: 26.1 KB
RevLine 
[2e44ac7]1/*
2    ##########################################################
3    #                                                        #
4    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
5    #   !!                                              !!   #
6    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
7    #   !!                                              !!   #
8    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
9    #                                                        #
10    ##########################################################
11*/
12
[9ee2756]13// NOTE: the following macros are defined in generate.py:
14//
[6aee3ab]15//  MAX_PD : the maximum number of dispersity loops allowed for this model,
16//      which will be at most modelinfo.MAX_PD.
[9ee2756]17//  NUM_PARS : the number of parameters in the parameter table
18//  NUM_VALUES : the number of values to skip at the start of the
19//      values array before you get to the dispersity values.
20//  PARAMETER_TABLE : list of parameter declarations used to create the
21//      ParameterTable type.
22//  KERNEL_NAME : model_Iq, model_Iqxy or model_Imagnetic.  This code is
23//      included three times, once for each kernel type.
[2c108a3]24//  MAGNETIC : defined when the magnetic kernel is being instantiated
25//  NUM_MAGNETIC : the number of magnetic parameters
[9ee2756]26//  MAGNETIC_PARS : a comma-separated list of indices to the sld
27//      parameters in the parameter table.
28//  CALL_VOLUME(table) : call the form volume function
29//  CALL_IQ(q, table) : call the Iq function for 1D calcs.
30//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data.
31//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes
32//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes
[108e70e]33//  CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models
[9ee2756]34//  INVALID(table) : test if the current point is feesible to calculate.  This
35//      will be defined in the kernel definition file.
[ff10479]36//  PROJECTION : equirectangular=1, sinusoidal=2
37//      see explore/jitter.py for definitions.
[767dca8]38
[03cac08]39#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
40#define _PAR_BLOCK_
[2e44ac7]41
42typedef struct {
[60eab2a]43#if MAX_PD > 0
[2c108a3]44    int32_t pd_par[MAX_PD];     // id of the nth dispersity variable
45    int32_t pd_length[MAX_PD];  // length of the nth dispersity weight vector
[0a7e5eb4]46    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
[5cf3c33]47    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
[60eab2a]48#endif // MAX_PD > 0
[bde38b5]49    int32_t num_eval;           // total number of voxels in hypercube
50    int32_t num_weights;        // total length of the weights vector
[5ff1b03]51    int32_t num_active;         // number of non-trivial pd loops
[8698a0d]52    int32_t theta_par;          // id of first orientation variable
[2e44ac7]53} ProblemDetails;
54
[bde38b5]55// Intel HD 4000 needs private arrays to be a multiple of 4 long
[2e44ac7]56typedef struct {
[56547a8]57    PARAMETER_TABLE
[bde38b5]58} ParameterTable;
59typedef union {
60    ParameterTable table;
61    double vector[4*((NUM_PARS+3)/4)];
[2e44ac7]62} ParameterBlock;
[9eb3632]63#endif // _PAR_BLOCK_
[03cac08]64
[9ee2756]65#if defined(MAGNETIC) && NUM_MAGNETIC > 0
66// ===== Helper functions for magnetism =====
[a4280bd]67
[32e3c9b]68// Return value restricted between low and high
69static double clip(double value, double low, double high)
70{
[b966a96]71  return (value < low ? low : (value > high ? high : value));
[32e3c9b]72}
73
74// Compute spin cross sections given in_spin and out_spin
75// To convert spin cross sections to sld b:
76//     uu * (sld - m_sigma_x);
77//     dd * (sld + m_sigma_x);
[2c108a3]78//     ud * (m_sigma_y - 1j*m_sigma_z);
79//     du * (m_sigma_y + 1j*m_sigma_z);
80// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag
[885753a]81static void set_spin_weights(double in_spin, double out_spin, double weight[6])
[32e3c9b]82{
[b966a96]83  in_spin = clip(in_spin, 0.0, 1.0);
84  out_spin = clip(out_spin, 0.0, 1.0);
[7c35fda]85  // Previous version of this function took the square root of the weights,
86  // under the assumption that
87  //
[885753a]88  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
[7c35fda]89  //
90  // However, since the weights are applied to the final intensity and
91  // are not interned inside the I(q) function, we want the full
92  // weight and not the square root.  Any function using
93  // set_spin_weights as part of calculating an amplitude will need to
94  // manually take that square root, but there is currently no such
95  // function.
[84e89cd]96  weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd
97  weight[1] = (1.0-in_spin) * out_spin;       // du
98  weight[2] = in_spin * (1.0-out_spin);       // ud
99  weight[3] = in_spin * out_spin;             // uu
[885753a]100  weight[4] = weight[1]; // du.imag
101  weight[5] = weight[2]; // ud.imag
[a4280bd]102}
103
[9ee2756]104// Compute the magnetic sld
[2c108a3]105static double mag_sld(
[885753a]106  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag
[2c108a3]107  const double qx, const double qy,
108  const double px, const double py,
109  const double sld,
110  const double mx, const double my, const double mz
111)
[a4280bd]112{
[2c108a3]113  if (xs < 4) {
[a4280bd]114    const double perp = qy*mx - qx*my;
[2c108a3]115    switch (xs) {
[aadec17]116      default: // keep compiler happy; condition ensures xs in [0,1,2,3]
[2c108a3]117      case 0: // uu => sld - D M_perpx
118          return sld - px*perp;
[885753a]119      case 1: // ud.real => -D M_perpy
[2c108a3]120          return py*perp;
[885753a]121      case 2: // du.real => -D M_perpy
[2c108a3]122          return py*perp;
[885753a]123      case 3: // dd => sld + D M_perpx
[2c108a3]124          return sld + px*perp;
125    }
126  } else {
127    if (xs== 4) {
[885753a]128      return -mz;  // du.imag => +D M_perpz
[2c108a3]129    } else { // index == 5
[885753a]130      return +mz;  // ud.imag => -D M_perpz
[2c108a3]131    }
132  }
[32e3c9b]133}
[9ee2756]134
[2c108a3]135
[9ee2756]136#endif
137
138// ===== Helper functions for orientation and jitter =====
139
[8698a0d]140// To change the definition of the angles, run explore/angles.py, which
141// uses sympy to generate the equations.
142
[9ee2756]143#if !defined(_QAC_SECTION) && defined(CALL_IQ_AC)
144#define _QAC_SECTION
145
146typedef struct {
147    double R31, R32;
148} QACRotation;
149
150// Fill in the rotation matrix R from the view angles (theta, phi) and the
151// jitter angles (dtheta, dphi).  This matrix can be applied to all of the
152// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qc]'
153static void
154qac_rotation(
155    QACRotation *rotation,
156    double theta, double phi,
157    double dtheta, double dphi)
[8698a0d]158{
159    double sin_theta, cos_theta;
160    double sin_phi, cos_phi;
161
[9ee2756]162    // reverse view matrix
[8698a0d]163    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
164    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
[9ee2756]165    const double V11 = cos_phi*cos_theta;
166    const double V12 = sin_phi*cos_theta;
167    const double V21 = -sin_phi;
168    const double V22 = cos_phi;
169    const double V31 = sin_theta*cos_phi;
170    const double V32 = sin_phi*sin_theta;
[8698a0d]171
[9ee2756]172    // reverse jitter matrix
173    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta);
174    SINCOS(dphi*M_PI_180, sin_phi, cos_phi);
175    const double J31 = sin_theta;
176    const double J32 = -sin_phi*cos_theta;
177    const double J33 = cos_phi*cos_theta;
[8698a0d]178
[9ee2756]179    // reverse matrix
180    rotation->R31 = J31*V11 + J32*V21 + J33*V31;
181    rotation->R32 = J31*V12 + J32*V22 + J33*V32;
[8698a0d]182}
183
[9ee2756]184// Apply the rotation matrix returned from qac_rotation to the point (qx,qy),
185// returning R*[qx,qy]' = [qa,qc]'
[d86f0fc]186static void
[9ee2756]187qac_apply(
[ec8d4ac]188    QACRotation *rotation,
[9ee2756]189    double qx, double qy,
190    double *qa_out, double *qc_out)
[8698a0d]191{
[ec8d4ac]192    const double dqc = rotation->R31*qx + rotation->R32*qy;
[9ee2756]193    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
194    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
[8698a0d]195
[9ee2756]196    *qa_out = dqa;
197    *qc_out = dqc;
[8698a0d]198}
[9ee2756]199#endif // _QAC_SECTION
200
201#if !defined(_QABC_SECTION) && defined(CALL_IQ_ABC)
202#define _QABC_SECTION
203
204typedef struct {
205    double R11, R12;
206    double R21, R22;
207    double R31, R32;
208} QABCRotation;
209
210// Fill in the rotation matrix R from the view angles (theta, phi, psi) and the
211// jitter angles (dtheta, dphi, dpsi).  This matrix can be applied to all of the
212// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qb,qc]'
[8698a0d]213static void
[9ee2756]214qabc_rotation(
215    QABCRotation *rotation,
[8698a0d]216    double theta, double phi, double psi,
[9ee2756]217    double dtheta, double dphi, double dpsi)
[8698a0d]218{
219    double sin_theta, cos_theta;
220    double sin_phi, cos_phi;
221    double sin_psi, cos_psi;
222
223    // reverse view matrix
224    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
225    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
226    SINCOS(psi*M_PI_180, sin_psi, cos_psi);
[9ee2756]227    const double V11 = -sin_phi*sin_psi + cos_phi*cos_psi*cos_theta;
228    const double V12 = sin_phi*cos_psi*cos_theta + sin_psi*cos_phi;
229    const double V21 = -sin_phi*cos_psi - sin_psi*cos_phi*cos_theta;
230    const double V22 = -sin_phi*sin_psi*cos_theta + cos_phi*cos_psi;
[8698a0d]231    const double V31 = sin_theta*cos_phi;
232    const double V32 = sin_phi*sin_theta;
233
234    // reverse jitter matrix
[9ee2756]235    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta);
236    SINCOS(dphi*M_PI_180, sin_phi, cos_phi);
237    SINCOS(dpsi*M_PI_180, sin_psi, cos_psi);
[8698a0d]238    const double J11 = cos_psi*cos_theta;
[9ee2756]239    const double J12 = sin_phi*sin_theta*cos_psi + sin_psi*cos_phi;
240    const double J13 = sin_phi*sin_psi - sin_theta*cos_phi*cos_psi;
241    const double J21 = -sin_psi*cos_theta;
242    const double J22 = -sin_phi*sin_psi*sin_theta + cos_phi*cos_psi;
243    const double J23 = sin_phi*cos_psi + sin_psi*sin_theta*cos_phi;
[8698a0d]244    const double J31 = sin_theta;
245    const double J32 = -sin_phi*cos_theta;
246    const double J33 = cos_phi*cos_theta;
247
248    // reverse matrix
[9ee2756]249    rotation->R11 = J11*V11 + J12*V21 + J13*V31;
250    rotation->R12 = J11*V12 + J12*V22 + J13*V32;
251    rotation->R21 = J21*V11 + J22*V21 + J23*V31;
252    rotation->R22 = J21*V12 + J22*V22 + J23*V32;
253    rotation->R31 = J31*V11 + J32*V21 + J33*V31;
254    rotation->R32 = J31*V12 + J32*V22 + J33*V32;
[8698a0d]255}
256
[9ee2756]257// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy),
258// returning R*[qx,qy]' = [qa,qb,qc]'
[d86f0fc]259static void
[9ee2756]260qabc_apply(
[ec8d4ac]261    QABCRotation *rotation,
[9ee2756]262    double qx, double qy,
263    double *qa_out, double *qb_out, double *qc_out)
[8698a0d]264{
[ec8d4ac]265    *qa_out = rotation->R11*qx + rotation->R12*qy;
266    *qb_out = rotation->R21*qx + rotation->R22*qy;
267    *qc_out = rotation->R31*qx + rotation->R32*qy;
[8698a0d]268}
269
[9ee2756]270#endif // _QABC_SECTION
271
272
273// ==================== KERNEL CODE ========================
[2e44ac7]274
[03cac08]275kernel
276void KERNEL_NAME(
[5cf3c33]277    int32_t nq,                 // number of q values
[2c108a3]278    const int32_t pd_start,     // where we are in the dispersity loop
279    const int32_t pd_stop,      // where we are stopping in the dispersity loop
[6e7ff6d]280    global const ProblemDetails *details,
[9eb3632]281    global const double *values,
[2e44ac7]282    global const double *q, // nq q values, with padding to boundary
[9eb3632]283    global double *result,  // nq+1 return values, again with padding
[2c108a3]284    const double cutoff     // cutoff in the dispersity weight product
[2e44ac7]285    )
286{
[9ee2756]287#ifdef USE_OPENCL
288  // who we are and what element we are working with
289  const int q_index = get_global_id(0);
290  if (q_index >= nq) return;
291#else
292  // Define q_index here so that debugging statements can be written to work
293  // for both OpenCL and DLL using:
294  //    if (q_index == 0) {printf(...);}
295  int q_index = 0;
296#endif
297
[6aee3ab]298  // ** Fill in the local values table **
299  // Storage for the current parameter values.
300  // These will be updated as we walk the dispersity mesh.
[9eb3632]301  ParameterBlock local_values;
[6aee3ab]302  //   values[0] is scale
303  //   values[1] is background
304  #ifdef USE_OPENMP
305  #pragma omp parallel for
306  #endif
307  for (int i=0; i < NUM_PARS; i++) {
308    local_values.vector[i] = values[2+i];
309    //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
310  }
311  //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
312  //if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop);
[2e44ac7]313
[6aee3ab]314  // ** Precompute magnatism values **
[a4280bd]315#if defined(MAGNETIC) && NUM_MAGNETIC>0
[bde38b5]316  // Location of the sld parameters in the parameter vector.
[9eb3632]317  // These parameters are updated with the effective sld due to magnetism.
318  const int32_t slds[] = { MAGNETIC_PARS };
[32e3c9b]319
320  // Interpret polarization cross section.
[a4280bd]321  //     up_frac_i = values[NUM_PARS+2];
322  //     up_frac_f = values[NUM_PARS+3];
323  //     up_angle = values[NUM_PARS+4];
[9ee2756]324  // TODO: could precompute more magnetism parameters before calling the kernel.
[885753a]325  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
[32e3c9b]326  double cos_mspin, sin_mspin;
[885753a]327  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights);
[a4280bd]328  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
[9eb3632]329#endif // MAGNETIC
[3044216]330
[6aee3ab]331  // ** Fill in the initial results **
[9ee2756]332  // If pd_start is zero that means that we are starting a new calculation,
333  // and must initialize the result to zero.  Otherwise, we are restarting
334  // the calculation from somewhere in the middle of the dispersity mesh,
335  // and we update the value rather than reset it. Similarly for the
336  // normalization factor, which is stored as the final value in the
337  // results vector (one past the number of q values).
338  //
339  // The code differs slightly between opencl and dll since opencl is only
340  // seeing one q value (stored in the variable "this_result") while the dll
341  // version must loop over all q.
342  #ifdef USE_OPENCL
343    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
344    double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
345  #else // !USE_OPENCL
346    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
347    if (pd_start == 0) {
348      #ifdef USE_OPENMP
349      #pragma omp parallel for
350      #endif
351      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
352    }
[6aee3ab]353    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
[9ee2756]354#endif // !USE_OPENCL
355
356
357// ====== macros to set up the parts of the loop =======
358/*
359Based on the level of the loop, uses C preprocessor magic to construct
360level-specific looping variables, including these from loop level 3:
361
362  int n3 : length of loop for mesh level 3
363  int i3 : current position in the loop for level 3, which is calculated
364       from a combination of pd_start, pd_stride[3] and pd_length[3].
365  int p3 : is the index into the parameter table for mesh level 3
366  double v3[] : pointer into dispersity array to values for loop 3
367  double w3[] : pointer into dispersity array to weights for loop 3
368  double weight3 : the product of weights from levels 3 and up, computed
369       as weight5*weight4*w3[i3].  Note that we need an outermost
370       value weight5 set to 1.0 for this to work properly.
371
372After expansion, the loop struction will look like the following:
373
374  // --- PD_INIT(4) ---
375  const int n4 = pd_length[4];
376  const int p4 = pd_par[4];
377  global const double *v4 = pd_value + pd_offset[4];
378  global const double *w4 = pd_weight + pd_offset[4];
379  int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start
380
381  // --- PD_INIT(3) ---
382  const int n3 = pd_length[3];
383  ...
384  int i3 = (pd_start/pd_stride[3])%n3;  // position in level 3 at pd_start
385
386  PD_INIT(2)
387  PD_INIT(1)
388  PD_INIT(0)
389
390  // --- PD_OUTERMOST_WEIGHT(5) ---
391  const double weight5 = 1.0;
392
393  // --- PD_OPEN(4,5) ---
394  while (i4 < n4) {
395    parameter[p4] = v4[i4];  // set the value for pd parameter 4 at this mesh point
396    const double weight4 = w4[i4] * weight5;
397
398    // from PD_OPEN(3,4)
399    while (i3 < n3) {
400      parameter[p3] = v3[i3];  // set the value for pd parameter 3 at this mesh point
401      const double weight3 = w3[i3] * weight4;
402
403      PD_OPEN(3,2)
404      PD_OPEN(2,1)
405      PD_OPEN(0,1)
406
[6aee3ab]407      // ... main loop body ...
408      APPLY_PROJECTION    // convert jitter values to spherical coords
409      BUILD_ROTATION      // construct the rotation matrix qxy => qabc
410      for each q
411          FETCH_Q         // set qx,qy from the q input vector
412          APPLY_ROTATION  // convert qx,qy to qa,qb,qc
413          CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
[9ee2756]414
415      ++step;  // increment counter representing position in dispersity mesh
416
417      PD_CLOSE(0)
418      PD_CLOSE(1)
419      PD_CLOSE(2)
420
421      // --- PD_CLOSE(3) ---
422      if (step >= pd_stop) break;
423      ++i3;
424    }
425    i3 = 0; // reset loop counter for next round through the loop
426
427    // --- PD_CLOSE(4) ---
428    if (step >= pd_stop) break;
429    ++i4;
[2e44ac7]430  }
[9ee2756]431  i4 = 0; // reset loop counter even though no more rounds through the loop
432
433*/
434
435
[6aee3ab]436// ** prepare inner loops **
[9ee2756]437
438// Depending on the shape type (radial, axial, triaxial), the variables
[6aee3ab]439// and calling parameters in the loop body will be slightly different.
440// Macros capture the differences in one spot so the rest of the code
441// is easier to read. The code below both declares variables for the
442// inner loop and defines the macros that use them.
[ff10479]443
[9ee2756]444#if defined(CALL_IQ)
445  // unoriented 1D
446  double qk;
[2c108a3]447  #define FETCH_Q() do { qk = q[q_index]; } while (0)
448  #define BUILD_ROTATION() do {} while(0)
449  #define APPLY_ROTATION() do {} while(0)
450  #define CALL_KERNEL() CALL_IQ(qk, local_values.table)
[9ee2756]451
452#elif defined(CALL_IQ_A)
453  // unoriented 2D
454  double qx, qy;
[2c108a3]455  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
456  #define BUILD_ROTATION() do {} while(0)
457  #define APPLY_ROTATION() do {} while(0)
458  #define CALL_KERNEL() CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table)
[9ee2756]459
460#elif defined(CALL_IQ_AC)
461  // oriented symmetric 2D
462  double qx, qy;
[2c108a3]463  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
[9ee2756]464  double qa, qc;
465  QACRotation rotation;
[ff10479]466  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
[767dca8]467  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
[ec8d4ac]468  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc)
[2c108a3]469  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
[9ee2756]470
471#elif defined(CALL_IQ_ABC)
472  // oriented asymmetric 2D
473  double qx, qy;
[2c108a3]474  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
[9ee2756]475  double qa, qb, qc;
476  QABCRotation rotation;
[ff10479]477  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
478  // psi and dpsi are only for IQ_ABC, so they are processed here.
479  const double psi = values[details->theta_par+4];
[6aee3ab]480  local_values.table.psi = 0.;
[ff10479]481  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
[ec8d4ac]482  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
[ff10479]483  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
[108e70e]484#elif defined(CALL_IQ_XY)
485  // direct call to qx,qy calculator
486  double qx, qy;
487  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
488  #define BUILD_ROTATION() do {} while(0)
489  #define APPLY_ROTATION() do {} while(0)
490  #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table)
[ff10479]491#endif
492
[2a7e20e]493// Define APPLY_PROJECTION depending on model symmetries. We do this outside
494// the previous if block so that we don't need to repeat the identical
495// logic in the IQ_AC and IQ_ABC branches.  This will become more important
496// if we implement more projections, or more complicated projections.
497#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
[ff10479]498  #define APPLY_PROJECTION() const double weight=weight0
[2a7e20e]499#elif defined(CALL_IQ_XY) // pass orientation to the model
[108e70e]500  // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc
501  // Need to plug the values for the orientation angles back into parameter
502  // table in case they were overridden by the orientation offset.  This
503  // means that orientation dispersity will not work for these models, but
504  // it was broken anyway, so no matter.  Still want to provide Iqxy in case
505  // the user model wants full control of orientation/magnetism.
506  #if defined(HAVE_PSI)
507    const double theta = values[details->theta_par+2];
508    const double phi = values[details->theta_par+3];
509    const double psi = values[details->theta_par+4];
510    double weight;
511    #define APPLY_PROJECTION() do { \
512      local_values.table.theta = theta; \
513      local_values.table.phi = phi; \
514      local_values.table.psi = psi; \
515      weight=weight0; \
516    } while (0)
517  #elif defined(HAVE_THETA)
518    const double theta = values[details->theta_par+2];
519    const double phi = values[details->theta_par+3];
520    double weight;
521    #define APPLY_PROJECTION() do { \
522      local_values.table.theta = theta; \
523      local_values.table.phi = phi; \
524      weight=weight0; \
525    } while (0)
526  #else
527    #define APPLY_PROJECTION() const double weight=weight0
528  #endif
[2a7e20e]529#else // apply jitter and view before calling the model
[9ee2756]530  // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
531  const double theta = values[details->theta_par+2];
532  const double phi = values[details->theta_par+3];
[6aee3ab]533  // Make sure jitter angle defaults to zero if there is no jitter distribution
534  local_values.table.theta = 0.;
535  local_values.table.phi = 0.;
[ff10479]536  // The "jitter" angles (dtheta, dphi, dpsi) are stored with the
537  // dispersity values and copied to the local parameter table as
538  // we go through the mesh.
539  double dtheta, dphi, weight;
[2a7e20e]540  #if PROJECTION == 1 // equirectangular
[ff10479]541    #define APPLY_PROJECTION() do { \
[767dca8]542      dtheta = local_values.table.theta; \
543      dphi = local_values.table.phi; \
[ff10479]544      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
[767dca8]545    } while (0)
[2a7e20e]546  #elif PROJECTION == 2 // sinusoidal
[ff10479]547    #define APPLY_PROJECTION() do { \
[767dca8]548      dtheta = local_values.table.theta; \
549      dphi = local_values.table.phi; \
[ff10479]550      weight = weight0; \
551      if (dtheta != 90.0) dphi /= cos(dtheta*M_PI_180); \
552      else if (dphi != 0.0) weight = 0.; \
553      if (fabs(dphi) >= 180.) weight = 0.; \
[767dca8]554    } while (0)
555  #endif
[2a7e20e]556#endif // done defining APPLY_PROJECTION
[9ee2756]557
[6aee3ab]558// ** define looping macros **
559
560// Define looping variables
561#define PD_INIT(_LOOP) \
562  const int n##_LOOP = details->pd_length[_LOOP]; \
563  const int p##_LOOP = details->pd_par[_LOOP]; \
564  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
565  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
566  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
567
568// Jump into the middle of the dispersity loop
569#define PD_OPEN(_LOOP,_OUTER) \
570  while (i##_LOOP < n##_LOOP) { \
571    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \
572    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER;
573
574// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD).
575#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0;
576#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n)
577
578// Close out the loop
579#define PD_CLOSE(_LOOP) \
580    if (step >= pd_stop) break; \
581    ++i##_LOOP; \
582  } \
583  i##_LOOP = 0;
[9eb3632]584
[9ee2756]585// ====== construct the loops =======
586
587// Pointers to the start of the dispersity and weight vectors, if needed.
[7b7da6b]588#if MAX_PD>0
[bde38b5]589  global const double *pd_value = values + NUM_VALUES;
590  global const double *pd_weight = pd_value + details->num_weights;
[7b7da6b]591#endif
[9eb3632]592
[9ee2756]593// The variable "step" is the current position in the dispersity loop.
594// It will be incremented each time a new point in the mesh is accumulated,
595// and used to test whether we have reached pd_stop.
596int step = pd_start;
597
[6aee3ab]598// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 ***
599
[9ee2756]600// define looping variables
[9eb3632]601#if MAX_PD>4
[9ee2756]602  PD_INIT(4)
[9eb3632]603#endif
604#if MAX_PD>3
[9ee2756]605  PD_INIT(3)
[9eb3632]606#endif
607#if MAX_PD>2
[9ee2756]608  PD_INIT(2)
[9eb3632]609#endif
610#if MAX_PD>1
[9ee2756]611  PD_INIT(1)
[9eb3632]612#endif
613#if MAX_PD>0
[9ee2756]614  PD_INIT(0)
[9eb3632]615#endif
[2e44ac7]616
[9ee2756]617// open nested loops
618PD_OUTERMOST_WEIGHT(MAX_PD)
[9eb3632]619#if MAX_PD>4
[9ee2756]620  PD_OPEN(4,5)
[9eb3632]621#endif
622#if MAX_PD>3
[9ee2756]623  PD_OPEN(3,4)
[9eb3632]624#endif
625#if MAX_PD>2
[9ee2756]626  PD_OPEN(2,3)
[9eb3632]627#endif
628#if MAX_PD>1
[9ee2756]629  PD_OPEN(1,2)
[9eb3632]630#endif
631#if MAX_PD>0
[9ee2756]632  PD_OPEN(0,1)
[9eb3632]633#endif
[5ff1b03]634
[9ee2756]635//if (q_index==0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n");}
636
637  // ====== loop body =======
638  #ifdef INVALID
639  if (!INVALID(local_values.table))
640  #endif
641  {
[ff10479]642     APPLY_PROJECTION();
643
[9ee2756]644    // Accumulate I(q)
645    // Note: weight==0 must always be excluded
[767dca8]646    if (weight > cutoff) {
647      pd_norm += weight * CALL_VOLUME(local_values.table);
[2c108a3]648      BUILD_ROTATION();
[9ee2756]649
650#ifndef USE_OPENCL
651      // DLL needs to explicitly loop over the q values.
652      #ifdef USE_OPENMP
653      #pragma omp parallel for
654      #endif
655      for (q_index=0; q_index<nq; q_index++)
656#endif // !USE_OPENCL
657      {
658
[2c108a3]659        FETCH_Q();
660        APPLY_ROTATION();
[9ee2756]661
662        // ======= COMPUTE SCATTERING ==========
663        #if defined(MAGNETIC) && NUM_MAGNETIC > 0
[2c108a3]664          // Compute the scattering from the magnetic cross sections.
665          double scattering = 0.0;
666          const double qsq = qx*qx + qy*qy;
667          if (qsq > 1.e-16) {
668            // TODO: what is the magnetic scattering at q=0
669            const double px = (qy*cos_mspin + qx*sin_mspin)/qsq;
670            const double py = (qy*sin_mspin - qx*cos_mspin)/qsq;
671
672            // loop over uu, ud real, du real, dd, ud imag, du imag
[aadec17]673            for (unsigned int xs=0; xs<6; xs++) {
[885753a]674              const double xs_weight = xs_weights[xs];
[2c108a3]675              if (xs_weight > 1.e-8) {
676                // Since the cross section weight is significant, set the slds
677                // to the effective slds for this cross section, call the
678                // kernel, and add according to weight.
679                for (int sk=0; sk<NUM_MAGNETIC; sk++) {
680                  const int32_t mag_index = NUM_PARS+5 + 3*sk;
681                  const int32_t sld_index = slds[sk];
682                  const double mx = values[mag_index];
683                  const double my = values[mag_index+1];
684                  const double mz = values[mag_index+2];
685                  local_values.vector[sld_index] =
686                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz);
[885753a]687//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n",
[dc6f601]688//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz);
[2c108a3]689                }
690                scattering += xs_weight * CALL_KERNEL();
[a4280bd]691              }
[9eb3632]692            }
[32e3c9b]693          }
[9ee2756]694        #else  // !MAGNETIC
[2c108a3]695          const double scattering = CALL_KERNEL();
[9ee2756]696        #endif // !MAGNETIC
[767dca8]697//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
[9ee2756]698
699        #ifdef USE_OPENCL
[767dca8]700          this_result += weight * scattering;
[9ee2756]701        #else // !USE_OPENCL
[767dca8]702          result[q_index] += weight * scattering;
[9ee2756]703        #endif // !USE_OPENCL
[3044216]704      }
[03cac08]705    }
[2e44ac7]706  }
[9ee2756]707
708// close nested loops
709++step;
710#if MAX_PD>0
711  PD_CLOSE(0)
[9eb3632]712#endif
713#if MAX_PD>1
[9ee2756]714  PD_CLOSE(1)
[9eb3632]715#endif
716#if MAX_PD>2
[9ee2756]717  PD_CLOSE(2)
[9eb3632]718#endif
719#if MAX_PD>3
[9ee2756]720  PD_CLOSE(3)
[9eb3632]721#endif
722#if MAX_PD>4
[9ee2756]723  PD_CLOSE(4)
[9eb3632]724#endif
[f2f67a6]725
[9ee2756]726// Remember the current result and the updated norm.
727#ifdef USE_OPENCL
728  result[q_index] = this_result;
729  if (q_index == 0) result[nq] = pd_norm;
730//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
731#else // !USE_OPENCL
[a738209]732  result[nq] = pd_norm;
[9ee2756]733//printf("res: %g/%g\n", result[0], pd_norm);
734#endif // !USE_OPENCL
[2c108a3]735
[6aee3ab]736// ** clear the macros in preparation for the next kernel **
[2c108a3]737#undef PD_INIT
738#undef PD_OPEN
739#undef PD_CLOSE
740#undef FETCH_Q
[ff10479]741#undef APPLY_PROJECTION
[2c108a3]742#undef BUILD_ROTATION
743#undef APPLY_ROTATION
744#undef CALL_KERNEL
[2e44ac7]745}
Note: See TracBrowser for help on using the repository browser.