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

   1 /* specfunc/coulomb_bound.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_gamma.h>
  27 #include <gsl/gsl_sf_pow_int.h>
  28 #include <gsl/gsl_sf_laguerre.h>
  29 #include <gsl/gsl_sf_coulomb.h>
  30
  31 #include "error.h"
  32 #include "check.h"
  33
  34 /* normalization for hydrogenic wave functions */
  35 static
  36 int
  37 R_norm(const int n, const int l, const double Z, gsl_sf_result * result)
  38 {
  39   double A   = 2.0*Z/n;
  40   double pre = sqrt(A*A*A /(2.0*n));
  41   gsl_sf_result ln_a, ln_b;
  42   gsl_sf_result ex;
  43   int stat_a = gsl_sf_lnfact_e(n+l, &ln_a);
  44   int stat_b = gsl_sf_lnfact_e(n-l-1, &ln_b);
  45   double diff_val = 0.5*(ln_b.val - ln_a.val);
  46   double diff_err = 0.5*(ln_b.err + ln_a.err) + GSL_DBL_EPSILON * fabs(diff_val);
  47   int stat_e = gsl_sf_exp_err_e(diff_val, diff_err, &ex);
  48   result->val  = pre * ex.val;
  49   result->err  = pre * ex.err;
  50   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  51   return GSL_ERROR_SELECT_3(stat_e, stat_a, stat_b);
  52 }
  53
  54
  55 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  56
  57 int
  58 gsl_sf_hydrogenicR_1_e(const double Z, const double r, gsl_sf_result * result)
  59 {
  60   if(Z > 0.0 && r >= 0.0) {
  61     double A = 2.0*Z;
  62     double norm = A*sqrt(Z);
  63     double ea = exp(-Z*r);
  64     result->val = norm*ea;
  65     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) * fabs(Z*r);
  66     CHECK_UNDERFLOW(result);
  67     return GSL_SUCCESS;
  68   }
  69   else {
  70     DOMAIN_ERROR(result);
  71   }
  72 }
  73
  74
  75 int
  76 gsl_sf_hydrogenicR_e(const int n, const int l,
  77                         const double Z, const double r,
  78                         gsl_sf_result * result)
  79 {
  80   if(n < 1 || l > n-1 || Z <= 0.0 || r < 0.0) {
  81     DOMAIN_ERROR(result);
  82   }
  83   else {
  84     double A = 2.0*Z/n;
  85     gsl_sf_result norm;
  86     int stat_norm = R_norm(n, l, Z, &norm);
  87     double rho = A*r;
  88     double ea = exp(-0.5*rho);
  89     double pp = gsl_sf_pow_int(rho, l);
  90     gsl_sf_result lag;
  91     int stat_lag = gsl_sf_laguerre_n_e(n-l-1, 2*l+1, rho, &lag);
  92     double W_val = norm.val * ea * pp;
  93     double W_err = norm.err * ea * pp;
  94     W_err += norm.val * ((0.5*rho + 1.0) * GSL_DBL_EPSILON) * ea * pp;
  95     W_err += norm.val * ea * ((l+1.0) * GSL_DBL_EPSILON) * pp;
  96     result->val  = W_val * lag.val;
  97     result->err  = W_val * lag.err + W_err * fabs(lag.val);
  98     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  99     if ((l == 0 || (r > 0 && l > 0)) && lag.val != 0.0
 100         && stat_lag == GSL_SUCCESS && stat_norm == GSL_SUCCESS) {
 101       CHECK_UNDERFLOW(result);
 102     };
 103     return GSL_ERROR_SELECT_2(stat_lag, stat_norm);
 104   }
 105 }
 106
 107 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 108
 109 #include "eval.h"
 110
 111 double gsl_sf_hydrogenicR_1(const double Z, const double r)
 112 {
 113   EVAL_RESULT(gsl_sf_hydrogenicR_1_e(Z, r, &result));
 114 }
 115
 116
 117 double gsl_sf_hydrogenicR(const int n, const int l, const double Z, const double r)
 118 {
 119   EVAL_RESULT(gsl_sf_hydrogenicR_e(n, l, Z, r, &result));
 120 }