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
}
}
虚数
最新推荐文章于 2024-01-26 09:26:04 发布