specialfunction.java

java 作图的程序

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);
     }

  /**