/*******************************************************************************
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SANSCurve
*
* %I
* Written by: Søren Kynde (kynde@nbi.dk)
* Rewritten by: Martin Cramer Pedersen (mcpe@nbi.dk)
* Date: October 17, 2012
* Version: $Revision$
* Origin: KU-Science
* Release: McStas 1.12c
*
* A component mimicking the scattering from a given I(q)-curve by using
* linear interpolation between the given points.
*
* %D
* 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 neutron.
*
* %P
* 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 neutrons).
* DetectorRadius          	: [m] 		Radius of the detector (for focusing the 
*										scattered neutrons).
* FileWithCurve				: []		Datafile with the given I(q).
*
* %E
*******************************************************************************/

DEFINE COMPONENT SANSCurve

DEFINITION PARAMETERS ()

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

OUTPUT PARAMETERS ()

DECLARE
%{
	double Absorption;
	double q;
	double* Q;
	double* I;
	double ForwardScattering;
	double Prefactor;
	int NumberOfDatapoints;
	double NumberDensity; 

	// 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** I) 
	{
		// Declarations
        FILE* File;
        
        int i = 0;
		int NumberOfDatapoints;

        char Line[256];
        
		double *IntensityArray;
		double *qArray;

		// Reading file
        if ((File = fopen(Filename, "r")) == 0) {
            printf("Cannot open file: %s...\n", Filename);
            exit(0);
        }

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

        *I = IntensityArray;
        *Q = qArray;

        return NumberOfDatapoints;
    }
%}

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


	// Initializing curve from file
    NumberOfDatapoints = LoadCurve(FileWithCurve, &Q, &I);
    ForwardScattering = I[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 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;
	double Slope;	
	double Offset;
	int i;

	// Computation
	Intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
    
	//fprintf(stderr, "Start of 'SANSCurve'... \n");

	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 < Q[0]) || (q > Q[NumberOfDatapoints - 1])) {
		    ABSORB;
		}
   
		// Find the first value of q in the curve larger than that of the neutron
		i = 1;

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

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

		Slope = (I[i] - I[i - 1]) / (Q[i] - Q[i - 1]);
		Offset = I[i] - Slope * Q[i];

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

		p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * Intensity * exp(- Absorption * (l + l1) / v);
		
		SCATTER;
	}
	
	//fprintf(stderr, "End of 'SANSCurve'... \n");
	
%}

MCDISPLAY
%{
	magnify("xyz");
	box(0, 0, 0, xwidth, yheight, zdepth);
%}

END