/*******************************************************************************
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SANSNanodiscsWithTags
*
* %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) - with
* histidine tag still on the belt proteins.
*
* %D
* A component simulating the scattering from a box-shaped, thin solution (water)
* of monodisperse phospholipid bilayer nanodiscs - with histidine tags still
* on the belt proteins.
*
* %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.
* RadiusOfGyrationForHisTag: [AA]    Radius of gyration for the his-tag.
* 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.
* VolumeOfOneHisTag: [AA^3]          Volume of one his-tag.
* 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.
* ScatteringLengthOfOneHisTag: [cm]  Scattering length of one his-tag.
* 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).
*
* %E
*******************************************************************************/

DEFINE COMPONENT SANSNanodiscsWithTags



SETTING PARAMETERS (AxisRatio = 1.4, NumberOfLipids = 130.0, AreaPerLipidHeadgroup = 65.0, HeightOfMSP = 24.0, RadiusOfGyrationForHisTag = 10.0,
VolumeOfOneMSP = 26296.5, VolumeOfHeadgroup = 319.0, VolumeOfCH2Tail = 818.8, VolumeOfCH3Tail = 108.6, VolumeOfOneHisTag = 2987.3,
ScatteringLengthOfOneMSP = 8.8E-10, ScatteringLengthOfHeadgroup = 7.05E-12, ScatteringLengthOfCH2Tail = -1.76E-12, ScatteringLengthOfCH3Tail = -9.15E-13,ScatteringLengthOfOneHisTag = 1.10E-10,
Roughness = 5.0, Concentration = 0.01, RhoSolvent = 6.4e-14, AbsorptionCrosssection = 0.0,
xwidth, yheight, zdepth, SampleToDetectorDistance, DetectorRadius)


DEPENDENCY " @GSLFLAGS@ "
NOACC

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

// Functions used for compution the formfactor of a cylinder
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 = gsl_sf_bessel_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;
	}

	// Functions used to compute the characteristics of the histidine tags
	double Debye(double q, double r)
	{
		// Variables needed in functions
		double ReturnValue;

		// Constants needed in function
		double Dummy1;
		double Dummy2;

		// Computation
		Dummy1 = (q * q) * (r * r);
		Dummy2 = 2 * ((exp(-Dummy1) + Dummy1 - 1) / (Dummy1 * Dummy1));

		ReturnValue = Dummy2;

		return ReturnValue;
	}

	double XiRodWithoutEndcaps(double q, double r, double l, double Alpha)
	{
		// Variables needed in function
		double ReturnValue;

		// Computation
		ReturnValue = j0(q * r * sin(Alpha)) * sin(q * l * cos(Alpha) / 2.f) / (q * l * cos(Alpha) / 2.f);

		return ReturnValue;
	}

	double PsiHammouda(double x)
	{
		// Variables used in function
		double ReturnValue;

		// Computation
		ReturnValue = (1.0 - exp(-pow(x, 2))) / pow(x, 2);

		return ReturnValue;
	}

	// Computation of the intensity from a given nanodisc
	double IntensityOfEmptyNanodiscsWithTags(
		double q,
		double MajorSemiAxis,
		double MinorSemiAxis,
		double ThicknessOfBelt,
		double HeightOfBelt,
		double HeightOfLipids,
		double HeightOfTails,
		double HeightOfCH3,
		double RadiusOfGyrationForHisTag,
		double VolumeOfOneHisTag,
		double DeltaRhoBelt,
		double DeltaRhoHead,
		double DeltaRhoCH2Tail,
		double DeltaRhoCH3Tail,
		double DeltaRhoTag)
	{
		// Declarations
		double Intensity;
		double IntensityPart;
		double IntensityDisc;

		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 AverageRadiusOfBelt = sqrt(pow(OuterMajorSemiAxis, 2) + pow(OuterMinorSemiAxis, 2));

		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;

		double Autocorrelation;
		double Crosscorrelation;
		double Disccorrelation;
		double FormfactorOfTags;

		// 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
				IntensityDisc = pow(AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2);

				// Add the histidine tags
		        FormfactorOfTags = XiRodWithoutEndcaps(q, AverageRadiusOfBelt + RadiusOfGyrationForHisTag, HeightOfBelt, Alpha);

		        Crosscorrelation = 2.0f * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * pow(FormfactorOfTags, 2) * pow(PsiHammouda(q * RadiusOfGyrationForHisTag), 2);

		        Autocorrelation  = 2.0f * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * Debye(q, RadiusOfGyrationForHisTag);

		        Disccorrelation  = 4.0f * (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail) *
									      DeltaRhoTag * VolumeOfOneHisTag * FormfactorOfTags * PsiHammouda(q * RadiusOfGyrationForHisTag);

		        IntensityPart = IntensityDisc + Crosscorrelation + Autocorrelation + Disccorrelation;

				// Return the final intensity
				Intensity += 2 / PI * sin(Alpha) * IntensityPart * AlphaStep * BetaStep;
			}
		}
		return Intensity;
	}
%}

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

	// Scattering lengths
	double RhoBelt;
	double RhoHead;
	double RhoCH2Tail;
	double RhoCH3Tail;
	double RhoTag;

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

	// Geometric properties
	double MajorSemiAxis;
	double MinorSemiAxis;
	double ThicknessOfBelt;

	double HeightOfBelt;
	double HeightOfLipids;
	double HeightOfTails;
	double HeightOfCH3;
%}

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;

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

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

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

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

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

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

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

	// 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 = exp(- pow(q * Roughness, 2)) * IntensityOfEmptyNanodiscsWithTags(q, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, HeightOfCH3,
																			 		 RadiusOfGyrationForHisTag, VolumeOfOneHisTag, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail, DeltaRhoTag);

		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