|
Packit |
67cb25 |
/* specfunc/synchrotron.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_exp.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_pow_int.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_synchrotron.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "chebyshev.h"
|
|
Packit |
67cb25 |
#include "cheb_eval.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double synchrotron1_data[13] = {
|
|
Packit |
67cb25 |
30.364682982501076273,
|
|
Packit |
67cb25 |
17.079395277408394574,
|
|
Packit |
67cb25 |
4.560132133545072889,
|
|
Packit |
67cb25 |
0.549281246730419979,
|
|
Packit |
67cb25 |
0.372976075069301172e-01,
|
|
Packit |
67cb25 |
0.161362430201041242e-02,
|
|
Packit |
67cb25 |
0.481916772120371e-04,
|
|
Packit |
67cb25 |
0.10512425288938e-05,
|
|
Packit |
67cb25 |
0.174638504670e-07,
|
|
Packit |
67cb25 |
0.22815486544e-09,
|
|
Packit |
67cb25 |
0.240443082e-11,
|
|
Packit |
67cb25 |
0.2086588e-13,
|
|
Packit |
67cb25 |
0.15167e-15
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series synchrotron1_cs = {
|
|
Packit |
67cb25 |
synchrotron1_data,
|
|
Packit |
67cb25 |
12,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
9
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double synchrotron2_data[12] = {
|
|
Packit |
67cb25 |
0.4490721623532660844,
|
|
Packit |
67cb25 |
0.898353677994187218e-01,
|
|
Packit |
67cb25 |
0.81044573772151290e-02,
|
|
Packit |
67cb25 |
0.4261716991089162e-03,
|
|
Packit |
67cb25 |
0.147609631270746e-04,
|
|
Packit |
67cb25 |
0.3628633615300e-06,
|
|
Packit |
67cb25 |
0.66634807498e-08,
|
|
Packit |
67cb25 |
0.949077166e-10,
|
|
Packit |
67cb25 |
0.1079125e-11,
|
|
Packit |
67cb25 |
0.10022e-13,
|
|
Packit |
67cb25 |
0.77e-16,
|
|
Packit |
67cb25 |
0.5e-18
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series synchrotron2_cs = {
|
|
Packit |
67cb25 |
synchrotron2_data,
|
|
Packit |
67cb25 |
11,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
7
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double synchrotron1a_data[23] = {
|
|
Packit |
67cb25 |
2.1329305161355000985,
|
|
Packit |
67cb25 |
0.741352864954200240e-01,
|
|
Packit |
67cb25 |
0.86968099909964198e-02,
|
|
Packit |
67cb25 |
0.11703826248775692e-02,
|
|
Packit |
67cb25 |
0.1645105798619192e-03,
|
|
Packit |
67cb25 |
0.240201021420640e-04,
|
|
Packit |
67cb25 |
0.35827756389389e-05,
|
|
Packit |
67cb25 |
0.5447747626984e-06,
|
|
Packit |
67cb25 |
0.838802856196e-07,
|
|
Packit |
67cb25 |
0.13069882684e-07,
|
|
Packit |
67cb25 |
0.2053099071e-08,
|
|
Packit |
67cb25 |
0.325187537e-09,
|
|
Packit |
67cb25 |
0.517914041e-10,
|
|
Packit |
67cb25 |
0.83002988e-11,
|
|
Packit |
67cb25 |
0.13352728e-11,
|
|
Packit |
67cb25 |
0.2159150e-12,
|
|
Packit |
67cb25 |
0.349967e-13,
|
|
Packit |
67cb25 |
0.56994e-14,
|
|
Packit |
67cb25 |
0.9291e-15,
|
|
Packit |
67cb25 |
0.152e-15,
|
|
Packit |
67cb25 |
0.249e-16,
|
|
Packit |
67cb25 |
0.41e-17,
|
|
Packit |
67cb25 |
0.7e-18
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series synchrotron1a_cs = {
|
|
Packit |
67cb25 |
synchrotron1a_data,
|
|
Packit |
67cb25 |
22,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
11
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double synchrotron21_data[13] = {
|
|
Packit |
67cb25 |
38.617839923843085480,
|
|
Packit |
67cb25 |
23.037715594963734597,
|
|
Packit |
67cb25 |
5.3802499868335705968,
|
|
Packit |
67cb25 |
0.6156793806995710776,
|
|
Packit |
67cb25 |
0.406688004668895584e-01,
|
|
Packit |
67cb25 |
0.17296274552648414e-02,
|
|
Packit |
67cb25 |
0.51061258836577e-04,
|
|
Packit |
67cb25 |
0.110459595022e-05,
|
|
Packit |
67cb25 |
0.18235530206e-07,
|
|
Packit |
67cb25 |
0.2370769803e-09,
|
|
Packit |
67cb25 |
0.24887296e-11,
|
|
Packit |
67cb25 |
0.21529e-13,
|
|
Packit |
67cb25 |
0.156e-15
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series synchrotron21_cs = {
|
|
Packit |
67cb25 |
synchrotron21_data,
|
|
Packit |
67cb25 |
12,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
9
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double synchrotron22_data[13] = {
|
|
Packit |
67cb25 |
7.9063148270660804288,
|
|
Packit |
67cb25 |
3.1353463612853425684,
|
|
Packit |
67cb25 |
0.4854879477453714538,
|
|
Packit |
67cb25 |
0.394816675827237234e-01,
|
|
Packit |
67cb25 |
0.19661622334808802e-02,
|
|
Packit |
67cb25 |
0.659078932293042e-04,
|
|
Packit |
67cb25 |
0.15857561349856e-05,
|
|
Packit |
67cb25 |
0.286865301123e-07,
|
|
Packit |
67cb25 |
0.4041202360e-09,
|
|
Packit |
67cb25 |
0.45568444e-11,
|
|
Packit |
67cb25 |
0.420459e-13,
|
|
Packit |
67cb25 |
0.3232e-15,
|
|
Packit |
67cb25 |
0.21e-17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series synchrotron22_cs = {
|
|
Packit |
67cb25 |
synchrotron22_data,
|
|
Packit |
67cb25 |
12,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
8
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double synchrotron2a_data[17] = {
|
|
Packit |
67cb25 |
2.020337094170713600,
|
|
Packit |
67cb25 |
0.10956237121807404e-01,
|
|
Packit |
67cb25 |
0.8542384730114676e-03,
|
|
Packit |
67cb25 |
0.723430242132822e-04,
|
|
Packit |
67cb25 |
0.63124427962699e-05,
|
|
Packit |
67cb25 |
0.5648193141174e-06,
|
|
Packit |
67cb25 |
0.512832480138e-07,
|
|
Packit |
67cb25 |
0.47196532914e-08,
|
|
Packit |
67cb25 |
0.4380744214e-09,
|
|
Packit |
67cb25 |
0.410268149e-10,
|
|
Packit |
67cb25 |
0.38623072e-11,
|
|
Packit |
67cb25 |
0.3661323e-12,
|
|
Packit |
67cb25 |
0.348023e-13,
|
|
Packit |
67cb25 |
0.33301e-14,
|
|
Packit |
67cb25 |
0.319e-15,
|
|
Packit |
67cb25 |
0.307e-16,
|
|
Packit |
67cb25 |
0.3e-17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series synchrotron2a_cs = {
|
|
Packit |
67cb25 |
synchrotron2a_data,
|
|
Packit |
67cb25 |
16,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
8
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_synchrotron_1_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*M_SQRT2 * GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
/* BJG: added first order correction term. The taylor series
|
|
Packit |
67cb25 |
is S1(x) = ((4pi)/(sqrt(3)gamma(1/3))) * (x/2)^(1/3)
|
|
Packit |
67cb25 |
* (1 - (gamma(1/3)/2)*(x/2)^2/3 + (3/4) * (x/2)^2 ....) */
|
|
Packit |
67cb25 |
double z = pow(x, 1.0/3.0);
|
|
Packit |
67cb25 |
double cf = 1 - 8.43812762813205e-01 * z * z;
|
|
Packit |
67cb25 |
result->val = 2.14952824153447863671 * z * cf;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double c0 = M_PI/M_SQRT3;
|
|
Packit |
67cb25 |
const double px = pow(x,1.0/3.0);
|
|
Packit |
67cb25 |
const double px11 = gsl_sf_pow_int(px,11);
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
gsl_sf_result result_c2;
|
|
Packit |
67cb25 |
cheb_eval_e(&synchrotron1_cs, t, &result_c1);
|
|
Packit |
67cb25 |
cheb_eval_e(&synchrotron2_cs, t, &result_c2);
|
|
Packit |
67cb25 |
result->val = px * result_c1.val - px11 * result_c2.val - c0 * x;
|
|
Packit |
67cb25 |
result->err = px * result_c1.err + px11 * result_c2.err + c0 * x * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -8.0*GSL_LOG_DBL_MIN/7.0) {
|
|
Packit |
67cb25 |
const double c0 = 0.2257913526447274323630976; /* log(sqrt(pi/2)) */
|
|
Packit |
67cb25 |
const double t = (12.0 - x) / (x + 4.0);
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
cheb_eval_e(&synchrotron1a_cs, t, &result_c1);
|
|
Packit |
67cb25 |
result->val = sqrt(x) * result_c1.val * exp(c0 - x);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * result->val * (fabs(c0-x)+1.0);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_synchrotron_2_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
/* BJG: added first order correction term. The taylor series
|
|
Packit |
67cb25 |
is S2(x) = ((2pi)/(sqrt(3)*gamma(1/3))) * (x/2)^(1/3)
|
|
Packit |
67cb25 |
* (1 - (gamma(1/3)/gamma(4/3))*(x/2)^(4/3) + (gamma(1/3)/gamma(4/3))*(x/2)^2...) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double z = pow(x, 1.0/3.0);
|
|
Packit |
67cb25 |
double cf = 1 - 1.17767156510235e+00 * z * x;
|
|
Packit |
67cb25 |
result->val = 1.07476412076723931836 * z * cf ;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double px = pow(x, 1.0/3.0);
|
|
Packit |
67cb25 |
const double px5 = gsl_sf_pow_int(px,5);
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result cheb1;
|
|
Packit |
67cb25 |
gsl_sf_result cheb2;
|
|
Packit |
67cb25 |
cheb_eval_e(&synchrotron21_cs, t, &cheb1);
|
|
Packit |
67cb25 |
cheb_eval_e(&synchrotron22_cs, t, &cheb2);
|
|
Packit |
67cb25 |
result->val = px * cheb1.val - px5 * cheb2.val;
|
|
Packit |
67cb25 |
result->err = px * cheb1.err + px5 * cheb2.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -8.0*GSL_LOG_DBL_MIN/7.0) {
|
|
Packit |
67cb25 |
const double c0 = 0.22579135264472743236; /* log(sqrt(pi/2)) */
|
|
Packit |
67cb25 |
const double t = (10.0 - x) / (x + 2.0);
|
|
Packit |
67cb25 |
gsl_sf_result cheb1;
|
|
Packit |
67cb25 |
cheb_eval_e(&synchrotron2a_cs, t, &cheb1);
|
|
Packit |
67cb25 |
result->val = sqrt(x) * exp(c0-x) * cheb1.val;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val * (fabs(c0-x)+1.0);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
UNDERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_synchrotron_1(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_synchrotron_1_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_synchrotron_2(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_synchrotron_2_e(x, &result));
|
|
Packit |
67cb25 |
}
|