/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SasView_fractal_core_shell
*
* %Identification
* Written by: Jose Robledo
* Based on sasmodels from SasView
* Origin: FZJ / DTU / ESS DMSC
*
*
* SasView fractal_core_shell model component as sample description.
*
* %Description
*
* SasView_fractal_core_shell component, generated from fractal_core_shell.c in sasmodels.
*
* Example: 
*  SasView_fractal_core_shell(radius, thickness, sld_core, sld_shell, sld_solvent, volfraction, fractal_dim, cor_length, 
*     model_scale=1.0, model_abs=0.0, xwidth=0.01, yheight=0.01, zdepth=0.005, R=0, 
*     int target_index=1, target_x=0, target_y=0, target_z=1,
*     focus_xw=0.5, focus_yh=0.5, focus_aw=0, focus_ah=0, focus_r=0, 
*     pd_radius=0.0, pd_thickness=0.0, pd_cor_length=0.0)
*
* %Parameters
* INPUT PARAMETERS:
* radius: [Ang] ([0.0, inf]) Sphere core radius.
* thickness: [Ang] ([0.0, inf]) Sphere shell thickness.
* sld_core: [1e-6/Ang^2] ([-inf, inf]) Sphere core scattering length density.
* sld_shell: [1e-6/Ang^2] ([-inf, inf]) Sphere shell scattering length density.
* sld_solvent: [1e-6/Ang^2] ([-inf, inf]) Solvent scattering length density.
* volfraction: [] ([0.0, inf]) Volume fraction of building block spheres.
* fractal_dim: [] ([0.0, 6.0]) Fractal dimension.
* cor_length: [Ang] ([0.0, inf]) Correlation length of fractal-like aggregates.
* Optional parameters:
* model_abs: [ ] Absorption cross section density at 2200 m/s.
* model_scale: [ ] Global scale factor for scattering kernel. For systems without inter-particle interference, the form factors can be related to the scattering intensity by the particle volume fraction.
* xwidth: [m] ([-inf, inf]) Horiz. dimension of sample, as a width.
* yheight: [m] ([-inf, inf]) vert . dimension of sample, as a height for cylinder/box
* zdepth: [m] ([-inf, inf]) depth of sample
* R: [m] Outer radius of sample in (x,z) plane for cylinder/sphere.
* target_x: [m] relative focus target position.
* target_y: [m] relative focus target position.
* target_z: [m] relative focus target position.
* target_index: [ ] Relative index of component to focus at, e.g. next is +1.
* focus_xw: [m] horiz. dimension of a rectangular area.
* focus_yh: [m], vert. dimension of a rectangular area.
* focus_aw: [deg], horiz. angular dimension of a rectangular area.
* focus_ah: [deg], vert. angular dimension of a rectangular area.
* focus_r: [m] case of circular focusing, focusing radius.
* pd_radius: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable.
* pd_thickness: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable.
* pd_cor_length: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable
*
* %Link
* %End
*******************************************************************************/
DEFINE COMPONENT SasView_fractal_core_shell

SETTING PARAMETERS (
        radius=60.0,
        thickness=10.0,
        sld_core=1.0,
        sld_shell=2.0,
        sld_solvent=3.0,
        volfraction=0.05,
        fractal_dim=2.0,
        cor_length=100.0,
        model_scale=1.0,
        model_abs=0.0,
        xwidth=0.01,
        yheight=0.01,
        zdepth=0.005,
        R=0,
        target_x=0,
        target_y=0,
        target_z=1,
        int target_index=1,
        focus_xw=0.5,
        focus_yh=0.5,
        focus_aw=0,
        focus_ah=0,
        focus_r=0,
        pd_radius=0.0,
        pd_thickness=0.0,
        pd_cor_length=0.0)


