/************************************************************

  elliptic_utils.c

  Various routines used for implementing
  the elliptic representation of the
  Boltzmann equation. Contains:

  get_X_gamm(Xs) = X/gamma(Xs)
  get_Pterm(Xs) = (1/2)(1-X/gamma)(Xs)
  get_Vterm(Xs) = (3(X/gamma)-1)/(2Xs)

  where Xs = X*X

  Created 2001 by Edward A. Richley

  This file is subject to the terms and conditions of the GNU General
  Public License. 

************************************************************/
#include <stdlib.h>
#include <math.h>

#define NG 50
#define EPS_MIN 1e-7
#define EPS_ANG 1e-9
#define EPS_NEWT 1e-11

static double get_g(double);
static double get_gamma(double);
double get_Pterm(double);
double get_Vterm(double);
double get_X_Gamma(double);
/************************************************************

  get_X_gamma(Xs)

  called with Xs = X*X

************************************************************/
double get_X_gamma(Xs)
double Xs;
{
  static int First_Time = 1;
  static double X_g[NG+1], X_gprime[NG+1],dN;
  static double a[NG+1], b[NG+1], c[NG+1], d[NG+1];
  int i,j;
  double xs, x, x_g, g, gprime;

  if (First_Time){
    dN = (double)NG;
    for(i=0; i<=NG; i++){
      xs = (double)i/dN;
      x = sqrt(xs);
      g = get_gamma(x);
      if (x < EPS_MIN){
	gprime = 3.0 ; 
	X_g[i] = .33333333333333333333333;
	X_gprime[i] = 0.8;
      } else {
	X_g[i] = x_g = x/g;
	gprime = (x < (1.0-EPS_MIN))?
	  ((1-g*g)/(1 + x*x - 2*x_g)):0.0;
	X_gprime[i] = .5*(1-x_g*gprime)/(g*x);
      }
      X_gprime[i] /= 2*dN;
    }
    for(i=0; i<NG; i++){
      j = i+1;
      a[i] = .25*(X_gprime[i]+X_gprime[j])+.25*(X_g[i]-X_g[j]);
      b[i] = .25*(X_gprime[j]-X_gprime[i]);
      c[i] = -.25*(X_gprime[i]+X_gprime[j])-.75*(X_g[i]-X_g[j]);
      d[i] = .5*(X_g[i]+X_g[j])+.25*(X_gprime[i]-X_gprime[j]);
    }
    First_Time = 0;
  }
  i = (int)(dN*Xs);
  if (i>=NG) {
    g = .5*(Xs+1);
    return(g);
  }
  j = i+1;
  x = (Xs+Xs)*dN-(double)(i+j);
  g = d[i]+x*(c[i]+x*(b[i]+x*a[i]));
  return(g);
}
/************************************************************

  get_Pterm(Xs)

  called with Xs = X*X

  Pterm = (1/2)(1-X/gamma)

************************************************************/
double get_Pterm(Xs)
double Xs;
{
  static int First_Time = 1;
  static double Pterm[NG+1], Ptermprime[NG+1],dN;
  static double a[NG+1], b[NG+1], c[NG+1], d[NG+1];
  int i,j;
  double x, xs, x_g, g, gprime;

  if (First_Time){
    dN = (double)NG;
    for(i=0; i<=NG; i++){
      xs = (double)i/dN;
      x = sqrt(xs);
      g = get_gamma(x);
      if (x < EPS_MIN){
	gprime = 3.0 ; 
	Pterm[i] = .33333333333333333333333;
	Ptermprime[i] = -0.4;
      } else {
	x_g = x/g;
	gprime = (x < (1.0-EPS_MIN))?
	  ((1-g*g)/(1 + x*x -2*x_g)):0.0;
	Pterm[i] = .5*(1-x_g);
	Ptermprime[i] = .25*(x_g*gprime-1)/(x*g);
      }
      Ptermprime[i] /= 2*dN;
    }
    for(i=0; i<NG; i++){
      j = i+1;
      a[i] = .25*(Ptermprime[i]+Ptermprime[j])+.25*(Pterm[i]-Pterm[j]);
      b[i] = .25*(Ptermprime[j]-Ptermprime[i]);
      c[i] = -.25*(Ptermprime[i]+Ptermprime[j])-.75*(Pterm[i]-Pterm[j]);
      d[i] = .5*(Pterm[i]+Pterm[j])+.25*(Ptermprime[i]-Ptermprime[j]);
    }
    First_Time = 0;
  }
  i = (int)(dN*Xs);
  if (i>=NG) {
    /*    g = .25*(1-Xs);*/
    g = 0.0;
    return(g);
  }
  j = i+1;
  x = (Xs+Xs)*dN-(double)(i+j);
  g = d[i]+x*(c[i]+x*(b[i]+x*a[i]));
  return(g);
}
/************************************************************

  get_Vterm(Xs)

  called with Xs = X*X

  Vterm = (3(X/gamma)-1)/(2X*X)

************************************************************/
double get_Vterm(Xs)
double Xs;
{
  static int First_Time = 1;
  static double Vterm[NG+1], Vtermprime[NG+1],dN;
  static double a[NG+1], b[NG+1], c[NG+1], d[NG+1];
  int i,j;
  double x, xs, x_g, g, gprime;

  if (First_Time){
    dN = (double)NG;
    for(i=0; i<=NG; i++){
      xs = (double)i/dN;
      x = sqrt(xs);
      g = get_gamma(x);
      if (x < EPS_MIN){
	gprime = 3.0 ; 
	Vterm[i] = 6.0/5.0;
	Vtermprime[i] = -18.0/(165.0);
      } else {
	x_g = x/g;
	gprime = (x < (1.0-EPS_MIN))?
	  ((1-g*g)/(1 + x*x - 2*x_g)): 0.0;
	Vterm[i] = (1.5/g-.5/x)/x;
	Vtermprime[i] = .5*(1/(xs*xs)-1.5*(1/x+gprime/g)/(g*xs));
      }
      Vtermprime[i] /= 2*dN;
    }
    for(i=0; i<NG; i++){
      j = i+1;
      a[i] = .25*(Vtermprime[i]+Vtermprime[j])+.25*(Vterm[i]-Vterm[j]);
      b[i] = .25*(Vtermprime[j]-Vtermprime[i]);
      c[i] = -.25*(Vtermprime[i]+Vtermprime[j])-.75*(Vterm[i]-Vterm[j]);
      d[i] = .5*(Vterm[i]+Vterm[j])+.25*(Vtermprime[i]-Vtermprime[j]);
    }
    First_Time = 0;
  }
  i = (int)(dN*Xs);
  if (i>=NG) {
    /*    g = 1.25-.25*Xs;*/
    g = 1.0;
    return(g);
  }
  j = i+1;
  x = (Xs+Xs)*dN-(double)(i+j);
  g = d[i]+x*(c[i]+x*(b[i]+x*a[i]));
  return(g);
}
/************************************************************

  get_gamma(X)

************************************************************/
static double get_gamma(X)
double X;
{
  double f, fp, g, del, eps, y;
  
  g = X;
  eps = 1.0;
  if (X < EPS_MIN) {return(3.0*X);}
  if (fabs(X - 1.0) < EPS_MIN) {return(1.0);}
  while(eps > EPS_NEWT){
    y = log((1+g)/(1-g));
    f = X-1/g+2.0/y;
    fp = -1/(g*g)-4.0/((1-g*g)*y*y);
    del = -f/fp;
    eps = (fabs(del/g));
    while ((g + del) >= 1.0) {
      del /= 2.0;
    }
    g = g + del;      
  }
  return(g);
}
/************************************************************

  get_g(X)

************************************************************/
static double get_g(X)
double X;
{
  static int First_Time = 1;
  static double gamma[NG+1], gprime[NG+1],dN;
  static double a[NG+1], b[NG+1], c[NG+1], d[NG+1];
  int i,j;
  double x, g, di;

  if (First_Time){
    dN = (double)NG;
    for(i=0; i<=NG; i++){
      x = (double)i/dN;
      gamma[i] = g = get_gamma(x);
      gprime[i] = (g < EPS_MIN)?
	3.0:
	(g < (1.0-EPS_MIN))?
	(g*(1-g*g)/(g-2*x+g*x)):
	0.0;
      gprime[i] /= 2*dN;
    }
    for(i=0; i<NG; i++){
      j = i+1;
      a[i] = .25*(gprime[i]+gprime[j])+.25*(gamma[i]-gamma[j]);
      b[i] = .25*(gprime[j]-gprime[i]);
      c[i] = -.25*(gprime[i]+gprime[j])-.75*(gamma[i]-gamma[j]);
      d[i] = .5*(gamma[i]+gamma[j])+.25*(gprime[i]-gprime[j]);
    }
    First_Time = 0;
  }
  i = (int)(dN*X);
  if (i >= NG){
    g = 1.0;
    return(g);
  }
  j = i+1;
  x = (X+X)*dN-(double)(i+j);
  g = d[i]+x*(c[i]+x*(b[i]+x*a[i]));
  return(g);
}