bessel_olver.c

多个常用的特殊函数的例子

/* specfunc/bessel_olver.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 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
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_airy.h>

#include "error.h"

#include "bessel.h"
#include "bessel_olver.h"

#include "chebyshev.h"
#include "cheb_eval.c"

/* fit for f(x) = zofmzeta((x+1)/2), 0 <= mzeta <= 1 */
static double zofmzeta_a_data[20] = {
  2.9332563730829348990,
  0.4896518224847036624,
  0.0228637617355380860,
 -0.0001715731377284693,
 -0.0000105927538148751,
  1.0595602530419e-6,
 -4.68016051691e-8,
  5.8310020e-12,
  1.766537581e-10,
 -1.45034640e-11,
  4.357772e-13,
  4.60971e-14,
 -2.57571e-14,
  2.26468e-14,
 -2.22053e-14,
  2.08593e-14,
 -1.84454e-14,
  1.50150e-14,
 -1.06506e-14,
  5.5375e-15
};
static cheb_series zofmzeta_a_cs = {
  zofmzeta_a_data,
  19,
  -1,1,
  8
};


/* fit for f(x) = zofmzeta((9x+11)/2), 1 <= mzeta <= 10 */
static double zofmzeta_b_data[30] = {
  22.40725276466303489,
  10.39808258825165581,
  1.092050144486018425,
 -0.071111274777921604,
  0.008990125336059704,
 -0.001201950338088875,
  0.000106686807968315,
  0.000017406491576830,
 -0.000014946669657805,
  6.189984487752e-6,
 -2.049466715178e-6,
  5.87189458020e-7,
 -1.46077514157e-7,
  2.9803936132e-8,
 -3.817692108e-9,
 -4.66980416e-10,
  5.83860334e-10,
 -2.78825299e-10,
  1.01682688e-10,
 -3.1209928e-11,
  8.111122e-12,
 -1.663986e-12,
  1.81364e-13,
  5.3414e-14,
 -4.7234e-14,
  2.1689e-14,
 -7.815e-15,
  2.371e-15,
 -6.04e-16,
  1.20e-16
};
static cheb_series zofmzeta_b_cs = {
  zofmzeta_b_data,
  29,
  -1,1,
  15
};


/* fit for f(x) = zofmzeta(mz(x))/mz(x)^(3/2),
 * mz(x) = (2/(x+1))^(2/3) 10
 * 10 <= mzeta <= Inf
 */
static double zofmzeta_c_data[11] = {
  1.3824761227122911500,
  0.0244856101686774245,
 -0.0000842866496282540,
  1.4656076569771e-6,
 -3.14874099476e-8,
  7.561134833e-10,
 -1.94531643e-11,
  5.245878e-13,
 -1.46380e-14,
  4.192e-16,
 -1.23e-17
};
static cheb_series zofmzeta_c_cs = {
  zofmzeta_c_data,
  10,
  -1,1,
  6
};


/* Invert [Abramowitz+Stegun, 9.3.39].
 * Assumes minus_zeta >= 0.
 */
double
gsl_sf_bessel_Olver_zofmzeta(double minus_zeta)
{
  if(minus_zeta < 1.0) {
    const double x = 2.0*minus_zeta - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&zofmzeta_a_cs, x, &c);
    return c.val;
  }
  else if(minus_zeta < 10.0) {
    const double x = (2.0*minus_zeta - 11.0)/9.0;
    gsl_sf_result c;
    cheb_eval_e(&zofmzeta_b_cs, x, &c);
    return c.val;
  }
  else {
    const double TEN_32 = 31.62277660168379332; /* 10^(3/2) */
    const double p = pow(minus_zeta, 3.0/2.0);
    const double x = 2.0*TEN_32/p - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&zofmzeta_c_cs, x, &c);
    return c.val * p;
  }
}


/* Chebyshev fit for f(x) = z(x)^6 A_3(z(x)),  z(x) = 22/(10(x+1)) */
static double A3_gt1_data[31] = {
  -0.123783199829515294670493131190,
   0.104636462534700704670877382304,
  -0.067500816575851826744877535903,
   0.035563362418888483652711005520,
  -0.0160738524035979408472979609051,
   0.0064497878252851092073278056238,
  -0.00235408261133449663958121821593,
   0.00079545702851302155411892534965,
  -0.00025214920745855079895784825637,
   0.00007574004596069392921153301833,
  -0.00002172917966339623434407978263,
   5.9914810727868915476543145465e-06,
  -1.5958781571808992162953719817e-06,
   4.1232986512903717525448312012e-07,
  -1.0369725993417659101913919101e-07,
   2.5457982304266541145999235022e-08,
  -6.1161715053791743082427422443e-09,
   1.4409346199138658887871461320e-09,
  -3.3350445956255561668232014995e-10,
   7.5950686572918996453336138108e-11,
  -1.7042296334409430377389900278e-11,
   3.7723525020626230919721640081e-12,
  -8.2460237635733980528416501227e-13,
   1.7816961527997797696251868875e-13,
  -3.8084101506541792942694560802e-14,
   8.0593669930916099079755351563e-15,
  -1.6896565961641739017452636964e-15,
   3.5115651805888443184822853595e-16,
  -7.2384771938569255638904297651e-17,
   1.4806598977677176106283840244e-17,
  -3.0069285750787303634897997963e-18
};
static cheb_series A3_gt1_cs = {
  A3_gt1_data,
  30,
  -1,1,
  17
};

