// by jcho

#include <math.h>

#include "libfunc.h"

#include <stdio.h>



//used in Si func

int factorial(int i) {

	int k, j;
	if (i<2){
		return 1;
	}

	k=1;

	for(j=1;j<i;j++) {
		k=k*(j+1);
	}

	return k;

}



// Used in pearl nec model

// Sine integral function: approximated within 1%!!!

// integral of sin(x)/x up to namx term nmax=6 looks the best.

double Si(double x)

{
	int i;
	int nmax=6;
	double out;
	long double power;
	double pi = 4.0*atan(1.0);

	if (x >= pi*6.2/4.0){
		double out_sin = 0.0;
		double out_cos = 0.0;
		out = pi/2.0;

		for (i=0; i<nmax-2; i+=1) {
			out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1);
			out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2);
		}

		out -= cos(x) * out_cos;
		out -= sin(x) * out_sin;
		return out;
	}

	out = 0.0;

	for (i=0; i<nmax; i+=1)	{
		if (i==0) {
			out += x;
			continue;
		}

		power = pow(x,(2 * i + 1));
		out += (double)pow(-1, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1));

		//printf ("Si=%g %g %d\n", x, out, i);
	}

	return out;
}



double sinc(double x)
{
  if (x==0.0){
    return 1.0;
  }
  return sin(x)/x;
}

//The code calculating magnetic sld below seems to be not used.
//Keeping it in order to check if it is not breaking anything.
// calculate magnetic sld and return total sld
// bn : contrast (not just sld of the layer)
// m0: max mag of M; mtheta: angle from x-z plane;
// mphi: angle (anti-clock-wise)of x-z projection(M) from x axis
// spinfraci: the fraction of UP among UP+Down (before sample)
// spinfracf: the fraction of UP among UP+Down (after sample and before detector)
// spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis
// Note: all angles are in degrees.
// void cal_msld(polar_sld *p_sld, int isangle, double qx, double qy, double bn,
// 				double m01, double mtheta1, double mphi1,
// 				double spinfraci, double spinfracf, double spintheta)
// {
// 	//locals
// 	double q_x = qx;
// 	double q_y = qy;
// 	double sld = bn;
// 	int is_angle = isangle;
// 	double pi = 4.0*atan(1.0);
// 	double s_theta = spintheta * pi/180.0;
// 	double m_max = m01;
// 	double m_phi = mphi1;
// 	double m_theta = mtheta1;
// 	double in_spin = spinfraci;
// 	double out_spin = spinfracf;
//
// 	double m_perp = 0.0;
// 	double m_perp_z = 0.0;
// 	double m_perp_y = 0.0;
// 	double m_perp_x = 0.0;
// 	double m_sigma_x = 0.0;
// 	double m_sigma_z = 0.0;
// 	double m_sigma_y = 0.0;
// 	//double b_m = 0.0;
// 	double q_angle = 0.0;
// 	double mx = 0.0;
// 	double my = 0.0;
// 	double mz = 0.0;
// 	double uu = sld;
// 	double dd = sld;
// 	double re_ud = 0.0;
// 	double im_ud = 0.0;
// 	double re_du = 0.0;
// 	double im_du = 0.0;
//  const double norm=outspin;


// The norm is needed to make sure that the scattering cross sections are
//correctly weighted, such that the sum of spin-resolved measurements adds up to
// the unpolarised or half-polarised scattering cross section. No intensity weighting
// needed on the incoming polariser side (assuming that a user), has normalised
// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively.

//if (out_spin < 0.5){norm=1-out_spin;}
//else{norm=out_spin;}

