/*******************************************************************************
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SANSEllipticCylinders
*
* %I
* Written by: Martin Cramer Pedersen (mcpe@nbi.dk)
* Date: October 17, 2012
* Origin: KU-Science
*
* A sample of monodisperse cylindrical particles with elliptic cross section in 
* solution.
*
* %D
* A component simulating the scattering from a box-shaped, thin solution
* of monodisperse, cylindrical particles with elliptic cross section.
*
* %P
* R1: [AA]                       First semiaxis of the cross section of the
*						elliptic cylinder.
* R2: [AA]                       Second semiaxis of the cross section of the
*						elliptic cylinder.
* Height: [AA]                   Height of the elliptic cylinder.
* 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 neutrons).
* DetectorRadius: [m]            Radius of the detector (for focusing the
* 						scattered neutrons).
* 
* %E
*******************************************************************************/

DEFINE COMPONENT SANSEllipticCylinders



SETTING PARAMETERS (R1 = 20.0, R2 = 40.0, Height = 100.0, Concentration = 0.01, DeltaRho = 1.0e-14, AbsorptionCrosssection = 0.0,
xwidth, yheight, zdepth, SampleToDetectorDistance, DetectorRadius)


DEPENDENCY " @GSLFLAGS@ "
NOACC

/*X-ray Parameters (x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p)*/
SHARE
%{
  #include <gsl/gsl_sf_bessel.h>
%}
DECLARE
%{
  // Declarations
  double Prefactor;
  double Absorption;
  double q;
  double NumberDensity;
%}

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);
  }

  Prefactor = NumberDensity * pow (PI * Height * R1 * R2, 2) * pow (DeltaRho, 2);

  Absorption = AbsorptionCrosssection;
%}

TRACE
%{
  // Declarations
  double t0;
  double t1;
  double l_full;
  double l;
  double l1;
  double Formfactor1;
  double Formfactor2;
  double Intensity;
  double SolidAngle;
  double qx;
  double qy;
  double qz;
  double v;
  double dt;
  double vx_i;
  double vy_i;
  double vz_i;
  double ProjectedRadius;
  char Intersect = 0;

  // Variables needed for integration over alpha
  int i;
  const int NumberOfStepsInAlpha = 30;
  double Alpha;
  const double AlphaMin = 0.0;
  const double AlphaMax = PI / 2.0;
  const double AlphaStep = (AlphaMax - AlphaMin) / (1.0 * NumberOfStepsInAlpha);

  // Variables needed in integration over beta
  int j;
  const int NumberOfStepsInBeta = 30;
  double Beta;
  const double BetaMin = 0.0;
  const double BetaMax = PI / 2.0;
  const double BetaStep = (BetaMax - BetaMin) / (1.0 * NumberOfStepsInBeta);

  // Computation
  Intersect = box_intersect (&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);

  if (Intersect) {

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

    // Compute properties of neutron
    v = sqrt (pow (vx, 2) + pow (vy, 2) + pow (vz, 2));
    l_full = v * (t1 - t0);
    dt = rand01 () * (t1 - t0) + t0;
    PROP_DT (dt);
    l = v * (dt - t0);

    // Store properties of incoming neutron
    vx_i = vx;
    vy_i = vy;
    vz_i = vz;

    // Generate new direction of neutron
    randvec_target_circle (&vx, &vy, &vz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius);

    NORM (vx, vy, vz);

    vx *= v;
    vy *= v;
    vz *= v;

    // Compute q
    qx = V2K * (vx_i - vx);
    qy = V2K * (vy_i - vy);
    qz = V2K * (vz_i - vz);

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

    // Compute scattering
    l1 = v * t1;

    Intensity = 0.0;

    for (i = 0; i < NumberOfStepsInAlpha; ++i) {
      Alpha = (i + 0.5) * AlphaStep;

      for (j = 0; j < NumberOfStepsInBeta; ++j) {
        Beta = (j + 0.5) * BetaStep;
        ProjectedRadius = sqrt (pow (R1 * sin (Beta), 2) + pow (R2 * cos (Beta), 2));

        Formfactor1 = gsl_sf_bessel_J1 (q * ProjectedRadius * sin (Alpha)) / (q * ProjectedRadius * sin (Alpha));
        Formfactor2 = sin (q * Height * cos (Alpha) / 2.0) / (q * Height * cos (Alpha) / 2.0);

        Intensity += 2 / PI * sin (Alpha) * Prefactor * pow (2 * Formfactor1 * Formfactor2, 2) * AlphaStep * BetaStep;
      }
    }

    p *= l_full * SolidAngle / (4.0 * PI) * Intensity * exp (-Absorption * (l + l1) / v);

    SCATTER;
  }
%}

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

END