source: sasmodels/sasmodels/kernel_iq.c @ 94f4543

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 94f4543 was 94f4543, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

update asymmetric orientation equations to be consistent between docs, jitter viewer and scattering pattern

  • Property mode set to 100644
File size: 17.0 KB
Line 
1
2/*
3    ##########################################################
4    #                                                        #
5    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
6    #   !!                                              !!   #
7    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
8    #   !!                                              !!   #
9    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
10    #                                                        #
11    ##########################################################
12*/
13
14#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
15#define _PAR_BLOCK_
16
17typedef struct {
18#if MAX_PD > 0
19    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
20    int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
21    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
22    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
23#endif // MAX_PD > 0
24    int32_t num_eval;           // total number of voxels in hypercube
25    int32_t num_weights;        // total length of the weights vector
26    int32_t num_active;         // number of non-trivial pd loops
27    int32_t theta_par;          // id of first orientation variable
28} ProblemDetails;
29
30// Intel HD 4000 needs private arrays to be a multiple of 4 long
31typedef struct {
32    PARAMETER_TABLE
33} ParameterTable;
34typedef union {
35    ParameterTable table;
36    double vector[4*((NUM_PARS+3)/4)];
37} ParameterBlock;
38#endif // _PAR_BLOCK_
39
40
41#if defined(MAGNETIC) && NUM_MAGNETIC>0
42
43// Return value restricted between low and high
44static double clip(double value, double low, double high)
45{
46  return (value < low ? low : (value > high ? high : value));
47}
48
49// Compute spin cross sections given in_spin and out_spin
50// To convert spin cross sections to sld b:
51//     uu * (sld - m_sigma_x);
52//     dd * (sld + m_sigma_x);
53//     ud * (m_sigma_y + 1j*m_sigma_z);
54//     du * (m_sigma_y - 1j*m_sigma_z);
55static void set_spins(double in_spin, double out_spin, double spins[4])
56{
57  in_spin = clip(in_spin, 0.0, 1.0);
58  out_spin = clip(out_spin, 0.0, 1.0);
59  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
60  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du
61  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud
62  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
63}
64
65static double mag_sld(double qx, double qy, double p,
66                       double mx, double my, double sld)
67{
68    const double perp = qy*mx - qx*my;
69    return sld + perp*p;
70}
71//#endif // MAGNETIC
72
73// TODO: way too hackish
74// For the 1D kernel, CALL_IQ_[A,AC,ABC] and MAGNETIC are not defined
75// so view_direct *IS NOT* included
76// For the 2D kernel, CALL_IQ_[A,AC,ABC] is defined but MAGNETIC is not
77// so view_direct *IS* included
78// For the magnetic kernel, CALL_IQ_[A,AC,ABC] is defined, but so is MAGNETIC
79// so view_direct *IS NOT* included
80#else // !MAGNETIC
81
82// ===== Implement jitter in orientation =====
83// To change the definition of the angles, run explore/angles.py, which
84// uses sympy to generate the equations.
85
86#if defined(CALL_IQ_AC) // oriented symmetric
87static double
88view_direct(double qx, double qy,
89             double theta, double phi,
90             ParameterTable table)
91{
92    double sin_theta, cos_theta;
93    double sin_phi, cos_phi;
94
95    // reverse view
96    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
97    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
98    const double qa = qx*cos_phi*cos_theta + qy*sin_phi*cos_theta;
99    const double qb = -qx*sin_phi + qy*cos_phi;
100    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta;
101
102    // reverse jitter after view
103    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta);
104    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi);
105    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta;
106
107    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
108    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
109
110    return CALL_IQ_AC(dqa, dqc, table);
111}
112
113#elif defined(CALL_IQ_ABC) // oriented asymmetric
114
115static double
116view_direct(double qx, double qy,
117             double theta, double phi, double psi,
118             ParameterTable table)
119{
120    double sin_theta, cos_theta;
121    double sin_phi, cos_phi;
122    double sin_psi, cos_psi;
123
124    // reverse view
125    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
126    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
127    SINCOS(psi*M_PI_180, sin_psi, cos_psi);
128    const double qa = qx*(-sin_phi*sin_psi + cos_phi*cos_psi*cos_theta) + qy*(sin_phi*cos_psi*cos_theta + sin_psi*cos_phi);
129    const double qb = qx*(-sin_phi*cos_psi - sin_psi*cos_phi*cos_theta) + qy*(-sin_phi*sin_psi*cos_theta + cos_phi*cos_psi);
130    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta;
131
132    // reverse jitter after view
133    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta);
134    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi);
135    SINCOS(table.psi*M_PI_180, sin_psi, cos_psi);
136    const double dqa = qa*cos_psi*cos_theta + qb*(sin_phi*sin_theta*cos_psi + sin_psi*cos_phi) + qc*(sin_phi*sin_psi - sin_theta*cos_phi*cos_psi);
137    const double dqb = -qa*sin_psi*cos_theta + qb*(-sin_phi*sin_psi*sin_theta + cos_phi*cos_psi) + qc*(sin_phi*cos_psi + sin_psi*sin_theta*cos_phi);
138    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta;
139
140    return CALL_IQ_ABC(dqa, dqb, dqc, table);
141}
142/* TODO: use precalculated jitter for faster 2D calcs on DLL.
143static void
144view_precalc(
145    double theta, double phi, double psi,
146    ParameterTable table,
147    double *R11, double *R12, double *R21,
148    double *R22, double *R31, double *R32)
149{
150    double sin_theta, cos_theta;
151    double sin_phi, cos_phi;
152    double sin_psi, cos_psi;
153
154    // reverse view matrix
155    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
156    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
157    SINCOS(psi*M_PI_180, sin_psi, cos_psi);
158    const double V11 = sin_phi*sin_psi + cos_phi*cos_psi*cos_theta;
159    const double V12 = sin_phi*cos_psi*cos_theta - sin_psi*cos_phi;
160    const double V21 = -sin_phi*cos_psi + sin_psi*cos_phi*cos_theta;
161    const double V22 = sin_phi*sin_psi*cos_theta + cos_phi*cos_psi;
162    const double V31 = sin_theta*cos_phi;
163    const double V32 = sin_phi*sin_theta;
164
165    // reverse jitter matrix
166    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta);
167    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi);
168    SINCOS(table.psi*M_PI_180, sin_psi, cos_psi);
169    const double J11 = cos_psi*cos_theta;
170    const double J12 = sin_phi*sin_theta*cos_psi - sin_psi*cos_phi;
171    const double J13 = -sin_phi*sin_psi - sin_theta*cos_phi*cos_psi;
172    const double J21 = sin_psi*cos_theta;
173    const double J22 = sin_phi*sin_psi*sin_theta + cos_phi*cos_psi;
174    const double J23 = sin_phi*cos_psi - sin_psi*sin_theta*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    *R11 = J11*V11 + J12*V21 + J13*V31;
181    *R12 = J11*V12 + J12*V22 + J13*V32;
182    *R21 = J21*V11 + J22*V21 + J23*V31;
183    *R22 = J21*V12 + J22*V22 + J23*V32;
184    *R31 = J31*V11 + J32*V21 + J33*V31;
185    *R32 = J31*V12 + J32*V22 + J33*V32;
186}
187
188static double
189view_apply(double qx, double qy,
190    double R11, double R12, double R21, double R22, double R31, double R32,
191    ParameterTable table)
192{
193    const double dqa = R11*qx + R12*qy;
194    const double dqb = R21*qx + R22*qy;
195    const double dqc = R31*qx + R32*qy;
196
197    CALL_IQ_ABC(dqa, dqb, dqc, table);
198}
199*/
200#endif
201
202#endif // !MAGNETIC
203
204kernel
205void KERNEL_NAME(
206    int32_t nq,                 // number of q values
207    const int32_t pd_start,     // where we are in the polydispersity loop
208    const int32_t pd_stop,      // where we are stopping in the polydispersity loop
209    global const ProblemDetails *details,
210    global const double *values,
211    global const double *q, // nq q values, with padding to boundary
212    global double *result,  // nq+1 return values, again with padding
213    const double cutoff     // cutoff in the polydispersity weight product
214    )
215{
216  // Storage for the current parameter values.  These will be updated as we
217  // walk the polydispersity cube.
218  ParameterBlock local_values;
219
220#if defined(MAGNETIC) && NUM_MAGNETIC>0
221  // Location of the sld parameters in the parameter vector.
222  // These parameters are updated with the effective sld due to magnetism.
223  #if NUM_MAGNETIC > 3
224  const int32_t slds[] = { MAGNETIC_PARS };
225  #endif
226
227  // TODO: could precompute these outside of the kernel.
228  // Interpret polarization cross section.
229  //     up_frac_i = values[NUM_PARS+2];
230  //     up_frac_f = values[NUM_PARS+3];
231  //     up_angle = values[NUM_PARS+4];
232  double spins[4];
233  double cos_mspin, sin_mspin;
234  set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins);
235  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
236#endif // MAGNETIC
237
238#if defined(CALL_IQ_AC) // oriented symmetric
239  const double theta = values[details->theta_par+2];
240  const double phi = values[details->theta_par+3];
241#elif defined(CALL_IQ_ABC) // oriented asymmetric
242  const double theta = values[details->theta_par+2];
243  const double phi = values[details->theta_par+3];
244  const double psi = values[details->theta_par+4];
245#endif
246
247  // Fill in the initial variables
248  //   values[0] is scale
249  //   values[1] is background
250  #ifdef USE_OPENMP
251  #pragma omp parallel for
252  #endif
253  for (int i=0; i < NUM_PARS; i++) {
254    local_values.vector[i] = values[2+i];
255//printf("p%d = %g\n",i, local_values.vector[i]);
256  }
257//printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
258//printf("start:%d  stop:%d\n", pd_start, pd_stop);
259
260  double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
261  if (pd_start == 0) {
262    #ifdef USE_OPENMP
263    #pragma omp parallel for
264    #endif
265    for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
266  }
267//printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
268
269#if MAX_PD>0
270  global const double *pd_value = values + NUM_VALUES;
271  global const double *pd_weight = pd_value + details->num_weights;
272#endif
273
274  // Jump into the middle of the polydispersity loop
275#if MAX_PD>4
276  int n4=details->pd_length[4];
277  int i4=(pd_start/details->pd_stride[4])%n4;
278  const int p4=details->pd_par[4];
279  global const double *v4 = pd_value + details->pd_offset[4];
280  global const double *w4 = pd_weight + details->pd_offset[4];
281#endif
282#if MAX_PD>3
283  int n3=details->pd_length[3];
284  int i3=(pd_start/details->pd_stride[3])%n3;
285  const int p3=details->pd_par[3];
286  global const double *v3 = pd_value + details->pd_offset[3];
287  global const double *w3 = pd_weight + details->pd_offset[3];
288//printf("offset %d: %d %d\n", 3, details->pd_offset[3], NUM_VALUES);
289#endif
290#if MAX_PD>2
291  int n2=details->pd_length[2];
292  int i2=(pd_start/details->pd_stride[2])%n2;
293  const int p2=details->pd_par[2];
294  global const double *v2 = pd_value + details->pd_offset[2];
295  global const double *w2 = pd_weight + details->pd_offset[2];
296#endif
297#if MAX_PD>1
298  int n1=details->pd_length[1];
299  int i1=(pd_start/details->pd_stride[1])%n1;
300  const int p1=details->pd_par[1];
301  global const double *v1 = pd_value + details->pd_offset[1];
302  global const double *w1 = pd_weight + details->pd_offset[1];
303#endif
304#if MAX_PD>0
305  int n0=details->pd_length[0];
306  int i0=(pd_start/details->pd_stride[0])%n0;
307  const int p0=details->pd_par[0];
308  global const double *v0 = pd_value + details->pd_offset[0];
309  global const double *w0 = pd_weight + details->pd_offset[0];
310//printf("w0:%p, values:%p, diff:%ld, %d\n",w0,values,(w0-values), NUM_VALUES);
311#endif
312
313
314  int step = pd_start;
315
316#if MAX_PD>4
317  const double weight5 = 1.0;
318  while (i4 < n4) {
319    local_values.vector[p4] = v4[i4];
320    double weight4 = w4[i4] * weight5;
321//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, local_values.vector[p4], weight4);
322#elif MAX_PD>3
323    const double weight4 = 1.0;
324#endif
325#if MAX_PD>3
326  while (i3 < n3) {
327    local_values.vector[p3] = v3[i3];
328    double weight3 = w3[i3] * weight4;
329//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, local_values.vector[p3], weight3);
330#elif MAX_PD>2
331    const double weight3 = 1.0;
332#endif
333#if MAX_PD>2
334  while (i2 < n2) {
335    local_values.vector[p2] = v2[i2];
336    double weight2 = w2[i2] * weight3;
337//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, local_values.vector[p2], weight2);
338#elif MAX_PD>1
339    const double weight2 = 1.0;
340#endif
341#if MAX_PD>1
342  while (i1 < n1) {
343    local_values.vector[p1] = v1[i1];
344    double weight1 = w1[i1] * weight2;
345//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, local_values.vector[p1], weight1);
346#elif MAX_PD>0
347    const double weight1 = 1.0;
348#endif
349#if MAX_PD>0
350  while(i0 < n0) {
351    local_values.vector[p0] = v0[i0];
352    double weight0 = w0[i0] * weight1;
353//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
354#else
355    const double weight0 = 1.0;
356#endif
357
358//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");
359//printf("sphcor: %g\n", spherical_correction);
360
361    #ifdef INVALID
362    if (!INVALID(local_values.table))
363    #endif
364    {
365      // Accumulate I(q)
366      // Note: weight==0 must always be excluded
367      if (weight0 > cutoff) {
368        pd_norm += weight0 * CALL_VOLUME(local_values.table);
369
370        #ifdef USE_OPENMP
371        #pragma omp parallel for
372        #endif
373        for (int q_index=0; q_index<nq; q_index++) {
374#if defined(MAGNETIC) && NUM_MAGNETIC > 0
375          const double qx = q[2*q_index];
376          const double qy = q[2*q_index+1];
377          const double qsq = qx*qx + qy*qy;
378
379          // Constant across orientation, polydispersity for given qx, qy
380          double scattering = 0.0;
381          // TODO: what is the magnetic scattering at q=0
382          if (qsq > 1.e-16) {
383            double p[4];  // dd, du, ud, uu
384            p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
385            p[3] = -p[0];
386            p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq;
387
388            for (int index=0; index<4; index++) {
389              const double xs = spins[index];
390              if (xs > 1.e-8) {
391                const int spin_flip = (index==1) || (index==2);
392                const double pk = p[index];
393                for (int axis=0; axis<=spin_flip; axis++) {
394                  #define M1 NUM_PARS+5
395                  #define M2 NUM_PARS+8
396                  #define M3 NUM_PARS+13
397                  #define SLD(_M_offset, _sld_offset) \
398                      local_values.vector[_sld_offset] = xs * (axis \
399                      ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
400                      : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
401                                (spin_flip ? 0.0 : values[_sld_offset+2])))
402                  #if NUM_MAGNETIC==1
403                      SLD(M1, MAGNETIC_PAR1);
404                  #elif NUM_MAGNETIC==2
405                      SLD(M1, MAGNETIC_PAR1);
406                      SLD(M2, MAGNETIC_PAR2);
407                  #elif NUM_MAGNETIC==3
408                      SLD(M1, MAGNETIC_PAR1);
409                      SLD(M2, MAGNETIC_PAR2);
410                      SLD(M3, MAGNETIC_PAR3);
411                  #else
412                  for (int sk=0; sk<NUM_MAGNETIC; sk++) {
413                      SLD(M1+3*sk, slds[sk]);
414                  }
415                  #endif
416#  if defined(CALL_IQ_A) // unoriented
417                  scattering += CALL_IQ_A(sqrt(qsq), local_values.table);
418#  elif defined(CALL_IQ_AC) // oriented symmetric
419                  scattering += view_direct(qx, qy, theta, phi, local_values.table);
420#  elif defined(CALL_IQ_ABC) // oriented asymmetric
421                  scattering += view_direct(qx, qy, theta, phi, psi, local_values.table);
422#  endif
423                }
424              }
425            }
426          }
427#elif defined(CALL_IQ) // 1d, not MAGNETIC
428          const double scattering = CALL_IQ(q[q_index], local_values.table);
429#else  // 2d data, not magnetic
430          const double qx = q[2*q_index];
431          const double qy = q[2*q_index+1];
432#  if defined(CALL_IQ_A) // unoriented
433          const double scattering = CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table);
434#  elif defined(CALL_IQ_AC) // oriented symmetric
435          const double scattering = view_direct(qx, qy, theta, phi, local_values.table);
436#  elif defined(CALL_IQ_ABC) // oriented asymmetric
437          const double scattering = view_direct(qx, qy, theta, phi, psi, local_values.table);
438#  endif
439#endif // !MAGNETIC
440//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0);
441          result[q_index] += weight0 * scattering;
442        }
443      }
444    }
445    ++step;
446#if MAX_PD>0
447    if (step >= pd_stop) break;
448    ++i0;
449  }
450  i0 = 0;
451#endif
452#if MAX_PD>1
453    if (step >= pd_stop) break;
454    ++i1;
455  }
456  i1 = 0;
457#endif
458#if MAX_PD>2
459    if (step >= pd_stop) break;
460    ++i2;
461  }
462  i2 = 0;
463#endif
464#if MAX_PD>3
465    if (step >= pd_stop) break;
466    ++i3;
467  }
468  i3 = 0;
469#endif
470#if MAX_PD>4
471    if (step >= pd_stop) break;
472    ++i4;
473  }
474  i4 = 0;
475#endif
476
477//printf("res: %g/%g\n", result[0], pd_norm);
478  // Remember the updated norm.
479  result[nq] = pd_norm;
480}
Note: See TracBrowser for help on using the repository browser.