source: sasmodels/sasmodels/kernel_iq.c @ 70530778

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 70530778 was 70530778, checked in by Paul Kienzle <pkienzle@…>, 14 months ago

avoid cancellation error causing NaN values to appear in 2D patterns

  • Property mode set to 100644
File size: 26.1 KB
Line 
1/*
2    ##########################################################
3    #                                                        #
4    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
5    #   !!                                              !!   #
6    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
7    #   !!                                              !!   #
8    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
9    #                                                        #
10    ##########################################################
11*/
12
13// NOTE: the following macros are defined in generate.py:
14//
15//  MAX_PD : the maximum number of dispersity loops allowed for this model,
16//      which will be at most modelinfo.MAX_PD.
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.
24//  MAGNETIC : defined when the magnetic kernel is being instantiated
25//  NUM_MAGNETIC : the number of magnetic parameters
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
33//  CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models
34//  INVALID(table) : test if the current point is feesible to calculate.  This
35//      will be defined in the kernel definition file.
36//  PROJECTION : equirectangular=1, sinusoidal=2
37//      see explore/jitter.py for definitions.
38
39#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
40#define _PAR_BLOCK_
41
42typedef struct {
43#if MAX_PD > 0
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
46    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
47    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
48#endif // MAX_PD > 0
49    int32_t num_eval;           // total number of voxels in hypercube
50    int32_t num_weights;        // total length of the weights vector
51    int32_t num_active;         // number of non-trivial pd loops
52    int32_t theta_par;          // id of first orientation variable
53} ProblemDetails;
54
55// Intel HD 4000 needs private arrays to be a multiple of 4 long
56typedef struct {
57    PARAMETER_TABLE
58} ParameterTable;
59typedef union {
60    ParameterTable table;
61    double vector[4*((NUM_PARS+3)/4)];
62} ParameterBlock;
63#endif // _PAR_BLOCK_
64
65#if defined(MAGNETIC) && NUM_MAGNETIC > 0
66// ===== Helper functions for magnetism =====
67
68// Return value restricted between low and high
69static double clip(double value, double low, double high)
70{
71  return (value < low ? low : (value > high ? high : value));
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);
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
81static void set_spin_weights(double in_spin, double out_spin, double weight[6])
82{
83  in_spin = clip(in_spin, 0.0, 1.0);
84  out_spin = clip(out_spin, 0.0, 1.0);
85  // Previous version of this function took the square root of the weights,
86  // under the assumption that
87  //
88  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
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.
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
100  weight[4] = weight[1]; // du.imag
101  weight[5] = weight[2]; // ud.imag
102}
103
104// Compute the magnetic sld
105static double mag_sld(
106  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag
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)
112{
113  if (xs < 4) {
114    const double perp = qy*mx - qx*my;
115    switch (xs) {
116      default: // keep compiler happy; condition ensures xs in [0,1,2,3]
117      case 0: // uu => sld - D M_perpx
118          return sld - px*perp;
119      case 1: // ud.real => -D M_perpy
120          return py*perp;
121      case 2: // du.real => -D M_perpy
122          return py*perp;
123      case 3: // dd => sld + D M_perpx
124          return sld + px*perp;
125    }
126  } else {
127    if (xs== 4) {
128      return -mz;  // du.imag => +D M_perpz
129    } else { // index == 5
130      return +mz;  // ud.imag => -D M_perpz
131    }
132  }
133}
134
135
136#endif
137
138// ===== Helper functions for orientation and jitter =====
139
140// To change the definition of the angles, run explore/angles.py, which
141// uses sympy to generate the equations.
142
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)
158{
159    double sin_theta, cos_theta;
160    double sin_phi, cos_phi;
161
162    // reverse view matrix
163    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
164    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
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;
171
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;
178
179    // reverse matrix
180    rotation->R31 = J31*V11 + J32*V21 + J33*V31;
181    rotation->R32 = J31*V12 + J32*V22 + J33*V32;
182}
183
184// Apply the rotation matrix returned from qac_rotation to the point (qx,qy),
185// returning R*[qx,qy]' = [qa,qc]'
186static void
187qac_apply(
188    QACRotation *rotation,
189    double qx, double qy,
190    double *qab_out, double *qc_out)
191{
192    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
193    const double dqc = rotation->R31*qx + rotation->R32*qy;
194    const double dqab_sq = -dqc*dqc + qx*qx + qy*qy;
195    //*qab_out = sqrt(fabs(dqab_sq));
196    *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0;
197    *qc_out = dqc;
198}
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]'
213static void
214qabc_rotation(
215    QABCRotation *rotation,
216    double theta, double phi, double psi,
217    double dtheta, double dphi, double dpsi)
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);
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;
231    const double V31 = sin_theta*cos_phi;
232    const double V32 = sin_phi*sin_theta;
233
234    // reverse jitter matrix
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);
238    const double J11 = cos_psi*cos_theta;
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;
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
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;
255}
256
257// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy),
258// returning R*[qx,qy]' = [qa,qb,qc]'
259static void
260qabc_apply(
261    QABCRotation *rotation,
262    double qx, double qy,
263    double *qa_out, double *qb_out, double *qc_out)
264{
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;
268}
269
270#endif // _QABC_SECTION
271
272
273// ==================== KERNEL CODE ========================
274
275kernel
276void KERNEL_NAME(
277    int32_t nq,                 // number of q values
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
280    global const ProblemDetails *details,
281    global const double *values,
282    global const double *q, // nq q values, with padding to boundary
283    global double *result,  // nq+1 return values, again with padding
284    const double cutoff     // cutoff in the dispersity weight product
285    )
286{
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
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.
301  ParameterBlock local_values;
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);
313
314  // ** Precompute magnatism values **
315#if defined(MAGNETIC) && NUM_MAGNETIC>0
316  // Location of the sld parameters in the parameter vector.
317  // These parameters are updated with the effective sld due to magnetism.
318  const int32_t slds[] = { MAGNETIC_PARS };
319
320  // Interpret polarization cross section.
321  //     up_frac_i = values[NUM_PARS+2];
322  //     up_frac_f = values[NUM_PARS+3];
323  //     up_angle = values[NUM_PARS+4];
324  // TODO: could precompute more magnetism parameters before calling the kernel.
325  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
326  double cos_mspin, sin_mspin;
327  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights);
328  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
329#endif // MAGNETIC
330
331  // ** Fill in the initial results **
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    }
353    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
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
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, ...)
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;
430  }
431  i4 = 0; // reset loop counter even though no more rounds through the loop
432
433*/
434
435
436// ** prepare inner loops **
437
438// Depending on the shape type (radial, axial, triaxial), the variables
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.
443
444#if defined(CALL_IQ)
445  // unoriented 1D
446  double qk;
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)
451
452#elif defined(CALL_IQ_A)
453  // unoriented 2D
454  double qx, qy;
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)
459
460#elif defined(CALL_IQ_AC)
461  // oriented symmetric 2D
462  double qx, qy;
463  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
464  double qa, qc;
465  QACRotation rotation;
466  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
467  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
468  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc)
469  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
470
471#elif defined(CALL_IQ_ABC)
472  // oriented asymmetric 2D
473  double qx, qy;
474  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
475  double qa, qb, qc;
476  QABCRotation rotation;
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];
480  local_values.table.psi = 0.;
481  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
482  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
483  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
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)
491#endif
492
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
498  #define APPLY_PROJECTION() const double weight=weight0
499#elif defined(CALL_IQ_XY) // pass orientation to the model
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
529#else // apply jitter and view before calling the model
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];
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.;
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;
540  #if PROJECTION == 1 // equirectangular
541    #define APPLY_PROJECTION() do { \
542      dtheta = local_values.table.theta; \
543      dphi = local_values.table.phi; \
544      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
545    } while (0)
546  #elif PROJECTION == 2 // sinusoidal
547    #define APPLY_PROJECTION() do { \
548      dtheta = local_values.table.theta; \
549      dphi = local_values.table.phi; \
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.; \
554    } while (0)
555  #endif
556#endif // done defining APPLY_PROJECTION
557
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;
584
585// ====== construct the loops =======
586
587// Pointers to the start of the dispersity and weight vectors, if needed.
588#if MAX_PD>0
589  global const double *pd_value = values + NUM_VALUES;
590  global const double *pd_weight = pd_value + details->num_weights;
591#endif
592
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
598// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 ***
599
600// define looping variables
601#if MAX_PD>4
602  PD_INIT(4)
603#endif
604#if MAX_PD>3
605  PD_INIT(3)
606#endif
607#if MAX_PD>2
608  PD_INIT(2)
609#endif
610#if MAX_PD>1
611  PD_INIT(1)
612#endif
613#if MAX_PD>0
614  PD_INIT(0)
615#endif
616
617// open nested loops
618PD_OUTERMOST_WEIGHT(MAX_PD)
619#if MAX_PD>4
620  PD_OPEN(4,5)
621#endif
622#if MAX_PD>3
623  PD_OPEN(3,4)
624#endif
625#if MAX_PD>2
626  PD_OPEN(2,3)
627#endif
628#if MAX_PD>1
629  PD_OPEN(1,2)
630#endif
631#if MAX_PD>0
632  PD_OPEN(0,1)
633#endif
634
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  {
642     APPLY_PROJECTION();
643
644    // Accumulate I(q)
645    // Note: weight==0 must always be excluded
646    if (weight > cutoff) {
647      pd_norm += weight * CALL_VOLUME(local_values.table);
648      BUILD_ROTATION();
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
659        FETCH_Q();
660        APPLY_ROTATION();
661
662        // ======= COMPUTE SCATTERING ==========
663        #if defined(MAGNETIC) && NUM_MAGNETIC > 0
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
673            for (unsigned int xs=0; xs<6; xs++) {
674              const double xs_weight = xs_weights[xs];
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);
687//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n",
688//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz);
689                }
690                scattering += xs_weight * CALL_KERNEL();
691              }
692            }
693          }
694        #else  // !MAGNETIC
695          const double scattering = CALL_KERNEL();
696        #endif // !MAGNETIC
697//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
698
699        #ifdef USE_OPENCL
700          this_result += weight * scattering;
701        #else // !USE_OPENCL
702          result[q_index] += weight * scattering;
703        #endif // !USE_OPENCL
704      }
705    }
706  }
707
708// close nested loops
709++step;
710#if MAX_PD>0
711  PD_CLOSE(0)
712#endif
713#if MAX_PD>1
714  PD_CLOSE(1)
715#endif
716#if MAX_PD>2
717  PD_CLOSE(2)
718#endif
719#if MAX_PD>3
720  PD_CLOSE(3)
721#endif
722#if MAX_PD>4
723  PD_CLOSE(4)
724#endif
725
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
732  result[nq] = pd_norm;
733//printf("res: %g/%g\n", result[0], pd_norm);
734#endif // !USE_OPENCL
735
736// ** clear the macros in preparation for the next kernel **
737#undef PD_INIT
738#undef PD_OPEN
739#undef PD_CLOSE
740#undef FETCH_Q
741#undef APPLY_PROJECTION
742#undef BUILD_ROTATION
743#undef APPLY_ROTATION
744#undef CALL_KERNEL
745}
Note: See TracBrowser for help on using the repository browser.