|
new file 100644
|
|
|
/** @file react_table.c
|
|
|
* @brief Handles the communication between the main code and the reaction tables.
|
|
|
*
|
|
|
* Rather than hardcoding a different function for each possible dependence
|
|
|
* of the reaction-rate on the electric field, these rate/field dependencies
|
|
|
* are written in a table. This allows reactions to be added/changed/removed
|
|
|
* without recompiling the source.
|
|
|
*/
|
|
|
|
|
|
/* Structure of the reaction table file between -----
|
|
|
|
|
|
------
|
|
|
Optional comments, maximum line length of 100 chars!
|
|
|
START
|
|
|
is_log
|
|
|
e_min
|
|
|
e_step
|
|
|
points
|
|
|
underflow
|
|
|
overflow
|
|
|
value(0)
|
|
|
value(1)
|
|
|
...
|
|
|
value(points-1)
|
|
|
-----
|
|
|
|
|
|
Anything before 'START' is ignored. 'START' denotes the start of the actual data
|
|
|
and is case-sensitive. 'is_log' equals one if the electric field values are spaced
|
|
|
logarithmically. 'e_min' denotes the lowest value of the electric field for
|
|
|
which a data-point exists. 'e_step' is the difference in fieldstrength between
|
|
|
consecutive data-points and 'points' is the number of data-points. 'value(x)' are
|
|
|
the actual data-points. Note that anything after 'points' points is ignored.
|
|
|
|
|
|
Underflow and overflow are reaction rates that are returned when the supplied field
|
|
|
strength is out of the bounds of the table.
|
|
|
|
|
|
IMPORTANT: electric field values are 10-logarithmic. So instead of 'e_min' actually
|
|
|
denotes log10(e_min)!
|
|
|
IMPORTANT: electric field values must be (logarithmically) equi-distant! */
|
|
|
|
|
|
/* Example file
|
|
|
|
|
|
-----
|
|
|
Sample reaction table. k(E) = E. 4 points are defined at E = 10^i, i \in {-2,-1,0,1}.
|
|
|
So the 'e_min' value is -2, 'e_step' is 1 and 'points' is 4. Overflow (100) and
|
|
|
underflow (0) values are added in case E-input values fall outside the bounds, in
|
|
|
this case E < 10^-2 or E > 10^1.
|
|
|
START
|
|
|
-2.0
|
|
|
1
|
|
|
4
|
|
|
0.0
|
|
|
100.0
|
|
|
0.01
|
|
|
0.1
|
|
|
1.0
|
|
|
10.0
|
|
|
This bit is ignored by the program, so can be used for comments. Though if i had anything
|
|
|
useful to say, it would've probably been better placed at the start of the file!
|
|
|
-----
|
|
|
*/
|
|
|
|
|
|
#include <stdlib.h>
|
|
|
#include <stdio.h>
|
|
|
#include <math.h>
|
|
|
|
|
|
#include "react_table.h"
|
|
|
|
|
|
/** @brief Reads a reaction rate table from file @a filename and stores the
|
|
|
* table in the form of a 'react_table' at location @a r.
|
|
|
*/
|
|
|
void
|
|
|
react_table_read(char *filename, react_table *r)
|
|
|
{
|
|
|
FILE *fp;
|
|
|
char *comments, buffer[101];
|
|
|
double e_min, e_step, underflow, overflow;
|
|
|
|
|
|
int cnt;
|
|
|
|
|
|
double values[MAX_TABLE_SIZE];
|
|
|
|
|
|
fp = fopen(filename,"r");
|
|
|
if (fp == NULL)
|
|
|
{
|
|
|
fprintf(stderr,"Unable to load reaction table in file: %s -- Shutting down\n",filename);
|
|
|
exit(1);
|
|
|
}
|
|
|
|
|
|
while (strncmp(buffer,"START",5) != 0)
|
|
|
fgets(buffer,100,fp);
|
|
|
|
|
|
fgets(buffer,100,fp);
|
|
|
r->e_min = atof(buffer);
|
|
|
fgets(buffer,100,fp);
|
|
|
r->e_step = atof(buffer);
|
|
|
fgets(buffer,100,fp);
|
|
|
r->steps = atoi(buffer) - 1;
|
|
|
fgets(buffer,100,fp);
|
|
|
r->underflow = atof(buffer);
|
|
|
fgets(buffer,100,fp);
|
|
|
r->overflow = atof(buffer);
|
|
|
|
|
|
for (cnt = 0; cnt <= r->steps; cnt++)
|
|
|
{
|
|
|
fgets(buffer,100,fp);
|
|
|
r->values[cnt] = atof(buffer);
|
|
|
}
|
|
|
|
|
|
fclose(fp);
|
|
|
}
|
|
|
|
|
|
/** @brief Computes an approximation for the reaction rate.
|
|
|
*
|
|
|
* Computes an approximation for the reaction rate at field-strength @a e
|
|
|
* by interpolating reaction rates from lookup table @a r. Stores the result
|
|
|
* in @a ra.
|
|
|
*/
|
|
|
void
|
|
|
react_table_lookup(react_table *r, double e, double *ra)
|
|
|
{
|
|
|
int pos;
|
|
|
double e1, e2, val1, val2, log_e, res;
|
|
|
|
|
|
// e == 0 is possible at initialization.
|
|
|
log_e = (e > 0) ? log10(e) : -1000.0;
|
|
|
|
|
|
/* If the supplied fieldstrength falls outside the boundaries of the table,
|
|
|
return predetermined under-/overflow values. In most cases the underflow(overflow)
|
|
|
will be equal to the lowest(highest) of the specified values in the table. */
|
|
|
if (log_e < r->e_min) { *ra = r->underflow; return; }
|
|
|
if (log_e > r->e_min + r->e_step * r->steps) { *ra = r->overflow; return; }
|
|
|
|
|
|
pos = floor((log_e - r->e_min) / r->e_step);
|
|
|
val1 = r->values[pos];
|
|
|
val2 = r->values[pos+1];
|
|
|
e1 = pow(10, r->e_min + r->e_step * (double) pos);
|
|
|
e2 = e1 * pow(10, r->e_step);
|
|
|
|
|
|
/* Linear interpolation of the 2 table values closest to e.
|
|
|
Note that the field-values in the table are logarithmic! */
|
|
|
|
|
|
res = val1 + ((val2 - val1) / (e2 - e1)) * (e - e1);
|
|
|
|
|
|
/* Dirty tricks with pointers to get the return value right. The normal 'return'
|
|
|
statement would give a completely different value on the "other side".
|
|
|
|
|
|
It's not pretty, but it works, so who cares? :) */
|
|
|
*ra = res;
|
|
|
}
|