gsl_complex
gsl_complex_sub (gsl_complex a, gsl_complex b)
{                               /* z=a-b */
  double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  double br = GSL_REAL (b), bi = GSL_IMAG (b);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, ar - br, ai - bi);
  return z;
}

gsl_complex
gsl_complex_add (gsl_complex a, gsl_complex b)
{                               /* z=a+b */
  double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  double br = GSL_REAL (b), bi = GSL_IMAG (b);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, ar + br, ai + bi);
  return z;
}

 /** Subtraction operator (vector) */
 vector<complex> vector<complex>::operator-(const vector<complex>& v) const
 {
     vector<complex> v1(_vector);
     gsl_complex z1;
     GSL_SET_COMPLEX(&z1, -1., 0.);
     if (gsl_blas_zaxpy(z1, v.as_gsl_type_ptr(), v1.as_gsl_type_ptr())) {
         std::cout << "\n Error in vector<complex> -" << std::endl;
         exit(EXIT_FAILURE);
     }
     return v1;
 }

/* ------------------------------------------------------ */
gsl_complex gsl_complex_min_real (gsl_complex a, gsl_complex b)
{
  gsl_complex z;
  double min;
	
  /* just consider real parts */
  min = GSL_REAL(a) < GSL_REAL(b) ? GSL_REAL(a) : GSL_REAL(b);
  GSL_SET_COMPLEX (&z, min, 0);
	
  return z;
}

gsl_complex
gsl_complex_cos (gsl_complex a)
{				/* z = cos(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);

  gsl_complex z;

  if (I == 0.0) 
    {
      /* avoid returing negative zero (-0.0) for the imaginary part  */

      GSL_SET_COMPLEX (&z, cos (R), 0.0);  
    } 
  else 
    {
      GSL_SET_COMPLEX (&z, cos (R) * cosh (I), sin (R) * sinh (-I));
    }

  return z;
}

gsl_complex
gsl_complex_pow_real (gsl_complex a, double b)
{				/* z=a^b */
  gsl_complex z;

  if (GSL_REAL (a) == 0 && GSL_IMAG (a) == 0)
    {
      GSL_SET_COMPLEX (&z, 0, 0);
    }
  else
    {
      double logr = gsl_complex_logabs (a);
      double theta = gsl_complex_arg (a);
      double rho = exp (logr * b);
      double beta = theta * b;
      GSL_SET_COMPLEX (&z, rho * cos (beta), rho * sin (beta));
    }

  return z;
}

/* ------------------------------------------------------ */
gsl_complex gsl_complex_max_real (gsl_complex a, gsl_complex b)
{
  gsl_complex z;
  double max;
	
  /* just consider real parts */
  max = GSL_REAL(a) > GSL_REAL(b) ? GSL_REAL(a) : GSL_REAL(b);
  GSL_SET_COMPLEX (&z, max, 0);
	
  return z;
}

complex_spinor complex_spinor::operator+(const complex_spinor & complex_spinor_param)const{

	gsl_complex aux_UP;
	gsl_complex aux_DOWN;

	GSL_SET_COMPLEX(&aux_UP, 0, 0);
	GSL_SET_COMPLEX(&aux_DOWN, 0, 0);
	
	int i;	
	int nSitios = this->_numSites;
	complex_spinor complex_spinor_res(nSitios);
	
	for(i=0; i<nSitios ; i++){
	aux_UP = gsl_complex_add(complex_spinor_param.complex_spinor_get(i, UP) , this->complex_spinor_get(i, UP)); 
  	aux_DOWN = gsl_complex_add(complex_spinor_param.complex_spinor_get(i, DOWN) , this->complex_spinor_get(i,DOWN));
	complex_spinor_res.complex_spinor_set(i, UP, aux_UP);
	complex_spinor_res.complex_spinor_set(i, DOWN, aux_DOWN);
	}
	return complex_spinor_res;
}

gsl_complex gsl_complex_pow(gsl_complex a, gsl_complex b)
{				/* z=a^b */
    gsl_complex z;

    if (GSL_REAL(a) == 0 && GSL_IMAG(a) == 0.0) {
	GSL_SET_COMPLEX(&z, 0.0, 0.0);
    } else {
	double logr = gsl_complex_logabs(a);
	double theta = gsl_complex_arg(a);

	double br = GSL_REAL(b), bi = GSL_IMAG(b);

	double rho = exp(logr * br - bi * theta);
	double beta = theta * br + bi * logr;

	GSL_SET_COMPLEX(&z, rho * cos(beta), rho * sin(beta));
    }

    return z;
}