SHARE %{
#ifdef OPENACC
#define USE_CUDA
#endif
%include "sas_kernel_header.c"

/* BEGIN Required header for SASmodel fractal_core_shell */
#define HAS_Iq
#define FORM_VOL

#ifndef SAS_HAVE_sas_3j1x_x
#define SAS_HAVE_sas_3j1x_x

#line 1 "sas_3j1x_x"
/**
* Spherical Bessel function 3*j1(x)/x
*
* Used for low q to avoid cancellation error.
* Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but
* in this case it is likely cancellation errors in the original expression
* using double precision that are the source.
*/
double sas_3j1x_x(double q);

// The choice of the number of terms in the series and the cutoff value for
// switching between series and direct calculation depends on the numeric
// precision.
//
// Point where direct calculation reaches machine precision:
//
//   single machine precision eps 3e-8 at qr=1.1 **
//   double machine precision eps 4e-16 at qr=1.1
//
// Point where Taylor series reaches machine precision (eps), where taylor
// series matches direct calculation (cross) and the error at that point:
//
//   prec   n eps  cross  error
//   single 3 0.28  0.4   6.2e-7
//   single 4 0.68  0.7   2.3e-7
//   single 5 1.18  1.2   7.5e-8
//   double 3 0.01  0.03  2.3e-13
//   double 4 0.06  0.1   3.1e-14
//   double 5 0.16  0.2   5.0e-15
//
// ** Note: relative error on single precision starts increase on the direct
// method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be
// safe for the sans range, with objects of 100 nm supported to a q of 0.1
// while maintaining 5 digits of precision.  For usans/sesans, the objects
// are larger but the q is smaller, so again it should be fine.
//
// See explore/sph_j1c.py for code to explore these ranges.

// Use 4th order series
#if FLOAT_SIZE>4
#define SPH_J1C_CUTOFF 0.1
#else
#define SPH_J1C_CUTOFF 0.7
#endif
#pragma acc routine seq
double sas_3j1x_x(double q)
{
    // 2017-05-18 PAK - support negative q
    if (fabs(q) < SPH_J1C_CUTOFF) {
        const double q2 = q*q;
        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.)))));
    } else {
        double sin_q, cos_q;
        SINCOS(q, sin_q, cos_q);
        return 3.0*(sin_q/q - cos_q)/(q*q);
    }
}


#endif // SAS_HAVE_sas_3j1x_x


#ifndef SAS_HAVE_sas_gamma
#define SAS_HAVE_sas_gamma

#line 1 "sas_gamma"
/*
The wrapper for gamma function from OpenCL and standard libraries
The OpenCL gamma function fails miserably on values lower than 1.0
while works fine on larger values.
We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute
to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1)
For t < 0, we use Gamma(t) = pi / ( Gamma(1 - t) * sin(pi * t) )
*/

#if defined(NEED_TGAMMA)
#pragma acc routine seq
static double cephes_stirf(double x)
{
	const double MAXSTIR=143.01608;
	const double SQTPI=2.50662827463100050242E0;
	double y, w, v;

	w = 1.0 / x;

	w = ((((
		7.87311395793093628397E-4*w
		- 2.29549961613378126380E-4)*w
		- 2.68132617805781232825E-3)*w
		+ 3.47222221605458667310E-3)*w
		+ 8.33333333333482257126E-2)*w
		+ 1.0;
	y = exp(x);
	if (x > MAXSTIR)
	{ /* Avoid overflow in pow() */
		v = pow(x, 0.5 * x - 0.25);
		y = v * (v / y);
	}
	else
	{
		y = pow(x, x - 0.5) / y;
	}
	y = SQTPI * y * w;
	return(y);
}

