tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Merge pull request #601 from diluculo/AiryBesselHankel 

Airy, Bessel, and Hankel functions
pull/587/merge
Christoph Ruegg 8 years ago
committed by GitHub
parent
b625985d48 03ef3b8a3a
commit
2be8c3a6ee
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 9575 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 92
      src/Numerics.Tests/SpecialFunctionsTests/AiryTests.cs
  2. 276
      src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
  3. 61
      src/Numerics.Tests/SpecialFunctionsTests/HankelTests.cs
  4. 138
      src/Numerics/SpecialFunctions/Airy.cs
  5. 7676
      src/Numerics/SpecialFunctions/Amos/AmosHelper.cs
  6. 1041
      src/Numerics/SpecialFunctions/Amos/AmosWrapper.cs
  7. 141
      src/Numerics/SpecialFunctions/Bessel.cs
  8. 72
      src/Numerics/SpecialFunctions/Hankel.cs
  9. 78
      src/Numerics/SpecialFunctions/SphericalBessel.cs

92
src/Numerics.Tests/SpecialFunctionsTests/AiryTests.cs
View File

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

276
src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs
View File

@ -0,0 +1,276 @@
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
#region Ratio of Bessel I : I(n + 1, z) / I(n, z)
[TestCase(0, 1, 1, 0.57495795977223997079, 0.35054769385125934266, 14)]
[TestCase(0, 10, 10, 0.97503660846868673171, 0.025655591609138262390, 14)]
[TestCase(0, 100, 100, 0.99750003174335729017, 0.0025062812447888833561, 12)]
[TestCase(0, 1E6, 1E6, 0.99999975000000000003, 2.5000006250003125000E-7, 8)] // cf. BesselI(0, 1e6 + 1e6 j) = 7.4E434290 - 6.9E434290 j and BesselI(1, 1e6 + 1e6 j) = 7.4E434290 - 6.9E434290 j
public void BesselIRatioExact(int n, double zr, double zi, double cyr, double cyi, int decimalPlaces)
{
var z = new Complex(zr, zi);
var actual = SpecialFunctions.BesselI(n + 1, z, true) / SpecialFunctions.BesselI(n, z, true);
AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces);
}
#endregion
#region Ratio of Bessel K : K(n + 1, z) / K(n, z)
[TestCase(0, 1E-6, 1E-6, 38803.844378426554737, -34562.243460439510287, 14)]
[TestCase(0, 1, 1, 1.2397012468882428039, -0.20950920530836317445, 14)]
[TestCase(0, 10, 10, 1.0249731393042359657, -0.024405853995209668702, 14)]
[TestCase(0, 100, 100, 1.0024999692332050665, -0.0024937812450508461147, 12)]
[TestCase(0, 1E6, 1E6, 1.0000002500000000000, -2.4999993750003125000E-7, 9)]
public void BesselKRatioExact(int n, double zr, double zi, double cyr, double cyi, int decimalPlaces)
{
var z = new Complex(zr, zi);
var actual = SpecialFunctions.BesselK(n + 1, z, true) / SpecialFunctions.BesselK(n, z, true);
AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces);
}
#endregion
}
}

61
src/Numerics.Tests/SpecialFunctionsTests/HankelTests.cs
View File

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

138
src/Numerics/SpecialFunctions/Airy.cs
View File