gsl_complex
gsl_complex_sqrt (gsl_complex a)
{				/* z=sqrt(a) */
  gsl_complex z;

  if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0)
    {
      GSL_SET_COMPLEX (&z, 0, 0);
    }
  else
    {
      double x = fabs (GSL_REAL (a));
      double y = fabs (GSL_IMAG (a));
      double w;

      if (x >= y)
	{
	  double t = y / x;
	  w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + t * t)));
	}
      else
	{
	  double t = x / y;
	  w = sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 + t * t)));
	}

      if (GSL_REAL (a) >= 0.0)
	{
	  double ai = GSL_IMAG (a);
	  GSL_SET_COMPLEX (&z, w, ai / (2.0 * w));
	}
      else
	{
	  double ai = GSL_IMAG (a);
	  double vi = (ai >= 0) ? w : -w;
	  GSL_SET_COMPLEX (&z, ai / (2.0 * vi), vi);
	}
    }

  return z;
}

static int
do_ccmul(lua_State *L, mMatComplex *a, double b_re, double b_im)
{
    mMatComplex *r = qlua_newMatComplex(L, a->l_size, a->r_size);
    gsl_complex z;

    gsl_matrix_complex_memcpy(r->m, a->m);
    GSL_SET_COMPLEX(&z, b_re, b_im);
    gsl_matrix_complex_scale(r->m, z);

    return 1;
}

 void
 complex::assign(const double& real=0., const double& imag=0.,
                 bool polar=false)
 {
   if (polar)
     _complex = gsl_complex_polar(real,imag);
   else
     {
       GSL_SET_COMPLEX(&_complex, real, imag);
       // *_complex = gsl_complex_rect(real,imag);
     }
 }

static int
cm_mul_cm(lua_State *L)
{
    mMatComplex *a = qlua_checkMatComplex(L, 1);
    mMatComplex *b = qlua_checkMatComplex(L, 2);
    int al = a->l_size;
    int ar = a->r_size;
    int bl = b->l_size;
    int br = b->r_size;
    mMatComplex *r = qlua_newMatComplex(L, al, br);
    gsl_complex z1, z0;

    if (ar != bl)
        return luaL_error(L, "matrix sizes mismatch in m * m");

    GSL_SET_COMPLEX(&z1, 1.0, 0.0);
    GSL_SET_COMPLEX(&z0, 0.0, 0.0);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, z1, a->m, b->m, z0, r->m);

    return 1;
}

void CModel::AllFacetHarmonics(int bw, double r, gsl_complex *coeff){
  for(int i=0; i<bw*bw; i++) GSL_SET_COMPLEX(&coeff[i], 0.0, 0.0);
  
  gsl_complex *temp = new gsl_complex[bw*bw];
  
  for(int i=0; i<f_no; i++){
    FacetHarmonics(i, bw, r, temp);
    for(int j=0; j<bw*bw; j++) coeff[j] = gsl_complex_add(coeff[j],temp[j]);
  }
	
  delete [] temp;
}

gsl_complex gsl_complex_arccos(gsl_complex a)
{				/* z = arccos(a) */
    double R = GSL_REAL(a), I = GSL_IMAG(a);
    gsl_complex z;

    if (I == 0) {
	z = gsl_complex_arccos_real(R);
    } else {
	double x = fabs(R), y = fabs(I);
	double r = hypot(x + 1, y), s = hypot(x - 1, y);
	double A = 0.5 * (r + s);
	double B = x / A;
	double y2 = y * y;

	double real, imag;

	const double A_crossover = 1.5, B_crossover = 0.6417;

	if (B <= B_crossover) {
	    real = acos(B);
	} else {
	    if (x <= 1) {
		double D =
		    0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));
		real = atan(sqrt(D) / x);
	    } else {
		double Apx = A + x;
		double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));
		real = atan((y * sqrt(D)) / x);
	    }
	}

	if (A <= A_crossover) {
	    double Am1;

	    if (x < 1) {
		Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));
	    } else {
		Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));
	    }

	    imag = log1p(Am1 + sqrt(Am1 * (A + 1)));
	} else {
	    imag = log(A + sqrt(A * A - 1));
	}

	GSL_SET_COMPLEX(&z, (R >= 0) ? real : M_PI - real,
			(I >= 0) ? -imag : imag);
    }

    return z;
}

