# Input is :
# First line = appropriate line from table29
# Next line(s) = temperature (if none, tabulate the ionization rate)
# Note: first exponential integral computed with polynomial approxomation
#       from J.Jin (http://jin.ece.uiuc.edu/routines/routines.html)

function E1(x) {	# Compute first exponential integral
    #print "#Compute E1(" x ")";
    if (x<0) return(0./0.);	# NaN
    if (x == 0) return(1./0.);	# Inf
    if (x < 1) return ( -log(x)+  \
       ((((1.07857e-3*x-9.76004e-3)*x+5.519968e-2)*x -0.24991055)*x \
        +0.99999193)*x-0.57721566);
    es1 = (((x+8.5733287401)*x+18.059016973)*x  +8.6347608925)*x+0.2677737343;
    es2 = (((x+9.5733223454)*x+25.6329561486)*x+21.0996530827)*x+3.9584969228;
    return(exp(-x)/x*es1/es2);
}

function title() {
    print "#======================";
    printf "#%9s %-10s# \n", " T[K]   ", "R(T)[cm3/s]"
    print "#======================";
}

BEGIN {
    k = 8.617e-5;	# Boltzmann in eV/K
    f = 2.
}

/^#/  { print; next }
/^ *$/ { next }

# First Line contains the ion, next lines contain a temperature
{ 
   if (Z<1) {		# First line = Z Ion Nspl I Tmin xi {rho}i
       Z = $1; Ion = $2; Ns = $3+0; I = $4+0; Tmin = $5+0; 
       for (i=1; i<=Ns; i++) X[i] = $(i+5)+0;
       for (i=1; i<=Ns; i++) Y[i] = $(i+5+Ns)*1.e-6;
       print "# Element=" Z "." Ion " I=" I "eV Tmin=" Tmin " Ns=" Ns;
       next
   }

   #print "#Input=" $0;

   if (nT == 0) title();
   nT += 1;		# Number of temperature(s) computed.
   T = $1+0;
   if (T < Tmin) { 
       print "#***T too small:", T " < " Tmin;
       next
   }
   t = k*T/I; 
   x = 1.0 - log(f)/log(t+f);

   #printf "#T=%-10.3e t=%-10.3e x=%-7.4f\n", T, t, x;

   # Find the point X[Is] w<hich is closest to x.
   Dx = 9.9; Is = 0;
   for (i=1; i<=Ns; i++) {
       dx = x - X[i]; if (dx < 0) dx = -dx;
       if (dx < Dx) { Is = i; Dx = dx; }
   }

   # Interpolate
   if (Is < Ns) Is1 = Is+1; else Is1 = Is-1;

   dx = (x - X[Is])/(X[Is1]-X[Is]);
    y = Y[Is] + (Y[Is1]-Y[Is])*dx;

    R = E1(1./t);
    R *= y / sqrt(t*I) / I ;

    printf "%-10.3e %-10.3e # x=%-7.4f y=%-7.4f from point#%-2d", T,R, \
           x,y*1.e6, Is;
    printf " (%7.4f,%7.4f)\n", X[Is], Y[Is]*1.e+6;
}

END { 
    if (nT < 1) {	# Compute from Tmin to 10+9
        title();
        for (i=1; i<=Ns; i++) {
	    x = X[i]; y = Y[i];
	    t = exp(log(f)/(1-x))-f;
	    T = t*(I/k);
	    R = E1(1./t);
	    R *= y / sqrt(t*I) / I ;
	    printf "%-10.3e %-10.3e # x=%-7.4f y=%-7.4f from point#%-2d\n", \
	       T, R, x, y*1.e6, i;
	}
    }
    else print "#Points computed:", nT 
}