/*******************************************************************************
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SANSNanodiscsFast
*
* %I
* Written by: Martin Cramer Pedersen (mcpe@nbi.dk)
* Date: October 17, 2012
* Origin: KU-Science
*
* A sample of monodisperse phospholipid bilayer nanodiscs in solution (water).
*
* %D
* A component very similar to EmptyNanodiscs.comp - however, the scattering
* profile is only computed once, and linear interpolation is the used to
* simulate the instrument.
*
* %P
* AxisRatio: []		                   Axis ratio of the bilayer patch.
* NumberOfLipids: []		             Number of lipids per nanodisc.
* AreaPerLipidHeadgroup: [AA^2]      Area per lipid headgroup - default is POPC.
* HeightOfMSP: [AA]                  Height of the belt protein - default is MSP1D1.
* VolumeOfOneMSP: [AA^3]             Volume of one belt protein - default is MSP1D1.
* 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.
* ScatteringLengthOfOneMSP: [cm]     Scattering length of one belt protein - default is MSP1D1.
* 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.
* Roughness: []		                   Factor used to smear the interfaces of the nanodisc.
* Concentration: [mM]                Concentration of sample.
* RhoSolvent: [cm/AA^3]              Scattering length density of solvent - default is D2O.
* 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).
* qMin: [AA^-1]                      Lowest q-value, for which a point is generated in the scattering profile
* qMax: [AA^-1]                      Highest q-value, for which a point is generated in the scattering profile
* NumberOfQBins: []		               Number of points generated in inital scattering profile.
*
* %E
*******************************************************************************/

DEFINE COMPONENT SANSNanodiscsFast



SETTING PARAMETERS (int NumberOfQBins = 200,
	AxisRatio = 1.4, NumberOfLipids = 130.0, AreaPerLipidHeadgroup = 65.0, HeightOfMSP = 24.0,
	VolumeOfOneMSP = 26296.5, VolumeOfHeadgroup = 319.0, VolumeOfCH2Tail = 818.8, VolumeOfCH3Tail = 108.6,
	ScatteringLengthOfOneMSP = 8.8E-10, ScatteringLengthOfHeadgroup = 7.05E-12, ScatteringLengthOfCH2Tail = -1.76E-12, ScatteringLengthOfCH3Tail = -9.15E-13,
	Roughness = 5.0, Concentration = 0.01, RhoSolvent = 6.4e-14, AbsorptionCrosssection = 0.0,
	xwidth, yheight, zdepth, SampleToDetectorDistance, DetectorRadius, qMin = 0.001, qMax = 1.0 )


DEPENDENCY " @GSLFLAGS@ "
NOACC

SHARE
%{
#include <gsl/gsl_sf_bessel.h>

// Functions used for compution the intensity from a given liposome
double FormfactorCylinder(double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta)
	{
		double const ProjectedRadius = sqrt(pow(MajorSemiAxis * sin(Beta), 2) + pow(MinorSemiAxis * cos(Beta), 2));

		double const Formfactor1 = j1(q * ProjectedRadius * sin(Alpha)) / (q * ProjectedRadius * sin(Alpha));
		double const Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0);

		return 2 * Formfactor1 * Formfactor2;
	}

	double IntensityOfEmptyNanodiscs(double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt,
									 double HeightOfBelt, double HeightOfLipids, double HeightOfTails, double HeightOfCH3,
									 double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail)
	{
		// Declaration
		double Intensity;
		double IntensityPart;

		double AmplitudeOfBelt;
		double AmplitudeOfHeads;
		double AmplitudeOfCH2Tail;
		double AmplitudeOfCH3Tail;

		double FormfactorOfBelt;
		double FormfactorOfHeads;
		double FormfactorOfCH2Tail;
		double FormfactorOfCH3Tail;

		const double OuterMajorSemiAxis = MajorSemiAxis + ThicknessOfBelt;
		const double OuterMinorSemiAxis = MinorSemiAxis + ThicknessOfBelt;

		const double VolumeOfBelt    = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis;
		const double VolumeOfHeads   = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis         - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis;
		const double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis          - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis;
		const double VolumeOfCH3Tail = PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis;

		// Variables needed for integration over alpha
		int i;
		const int NumberOfStepsInAlpha = 50;
		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 = 50;
		double Beta;
		const double BetaMin = 0.0;
		const double BetaMax = PI / 2.0;
		const double BetaStep = (BetaMax - BetaMin) / (1.0 * NumberOfStepsInBeta);

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

				// Compute formfactors
				FormfactorOfBelt = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis * FormfactorCylinder(q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) -
								    PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) /
								   (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis);

				FormfactorOfHeads = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) -
								     PI * HeightOfTails  * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails , Alpha, Beta)) /
								    (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis);

				FormfactorOfCH2Tail = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) -
								       PI * HeightOfCH3   * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3  , Alpha, Beta)) /
								      (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis);

				FormfactorOfCH3Tail = FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta);

				// Compute amplitudes
				AmplitudeOfBelt    = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt;
				AmplitudeOfHeads   = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads;
				AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail;
				AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail;

				// Perform integration
				IntensityPart = pow(AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2);

				Intensity += 2.0 / PI * sin(Alpha) * IntensityPart * AlphaStep * BetaStep;
			}
		}

		return Intensity;
	}