static VALUE Vector_complex_set_all(VALUE obj, VALUE va) {
  gsl_vector_complex * v;
  gsl_complex z;
  double d0, d1;
  d0 = NUM2DBL(rb_ary_entry(va,0));
  d1 = NUM2DBL(rb_ary_entry(va,1));
  GSL_SET_COMPLEX(&z, d0, d1);
  Data_Get_Struct(obj, gsl_vector_complex, v);

  gsl_vector_complex_set_all(v, z);

  return obj;
}

gsl_complex gsl_complex_arccot(gsl_complex a)
{				/* z = arccot(a) */
    gsl_complex z;

    if (GSL_REAL(a) == 0.0 && GSL_IMAG(a) == 0.0) {
	GSL_SET_COMPLEX(&z, M_PI_2, 0);
    } else {
	z = gsl_complex_inverse(a);
	z = gsl_complex_arctan(z);
    }

    return z;
}

//go through all sets of Bethe rapidities and calculate g1 for all identical sets. return value is sum of all these terms.
gsl_complex exact_integrationDE(double orteins, double ortzwei, double** kEp_plusminus, int anzahl_sets, double** coeffoverlap)
{

  //variables to make loop easily readable, could delete from final version
  double c, Energie, Energieprime, normierung, normierungprime;
  double k[N];
  double kprime[N];//corresponds to primed Bethe rapidity set in papers
  gsl_complex overlap, overlapprime, integral;
  GSL_SET_COMPLEX(&integral, 0.0, 0.0);//sum of all values, technically no needed but good for quick check of initial shape of g1
  
  int zaehler_integrale = 0;//counter of total sets of integrals
  
  double Nminusone_fak = 1.0;//(N-1) factorial, see definition of integrals on restricted spatial domain
  for(int j=1; j<=(N-1); j++)
    Nminusone_fak = Nminusone_fak * ((double) j);


  for(int zaehler_eigenstate = 0 ; zaehler_eigenstate < anzahl_sets; zaehler_eigenstate++){//loop over all Bethe rapidity sets
    
    c = kEp_plusminus[zaehler_eigenstate][0];//interaction strength
    for(int bb=0; bb<N; bb++)
      k[bb]=kEp_plusminus[zaehler_eigenstate][bb+1];//Bethe rapidities
    Energie = kEp_plusminus[zaehler_eigenstate][N+1];//energy
    normierung=kEp_plusminus[zaehler_eigenstate][N+6];//norm    
    GSL_SET_COMPLEX(&overlap, coeffoverlap[zaehler_eigenstate][0], coeffoverlap[zaehler_eigenstate][1]);//overlap of set with initial state

      //calculate g1_DE for this Bethe set
    gsl_complex integralwert = integraleigenfunction(k, k, c, normierung, normierung, orteins, ortzwei);
    
    integralwert = gsl_complex_mul_real(integralwert, (Nminusone_fak*Laenge));//integral evaluated in Rp only -> missing factor of (N-1)!. Laenge (=system length): definition of g1 has prefactor of N!/(n*(N-1)!) = Laenge

    integral = gsl_complex_add(integral, gsl_complex_mul(gsl_complex_mul(overlap, gsl_complex_conjugate(overlap)), integralwert));
     
  }//end for loop zaehler_eigenstate 

  
  return integral;
};

gsl_complex
gsl_complex_arccos_real (double a)
{				/* z = arccos(a) */
  gsl_complex z;

  if (fabs (a) <= 1.0)
    {
      GSL_SET_COMPLEX (&z, acos (a), 0);
    }
  else
    {
      if (a < 0.0)
	{
	  GSL_SET_COMPLEX (&z, M_PI, -acosh (-a));
	}
      else
	{
	  GSL_SET_COMPLEX (&z, 0, acosh (a));
	}
    }

  return z;
}

gsl_complex
gsl_complex_arccsc_real (double a)
{				/* z = arccsc(a) */
  gsl_complex z;

  if (a <= -1.0 || a >= 1.0)
    {
      GSL_SET_COMPLEX (&z, asin (1 / a), 0.0);
    }
  else
    {
      if (a >= 0.0)
	{
	  GSL_SET_COMPLEX (&z, M_PI_2, -acosh (1 / a));
	}
      else
	{
	  GSL_SET_COMPLEX (&z, -M_PI_2, -acosh (-1 / a));
	}
    }

  return z;
}