Core code analysis: How GISAS intensity is computed

last updated 26feb26

Generic simulation:

Simulation::simulate() -> Datafield
  prepareSimulation() virtual
  re_sample = ReSample::make(..)
  m_cache = vector<double>(nOutChannels(), 0)
  for (i < n_combinations): // loop over distributions
    runSingleSimulation(re_sample, batch...)
  return packResult() // m_cache -> Datafield

Simulation::runSingleSimulation(re_sample, batch...) -> void
  for threads:
    for i_pixel in thread:
      runComputation(re_sample, i_pixel, wgt) // pure virtual

Simulate scattering intensity:

ScatteringSimulation::runComputation(re_sample, i_pixel, weight) -> void
  // Computes relection and scattering towards m_pixel[i], and stores it in m_cache[i].
  // Accounts for beam intensity, footprint, and given weight (from distributions).
  DiffuseElement ele(m_beam, m_detector, m_pixel)
  m_cache[i] += weight * ... * Evaluate::scatteredContribution(re_sample, ele, ...)

Evaluate::scatteredContribution(re_sample, ele, ...) -> double
  // implemented in Sim/Evaluate/EvaluateScattering.cpp
  return sum over s in re_sample.structs() of:
    s->coverage() * s->numberDensity() * s->scalarScattering(ele)) // or polarized
    // with (area|volume)NumberDensity() for Struct2D|3D

// For all 2D structures (Sample/Aggregate/Struct2D.cpp):
Struct2D::scalarScattering(ele) -> double
  [mean, variance] := particle()->scatteringMoments(ele)
  return variance + structureFactor(ele.meanQ()) * |mean|²
  // where:
  //   mean     = <F>              (average form factor amplitude)
  //   variance = <|F|²> - |<F>|²  (incoherent contribution)
  // structureFactor() returns 1.0 in base class Struct2D (random assemblies)
  // structureFactor() is overridden in Ordered2D subclasses

For a simple particle, variance = 0 and the result is S(q) · |F|². For mixtures or compounds, the variance accounts for interference between different particle contributions.

Scattering moments for different particle types

Instances of IParticle::scatteringMoments:

// Simple particle (Sample/Particle/Particle.cpp):
Particle::scatteringMoments(ele) -> ScatteringMoments
  return {dwbaFF(ele), 0}   // mean = F, variance = 0

// Random mixture of particles (Sample/Particle/Mixture.cpp):
Mixture::scatteringMoments(ele) -> ScatteringMoments
  F, I = 0, 0
  for i <= n_particles:
    w = relative weight of particle i
    [f, v] = particle[i]->scatteringMoments(ele)
    F += w · f
    I += w · (v + |f|²)
  return {F, I - |F|²}

// Compound particle (Sample/Particle/Compound.cpp):
Compound::scatteringMoments(ele) -> ScatteringMoments
  // Sums contributions from constituent particles (fixed relative positions)
  F, V = 0, 0
  for p in particles:
    [f, v] = p->scatteringMoments(ele)
    F += f
    V += v
  return {F, V}

Polarized scattering

For polarized scattering, amplitudes become 2×2 spin matrices:

// Sample/Aggregate/Struct2D.cpp
Struct2D::polarizedScattering(ele) -> double
  [mean_amp, variance] := particle()->polScatteringMoments(ele)
  coherent := mean_amp · polarizer · mean_amp
  amplitude_matrix := analyzer · coherent
  intensity_matrix := analyzer · (variance + coherent)
  amplitude_trace := |trace(amplitude_matrix)|
  intensity_trace := |trace(intensity_matrix)|
  return intensity_trace + (structureFactor(q) - 1) · amplitude_trace

Here polScatteringMoments returns:

  • mean = <M> (average spin matrix amplitude)
  • variance = <M · P_in · M†> - <M> · P_in · <M>†

DWBA computation

// DWBA computation (Sample/Scattering/Scatterer.cpp)
// Scatterer is base class of IParticle and Struct3D
Scatterer::dwbaFF(ele) -> complex
  // For pure Born approximation (only_outgoing):
  return theFF(wavevectors)
  // For full DWBA, compute 4 terms (T_in, R_in) x (T_out, R_out):
  return T_in * T_out * theFF(q_TT)   // direct scattering
       + R_in * T_out * theFF(q_RT)   // reflection then scattering
       + T_in * R_out * theFF(q_TR)   // scattering then reflection
       + R_in * R_out * theFF(q_RR)   // reflection, scattering, reflection

IParticle::theFF(WavevectorInfo wavevectors) -> complex
  // applies position phase and rotation
  return e^(i*q*position) * formfactor->theFF(rotated wavevectors)

Formfactor::theFF(WavevectorInfo wavevectors) -> complex
  returns formfactor(wavevectors.getQ())

Formfactor::formfactor(C3 q) -> complex, virtual = 0

Structure factor

Paracrystal2D:

Ordered2D::structureFactor(q, outer_iff)
  returns iff_no_inner(q, outer_iff)

Ordered2D::iff_no_inner(q, outer_iff)
  returns DWfactor(q) * (iff_without_dw(q) * outer_iff - 1.0) + 1.0

ICrystal2D::iff_without_dw(q)
  performs integration or simply
  returns interferenceForXi(m_lattice->rotationAngle(), q.x(), q.y())

Paracrystal2D::interferenceForXi(xi, qx, qy)
  returns interference1D(qx, qy, xi, 0)
    * interference1D(qx, qy, xi + m_lattice->latticeAngle(), 1)

RadialParacrystal:

RadialParacrystal::iff_without_dw(q)
  directly performs Hosemann computation