虚数

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace Scientific
{
    public struct cplx
    {
        public double x;
        public double y;
        #region [ 构造函数 ]
        public cplx(double a, double b)
        {
            x = a;
            y = b;
        }
        #endregion
        #region [ Operator  + - x /  ]
        public static cplx operator +(cplx A, cplx B)
        {
            return new cplx(A.x + B.x, A.y + B.y);
        }
        public static cplx operator +(double a, cplx B)
        {
            return new cplx(a + B.x, B.y);
        }
        public static cplx operator +(cplx B, double a)
        {
            return new cplx(a + B.x, B.y);
        }
        public static cplx operator -(cplx A, cplx B)
        {
            return new cplx(A.x - B.x, A.y - B.y);
        }
        public static cplx operator -(cplx A, double b)
        {
            return new cplx(A.x - b, A.y);
        }
        public static cplx operator -(double a, cplx B)
        {
            return new cplx(a - B.x, 0 - B.y);
        }
        public static cplx operator *(cplx A, cplx B)
        {
            return new cplx((A.x * B.x) - (A.y * B.y), (A.y * B.x) + (A.x * B.y));
        }
        public static cplx operator *(double a, cplx B)
        {
            return new cplx(a * B.x, a * B.y);
        }
        public static cplx operator *(cplx B, double a)
        {
            return new cplx(a * B.x, a * B.y);
        }
        public static cplx operator /(cplx A, cplx B)
        {
            cplx cp = new cplx(1, 0);   //A=0;B=0
            double denominator = B.x * B.x + B.y * B.y;
            if (denominator > 1e-300)   //B!=0
            {
                cp.x = (A.x * B.x + A.y * B.y) / denominator;
                cp.y = (A.y * B.x - A.x * B.y) / denominator;
            }
            else                       //B=0
            {
                if (A.x * A.x + A.y * A.y > 1e-300) //B=0; A!=0
                {
                    cp.x = 1e300;
                    cp.y = 1e300;
                }
            }
            return cp;
        }
        public static cplx operator /(cplx A, double b)
        {
            return A / (new cplx(b, 0));
        }
        public static cplx operator /(double b, cplx A)
        {
            return (new cplx(b, 0)) / A;
        }
        #endregion
        #region  [ Math Functions ]
        public cplx Pow(double n)
        {
            cplx mp = ToMagPhase();  //Convert to Mag and Phase
            mp.x = Math.Pow(10, mp.x / 20);
            double MAG = Math.Pow(mp.x, n);            //Power MAG
            double ANG = mp.y * n;                     //Power ANG
            cplx cp = new cplx(1, 0);
            cp.x = MAG * Math.Cos(ANG);       //Real
            cp.y = MAG * Math.Sin(ANG);       //Imaginary
            return cp;
        }
        public static cplx Exp_j(double N)
        {
            double a = Math.Cos(N);
            double b = Math.Sin(N);
            return new cplx(a, b);
        }
        public cplx ToMagPhase()
        {
            cplx cp = new cplx(1, 0);
            cp.x = 20 * Math.Log10(Math.Sqrt(x * x + y * y));     //Mag
            cp.y = Math.Atan2(y, x);              //Phase in Arc
            return cp;
        }
        public cplx ToRealIm()
        {
            cplx cp = new cplx(1, 0);
            cp.x = x * Math.Cos(y);       //Real
            cp.y = x * Math.Sin(y);       //Imaginary
            return cp;
        }
        #endregion
        #region [ ToString ]
        public override string ToString()
        {
            return x.ToString() + " +j (" + y.ToString() + ")";
        }
        #endregion
    }
    public class Calc
    {
        public Calc()
        {
        }
        #region [ Base Functions ]
        public static double Calc_Arc(double Deg)
        {
            double Arc = Deg * (Math.PI / 180);
            return Arc;
        }
        public static double Calc_Deg(double Arc)
        {
            double Deg = Arc * 180 / Math.PI;
            return Deg;
        }
        public static double Calc_10Log(double abs)
        {
            return 10D * Math.Log10(Math.Abs(abs));
        }
        public static double Calc_20Log(double abs)
        {
            return 20D * Math.Log10(Math.Abs(abs));
        }
        public static double Calc_10Exp(double rl)
        {
            return Math.Pow(10, rl / 10D);
        }
        public static double Calc_20Exp(double rl)
        {
            return Math.Pow(10, rl / 20D);
        }
        #endregion
        #region [ Bessel ]
        private static void besselasympt0(double x, ref double pzero, ref double qzero)
        {
            double xsq = 0;
            double p2 = 0;
            double q2 = 0;
            double p3 = 0;
            double q3 = 0;
            xsq = 64.0 / (x * x);
            p2 = 0.0;
            p2 = 2485.271928957404011288128951 + xsq * p2;
            p2 = 153982.6532623911470917825993 + xsq * p2;
            p2 = 2016135.283049983642487182349 + xsq * p2;
            p2 = 8413041.456550439208464315611 + xsq * p2;
            p2 = 12332384.76817638145232406055 + xsq * p2;
            p2 = 5393485.083869438325262122897 + xsq * p2;
            q2 = 1.0;
            q2 = 2615.700736920839685159081813 + xsq * q2;
            q2 = 156001.7276940030940592769933 + xsq * q2;
            q2 = 2025066.801570134013891035236 + xsq * q2;
            q2 = 8426449.050629797331554404810 + xsq * q2;
            q2 = 12338310.22786324960844856182 + xsq * q2;
            q2 = 5393485.083869438325560444960 + xsq * q2;
            p3 = -0.0;
            p3 = -4.887199395841261531199129300 + xsq * p3;
            p3 = -226.2630641933704113967255053 + xsq * p3;
            p3 = -2365.956170779108192723612816 + xsq * p3;
            p3 = -8239.066313485606568803548860 + xsq * p3;
            p3 = -10381.41698748464093880530341 + xsq * p3;
            p3 = -3984.617357595222463506790588 + xsq * p3;
            q3 = 1.0;
            q3 = 408.7714673983499223402830260 + xsq * q3;
            q3 = 15704.89191515395519392882766 + xsq * q3;
            q3 = 156021.3206679291652539287109 + xsq * q3;
            q3 = 533291.3634216897168722255057 + xsq * q3;
            q3 = 666745.4239319826986004038103 + xsq * q3;
            q3 = 255015.5108860942382983170882 + xsq * q3;
            pzero = p2 / q2;
            qzero = 8 * p3 / q3 / x;
        }
        private static void besselasympt1(double x, ref double pzero, ref double qzero)
        {
            double xsq = 0;
            double p2 = 0;
            double q2 = 0;
            double p3 = 0;
            double q3 = 0;
            xsq = 64.0 / (x * x);
            p2 = -1611.616644324610116477412898;
            p2 = -109824.0554345934672737413139 + xsq * p2;
            p2 = -1523529.351181137383255105722 + xsq * p2;
            p2 = -6603373.248364939109255245434 + xsq * p2;
            p2 = -9942246.505077641195658377899 + xsq * p2;
            p2 = -4435757.816794127857114720794 + xsq * p2;
            q2 = 1.0;
            q2 = -1455.009440190496182453565068 + xsq * q2;
            q2 = -107263.8599110382011903063867 + xsq * q2;
            q2 = -1511809.506634160881644546358 + xsq * q2;
            q2 = -6585339.479723087072826915069 + xsq * q2;
            q2 = -9934124.389934585658967556309 + xsq * q2;
            q2 = -4435757.816794127856828016962 + xsq * q2;
            p3 = 35.26513384663603218592175580;
            p3 = 1706.375429020768002061283546 + xsq * p3;
            p3 = 18494.26287322386679652009819 + xsq * p3;
            p3 = 66178.83658127083517939992166 + xsq * p3;
            p3 = 85145.16067533570196555001171 + xsq * p3;
            p3 = 33220.91340985722351859704442 + xsq * p3;
            q3 = 1.0;
            q3 = 863.8367769604990967475517183 + xsq * q3;
            q3 = 37890.22974577220264142952256 + xsq * q3;
            q3 = 400294.4358226697511708610813 + xsq * q3;
            q3 = 1419460.669603720892855755253 + xsq * q3;
            q3 = 1819458.042243997298924553839 + xsq * q3;
            q3 = 708712.8194102874357377502472 + xsq * q3;
            pzero = p2 / q2;
            qzero = 8 * p3 / q3 / x;
        }
        public static double besselj0(double x)
        {
            double result = 0;
            double xsq = 0;
            double nn = 0;
            double pzero = 0;
            double qzero = 0;
            double p1 = 0;
            double q1 = 0;
            if (x < 0)
            {
                x = -x;
            }
            if (x > 8.0)
            {
                besselasympt0(x, ref pzero, ref qzero);
                nn = x - Math.PI / 4;
                result = Math.Sqrt(2 / Math.PI / x) * (pzero * Math.Cos(nn) - qzero * Math.Sin(nn));
                return result;
            }
            xsq = x * x;
            p1 = 26857.86856980014981415848441;
            p1 = -40504123.71833132706360663322 + xsq * p1;
            p1 = 25071582855.36881945555156435 + xsq * p1;
            p1 = -8085222034853.793871199468171 + xsq * p1;
            p1 = 1434354939140344.111664316553 + xsq * p1;
            p1 = -136762035308817138.6865416609 + xsq * p1;
            p1 = 6382059341072356562.289432465 + xsq * p1;
            p1 = -117915762910761053603.8440800 + xsq * p1;
            p1 = 493378725179413356181.6813446 + xsq * p1;
            q1 = 1.0;
            q1 = 1363.063652328970604442810507 + xsq * q1;
            q1 = 1114636.098462985378182402543 + xsq * q1;
            q1 = 669998767.2982239671814028660 + xsq * q1;
            q1 = 312304311494.1213172572469442 + xsq * q1;
            q1 = 112775673967979.8507056031594 + xsq * q1;
            q1 = 30246356167094626.98627330784 + xsq * q1;
            q1 = 5428918384092285160.200195092 + xsq * q1;
            q1 = 493378725179413356211.3278438 + xsq * q1;
            result = p1 / q1;
            return result;
        }
        public static double besselj1(double x)
        {
            double result = 0;
            double s = 0;
            double xsq = 0;
            double nn = 0;
            double pzero = 0;
            double qzero = 0;
            double p1 = 0;
            double q1 = 0;
            s = Math.Sign(x);
            if (x < 0)
            {
                x = -x;
            }
            if (x > 8.0)
            {
                besselasympt1(x, ref pzero, ref qzero);
                nn = x - 3 * Math.PI / 4;
                result = Math.Sqrt(2 / Math.PI / x) * (pzero * Math.Cos(nn) - qzero * Math.Sin(nn));
                if (s < 0)
                {
                    result = -result;
                }
                return result;
            }
            xsq = x * x;
            p1 = 2701.122710892323414856790990;
            p1 = -4695753.530642995859767162166 + xsq * p1;
            p1 = 3413234182.301700539091292655 + xsq * p1;
            p1 = -1322983480332.126453125473247 + xsq * p1;
            p1 = 290879526383477.5409737601689 + xsq * p1;
            p1 = -35888175699101060.50743641413 + xsq * p1;
            p1 = 2316433580634002297.931815435 + xsq * p1;
            p1 = -66721065689249162980.20941484 + xsq * p1;
            p1 = 581199354001606143928.050809 + xsq * p1;
            q1 = 1.0;
            q1 = 1606.931573481487801970916749 + xsq * q1;
            q1 = 1501793.594998585505921097578 + xsq * q1;
            q1 = 1013863514.358673989967045588 + xsq * q1;
            q1 = 524371026216.7649715406728642 + xsq * q1;
            q1 = 208166122130760.7351240184229 + xsq * q1;
            q1 = 60920613989175217.46105196863 + xsq * q1;
            q1 = 11857707121903209998.37113348 + xsq * q1;
            q1 = 1162398708003212287858.529400 + xsq * q1;
            result = s * x * p1 / q1;
            return result;
        }
        public static double besseljn(int n, double x)
        {
            double result = 0;
            double pkm2 = 0;
            double pkm1 = 0;
            double pk = 0;
            double xk = 0;
            double r = 0;
            double ans = 0;
            int k = 0;
            int sg = 0;
            if (n < 0)
            {
                n = -n;
                if (n % 2 == 0)
                {
                    sg = 1;
                }
                else
                {
                    sg = -1;
                }
            }
            else
            {
                sg = 1;
            }
            if (x < 0)
            {
                if (n % 2 != 0)
                {
                    sg = -sg;
                }
                x = -x;
            }
            if (n == 0)
            {
                result = sg * besselj0(x);
                return result;
            }
            if (n == 1)
            {
                result = sg * besselj1(x);
                return result;
            }
            if (n == 2)
            {
                if (x == 0)
                {
                    result = 0;
                }
                else
                {
                    result = sg * (2.0 * besselj1(x) / x - besselj0(x));
                }
                return result;
            }
            if (x < 5E-16)
            {
                result = 0;
                return result;
            }
            k = 53;
            pk = 2 * (n + k);
            ans = pk;
            xk = x * x;
            do
            {
                pk = pk - 2.0;
                ans = pk - xk / ans;
                k = k - 1;
            }
            while (k != 0);
            ans = x / ans;
            pk = 1.0;
            pkm1 = 1.0 / ans;
            k = n - 1;
            r = 2 * k;
            do
            {
                pkm2 = (pkm1 * r - pk * x) / x;
                pk = pkm1;
                pkm1 = pkm2;
                r = r - 2.0;
                k = k - 1;
            }
            while (k != 0);
            if (Math.Abs(pk) > Math.Abs(pkm1))
            {
                ans = besselj1(x) / pk;
            }
            else
            {
                ans = besselj0(x) / pkm1;
            }
            result = sg * ans;
            return result;
        }
        public static double Jn(int n, double x)
        {
            double y = 0D;
            y = besseljn(n, x);
            return y;
        }
        public static double Jd(int n, double x)
        {
            double y = 0D;
            y = (Jn(n - 1, x) - Jn(n + 1, x)) / 2;
            return y;
        }
        #endregion
    }
}