/*******************************************************************************
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: SAXSLiposomes
*
* %Identification
* Written by: Martin Cramer Pedersen (mcpe@nbi.dk)
* Date: May 2, 2012
* Origin: KU-Science
*
* A sample of polydisperse liposomes in solution (water).
*
* %Description
* A component simulating the scattering from a box-shaped, thin solution (water)
* of liposomes described by a pentuple-shell model.
*
* Example: SAXSLiposomes( xwidth = 0.01, yheight = 0.01, zdepth = 0.01, SampleToDetectorDistance = 0.48, DetectorRadius = 0.1 )
*
* %Parameters
* Radius: [AA]      Average thickness of the liposomes.
* Thickness: [AA]   Thickness of the bilayer.
* SigmaRadius: []   Relative Gaussian deviation of the radius in the distribution of liposomes.
* nRadius: []       Number of bins in Radius for polydisperse distribution.
* VolumeOfHeadgroup: [AA^3]   Volume of one lipid headgroup - default is POPC.
* VolumeOfCH2Tail: [AA^3]     Volume of the CH2-chains of one lipid - default is POPC.
* VolumeOfCH3Tail: [AA^3]     Volume of the CH3-tails of one lipid - default is POPC.
* ScatteringLengthOfHeadgroup: [cm] Scattering length of one lipid headgroup - default is POPC.
* ScatteringLengthOfCH2Tail: [cm]   Scattering length of the CH2-chains of one lipid - default is POPC.
* ScatteringLengthOfCH3Tail: [cm]   Scattering length of the CH3-tails of one lipid - default is POPC.
* 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).
*
* %End
*******************************************************************************/

DEFINE COMPONENT SAXSLiposomes



SETTING PARAMETERS (Radius = 800.0, Thickness = 38.89, SigmaRadius = 0.20, nRadius=100,
		    VolumeOfHeadgroup = 319.0, VolumeOfCH2Tail = 818.8, VolumeOfCH3Tail = 108.6,
		    ScatteringLengthOfHeadgroup = 4.62E-11, ScatteringLengthOfCH2Tail = 6.71E-11, ScatteringLengthOfCH3Tail = 5.08E-12,
		    Concentration = 0.01, AbsorptionCrosssection = 0.0,
		    xwidth, yheight, zdepth, SampleToDetectorDistance, DetectorRadius)



/*X-ray Parameters (x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p)*/
SHARE
%{
  // Functions used for compution the intensity from a given liposome
#pragma acc routine
  double FormfactorSphere(double q, double R)
  {
    return 3 * (sin(q * R) - q * R * cos(q * R)) / pow(q * R, 3);
  }

#pragma acc routine
  double IntensityOfLiposome(double q, double R, double ThicknessHead, double ThicknessTail, double ThicknessCH3, double DeltaRhoHead, double DeltaRhoCH2, double DeltaRhoCH3)
  {
    double RHeadOut,RTailOut,RCH3Out;
    double RCH3In,RTailIn,RHeadIn;
    double VolumeHeadOut, VolumeTailOut, VolumeCH3,VolumeTailIn,VolumeHeadIn;
    double AmplitudeHeadOut,AmplitudeTailOut,AmplitudeCH3,AmplitudeTailIn,AmplitudeHeadIn;
    double Intensity;

    RHeadOut = R + ThicknessHead + ThicknessTail + ThicknessCH3;
    RTailOut = R + ThicknessTail + ThicknessCH3;
    RCH3Out  = R + ThicknessCH3;
    RCH3In   = R - ThicknessCH3;
    RTailIn  = R - ThicknessTail - ThicknessCH3;
    RHeadIn  = R - ThicknessHead - ThicknessTail - ThicknessCH3;

    VolumeHeadOut = 4.0 / 3.0 * PI * (pow(RHeadOut, 3) - pow(RTailOut, 3));
    VolumeTailOut = 4.0 / 3.0 * PI * (pow(RTailOut, 3) - pow(RCH3Out, 3));
    VolumeCH3     = 4.0 / 3.0 * PI * (pow(RCH3Out, 3)  - pow(RCH3In, 3));
    VolumeTailIn  = 4.0 / 3.0 * PI * (pow(RCH3In, 3)   - pow(RTailIn, 3));
    VolumeHeadIn  = 4.0 / 3.0 * PI * (pow(RTailIn, 3)  - pow(RHeadIn, 3));

    AmplitudeHeadOut = DeltaRhoHead * VolumeHeadOut *
      (pow(RHeadOut, 3) * FormfactorSphere(q, RHeadOut) - pow(RTailOut, 3) * FormfactorSphere(q, RTailOut)) / (pow(RHeadOut, 3) - pow(RTailOut, 3));

    AmplitudeTailOut = DeltaRhoCH2  * VolumeTailOut *
      (pow(RTailOut, 3) * FormfactorSphere(q, RTailOut) - pow(RCH3Out, 3) * FormfactorSphere(q, RCH3Out)) / (pow(RTailOut, 3) - pow(RCH3Out, 3));

    AmplitudeCH3 = DeltaRhoCH3  * VolumeCH3 *
      (pow(RCH3Out, 3) * FormfactorSphere(q, RCH3Out) - pow(RCH3In, 3) * FormfactorSphere(q, RCH3In)) / (pow(RCH3Out, 3) - pow(RCH3In, 3));

    AmplitudeTailIn  = DeltaRhoCH2  * VolumeTailIn *
      (pow(RCH3In, 3) * FormfactorSphere(q, RCH3In) - pow(RTailIn, 3) * FormfactorSphere(q, RTailIn)) / (pow(RCH3In, 3) - pow(RTailIn, 3));

    AmplitudeHeadIn = DeltaRhoHead * VolumeHeadIn *
      (pow(RTailIn, 3) * FormfactorSphere(q, RTailIn) - pow(RHeadIn, 3) * FormfactorSphere(q, RHeadIn)) / (pow(RTailIn, 3) - pow(RHeadIn, 3));

    Intensity = pow(AmplitudeHeadOut + AmplitudeTailOut + AmplitudeCH3 + AmplitudeTailIn + AmplitudeHeadIn, 2);

    return Intensity;
  }
%}

