source: sasmodels/sasmodels/kernel_iq.c @ a34b811

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since a34b811 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

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