package graph;import java.lang.Math;import java.lang.ArithmeticException;/*******************************************************************************    Class  SpecialFunction            ******************************************************************************    Copyright (C) 1996 Leigh Brookshaw****    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.******************************************************************************    This class is an extension of java.lang.Math. It includes a number**    of special functions not found in the Math class.***************************************************************************//** * This class contains physical constants and special functions not found * in the java.lang.Math class. * Like the java.lang.Math class this class is final and cannot be * subclassed. * All physical constants are in cgs units. * <P> * <B>NOTE:</B> These special functions do not necessarily use the fastest * or most accurate algorithms. * * @version $Revision: 1.7 $, $Date: 1996/08/19 06:04:15 $ * @author Leigh Brookshaw  */public final class SpecialFunction extends Object {  /*  ** machine constants  */    private static final double MACHEP =  1.11022302462515654042E-16;    private static final double MAXLOG =  7.09782712893383996732E2;    private static final double MINLOG = -7.451332191019412076235E2;    private static final double MAXGAM = 171.624376956302725;    private static final double SQTPI  =  2.50662827463100050242E0;    private static final double SQRTH  =  7.07106781186547524401E-1;    private static final double LOGPI  =  1.14472988584940017414;  /*  ** Physical Constants in cgs Units  */  /**   * Boltzman Constant. Units erg/deg(K)   */    public static final double BOLTZMAN     = 1.3807e-16;  /**   * Elementary Charge. Units statcoulomb    */    public static final double ECHARGE      = 4.8032e-10;  /**   * Electron Mass. Units g   */    public static final double EMASS        = 9.1095e-28;  /**   * Proton Mass. Units g   */    public static final double PMASS        = 1.6726e-24;  /**   * Gravitational Constant. Units dyne-cm^2/g^2   */    public static final double GRAV         = 6.6720e-08;  /**   * Planck constant. Units erg-sec   */    public static final double PLANCK       = 6.6262e-27;  /**   * Speed of Light in a Vacuum. Units cm/sec   */    public static final double LIGHTSPEED   = 2.9979e10;  /**   * Stefan-Boltzman Constant. Units erg/cm^2-sec-deg^4   */    public static final double STEFANBOLTZ  = 5.6703e-5;  /**   * Avogadro Number. Units  1/mol   */    public static final double AVOGADRO     = 6.0220e23;  /**   * Gas Constant. Units erg/deg-mol   */    public static final double GASCONSTANT  = 8.3144e07;  /**   * Gravitational Acceleration at the Earths surface. Units cm/sec^2   */    public static final double GRAVACC      = 980.67;  /**   * Solar Mass. Units g   */    public static final double SOLARMASS    = 1.99e33;  /**   * Solar Radius. Units cm   */    public static final double SOLARRADIUS  = 6.96e10;  /**   * Solar Luminosity. Units erg/sec   */    public static final double SOLARLUM     = 3.90e33;  /**   * Solar Flux. Units erg/cm^2-sec   */    public static final double SOLARFLUX    = 6.41e10;  /**   * Astronomical Unit (radius of the Earth's orbit). Units cm   */    public static final double AU           = 1.50e13;    /**     * Don't let anyone instantiate this class.     */    private SpecialFunction() {}  /*  ** Function Methods  */  /**   * @param x a double value   * @return The log<sub>10</sub>   */    static public double log10(double x) throws ArithmeticException {         if( x <= 0.0 ) throw new ArithmeticException("range exception");         return Math.log(x)/2.30258509299404568401;    }  /**   * @param x a double value   * @return the hyperbolic cosine of the argument   */    static public double cosh(double x) throws ArithmeticException {      double a;      a = x;      if( a < 0.0 ) a = Math.abs(x);      a = Math.exp(a);      return 0.5*(a+1/a);    }  /**   * @param x a double value   * @return the hyperbolic sine of the argument   */    static public double sinh(double x) throws ArithmeticException {      double a;      if(x == 0.0) return x;      a = x;      if( a < 0.0 ) a = Math.abs(x);      a = Math.exp(a);      if( x < 0.0 )  return -0.5*(a-1/a);      else           return  0.5*(a-1/a);    }  /**   * @param x a double value   * @return the hyperbolic tangent of the argument   */    static public double tanh(double x) throws ArithmeticException {      double a;      if( x == 0.0 ) return x;      a = x;      if( a < 0.0 ) a = Math.abs(x);      a = Math.exp(2.0*a);      if(x < 0.0 ) return -( 1.0-2.0/(a+1.0) );      else         return  ( 1.0-2.0/(a+1.0) );    }  /**   * @param x a double value   * @return the hyperbolic arc cosine of the argument   */    static public double acosh(double x) throws ArithmeticException {      if( x < 1.0 ) throw new ArithmeticException("range exception");      return Math.log( x + Math.sqrt(x*x-1));    }  /**   * @param x a double value   * @return the hyperbolic arc sine of the argument   */    static public double asinh(double xx) throws ArithmeticException {      double x;      int sign;      if(xx == 0.0) return xx;      if( xx < 0.0 ) {                      sign = -1;                      x = -xx;      } else {                      sign = 1;                      x = xx;      }      return sign*Math.log( x + Math.sqrt(x*x+1));    }  /**   * @param x a double value   * @return the hyperbolic arc tangent of the argument   */    static public double atanh(double x) throws ArithmeticException {      if( x > 1.0 || x < -1.0 ) throw                          new ArithmeticException("range exception");      return 0.5 * Math.log( (1.0+x)/(1.0-x) );    }  /**   * @param x a double value   * @return the Bessel function of order 0 of the argument.   */    static public double j0(double x) throws ArithmeticException {        double ax;        if( (ax=Math.abs(x)) < 8.0 ) {           double y=x*x;           double ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7                       +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));           double ans2=57568490411.0+y*(1029532985.0+y*(9494680.718                       +y*(59272.64853+y*(267.8532712+y*1.0))));           return ans1/ans2;        } else {           double z=8.0/ax;           double y=z*z;           double xx=ax-0.785398164;           double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4                       +y*(-0.2073370639e-5+y*0.2093887211e-6)));           double ans2 = -0.1562499995e-1+y*(0.1430488765e-3                       +y*(-0.6911147651e-5+y*(0.7621095161e-6                       -y*0.934935152e-7)));                      return Math.sqrt(0.636619772/ax)*                  (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);	}    }  /**   * @param x a double value   * @return the Bessel function of order 1 of the argument.   */    static public double j1(double x) throws ArithmeticException {      double ax;      double y;      double ans1, ans2;      if ((ax=Math.abs(x)) < 8.0) {         y=x*x;         ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1               +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));         ans2=144725228442.0+y*(2300535178.0+y*(18583304.74               +y*(99447.43394+y*(376.9991397+y*1.0))));         return ans1/ans2;       } else {         double z=8.0/ax;         double xx=ax-2.356194491;         y=z*z;         ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4              +y*(0.2457520174e-5+y*(-0.240337019e-6))));         ans2=0.04687499995+y*(-0.2002690873e-3              +y*(0.8449199096e-5+y*(-0.88228987e-6              +y*0.105787412e-6)));         double ans=Math.sqrt(0.636619772/ax)*                   (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);         if (x < 0.0) ans = -ans;         return ans;       }    }  /**   * @param n integer order   * @param x a double value   * @return the Bessel function of order n of the argument.   */    static public double jn(int n, double x) throws ArithmeticException {       int j,m;       double ax,bj,bjm,bjp,sum,tox,ans;       boolean jsum;       double ACC   = 40.0;       double BIGNO = 1.0e+10;       double BIGNI = 1.0e-10;       if(n == 0) return j0(x);       if(n == 1) return j1(x);       ax=Math.abs(x);       if(ax == 0.0)  return 0.0;       else        if (ax > (double)n) {         tox=2.0/ax;         bjm=j0(ax);         bj=j1(ax);         for (j=1;j<n;j++) {            bjp=j*tox*bj-bjm;            bjm=bj;            bj=bjp;         }         ans=bj;       } else {         tox=2.0/ax;         m=2*((n+(int)Math.sqrt(ACC*n))/2);         jsum=false;         bjp=ans=sum=0.0;         bj=1.0;         for (j=m;j>0;j--) {            bjm=j*tox*bj-bjp;            bjp=bj;            bj=bjm;            if (Math.abs(bj) > BIGNO) {               bj *= BIGNI;               bjp *= BIGNI;               ans *= BIGNI;               sum *= BIGNI;            }            if (jsum) sum += bj;            jsum=!jsum;            if (j == n) ans=bjp;          }          sum=2.0*sum-bj;          ans /= sum;       }       return  x < 0.0 && n%2 == 1 ? -ans : ans;   }  /**   * @param x a double value   * @return the Bessel function of the second kind,    *          of order 0 of the argument.   */   static public double y0(double x) throws ArithmeticException {      if (x < 8.0) {         double y=x*x;         double ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6                        +y*(10879881.29+y*(-86327.92757+y*228.4622733))));         double ans2=40076544269.0+y*(745249964.8+y*(7189466.438                        +y*(47447.26470+y*(226.1030244+y*1.0))));         return (ans1/ans2)+0.636619772*j0(x)*Math.log(x);      } else {         double z=8.0/x;         double y=z*z;         double xx=x-0.785398164;         double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4                     +y*(-0.2073370639e-5+y*0.2093887211e-6)));         double ans2 = -0.1562499995e-1+y*(0.1430488765e-3                      +y*(-0.6911147651e-5+y*(0.7621095161e-6                      +y*(-0.934945152e-7))));         return Math.sqrt(0.636619772/x)*                (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);      }   }  /**   * @param x a double value   * @return the Bessel function of the second kind,   *  of order 1 of the argument.   */   static public double y1(double x) throws ArithmeticException {         if (x < 8.0) {         double y=x*x;         double ans1=x*(-0.4900604943e13+y*(0.1275274390e13                     +y*(-0.5153438139e11+y*(0.7349264551e9                     +y*(-0.4237922726e7+y*0.8511937935e4)))));         double ans2=0.2499580570e14+y*(0.4244419664e12                     +y*(0.3733650367e10+y*(0.2245904002e8                     +y*(0.1020426050e6+y*(0.3549632885e3+y)))));         return (ans1/ans2)+0.636619772*(j1(x)*Math.log(x)-1.0/x);      } else {         double z=8.0/x;         double y=z*z;         double xx=x-2.356194491;         double ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4                     +y*(0.2457520174e-5+y*(-0.240337019e-6))));         double ans2=0.04687499995+y*(-0.2002690873e-3                     +y*(0.8449199096e-5+y*(-0.88228987e-6                     +y*0.105787412e-6)));         return Math.sqrt(0.636619772/x)*                (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);      }   }  /**   * @param n integer order   * @param x a double value   * @return the Bessel function of the second kind,   *    of order n of the argument.   */    static public double yn(int n, double x) throws ArithmeticException {       double by,bym,byp,tox;       if(n == 0) return y0(x);       if(n == 1) return y1(x);       tox=2.0/x;       by=y1(x);       bym=y0(x);       for (int j=1;j<n;j++) {         byp=j*tox*by-bym;         bym=by;         by=byp;       }       return by;    }  /**   * @param x a double value   * @return the factorial of the argument   */     static public double fac(double x) throws ArithmeticException {        double d = Math.abs(x);        if(Math.floor(d) == d) return (double)fac( (int)x );        else                   return gamma(x+1.0);     }  /**