DECLARE
%{
	double Absorption;
	double NumberDensity;

	double RMin;
	double RMax;
	double RStep;

	// Scattering lengths
	double DeltaRhoHead;
	double DeltaRhoCH2Tail;
	double DeltaRhoCH3Tail;

	// Thickness
	double ThicknessOfHead;
	double ThicknessOfCH2Tail;
	double ThicknessOfCH3Tail;
%}

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

	Absorption = AbsorptionCrosssection;


	RMin = Radius - 3.0 * SigmaRadius * Radius;

	if (RMin < Thickness / 2.0) {
		RMin = Thickness / 2.0;
	}

	RMax = Radius + 3.0 * SigmaRadius * Radius;

	RStep = (RMax - RMin) / (1.0f * nRadius);

	// Molecular properties of liposomes
	const double ScatteringLengthOfWater = 2.82E-12;
	const double VolumeOfWater = 30.0;

	double RhoWater   = ScatteringLengthOfWater     / VolumeOfWater;
	double RhoHead    = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup;
	double RhoCH2Tail = ScatteringLengthOfCH2Tail   / VolumeOfCH2Tail;
	double RhoCH3Tail = ScatteringLengthOfCH3Tail   / VolumeOfCH3Tail;

	DeltaRhoHead    = RhoHead    - RhoWater;
	DeltaRhoCH2Tail = RhoCH2Tail - RhoWater;
	DeltaRhoCH3Tail = RhoCH3Tail - RhoWater;

	ThicknessOfHead    = Thickness * VolumeOfHeadgroup / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail);
	ThicknessOfCH2Tail = Thickness * VolumeOfCH2Tail   / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail);
	ThicknessOfCH3Tail = Thickness * VolumeOfCH3Tail   / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail);
%}

TRACE
%{
	// Declarations	
	double l0; 
	double l1; 
	double l_full;
	double l;
	double l_1;
	double Intensity;
	double Weight1;
	double Weight2;
	double IntensityPart;
	double SolidAngle;
	double qx; 
	double qy; 
	double qz;
	double k;
	double q;
	double dl;
	double kx_i;
	double ky_i;
	double kz_i;
	int Intersect = 0;
	double R;
	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(qx*qx+qy*qy+qz*qz);
		// Compute scattering
		Intensity = 0.0;
		Weight1 = 1.0 / (SigmaRadius * Radius * sqrt(2.0 * PI));
		for (i = 0; i < nRadius; i++) {
		  R = RMin + RStep* (i + 0.5);
		  IntensityPart = IntensityOfLiposome(q, R, ThicknessOfHead, ThicknessOfCH2Tail, ThicknessOfCH3Tail, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail);
		  Weight2 = exp(- pow((R - Radius) / (sqrt(2.0) * SigmaRadius * Radius), 2));

		  Intensity += Weight1 * Weight2 * IntensityPart * RStep;
		}

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

		SCATTER;
	}
%}

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

END