@ -0,0 +1,138 @@
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(z).
/// <p/>
/// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z).
/// </summary>
/// <param name="z">The value to compute the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static Complex AiryAi(Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCairy(z) : amos.Cairy(z);
}
/// <summary>
/// Airy function Ai(z).
/// <p/>
/// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z).
/// </summary>
/// <param name="z">The value to compute the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static double AiryAi(double z, bool expScaled = false)
{
if (expScaled)
{
var amos = new AmosWrapper();
return amos.ScaledCairy(z);
}
else
{
return AiryAi(new Complex(z, 0), expScaled).Real;
}
}
/// <summary>
/// Derivative of the Airy function Ai.
/// <p/>
/// If expScaled is true, returns Exp(zta) * d/dz Ai(z), where zta = (2/3) * z * Sqrt(z).
/// </summary>
/// <param name="z">The value to compute the derivative of the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static Complex AiryAiPrime(Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCairyPrime(z) : amos.CairyPrime(z);
}
/// <summary>
/// Derivative of the Airy function Ai.
/// <p/>
/// If expScaled is true, returns Exp(zta) * d/dz Ai(z), where zta = (2/3) * z * Sqrt(z).
/// </summary>
/// <param name="z">The value to compute the derivative of the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static double AiryAiPrime(double z, bool expScaled = false)
{
if (expScaled)
{
var amos = new AmosWrapper();
return amos.ScaledCairyPrime(z);
}
else
{
return AiryAiPrime(new Complex(z, 0), expScaled).Real;
}
}
/// <summary>
/// Airy function Bi(z).
/// <p/>
/// If expScaled is true, returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
/// </summary>
/// <param name="z">The value to compute the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static Complex AiryBi(Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbiry(z) : amos.Cbiry(z);
}
/// <summary>
/// Airy function Bi(x).
/// <p/>
/// If expScaled is true, returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
/// </summary>
/// <param name="z">The value to compute the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static double AiryBi(double z, bool expScaled = false)
{
return AiryBi(new Complex(z, 0), expScaled).Real;
}
/// <summary>
/// Derivative of the Airy function Bi(z).
/// <p/>
/// If expScaled is true, returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
/// </summary>
/// <param name="z">The value to compute the derivative of the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static Complex AiryBiPrime(Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbiryPrime(z) : amos.CbiryPrime(z);
}
/// <summary>
/// Derivative of the Airy function Bi(z).
/// <p/>
/// If expScaled is true, returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real).
/// </summary>
/// <param name="z">The value to compute the derivative of the Airy function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Airy function</param>
/// <returns></returns>
public static double AiryBiPrime(double z, bool expScaled = false)
{
return AiryBiPrime(new Complex(z, 0), expScaled).Real;
}
}
}

7676
src/Numerics/SpecialFunctions/Amos/AmosHelper.cs
View File

File diff suppressed because it is too large

1041
src/Numerics/SpecialFunctions/Amos/AmosWrapper.cs
View File

File diff suppressed because it is too large

141
src/Numerics/SpecialFunctions/Bessel.cs
View File

@ -0,0 +1,141 @@
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, J(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * J(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static Complex BesselJ(double v, Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesj(v, z) : amos.Cbesj(v, z);
}
/// <summary>
/// Bessel function of the first kind, J(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * J(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static double BesselJ(double v, double z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesj(v, z) : amos.Cbesj(v, z);
}
/// <summary>
/// Bessel function of the second kind, Y(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static Complex BesselY(double v, Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesy(v, z) : amos.Cbesy(v, z);
}
/// <summary>
/// Bessel function of the second kind, Y(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static double BesselY(double v, double z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesy(v, z) : amos.Cbesy(v, z);
}
/// <summary>
/// Modified Bessel function of the first kind, I(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(x)) * I(v, z) where x = z.Real.
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static Complex BesselI(double v, Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesi(v, z) : amos.Cbesi(v, z);
}
/// <summary>
/// Modified Bessel function of the first kind, I(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(x)) * I(v, z) where x = z.Real.
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static double BesselI(double v, double z, bool expScaled = false)
{
if (expScaled)
{
var amos = new AmosWrapper();
return amos.ScaledCbesi(v, z);
}
else
{
return BesselI(v, new Complex(z, 0), expScaled).Real;
}
}
/// <summary>
/// Modified Bessel function of the second kind, K(v, z).
/// <p/>
/// If expScaled is true, returns Exp(z) * K(v, z).
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static Complex BesselK(double v, Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesk(v, z) : amos.Cbesk(v, z);
}
/// <summary>
/// Modified Bessel function of the second kind, K(v, z).
/// <p/>
/// If expScaled is true, returns Exp(z) * K(v, z).
/// </summary>
/// <param name="v">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Bessel function</param>
/// <returns></returns>
public static double BesselK(double v, double z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesk(v, z) : amos.Cbesk(v, z);
}
}
}