#pragma acc routine seq
static double tgamma(double x) {
	double p, q, z;
	int sgngam;
	int i;

	sgngam = 1;
	if (isnan(x))
		return(x);
	q = fabs(x);

	if (q > 33.0)
	{
		if (x < 0.0)
		{
			p = floor(q);
			if (p == q)
			{
				return (NAN);
			}
			i = p;
			if ((i & 1) == 0)
				sgngam = -1;
			z = q - p;
			if (z > 0.5)
			{
				p += 1.0;
				z = q - p;
			}
			z = q * sin(M_PI * z);
			if (z == 0.0)
			{
				return(NAN);
			}
			z = fabs(z);
			z = M_PI / (z * cephes_stirf(q));
		}
		else
		{
			z = cephes_stirf(x);
		}
		return(sgngam * z);
	}

	z = 1.0;
	while (x >= 3.0)
	{
		x -= 1.0;
		z *= x;
	}

	while (x < 0.0)
	{
		if (x > -1.E-9)
			goto small;
		z /= x;
		x += 1.0;
	}

	while (x < 2.0)
	{
		if (x < 1.e-9)
			goto small;
		z /= x;
		x += 1.0;
	}

	if (x == 2.0)
		return(z);

	x -= 2.0;
	p = (((((
		+1.60119522476751861407E-4*x
		+ 1.19135147006586384913E-3)*x
		+ 1.04213797561761569935E-2)*x
		+ 4.76367800457137231464E-2)*x
		+ 2.07448227648435975150E-1)*x
		+ 4.94214826801497100753E-1)*x
		+ 9.99999999999999996796E-1;
	q = ((((((
		-2.31581873324120129819E-5*x
		+ 5.39605580493303397842E-4)*x
		- 4.45641913851797240494E-3)*x
		+ 1.18139785222060435552E-2)*x
		+ 3.58236398605498653373E-2)*x
		- 2.34591795718243348568E-1)*x
		+ 7.14304917030273074085E-2)*x
		+ 1.00000000000000000320E0;
	return(z * p / q);

small:
	if (x == 0.0)
	{
		return (NAN);
	}
	else
		return(z / ((1.0 + 0.5772156649015329 * x) * x));
}
#endif // NEED_TGAMMA

#pragma acc routine seq
inline double sas_gamma(double x)
{
    // Note: the builtin tgamma can give slow and unreliable results for x<1.
    // The following transform extends it to zero and to negative values.
    // It should return NaN for zero and negative integers but doesn't.
    // The accuracy is okay but not wonderful for negative numbers, maybe
    // one or two digits lost in the calculation. If higher accuracy is
    // needed, you could test the following loop:
    //    double norm = 1.;
    //    while (x<1.) { norm*=x; x+=1.; }
    //    return tgamma(x)/norm;
    return (x<0. ? M_PI/tgamma(1.-x)/sin(M_PI*x) : tgamma(x+1)/x);
}


#endif // SAS_HAVE_sas_gamma


#ifndef SAS_HAVE_core_shell
#define SAS_HAVE_core_shell

#line 1 "core_shell"
/*******************************************************************

core_shell_kernel

Form factor used in core_shell and fractal_core_shell

********************************************************************/
#pragma acc routine seq
static
double core_shell_fq(double q,
                         double radius,
                         double thickness,
                         double core_sld,
                         double shell_sld,
                         double solvent_sld)
{
    // Core first, then add in shell
    const double core_qr = q * radius;
    const double core_contrast = core_sld - shell_sld;
    const double core_bes = sas_3j1x_x(core_qr);
    const double core_volume = M_4PI_3 * cube(radius);
    double f = core_volume * core_bes * core_contrast;

    // Now the shell
    const double shell_qr = q * (radius + thickness);
    const double shell_contrast = shell_sld - solvent_sld;
    const double shell_bes = sas_3j1x_x(shell_qr);
    const double shell_volume = M_4PI_3 * cube(radius + thickness);
    f += shell_volume * shell_bes * shell_contrast;
    return f;
}

// Deprecated function: use core_shell_fq instead.
#pragma acc routine seq
static
double core_shell_kernel(double q,
                         double radius,
                         double thickness,
                         double core_sld,
                         double shell_sld,
                         double solvent_sld)
{
    const double fq = core_shell_fq(q, radius, thickness, core_sld, shell_sld, solvent_sld);
    return 1.0e-4 * fq*fq;
}


#endif // SAS_HAVE_core_shell


#ifndef SAS_HAVE_fractal_sq
#define SAS_HAVE_fractal_sq

