/*******************************************************************************
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: SAXSShells
*
* %Identification
* Written by: Martin Cramer Pedersen (mcpe@nbi.dk)
* Date: May 11, 2012
* Origin: KU-Science
* Release: McXtrace 1.0
*
* A sample of monodisperse shell-like particles in solution.
*
* %Description
* A simple component simulating the scattering from a box-shaped, thin solution
* of monodisperse, shell-like particles.
*
* Example: Sample1 = SAXSShells( xwidth = 0.01, yheight = 0.01, zdepth = 0.01, SampleToDetectorDistance = 0.5, DetectorRadius = 0.1, R = 50.0, Thickness = 20.0 )
*
* %Parameters
* R: [AA]          Average radius of the particles.
* Thickness: [AA]  Thickness of the shell - so that the outer radius is R + Thickness and the inner is R - Thickness.
* Concentration: [mM]  Concentration of sample.
* DeltaRho: [cm/AA^3]  Excess scattering length density of the particles.
* AbsorptionCrosssection: [1/m]  Absorption cross section of the sample.
* xwidth: [m]      Dimension of component in the x-direction.
* yheight: [m]     Dimension of component in the y-direction.
* zdepth: [m]      Dimension of component in the z-direction.
* SampleToDetectorDistance: [m]   Distance from sample to detector (for focusing the scattered x-rays).
* DetectorRadius: [m]  Radius of the detector (for focusing the scattered x-rays).
*
* %End
*******************************************************************************/

DEFINE COMPONENT SAXSShells



SETTING PARAMETERS (R = 100.0, Thickness = 5.0, Concentration = 0.01, DeltaRho = 1.0e-14, AbsorptionCrosssection = 0.0,
		    		xwidth, yheight, zdepth, SampleToDetectorDistance, DetectorRadius)



/*X-ray PARAMETERS (x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p)*/

DECLARE
%{
  double Prefactor;
  double Absorption;
  double q;
  double NumberDensity;

  double RBig;
  double RSmall;

  double VolumeBigSphere;
  double VolumeSmallSphere;
  double Volume;
%}

INITIALIZE
%{
  // Rescale concentration into number of aggregates per m^3 times 10^-4
  NumberDensity = Concentration * 6.02214129e19;

  // Computations
  if (!xwidth || !yheight || !zdepth) {
    printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP);
  }

  if (Thickness >= R) {
    printf ("%s: Thickness of shell larger than radius of shell!\n", NAME_CURRENT_COMP);
  }

  RBig = R + Thickness / 2.0;
  RSmall = R - Thickness / 2.0;

  VolumeBigSphere = 4.0 / 3.0 * PI * pow (RBig, 3);
  VolumeSmallSphere = 4.0 / 3.0 * PI * pow (RSmall, 3);

  Volume = VolumeBigSphere - VolumeSmallSphere;

  Prefactor = NumberDensity * pow (Volume, 2) * pow (DeltaRho, 2);

  Absorption = AbsorptionCrosssection;
%}

TRACE
%{
  // Declarations
  double l0;
  double l1;
  double l_full;
  double l;
  double l_1;
  double FormfactorBigSphere;
  double FormfactorSmallSphere;
  double Formfactor;
  double SolidAngle;
  double qx;
  double qy;
  double qz;
  double k;
  double dl;
  double kx_i;
  double ky_i;
  double kz_i;
  char Intersect = 0;

  // Computation
  Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);

  if (Intersect) {
    if (l0 < 0.0) {
      fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP);
      ABSORB;
    }

    // Compute properties of photon
    k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2));
    l_full = l1 - l0;
    dl = rand01 () * (l1 - l0) + l0;
    PROP_DL (dl);
    l = dl - l0;

    // Store properties of incoming photon
    kx_i = kx;
    ky_i = ky;
    kz_i = kz;

    // Generate new direction of photon
    randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius);

    NORM (kx, ky, kz);

    kx *= k;
    ky *= k;
    kz *= k;

    // Compute q
    qx = kx_i - kx;
    qy = ky_i - ky;
    qz = kz_i - kz;

    q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2));

    // Compute scattering
    FormfactorBigSphere = 3.0 * (sin (q * RBig) - q * RBig * cos (q * RBig)) / pow (q * RBig, 3);
    FormfactorSmallSphere = 3.0 * (sin (q * RSmall) - q * RSmall * cos (q * RSmall)) / pow (q * RSmall, 3);
    Formfactor = (FormfactorBigSphere * VolumeBigSphere - FormfactorSmallSphere * VolumeSmallSphere) / (VolumeBigSphere - VolumeSmallSphere);

    p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * pow (Formfactor, 2) * exp (-Absorption * (l + l1));

    SCATTER;
  }
%}

MCDISPLAY
%{
  box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0);
%}

END