72
src/Numerics/SpecialFunctions/Hankel.cs
View File

@ -0,0 +1,72 @@
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, H1(n, z).
/// <p/>
/// If expScaled is true, returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1).
/// </summary>
/// <param name="n">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
/// <returns></returns>
public static Complex HankelH1(double n, Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesh1(n, z) : amos.Cbesh1(n, z);
}
/// <summary>
/// Hankel function of the first kind, H1(n, z).
/// <p/>
/// If expScaled is true, returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1).
/// </summary>
/// <param name="n">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
/// <returns></returns>
public static double HankelH1(double n, double z, bool expScaled = false)
{
return HankelH1(n, new Complex(z, 0), expScaled).Real;
}
/// <summary>
/// Hankel function of the second kind, H2(n, z).
/// <p/>
/// If expScaled is true, returns Exp(z * j) * H2(n, z) where j = Sqrt(-1).
/// </summary>
/// <param name="n">The order of the Hankel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
/// <returns></returns>
public static Complex HankelH2(double n, Complex z, bool expScaled = false)
{
var amos = new AmosWrapper();
return (expScaled) ? amos.ScaledCbesh2(n, z) :amos.Cbesh2(n, z);
}
/// <summary>
/// Hankel function of the second kind, H2(n, z).
/// <p/>
/// If expScaled is true, returns Exp(z * j) * H2(n, z) where j = Sqrt(-1).
/// </summary>
/// <param name="n">The order of the Bessel function</param>
/// <param name="z">The value to compute the Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled Hankel function</param>
/// <returns></returns>
public static double HankelH2(double n, double z, bool expScaled = false)
{
return HankelH2(n, new Complex(z, 0), expScaled).Real;
}
}
}

78
src/Numerics/SpecialFunctions/SphericalBessel.cs
View File

@ -0,0 +1,78 @@
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, j(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * j(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the spherical Bessel function</param>
/// <param name="z">The value to compute the spherical Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
/// <returns></returns>
public static Complex SphericalBesselJ(double v, Complex z, bool expScaled = false)
{
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
return rthpi * BesselJ(v + 0.5, z, expScaled) / Complex.Sqrt(z);
}
/// <summary>
/// Spherical Bessel function of the first kind, j(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * j(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the spherical Bessel function</param>
/// <param name="z">The value to compute the spherical Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
/// <returns></returns>
public static double SphericalBesselJ(double v, double z, bool expScaled = false)
{
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
return rthpi * BesselJ(v + 0.5, z, expScaled) / Math.Sqrt(z);
}
/// <summary>
/// Spherical Bessel function of the second kind, y(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * y(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the spherical Bessel function</param>
/// <param name="z">The value to compute the spherical Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
/// <returns></returns>
public static Complex SphericalBesselY(double v, Complex z, bool expScaled = false)
{
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
return rthpi * BesselY(v + 0.5, z, expScaled) / Complex.Sqrt(z);
}
/// <summary>
/// Spherical Bessel function of the second kind, y(v, z).
/// <p/>
/// If expScaled is true, returns Exp(-Abs(y)) * y(v, z) where y = z.Imaginary.
/// </summary>
/// <param name="v">The order of the spherical Bessel function</param>
/// <param name="z">The value to compute the spherical Bessel function of.</param>
/// <param name="expScaled">If true, returns exponentially-scaled spherical Bessel function</param>
/// <returns></returns>
public static double SphericalBesselY(double v, double z, bool expScaled = false)
{
const double rthpi = 1.2533141373155002512; //sqrt(pi/2)
return rthpi * BesselY(v + 0.5, z, expScaled) / Math.Sqrt(z);
}
}
}
Write Preview
Loading…
Cancel
Save