forked from tsai/mathnet-numerics
diluculo
8 years ago
9 changed files with 9075 additions and 0 deletions
@ -0,0 +1,92 @@ |
|||
using NUnit.Framework; |
|||
using System.Numerics; |
|||
|
|||
namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests |
|||
{ |
|||
/// <summary>
|
|||
/// Airy functions tests.
|
|||
/// </summary>
|
|||
[TestFixture, Category("Functions")] |
|||
public class AiryTests |
|||
{ |
|||
[TestCase(0.0, 0.0, 0.35502805388781723926, 0.0, 14)] |
|||
[TestCase(1.0, 0.0, 0.13529241631288141552, 0.0, 14)] |
|||
[TestCase(-1.0, 0.0, 0.53556088329235211880, 0.0, 14)] |
|||
[TestCase(0.0, 1.0, 0.33149330543214118898, -0.31744985896844377348, 14)] |
|||
[TestCase(0.0, -1.0, 0.33149330543214118898, 0.31744985896844377348, 14)] |
|||
[TestCase(1.0, 1.0, 0.060458308371838149197, -0.15188956587718140235, 13)] |
|||
[TestCase(-1.0, -1.0, 0.82211742655527259396, 0.11996634266442434389, 14)] |
|||
[TestCase(10.0, 5.0, -7.0165968356580205348E-10, 2.7637938765570892385E-10, 14)] |
|||
[TestCase(-10.0, -5.0, 1.2329510105552886439E6, 502435.45036313325712, 13)] |
|||
[TestCase(100.0, 100.0, 2.9099582462207032076E-188, 2.3530135917061787560E-188, 14)] |
|||
[TestCase(32.0, -64.0, 5.3014568355995704254E14, 2.0523039181737934724E13, 11)] |
|||
public void AiryAiApprox(double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.AiryAi(new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
[TestCase(0.0, 0.0, -0.25881940379280679841, 0.0, 14)] |
|||
[TestCase(1.0, 0.0, -0.15914744129679321279, 0.0, 14)] |
|||
[TestCase(-1.0, 0.0, -0.010160567116645209395, 0.0, 13)] |
|||
[TestCase(0.0, 1.0, -0.43249265984180709931, 0.098047856229243232384, 14)] |
|||
[TestCase(0.0, -1.0, -0.43249265984180709931, -0.098047856229243232384, 14)] |
|||
[TestCase(1.0, 1.0, -0.13062795349964751771, 0.16306759644932391574, 13)] |
|||
[TestCase(-1.0, -1.0, -0.37906047922683349624, 0.60450013086224607164, 14)] |
|||
[TestCase(10.0, 5.0, 2.5069533633744176669E-9, -3.7264734152221837760E-10, 13)] |
|||
[TestCase(-10.0, -5.0, -2.5522092133303029354E6, 3.6244486331968648174E6, 13)] |
|||
[TestCase(100.0, 100.0, -2.1269500702792576318E-187, -3.9094417620496355697E-187, 14)] |
|||
[TestCase(32.0, -64.0, -3.9067669879175239067E15, 2.2082697395233711387E15, 12)] |
|||
public void AiryAiPrimeApprox(double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.AiryAiPrime(new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
[TestCase(0.0, 0.0, 0.61492662744600073515, 0.0, 14)] |
|||
[TestCase(1.0, 0.0, 1.2074235949528712594, 0.0, 14)] |
|||
[TestCase(-1.0, 0.0, 0.10399738949694461189, 0.0, 14)] |
|||
[TestCase(0.0, 1.0, 0.64885820833039494458, 0.34495863476804837025, 14)] |
|||
[TestCase(0.0, -1.0, 0.64885820833039494458, -0.34495863476804837025, 14)] |
|||
[TestCase(1.0, 1.0, 0.71665807338276843179, 0.61988929040084476435, 13)] |
|||
[TestCase(-1.0, -1.0, 0.21429040153487357398, -0.67391692372270520960, 13)] |
|||
[TestCase(10.0, 5.0, -6.2471332619252646778E7, -9.0137631732829873100E6, 13)] |
|||
[TestCase(-10.0, -5.0, 502435.45036315399060, -1.2329510105552595201E6, 13)] |
|||
[TestCase(100.0, 100.0, 1.7086751714463652039E185, -3.1416590020830804578E185, 11)] |
|||
[TestCase(32.0, -64.0, 2.0523039181737934724E13, -5.3014568355995704254E14, 11)] |
|||
public void AiryBiApprox(double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.AiryBi(new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
[TestCase(0.0, 0.0, 0.44828835735382635791, 0.0, 14)] |
|||
[TestCase(1.0, 0.0, 0.93243593339277563296, 0.0, 14)] |
|||
[TestCase(-1.0, 0.0, 0.59237562642279235082, 0.0, 14)] |
|||
[TestCase(0.0, 1.0, 0.13502664671081897270, -0.12883738678125487904, 14)] |
|||
[TestCase(0.0, -1.0, 0.13502664671081897270, 0.12883738678125487904, 14)] |
|||
[TestCase(1.0, 1.0, 0.075662844174965992918, 0.78370099878545527505, 13)] |
|||
[TestCase(-1.0, -1.0, 0.83447348852278263690, 0.34652606326682852855, 14)] |
|||
[TestCase(10.0, 5.0, -1.9502117576621140236E8, -7.7790596402473222951E7, 14)] |
|||
[TestCase(-10.0, -5.0, 3.6244486331969762244E6, 2.5522092133302581994E6, 13)] |
|||
[TestCase(100.0, 100.0, 3.3072107798533888664E186, -2.6734837736904900245E186, 12)] |
|||
[TestCase(32.0, -64.0, 2.2082697395233711387E15, 3.9067669879175239067E15, 12)] |
|||
public void AiryBiPrimeApprox(double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.AiryBiPrime(new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,245 @@ |
|||
using MathNet.Numerics.UnitTests; |
|||
using NUnit.Framework; |
|||
using System; |
|||
using System.Numerics; |
|||
|
|||
namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests |
|||
{ |
|||
/// <summary>
|
|||
/// Bessel functions tests.
|
|||
/// </summary>
|
|||
[TestFixture, Category("Functions")] |
|||
public class BesselTests |
|||
{ |
|||
#region Bessel J
|
|||
|
|||
[Test] |
|||
public void BesselJ0Approx([Range(-3, 3, 0.25)] double x) |
|||
{ |
|||
// Approx by Abramowitz/Stegun 9.4.1
|
|||
Assert.AreEqual(Evaluate.Polynomial(x / 3.0, 1.0, 0.0, -2.2499997, 0.0, 1.2656208, 0.0, -0.3163866, 0.0, 0.0444479, 0.0, -0.0039444, 0.0, 0.0002100), SpecialFunctions.BesselJ(0, x), 1e-7); |
|||
} |
|||
|
|||
[TestCase(0, 0.0, 0.0, 1.0000000000000000000, 0.0000000000000000000, 14)] |
|||
[TestCase(0, 0.25, 0.0, 0.98443592929585270492, 0.0000000000000000000, 14)] |
|||
[TestCase(0, 0.5, 0.0, 0.93846980724081290423, 0.0000000000000000000, 14)] |
|||
[TestCase(0, 1.0, 0.0, 0.76519768655796655145, 0.0000000000000000000, 14)] |
|||
[TestCase(0, -1.0, 0.0, 0.76519768655796655145, 0.0000000000000000000, 14)] |
|||
[TestCase(0, 0.0, 1.0, 1.2660658777520083356, 0.0000000000000000000, 14)] |
|||
[TestCase(0, 0.0, -1.0, 1.2660658777520083356, 0.0000000000000000000, 14)] |
|||
[TestCase(0, 1.0, 1.0, 0.93760847680602927660, -0.49652994760912213217, 14)] |
|||
[TestCase(0, -1.0, -1.0, 0.93760847680602927660, -0.49652994760912213217, 14)] |
|||
[TestCase(0, 10.0, 5.0, -17.789591129450371518, 0.20071161672120485098, 13)] |
|||
[TestCase(0, -10.0, -5.0, -17.789591129450371518, 0.20071161672120485098, 13)] |
|||
[TestCase(0.001, 1.0, 1.0, 0.93830680591878982122, -0.49541402042792734422, 14)] |
|||
[TestCase(1, 0.0, -10.0, 0.0000000000000000000, -2670.9883037012546543, 12)] |
|||
[TestCase(1, 10.0, 5.0, -0.92143143564032744340, -17.436439610477507252, 13)] |
|||
[TestCase(-1, 10.0, 5.0, 0.92143143564032744339, 17.436439610477507252, 13)] |
|||
[TestCase(1, 100.0, 100.0, -7.1708320564209428560E41, 5.4401809306071216972E41, 14)] |
|||
[TestCase(2, 4.0, 0.0, 0.36412814585207280421, 0.0000000000000000000, 14)] |
|||
[TestCase(2, 32.0, -64.0, -2.6838844111822281863E26, -1.0228528267491243961E26, 14)] |
|||
[TestCase(100, 1.0, 1.0, -9.5168076700036928783E-174, 4.7113294059604970982E-176, 14)] |
|||
public void ComplexBesselJnExact(double v, double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.BesselJ(v, new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region Bessel Y
|
|||
|
|||
[Test] |
|||
public void BesselY0Approx([Range(0.25, 3, 0.25)] double x) |
|||
{ |
|||
// Approx by Abramowitz/Stegun 9.4.2
|
|||
Assert.AreEqual( |
|||
Evaluate.Polynomial(x / 3.0, 2.0 / Math.PI * Math.Log(x / 2.0) * SpecialFunctions.BesselJ(0, x) + 0.36746691, 0.0, 0.60559366, 0.0, -0.74350384, 0.0, 0.25300117, 0.0, -0.04261214, 0.0, 0.00427916, 0.0, -0.00024846), |
|||
SpecialFunctions.BesselY(0, x), 1e-7); |
|||
} |
|||
|
|||
[TestCase(0, 0.0, 0.0, double.NegativeInfinity, 0.0, 14)] |
|||
[TestCase(0, 1.0, 0.0, 0.088256964215676957983, 0.0, 14)] |
|||
[TestCase(0, 0.0, 1.0, -0.26803248203398854876, 1.2660658777520083356, 14)] |
|||
[TestCase(0, 1.0, 1.0, 0.44547448893603251403, 0.71015858200373452118, 14)] |
|||
[TestCase(0, 10.0, 5.0, -0.20001363358737298600, -17.788152066750989245, 13)] |
|||
[TestCase(0, -10.0, -5.0, 0.20140959985503671596, 17.791030192149753792, 13)] |
|||
[TestCase(0.001, 1.0, 1.0, 0.44400137479127562106, 0.71093732860971252181, 14)] |
|||
[TestCase(1, 10.0, 5.0, 17.437935555053368400, -0.92208733599150420146, 13)] |
|||
[TestCase(-1, 10.0, 5.0, -17.437935555053368400, 0.92208733599150420146, 13)] |
|||
[TestCase(1, 100.0, 100.0, -5.4401809306071216972E41, -7.1708320564209428560E41, 14)] |
|||
[TestCase(2, 4.0, 0.0, 0.21590359460361499453, 0.0, 14)] |
|||
[TestCase(2, 32.0, -64.0, -1.0228528267491243961E26, 2.6838844111822281863E26, 14)] |
|||
[TestCase(100, 1.0, 1.0, 3.3446291442187872047E170, 1.6892209981554826575E168, 10)] |
|||
public void ComplexBesselYnExact(double v, double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.BesselY(v, new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region Modified Bessel I
|
|||
|
|||
[Test] |
|||
public void BesselI0Approx([Range(-3.75, 3.75, 0.25)] double x) |
|||
{ |
|||
// Approx by Abramowitz/Stegun 9.8.1
|
|||
Assert.AreEqual(Evaluate.Polynomial(x / 3.75, 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813), SpecialFunctions.BesselI(0, x), 1e-7); |
|||
} |
|||
|
|||
[TestCase(0.0, 1.0)] |
|||
[TestCase(0.005, 1.000006250009766)] |
|||
[TestCase(0.5, 1.063483370741324)] |
|||
[TestCase(1.5, 1.646723189772891)] |
|||
[TestCase(10.0, 2815.716628466254)] |
|||
[TestCase(100.0, 1.073751707131074e+42)] |
|||
[TestCase(-0.005, 1.000006250009766)] |
|||
[TestCase(-10.0, 2815.716628466254)] |
|||
public void BesselI0Exact(double x, double expected) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative(expected, SpecialFunctions.BesselI(0, x), 14); |
|||
} |
|||
|
|||
[Test] |
|||
public void BesselI1Approx([Range(-3.75, 3.75, 0.25)] double x) |
|||
{ |
|||
// Approx by Abramowitz/Stegun 9.8.3
|
|||
Assert.AreEqual(Evaluate.Polynomial(x / 3.75, 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411) * x, SpecialFunctions.BesselI(1, x), 1e-8); |
|||
} |
|||
|
|||
[TestCase(0.0, 0.0)] |
|||
[TestCase(0.005, 0.002500007812508138)] |
|||
[TestCase(0.5, 0.2578943053908963)] |
|||
[TestCase(1.5, 0.9816664285779076)] |
|||
[TestCase(10.0, 2670.988303701255)] |
|||
[TestCase(100.0, 1.068369390338162e+42)] |
|||
[TestCase(-0.005, -0.002500007812508138)] |
|||
[TestCase(-10.0, -2670.988303701255)] |
|||
public void BesselI1Exact(double x, double expected) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative(expected, SpecialFunctions.BesselI(1, x), 14); |
|||
} |
|||
|
|||
[TestCase(0, +10, 2815.7166284662544715)] |
|||
[TestCase(0, -10, 2815.7166284662544715)] |
|||
[TestCase(1, +10, 2670.9883037012546543)] |
|||
[TestCase(1, -10, -2670.9883037012546543)] |
|||
[TestCase(2, +10, 2281.5189677260035406)] |
|||
[TestCase(2, -10, 2281.5189677260035406)] |
|||
[TestCase(3, +10, 1758.3807166108532381)] |
|||
[TestCase(3, -10, -1758.3807166108532381)] |
|||
[TestCase(4, +10, 1226.4905377594915977)] |
|||
[TestCase(4, -10, 1226.4905377594915977)] |
|||
[TestCase(5, +10, 777.18828640325995991)] |
|||
[TestCase(5, -10, -777.18828640325995991)] |
|||
public void BesselInExact(int n, double x, double expected) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative(expected, SpecialFunctions.BesselI(n, x), 14); |
|||
} |
|||
|
|||
[TestCase(0, 0.0, 0.0, 1.0, 0.0, 14)] |
|||
[TestCase(0, 1.0, 0.0, 1.2660658777520083356, 0.0, 14)] |
|||
[TestCase(0, -1.0, 0.0, 1.2660658777520083356, 0.0, 14)] |
|||
[TestCase(0, 0.0, 1.0, 0.7651976865579665514497, 0.0, 14)] |
|||
[TestCase(0, 0.0, -1.0, 0.7651976865579665514497, 0.0, 14)] |
|||
[TestCase(0, 1.0, 1.0, 0.93760847680602927660, 0.49652994760912213217, 14)] |
|||
[TestCase(0, 10.0, 5.0, 133.59268491136563443, -2651.8719709592179510, 14)] |
|||
[TestCase(0, -10.0, -5.0, 133.59268491136563443, -2651.8719709592179510, 14)] |
|||
[TestCase(0.001, 1.0, 1.0, 0.93752745412589393347, 0.49688729751353492126, 14)] |
|||
[TestCase(1, 10.0, 5.0, 183.59643950073504037, -2541.3865266173166975, 14)] |
|||
[TestCase(-1, 10.0, 5.0, 183.59643950073504037, -2541.3865266173166975, 14)] |
|||
[TestCase(1, 100.0, 100.0, 5.4401809306071216972E41, -7.1708320564209428560E41, 14)] |
|||
[TestCase(2, 4.0, 0.0, 6.4221893752841055416, 0.0, 14)] |
|||
[TestCase(2, 32.0, -64.0, 2.9566228902322313083E12, -2.1928744176387872855E12, 14)] |
|||
[TestCase(100, 1.0, 1.0, -9.5168076700036928783E-174, -4.7113294059604970982E-176, 14)] |
|||
[TestCase(100, -1.0, -1.0, -9.5168076700036928783E-174, -4.7113294059604970982E-176, 14)] |
|||
public void ComplexBesselInExact(double v, double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.BesselI(v, new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region Modified Bessel K
|
|||
|
|||
[Test] |
|||
public void BesselK0Approx([Range(0.20, 2.0, 0.20)] double x) |
|||
{ |
|||
// Approx by Abramowitz/Stegun 9.8.5
|
|||
Assert.AreEqual(Evaluate.Polynomial(x / 2.0, -Math.Log(x / 2.0) * SpecialFunctions.BesselI(0, x) - 0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740), SpecialFunctions.BesselK(0, x), 1e-8); |
|||
} |
|||
|
|||
[TestCase(1e-10, 23.14178244559887)] |
|||
[TestCase(1e-5, 11.62885698094436)] |
|||
[TestCase(0.005, 5.414288971329485)] |
|||
[TestCase(0.5, 0.9244190712276659)] |
|||
[TestCase(1.5, 0.2138055626475257)] |
|||
[TestCase(10.0, 0.00001778006231616765)] |
|||
[TestCase(100.0, 4.656628229175902e-45)] |
|||
public void BesselK0Exact(double x, double expected) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative(expected, SpecialFunctions.BesselK(0, x), 14); |
|||
} |
|||
|
|||
[Test] |
|||
public void BesselK1Approx([Range(0.20, 2.0, 0.20)] double x) |
|||
{ |
|||
// Approx by Abramowitz/Stegun 9.8.7
|
|||
Assert.AreEqual(Evaluate.Polynomial(x / 2.0, x * Math.Log(x / 2.0) * SpecialFunctions.BesselI(1, x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686), SpecialFunctions.BesselK(1, x) * x, 1e-8); |
|||
} |
|||
|
|||
[TestCase(1e-10, 1.0e+10)] |
|||
[TestCase(1e-5, 99999.99993935572)] |
|||
[TestCase(0.005, 199.9852143257300)] |
|||
[TestCase(0.5, 1.656441120003301)] |
|||
[TestCase(1.5, 0.2773878004568438)] |
|||
[TestCase(10.0, 0.00001864877345382558)] |
|||
[TestCase(100.0, 4.679853735636909e-45)] |
|||
public void BesselK1Exact(double x, double expected) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative(expected, SpecialFunctions.BesselK(1, x), 14); |
|||
} |
|||
|
|||
[TestCase(0, 0.0, 0.0, double.PositiveInfinity, 0.0, 14)] |
|||
[TestCase(0, 1.0, 0.0, 0.42102443824070833334, 0.0, 14)] |
|||
[TestCase(0, -1.0, 0.0, 0.42102443824070833334, -3.9774632605064226373, 14)] |
|||
[TestCase(0, 0.0, 1.0, -0.13863371520405399968, -1.2019697153172064991, 14)] |
|||
[TestCase(0, 0.0, -1.0, -0.13863371520405399968, 1.2019697153172064991, 14)] |
|||
[TestCase(0, 1.0, 1.0, 0.080197726946517818727, -0.35727745928533025061, 13)] |
|||
[TestCase(0, -1.0, -1.0, -1.4796971087496251933, 2.5883064433920073708, 14)] |
|||
[TestCase(0, 10.0, 5.0, 8.2975524594259100800E-6, 0.000014668532910602615041, 13)] |
|||
[TestCase(0, -10.0, -5.0, 8331.1015105237170946, 419.69381215941520621, 8)] |
|||
[TestCase(0.001, 1.0, 1.0, 0.080197683816936985794, -0.35727755473115169498, 14)] |
|||
[TestCase(1, 10.0, 5.0, 8.9074008179539983668E-6, 0.000015086785431163546150, 13)] |
|||
[TestCase(-1, 10.0, 5.0, 8.9074008179539983668E-6, 0.000015086785431163546150, 13)] |
|||
[TestCase(1, 100.0, 100.0, 3.8914899953028699000E-45, 5.3411641462767852385E-46, 14)] |
|||
[TestCase(-1, 100.0, 100.0, 3.8914899953028699000E-45, 5.3411641462767852385E-46, 14)] |
|||
[TestCase(2, 4.0, 0.0, 0.017401425529487240005, 0.0, 14)] |
|||
[TestCase(2, 32.0, -64.0, -3.2911230526496297606E-16, 1.8699562198185709108E-15, 14)] |
|||
[TestCase(70, -32.0, 0.0, 1.1846697842374187940E12, -1.7227016660469842981E-14, 14)] |
|||
[TestCase(70, 0.0, -1000.0, -0.019114301898634193541, -0.034775023193538639117, 13)] |
|||
[TestCase(70, -100.0, -100.0, 1.0618924726225417382E37, 1.3273296182243209319E36, 11)] |
|||
[TestCase(100, 1.0, 1.0, -5.2537311742300294807E170, 2.6534221390474409958E168, 10)] |
|||
public void ComplexBesselKnExact(double v, double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.BesselK(v, new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
#endregion
|
|||
} |
|||
} |
|||
@ -0,0 +1,61 @@ |
|||
using MathNet.Numerics.UnitTests; |
|||
using NUnit.Framework; |
|||
using System.Numerics; |
|||
|
|||
namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests |
|||
{ |
|||
/// <summary>
|
|||
/// Hankel functions tests.
|
|||
/// </summary>
|
|||
[TestFixture, Category("Functions")] |
|||
public class HankelTests |
|||
{ |
|||
[TestCase(0, 1.0, 0.0, 0.76519768655796655145, 0.088256964215676957983, 14)] |
|||
[TestCase(0, -1.0, 0.0, -0.76519768655796655145, 0.088256964215676957983, 14)] |
|||
[TestCase(0, 0.0, 1.0, 0.0, -0.26803248203398854876, 14)] |
|||
[TestCase(0, 0.0, -1.0, 2.5321317555040166712, -0.26803248203398854876, 14)] |
|||
[TestCase(0, 1.0, 1.0, 0.22744989480229475542, -0.051055458673089618135, 14)] |
|||
[TestCase(0, -1.0, -1.0, 2.1026668484143533086, -1.0441153538913338825, 14)] |
|||
[TestCase(0, 10.0, 5.0, -0.0014390626993822736024, 0.00069798313383186497809, 13)] |
|||
[TestCase(0, -10.0, -5.0, -35.580621321600125310, 0.40212121657624156694, 13)] |
|||
[TestCase(0.001, 1.0, 1.0, 0.22736947730907729941, -0.051412645636651723164, 14)] |
|||
[TestCase(1, 10.0, 5.0, 0.00065590035117675806376, 0.0014959445758611476093, 13)] |
|||
[TestCase(-1, 10.0, 5.0, -0.00065590035117675806376, -0.0014959445758611476093, 13)] |
|||
[TestCase(1, 100.0, 100.0, -2.4773994749804332257E-45, 3.4002907029806139579E-46, 14)] |
|||
[TestCase(2, 4.0, 0.0, 0.36412814585207280421, 0.21590359460361499453, 14)] |
|||
[TestCase(2, 32.0, -64.0, -5.3677688223644563726E26, -2.0457056534982487923E26, 14)] |
|||
[TestCase(100, 1.0, 1.0, -1.6892209981554826575E168, 3.3446291442187872047E170, 11)] |
|||
public void HankelH1Approx(double v, double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.HankelH1(v, new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
|
|||
[TestCase(0, 1.0, 0.0, 0.76519768655796655145, -0.088256964215676957983, 14)] |
|||
[TestCase(0, -1.0, 0.0, 2.2955930596738996543, -0.088256964215676957983, 14)] |
|||
[TestCase(0, 0.0, 1.0, 2.5321317555040166712, 0.26803248203398854876, 14)] |
|||
[TestCase(0, 0.0, -1.0, 0.0, 0.26803248203398854876, 14)] |
|||
[TestCase(0, 1.0, 1.0, 1.6477670588097637978, -0.94200443654515464620, 14)] |
|||
[TestCase(0, -1.0, -1.0, -0.22744989480229475542, 0.051055458673089618135, 14)] |
|||
[TestCase(0, 10.0, 5.0, -35.577743196201360763, 0.40072525030857783699, 13)] |
|||
[TestCase(0, -10.0, -5.0, 0.0014390626993822736024, -0.00069798313383186497809, 13)] |
|||
[TestCase(0.001, 1.0, 1.0, 1.6492441345285023430, -0.93941539521920296528, 14)] |
|||
[TestCase(1, 10.0, 5.0, -1.8435187716318316449, -34.874375165530875652, 13)] |
|||
[TestCase(-1, 10.0, 5.0, 1.8435187716318316449, 34.874375165530875652, 13)] |
|||
[TestCase(1, 100.0, 100.0, -1.4341664112841885712E42, 1.0880361861214243394E42, 14)] |
|||
[TestCase(2, 4.0, 0.0, 0.36412814585207280421, -0.21590359460361499453, 14)] |
|||
[TestCase(2, 32.0, -64.0, -1.1400282228790383437E-29, -1.0479217248690493285E-29, 14)] |
|||
[TestCase(100, 1.0, 1.0, 1.6892209981554826575E168, -3.3446291442187872047E170, 10)] |
|||
public void HankelH2Approx(double v, double zr, double zi, double cyr, double cyi, int decimalPlaces) |
|||
{ |
|||
AssertHelpers.AlmostEqualRelative( |
|||
new Complex(cyr, cyi), |
|||
SpecialFunctions.HankelH2(v, new Complex(zr, zi)), |
|||
decimalPlaces |
|||
); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,98 @@ |
|||
using System; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using System.Numerics; |
|||
using System.Text; |
|||
|
|||
namespace MathNet.Numerics |
|||
{ |
|||
/// <summary>
|
|||
/// This partial implementation of the SpecialFunctions class contains all methods related to the Airy functions.
|
|||
/// </summary>
|
|||
public static partial class SpecialFunctions |
|||
{ |
|||
/// <summary>
|
|||
/// Airy function Ai
|
|||
/// </summary>
|
|||
/// <param name="z">The value to compute the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex AiryAi(Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cairy(z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Airy function Ai
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double AiryAi(double x) |
|||
{ |
|||
return AiryAi(new Complex(x, 0)).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Derivative of the Airy function Ai
|
|||
/// </summary>
|
|||
/// <param name="z">The value to compute the derivative of the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex AiryAiPrime(Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.CairyPrime(z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Derivative of the Airy function Ai
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the derivative of the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double AiryAiPrime(double x) |
|||
{ |
|||
return AiryAiPrime(new Complex(x, 0)).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Airy function Bi(z)
|
|||
/// </summary>
|
|||
/// <param name="z">The value to compute the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex AiryBi(Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cbiry(z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Airy function Bi(x)
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double AiryBi(double x) |
|||
{ |
|||
return AiryBi(new Complex(x, 0)).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Derivative of the Airy function Bi(z)
|
|||
/// </summary>
|
|||
/// <param name="z">The value to compute the derivative of the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex AiryBiPrime(Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.CbiryPrime(z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Derivative of the Airy function Bi(z)
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the derivative of the Airy function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double AiryBiPrime(double x) |
|||
{ |
|||
return AiryBiPrime(new Complex(x, 0)).Real; |
|||
} |
|||
} |
|||
} |
|||
File diff suppressed because it is too large
@ -0,0 +1,664 @@ |
|||
using System; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using System.Numerics; |
|||
using System.Text; |
|||
|
|||
namespace MathNet.Numerics |
|||
{ |
|||
// References:
|
|||
// [1] https://github.com/scipy/scipy/blob/master/scipy/special/amos_wrappers.c
|
|||
public static partial class SpecialFunctions |
|||
{ |
|||
private class AmosWrapper |
|||
{ |
|||
#region AiryAi
|
|||
|
|||
public Complex Cairy(Complex z) |
|||
{ |
|||
int id = 0; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double air = double.NaN; |
|||
double aii = double.NaN; |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr); |
|||
return new Complex(air, aii); |
|||
} |
|||
|
|||
public Complex CairyPrime(Complex z) |
|||
{ |
|||
int id = 1; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double air = double.NaN; |
|||
double aii = double.NaN; |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr); |
|||
return new Complex(air, aii); |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region AiryBi
|
|||
|
|||
public Complex Cbiry(Complex z) |
|||
{ |
|||
int id = 0; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double bir = double.NaN; |
|||
double bii = double.NaN; |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbiry(z.Real, z.Imaginary, id, kode, ref bir, ref bii, ref nz, ref ierr); |
|||
return new Complex(bir, bii); |
|||
} |
|||
|
|||
public Complex CbiryPrime(Complex z) |
|||
{ |
|||
int id = 1; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double bipr = double.NaN; |
|||
double bipi = double.NaN; |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbiry(z.Real, z.Imaginary, id, kode, ref bipr, ref bipi, ref nz, ref ierr); |
|||
return new Complex(bipr, bipi); |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region BesselJ
|
|||
|
|||
public Complex Cbesj(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); |
|||
} |
|||
|
|||
int sign = 1; |
|||
if (v < 0) |
|||
{ |
|||
v = -v; |
|||
sign = -1; |
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
double[] cyjr = new double[n]; |
|||
double[] cyji = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyjr[i] = double.NaN; |
|||
cyji[i] = double.NaN; |
|||
} |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr); |
|||
Complex cyj = new Complex(cyjr[0], cyji[0]); |
|||
|
|||
if (ierr == 2) |
|||
{ |
|||
//overflow
|
|||
cyj = CbesjScaled(v, z); |
|||
cyj = new Complex(cyj.Real * double.PositiveInfinity, cyj.Imaginary * double.PositiveInfinity); |
|||
} |
|||
|
|||
if (sign == -1) |
|||
{ |
|||
if (!ReflectJY(ref cyj, v)) |
|||
{ |
|||
double[] cyyr = new double[n]; |
|||
double[] cyyi = new double[n]; |
|||
double[] cwrkr = new double[n]; |
|||
double[] cwrki = new double[n]; |
|||
|
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyyr[i] = double.NaN; |
|||
cyyi[i] = double.NaN; |
|||
cwrkr[i] = double.NaN; |
|||
cwrki[i] = double.NaN; |
|||
} |
|||
|
|||
amos.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cwrkr, cwrki, ref ierr); |
|||
Complex cyy = new Complex(cyyr[0], cyyi[0]); |
|||
|
|||
cyj = RotateJY(cyj, cyy, v); |
|||
} |
|||
} |
|||
return cyj; |
|||
} |
|||
|
|||
private Complex CbesjScaled(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); |
|||
} |
|||
|
|||
int sign = 1; |
|||
if (v < 0) |
|||
{ |
|||
v = -v; |
|||
sign = -1; |
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 2; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double[] cyjr = new double[n]; |
|||
double[] cyji = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyjr[i] = double.NaN; |
|||
cyji[i] = double.NaN; |
|||
} |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr); |
|||
Complex cyj = new Complex(cyjr[0], cyji[0]); |
|||
|
|||
if (sign == -1) |
|||
{ |
|||
if (!ReflectJY(ref cyj, v)) |
|||
{ |
|||
double[] cyyr = new double[n]; |
|||
double[] cyyi = new double[n]; |
|||
double[] cworkr = new double[n]; |
|||
double[] cworki = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyyr[i] = double.NaN; |
|||
cyyi[i] = double.NaN; |
|||
cworkr[i] = double.NaN; |
|||
cworki[i] = double.NaN; |
|||
} |
|||
|
|||
amos.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cworkr, cworki, ref ierr); |
|||
Complex cyy = new Complex(cyyr[0], cyyi[0]); |
|||
|
|||
cyj = RotateJY(cyj, cyy, v); |
|||
} |
|||
} |
|||
return cyj; |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region BesselY
|
|||
|
|||
public Complex Cbesy(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); |
|||
} |
|||
|
|||
int sign = 1; |
|||
if (v < 0) |
|||
{ |
|||
v = -v; |
|||
sign = -1; |
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
Complex cyy; |
|||
|
|||
var amos = new AmosHelper(); |
|||
if (z.Real == 0 && z.Imaginary == 0) |
|||
{ |
|||
//overflow
|
|||
cyy = new Complex(double.NegativeInfinity, 0); |
|||
} |
|||
else |
|||
{ |
|||
double[] cyyr = new double[n]; |
|||
double[] cyyi = new double[n]; |
|||
double[] cworkr = new double[n]; |
|||
double[] cworki = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyyr[i] = double.NaN; |
|||
cyyi[i] = double.NaN; |
|||
cworkr[i] = double.NaN; |
|||
cworki[i] = double.NaN; |
|||
} |
|||
|
|||
amos.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cworkr, cworki, ref ierr); |
|||
cyy = new Complex(cyyr[0], cyyi[0]); |
|||
|
|||
if (ierr == 2) |
|||
{ |
|||
if (z.Real >= 0 && z.Imaginary == 0) |
|||
{ |
|||
//overflow
|
|||
cyy = new Complex(double.NegativeInfinity, 0); |
|||
} |
|||
} |
|||
} |
|||
|
|||
if (sign == -1) |
|||
{ |
|||
if (!ReflectJY(ref cyy, v)) |
|||
{ |
|||
double[] cyjr = new double[n]; |
|||
double[] cyji = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyjr[i] = double.NaN; |
|||
cyji[i] = double.NaN; |
|||
} |
|||
|
|||
amos.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr); |
|||
Complex cyj = new Complex(cyjr[0], cyji[0]); |
|||
|
|||
cyy = RotateJY(cyy, cyj, -v); |
|||
} |
|||
} |
|||
return cyy; |
|||
} |
|||
|
|||
public double CbesyReal(double v, double x) |
|||
{ |
|||
if (x < 0.0) |
|||
{ |
|||
return double.NaN; |
|||
} |
|||
|
|||
Complex z = new Complex(x, 0.0); |
|||
return Cbesy(v, z).Real; |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region BesselI
|
|||
|
|||
public Complex Cbesi(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); |
|||
} |
|||
|
|||
int sign = 1; |
|||
if (v < 0) |
|||
{ |
|||
v = -v; |
|||
sign = -1; |
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double[] cyir = new double[n]; |
|||
double[] cyii = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyir[i] = double.NaN; |
|||
cyii[i] = double.NaN; |
|||
} |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbesi(z.Real, z.Imaginary, v, kode, n, cyir, cyii, ref nz, ref ierr); |
|||
Complex cyi = new Complex(cyir[0], cyii[0]); |
|||
|
|||
if (ierr == 2) |
|||
{ |
|||
//overflow
|
|||
if (z.Imaginary == 0 && (z.Real >= 0 || v == Math.Floor(v))) |
|||
{ |
|||
if (z.Real < 0 && v / 2 != Math.Floor(v / 2)) |
|||
cyi = new Complex(double.NegativeInfinity, 0); |
|||
else |
|||
cyi = new Complex(double.PositiveInfinity, 0); |
|||
} |
|||
else |
|||
{ |
|||
cyi = CbesiScaled(v * sign, z); |
|||
cyi = new Complex(cyi.Real * double.PositiveInfinity, cyi.Imaginary * double.PositiveInfinity); |
|||
} |
|||
} |
|||
|
|||
if (sign == -1) |
|||
{ |
|||
if (!ReflectI(cyi, v)) |
|||
{ |
|||
double[] cykr = new double[n]; |
|||
double[] cyki = new double[n]; |
|||
amos.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr); |
|||
Complex cyk = new Complex(cykr[0], cyki[0]); |
|||
|
|||
cyi = RotateI(cyi, cyk, v); |
|||
} |
|||
} |
|||
|
|||
return cyi; |
|||
} |
|||
|
|||
private Complex CbesiScaled(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); |
|||
} |
|||
|
|||
int sign = 1; |
|||
if (v < 0) |
|||
{ |
|||
v = -v; |
|||
sign = -1; |
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 2; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double[] cyir = new double[n]; |
|||
double[] cyii = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyir[i] = double.NaN; |
|||
cyii[i] = double.NaN; |
|||
} |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbesi(z.Real, z.Imaginary, v, kode, n, cyir, cyii, ref nz, ref ierr); |
|||
Complex cyi = new Complex(cyir[0], cyii[0]); |
|||
|
|||
if (sign == -1) |
|||
{ |
|||
if (!ReflectI(cyi, v)) |
|||
{ |
|||
double[] cykr = new double[n]; |
|||
double[] cyki = new double[n]; |
|||
amos.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr); |
|||
Complex cyk = new Complex(cykr[0], cyki[0]); |
|||
|
|||
//adjust scaling to match zbesi
|
|||
cyk = Rotate(cyk, -z.Imaginary / Math.PI); |
|||
if (z.Real > 0) |
|||
{ |
|||
cyk = new Complex(cyk.Real * Math.Exp(-2 * z.Real), cyk.Imaginary * Math.Exp(-2 * z.Real)); |
|||
} |
|||
//v -> -v
|
|||
cyi = RotateI(cyi, cyk, v); |
|||
} |
|||
} |
|||
|
|||
return cyi; |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region BesselK
|
|||
|
|||
public Complex Cbesk(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); |
|||
} |
|||
if (v < 0) |
|||
{ |
|||
//K_v == K_{-v} even for non-integer v
|
|||
v = -v; |
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double[] cykr = new double[n]; |
|||
double[] cyki = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cykr[i] = double.NaN; |
|||
cyki[i] = double.NaN; |
|||
} |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr); |
|||
Complex cyk = new Complex(cykr[0], cyki[0]); |
|||
|
|||
if (ierr == 1) |
|||
{ |
|||
if (z.Real == 0.0 && z.Imaginary == 0.0) |
|||
{ |
|||
cyk = new Complex(double.PositiveInfinity, 0); |
|||
} |
|||
} |
|||
else if (ierr == 2) |
|||
{ |
|||
if (z.Real >= 0 && z.Imaginary == 0) |
|||
{ |
|||
//overflow
|
|||
cyk = new Complex(double.PositiveInfinity, 0); |
|||
} |
|||
} |
|||
|
|||
return cyk; |
|||
} |
|||
|
|||
public double CbeskReal(double v, double z) |
|||
{ |
|||
if (z < 0) |
|||
{ |
|||
return double.NaN; |
|||
} |
|||
else if (z == 0) |
|||
{ |
|||
return double.PositiveInfinity; |
|||
} |
|||
else if (z > 710 * (1 + Math.Abs(v))) |
|||
{ |
|||
// Underflow. See uniform expansion https://dlmf.nist.gov/10.41
|
|||
// This condition is not a strict bound (it can underflow earlier),
|
|||
// rather, we are here working around a restriction in AMOS.
|
|||
|
|||
return 0; |
|||
} |
|||
else |
|||
{ |
|||
Complex w = new Complex(z, 0); |
|||
return Cbesk(v, w).Real; |
|||
} |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region HankelH1
|
|||
|
|||
public Complex Cbesh1(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); ; |
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 1; |
|||
int m = 1; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double[] cyhr = new double[n]; |
|||
double[] cyhi = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyhr[i] = double.NaN; |
|||
cyhi[i] = double.NaN; |
|||
} |
|||
|
|||
int sign = 1; |
|||
if (v < 0) |
|||
{ |
|||
v = -v; |
|||
sign = -1; |
|||
} |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr); |
|||
Complex cyh = new Complex(cyhr[0], cyhi[0]); |
|||
|
|||
if (sign == -1) |
|||
{ |
|||
cyh = Rotate(cyh, v); |
|||
} |
|||
|
|||
return cyh; |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region HankelH2
|
|||
|
|||
public Complex Cbesh2(double v, Complex z) |
|||
{ |
|||
if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); |
|||
} |
|||
|
|||
if (v == 0 && z.Real == 0 && z.Imaginary == 0) |
|||
{ |
|||
return new Complex(double.NaN, double.NaN); // ComplexInfinity
|
|||
} |
|||
|
|||
int n = 1; |
|||
int kode = 1; |
|||
int m = 2; |
|||
int nz = 0; |
|||
int ierr = 0; |
|||
|
|||
double[] cyhr = new double[n]; |
|||
double[] cyhi = new double[n]; |
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
cyhr[i] = double.NaN; |
|||
cyhi[i] = double.NaN; |
|||
} |
|||
|
|||
int sign = 1; |
|||
if (v < 0) |
|||
{ |
|||
v = -v; |
|||
sign = -1; |
|||
} |
|||
|
|||
var amos = new AmosHelper(); |
|||
amos.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr); |
|||
Complex cyh = new Complex(cyhr[0], cyhi[0]); |
|||
|
|||
if (sign == -1) |
|||
{ |
|||
cyh = Rotate(cyh, -v); |
|||
} |
|||
return cyh; |
|||
} |
|||
|
|||
#endregion
|
|||
|
|||
#region utilities
|
|||
|
|||
private double SinPi(double x) |
|||
{ |
|||
if (Math.Floor(x) == x && Math.Abs(x) < 1.0e14) |
|||
{ |
|||
//Return 0 when at exact zero, as long as the floating point number is
|
|||
//small enough to distinguish integer points from other points.
|
|||
|
|||
return 0; |
|||
} |
|||
return Math.Sin(Math.PI * x); |
|||
} |
|||
|
|||
private double CosPi(double x) |
|||
{ |
|||
if (Math.Floor(x + 0.5) == x + 0.5 && Math.Abs(x) < 1.0E14) |
|||
{ |
|||
//Return 0 when at exact zero, as long as the floating point number is
|
|||
//small enough to distinguish integer points from other points.
|
|||
|
|||
return 0; |
|||
} |
|||
return Math.Cos(Math.PI * x); |
|||
} |
|||
|
|||
private Complex Rotate(Complex z, double v) |
|||
{ |
|||
double c = CosPi(v); |
|||
double s = SinPi(v); |
|||
return new Complex(z.Real * c - z.Imaginary * s, z.Real * s + z.Imaginary * c); |
|||
} |
|||
|
|||
private Complex RotateJY(Complex j, Complex y, double v) |
|||
{ |
|||
double c = CosPi(v); |
|||
double s = SinPi(v); |
|||
return new Complex(j.Real * c - y.Real * s, j.Imaginary * c - y.Imaginary * s); |
|||
} |
|||
|
|||
private bool ReflectJY(ref Complex jy, double v) |
|||
{ |
|||
//NB: Y_v may be huge near negative integers -- so handle exact
|
|||
// integers carefully
|
|||
|
|||
if (v != Math.Floor(v)) |
|||
{ |
|||
return false; |
|||
} |
|||
|
|||
int i = (int)(v - 16384.0 * Math.Floor(v / 16384.0)); |
|||
if (i % 2 == 1) |
|||
{ |
|||
jy = new Complex(-jy.Real, -jy.Imaginary); |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
|
|||
private bool ReflectI(Complex ik, double v) |
|||
{ |
|||
if (v != Math.Floor(v)) |
|||
{ |
|||
return false; |
|||
} |
|||
|
|||
return true; //I is symmetric for integer v
|
|||
} |
|||
|
|||
private Complex RotateI(Complex i, Complex k, double v) |
|||
{ |
|||
double s = Math.Sin(v * Math.PI) * (2.0 / Math.PI); |
|||
return new Complex(i.Real + s * k.Real, i.Imaginary + s * k.Imaginary); |
|||
} |
|||
|
|||
#endregion
|
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,108 @@ |
|||
using System; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using System.Numerics; |
|||
using System.Text; |
|||
|
|||
namespace MathNet.Numerics |
|||
{ |
|||
/// <summary>
|
|||
/// This partial implementation of the SpecialFunctions class contains all methods related to the Bessel functions.
|
|||
/// </summary>
|
|||
public static partial class SpecialFunctions |
|||
{ |
|||
/// <summary>
|
|||
/// Bessel function of the first kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="z">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex BesselJ(double n, Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cbesj(n, z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Bessel function of the first kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="x">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double BesselJ(double n, double x) |
|||
{ |
|||
return BesselJ(n, new Complex(x, 0)).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Bessel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="z">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex BesselY(double n, Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cbesy(n, z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Bessel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="x">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double BesselY(double n, double x) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.CbesyReal(n, x); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Modified Bessel function of the first kind, of order n
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="z">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex BesselI(double n, Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cbesi(n, z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Modified Bessel function of the first kind, of order n
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="x">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double BesselI(double n, double x) |
|||
{ |
|||
return BesselI(n, new Complex(x, 0)).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Modified Bessel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="z">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex BesselK(double n, Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cbesk(n, z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Modified Bessel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="x">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double BesselK(double n, double x) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.CbeskReal(n, x); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,60 @@ |
|||
using System; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using System.Numerics; |
|||
using System.Text; |
|||
|
|||
namespace MathNet.Numerics |
|||
{ |
|||
/// <summary>
|
|||
/// This partial implementation of the SpecialFunctions class contains all methods related to the Hankel function.
|
|||
/// </summary>
|
|||
public static partial class SpecialFunctions |
|||
{ |
|||
/// <summary>
|
|||
/// Hankel function of the first kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="z">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex HankelH1(double n, Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cbesh1(n, z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Hankel function of the first kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="x">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double HankelH1(double n, double x) |
|||
{ |
|||
return HankelH1(n, new Complex(x, 0)).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Hankel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Hankel function</param>
|
|||
/// <param name="z">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex HankelH2(double n, Complex z) |
|||
{ |
|||
var amos = new AmosWrapper(); |
|||
return amos.Cbesh2(n, z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Hankel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the Bessel function</param>
|
|||
/// <param name="x">The value to compute the Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double HankelH2(double n, double x) |
|||
{ |
|||
return HankelH2(n, new Complex(x, 0)).Real; |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,66 @@ |
|||
using System; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using System.Numerics; |
|||
using System.Text; |
|||
|
|||
namespace MathNet.Numerics |
|||
{ |
|||
/// <summary>
|
|||
/// This partial implementation of the SpecialFunctions class contains all methods related to the spherical Bessel functions.
|
|||
/// </summary>
|
|||
public static partial class SpecialFunctions |
|||
{ |
|||
/// <summary>
|
|||
/// Spherical Bessel function of the first kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the spherical Bessel function</param>
|
|||
/// <param name="z">The value to compute the spherical Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex SphericalBesselJ(double n, Complex z) |
|||
{ |
|||
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
|
|||
|
|||
return rthpi * BesselJ(n + 0.5, z) / Complex.Sqrt(z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Spherical Bessel function of the first kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the spherical Bessel function</param>
|
|||
/// <param name="x">The value to compute the spherical Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double SphericalBesselJ(double n, double x) |
|||
{ |
|||
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
|
|||
|
|||
return rthpi * BesselJ(n + 0.5, x) / Math.Sqrt(x); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Spherical Bessel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the spherical Bessel function</param>
|
|||
/// <param name="z">The value to compute the spherical Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static Complex SphericalBesselY(double n, Complex z) |
|||
{ |
|||
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
|
|||
|
|||
return rthpi * BesselY(n + 0.5, z) / Complex.Sqrt(z); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Spherical Bessel function of the second kind
|
|||
/// </summary>
|
|||
/// <param name="n">The order of the spherical Bessel function</param>
|
|||
/// <param name="x">The value to compute the spherical Bessel function of.</param>
|
|||
/// <returns></returns>
|
|||
public static double SphericalBesselY(double n, double x) |
|||
{ |
|||
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
|
|||
|
|||
return rthpi * BesselY(n + 0.5, x) / Math.Sqrt(x); |
|||
} |
|||
} |
|||
} |
|||
Loading…
Reference in new issue