#line 1 "fractal_sq"
#pragma acc routine seq
static double
fractal_sq(double q, double radius, double fractal_dim, double cor_length)
{
    //calculate S(q),  using Teixeira, Eq(15)
    // mathematica query to check limiting conditions:
    //    lim x->0 of [ x gamma(x-1) sin(arctan(q c (x-1))) (q r)^(-x) (1 + 1/(q c)^2)^((1-x)/2) ]
    // Note: gamma(x) may be unreliable for x<0, so the gamma(D-1) is risky.
    // We instead transform D*gamma(D-1) into gamma(D+1)/(D-1).
    double term;
    if (q == 0.) {
        const double D = fractal_dim;
        term = pow(cor_length/radius, D)*sas_gamma(D+1.);
    } else if (fractal_dim == 0.) {
        term = 1.0;
    } else if (fractal_dim == 1.) {
        term = atan(q*cor_length)/(q*radius);
    } else {
        // q>0, D>0
        const double D = fractal_dim;
        const double Dm1 = fractal_dim - 1.0;
        // Note: for large Dm1, sin(Dm1*atan(q*cor_length) can go negative
        const double t1 = sas_gamma(D+1.)/Dm1*sin(Dm1*atan(q*cor_length));
        const double t2 = pow(q*radius, -D);
        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1);
        term = t1 * t2 * t3;
    }
    return 1.0 + term;
}


#endif // SAS_HAVE_fractal_sq


#ifndef SAS_HAVE_fractal_core_shell
#define SAS_HAVE_fractal_core_shell

#line 1 "fractal_core_shell"
static double
form_volume_fractal_core_shell(double radius, double thickness)
{
    return M_4PI_3 * cube(radius + thickness);
}

static double
Iq_fractal_core_shell(double q,
   double radius,
   double thickness,
   double core_sld,
   double shell_sld,
   double solvent_sld,
   double volfraction,
   double fractal_dim,
   double cor_length)
{
    //The radius for the building block of the core shell particle that is
    //needed by the Teixeira fractal S(q) is the radius of the whole particle.
    const double cs_radius = radius + thickness;
    const double sq = fractal_sq(q, cs_radius, fractal_dim, cor_length);
    const double fq = core_shell_fq(q, radius, thickness,
                                    core_sld, shell_sld, solvent_sld);

    return 1.0e-4 * volfraction * sq * fq * fq;
}


#endif // SAS_HAVE_fractal_core_shell



/* END Required header for SASmodel fractal_core_shell */
%}
    DECLARE
%{
  double shape;
  double my_a_v;
%}

INITIALIZE
%{
shape=-1;  /* -1:no shape, 0:cyl, 1:box, 2:sphere  */
if (xwidth && yheight && zdepth)
    shape=1;
  else if (R > 0 && yheight)
    shape=0;
  else if (R > 0 && !yheight)
    shape=2;
  if (shape < 0)
    exit(fprintf(stderr, "SasView_model: %s: sample has invalid dimensions.\n"
                         "ERROR     Please check parameter values.\n", NAME_CURRENT_COMP));

  /* now compute target coords if a component index is supplied */
  if (!target_index && !target_x && !target_y && !target_z) target_index=1;
  if (target_index)
  {
    Coords ToTarget;
    ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
    ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
    coords_get(ToTarget, &target_x, &target_y, &target_z);
  }

  if (!(target_x || target_y || target_z)) {
    printf("SasView_model: %s: The target is not defined. Using direct beam (Z-axis).\n",
      NAME_CURRENT_COMP);
    target_z=1;
  }

  my_a_v = model_abs*2200*100; /* Is not yet divided by v. 100: Convert barns -> fm^2 */

%}


