/*******************************************************************************
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: SAXSCurve
*
* %Identification
* Written by: Martin Cramer Pedersen (mcpe@nbi.dk)
* Based on a SANS-component in McStas by Søren Kynde
* Date: May 2, 2012
* Origin: KU-Science
*
* A component mimicking the scattering from a given I(q)-curve by using
* linear interpolation between the given points.
*
* %Description
* A box-shaped component simulating the scattering from any given I(q)-curve.
* The component uses linear interpolation to generate the points necessary to
* compute the scattering of any given photon.
*
* %Parameters
* DeltaRho:       [cm/AA^3]	Excess scattering length density of the particles.
* Volume:         [AA^3]	Volume of the particles.
* Concentration:  [mM]		Concentration of sample.
* 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).
* FileWithCurve:  [str]		Datafile with the given I(q).
*
* %End
*******************************************************************************/

DEFINE COMPONENT SAXSCurve



SETTING PARAMETERS (DeltaRho = 1.0e-14, Volume = 10000.0, 
    Concentration = 0.01, AbsorptionCrosssection = 0.0,
    xwidth, yheight, zdepth, 
    SampleToDetectorDistance, DetectorRadius,
    string FileWithCurve = "Curve.dat")


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


SHARE
%{
  // Function used to determine the number of datapoints in the input file
  int CountLines(FILE* File)
  {
    // Declarations        
    double Dummy1;
    double Dummy2;
    char Line[256];
    int NumberOfDatapoints = 0;

    // I/O
    while (fgets(Line, sizeof(Line), File) != NULL) {

      if (sscanf(Line, "%lf %lf", &Dummy1, &Dummy2) == 2) {
        ++NumberOfDatapoints;
      }
    }

    return NumberOfDatapoints;
  }

  // Function used to extract the scattering profile from a given curve
  int LoadCurve(char Filename[], double** Q, double** II)
  {
    // Declarations
    FILE* File;

    int i = 0;
    int NumberOfDatapoints;

    char Line[256];

    double *IntensityArray;
    double *qArray;

    // Reading file
    if ((File = fopen(Filename, "r")) == 0) {
      fprintf(stderr,"ERROR: (%s): Cannot open file: %s...\n", "SAXSCurve.comp",Filename);
      exit(-1);
    }

    NumberOfDatapoints = CountLines(File);

    qArray         = (double *) calloc(NumberOfDatapoints, sizeof(double));
    IntensityArray = (double *) calloc(NumberOfDatapoints, sizeof(double));

    rewind(File);

    while (i < NumberOfDatapoints && fgets(Line, sizeof(Line), File) != NULL) {
      if (sscanf(Line, "%lf %lf", &qArray[i], &IntensityArray[i]) == 2) {
        ++i;
      }
    }

    *II = IntensityArray;
    *Q = qArray;

    return NumberOfDatapoints;
  }
%}


DECLARE
%{
  double Absorption;
  double* QQ;
  double* II;
  double ForwardScattering;
  double Prefactor;
  int NumberOfDatapoints;
  double NumberDensity; 
%}

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

  // Initializing curve from file
  NumberOfDatapoints = LoadCurve(FileWithCurve, &QQ, &II);
  ForwardScattering = II[0];

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

  Absorption = AbsorptionCrosssection;
  Prefactor = NumberDensity * pow(DeltaRho * Volume, 2);
%}

TRACE
%{
  // Declarations	
  double l0; 
  double l1; 
  double l_full;
  double l;
  double l_1;
  double Intensity;
  double Weight;
  double IntensityPart;
  double SolidAngle;
  double q;
  double qx; 
  double qy; 
  double qz;
  double k;
  double dl;
  double kx_i;
  double ky_i;
  double kz_i;
  char Intersect = 0;
  double Slope;	
  double Offset;
  int i;

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

    // Discard photon, if q is out of range
    if ((q < QQ[0]) || (q > QQ[NumberOfDatapoints - 1])) {
      ABSORB;
    }

    // Find the first value of q in the curve larger than that of the photon
    i = 1;

    while (q > QQ[i]) {
      ++i;
    }

    // Do a linear interpolation
    Slope = (II[i] - II[i - 1]) / (QQ[i] - QQ[i - 1]);
    Offset = II[i] - Slope * QQ[i];

    Intensity = (Slope * q + Offset) / ForwardScattering;

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

    SCATTER;
  }
%}

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

END