htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/synchrotron.c

   1 /* specfunc/synchrotron.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_exp.h>
  26 #include <gsl/gsl_sf_pow_int.h>
  27 #include <gsl/gsl_sf_synchrotron.h>
  28
  29 #include "error.h"
  30
  31 #include "chebyshev.h"
  32 #include "cheb_eval.c"
  33
  34 static double synchrotron1_data[13] = {
  35   30.364682982501076273,
  36   17.079395277408394574,
  37    4.560132133545072889,
  38    0.549281246730419979,
  39    0.372976075069301172e-01,
  40    0.161362430201041242e-02,
  41    0.481916772120371e-04,
  42    0.10512425288938e-05,
  43    0.174638504670e-07,
  44    0.22815486544e-09,
  45    0.240443082e-11,
  46    0.2086588e-13,
  47    0.15167e-15
  48 };
  49 static cheb_series synchrotron1_cs = {
  50   synchrotron1_data,
  51   12,
  52   -1.0, 1.0,
  53   9
  54 };
  55
  56 static double synchrotron2_data[12] = {
  57   0.4490721623532660844,
  58   0.898353677994187218e-01,
  59   0.81044573772151290e-02,
  60   0.4261716991089162e-03,
  61   0.147609631270746e-04,
  62   0.3628633615300e-06,
  63   0.66634807498e-08,
  64   0.949077166e-10,
  65   0.1079125e-11,
  66   0.10022e-13,
  67   0.77e-16,
  68   0.5e-18
  69 };
  70 static cheb_series synchrotron2_cs = {
  71   synchrotron2_data,
  72   11,
  73   -1.0, 1.0,
  74   7
  75 };
  76
  77 static double synchrotron1a_data[23] = {
  78   2.1329305161355000985,
  79   0.741352864954200240e-01,
  80   0.86968099909964198e-02,
  81   0.11703826248775692e-02,
  82   0.1645105798619192e-03,
  83   0.240201021420640e-04,
  84   0.35827756389389e-05,
  85   0.5447747626984e-06,
  86   0.838802856196e-07,
  87   0.13069882684e-07,
  88   0.2053099071e-08,
  89   0.325187537e-09,
  90   0.517914041e-10,
  91   0.83002988e-11,
  92   0.13352728e-11,
  93   0.2159150e-12,
  94   0.349967e-13,
  95   0.56994e-14,
  96   0.9291e-15,
  97   0.152e-15,
  98   0.249e-16,
  99   0.41e-17,
 100   0.7e-18
 101 };
 102 static cheb_series synchrotron1a_cs = {
 103   synchrotron1a_data,
 104   22,
 105   -1.0, 1.0,
 106   11
 107 };
 108
 109 static double synchrotron21_data[13] = {
 110   38.617839923843085480,
 111   23.037715594963734597,
 112   5.3802499868335705968,
 113   0.6156793806995710776,
 114   0.406688004668895584e-01,
 115   0.17296274552648414e-02,
 116   0.51061258836577e-04,
 117   0.110459595022e-05,
 118   0.18235530206e-07,
 119   0.2370769803e-09,
 120   0.24887296e-11,
 121   0.21529e-13,
 122   0.156e-15
 123 };
 124 static cheb_series synchrotron21_cs = {
 125   synchrotron21_data,
 126   12,
 127   -1.0, 1.0,
 128   9
 129 };
 130
 131 static double synchrotron22_data[13] = {
 132    7.9063148270660804288,
 133    3.1353463612853425684,
 134    0.4854879477453714538,
 135    0.394816675827237234e-01,
 136    0.19661622334808802e-02,
 137    0.659078932293042e-04,
 138    0.15857561349856e-05,
 139    0.286865301123e-07,
 140    0.4041202360e-09,
 141    0.45568444e-11,
 142    0.420459e-13,
 143    0.3232e-15,
 144    0.21e-17
 145 };
 146 static cheb_series synchrotron22_cs = {
 147   synchrotron22_data,
 148   12,
 149   -1.0, 1.0,
 150   8
 151 };
 152
 153 static double synchrotron2a_data[17] = {
 154   2.020337094170713600,
 155   0.10956237121807404e-01,
 156   0.8542384730114676e-03,
 157   0.723430242132822e-04,
 158   0.63124427962699e-05,
 159   0.5648193141174e-06,
 160   0.512832480138e-07,
 161   0.47196532914e-08,
 162   0.4380744214e-09,
 163   0.410268149e-10,
 164   0.38623072e-11,
 165   0.3661323e-12,
 166   0.348023e-13,
 167   0.33301e-14,
 168   0.319e-15,
 169   0.307e-16,
 170   0.3e-17
 171 };
 172 static cheb_series synchrotron2a_cs = {
 173   synchrotron2a_data,
 174   16,
 175   -1.0, 1.0,
 176   8
 177 };
 178
 179
 180 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 181
 182 int gsl_sf_synchrotron_1_e(const double x, gsl_sf_result * result)
 183 {
 184   /* CHECK_POINTER(result) */
 185
 186   if(x < 0.0) {
 187     DOMAIN_ERROR(result);
 188   }
 189   else if(x < 2.0*M_SQRT2 * GSL_SQRT_DBL_EPSILON) {
 190     /* BJG: added first order correction term.  The taylor series
 191        is  S1(x) = ((4pi)/(sqrt(3)gamma(1/3))) * (x/2)^(1/3)
 192        * (1 - (gamma(1/3)/2)*(x/2)^2/3 + (3/4) * (x/2)^2 ....) */
 193     double z = pow(x, 1.0/3.0);
 194     double cf = 1 - 8.43812762813205e-01 * z * z;
 195     result->val = 2.14952824153447863671 * z * cf;
 196     result->err = GSL_DBL_EPSILON * result->val;
 197     return GSL_SUCCESS;
 198   }
 199   else if(x <= 4.0) {
 200     const double c0   = M_PI/M_SQRT3;
 201     const double px   = pow(x,1.0/3.0);
 202     const double px11 = gsl_sf_pow_int(px,11);
 203     const double t = x*x/8.0 - 1.0;
 204     gsl_sf_result result_c1;
 205     gsl_sf_result result_c2;
 206     cheb_eval_e(&synchrotron1_cs, t, &result_c1);
 207     cheb_eval_e(&synchrotron2_cs, t, &result_c2);
 208     result->val  = px * result_c1.val - px11 * result_c2.val - c0 * x;
 209     result->err  = px * result_c1.err + px11 * result_c2.err + c0 * x * GSL_DBL_EPSILON;
 210     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 211     return GSL_SUCCESS;
 212   }
 213   else if(x < -8.0*GSL_LOG_DBL_MIN/7.0) {
 214     const double c0 = 0.2257913526447274323630976; /* log(sqrt(pi/2)) */
 215     const double t = (12.0 - x) / (x + 4.0);
 216     gsl_sf_result result_c1;
 217     cheb_eval_e(&synchrotron1a_cs, t, &result_c1);
 218     result->val = sqrt(x) * result_c1.val * exp(c0 - x);
 219     result->err = 2.0 * GSL_DBL_EPSILON * result->val * (fabs(c0-x)+1.0);
 220     return GSL_SUCCESS;
 221   }
 222   else {
 223     UNDERFLOW_ERROR(result);
 224   }
 225 }
 226
 227
 228 int gsl_sf_synchrotron_2_e(const double x, gsl_sf_result * result)
 229 {
 230   /* CHECK_POINTER(result) */
 231
 232   if(x < 0.0) {
 233     DOMAIN_ERROR(result);
 234   }
 235   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
 236     /* BJG: added first order correction term.  The taylor series
 237        is  S2(x) = ((2pi)/(sqrt(3)*gamma(1/3))) * (x/2)^(1/3)
 238        * (1 - (gamma(1/3)/gamma(4/3))*(x/2)^(4/3) + (gamma(1/3)/gamma(4/3))*(x/2)^2...) */
 239
 240     double z = pow(x, 1.0/3.0);
 241     double cf = 1 - 1.17767156510235e+00 * z * x;
 242     result->val = 1.07476412076723931836 * z * cf ;
 243     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
 244     return GSL_SUCCESS;
 245   }
 246   else if(x <= 4.0) {
 247     const double px  = pow(x, 1.0/3.0);
 248     const double px5 = gsl_sf_pow_int(px,5);
 249     const double t = x*x/8.0 - 1.0;
 250     gsl_sf_result cheb1;
 251     gsl_sf_result cheb2;
 252     cheb_eval_e(&synchrotron21_cs, t, &cheb1);
 253     cheb_eval_e(&synchrotron22_cs, t, &cheb2);
 254     result->val  = px * cheb1.val - px5 * cheb2.val;
 255     result->err  = px * cheb1.err + px5 * cheb2.err;
 256     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 257     return GSL_SUCCESS;
 258   }
 259   else if(x < -8.0*GSL_LOG_DBL_MIN/7.0) {
 260     const double c0 = 0.22579135264472743236;   /* log(sqrt(pi/2)) */
 261     const double t  = (10.0 - x) / (x + 2.0);
 262     gsl_sf_result cheb1;
 263     cheb_eval_e(&synchrotron2a_cs, t, &cheb1);
 264     result->val = sqrt(x) * exp(c0-x) * cheb1.val;
 265     result->err = GSL_DBL_EPSILON * result->val * (fabs(c0-x)+1.0);
 266     return GSL_SUCCESS;
 267   }
 268   else {
 269     UNDERFLOW_ERROR(result);
 270   }
 271 }
 272
 273 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 274
 275 #include "eval.h"
 276
 277 double gsl_sf_synchrotron_1(const double x)
 278 {
 279   EVAL_RESULT(gsl_sf_synchrotron_1_e(x, &result));
 280 }
 281
 282 double gsl_sf_synchrotron_2(const double x)
 283 {
 284   EVAL_RESULT(gsl_sf_synchrotron_2_e(x, &result));
 285 }