TRACE
%{
  double t0, t1, v, l_full, l, l_1, dt, d_phi, my_s;
  double aim_x=0, aim_y=0, aim_z=1, axis_x, axis_y, axis_z;
  double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z;
  double f, solid_angle, vx_i, vy_i, vz_i, q, qx, qy, qz;
  char intersect=0;

  /* Intersection neutron trajectory / sample (sample surface) */
  if (shape == 0){
    intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, R, yheight);}
  else if (shape == 1){
    intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);}
  else if (shape == 2){
    intersect = sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, R);}
  if(intersect)
  {
    if(t0 < 0)
      ABSORB;

    /* Neutron enters at t=t0. */
    v = sqrt(vx*vx + vy*vy + vz*vz);
    l_full = v * (t1 - t0);          /* Length of full path through sample */
    dt = rand01()*(t1 - t0) + t0;    /* Time of scattering */
    PROP_DT(dt);                     /* Point of scattering */
    l = v*(dt-t0);                   /* Penetration in sample */

    vx_i=vx;
    vy_i=vy;
    vz_i=vz;
    if ((target_x || target_y || target_z)) {
      aim_x = target_x-x;            /* Vector pointing at target (anal./det.) */
      aim_y = target_y-y;
      aim_z = target_z-z;
    }
    if(focus_aw && focus_ah) {
      randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
        aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
    } else if(focus_xw && focus_yh) {
      randvec_target_rect(&vx, &vy, &vz, &solid_angle,
        aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP);
    } else {
      randvec_target_circle(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r);
    }
    NORM(vx, vy, vz);
    vx *= v;
    vy *= v;
    vz *= v;
    qx = V2K*(vx_i-vx);
    qy = V2K*(vy_i-vy);
    qz = V2K*(vz_i-vz);
    q = sqrt(qx*qx+qy*qy+qz*qz);
    
    double trace_radius=radius;
    double trace_thickness=thickness;
    double trace_cor_length=cor_length;
    if ( pd_radius!=0.0 || pd_thickness!=0.0 || pd_cor_length!=0.0 ){
    trace_radius = (randnorm()*pd_radius+1.0)*radius;
    trace_thickness = (randnorm()*pd_thickness+1.0)*thickness;
    trace_cor_length = (randnorm()*pd_cor_length+1.0)*cor_length;
    }

        


    // Sample dependent. Retrieved from SasView./////////////////////
    float Iq_out;
    Iq_out = 1;

    Iq_out = Iq_fractal_core_shell(q, trace_radius, trace_thickness, sld_core, sld_shell, sld_solvent, volfraction, fractal_dim, trace_cor_length);


    float vol;
    vol = 1;

    // Scale by 1.0E2 [SasView: 1/cm  ->   McStas: 1/m]
    Iq_out = model_scale*Iq_out / vol * 1.0E2;

    l_1 = v*t1;
    p *= l_full*solid_angle/(4*PI)*Iq_out*exp(-my_a_v*(l+l_1)/v);
    SCATTER;
  }
%}

MCDISPLAY
%{

  if (shape == 0) {	/* cylinder */
    circle("xz", 0,  yheight/2.0, 0, R);
    circle("xz", 0, -yheight/2.0, 0, R);
    line(-R, -yheight/2.0, 0, -R, +yheight/2.0, 0);
    line(+R, -yheight/2.0, 0, +R, +yheight/2.0, 0);
    line(0, -yheight/2.0, -R, 0, +yheight/2.0, -R);
    line(0, -yheight/2.0, +R, 0, +yheight/2.0, +R);
  }
  else if (shape == 1) { 	/* box */
    double xmin = -0.5*xwidth;
    double xmax =  0.5*xwidth;
    double ymin = -0.5*yheight;
    double ymax =  0.5*yheight;
    double zmin = -0.5*zdepth;
    double zmax =  0.5*zdepth;
    multiline(5, xmin, ymin, zmin,
                 xmax, ymin, zmin,
                 xmax, ymax, zmin,
                 xmin, ymax, zmin,
                 xmin, ymin, zmin);
    multiline(5, xmin, ymin, zmax,
                 xmax, ymin, zmax,
                 xmax, ymax, zmax,
                 xmin, ymax, zmax,
                 xmin, ymin, zmax);
    line(xmin, ymin, zmin, xmin, ymin, zmax);
    line(xmax, ymin, zmin, xmax, ymin, zmax);
    line(xmin, ymax, zmin, xmin, ymax, zmax);
    line(xmax, ymax, zmin, xmax, ymax, zmax);
  }
  else if (shape == 2) {	/* sphere */
    circle("xy", 0,  0.0, 0, R);
    circle("xz", 0,  0.0, 0, R);
    circle("yz", 0,  0.0, 0, R);
  }
%}
END