/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SasView_superball
*
* %Identification
* Written by: Jose Robledo
* Based on sasmodels from SasView
* Origin: FZJ / DTU / ESS DMSC
*
*
* SasView superball model component as sample description.
*
* %Description
*
* SasView_superball component, generated from superball.c in sasmodels.
*
* Example: 
*  SasView_superball(sld, sld_solvent, length_a, exponent_p, 
*     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_length_a=0.0)
*
* %Parameters
* INPUT PARAMETERS:
* sld: [1e-6/Ang^2] ([-inf, inf]) Superball scattering length density.
* sld_solvent: [1e-6/Ang^2] ([-inf, inf]) Solvent scattering length density.
* length_a: [Ang] ([0, inf]) Cube edge length of the superball.
* exponent_p: [] ([0, inf]) Exponent describing the roundness of the superball.
* 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_length_a: [] (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_superball

SETTING PARAMETERS (
        sld=4,
        sld_solvent=1,
        length_a=50,
        exponent_p=2.5,
        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_length_a=0.0)


SHARE %{
%include "sas_kernel_header.c"

/* BEGIN Required header for SASmodel superball */
#define HAS_Iqabc
#define HAS_FQ
#define FORM_VOL

#ifndef SAS_HAVE_gauss20
#define SAS_HAVE_gauss20

#line 1 "gauss20"
// Created by Andrew Jackson on 4/23/07

 #ifdef GAUSS_N
 # undef GAUSS_N
 # undef GAUSS_Z
 # undef GAUSS_W
 #endif
 #define GAUSS_N 20
 #define GAUSS_Z Gauss20Z
 #define GAUSS_W Gauss20Wt

// Gaussians
constant double Gauss20Wt[20]={
	.0176140071391521,
	.0406014298003869,
	.0626720483341091,
	.0832767415767047,
	.10193011981724,
	.118194531961518,
	.131688638449177,
	.142096109318382,
	.149172986472604,
	.152753387130726,
	.152753387130726,
	.149172986472604,
	.142096109318382,
	.131688638449177,
	.118194531961518,
	.10193011981724,
	.0832767415767047,
	.0626720483341091,
	.0406014298003869,
	.0176140071391521
};

constant double Gauss20Z[20]={
	-.993128599185095,
	-.963971927277914,
	-.912234428251326,
	-.839116971822219,
	-.746331906460151,
	-.636053680726515,
	-.510867001950827,
	-.37370608871542,
	-.227785851141645,
	-.076526521133497,
	.0765265211334973,
	.227785851141645,
	.37370608871542,
	.510867001950827,
	.636053680726515,
	.746331906460151,
	.839116971822219,
	.912234428251326,
	.963971927277914,
	.993128599185095
};

#pragma acc declare copyin( Gauss20Wt[0:20], Gauss20Z[0:20] )

#endif // SAS_HAVE_gauss20


#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_superball
#define SAS_HAVE_superball

#line 1 "superball"
static double
form_volume_superball(double length_a, double exponent_p)
{
  double g1 = sas_gamma(1.0 / (2.0 * exponent_p));
  double g3 = sas_gamma(3.0 / (2.0 * exponent_p));
  return cube(length_a) / 12.0 / square(exponent_p) * cube(g1) / g3;
}

static double
radius_from_excluded_volume_superball(double length_a, double exponent_p)
{
  double g1 = sas_gamma(1.0 / (2.0 * exponent_p));
  double g3 = sas_gamma(3.0 / (2.0 * exponent_p));
  double g5 = sas_gamma(5.0 / (2.0 * exponent_p));

  return length_a * g3 * sqrt(3.0 / 10.0 / g1 / g5);
}

static double

radius_effective_superball(int mode, double length_a, double exponent_p)

{
  switch (mode)
  {
  default:
  case 1: // radius of gyration
    return radius_from_excluded_volume_superball(length_a, exponent_p);
  case 2: // equivalent volume sphere
    return cbrt(form_volume_superball(length_a, exponent_p) / M_4PI_3);
  case 3: // half length_a
    return 0.5 * length_a;
  }
}

static double oriented_superball(
    double qx,
    double qy,
    double qz,
    double length_a,
    double exponent_p)
{
  // oriented superball form factor

  // outer integral for x
  const double radius = length_a / 2.0; // superball radius
  const double inverse_2p = 1.0 / (2.0 * exponent_p);

  double outer_integral = 0.0; //initialize integral

  for (int i_x = 0; i_x < GAUSS_N; i_x++)
  {
    const double x = 0.5 * (GAUSS_Z[i_x] + 1.0); // integrate 0, 1
    const double x2p = pow(x, 2.0 * exponent_p);
    const double gamma = pow(1.0 - x2p, inverse_2p);

    // inner integral for y
    double inner_integral = 0.0; //initialize integral
    for (int i_y = 0; i_y < GAUSS_N; i_y++)
    {
      const double y = 0.5 * gamma * (GAUSS_Z[i_y] + 1.0); // integrate 0, gamma
      const double y2p = pow(y, 2.0 * exponent_p);
      const double zeta = pow(1.0 - x2p - y2p, inverse_2p);
      const double cos1 = cos(radius * qy * y);
      const double sinc2 = qz == 0 ? radius * zeta : sin(radius * qz * zeta) / qz;
      const double fq = cos1 * sinc2;
      inner_integral += GAUSS_W[i_y] * fq;
    }

    const double co = cos(radius * qx * x);

    // integration factor for -1,1 quadrature to 0, gamma: gamma/2
    const double integration_factor = 0.5 * gamma;

    // Eq. 21 in [Dresen2021]
    outer_integral += GAUSS_W[i_x] * integration_factor * inner_integral * co * 2.0 * square(length_a);

  }
// Needed to normalise the oriented form factor, but would be reverted later with s = SLD contrast * volume
// outer_integral /= form_volume_superball(length_a, exponent_p); 

  // integration factor for -1,1 quadrature to 0, 1: 1/2
  return 0.5 * outer_integral;
}

static void
Fq_superball(double q,
   double *F1,
   double *F2,
   double sld,
   double solvent_sld,
   double length_a,
   double exponent_p)
{

  // translate a point in [-1,1] to a point in [0, pi/2]
  const double zm = M_PI_4;
  const double zb = M_PI_4;

  double orient_averaged_outer_total_F1 = 0.0; //initialize integral
  double orient_averaged_outer_total_F2 = 0.0; //initialize integral
  // phi integral
  for (int i_phi = 0; i_phi < GAUSS_N; i_phi++)
  {

    const double phi = GAUSS_Z[i_phi]*zm +zb; // integrate 0 .. pi/2

    double sin_phi, cos_phi;
    SINCOS(phi, sin_phi, cos_phi);

    double orient_averaged_inner_total_F1 = 0.0; //initialize integral
    double orient_averaged_inner_total_F2 = 0.0; //initialize integral
    // theta integral
    for (int i_theta = 0; i_theta < GAUSS_N; i_theta++)
    {

      const double cos_theta = GAUSS_Z[i_theta]*0.5 + 0.5; // integrate 0, 1
      const double sin_theta = sqrt( 1.0 - square(cos_theta) );


      const double qx = q * cos_phi * sin_theta;
      const double qy = q * sin_phi * sin_theta;
      const double qz = q * cos_theta;

      const double f_oriented = oriented_superball(qx, qy, qz, length_a, exponent_p);


      orient_averaged_inner_total_F1 += GAUSS_W[i_theta] * f_oriented;
      orient_averaged_inner_total_F2 += GAUSS_W[i_theta] * square(f_oriented);

    }
    orient_averaged_outer_total_F1 += GAUSS_W[i_phi] * orient_averaged_inner_total_F1;
    orient_averaged_outer_total_F2 += GAUSS_W[i_phi] * orient_averaged_inner_total_F2;
  }


  // integration factors for phi and theta integral, divided by solid angle of pi/2
  orient_averaged_outer_total_F1 *= 0.25;
  orient_averaged_outer_total_F2 *= 0.25;
  // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
  const double s =  (sld - solvent_sld) ;

  *F1 = 1.0e-2 * s * orient_averaged_outer_total_F1;
  *F2 = 1.0e-4 * s * s * orient_averaged_outer_total_F2;
}

static double
Iqabc_superball(double qa, double qb, double qc,
      double sld,
      double solvent_sld,
      double length_a,
      double exponent_p)
{
  const double f_oriented = oriented_superball(qa, qb, qc, length_a, exponent_p);

  const double s = (sld - solvent_sld); 


  const double form = square(s * f_oriented);
  // Square and convert from [1e-12 A-1] to [cm-1]
  return 1.0e-4 * form;
}


#endif // SAS_HAVE_superball



/* END Required header for SASmodel superball */
%}
    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_length_a=length_a;
    if ( pd_length_a!=0.0 ){
    trace_length_a = (randnorm()*pd_length_a+1.0)*length_a;
    }

        


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

    double F1=0.0, F2=0.0;
    Fq_superball(q, &F1, &F2, sld, sld_solvent, trace_length_a, exponent_p);
    Iq_out = F2;


    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