// 	//No mag means no further calculation
// 	if (isangle>0) {
// 		if (m_max < 1.0e-32){
// 			uu = sqrt(in_spin * out_spin/ norm) * uu ;
// 			dd = sqrt((1.0 - in_spin) * (1.0 - out_spin)/ norm) * dd ;
// 		}
// 	}
// 	else if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){
// 			uu = sqrt(in_spin * out_spin/ norm) * uu;
// 			dd = sqrt((1.0 - in_spin) * (1.0 - out_spin)/ norm) * dd;
// 	} else {
//
// 		//These are needed because of the precision of inputs
// 		if (in_spin < 0.0) in_spin = 0.0;
// 		if (in_spin > 1.0) in_spin = 1.0;
// 		if (out_spin < 0.0) out_spin = 0.0;
// 		if (out_spin > 1.0) out_spin = 1.0;
//
// 		if (q_x == 0.0) q_angle = pi / 2.0;
// 		else q_angle = atan(q_y/q_x);
// 		if (q_y < 0.0 && q_x < 0.0) q_angle -= pi;
// 		else if (q_y > 0.0 && q_x < 0.0) q_angle += pi;
//
// 		q_angle = pi/2.0 - q_angle;
// 		if (q_angle > pi) q_angle -= 2.0 * pi;
// 		else if (q_angle < -pi) q_angle += 2.0 * pi;
//
// 		if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){
// 			m_perp = 0.0;
// 			}
// 		else {
// 			m_perp = m_max;
// 			}
// 		if (is_angle > 0){
// 			m_phi *= pi/180.0;
// 			m_theta *= pi/180.0;
// 			mx = m_perp * cos(m_theta) * cos(m_phi);
// 			my = m_perp * sin(m_theta);
// 			mz = -(m_perp * cos(m_theta) * sin(m_phi));
// 		}
// 		else{
// 			mx = m_perp;
// 			my = m_phi;
// 			mz = m_theta;
// 		}
// 		//ToDo: simplify these steps
// 		// m_perp1 -m_perp2
// 		m_perp_x = (mx) *  cos(q_angle);
// 		m_perp_x -= (my) * sin(q_angle);
// 		m_perp_y = m_perp_x;
// 		m_perp_x *= cos(-q_angle);
// 		m_perp_y *= sin(-q_angle);
// 		m_perp_z = mz;
//
// 		m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta));
// 		m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta));
// 		m_sigma_z = (m_perp_z);
//
// 		//Find b
// 		uu += m_sigma_x;
// 		dd -= m_sigma_x;
// 		re_ud = m_sigma_y;
// 		re_du = m_sigma_y;
// 		im_ud = m_sigma_z;
// 		im_du = -m_sigma_z;
//
// 		uu = sqrt((in_spin) * (out_spin)/ norm) * uu;
// 		dd = sqrt((1.0 - in_spin) * (1.0 - out_spin)/ norm) * dd;
//
// 		re_ud = sqrt(in_spin * (1.0 - out_spin)/ norm) * re_ud;
// 		im_ud = sqrt(in_spin * (1.0 - out_spin)/ norm) * im_ud;
// 		re_du = sqrt((1.0 - in_spin) * out_spin/ norm) * re_du;
// 		im_du = sqrt((1.0 - in_spin) * out_spin/ norm) * im_du;
// 	}
// 	p_sld->uu = uu;
// 	p_sld->dd = dd;
// 	p_sld->re_ud = re_ud;
// 	p_sld->im_ud = im_ud;
// 	p_sld->re_du = re_du;
// 	p_sld->im_du = im_du;
// }


/** Modifications below by kieranrcampbell@gmail.com
    Institut Laue-Langevin, July 2012
**/

/**
Wojtek's comment Mar 22 2016: The remaing code can mostly likely be deleated
Keeping it in order to check if it is not breaking anything
**/

/*
#define ITMAX 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30

void gser(float *gamser, float a, float x, float *gln) {
  int n;
  float sum,del,ap;

  *gln = lgamma(a);
  if(x <= 0.0) {
    if (x < 0.0) printf("Error: x less than 0 in routine gser");
    *gamser = 0.0;
    return;
  } else {
    ap = a;
    del = sum = 1.0/a;

    for(n=1;n<=ITMAX;n++) {
      ++ap;
      del *= x/ap;
      sum += del;
      if(fabs(del) < fabs(sum)*EPS) {
	*gamser = sum * exp(-x + a * log(x) - (*gln));
	return;
      }
    }
    printf("a too large, ITMAX too small in routine gser");
    return;

  }

}
*/
/**
   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction
   representation
**/
/*
void gcf(float *gammcf, float a, float x, float *gln) {
  int i;
  float an,b,c,d,del,h;

  *gln = lgamma(a);
  b = x+1.0-a;
  c = 1.0/FPMIN;
  d = 1.0/b;
  h=d;
  for (i=1;i <= ITMAX; i++) {
    an = -i*(i-a);
    b += 2.0;
    d = an*d + b;
    if (fabs(d) < FPMIN) d = FPMIN;
    c = b+an/c;
    if (fabs(c) < FPMIN) c = FPMIN;
    d = 1.0/d;
    del = d*c;
    h += del;
    if (fabs(del-1.0) < EPS) break;
  }
  if (i > ITMAX) printf("a too large, ITMAX too small in gcf");
  *gammcf = exp(-x+a*log(x)-(*gln))*h;
  return;
}
*/
/**
   Represents incomplete error function, P(a,x)
**/
/*
float gammp(float a, float x) {
  float gamser,gammcf,gln;
  if(x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammp");
  if (x < (a+1.0)) {
    gser(&gamser,a,x,&gln);
    return gamser;
  } else {
    gcf(&gammcf,a,x,&gln);
    return 1.0 - gammcf;
  }
}
*/
/**
    Implementation of the error function, erf(x)
**/
/*
float erff(float x) {
  return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}
*/