%}


DECLARE
%{
// Declarations
double Absorption;
double NumberDensity;

// Curve used in linear
double *qArray;
double *IArray;
%}

INITIALIZE
%{
// Functions used for compution the intensity from a given liposome

>>>>>>> b6fb3daa6... move function declarations to SHARE

	// Declarations
	double RhoBelt;
	double RhoHead;
	double RhoCH2Tail;
	double RhoCH3Tail;

	double DeltaRhoHead;
	double DeltaRhoBelt;
	double DeltaRhoCH2Tail;
	double DeltaRhoCH3Tail;

	double MajorSemiAxis;
	double MinorSemiAxis;
	double ThicknessOfBelt;

	double HeightOfBelt;
	double HeightOfLipids;
	double HeightOfTails;
	double HeightOfCH3;

	double AreaOfLipids;
	// Declarations
	double Absorption;
	double q;
	double NumberDensity;

	// Curve used in linear
	double *qArray;
	double *IArray;


INITIALIZE
%{
	const double ScatteringLengthOfWater = 2.82E-12;
	const double VolumeOfWater = 30.0;

	const double qStep = (qMax - qMin) / (1.0 * NumberOfQBins);
	int i;
	double qDummy;

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

	// Scattering properties of different components
	RhoBelt    = ScatteringLengthOfOneMSP    / VolumeOfOneMSP;
	RhoHead    = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup;
	RhoCH2Tail = ScatteringLengthOfCH2Tail   / VolumeOfCH2Tail;
	RhoCH3Tail = ScatteringLengthOfCH3Tail   / VolumeOfCH3Tail;

	DeltaRhoBelt	= RhoBelt    - RhoSolvent;
	DeltaRhoHead    = RhoHead    - RhoSolvent;
	DeltaRhoCH2Tail = RhoCH2Tail - RhoSolvent;
	DeltaRhoCH3Tail = RhoCH3Tail - RhoSolvent;

	// Geometric properties of different components
	AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0;

	MinorSemiAxis = sqrt(AreaOfLipids / (PI * AxisRatio));
	MajorSemiAxis = MinorSemiAxis * AxisRatio;

	HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup;
	HeightOfTails  = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup;
	HeightOfCH3    = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup;

	ThicknessOfBelt = sqrt(pow(MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0;

	// Compute scattering from nanodiscs in predecided points
	qArray = (double *) calloc(NumberOfQBins, sizeof(double));
	IArray = (double *) calloc(NumberOfQBins, sizeof(double));

	for (i = 0; i < NumberOfQBins; ++i) {
		qDummy = qMin + (i + 0.5) * qStep;

		qArray[i] = qDummy;
	 	IArray[i] = IntensityOfEmptyNanodiscs(qDummy, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids,
											  HeightOfTails, HeightOfCH3, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail);
	}
%}

TRACE
%{
	// Declarations
	double t0;
	double t1;
	double l_full;
	double l;
	double l1;
	double Intensity;
	double Weight;
	double IntensityPart;
	double SolidAngle;
	double q;
	double qx;
	double qy;
	double qz;
	double v;
	double dt;
	double vx_i;
	double vy_i;
	double vz_i;
	char Intersect = 0;
	int i;
	double Slope;
	double Offset;

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

		// Discard neutron, if q is out of range
		if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) {
		    ABSORB;
		}

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

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

		// Do a linear interpolation
		l1 = v * t1;

		Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]);
		Offset = IArray[i] - Slope * qArray[i];

		Intensity = (Slope * q + Offset) * exp(- pow(q * Roughness, 2));

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

		SCATTER;
	}
%}

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

END