ein_tab.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "glafic.h"

#define NUM_LNX 801
#define LNX_MIN -23.0
#define DLNX 0.0575 
#define NUM_ALP 99
#define ALP_MIN 0.02
#define DALP 0.01
#define NUM_FLAG 13

double intpol_ein(double tab[NUM_ALP][NUM_LNX], double lnx, double alp);
double intpol_ein_lin(double tab[NUM_ALP][NUM_LNX], double lnx, double alp);

static double tab_alp[NUM_ALP];
static double tab_lnx[NUM_LNX];
static double tab_lnkappa_ein[NUM_ALP][NUM_LNX];
static double tab_lndkappa_ein[NUM_ALP][NUM_LNX];
static double tab_lndphi_ein[NUM_ALP][NUM_LNX];

/*--------------------------------------------------------------
  make tables
*/

void ein_maketable(void)
{
  int i, j, f;
  double x, lx, alp;
  static int flag;

 if(flag !=  NUM_FLAG){
    flag = NUM_FLAG;

    fprintf(stderr, " making tables for ein... [%d x %d = %d points]\n\n", NUM_ALP, NUM_LNX, NUM_ALP*NUM_LNX);

    f = ein_usetab;
    ein_usetab = 0;

    for(i=0;i<NUM_ALP;i++){
      tab_alp[i] = (ALP_MIN) + ((double)i) * (DALP);
    }
    
    for(i=0;i<NUM_LNX;i++){
      tab_lnx[i] = (LNX_MIN) + ((double)i) * (DLNX);
    }
    
    for(i=0;i<NUM_ALP;i++){
      alp = tab_alp[i];
      b_func_ein(1.0, 1.0, alp); /*  set alpha value  */
      for(j=0;j<NUM_LNX;j++){
	lx = tab_lnx[j];
	x = exp(lx);
	tab_lnkappa_ein[i][j] = log(kappa_ein_dl(x) + OFFSET_LOG);
	tab_lndkappa_ein[i][j] = log((-1.0) * dkappa_ein_dl(x) + OFFSET_LOG);
	tab_lndphi_ein[i][j] = log(dphi_ein_dl(x) + OFFSET_LOG);
      }
    }

    ein_usetab = f;
  }

  return;
}

/*--------------------------------------------------------------
  estimate values from the tables
*/

double kappa_ein_dl_tab(double x, double alpha)
{
  double lx, ll;

  lx = log(x);
  ll = intpol_ein(tab_lnkappa_ein, lx, alpha);

  return exp(ll);
}

double dkappa_ein_dl_tab(double x, double alpha)
{
  double lx, ll;

  lx = log(x);
  ll = intpol_ein(tab_lndkappa_ein, lx, alpha);

  return (-1.0) * exp(ll);
}

double dphi_ein_dl_tab(double x, double alpha)
{
  double lx, ll;

  lx = log(x);
  ll = intpol_ein(tab_lndphi_ein, lx, alpha);

  return exp(ll);
}

double intpol_ein(double tab[NUM_ALP][NUM_LNX], double lnx, double alp)
{
  return intpol_ein_lin(tab, lnx, alp);
  /* return intpol_ein_bcu(tab, lnx, alp); */
}

double intpol_ein_lin(double tab[NUM_ALP][NUM_LNX], double lnx, double alp)
{
  int i, j;
  double t, u;

  i = (int)((alp-(ALP_MIN)) / DALP);
  j = (int)((lnx-(LNX_MIN)) / DLNX);

  if(i<0){ i = 0; alp = tab_alp[i]; }
  if(j<0){ j = 0; lnx = tab_lnx[j]; }
  if(i>=NUM_ALP){ i = NUM_ALP-1; alp = tab_alp[i]; }
  if(j>=NUM_LNX){ j = NUM_LNX-1; lnx = tab_lnx[j]; }
  
  t = (alp - tab_alp[i]) / (tab_alp[i + 1] - tab_alp[i]);
  u = (lnx - tab_lnx[j]) / (tab_lnx[j + 1] - tab_lnx[j]);
  
  return (1.0 - t) * (1.0 - u) * tab[i][j] + t * (1.0 - u) * tab[i + 1][j] + t * u * tab[i + 1][j + 1] + (1.0 - t) * u * tab[i][j + 1];
}

#undef NUM_LNX 
#undef LNX_MIN 
#undef DLNX 
#undef NUM_ALP 
#undef ALP_MIN 
#undef DALP 
#undef NUM_FLAG