1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-22 18:28:51 -05:00

Merge pull request #29 from hrosiak/eu

Eu
This commit is contained in:
Andrej Prochazka 2018-03-22 00:50:01 +01:00 committed by GitHub
commit 3394aeb1a2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 446 additions and 17 deletions

View File

@ -2,6 +2,7 @@
#include <algorithm>
#include <cassert>
#include "catima/calculations.h"
#include "catima/build_config.h"
#include "catima/constants.h"
#include "catima/data_ionisation_potential.h"
#include "catima/data_atima.h"

View File

@ -27,11 +27,10 @@
#include "catima/structures.h"
#include "catima/calculations.h"
#include "catima/material_database.h"
#include "catima/storage.h"
namespace catima{
class Interpolator;
/**
* calculate dEdx for projectile-Material combination
* @param p - Projectile

258
spline.cpp Normal file
View File

@ -0,0 +1,258 @@
/*
* This is modification of Tino Kluge tk spline
* https://github.com/ttk592/spline/
*
* the modification is in LU caclulation,
* optimized for tridiagonal matrices
*/
#include "spline.h"
namespace catima{
band_matrix::band_matrix(int dim)
{
resize(dim);
}
void band_matrix::resize(int dim)
{
assert(dim>0);
a.resize(dim);
d.resize(dim);
c.resize(dim);
}
int band_matrix::dim() const
{
return d.size();
}
// defines the new operator (), so that we can access the elements
// by A(i,j), index going from i=0,...,dim()-1
double & band_matrix::operator () (int i, int j)
{
int k=j-i; // what band is the entry
assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
assert(k<2 && k>-2);
if(k>0)return c[i];
else if(k==0) return d[i];
else return a[i];
}
double band_matrix::operator () (int i, int j) const
{
int k=j-i; // what band is the entry
assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
if(k>0)return c[i];
else if(k==0) return d[i];
else return a[i];
}
std::vector<double> band_matrix::trig_solve(const std::vector<double>& b) const
{
assert( this->dim()==(int)b.size() );
std::vector<double> x(this->dim());
std::vector<double> g(this->dim());
int j_stop;
double sum;
assert(d[0]!=0.0);
x[0] = b[0]/d[0];
double bet = d[0];
for(int j=1;j<this->dim();j++){
g[j] = c[j-1]/bet;
bet = d[j] - (a[j]*g[j]);
assert(bet != 0.0);
x[j] = (b[j]-a[j]*x[j-1])/bet;
}
for(int j=this->dim()-2;j>=0;j--){
x[j] -= g[j+1]*x[j+1];
}
return x;
}
// spline implementation
// -----------------------
void spline::set_boundary(spline::bd_type left, double left_value,
spline::bd_type right, double right_value,
bool force_linear_extrapolation)
{
assert(n==0); // set_points() must not have happened yet
m_left=left;
m_right=right;
m_left_value=left_value;
m_right_value=right_value;
m_force_linear_extrapolation=force_linear_extrapolation;
}
void spline::set_points(const double *x,
const double *y,
const size_t num
)
{
assert(num>2);
m_x=x;
m_y=y;
n=num;
// TODO: maybe sort x and y, rather than returning an error
for(int i=0; i<n-1; i++) {
assert(m_x[i]<m_x[i+1]);
}
// setting up the matrix and right hand side of the equation system
// for the parameters b[]
band_matrix A(n);
std::vector<double> rhs(n);
for(int i=1; i<n-1; i++) {
A(i,i-1)=1.0/3.0*(x[i]-x[i-1]);
A(i,i)=2.0/3.0*(x[i+1]-x[i-1]);
A(i,i+1)=1.0/3.0*(x[i+1]-x[i]);
rhs[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
}
// boundary conditions
if(m_left == spline::bd_type::second_deriv) {
// 2*b[0] = f''
A(0,0)=2.0;
A(0,1)=0.0;
rhs[0]=m_left_value;
} else{
// c[0] = f', needs to be re-expressed in terms of b:
// (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
A(0,0)=2.0*(x[1]-x[0]);
A(0,1)=1.0*(x[1]-x[0]);
rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-m_left_value);
}
if(m_right == spline::bd_type::second_deriv) {
// 2*b[n-1] = f''
A(n-1,n-1)=2.0;
A(n-1,n-2)=0.0;
rhs[n-1]=m_right_value;
} else{
// c[n-1] = f', needs to be re-expressed in terms of b:
// (b[n-2]+2b[n-1])(x[n-1]-x[n-2])
// = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
A(n-1,n-1)=2.0*(x[n-1]-x[n-2]);
A(n-1,n-2)=1.0*(x[n-1]-x[n-2]);
rhs[n-1]=3.0*(m_right_value-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
// solve the equation system to obtain the parameters b[]
//m_b=A.lu_solve(rhs);
m_b=A.trig_solve(rhs);
// calculate parameters a[] and c[] based on b[]
m_a.resize(n);
m_c.resize(n);
for(int i=0; i<n-1; i++) {
m_a[i]=1.0/3.0*(m_b[i+1]-m_b[i])/(x[i+1]-x[i]);
m_c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
- 1.0/3.0*(2.0*m_b[i]+m_b[i+1])*(x[i+1]-x[i]);
}
// for left extrapolation coefficients
m_b0 = (m_force_linear_extrapolation==false) ? m_b[0] : 0.0;
m_c0 = m_c[0];
// for the right extrapolation coefficients
// f_{n-1}(x) = b*(x-x_{n-1})^2 + c*(x-x_{n-1}) + y_{n-1}
double h=x[n-1]-x[n-2];
// m_b[n-1] is determined by the boundary condition
m_a[n-1]=0.0;
m_c[n-1]=3.0*m_a[n-2]*h*h+2.0*m_b[n-2]*h+m_c[n-2]; // = f'_{n-2}(x_{n-1})
if(m_force_linear_extrapolation==true)
m_b[n-1]=0.0;
}
double spline::operator() (double x) const
{
assert(n>2);
// find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
auto it=std::lower_bound(m_x,m_x+n,x);
//int idx=std::max( int(it-m_x)-1, 0);
if(it!=m_x)it--;
int idx = std::distance(m_x,it);
double mx = *it;
double h=x-mx;
double interpol;
if(x<m_x[0]) {
// extrapolation to the left
interpol=(m_b0*h + m_c0)*h + m_y[0];
} else if(x>m_x[n-1]) {
// extrapolation to the right
interpol=(m_b[n-1]*h + m_c[n-1])*h + m_y[n-1];
} else {
// interpolation
interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
}
return interpol;
}
double spline::deriv(int order, double x) const
{
assert(order>0);
// find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
auto it=std::lower_bound(m_x,m_x+n,x);
//int idx=std::max( int(it-m_x)-1, 0);
if(it!=m_x)it--;
int idx = std::distance(m_x,it);
double mx = *it;
double h=x-mx;
double interpol;
if(x<m_x[0]) {
// extrapolation to the left
switch(order) {
case 1:
interpol=2.0*m_b0*h + m_c0;
break;
case 2:
interpol=2.0*m_b0*h;
break;
default:
interpol=0.0;
break;
}
} else if(x>m_x[n-1]) {
// extrapolation to the right
switch(order) {
case 1:
interpol=2.0*m_b[n-1]*h + m_c[n-1];
break;
case 2:
interpol=2.0*m_b[n-1];
break;
default:
interpol=0.0;
break;
}
} else {
// interpolation
switch(order) {
case 1:
interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx];
break;
case 2:
interpol=6.0*m_a[idx]*h + 2.0*m_b[idx];
break;
case 3:
interpol=6.0*m_a[idx];
break;
default:
interpol=0.0;
break;
}
}
return interpol;
}
} // namespace tk

93
spline.h Normal file
View File

@ -0,0 +1,93 @@
/*
*
* This is modification of Tino Kluge tk spline
* calculation is optimized for tridiagonal matrices
*
* Copyright(C) 2017
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef TK_SPLINE_H
#define TK_SPLINE_H
#include <cstdio>
#include <cassert>
#include <vector>
#include <algorithm>
namespace catima
{
// band matrix solver
class band_matrix
{
private:
std::vector<double> a;
std::vector<double> d;
std::vector<double> c;
std::vector<double> save;
public:
band_matrix() {}; // constructor
band_matrix(int dim); // constructor
~band_matrix() {}; // destructor
void resize(int dim); // init with dim,n_u,n_l
int dim() const; // matrix dimension
// access operator
double & operator () (int i, int j); // write
double operator () (int i, int j) const; // read
// we can store an additional diogonal (in m_lower)
std::vector<double> trig_solve(const std::vector<double>& b) const;
};
// spline interpolation
class spline
{
public:
enum class bd_type {
first_deriv = 1,
second_deriv = 2
};
private:
const double *m_x, *m_y; // x,y coordinates of points
size_t n=0;
// interpolation parameters
// f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i
std::vector<double> m_a,m_b,m_c; // spline coefficients
double m_b0, m_c0; // for left extrapol
bd_type m_left = bd_type::second_deriv;
bd_type m_right = bd_type::second_deriv;
double m_left_value = 0.0;
double m_right_value = 0.0;
bool m_force_linear_extrapolation = false;
public:
// set default boundary condition to be zero curvature at both ends
spline(){}
// optional, but if called it has to come be before set_points()
void set_boundary(bd_type left, double left_value,
bd_type right, double right_value,
bool force_linear_extrapolation=false);
void set_points(const double *x,
const double *y,
const size_t num);
double operator() (double x) const;
double deriv(int order, double x) const;
};
} // namespace tk
#endif /* TK_SPLINE_H */

View File

@ -1,6 +1,22 @@
/*
* Copyright(C) 2017
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <iostream>
#include "storage.h"
#include "catima/catima.h"
namespace catima {
Data _storage;
EnergyTable<max_datapoints> energy_table(logEmin,logEmax);
@ -58,7 +74,7 @@ DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
}
//////////// Interpolator ////////////////////////////////
Interpolator::Interpolator(const double *x, const double *y, int num,interpolation_t type){
InterpolatorGSL::InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type){
acc = gsl_interp_accel_alloc ();
if(type==cspline)
@ -71,7 +87,7 @@ Interpolator::Interpolator(const double *x, const double *y, int num,interpolati
max= x[num-1];
}
Interpolator::Interpolator(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type){
InterpolatorGSL::InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type){
//Interpolator(x.data(),y.data(),x.size());
acc = gsl_interp_accel_alloc ();
if(type==cspline)
@ -84,21 +100,44 @@ Interpolator::Interpolator(const std::vector<double>& x, const std::vector<doubl
max= x[x.size()-1];
}
Interpolator::~Interpolator(){
InterpolatorGSL::~InterpolatorGSL(){
gsl_interp_accel_free (acc);
gsl_spline_free (spline);
}
double Interpolator::eval(double x){
double InterpolatorGSL::eval(double x){
if(x<min)x=min;
if(x>max)x=max;
return gsl_spline_eval(spline, x, acc);
}
double Interpolator::derivative(double x){
double InterpolatorGSL::derivative(double x){
if(x<min)x=min;
if(x>max)x=max;
return gsl_spline_eval_deriv (spline, x, acc);
}
//////////// Interpolator2 ////////////////////////////////
#ifdef BUILTIN_SPLINE
Interpolator2::Interpolator2(const double *x, const double *y, int num){
ss.set_points(x,y,num);
min= x[0];
max= x[num-1];
}
double Interpolator2::eval(double x){
if(x<min)x=min;
if(x>max)x=max;
return ss(x);
}
double Interpolator2::derivative(double x){
if(x<min)x=min;
if(x>max)x=max;
return ss.deriv(1,x);
}
#endif
}

View File

@ -19,10 +19,17 @@
#include <vector>
#include <iterator>
#include <cmath>
//#include <unordered_set>
#include <gsl/gsl_spline.h>
#include "catima/build_config.h"
#include "catima/constants.h"
#include "catima/catima.h"
#include "catima/structures.h"
#include "catima/config.h"
#ifdef BUILTIN_SPLINE
#include "catima/spline.h"
#endif
namespace catima{
@ -49,20 +56,20 @@ namespace catima{
template<int N>
int EnergyTable_index(const EnergyTable<N> &table, double val){
double lxval = log(val)/M_LN10;
if(val<table.values[0] || val>table.values[table.num-1])return -1;
int i = (int)lxval/table.step;
double lxval = (log(val/table.values[0])/M_LN10);
int i = (int)std::floor(lxval/table.step);
return i;
}
template<int N>
double EnergyTable_interpolate(const EnergyTable<N> &table, double xval, double *y){
double r;
double lxval = log(xval)/M_LN10;
if(xval<table.values[0] || xval>table.values[table.num-1])return 0.0;
if(xval==table.values[table.num-1])return y[table.num-1];
int i = (int)(lxval/table.step);
double linstep = table.values[i+1] - table.values[i];
double lxval = (log(xval/table.values[0])/M_LN10);
int i = (int)std::floor(lxval/table.step);
double linstep = table.values[i+1] - table.values[i];
double x = 1.0 - ((xval - table.values[i])/linstep);
r = (x*y[i]) + ((1-x)*y[i+1]);
return r;
@ -113,11 +120,11 @@ namespace catima{
};
/// Interpolation class, to store interpolated values
class Interpolator{
class InterpolatorGSL{
public:
Interpolator(const double *x, const double *y, int num,interpolation_t type=cspline);
Interpolator(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type=cspline);
~Interpolator();
InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type=cspline);
InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type=cspline);
~InterpolatorGSL();
double operator()(double x){return eval(x);};
double eval(double x);
double derivative(double x);
@ -131,6 +138,22 @@ namespace catima{
gsl_spline *spline;
};
#ifdef BUILTIN_SPLINE
class Interpolator2{
public:
Interpolator2(const double *x, const double *y, int num);
double operator()(double x){return eval(x);};
double eval(double x);
double derivative(double x);
double get_min(){return min;};
double get_max(){return max;};
private:
double min=0;
double max=0;
spline ss;
};
#endif
extern Data _storage;
inline DataPoint& get_data(const Projectile &p, const Material &t, const Config &c=default_config){
@ -138,6 +161,9 @@ namespace catima{
}
bool operator==(const DataPoint &a, const DataPoint &b);
using InterpolatorLinear = InterpolatorGSL;
using Interpolator = InterpolatorGSL;
}
#endif

View File

@ -121,6 +121,19 @@ const lest::test specification[] =
EXPECT(catima::energy_table.values[0]==exp(M_LN10*(catima::logEmin)));
EXPECT(catima::energy_table.values[1]==exp(M_LN10*(catima::logEmin+step)));
EXPECT(catima::energy_table.values[catima::max_datapoints-1]==approx(exp(M_LN10*(catima::logEmax))).epsilon(1e-6));
},
CASE("indexing"){
double step = catima::energy_table.step;
double val, dif;
EXPECT(EnergyTable_index(catima::energy_table, 0.0)==-1);
for(int i:{5,10,100,498}){
val = catima::energy_table.values[i];
dif = catima::energy_table.values[i+1] - val;
EXPECT(EnergyTable_index(catima::energy_table, val)==i);
EXPECT(EnergyTable_index(catima::energy_table, val+0.5*dif)==i);
}
}
};