home *** CD-ROM | disk | FTP | other *** search
/ System Booster / System Booster.iso / Archives / ForCLI / Eval113.lha / Eval / source / emath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-04-14  |  19.0 KB  |  879 lines
  1. /*
  2. **
  3. ** EMATH.C     A module to be included with EVAL.C that supplies
  4. **             math functions not included in ANSI C.
  5. **
  6. ** To add/change the functions in Eval, see funcs.c.
  7. **
  8. ** NOTICE about this part of the Eval source code:
  9. ** ------------------------------------------------------------------------
  10. ** If you add source code to this module that is copied from "Numerical
  11. ** Recipes in C" or some other book of C math functions, make sure you
  12. ** clearly read the license agreement from that book.  In particular,
  13. ** "Numerical Recipes in C" forbids you to distribute their source code.
  14. ** Because my license agreement requires free distribution of ALL source
  15. ** code, you may not add code from a book like "Numerical Recipes in C" and
  16. ** then redistribute any of the source or executables.  The upshot of all
  17. ** this is that you must write your own versions of whatever math functions
  18. ** you add to this code if you plan to redistribute it.
  19. **
  20. **
  21. ** Eval is a floating point expression evaluator.
  22. ** This file last updated in version 1.12
  23. ** For the version number, see eval.h
  24. ** Copyright (C) 1993  Will Menninger
  25. **
  26. ** This program is free software; you can redistribute it and/or modify it
  27. ** under the terms of the GNU General Public License as published by the
  28. ** Free Software Foundation; either version 2 of the License, or any
  29. ** later version.
  30. **
  31. ** This program is distributed in the hope that it will be useful, but
  32. ** WITHOUT ANY WARRANTY; without even the implied warranty of
  33. ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  34. ** General Public License for more details.
  35. **
  36. ** You should have received a copy of the GNU General Public License along
  37. ** with this program; if not, write to the Free Software Foundation, Inc.,
  38. ** 675 Mass Ave, Cambridge, MA 02139, USA.
  39. **
  40. ** The author until 9/93 can be contacted at:
  41. ** e-mail:     willus@ilm.pfc.mit.edu
  42. ** U.S. mail:  Will Menninger, 45 River St., #2, Boston, MA 02108-1124
  43. **
  44. */
  45.  
  46. #include "eval.h"
  47.  
  48.  
  49. #define ACCURACY    ((double)1.0e-10)
  50. #define SMALL_PI    3.14
  51. #define MAXRESULT   (DBL_MAX/1.e3)
  52. #define ACC         40.0
  53. #define EXTR        1.0e10
  54. #define MAX         1e14
  55. #define FACTOR      2
  56. #define FIRSTGUESS  .11
  57. #define TINYNUM     1e-6
  58. #define NITER       30
  59.  
  60. static double x0;
  61.  
  62. static double bessi0(double x);
  63. static double bessi1(double x);
  64. static double bessj0(double x);
  65. static double bessj1(double x);
  66. static double bessk0(double x);
  67. static double bessk1(double x);
  68. static double hypcos(double x);
  69. static double hypsin(double x);
  70. static double hyptan(double x);
  71.  
  72. /*
  73. **
  74. **  acosh      Calculates arc-hyperbolic-cosine of a value.
  75. **             Arguments must be >=1.  If not, -1 is returned.
  76. **             Positive values are always returned if the
  77. **             argument is within the correct domain.
  78. **             Abs(Maximum value returned) = 1e14
  79. **
  80. */
  81. double acosh(double x)
  82.  
  83.    {
  84.    int     i;
  85.    double  b1,b2;
  86.  
  87.    if (x<1)
  88.        return(-1.);
  89.    if (x==1.)
  90.        return(0.);
  91.    x0=x;
  92.    if (cosh(.11)>x)
  93.        {
  94.        b1=0.;
  95.        b2=.11;
  96.        }
  97.    else
  98.        {
  99.        for (b2=.1;cosh(b2)<x && b2<MAX;b2*=FACTOR);
  100.        if (b2>MAX)
  101.            return(MAX);
  102.        b1=b2/FACTOR;
  103.        }
  104.    return(root(hypcos,b1,b2,ACCURACY,&i));
  105.    }
  106.  
  107.  
  108. static double hypcos(double x)
  109.  
  110.    {
  111.    return(cosh(x)-x0);
  112.    }
  113.  
  114.  
  115. /*
  116. **
  117. **  asinh      Calculates arc-hyperbolic-sine of a value.
  118. **
  119. **             Abs(Maximum value returned) = 1e14
  120. **
  121. */
  122. double asinh(double x)
  123.  
  124.    {
  125.    int     i,sign;
  126.    double  b1,b2,acc;
  127.  
  128.    if (x==0.)
  129.        return(0.);
  130.    if (x>0.)
  131.        sign=1;
  132.    else
  133.        {
  134.        sign=-1;
  135.        x=-x;
  136.        }
  137.    x0=x;
  138.    acc=ACCURACY;
  139.    if (sinh(FIRSTGUESS)>x)
  140.        {
  141.        b1=0.;
  142.        b2=FIRSTGUESS;
  143.        }
  144.    else
  145.        {
  146.        for (b2=.1;sinh(b2)<x && b2<MAX;b2*=FACTOR,acc*=FACTOR);
  147.        if (b2>MAX)
  148.            return(MAX*sign);
  149.        b1=b2/FACTOR;
  150.        }
  151.    return(root(hypsin,b1,b2,acc,&i)*sign);
  152.    }
  153.  
  154.  
  155. static double hypsin(double x)
  156.  
  157.    {
  158.    return(sinh(x)-x0);
  159.    }
  160.  
  161.  
  162. /*
  163. **
  164. **  atanh      Calculates arc-hyperbolic-tangent of a value.
  165. **             Value must be between but not including -1 and 1
  166. **             or else 0 is returned.
  167. **
  168. **             Abs(Maximum value returned) = 1e14
  169. **
  170. */
  171. double atanh(double x)
  172.  
  173.    {
  174.    int     i,sign;
  175.    double  b1,b2;
  176.  
  177.    if (x<=-1 || x>=1)
  178.        return(0.);
  179.    if (x==0.)
  180.        return(0.);
  181.    if (x>0.)
  182.        sign=1;
  183.    else
  184.        {
  185.        sign=-1;
  186.        x=-x;
  187.        }
  188.    x0=x;
  189.    if (tanh(.11)>x)
  190.        {
  191.        b1=0.;
  192.        b2=.11;
  193.        }
  194.    else
  195.        {
  196.        for (b2=.1;tanh(b2)<x && b2<MAX;b2*=FACTOR);
  197.        if (b2>MAX)
  198.            return(MAX*sign);
  199.        b1=b2/FACTOR;
  200.        }
  201.    return(root(hyptan,b1,b2,ACCURACY,&i)*sign);
  202.    }
  203.  
  204.  
  205. static double hyptan(double x)
  206.  
  207.    {
  208.    return(tanh(x)-x0);
  209.    }
  210.  
  211.  
  212. /*
  213. **
  214. ** bessi    Calculates the value of Im(x) using a downward recurrence
  215. **          (Miller's) formula.  Calls bessi0 to normalize.
  216. **
  217. ** The idea for Miller's algorithm, and where to start the series,
  218. ** was taken from Numerical Recipes in C, but this function is NOT
  219. ** a copy of their function, and, I believe, does not infringe on
  220. ** their copyright.
  221. **
  222. ** Properties:  Im(x) is a solution to the modified (hyperbolic)
  223. **              Bessel equation.
  224. **              Im(x) is monotonically increasing and diverges
  225. **              at infinity.
  226. **
  227. **              I(-m)(x) = Im(x)
  228. **              Im(-x)   = (-1)^m*Im(x)
  229. **              I(m+1)(x) = -m*2/x*Im(x) + I(m-1)(x)
  230. **              I(m-1)(x) = m*2/x*Im(x) + I(m+1)(x)
  231. **
  232. */
  233. double bessi(int m,double x)
  234.  
  235.     {
  236.     int     i,negx;
  237.     double  c,ans,nextterm,seed0,seed1;
  238.  
  239.     if (m<0)
  240.         m=-m;
  241.     if (x==0.)
  242.         return(0.);
  243.     if (m==0)
  244.         return(bessi0(x));
  245.     if (m==1)
  246.         return(bessi1(x));
  247.     negx=(x<0);
  248.     x=fabs(x);
  249.     if (x==0.)
  250.         return(0.);
  251.     c=2./x;
  252.     seed1=1.0;  /* Initial seeds don't matter.  Series converges. */
  253.     seed0=0.0;
  254.     nextterm=0.;
  255.     ans=0.;
  256.     for (i=(int)(m+sqrt(64.*m))&(~1);i>0;i--)
  257.         {
  258.         nextterm=i*c*seed1+seed0;
  259.         if (fabs(nextterm)>EXTR)
  260.             {
  261.             nextterm/=EXTR;
  262.             ans/=EXTR;
  263.             seed1/=EXTR;
  264.             }
  265.         if (i==m)
  266.             ans=seed1;
  267.         seed0=seed1;
  268.         seed1=nextterm;
  269.         }
  270.     ans*=bessi0(x)/nextterm;
  271.     return(negx && (m&1) ? -ans : ans);
  272.     }
  273.  
  274.  
  275. /*
  276. ** bessi0   Calculate I0(x) from Abromowitz and Stegun, p. 378.
  277. **          Accuracy is better than 2e-7.
  278. */
  279. static double bessi0(double x)
  280.  
  281.     {
  282.     double  t,t2,t3,t4,t8,ans;
  283.  
  284.     x=fabs(x);
  285.     if (x<3.75)
  286.         {
  287.         /* Eq 9.8.1 */
  288.         t=(x/3.75);
  289.         t2=t*t;
  290.         t4=t2*t2;
  291.         t8=t4*t4;
  292.         ans=1+3.5156229*t2+3.0899424*t4+1.2067492*t4*t2+.2659732*t8
  293.              +.0360768*t8*t2+.0045813*t4*t8;
  294.         }
  295.     else
  296.         {
  297.         /* Eq 9.8.2 */
  298.         t=3.75/x;
  299.         t2=t*t;
  300.         t3=t2*t;
  301.         t4=t2*t2;
  302.         ans=.39894228+.01328592*t+.00225319*t2-.00157565*t3
  303.             +.00916281*t4-.02057706*t*t4+.02635537*t2*t4
  304.             -.01647633*t3*t4+.00392377*t4*t4;
  305.         ans=exp(x)*ans/sqrt(x);
  306.         }
  307.     return(ans);
  308.     }
  309.  
  310.  
  311. /*
  312. ** bessi1   Calculate I1(x) from Abromowitz and Stegun, p. 378.
  313. **          Accuracy is better than 3e-7.
  314. */
  315. static double bessi1(double x)
  316.  
  317.     {
  318.     int     negx;
  319.     double  t,t2,t3,t4,t8,ans;
  320.  
  321.     negx=(x<0.);
  322.     x=fabs(x);
  323.     if (x<3.75)
  324.         {
  325.         /* Eq 9.8.3 */
  326.         t=(x/3.75);
  327.         t2=t*t;
  328.         t4=t2*t2;
  329.         t8=t4*t4;
  330.         ans=.5+.87890594*t2+.51498869*t4+.15084934*t4*t2+.02658733*t8
  331.               +.00301532*t8*t2+.00032411*t8*t4;
  332.         ans*=x;
  333.         }
  334.     else
  335.         {
  336.         /* Eq 9.8.4 */
  337.         t=3.75/x;
  338.         t2=t*t;
  339.         t3=t2*t;
  340.         t4=t2*t2;
  341.         ans=.39894228-.03988024*t-.00362018*t2+.00163801*t3
  342.             -.01031555*t4+.02282967*t*t4-.02895312*t4*t2
  343.             +.01787654*t4*t3-.00420059*t4*t4;
  344.         ans=exp(x)*ans/sqrt(x);
  345.         }
  346.     return(negx ? -ans : ans);
  347.     }
  348.  
  349.  
  350. /*
  351. **
  352. ** bessj    Calculates the value of Jm(x) using the bessj0 and
  353. **          bessj0 functions and the Bessel function recurrence relation.
  354. **
  355. ** The idea for Miller's algorithm, when to use it, and where to start the
  356. ** series, was taken from Numerical Recipes in C, but this function is NOT
  357. ** a copy of their function, and, I believe, does not infringe on their
  358. ** copyright.
  359. **
  360. ** Properties:  Jm(x) is a solution to the Bessel function.
  361. **              It is oscillatory in nature and does not diverge at
  362. **              any point.
  363. **
  364. **              J(-m)(x) = (-1)^m*Jm(x)
  365. **              Jm(-x)   = (-1)^m*Jm(x)
  366. **              J(m+1)(x) = m*2/x*Jm(x) - J(m-1)(x)
  367. **              J(m-1)(x) = m*2/x*Jm(x) - J(m+1)(x)
  368. **
  369. **
  370. */
  371. double bessj(int m,double x)
  372.  
  373.     {
  374.     int     i,negm,negx;
  375.     double  c,j0,j1,ans,nextterm,seed1,seed0;
  376.  
  377.     negm=(m<0);
  378.     if (m<0)
  379.         m=-m;
  380.     if (!(m&1))
  381.         negm=0;
  382.     if (m==0)
  383.         return(bessj0(x));
  384.     if (m==1)
  385.         return(negm ? -bessj1(x) : bessj1(x));
  386.     negx=(x<0);
  387.     x=fabs(x);
  388.     if (x==0.)
  389.         return(0.);
  390.     c=2./x;
  391.     if (x>m)
  392.         {
  393.         /* Use forward recurrence relation */
  394.         j0=bessj0(x);
  395.         j1=bessj1(x);
  396.         for (i=1;i<m;i++)
  397.             {
  398.             ans=i*c*j1-j0;
  399.             j0=j1;
  400.             j1=ans;
  401.             }
  402.         }
  403.     else
  404.         {
  405.         /* Use Miller's algorithm of downward recurrence for x<m */
  406.         seed1=1.; /* Start with any seed value--series will converge */
  407.         seed0=0.;
  408.         ans=0.;
  409.         for (i=(int)(m+sqrt(64.*m))&(~1);i>0;i--)
  410.             {
  411.             nextterm=i*c*seed1-seed0;
  412.             if (fabs(nextterm)>EXTR)
  413.                 {
  414.                 nextterm/=EXTR;
  415.                 ans/=EXTR;
  416.                 seed1/=EXTR;
  417.                 }
  418.             if (i==m)
  419.                 ans=seed1;
  420.             seed0=seed1;
  421.             seed1=nextterm;
  422.             }
  423.         /* Normalize with bessj0 */
  424.         ans*=(bessj0(x)/nextterm);
  425.         }
  426.     if (negx && (m&1))
  427.         ans=-ans;
  428.     return(negm ? -ans : ans);
  429.     }
  430.  
  431.  
  432. /*
  433. ** bessj0   Calculate J0(x) from Abromowitz and Stegun, p. 369.
  434. **          Accuracy is better than 5e-8.
  435. */
  436. static double bessj0(double x)
  437.  
  438.     {
  439.     double  ans;
  440.     double  x2,x4,x8;
  441.     double  f0,theta0,thox,thox2,thox4;
  442.  
  443.     x=fabs(x);
  444.     if (x<3)
  445.         {
  446.         /* Eq 9.4.1 */
  447.  
  448.         x2=(x/3)*(x/3);
  449.         x4=x2*x2;
  450.         x8=x4*x4;
  451.         ans=1-2.2499997*x2+1.2656208*x4-.3163866*x2*x4
  452.              +.04444479*x8-.0039444*x8*x2+.0002100*x8*x4;
  453.         }
  454.     else
  455.         {
  456.         /* Eq 9.4.3 */
  457.  
  458.         thox=3./x;
  459.         thox2=thox*thox;
  460.         thox4=thox2*thox2;
  461.         f0=.79788456-.00000077*thox-.00552740*thox2-.00009512*thox2*thox
  462.            +.00137237*thox4-.00072805*thox4*thox+.00014476*thox4*thox2;
  463.         theta0=x-.78539816-.04166397*thox-.00003954*thox2+.00262573*thox2*thox
  464.                -.00054125*thox4-.00029333*thox4*thox+.00013558*thox2*thox4;
  465.         ans=f0*cos(theta0)/sqrt(x);
  466.         }
  467.     return(ans);
  468.     }
  469.  
  470.  
  471. /*
  472. ** bessj1   Calculate J1(x) from Abromowitz and Stegun, p. 370.
  473. **          Accuracy is better than 5e-8.
  474. */
  475. static double bessj1(double x)
  476.  
  477.     {
  478.     int     negx;
  479.     double  ans;
  480.     double  x2,x4,x8;
  481.     double  f1,theta1,thox,thox2,thox4;
  482.  
  483.     negx=(x<0);
  484.     x=fabs(x);
  485.     if (x<3)
  486.         {
  487.         /* Eq 9.4.4 */
  488.  
  489.         x2=(x/3)*(x/3);
  490.         x4=x2*x2;
  491.         x8=x4*x4;
  492.         ans=0.5-.56249985*x2+.21093573*x4-.03954289*x4*x2+.00443319*x8
  493.             -.00031761*x8*x2+.00001109*x4*x8;
  494.         ans*=x;
  495.         }
  496.     else
  497.         {
  498.         /* Eq 9.4.6 */
  499.  
  500.         thox=3./x;
  501.         thox2=thox*thox;
  502.         thox4=thox2*thox2;
  503.         f1=.79788456+.00000156*thox+.01659667*thox2+.00017105*thox*thox2
  504.            -.00249511*thox4+.00113653*thox4*thox-.00020033*thox2*thox4;
  505.         theta1=x-2.35619449+.12499612*thox+.00005650*thox2-.00637879*thox*thox2
  506.                +.00074348*thox4+.00079824*thox4*thox-.00029166*thox4*thox2;
  507.         ans=f1*cos(theta1)/sqrt(x);
  508.         }
  509.     return(negx ? -ans : ans);
  510.     }
  511.  
  512.  
  513. /*
  514. **
  515. ** bessk   Calculates the value of Km(x) using Bessel function
  516. **         recurrence relations.
  517. **
  518. ** Properties:  Km(x) is a solution to the modified
  519. **              (hyperbolic) Bessel equation.  It diverges
  520. **              at zero and is monotonically decreasing.
  521. **
  522. **              Km(x) can also be written as:
  523. **
  524. **                           pi   I(-v)(x) - Iv(x)
  525. **              lim (v->m)  ---- ------------------
  526. **                            2       sin(v*x)
  527. **
  528. **              (only when m is an integer)
  529. **
  530. **
  531. **              K(-m)(x) = Km(x)
  532. **              Km(-x)   = (-1)^m*Km(x)
  533. **              K(m+1)(x) = m*2/x*Km(x) + K(m-1)(x)
  534. **
  535. */
  536. double bessk(int m,double x)
  537.  
  538.     {
  539.     int     i,negx;
  540.     double  c,k0,k1,ans;
  541.  
  542.     if (m<0)
  543.         m=-m;
  544.     if (m==0)
  545.         return(bessk0(x));
  546.     if (m==1)
  547.         return(bessk1(x));
  548.     negx=(x<0.);
  549.     x=fabs(x);
  550.     if (x==0.)
  551.         return(DBL_MAX);
  552.     c=2./x;
  553.     k0=bessk0(x);
  554.     k1=bessk1(x);
  555.     for (i=1;i<m;i++)
  556.         {
  557.         ans=c*i*k1+k0;
  558.         k0=k1;
  559.         k1=ans;
  560.         }
  561.     return(negx && (m&1) ? -ans : ans);
  562.     }
  563.  
  564.  
  565. /*
  566. ** bessk0   Calculate K0(x) from Abromowitz and Stegun, p. 379.
  567. **          Accuracy is better than 2e-7.
  568. */
  569. static double bessk0(double x)
  570.  
  571.     {
  572.     double  ans,t,t2,t4,t8;
  573.  
  574.     x=fabs(x);
  575.     if (x==0.)
  576.         return(MAXRESULT);
  577.     if (x<2)
  578.         {
  579.         /* Eq 9.8.5 */
  580.         t=x/2.;
  581.         t2=t*t;
  582.         t4=t2*t2;
  583.         t8=t4*t4;
  584.         ans=-.57721566+.42278420*t2+.23069756*t4+.03488590*t4*t2
  585.             +.00262698*t8+.00010750*t8*t2+.00000740*t8*t4;
  586.         ans-=log(x/2)*bessi0(x);
  587.         }
  588.     else
  589.         {
  590.         /* Eq 9.8.6 */
  591.         t=2./x;
  592.         t2=t*t;
  593.         t4=t2*t2;
  594.         ans=1.25331414-.07832358*t+.02189568*t2-.01062446*t2*t
  595.             +.00587872*t4-.00251540*t4*t+.00053208*t4*t2;
  596.         ans*=exp(-x)/sqrt(x);
  597.         }
  598.     return(ans);
  599.     }
  600.  
  601.  
  602. /*
  603. ** bessk1   Calculate K1(x) from Abromowitz and Stegun, p. 379.
  604. **          Accuracy is better than 3e-7.
  605. */
  606. static double bessk1(double x)
  607.  
  608.     {
  609.     int     negx;
  610.     double  ans,t,t2,t4,t8;
  611.  
  612.     negx=(x<0.);
  613.     x=fabs(x);
  614.     if (x==0.)
  615.         return(MAXRESULT);
  616.     if (x<2)
  617.         {
  618.         /* Eq 9.8.7 */
  619.         t=x/2.;
  620.         t2=t*t;
  621.         t4=t2*t2;
  622.         t8=t4*t4;
  623.         ans=1+.15443144*t2-.67278579*t4-.18156897*t4*t2
  624.              -.01919402*t8-.00110404*t8*t2-.00004686*t8*t4;
  625.         ans=(ans/x)+log(x/2)*bessi1(x);
  626.         }
  627.     else
  628.         {
  629.         /* Eq 9.8.8 */
  630.         t=2./x;
  631.         t2=t*t;
  632.         t4=t2*t2;
  633.         ans=1.25331414+.23498619*t-.03655620*t2+.01504268*t2*t
  634.             -.00780353*t4+.00325614*t4*t-.00068245*t4*t2;
  635.         ans*=exp(-x)/sqrt(x);
  636.         }
  637.     return(negx ? -ans : ans);
  638.     }
  639.  
  640.  
  641. /*
  642. **
  643. **  dbessi      Calculates the value of Im'(x) using bessi and
  644. **              Bessel function recurrence relations.
  645. **  Im'(x) = m/x*Im(x) + I(m+1)(x)
  646. **
  647. */
  648. double dbessi(int m,double x)
  649.  
  650.     {
  651.     if (fabs(x)<TINYNUM)
  652.         x=x<0 ? -TINYNUM : TINYNUM;
  653.     return(m*bessi(m,x)/x+bessi(m+1,x));
  654.     }
  655.  
  656.  
  657. /*
  658. **
  659. **  dbessj      Calculates the value of Jm'(x) using bessj and
  660. **              Bessel function recurrence relations.
  661. **  Jm'(x) = m/x*Jm(x) - J(m+1)(x)
  662. **
  663. */
  664. double dbessj(int m,double x)
  665.  
  666.     {
  667.     if (fabs(x)<TINYNUM)
  668.         x=x<0 ? -TINYNUM : TINYNUM;
  669.     return(m*bessj(m,x)/x-bessj(m+1,x));
  670.     }
  671.  
  672.  
  673. /*
  674. **
  675. **  dbessk      Calculates the value of Km'(x) using bessk and
  676. **              Bessel function recurrence relations.
  677. **  Km'(x) = m/x*Km(x) - K(m+1)(x)
  678. **
  679. */
  680. double dbessk(int m,double x)
  681.  
  682.     {
  683.     if (fabs(x)<TINYNUM)
  684.         x=x<0 ? -TINYNUM : TINYNUM;
  685.     return(m*bessk(m,x)/x-bessk(m+1,x));
  686.     }
  687.  
  688.  
  689. /*
  690. **
  691. ** djroot       Calculates the nth non-zero root of the mth
  692. **              order bessel function's derivative, Jm'(x)
  693. **
  694. */
  695. double djroot(int m,int n)
  696.  
  697.    {
  698.    double  f0,f1,step,x0,x1;
  699.    int     flag,nroot;
  700.  
  701.    if (n<1)
  702.        return(0.);
  703.    x0=(double)m+.5;
  704.    step= n==1 ? SMALL_PI/4. : SMALL_PI;
  705.    x1=x0+step;
  706.    f0=dbessj(m,x0);
  707.    f1=dbessj(m,x1);
  708.    for (nroot=0 ;1; x0=x1,f0=f1,x1+=step,f1=dbessj(m,x1))
  709.        {
  710.        if (SGN(f0)*SGN(f1)<0)
  711.            {
  712.            nroot++;
  713.            if (nroot==n)
  714.                break;
  715.            if (nroot==n-1)
  716.                step/=4.;
  717.            }
  718.        }
  719.    return(root2(dbessj,m,x0,x1,ACCURACY,&flag));
  720.    }
  721.  
  722.  
  723. /*
  724. ** factorial
  725. */
  726. double factorial(int m)
  727.  
  728.     {
  729.     int     i;
  730.     double  ans;
  731.  
  732.     if (m<2)
  733.         return(1.);
  734.     ans=1.0;
  735.     for (i=2;i<=m;i++)
  736.         ans*=i;
  737.     return(ans);
  738.     }
  739.  
  740.  
  741. /*
  742. **
  743. **  jroot       Calculates the nth non-zero root of the mth
  744. **              order bessel function, Jm(x)
  745. **
  746. */
  747. double jroot(int m,int n)
  748.  
  749.    {
  750.    double  f0,f1,step,x0,x1;
  751.    int     flag,nroot;
  752.  
  753.    if (n<1)
  754.        return(0.);
  755.    x0=(double)m+2.;
  756.    step= n==1 ? SMALL_PI/4. : SMALL_PI;
  757.    x1=x0+step;
  758.    f0=bessj(m,x0);
  759.    f1=bessj(m,x1);
  760.    for (nroot=0 ;1; x0=x1,f0=f1,x1+=step,f1=bessj(m,x1))
  761.        {
  762.        if (SGN(f0)*SGN(f1)<0)
  763.            {
  764.            nroot++;
  765.            if (nroot==n)
  766.                break;
  767.            if (nroot==n-1)
  768.                step/=4.;
  769.            }
  770.        }
  771.    return(root2(bessj,m,x0,x1,ACCURACY,&flag));
  772.    }
  773.  
  774.  
  775. /*
  776. **
  777. ** root     Finds root of function F(x) between x0 and x1
  778. **          to precision: F(root) = 0. +/- xacc
  779. **
  780. ** Returns:  The root of the function.
  781. **
  782. ** Also returns status in "flag" variable.
  783. ** flag=0 if root found.
  784. ** flag=1 if root not found.  This could be due to NITER iterations being
  785. **        exceeded (a diverging solution pehaps?) or x0 and x1 not bracketing
  786. **        the root.
  787. **
  788. */
  789. double root(double(*F)(double),double x0,double x1,double xacc,int *flag)
  790.  
  791.    {
  792.    int     n;
  793.    double  f2,f3,f4,x2,x3,x4;
  794.  
  795.    (*flag)=1;
  796.    if (x0==x1)
  797.        return(0.);
  798.    x2=x0;
  799.    x3=x1;
  800.    f2=(*F)(x2);
  801.    f3=(*F)(x3);
  802.    if (SGN(f2)*SGN(f3)>0)
  803.        return(0.);
  804.    for (n=1;n<=NITER;n++)
  805.        {
  806.        if (n<=4)
  807.            x4=(x2+x3)/2.;
  808.        else
  809.            x4=(f3*x2-f2*x3)/(f3-f2);
  810.        f4=(*F)(x4);
  811.        if (ABS(f4)<xacc)
  812.            {
  813.            (*flag)=0;
  814.            return(x4);
  815.            }
  816.        if (SGN(f2)*SGN(f4)<0)
  817.            {
  818.            if (n>4 && ABS(f4)>ABS(f3))
  819.                return(0.);
  820.            f3=f4;
  821.            x3=x4;
  822.            continue;
  823.            }
  824.        if (n>4 && ABS(f4)>ABS(f2))
  825.            return(0.);
  826.        f2=f4;
  827.        x2=x4;
  828.        }
  829.    return(0.);
  830.    }
  831.  
  832.  
  833. /*
  834. **  root2   Finds root of function F(m,x).  Works the same way as root.
  835. */
  836. double root2(double(*F)(int,double),int m,double x0,double x1,double xacc,
  837.              int *flag)
  838.  
  839.    {
  840.    int     n;
  841.    double  f2,f3,f4,x2,x3,x4;
  842.  
  843.    (*flag)=1;
  844.    if (x0==x1)
  845.        return(0.);
  846.    x2=x0;
  847.    x3=x1;
  848.    f2=(*F)(m,x2);
  849.    f3=(*F)(m,x3);
  850.    if (SGN(f2)*SGN(f3)>0)
  851.        return(0.);
  852.    for (n=1;n<=NITER;n++)
  853.        {
  854.        if (n<=4)
  855.            x4=(x2+x3)/2.;
  856.        else
  857.            x4=(f3*x2-f2*x3)/(f3-f2);
  858.        f4=(*F)(m,x4);
  859.        if (ABS(f4)<xacc)
  860.            {
  861.            (*flag)=0;
  862.            return(x4);
  863.            }
  864.        if (SGN(f2)*SGN(f4)<0)
  865.            {
  866.            if (n>4 && ABS(f4)>ABS(f3))
  867.                return(0.);
  868.            f3=f4;
  869.            x3=x4;
  870.            continue;
  871.            }
  872.        if (n>4 && ABS(f4)>ABS(f2))
  873.            return(0.);
  874.        f2=f4;
  875.        x2=x4;
  876.        }
  877.    return(0.);
  878.    }
  879.