/* chebyshev expansion for f(x) = z(x)^8 A_4(z(x)), z(x) = 12/(5(x+1)) */
static double A4_gt1_data[30] = {
  1.15309329391198493586724229008,
 -1.01812701728669338904729927846,
  0.71964022270555684403652781941,
 -0.42359963977172689685150061355,
  0.215024488759339557817435404261,
 -0.096751915348145944032096342479,
  0.039413982058824310099856035361,
 -0.014775225692561697963781115014,
  0.005162114514159370516947823271,
 -0.00169783446445524322560925166335,
  0.00052995667873006847211519193478,
 -0.00015802027574996477115667974856,
  0.000045254366680989687988902825193,
 -0.000012503722965474638015488600967,
  3.3457656998119148699124716204e-06,
 -8.6981575241150758412492331833e-07,
  2.2030895484325645640823940625e-07,
 -5.4493369492600677068285936533e-08,
  1.3190457281724829107139385556e-08,
 -3.1301560183377379158951191769e-09,
  7.2937802527123344842593076131e-10,
 -1.6712080137945140407348940109e-10,
  3.7700053248213600430503521194e-11,
 -8.3824538848817227637828899571e-12,
  1.8388741910049766865274037194e-12,
 -3.9835919980753778560117573063e-13,
  8.5288827136546615604290389711e-14,
 -1.8060227869114416998653266836e-14,
  3.7849342199690728470461022877e-15,
 -7.8552867468122209577151823365e-16
};
static cheb_series A4_gt1_cs = {
  A4_gt1_data,
  17, /* 29, */
  -1, 1,
  17
};

/* Chebyshev fit for f(x) = z(x)^3 B_2(z(x)), z(x) = 12/(5(x+1)) */
static double B2_gt1_data[40] = {
  0.00118587147272683864479328868589,
  0.00034820459990648274622193981840,
 -0.00030411304425639768103075864567,
  0.00002812066284012343531484682886,
  0.00004493525295901613184489898748,
 -0.00003037629997093072196779489677,
  0.00001125979647123875721949743970,
 -2.4832533969517775991951008218e-06,
 -9.9003813640537799587086928278e-08,
  4.9259859656183110299492296029e-07,
 -3.7644120964426705960749504975e-07,
  2.2887828521334625189639122509e-07,
 -1.3202687370822203731489855050e-07,
  7.7019669092537400811434860763e-08,
 -4.6589706973010511603890144294e-08,
  2.9396476233013923711978522963e-08,
 -1.9293230611988282919101954538e-08,
  1.3099107013728717842406906896e-08,
 -9.1509111940885962831104149355e-09,
  6.5483472971925614347299375295e-09,
 -4.7831253582139967461241674569e-09,
  3.5562625457426178152760148639e-09,
 -2.6853389444008414186916562103e-09,
  2.0554738667134200145781857289e-09,
 -1.5923172019517426277886522758e-09,
  1.2465923213464381457319481498e-09,
 -9.8494846881180588507969988989e-10,
  7.8438674499372126663957464312e-10,
 -6.2877567918342950225937136855e-10,
  5.0662318868755257959686944117e-10,
 -4.0962270881243451160378710952e-10,
  3.3168684677374908553161911299e-10,
 -2.6829406619847450633596163305e-10,
  2.1603988122184568375561077873e-10,
 -1.7232373309560278402012124481e-10,
  1.3512709089611470626617830434e-10,
 -1.0285354732538663013167579792e-10,
  7.4211345443901713467637018423e-11,
 -4.8124980266864320351456993068e-11,
  2.3666534694476306077416831958e-11
};
static cheb_series B2_gt1_cs = {
  B2_gt1_data,
  39,
  -1, 1,
  30
};


/* Chebyshev fit for f(x) = z(x)^6 B_3(z(x)), z(x) = 12/(5(x+1)) */
static double B3_gt1_data[30] = {
 -0.0102445379362695740863663926486,
  0.0036618484329295342954730801917,
  0.0026154252498599303282569321117,
 -0.0036187389410353156728771706336,
  0.0021878564157692275944613452462,
 -0.0008219952303590803584426516821,
  0.0001281773889155631494321316520,
  0.0001000944653368032985720548637,
 -0.0001288293344663774273453147788,
  0.00010136264202696513867821487205,
 -0.00007000275849659556221916572733,
  0.00004694886396757430431607955146,
 -0.00003190003869717837686356945696,
  0.00002231453668447775219665947479,
 -0.00001611102197712439539300336438,
  0.00001196634424990735214466633513,
 -9.0986920398931223804111374679e-06,
  7.0492613694235423068926562567e-06,
 -5.5425216624642184684300615394e-06,
  4.4071884714230296614449244106e-06,
 -3.5328595506791663127928952625e-06,
  2.84594975572077091520522824686e-06,
 -2.29592697828824392391071619788e-06,
  1.84714740375289956396370322228e-06,
 -1.47383331248116454652025598620e-06,
  1.15687781098593231076084710267e-06,
 -8.8174688524627071175315084910e-07,
  6.3705856964426840441434605593e-07,
 -4.1358791499961929237755474814e-07,
  2.0354151158738819867477996807e-07
};
static cheb_series B3_gt1_cs = {
  B3_gt1_data,
  29,