diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..412eeda --- /dev/null +++ b/.gitattributes @@ -0,0 +1,22 @@ +# Auto detect text files and perform LF normalization +* text=auto + +# Custom for Visual Studio +*.cs diff=csharp +*.sln merge=union +*.csproj merge=union +*.vbproj merge=union +*.fsproj merge=union +*.dbproj merge=union + +# Standard to msysgit +*.doc diff=astextplain +*.DOC diff=astextplain +*.docx diff=astextplain +*.DOCX diff=astextplain +*.dot diff=astextplain +*.DOT diff=astextplain +*.pdf diff=astextplain +*.PDF diff=astextplain +*.rtf diff=astextplain +*.RTF diff=astextplain diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..05651d4 --- /dev/null +++ b/.gitignore @@ -0,0 +1,74 @@ +############ +## Windows +############ + +# Windows image file caches +Thumbs.db +*.o + +# Folder config file +Desktop.ini + +################# +## Visual Studio +################# + +## Ignore Visual Studio temporary files, build results, and +## files generated by popular Visual Studio add-ons. +# User-specific files +**/Debug +**/Release +win32/packages/ +**/.vs +*.tlog +*.obj +*.log +*.filters +*.lastbuildstate +*.manifest +*.cache +*.users +*.user +*.res +*.opendb +*.db +*.unsuccessfulbuild +*.ipch +*.pdb +*.exp +*.ilk +*.idb +*.opensdf +*.sdf +*.u2d +*.suo +*.tiff +*.txt +*.aps + +################ +# Linux and OSX +############### +\.DS* +*.so +*.d + + +# Created by https://www.gitignore.io + +### Xcode ### +build/ +*.pbxuser +!default.pbxuser +*.mode1v3 +!default.mode1v3 +*.mode2v3 +!default.mode2v3 +*.perspectivev3 +!default.perspectivev3 +xcuserdata +*.xccheckout +*.moved-aside +DerivedData +*.xcuserstate + diff --git a/LFMF_Examples.csv b/LFMF_Examples.csv new file mode 100644 index 0000000..0d036f0 --- /dev/null +++ b/LFMF_Examples.csv @@ -0,0 +1,6 @@ +h_tx__meter,h_rx__meter,f__mhz,P_tx__watt,N_s,d__km,epsilon,sigma,pol,rtn,A_btl__db,E_dBuVm,P_rx__dbm,method +0,0,0.01,1000,301,1000,15,0.005,0,0,184.5,-82.5,-114.9,1 +0,0,1,5000,301,5000,15,0.005,1,0,536.5,-387.5,-459.9,1 +5.5,1.5,10,500,315,15,15,0.005,0,0,151.3,7.6,-84.8,0 +1,1,0.45,1000,315,3000,15,0.005,1,0,264.3,-129.3,-194.8,1 +10,10,30,10000,315,5000,15,0.005,0,0,1574.9,-1393.3,-1495.3,1 diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..12746a1 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,14 @@ +SOFTWARE DISCLAIMER / RELEASE + +This software was developed by employees of the National Telecommunications and Information Administration (NTIA), an agency of the Federal Government and is provided to you as a public service. Pursuant to Title 15 United States Code Section 105, works of NTIA employees are not subject to copyright protection within the United States. + +The software is provided by NTIA “AS IS.” NTIA MAKES NO WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR STATUTORY, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NON-INFRINGEMENT AND DATA ACCURACY. NTIA does not warrant or make any representations regarding the use of the software or the results thereof, including but not limited to the correctness, accuracy, reliability or usefulness of the software. + +To the extent that NTIA holds rights in countries other than the United States, you are hereby granted the non-exclusive irrevocable and unconditional right to print, publish, prepare derivative works and distribute the NTIA software, in any medium, or authorize others to do so on your behalf, on a royalty-free basis throughout the World. + +You may improve, modify, and create derivative works of the software or any portion of the software, and you may copy and distribute such modifications or works. Modified works should carry a notice stating that you changed the software and should note the date and nature of any such change. + +You are solely responsible for determining the appropriateness of using and distributing the software and you assume all risks associated with its use, including but not limited to the risks and costs of program errors, compliance with applicable laws, damage to or loss of data, programs or equipment, and the unavailability or interruption of operation. This software is not intended to be used in any situation where a failure could cause risk of injury or damage to property. + +Please provide appropriate acknowledgments of NTIA’s creation of the software in any copies or derivative works of this software. + diff --git a/README.md b/README.md index 56c3c73..bcd516a 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,77 @@ -# LFMF -Low Frequency / Medium Frequency (LF/MF) Propagation Model +# Low Frequency / Medium Frequency (LF/MF) Propagation Model + +This code repository contains an the NTIA/ITS implementation of the Low Frequency / Medium Frequency (LF/MF) Propagation Model. LF/MF predicts basic transmission loss in the frequency range 0.01 - 30 MHz for propagation paths over a smooth Earth and antenna heights less than 50 meters. + +## Inputs ## + +| Variable | Type | Units | Limits | Description | +|-------------------|--------|-------|--------------|--------------| +| `h_tx__meter` | double | meter | 0 <= `h_tx__meter` <= 50 | Height of the transmitter | +| `h_rx__meter` | double | meter | 0 <= `h_rx__meter` <= 50 | Height of the receiver | +| `f__mhz` | double | MHz | 0.01 <= `f__mhz` <= 30 | Frequency | +| `P_tx__watt` | double | Watt | 0 < `P_tx__watt` | Transmitter power | +| `N_s` | double | N-Units | 250 <= `N_s` <= 400 | Surface refractivity | +| `d__km` | double | km | 0 < `d__km` | Path distance | +| `epsilon` | double | | 1 <= `epsilon` | Relative permittivity | +| `sigma` | double | S/m | 0 < `sigma` | Conductivity | +| `pol` | int | | | Polarization | + +## Outputs ## + +Outputs to LFMF are contained within a defined `Result` structure. + +| Variable | Type | Units | Description | +|---------------|--------|-------|-------------| +| `A_btl__db` | double | dB | Basic transmission loss | +| `E_dBuVm` | double | dB(uV/m) | Electrice field strength | +| `P_rx__dbm` | double | dBm | Received power | +| `method` | int | | Solution method | + +## Return Codes ## + +Possible return codes, including the corresponding defined constant name as defined in `LFMF.h`: + +| Value | Const Name | Description | +| ------|----------------------------------|--------------| +| 0 | `SUCCESS` | Successful execution | +| 1000 | `ERROR__TX_TERMINAL_HEIGHT` | TX terminal height is out of range | +| 1001 | `ERROR__RX_TERMINAL_HEIGHT` | RX terminal height is out of range | +| 1002 | `ERROR__FREQUENCY` | Frequency is out of range | +| 1003 | `ERROR__TX_POWER` | Transmit power is out of range | +| 1004 | `ERROR__SURFACE_REFRACTIVITY` | Surface refractivity is out of range | +| 1005 | `ERROR__PATH_DISTANCE` | Path distance is out of range | +| 1006 | `ERROR__EPSILON` | Epsilon is out of range | +| 1007 | `ERROR__SIGMA` | Sigma is out of range | +| 1008 | `ERROR__POLARIZATION` | Invalid value for polarization | + +## Example Values ## + +A set of example inputs and outputs are provided for testing purposes. This is not a comprehensive validation test set. The test set can be found in [LFMF_Examples.csv](LFMF_Examples.csv). + +## Notes on Code Style ## + +* In general, variables follow the naming convention in which a single underscore denotes a subscript (pseudo-LaTeX format), where a double underscore is followed by the units, i.e. h_tx__meter. +* Variables are named to match their corresponding mathematical variables from their publication text. +* Wherever possible, equation numbers are provided. + +## Configure and Build ## + +### C++ Software + +The software is designed to be built into a DLL (or corresponding library for non-Windows systems). The source code can be built for any OS that supports the standard C++ libraries. A Visual Studio 2019 project file is provided for Windows users to support the build process and configuration. + +### C#/.NET Wrapper Software + +The .NET support of LFMF consists of a simple pass-through wrapper around the native DLL. It is compiled to target .NET Framework 4.7.2. Distribution and updates are provided through the published NuGet package. + +## References ## + +* Bremmer, H. "Terrestrial Radio Waves" _Elsevier_, 1949. +* DeMinco, N. "Medium Frequency Propagation Prediction Techniques and Antenna Modeling for Intelligent Transportation Systems (ITS) Broadcast Applications", [_NTIA Report 99-368_](https://www.its.bldrdoc.gov/publications/2399.aspx), August 1999 +* DeMinco, N. "Ground-wave Analysis Model For MF Broadcast System", [_NTIA Report 86-203_](https://www.its.bldrdoc.gov/publications/2226.aspx), September 1986 +* Sommerfeld, A. "The propagation of waves in wireless telegraphy", _Ann. Phys._, 1909, 28, p.665 +* Wait, J. "Radiation From a Vertical Antenna Over a Curved Stratified Ground", _Journal of Research of the National Bureau of Standards_. Vol 56, No. 4, April 1956. Research Paper 2671 + +## Contact ## + +For questions, contact Nick DeMinco, (303) 497-3453, NDeminco@ntia.gov \ No newline at end of file diff --git a/dotnet/ITS.Propagation.LFMF/ITS.Propagation.LFMF.csproj b/dotnet/ITS.Propagation.LFMF/ITS.Propagation.LFMF.csproj new file mode 100644 index 0000000..d12cec5 --- /dev/null +++ b/dotnet/ITS.Propagation.LFMF/ITS.Propagation.LFMF.csproj @@ -0,0 +1,51 @@ + + + + + Debug + AnyCPU + {1C761232-C4C0-4403-AA3A-7160DCDF48C8} + Library + Properties + ITS.Propagation + ITS.Propagation.LFMF + v4.7.2 + 512 + true + + + true + full + false + bin\Debug\ + DEBUG;TRACE + prompt + 4 + bin\Debug\ITS.Propagation.LFMF.xml + + + pdbonly + true + bin\Release\ + TRACE + prompt + 4 + bin\Release\ITS.Propagation.LFMF.xml + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/dotnet/ITS.Propagation.LFMF/ITS.Propagation.LFMF.sln b/dotnet/ITS.Propagation.LFMF/ITS.Propagation.LFMF.sln new file mode 100644 index 0000000..46e729b --- /dev/null +++ b/dotnet/ITS.Propagation.LFMF/ITS.Propagation.LFMF.sln @@ -0,0 +1,25 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +# Visual Studio Version 16 +VisualStudioVersion = 16.0.30320.27 +MinimumVisualStudioVersion = 10.0.40219.1 +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "ITS.Propagation.LFMF", "ITS.Propagation.LFMF.csproj", "{1C761232-C4C0-4403-AA3A-7160DCDF48C8}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Any CPU = Debug|Any CPU + Release|Any CPU = Release|Any CPU + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {1C761232-C4C0-4403-AA3A-7160DCDF48C8}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {1C761232-C4C0-4403-AA3A-7160DCDF48C8}.Debug|Any CPU.Build.0 = Debug|Any CPU + {1C761232-C4C0-4403-AA3A-7160DCDF48C8}.Release|Any CPU.ActiveCfg = Release|Any CPU + {1C761232-C4C0-4403-AA3A-7160DCDF48C8}.Release|Any CPU.Build.0 = Release|Any CPU + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection + GlobalSection(ExtensibilityGlobals) = postSolution + SolutionGuid = {ECF2AF6A-510E-4A3F-BAA6-51FDC7F311C6} + EndGlobalSection +EndGlobal diff --git a/dotnet/ITS.Propagation.LFMF/LFMF.cs b/dotnet/ITS.Propagation.LFMF/LFMF.cs new file mode 100644 index 0000000..e72c476 --- /dev/null +++ b/dotnet/ITS.Propagation.LFMF/LFMF.cs @@ -0,0 +1,82 @@ +using System; +using System.Runtime.InteropServices; + +namespace ITS.Propagation +{ + /// + /// The Low Frequency / Medium Frequency (LF/MF) Propagation Model + /// + public static partial class LFMF + { + #region 32-Bit P/Invoke Definitions + + [DllImport("LFMF_x86.dll", CallingConvention = CallingConvention.StdCall, CharSet = CharSet.Ansi, EntryPoint = "LFMF")] + private static extern int LFMF_x86(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol, out Result result); + + #endregion + + #region 64-Bit P/Invoke Definitions + + [DllImport("LFMF_x64.dll", CallingConvention = CallingConvention.StdCall, CharSet = CharSet.Ansi, EntryPoint = "LFMF")] + private static extern int LFMF_x64(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol, out Result result); + + #endregion + + private delegate int LFMF_Delegate(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol, out Result result); + + private static LFMF_Delegate LFMF_Invoke; + + static LFMF() + { + if (Environment.Is64BitProcess) + LFMF_Invoke = LFMF_x64; + else + LFMF_Invoke = LFMF_x86; + } + + /// + /// Compute the LFMF propagation prediction + /// + /// Transmitter height, in meters + /// Receiver height, in meters + /// Frequency, in MHz + /// Transmit power, in Watts + /// Surface refractivity, in N-Units + /// Path distance, in km + /// Relative permittivity + /// Conductivity + /// Polarization + /// Prediction result + /// Error code + public static int Invoke(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, Polarization pol, out Result result) + { + return LFMF_Invoke(h_tx__meter, h_rx__meter, f__mhz, P_tx__watt, N_s, + d__km, epsilon, sigma, (int)pol, out result); + } + + /// + /// Compute the LFMF propagation prediction + /// + /// Transmitter height, in meters + /// Receiver height, in meters + /// Frequency, in MHz + /// Transmit power, in Watts + /// Surface refractivity, in N-Units + /// Path distance, in km + /// Relative permittivity + /// Conductivity + /// Polarization + /// Prediction result + /// Error code + public static int Invoke(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol, out Result result) + { + return LFMF_Invoke(h_tx__meter, h_rx__meter, f__mhz, P_tx__watt, N_s, + d__km, epsilon, sigma, pol, out result); + } + } +} diff --git a/dotnet/ITS.Propagation.LFMF/Properties/AssemblyInfo.cs b/dotnet/ITS.Propagation.LFMF/Properties/AssemblyInfo.cs new file mode 100644 index 0000000..57259c0 --- /dev/null +++ b/dotnet/ITS.Propagation.LFMF/Properties/AssemblyInfo.cs @@ -0,0 +1,36 @@ +using System.Reflection; +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; + +// General Information about an assembly is controlled through the following +// set of attributes. Change these attribute values to modify the information +// associated with an assembly. +[assembly: AssemblyTitle("ITS.Propagation.LFMF")] +[assembly: AssemblyDescription("")] +[assembly: AssemblyConfiguration("")] +[assembly: AssemblyCompany("")] +[assembly: AssemblyProduct("ITS.Propagation.LFMF")] +[assembly: AssemblyCopyright("")] +[assembly: AssemblyTrademark("")] +[assembly: AssemblyCulture("")] + +// Setting ComVisible to false makes the types in this assembly not visible +// to COM components. If you need to access a type in this assembly from +// COM, set the ComVisible attribute to true on that type. +[assembly: ComVisible(false)] + +// The following GUID is for the ID of the typelib if this project is exposed to COM +[assembly: Guid("1c761232-c4c0-4403-aa3a-7160dcdf48c8")] + +// Version information for an assembly consists of the following four values: +// +// Major Version +// Minor Version +// Build Number +// Revision +// +// You can specify all the values or you can default the Build and Revision Numbers +// by using the '*' as shown below: +// [assembly: AssemblyVersion("1.0.*")] +[assembly: AssemblyVersion("1.0.0.0")] +[assembly: AssemblyFileVersion("1.0.0.0")] diff --git a/dotnet/ITS.Propagation.LFMF/Result.cs b/dotnet/ITS.Propagation.LFMF/Result.cs new file mode 100644 index 0000000..c083281 --- /dev/null +++ b/dotnet/ITS.Propagation.LFMF/Result.cs @@ -0,0 +1,67 @@ +using System; +using System.Runtime.InteropServices; + +namespace ITS.Propagation +{ + public static partial class LFMF + { + /// + /// Polarization + /// + public enum Polarization : int + { + /// + /// Horizontal polarization + /// + Horizontal = 0, + + /// + /// Vertical polarization + /// + Vertical = 1 + } + + /// + /// Solution method used + /// + public enum Method : int + { + /// + /// Flat earth with curve correction + /// + FlatEarthCurveCorrection = 0, + + /// + /// Residue series + /// + ResidueSeries = 1 + } + + /// + /// LF/MF prediction results + /// + [StructLayout(LayoutKind.Sequential)] + public struct Result + { + /// + /// Basic transmission loss, in dB + /// + public double A_btl__db; + + /// + /// Electic field strength, in db(uV/m) + /// + public double E_dBuVm; + + /// + /// Received power, in dBm + /// + public double P_rx__dbm; + + /// + /// Solution method used + /// + public Method Method; + } + } +} diff --git a/dotnet/nuget/LFMF.nuspec b/dotnet/nuget/LFMF.nuspec new file mode 100644 index 0000000..99f0d27 --- /dev/null +++ b/dotnet/nuget/LFMF.nuspec @@ -0,0 +1,30 @@ + + + + LFMF + 1.0.0 + The Institute for Telecommunication Sciences + The Institute for Telecommunication Sciences + LICENSE.md + https://github.com/NTIA/LFMF + images\itslogo.png + true + Low Frequency / Medium Frequency (LF/MF) Propagation Model + None + Propagation + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/dotnet/nuget/build/net472/LFMF.targets b/dotnet/nuget/build/net472/LFMF.targets new file mode 100644 index 0000000..d4c1b1d --- /dev/null +++ b/dotnet/nuget/build/net472/LFMF.targets @@ -0,0 +1,16 @@ + + + + + PreserveNewest + LFMF_x86.dll + + + + + + PreserveNewest + LFMF_x64.dll + + + \ No newline at end of file diff --git a/dotnet/nuget/images/itslogo.png b/dotnet/nuget/images/itslogo.png new file mode 100644 index 0000000..3d2f319 Binary files /dev/null and b/dotnet/nuget/images/itslogo.png differ diff --git a/include/LFMF.h b/include/LFMF.h new file mode 100644 index 0000000..f86066a --- /dev/null +++ b/include/LFMF.h @@ -0,0 +1,89 @@ + +#include + +// Export the DLL functions as "C" and not C++ +#define DLLEXPORT extern "C" __declspec(dllexport) +#define MAX(x, y) (((x) > (y)) ? (x) : (y)) +#define MIN(x, y) (((x) < (y)) ? (x) : (y)) +#define DIM(x, y) (((x) > (y)) ? (x - y) : (0)) + +using std::complex; + +#define PI 3.1415926535897932384 +#define epsilon_0 8.854187817e-12 // Vacuum permittivity (F/m) +#define a_0__km 6370 // Earth radius, in km +#define C 299792458.0 // Speed of light (m/s) +#define SQRTPI sqrt(PI) +#define THIRD 1.0 / 3.0 +#define FALSE 0 +#define TRUE 1 +#define D2R PI/180.0 +#define R2D 180.0/PI +#define ETA 119.9169832*PI // Intrinsic impedance of free space (ohms) + +#define POLARIZATION__HORIZONTAL 0 +#define POLARIZATION__VERTICAL 1 + +#define METHOD__FLAT_EARTH_CURVE 0 +#define METHOD__RESIUDE_SERIES 1 + +#define WERF_TERMS 30 // Number of terms for the complex error function + +#define YES 1 // Find the derivative i.e., Ai'() or Bi'() +#define NO 0 // Find Ai() or Bi() +// kind +#define AIRY 1 // Find the Airy Function +#define AIRYD 2 // Find the Airy function Derivative +#define BAIRY 3 // Find the Bairy Function +#define BAIRYD 4 // Find the Bairy function Derivative +#define WTWO 5 // find Hufford Wi(2) or Wait W2 +#define DWTWO 6 // find Hufford Wi'(2) or Wait W2' +#define WONE 7 // find Hufford Wi(1) or Wait W1 +#define DWONE 8 // find Hufford Wi'(1) or Wait W1' +// scaling +#define HUFFORD 9 // Use Hufford scaling +#define WAIT 10 // Use Wait scaling +#define NONE 11 // No Scaling + +// Return codes +#define SUCCESS 0 + +// Error codes +#define ERROR__TX_TERMINAL_HEIGHT 1000 // TX terminal height is out of range +#define ERROR__RX_TERMINAL_HEIGHT 1001 // RX terminal height is out of range +#define ERROR__FREQUENCY 1002 // Frequency is out of range +#define ERROR__TX_POWER 1003 // Transmit power is out of range +#define ERROR__SURFACE_REFRACTIVITY 1004 // Surface refractivity is out of range +#define ERROR__PATH_DISTANCE 1005 // Path distance is out of range +#define ERROR__EPSILON 1006 // Epsilon is out of range +#define ERROR__SIGMA 1007 // Sigma is out of range +#define ERROR__POLARIZATION 1008 // Invalid value for polarization + +////////////////////////////////////// +// Data Structures + +struct Result +{ + double A_btl__db; + double E_dBuVm; + double P_rx__dbm; + + int method; +}; + +////////////////////////////////////// +// Main LFMF Function + +DLLEXPORT int LFMF(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol, Result *result); + +////////////////////////////////////// +// Helper Functions + +double FlatEarthCurveCorrection(complex delta, complex q, double h_1__km, double h_2__km, double d, double k); +double ResidueSeries(double d, double k, double h_1__km, double h_2__km, double nu, double theta, complex q); +complex werf(complex qi); +complex Airy(complex Z, int kind, int scaling); +complex WiRoot(int i, complex *DWi, complex q, complex *Wi, int kind, int scaling); +int ValidateInput(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol); diff --git a/src/Airy.cpp b/src/Airy.cpp new file mode 100644 index 0000000..37355a1 --- /dev/null +++ b/src/Airy.cpp @@ -0,0 +1,818 @@ +#include "..\include\LFMF.h" + +/*============================================================================= + | + | Description: This routine finds the Airy, Bariy, Wi(1) and Wi(2) + | functions and their derivatives for a complex input argument + | from a shifted Taylor series or by asymptotic approximation + | depending of the location of the input argument. + | + | This routine determines the so-called "Airy Functions of the + | third kind" Wi(1) and Wi(2) that are found in equation 38 + | of NTIA Report 87-219 "A General Theory of Radio + | Propagation through a Stratified Atmosphere", George + | Hufford, July 1987 + | + | The Airy function that appeared in the original GWINT and + | GWRES had the switches all mangled from what George Hufford + | had in mind this routine has the corrected switches. Please + | see the Airy function code that appears in the appendix of + | OT/ITS RR 11 "A Wave Hop Propagation Program for an + | Anisotropic Ionosphere" L. A. Berry and J. E. Herman + | April 1971 + | + | Input: Z - Input argument + | kind - Switch that indicates what type of Airy + | function to solve for + | scaling - + | + | Outputs: [None] + | + | Returns: Ai - The desired Airy function calculated at Z + | + | Note: A note on scaling the output from this program + | + | There is a definitional problem with the Airy function + | which is inevitable relative to how it was defined in the + | original LFMF code originated with the Hufford's AIRY + | subroutine. + | + | Using the scaling equal to HUFFORD in this program follows + | the definitions of Wi(1) and Wi(2) as defined by Hufford + | (87-219) + | + | Using the scaling equal to WAIT in this program uses the + | definitions of W1 and W2 defined in Deminco (99-368) and + | in the original LFMF code following Berry via Wait. + | + | The two solutions differ by a constant. As Hufford notes + | concerning Wi(1) and Wi(2) in 87-219 + | + | "Except for multiplicative constants they correspond to + | what Fock (1965) calls w1 and w2 and to what Wait (1962) + | calls w2 and w1" + | + | The following are the multiplicative constants that allow + | for the translation between Hufford Wi(2) and Wi(1) with + | Wait W1 and W2, respectively. These are given here as a + | reference if this function is used for programs other + | than LFMF. + | + | // Wait + | complex WW2 = complex( sqrt(3.0*PI), sqrt(PI)); + | complex WDW2 = complex(-1.0*sqrt(3.0*PI), sqrt(PI)); + | complex WW1 = complex( sqrt(3.0*PI), -1.0*sqrt(PI)); + | complex WDW1 = complex(-1.0*sqrt(3.0*PI), -1.0*sqrt(PI)); + | + | // Hufford + | complex HW2 = 2.0*complex(cos( PI/3.0), sin( PI/3.0)); + | complex HDW2 = 2.0*complex(cos(-PI/3.0), sin(-PI/3.0)); + | complex HW1 = 2.0*complex(cos(-PI/3.0), sin(-PI/3.0)); + | complex HDW1 = 2.0*complex(cos( PI/3.0), sin( PI/3.0)); + | + | // (Multiplicative constant) * Huffords Wi'(1) = Wait W1' + | // So the multiplicative constants are: + | complex uDW2 = WDW2/HDW1; // uDW2 = complex(0.0, sqrt(PI)) + | complex uW2 = WW2/HW1; // uW2 = complex(0.0, sqrt(PI)) + | complex uDW1 = WDW1/HDW2; // uDW1 = complex(0.0, -sqrt(PI)) + | complex uW1 = WW1/HW2; // uW1 = complex(0.0, -sqrt(PI)) + | + | To make the solutions that are generated by this program + | for the Hufford Airy functions of the "3rd kind" abundantly + | clear please examine the following examples. + | + | For Z = 8.0 + 8.0 i the Asymptotic Solution is used + | + | Ai( 8.0 + 8.0 i) = 6.576933e-007 + 9.312331e-006 i + | Ai'(8.0 + 8.0 i) = 9.79016e-006 + -2.992170e-005 i + | Bi( 8.0 + 8.0 i) = -1.605154e+003 + -4.807200e+003 i + | Bi'(8.0 + 8.0 i) = 1301.23 + -16956 i + | Wi(1)(8.0 + 8.0 i) = -4.807200e+003 + 1.605154e+003 i + | Wi(2)(8.0 + 8.0 i) = 4.807200e+003 + -1.605154e+003 i + | Ai(z) - j*Bi(z) = -4.807200e+003 + 1.605154e+003 i + | Ai(z) + j*Bi(z) = 4.807200e+003 + -1.605154e+003 i + | + | For Z = 1.0 - 2.0 i the Taylor series with a shifted + | center of expansion solution used. + | + | Ai( 1.0 - 2.0 i) = -2.193862e-001 + 1.753859e-001 i + | Ai'(1.0 - 2.0 i) = 0.170445 + -0.387622 i + | Bi( 1.0 - 2.0 i) = 4.882205e-002 + -1.332740e-001 i + | Bi'(1.0 - 2.0 i) = -0.857239 + -0.495506 i + | Wi(1)(1.0 - 2.0 i) = -3.526603e-001 + 1.265638e-001 i + | Wi(2)(1.0 - 2.0 i) = -8.611221e-002 + 2.242079e-001 i + | Ai(z) - j*Bi(z) = -3.526603e-001 + 1.265639e-001 i + | Ai(z) + j*Bi(z) = -8.611221e-002 + 2.242080e-001 i + | + *===========================================================================*/ +complex Airy(complex Z, int kind, int scaling) +{ + // NQTT, ASLT data + int NQTT[15] = { 1,3,7,12,17,23,29,35,41,47,53,59,64,68,71 }; // Centers of Expansion of Taylor series on real axis indices into the + // AV, APV, BV and BPV arrays + int N; // Index into NQTT[] array + int NQ8; // Index that indicates the radius of convergence of the Taylor series solution + int CoERealidx; // Center of Expansion of the Taylor Series real index + int CoEImagidx; // Center of Expansion of the Taylor Series imaginary index + int cnt; // loop counter for the Taylor series calculation + int derivative; // index for derivative + + bool reflection; // Flag to indicate that the answer needs to be flipped over since this routine only finds solutions in quadrant 1 and 2 + + complex A[2], ZT, B0, B1, B2, B3, AN, U, AIRYW, ZA, ZB, ZE, ZR, V, ZV, ZU, PHZU; // Temps + complex CoE; // Center of Expansion of the Taylor series + complex Ai; // Ai is either Ai(at the center of expansion of the Taylor series) or Bi( at the center of expansion of the Taylor series ) + complex Aip; // Aip is the derivative of the above + complex sum1; // Temp Sum for the asymptotic solution + complex sum2; // Temp Sum for the asymptotic solution + complex ZB2, ZB1; + + double one; // Used in the calculation of the asymptotic solution is either -1 or 1 + + // terms for asymptotic series. second column is for derivative + int SIZE_OF_ASV = 15; + double ASV[15][2] = { {0.5989251E+5, -0.6133571E+5}, + {0.9207207E+4, -0.9446355E+4}, + {0.1533169E+4, -0.1576357E+4}, + {0.2784651E+3, -0.2870332E+3}, + {0.5562279E+2, -0.5750830E+2}, + {0.1234157E+2, -0.1280729E+2}, + {0.3079453E+1, -0.3210494E+1}, + {0.8776670E+0, -0.9204800E+0}, + {0.2915914E+0, -0.3082538E+0}, + {0.1160991E+0, -0.1241059E+0}, + {0.5764919E-1, -0.6266216E-1}, + {0.3799306E-1, -0.4246283E-1}, + {0.3713349E-1, -0.4388503E-1}, + {0.6944444E-1, -0.9722222E-1}, + {0.1000000E+1, 0.1000000E+1} }; + + // Complex array of the value Ai(a) Airy function + // where a is the complex location of the center of expansion of the Taylor series + complex AV[70]; + + // Complex array of the value Ai'(a) derivative of the Airy function + // where a is the complex location of the center of expansion of the Taylor series + // presumably for Ai'[z] + complex APV[70]; + + ////////////////////////////////////////////////////////////////////////// + // Initialize the center of expansion arrays. // + ////////////////////////////////////////////////////////////////////////// + // The array AV[] (and BV[]) is the Airy function for Ai(a) to shift // + // the Taylor series from the origin to the point a Thus the series is // + // f(z -a) // + // Why George Hufford choose these particular locations for the // + // of the centers of expansion is unknown. The centers of expansion are // + // included here to remove any ambiguity in the method. // + ////////////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////// Center of expansion + /////////////////////////////////////////////////// ( Real, Imaginary) + AV[0] = complex(-3.2914520e-001, +0.0000000e+000);// (-6,0) + AV[1] = complex(-2.6780040e+000, +1.4774590e+000);// (-6,1/sin(pi/3)) + AV[2] = complex(+3.5076100e-001, +0.0000000e+000);// (-5,0) + AV[3] = complex(+2.4122260e+000, +6.9865120e-001);// (-5,1/sin(pi/3)) + AV[4] = complex(+3.3635530e+001, -3.4600960e+000);// (-5,2/sin(pi/3)) + AV[5] = complex(+3.4449740e+002, -3.3690890e+002);// (-5,3/sin(pi/3)) + AV[6] = complex(-7.0265530e-002, +0.0000000e+000);// (-4,0) + AV[7] = complex(-5.4818220e-001, -1.9207370e+000);// (-4,1/sin(pi/3)) + AV[8] = complex(-1.3383400e+001, -1.6022590e+001);// (-4,2/sin(pi/3)) + AV[9] = complex(-2.2967800e+002, -3.2072450e+001);// (-4,3/sin(pi/3)) + AV[10] = complex(-1.8040780e+003, +2.1917680e+003);// (-4,4/sin(pi/3)) + AV[11] = complex(-3.7881430e-001, +0.0000000e+000);// (-3,0) + AV[12] = complex(-1.3491840e+000, +8.4969080e-001);// (-3,1/sin(pi/3)) + AV[13] = complex(-6.0453340e+000, +1.0623180e+001);// (-3,2/sin(pi/3)) + AV[14] = complex(+3.1169620e+001, +9.8813520e+001);// (-3,3/sin(pi/3)) + AV[15] = complex(+9.8925350e+002, +1.3905290e+002);// (-3,4/sin(pi/3)) + AV[16] = complex(+2.2740740e-001, +0.0000000e+000);// (-2,0) + AV[17] = complex(+7.1857400e-001, +9.7809090e-001);// (-2,1/sin(pi/3)) + AV[18] = complex(+6.0621090e+000, +2.7203010e+000);// (-2,2/sin(pi/3)) + AV[19] = complex(+3.6307080e+001, -2.0961360e+001);// (-2,3/sin(pi/3)) + AV[20] = complex(-6.7139790e+001, -3.0904640e+002);// (-2,4/sin(pi/3)) + AV[21] = complex(-2.8001650e+003, +4.6649370e+002);// (-2,5/sin(pi/3)) + AV[22] = complex(+5.3556090e-001, +0.0000000e+000);// (-1,0) + AV[23] = complex(+9.2407370e-001, -1.9106560e-001);// (-1,1/sin(pi/3)) + AV[24] = complex(+1.8716190e+000, -2.5743310e+000);// (-1,2/sin(pi/3)) + AV[25] = complex(-7.2188440e+000, -1.2924200e+001);// (-1,3/sin(pi/3)) + AV[26] = complex(-8.1787380e+001, +3.2087010e+001);// (-1,4/sin(pi/3)) + AV[27] = complex(+2.9933950e+002, +5.6922180e+002);// (-1,5/sin(pi/3)) + AV[28] = complex(+3.5502810e-001, +0.0000000e+000);// ( 0,0) + AV[29] = complex(+3.1203440e-001, -3.8845390e-001);// ( 0,1/sin(pi/3)) + AV[30] = complex(-5.2840000e-001, -1.0976410e+000);// ( 0,2/sin(pi/3)) + AV[31] = complex(-4.2009350e+000, +1.1940150e+000);// ( 0,3/sin(pi/3)) + AV[32] = complex(+7.1858830e+000, +1.9600910e+001);// ( 0,4/sin(pi/3)) + AV[33] = complex(+1.0129120e+002, -7.5951230e+001);// ( 0,5/sin(pi/3)) + AV[34] = complex(+1.3529240e-001, +0.0000000e+000);// ( 1,0) + AV[35] = complex(+3.2618480e-002, -1.7084870e-001);// ( 1,1/sin(pi/3)) + AV[36] = complex(-3.4215380e-001, -8.9067650e-002);// ( 1,2/sin(pi/3)) + AV[37] = complex(-1.4509640e-001, +1.0328020e+000);// ( 1,3/sin(pi/3)) + AV[38] = complex(+4.1001970e+000, -6.8936910e-001);// ( 1,4/sin(pi/3)) + AV[39] = complex(-1.3030120e+001, -1.6910540e+001);// ( 1,5/sin(pi/3)) + AV[40] = complex(+3.4924130e-002, +0.0000000e+000);// ( 2,0) + AV[41] = complex(-8.4464730e-003, -4.2045150e-002);// ( 2,1/sin(pi/3)) + AV[42] = complex(-6.9313270e-002, +3.5364800e-002);// ( 2,2/sin(pi/3)) + AV[43] = complex(+1.5227620e-001, +1.2848450e-001);// ( 2,3/sin(pi/3)) + AV[44] = complex(+1.0681370e-001, -6.7766150e-001);// ( 2,4/sin(pi/3)) + AV[45] = complex(-2.6193430e+000, +1.5699860e+000);// ( 2,5/sin(pi/3)) + AV[46] = complex(+6.5911390e-003, +0.0000000e+000);// ( 3,0) + AV[47] = complex(-3.9443990e-003, -6.8060110e-003);// ( 3,1/sin(pi/3)) + AV[48] = complex(-5.9820130e-003, +1.1799010e-002);// ( 3,2/sin(pi/3)) + AV[49] = complex(+2.9922500e-002, -5.9772930e-003);// ( 3,3/sin(pi/3)) + AV[50] = complex(-7.7464130e-002, -5.2292400e-002);// ( 3,4/sin(pi/3)) + AV[51] = complex(+1.1276590e-001, +3.5112440e-001);// ( 3,5/sin(pi/3)) + AV[52] = complex(+9.5156390e-004, +0.0000000e+000);// ( 4,0) + AV[53] = complex(-8.0843000e-004, -7.6590130e-004);// ( 4,1/sin(pi/3)) + AV[54] = complex(+1.6147820e-004, +1.7661760e-003);// ( 4,2/sin(pi/3)) + AV[55] = complex(+2.0138720e-003, -3.1976720e-003);// ( 4,3/sin(pi/3)) + AV[56] = complex(-9.5086780e-003, +4.5377830e-003);// ( 4,4/sin(pi/3)) + AV[57] = complex(+3.7560190e-002, +5.7361920e-004);// ( 4,5/sin(pi/3)) + AV[58] = complex(+1.0834440e-004, +0.0000000e+000);// ( 5,0) + AV[59] = complex(-1.0968610e-004, -5.9902330e-005);// ( 5,1/sin(pi/3)) + AV[60] = complex(+1.0778190e-004, +1.5771600e-004);// ( 5,2/sin(pi/3)) + AV[61] = complex(-6.8980940e-005, -3.7626460e-004);// ( 5,3/sin(pi/3)) + AV[62] = complex(-1.6166130e-004, +9.7457770e-004);// ( 5,4/sin(pi/3)) + AV[63] = complex(+9.9476940e-006, +0.0000000e+000);// ( 6,0) + AV[64] = complex(-1.0956820e-005, -2.9508800e-006);// ( 6,1/sin(pi/3)) + AV[65] = complex(+1.4709070e-005, +8.1042090e-006);// ( 6,2/sin(pi/3)) + AV[66] = complex(-2.4446020e-005, -2.0638140e-005);// ( 6,3/sin(pi/3)) + AV[67] = complex(+7.4921290e-007, +0.0000000e+000);// ( 7,0) + AV[68] = complex(-8.4619070e-007, -3.6807340e-008);// ( 7,1/sin(pi/3)) + AV[69] = complex(+1.2183960e-006, +8.3589200e-008);// ( 7,2/sin(pi/3)) + + ////////////////////////////////////////////////////////////////////////// + // This array APV[] is the derivative of the Airy function for Ai'(a) // + ////////////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////// Center of expansion + /////////////////////////////////////////////////// ( Real, Imaginary) + APV[0] = complex(+3.4593550e-001, +0.0000000e+000);// (-6,0) + APV[1] = complex(+4.1708880e+000, +6.2414440e+000);// (-6,1/sin(pi/3)) + APV[2] = complex(+3.2719280e-001, +0.0000000e+000);// (-5,0) + APV[3] = complex(+1.0828740e+000, -5.4928300e+000);// (-5,1/sin(pi/3)) + APV[4] = complex(-2.3363520e+001, -7.4901850e+001);// (-5,2/sin(pi/3)) + APV[5] = complex(-1.0264880e+003, -5.6707940e+002);// (-5,3/sin(pi/3)) + APV[6] = complex(-7.9062860e-001, +0.0000000e+000);// (-4,0) + APV[7] = complex(-3.8085830e+000, +1.5129610e+000);// (-4,1/sin(pi/3)) + APV[8] = complex(-2.6086380e+001, +3.5540710e+001);// (-4,2/sin(pi/3)) + APV[9] = complex(+1.0761840e+002, +5.1239940e+002);// (-4,3/sin(pi/3)) + APV[10] = complex(+6.6597800e+003, +1.8096190e+003);// (-4,4/sin(pi/3)) + APV[11] = complex(+3.1458380e-001, +0.0000000e+000);// (-3,0) + APV[12] = complex(+1.8715430e+000, +2.0544840e+000);// (-3,1/sin(pi/3)) + APV[13] = complex(+2.2591740e+001, +4.8563000e+000);// (-3,2/sin(pi/3)) + APV[14] = complex(+1.6163000e+002, -1.4335600e+002);// (-3,3/sin(pi/3)) + APV[15] = complex(-8.0047160e+002, -2.1527450e+003);// (-3,4/sin(pi/3)) + APV[16] = complex(+6.1825900e-001, +0.0000000e+000);// (-2,0) + APV[17] = complex(+1.3019600e+000, -1.2290770e+000);// (-2,1/sin(pi/3)) + APV[18] = complex(+1.5036120e-001, -1.1008090e+001);// (-2,2/sin(pi/3)) + APV[19] = complex(-7.0116800e+001, -4.0480820e+001);// (-2,3/sin(pi/3)) + APV[20] = complex(-4.8317170e+002, +4.9692760e+002);// (-2,4/sin(pi/3)) + APV[21] = complex(+4.8970660e+003, +4.8627290e+003);// (-2,5/sin(pi/3)) + APV[22] = complex(-1.0160570e-002, +0.0000000e+000);// (-1,0) + APV[23] = complex(-5.4826640e-001, -7.1365290e-001);// (-1,1/sin(pi/3)) + APV[24] = complex(-4.6749130e+000, -1.1924250e-001);// (-1,2/sin(pi/3)) + APV[25] = complex(-1.0536400e+001, +2.4943710e+001);// (-1,3/sin(pi/3)) + APV[26] = complex(+1.6333770e+002, +9.0394910e+001);// (-1,4/sin(pi/3)) + APV[27] = complex(+5.6449460e+002, -1.4248320e+003);// (-1,5/sin(pi/3)) + APV[28] = complex(-2.5881940e-001, +0.0000000e+000);// ( 0,0) + APV[29] = complex(-4.8620750e-001, +1.5689920e-001);// ( 0,1/sin(pi/3)) + APV[30] = complex(-4.7348130e-001, +1.7093440e+000);// ( 0,2/sin(pi/3)) + APV[31] = complex(+7.0373840e+000, +3.6281820e+000);// ( 0,3/sin(pi/3)) + APV[32] = complex(+1.7739590e+001, -4.0360420e+001);// ( 0,4/sin(pi/3)) + APV[33] = complex(-2.9791510e+002, -3.8408890e+001);// ( 0,5/sin(pi/3)) + APV[34] = complex(-1.5914740e-001, +0.0000000e+000);// ( 1,0) + APV[35] = complex(-1.1340420e-001, +1.9730500e-001);// ( 1,1/sin(pi/3)) + APV[36] = complex(+4.0126210e-001, +3.9223000e-001);// ( 1,2/sin(pi/3)) + APV[37] = complex(+1.3348650e+000, -1.4377270e+000);// ( 1,3/sin(pi/3)) + APV[38] = complex(-7.9022490e+000, -4.2063640e+000);// ( 1,4/sin(pi/3)) + APV[39] = complex(-1.3892750e+000, +5.1229420e+001);// ( 1,5/sin(pi/3)) + APV[40] = complex(-5.3090380e-002, +0.0000000e+000);// ( 2,0) + APV[41] = complex(-1.6832970e-003, +6.8366970e-002);// ( 2,1/sin(pi/3)) + APV[42] = complex(+1.3789400e-001, -1.1613800e-002);// ( 2,2/sin(pi/3)) + APV[43] = complex(-1.4713730e-001, -3.7151990e-001);// ( 2,3/sin(pi/3)) + APV[44] = complex(-1.0070200e+000, +1.1591350e+000);// ( 2,4/sin(pi/3)) + APV[45] = complex(+7.5045050e+000, +4.6913120e-001);// ( 2,5/sin(pi/3)) + APV[46] = complex(-1.1912980e-002, +0.0000000e+000);// ( 3,0) + APV[47] = complex(+5.1468570e-003, +1.3660890e-002);// ( 3,1/sin(pi/3)) + APV[48] = complex(+1.8309710e-002, -1.8808590e-002);// ( 3,2/sin(pi/3)) + APV[49] = complex(-6.4461590e-002, -1.3611790e-002);// ( 3,3/sin(pi/3)) + APV[50] = complex(+1.0516240e-001, +1.9313050e-001);// ( 3,4/sin(pi/3)) + APV[51] = complex(+2.0520050e-001, -9.1772620e-001);// ( 3,5/sin(pi/3)) + APV[52] = complex(-1.9586410e-003, +0.0000000e+000);// ( 4,0) + APV[53] = complex(+1.4695650e-003, +1.8086380e-003);// ( 4,1/sin(pi/3)) + APV[54] = complex(+5.9709950e-004, -3.8332700e-003);// ( 4,2/sin(pi/3)) + APV[55] = complex(-6.8910890e-003, +5.4467430e-003);// ( 4,3/sin(pi/3)) + APV[56] = complex(+2.6167930e-002, -8.4092000e-004);// ( 4,4/sin(pi/3)) + APV[57] = complex(-8.8284470e-002, -4.6475310e-002);// ( 4,5/sin(pi/3)) + APV[58] = complex(-2.4741390e-004, +0.0000000e+000);// ( 5,0) + APV[59] = complex(+2.3707840e-004, +1.6461110e-004);// ( 5,1/sin(pi/3)) + APV[60] = complex(-1.7465570e-004, -4.2026780e-004);// ( 5,2/sin(pi/3)) + APV[61] = complex(-1.0394520e-004, +9.4761840e-004);// ( 5,3/sin(pi/3)) + APV[62] = complex(+1.3004110e-003, -2.2446660e-003);// ( 5,4/sin(pi/3)) + APV[63] = complex(-2.4765200e-005, +0.0000000e+000);// ( 6,0) + APV[64] = complex(+2.6714870e-005, +9.8691570e-006);// ( 6,1/sin(pi/3)) + APV[65] = complex(-3.3539770e-005, -2.7113280e-005);// ( 6,2/sin(pi/3)) + APV[66] = complex(+4.9197840e-005, +6.9349090e-005);// ( 6,3/sin(pi/3)) + APV[67] = complex(-2.0081510e-006, +0.0000000e+000);// ( 7,0) + APV[68] = complex(+2.2671240e-006, +2.7848510e-007);// ( 7,1/sin(pi/3)) + APV[69] = complex(-3.2692130e-006, -7.3943490e-007);// ( 7,2/sin(pi/3)) + + ///////////////////////////////////////////////////////////////////////// + + // Complex array of the value Bi(a) Airy function + // where a is the complex location of the center of expansion of the Taylor series + complex BV[70]; + + // Complex array of the value Bi'(a) derivative of the Airy function + // where a is the complex location of the center of expansion of the Taylor series + // presumably for Ai'[z] + complex BPV[70]; + + BV[0] = complex(-1.466984e-001, -9.813078e-017);// (-6,0) + BV[1] = complex(-1.489391e+000, -2.660635e+000);// (-6,1/sin(pi/3)) + BV[2] = complex(-1.383691e-001, +0.000000e+000);// (-5,0) + BV[3] = complex(-7.034482e-001, +2.384547e+000);// (-5,1/sin(pi/3)) + BV[4] = complex(+3.460723e+000, +3.363363e+001);// (-5,2/sin(pi/3)) + BV[5] = complex(+3.369090e+002, +3.444973e+002);// (-5,3/sin(pi/3)) + BV[6] = complex(+3.922347e-001, -1.041605e-016);// (-4,0) + BV[7] = complex(+1.956219e+000, -5.327226e-001);// (-4,1/sin(pi/3)) + BV[8] = complex(+1.602464e+001, -1.338050e+001);// (-4,2/sin(pi/3)) + BV[9] = complex(+3.207239e+001, -2.296777e+002);// (-4,3/sin(pi/3)) + BV[10] = complex(-2.191768e+003, -1.804078e+003);// (-4,4/sin(pi/3)) + BV[11] = complex(-1.982896e-001, +4.440892e-016);// (-3,0) + BV[12] = complex(-8.880754e-001, -1.308713e+000);// (-3,1/sin(pi/3)) + BV[13] = complex(-1.062975e+001, -6.044056e+000);// (-3,2/sin(pi/3)) + BV[14] = complex(-9.881405e+001, +3.116914e+001);// (-3,3/sin(pi/3)) + BV[15] = complex(-1.390528e+002, +9.892534e+002);// (-3,4/sin(pi/3)) + BV[16] = complex(-4.123026e-001, +1.451806e-016);// (-2,0) + BV[17] = complex(-1.034766e+000, +6.541962e-001);// (-2,1/sin(pi/3)) + BV[18] = complex(-2.720266e+000, +6.048328e+000);// (-2,2/sin(pi/3)) + BV[19] = complex(+2.096300e+001, +3.630613e+001);// (-2,3/sin(pi/3)) + BV[20] = complex(+3.090465e+002, -6.713963e+001);// (-2,4/sin(pi/3)) + BV[21] = complex(-4.664937e+002, -2.800165e+003);// (-2,5/sin(pi/3)) + BV[22] = complex(+1.039974e-001, +0.000000e+000);// (-1,0) + BV[23] = complex(+2.797458e-001, +8.086491e-001);// (-1,1/sin(pi/3)) + BV[24] = complex(+2.606133e+000, +1.870297e+000);// (-1,2/sin(pi/3)) + BV[25] = complex(+1.292648e+001, -7.213647e+000);// (-1,3/sin(pi/3)) + BV[26] = complex(-3.208774e+001, -8.178697e+001);// (-1,4/sin(pi/3)) + BV[27] = complex(-5.692218e+002, +2.993394e+002);// (-1,5/sin(pi/3)) + BV[28] = complex(+6.149266e-001, +0.000000e+000);// ( 0,0) + BV[29] = complex(+6.732023e-001, +3.575876e-001);// ( 0,1/sin(pi/3)) + BV[30] = complex(+1.125057e+000, -4.471292e-001);// ( 0,2/sin(pi/3)) + BV[31] = complex(-1.211148e+000, -4.191469e+000);// ( 0,3/sin(pi/3)) + BV[32] = complex(-1.960240e+001, +7.182663e+000);// ( 0,4/sin(pi/3)) + BV[33] = complex(+7.595175e+001, +1.012911e+002);// ( 0,5/sin(pi/3)) + BV[34] = complex(+1.207424e+000, +0.000000e+000);// ( 1,0) + BV[35] = complex(+5.951440e-001, +6.156664e-001);// ( 1,1/sin(pi/3)) + BV[36] = complex(-1.002325e-001, -1.338228e-001);// ( 1,2/sin(pi/3)) + BV[37] = complex(-1.089323e+000, -2.019524e-001);// ( 1,3/sin(pi/3)) + BV[38] = complex(+7.047139e-001, +4.091592e+000);// ( 1,4/sin(pi/3)) + BV[39] = complex(+1.691067e+001, -1.302705e+001);// ( 1,5/sin(pi/3)) + BV[40] = complex(+3.298095e+000, +0.000000e+000);// ( 2,0) + BV[41] = complex(+2.244706e-001, +2.421124e+000);// ( 2,1/sin(pi/3)) + BV[42] = complex(-1.199515e+000, -1.167656e-001);// ( 2,2/sin(pi/3)) + BV[43] = complex(+6.781072e-003, -2.225418e-001);// ( 2,3/sin(pi/3)) + BV[44] = complex(+7.470822e-001, +1.832986e-001);// ( 2,4/sin(pi/3)) + BV[45] = complex(-1.590993e+000, -2.617694e+000);// ( 2,5/sin(pi/3)) + BV[46] = complex(+1.403733e+001, +0.000000e+000);// ( 3,0) + BV[47] = complex(-3.731398e+000, +1.066394e+001);// ( 3,1/sin(pi/3)) + BV[48] = complex(-4.440986e+000, -4.309647e+000);// ( 3,2/sin(pi/3)) + BV[49] = complex(+2.373933e+000, -5.300179e-001);// ( 3,3/sin(pi/3)) + BV[50] = complex(-2.821481e-001, +5.657373e-001);// ( 3,4/sin(pi/3)) + BV[51] = complex(-3.904913e-001, -5.168316e-002);// ( 3,5/sin(pi/3)) + BV[52] = complex(+8.384707e+001, +0.000000e+000);// ( 4,0) + BV[53] = complex(-4.356467e+001, +5.497027e+001);// ( 4,1/sin(pi/3)) + BV[54] = complex(-7.156364e+000, -4.113550e+001);// ( 4,2/sin(pi/3)) + BV[55] = complex(+1.455852e+001, +1.109071e+001);// ( 4,3/sin(pi/3)) + BV[56] = complex(-6.111359e+000, -1.094609e-001);// ( 4,4/sin(pi/3)) + BV[57] = complex(+1.403434e+000, -7.255043e-001);// ( 4,5/sin(pi/3)) + BV[58] = complex(+6.577920e+002, +0.000000e+000);// ( 5,0) + BV[59] = complex(-4.598656e+002, +3.242259e+002);// ( 5,1/sin(pi/3)) + BV[60] = complex(+1.324505e+002, -3.294705e+002);// ( 5,2/sin(pi/3)) + BV[61] = complex(+2.057579e+001, +1.674034e+002);// ( 5,3/sin(pi/3)) + BV[62] = complex(-3.161505e+001, -5.302141e+001);// ( 5,4/sin(pi/3)) + BV[63] = complex(+6.536446e+003, +0.000000e+000);// ( 6,0) + BV[64] = complex(-5.316522e+003, +1.992175e+003);// ( 6,1/sin(pi/3)) + BV[65] = complex(+2.888529e+003, -2.373473e+003);// ( 6,2/sin(pi/3)) + BV[66] = complex(-1.078657e+003, +1.551931e+003);// ( 6,3/sin(pi/3)) + BV[67] = complex(+8.032779e+004, +0.000000e+000);// ( 7,0) + BV[68] = complex(-7.001987e+004, +8.828416e+003);// ( 7,1/sin(pi/3)) + BV[69] = complex(+4.676699e+004, -1.085924e+004);// ( 7,2/sin(pi/3)) + + ////////////////////////////////////////////////////////////// + + BPV[0] = complex(-8.128988e-001, +3.365185e-016);// (-6,0) + BPV[1] = complex(-6.287609e+000, +4.146176e+000);// (-6,1/sin(pi/3)) + BPV[2] = complex(+7.784118e-001, -3.004629e-016);// (-5,0) + BPV[3] = complex(+5.554036e+000, +1.063645e+000);// (-5,1/sin(pi/3)) + BPV[4] = complex(+7.490659e+001, -2.336310e+001);// (-5,2/sin(pi/3)) + BPV[5] = complex(+5.670796e+002, -1.026488e+003);// (-5,3/sin(pi/3)) + BPV[6] = complex(-1.166706e-001, +2.371654e-015);// (-4,0) + BPV[7] = complex(-1.532413e+000, -3.730947e+000);// (-4,1/sin(pi/3)) + BPV[8] = complex(-3.554558e+001, -2.608034e+001);// (-4,2/sin(pi/3)) + BPV[9] = complex(-5.124001e+002, +1.076185e+002);// (-4,3/sin(pi/3)) + BPV[10] = complex(-1.809619e+003, +6.659780e+003);// (-4,4/sin(pi/3)) + BPV[11] = complex(-6.756112e-001, -2.403703e-017);// (-3,0) + BPV[12] = complex(-2.142202e+000, +1.818610e+000);// (-3,1/sin(pi/3)) + BPV[13] = complex(-4.863137e+000, +2.258023e+001);// (-3,2/sin(pi/3)) + BPV[14] = complex(+1.433564e+002, +1.616285e+002);// (-3,3/sin(pi/3)) + BPV[15] = complex(+2.152746e+003, -8.004716e+002);// (-3,4/sin(pi/3)) + BPV[16] = complex(+2.787952e-001, +0.000000e+000);// (-2,0) + BPV[17] = complex(+1.300360e+000, +1.185229e+000);// (-2,1/sin(pi/3)) + BPV[18] = complex(+1.103082e+001, +1.397575e-001);// (-2,2/sin(pi/3)) + BPV[19] = complex(+4.048422e+001, -7.011484e+001);// (-2,3/sin(pi/3)) + BPV[20] = complex(-4.969277e+002, -4.831712e+002);// (-2,4/sin(pi/3)) + BPV[21] = complex(-4.862729e+003, +4.897066e+003);// (-2,5/sin(pi/3)) + BPV[22] = complex(+5.923756e-001, +0.000000e+000);// (-1,0) + BPV[23] = complex(+9.080502e-001, -5.080757e-001);// (-1,1/sin(pi/3)) + BPV[24] = complex(+1.499485e-001, -4.631403e+000);// (-1,2/sin(pi/3)) + BPV[25] = complex(-2.494926e+001, -1.052676e+001);// (-1,3/sin(pi/3)) + BPV[26] = complex(-9.039663e+001, +1.633370e+002);// (-1,4/sin(pi/3)) + BPV[27] = complex(+1.424833e+003, +5.644943e+002);// (-1,5/sin(pi/3)) + BPV[28] = complex(+4.482884e-001, +0.000000e+000);// ( 0,0) + BPV[29] = complex(+2.493288e-002, -1.876446e-001);// ( 0,1/sin(pi/3)) + BPV[30] = complex(-1.774795e+000, -3.533830e-001);// ( 0,2/sin(pi/3)) + BPV[31] = complex(-3.663891e+000, +7.026174e+000);// ( 0,3/sin(pi/3)) + BPV[32] = complex(+4.036322e+001, +1.773235e+001);// ( 0,4/sin(pi/3)) + BPV[33] = complex(+3.840990e+001, -2.979143e+002);// ( 0,5/sin(pi/3)) + BPV[34] = complex(+9.324359e-001, +0.000000e+000);// ( 1,0) + BPV[35] = complex(-1.293870e-001, +7.817697e-001);// ( 1,1/sin(pi/3)) + BPV[36] = complex(-8.385825e-001, +4.901385e-001);// ( 1,2/sin(pi/3)) + BPV[37] = complex(+1.421331e+000, +1.181168e+000);// ( 1,3/sin(pi/3)) + BPV[38] = complex(+4.244380e+000, -7.895016e+000);// ( 1,4/sin(pi/3)) + BPV[39] = complex(-5.123410e+001, -1.383387e+000);// ( 1,5/sin(pi/3)) + BPV[40] = complex(+4.100682e+000, +0.000000e+000);// ( 2,0) + BPV[41] = complex(-9.576171e-001, +3.432468e+000);// ( 2,1/sin(pi/3)) + BPV[42] = complex(-1.747487e+000, -8.602854e-001);// ( 2,2/sin(pi/3)) + BPV[43] = complex(+9.978890e-001, -6.434913e-001);// ( 2,3/sin(pi/3)) + BPV[44] = complex(-1.127841e+000, -7.762214e-001);// ( 2,4/sin(pi/3)) + BPV[45] = complex(-5.136144e-001, +7.476880e+000);// ( 2,5/sin(pi/3)) + BPV[46] = complex(+2.292221e+001, +0.000000e+000);// ( 3,0) + BPV[47] = complex(-1.021000e+001, +1.662556e+001);// ( 3,1/sin(pi/3)) + BPV[48] = complex(-5.018884e+000, -1.067168e+001);// ( 3,2/sin(pi/3)) + BPV[49] = complex(+5.067979e+000, +1.074279e+000);// ( 3,3/sin(pi/3)) + BPV[50] = complex(-1.620678e+000, +1.029461e+000);// ( 3,4/sin(pi/3)) + BPV[51] = complex(+1.055970e+000, -2.041230e-001);// ( 3,5/sin(pi/3)) + BPV[52] = complex(+1.619267e+002, +0.000000e+000);// ( 4,0) + BPV[53] = complex(-1.021827e+002, +9.434616e+001);// ( 4,1/sin(pi/3)) + BPV[54] = complex(+9.638391e+000, -8.764529e+001);// ( 4,2/sin(pi/3)) + BPV[55] = complex(+2.157904e+001, +3.568852e+001);// ( 4,3/sin(pi/3)) + BPV[56] = complex(-1.346647e+001, -6.665778e+000);// ( 4,4/sin(pi/3)) + BPV[57] = complex(+4.276735e+000, -9.657825e-002);// ( 4,5/sin(pi/3)) + BPV[58] = complex(+1.435819e+003, +0.000000e+000);// ( 5,0) + BPV[59] = complex(-1.099383e+003, +5.897525e+002);// ( 5,1/sin(pi/3)) + BPV[60] = complex(+4.709887e+002, -6.717565e+002);// ( 5,2/sin(pi/3)) + BPV[61] = complex(-7.965464e+001, +4.040825e+002);// ( 5,3/sin(pi/3)) + BPV[62] = complex(-2.418965e+001, -1.582957e+002);// ( 5,4/sin(pi/3)) + BPV[63] = complex(+1.572560e+004, +0.000000e+000);// ( 6,0) + BPV[64] = complex(-1.334470e+004, +3.525428e+003);// ( 6,1/sin(pi/3)) + BPV[65] = complex(+8.229089e+003, -4.446372e+003);// ( 6,2/sin(pi/3)) + BPV[66] = complex(-3.795705e+003, +3.141147e+003);// ( 6,3/sin(pi/3)) + BPV[67] = complex(+2.095527e+005, +0.000000e+000);// ( 7,0) + BPV[68] = complex(-1.853403e+005, +7.452520e+003);// ( 7,1/sin(pi/3)) + BPV[69] = complex(+1.286235e+005, -8.069218e+003);// ( 7,2/sin(pi/3)) + + //////////////////////////////////////// + // Validate input - kind & derivative // + //////////////////////////////////////// + if ((kind != AIRY) && + (kind != BAIRY) && + (kind != WONE) && + (kind != DWONE) && + (kind != WTWO) && + (kind != DWTWO)) { + return complex(0, 0); // Airy Error: Invalid kind value + }; + + if ((scaling != NONE) && (scaling != HUFFORD) && (scaling != WAIT)) { + return complex(0, 0); // Airy Error: Invalid scaling value + }; + + // Set the derivative flag + if (kind == DWTWO || kind == DWONE || kind == AIRYD || kind == BAIRYD) { + derivative = YES; + } + else { + derivative = NO; + }; + + // Now do something productive with numbers... + + // In the original version of this by George Hufford it calculated Ai(Z), Ai'(Z) and the + // Airy functions of the "3rd kind" Wi(1)(Z) and Wi(2)(Z) (see Eqn 38 in NTIA Report 87-219) + // Note: Theta in 87-219 is Z in this program. + // The input switch had three values, here we are going to have four so that the + // Bairy function doesn't feel left out. + + // The following scales the input parameter Z depending on what the user is trying to do. + // If the user is trying to find just the Ai(Z), Ai'(Z), Bi(Z) or Bi'(Z) there is no scaling. + if (kind == AIRY || kind == BAIRY || kind == AIRYD || kind == BAIRYD) { + // For Ai(Z) and Bi(Z) No translation in the complex plane + U = complex(1.0, 0.0); + } + // Note that W1 Wait = Wi(2) Hufford and W2 Wait = Wi(1) Hufford + // So the following inequalities keep this all straight + else if (((kind == DWONE || kind == WONE) && scaling == HUFFORD) + || + ((kind == DWTWO || kind == WTWO) && scaling == WAIT)) { + // This corresponds to Wi(1)(Z) in Eqn 38 Hufford NTIA Report 87-219 + // or Wait W2 + U = complex(cos(2.0*PI / 3.0), sin(2.0*PI / 3.0)); + } + else if (((kind == DWTWO || kind == WTWO) && scaling == HUFFORD) + || + ((kind == DWONE || kind == WONE) && scaling == WAIT)) { + // This corresponds to Wi(2)(Z) in Eqn 38 Hufford NTIA Report 87-219 + // or Wait W1 + U = complex(cos(-2.0*PI / 3.0), sin(-2.0*PI / 3.0)); + }; + + // Translate the input parameter + ZU = Z * U; + + // We will be only calculating for quadrant 1 and 2. If the desired value is in 3 or 4 we + // will have to flip it over after the calculation + reflection = false; + if (ZU.imag() <= 0) { + reflection = true; // reflection = true means Z.imag() <= 0, use reflection formula to get result + ZU = complex(ZU.real(), -ZU.imag()); + }; + + // Begin the calculation to determine if + // a) the shifted Taylor series will be used or + // b) the Asymptotic approximation is used. + // A shifted Taylor series is necessary because the Taylor series is defined with a center of expansion + // at the origin is a poor approximation to the true value of the (B)Airy function. + + // NOTE: The condition for which method is used is dependant on that value of Z and not the + // transformed version of ZU, which is the shifted Z by exp(+-j*2*pi/3). + // For the calculation of Ai(Z) and Bi(Z) Z is not shifted. + // For the calculation of Wi(1)(Z) and Wi(2)(Z) the value of Ai(ZU) is found. + + // Initialize the indexes for the center of expansion + // Note these are used in the if statements below as flags + N = 0; + NQ8 = 0; + + // If Z is small, use Taylor's series at various centers of expansion chosen by George Hufford + // If Z is large, use Asymptotic series NIST DLMF 9.4.5 - 9.4.8 + + // The following inequality is formed from the implicit arguments for the AV[], BV[], BPV[] and APV[] + // The inequality makes sure that the center of expansion for the Taylor series solution is not + // exceeded. + // (ZU.real() >= -6.5) -6.5 is 0.5 is the real value of the center of expansion in the array + // (ZU.real() <= 7.5) 7.5 is 0.5 is the real value of the center of expansion in the array + // (ZU.imag() <= 6.35) 6.35 is 5.5/sin(PI/3) which is 0.5 past 5/sin(PI/3) + if ((ZU.real() >= -6.5) && (ZU.real() <= 7.5) && (ZU.imag() <= 6.35)) { + + // choose center of expansion of the Taylor series + CoERealidx = (int)(ZU.real() + _copysign(0.5, ZU.real())); + CoEImagidx = (int)(sin(PI / 3.0) * (ZU.imag() + 0.5)); // sin(60)*(Z.imag()+0.5) + + N = NQTT[CoERealidx + 6] + CoEImagidx; // N is index of center of expansion + + // Check to see if N is out of bounds + if (N >= 70) { // Stop if the index N reaches the limit of array AV[] which is 70 + //printf("Airy() Error: Z is too large\n"); + return complex(0, 0); + }; + + NQ8 = NQTT[CoERealidx + 7]; // The next real center of expansion or what is know here ... + // as the area of the Taylor's series + + // if Z is inside Taylor's series area, continue. Otherwise, go to asymptotic series + if (N < NQ8) { + + /////////////////////////////////////////// + // Compute the function by Taylor Series // + /////////////////////////////////////////// + + // sum Taylor's series around nearest expansion point + // The arrays AV[] and APV[] are incremented in the complex domain by 1/sin(PI/3) + CoE = complex((double)CoERealidx, (double)CoEImagidx / sin(PI / 3.0)); + + // Translate the input parameter to the new location + ZU = ZU - CoE; + + // Calculate the first term of the Taylor Series + // To do this we need to find the Airy or Bairy function at the center of + // expansion, CoE, that has been precalculated in the arrays above. + if (kind == BAIRY || kind == BAIRYD) { + Ai = BV[N - 1]; // Bi(CoE) + Aip = BPV[N - 1]; // Bi'(CoE) + } + else { // All other cases use the Coe for Ai(z) + Ai = AV[N - 1]; // Ai(CoE) + Aip = APV[N - 1]; // Ai'(CoE) + }; + + // Find the first elements of the Taylor series + // Translation + B1 = Ai; // B1 is first term for function Ai(a) + B3 = B1 * CoE*ZU; // B3 is second term for derivative Ai(a)*a*(z-a) + A[1] = Aip; // A is first term for derivation Ai'(a) + B2 = A[1] * ZU; // B2 is second term for function Ai'(a)(z-a) + A[0] = B2 + B1; // A[0] is the sum of Ai() or Bi() Ai'(a)(z-a) + Ai(a) + A[1] = A[1] + B3; // A[1] is the sum of Ai'() or Bi'() Ai'(a) + Ai(a)*a*(z-a) + + AN = 1.0; + + // Initialize the counter + cnt = 0; + + // compute terms of series and sum until convergence + do { + do { + AN = AN + 1.0; + B3 = B3 * ZU / AN; + A[0] = B3 + A[0]; + B0 = B1; + B1 = B2; + B2 = B3; + B3 = (CoE*B1 + ZU * B0)*ZU / AN; + A[1] = B3 + A[1]; + + } while ((abs(B2) > (0.5E-7 * abs(A[0]))) + || + (abs(B3) > (0.5E-7 * abs(A[1])))); // Has the convergence criteria been met? + cnt++; + + } while (cnt < 3); // require that the loop be executed 3 times + } + + }; // if ((ZU.real() >= -6.5) && (ZU.real() <= 7.5) && (ZU.imag() <= 6.35)) + +// Determine if the data for the center of expansion is exceeded + if (((ZU.real() < -6.5) || (ZU.real() > 7.5) || (ZU.imag() > 6.35)) + || + (N >= NQ8) + || + (Z.real() == 0 && Z.imag() == 0)) { + + + /////////////////////////////////////////////// + // Compute the function by Asymptotic Series // + /////////////////////////////////////////////// + + /////////////////////////////////////////////////////// + // Please see // + // "On the Asymptotic Expansion of Airy's Integral" // + // E. T. Copson, Cambridge Press, November, 1962 // + // for details on why a second series is necessary // + // The equations found in the reference above appear // + // to be what Hufford used in the creation of this // + // algorithm // + /////////////////////////////////////////////////////// + + // Find intermediate values + ZA = sqrt(ZU); // zeta^(1/2) + ZT = (2.0 / 3.0)*ZU*ZA; // NIST DLMF 9.7.1 => -(2/3)zeta^(3/2) + + if (kind == BAIRY || kind == BAIRYD) { + one = 1.0; // Terms for the Bairy sum do not alternate sign + } + else { + one = -1.0; // All other functions use the Airy whose sum alternates sign + } + + // Compute the asymptotic series either sum over k for (u_k*zeta^-k) or sum over k for (v_k*zeta^-k) + // Which is used depends on M => M = 0 use u_k M = 1 use v_k + // By doing this backward you don't have to do multiple powers zeta^-1 + // Note the coefficients are backward so the for loop will be forward + sum1 = complex(0.0, 0.0); // Initialize the sum + for (int i = 0; i < 14; i++) { + sum1 = (pow(one, i)*ASV[i][derivative] + sum1) / ZT; + }; + sum1 = ASV[SIZE_OF_ASV - 1][derivative] + sum1; // Add the first element that is a function of zeta^0 + + // Now determine if a second series is necessary + // If it is not set the second sum to zero + + //////////////////////////////////////////////////////////////////////////////////// + // Historic Note: + // Hufford originally used the following inequality in AIRY() + // (See OT/ITS RR 11) IF(XT(2) .GT. 0. .AND. XT(l) .LT. ll.8595) LG=4 + // to determine if a second series is required for a reasonable level of accuracy. + // The C translation for the variables defined here + // would be if(ZT.imag() <= 0.0) && (ZT.real() >= -ll.8595) + // In the LFMF code the inequality for the same purpose is (translated to C) + // if((ZT.imag() <= 0.0) && (ZT.real() >= -8.4056)) + // Since ZT = (2/3)ZU^(3/2) and for ZT.imag() = 0.0 and ZT.real() = -8.4056 + // ZU = -2.70859033 + j*4.69141606 which has an angle of 2*PI/3 + // Similarly for Hufford's original code + // ZU = -3.40728475 + j*5.90159031 which has an angle of 2*PI/3 + // Thus we could replace the following with the inequality + // if(ZU.arg() > 2.0*PI/3.0) + ////////////////////////////////////////////////////////////////////////////////////// + + // From Copson the F(z) solution is only valid for phase(z) <= PI/3.0 + // While the F(z) + i*G(z) solution is necessary for phase(z) > PI/3.0 + if (abs(arg(ZU)) > PI / 3.0) { + sum2 = complex(0.0, 0.0); // Initialize the second sum + for (int i = 0; i < 14; i++) { + sum2 = (ASV[i][derivative] + sum2) / ZT; + }; + sum2 = ASV[SIZE_OF_ASV - 1][derivative] + sum2; // Add the first element that is a function of zeta^0 + } + else { // Only one series is necessary for accuracy + sum2 = complex(0.0, 0.0); + }; + + // Now do the final function that leads the sum depending on what the user wants. + // The leading function has to be taken apart so that it can be assembled as necessary for + // the possible two parts of the sum + if (kind == BAIRY || kind == BAIRYD) { + if (derivative == NO) { + ZB = 1.0 / (sqrt(ZA)); // NIST DLMF 9.7.8 + } + else if (derivative == YES) { + ZB = sqrt(ZA); // NIST DLMF 9.7.7 + }; + ZB1 = ZB * exp(ZT) / sqrt(PI); // For Bairy multiply by e^(zeta)/sqrt(PI) + ZB2 = ZB * 1.0 / (exp(ZT)*sqrt(PI)); + + } + else { // All other kind use Airy + if (derivative == NO) { + ZB = 1.0 / sqrt(ZA); // NIST DLMF 9.7.6 + } + else if (derivative == YES) { + ZB = -1.0*sqrt(ZA); // NIST DLMF 9.7.5 + }; + ZB1 = ZB * 1.0 / (2.0*exp(ZT)*sqrt(PI)); // For Airy multiply be e^(-zeta)/(2.0*sqrt(PI)) + ZB2 = ZB * exp(ZT) / (2.0*sqrt(PI)); + }; + + + // Multiply by the leading coefficient to get the results for NIST DLMF 9.7.5 - 9.7.8 + if (derivative == YES) { + A[derivative] = ZB1 * sum1 - complex(0.0, 1.0)*ZB2*sum2; + } + else if (derivative == NO) { + A[derivative] = ZB1 * sum1 + complex(0.0, 1.0)*ZB2*sum2; + }; + + }; // if (( Z.real() < 6.5 || Z.real() > 7.5 || Z.imag() > 6.35 || N > NQ8 || ((Z.real() == 0 && Z.imag() == 0)))) + + + ////////////////////////////////////////////// + // End of the Asymptotic Series Calculation // + ////////////////////////////////////////////// + + // Store the desired quantity + Ai = A[derivative]; + + // Final Transform to get the desired function + // Was the input parameter in quadrant 3 or 4? + // If it was we have to take the conjugate of the calculation result + if (reflection != false) { + Ai = complex(Ai.real(), -Ai.imag()); + }; + + // The final scaling factor is a function of the kind, derivative and scaling flags + if (scaling == NONE) { + // The number from the Taylor series or asymptotic calculation + // does not need to multiplied by anything + U = complex(1.0, 0.0); + } + // Hufford Wi(1) and Wi'(1) + else if ((kind == WONE || kind == DWONE) && (scaling == HUFFORD)) { + if (derivative == NO) { + U = 2.0*complex(cos(-PI / 3.0), sin(-PI / 3.0)); + } + else if (derivative == YES) { + U = 2.0*complex(cos(PI / 3.0), sin(PI / 3.0)); + }; + } + // Hufford Wi(2) and Wi'(2) + else if ((kind == WTWO || kind == DWTWO) && (scaling == HUFFORD)) { + if (derivative == NO) { + U = 2.0*complex(cos(PI / 3.0), sin(PI / 3.0)); + } + else if (derivative == YES) { + U = 2.0*complex(cos(-PI / 3.0), sin(-PI / 3.0)); + }; + } + // Wait W1 and W1' + else if ((kind == WONE || kind == DWONE) && (scaling == WAIT)) { + if (derivative == NO) { + U = complex(sqrt(3.0*PI), -1.0*sqrt(PI)); + } + else if (derivative == YES) { + U = complex(-1.0*sqrt(3.0*PI), -1.0*sqrt(PI)); + }; + } + // Wait W2 and W2' + else if ((kind == WTWO || kind == DWTWO) && (scaling == WAIT)) { + if (derivative == NO) { + U = complex(sqrt(3.0*PI), sqrt(PI)); + } + else if (derivative == YES) { + U = complex(-1.0*sqrt(3.0*PI), sqrt(PI)); + }; + }; + + // Scale the return value + Ai = Ai * U; + + return Ai; + +}; \ No newline at end of file diff --git a/src/FlatEarthCurveCorrection.cpp b/src/FlatEarthCurveCorrection.cpp new file mode 100644 index 0000000..3e6a94c --- /dev/null +++ b/src/FlatEarthCurveCorrection.cpp @@ -0,0 +1,56 @@ +#include "..\include\LFMF.h" + +/*============================================================================= + | + | Description: Calculates the groundwave field strength using the flat Earth + | approximation with curvature correction. + | + | References: 99-368 "Medium Frequency Propagation + | Prediction Techniques and Antenna Modeling for + | Intelligent Transportation Systems (ITS) Broadcast + | Applications", Nicholas DeMinco. Eq (31) + | J. Wait, "Radiation From a Vertical Antenna Over a Curved + | Stratified Ground", Journal of Research of the National + | Bureau of Standards Vol 56, No. 4, April 1956 + | Research Paper 2671 + | + | Input: Delta - Surface impedance + | q - Intermediate value -j*nu*delta + | h_1__km - Height of the higher antenna, in km + | h_2__km - Height of the lower antenna, in km + | d__km - Path distance, in km + | k - Wavenumber, in rad/km + | + | Returns: E_gw - Normalized field strength in mV/m + | + *===========================================================================*/ +double FlatEarthCurveCorrection(complex Delta, complex q, double h_1__km, double h_2__km, double d__km, double k) +{ + complex j = complex(0.0, 1.0); + + // In order for the werf() function to be used both here and in gwfe() + // the argument, qi, has to be defined correctly. The following is how + // it is done in the original GWFEC.FOR + complex qi = (-0.5 + j * 0.5)*sqrt(k*d__km)*Delta; + complex p = qi * qi; + + complex q3 = pow(q, 3); + complex q6 = pow(q, 6); + + // Find F(p) Eqn (32) NTIA Report 99-368 + complex Fofp = 1.0 + sqrt(PI)*j*qi*werf(qi); + + // Calculate f(x) which is the normalized electric field, E_ratio; Eqn (31) NTIA Report 99-368 + complex fofx = Fofp + (1.0 - j * sqrt(PI * p) - (1.0 + 2.0 * p) * Fofp) / (4.0 * q3) + (1.0 - j * sqrt(PI * p) * (1.0 - p) - 2.0 * p + 5.0 * p * p / 6.0 + (p * p / 2.0 - 1.0) * Fofp) / (4.0 * q6); + + // Now find the final normalized field strength from f(x) and the height-gain function for each antenna + // A height-gain function for an antenna is expressed as two terms of a Taylor series + // (See DeMinco NTIA Report 99-368 Aug 1999 + // "Medium Frequency Propagation Prediction Techniques and + // Antenna Modeling for Intelligent Transportation Systems (ITS) Broadcast Applications" + // Equation 36) + + double E_gw__mVm = abs(fofx*(1.0 + j * k*h_2__km*Delta)*(1.0 + j * k*h_1__km*Delta)); + + return E_gw__mVm; +}; \ No newline at end of file diff --git a/src/LFMF.cpp b/src/LFMF.cpp new file mode 100644 index 0000000..d76e6ce --- /dev/null +++ b/src/LFMF.cpp @@ -0,0 +1,105 @@ +#include "..\include\LFMF.h" + +/*============================================================================= + | + | Description: Compute the LFMF propagation prediction + | + | Input: h_tx__meter - Height of the transmitter, in meter + | h_rx__meter - Height of the receiver, in meter + | f__mhz - Frequency, in MHz + | P_tx__watt - Transmitter power, in Watts + | N_s - Surface refractivity, in N-Units + | d__km - Path distance, in km + | epsilon - Relative permittivity + | sigma - Conductivity + | pol - Polarization + | + 0 : POLARIZATION__HORIZONTAL + | + 1 : POLARIZATION__VERTICAL + | + | Outputs: result - Result structure + | + | Returns: error - Error code + | + *===========================================================================*/ +int LFMF(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol, Result *result) +{ + int rtn = ValidateInput(h_tx__meter, h_rx__meter, f__mhz, P_tx__watt, N_s, + d__km, epsilon, sigma, pol); + if (rtn != SUCCESS) + return rtn; + + // Create the complex value j since this was written by electrical engineers + complex j = complex(0.0, 1.0); + + double f__hz = f__mhz * 1e6; + double lambda__meter = C / f__hz; // wavelength, in meters + + double h_1__km = MIN(h_tx__meter, h_rx__meter) / 1000; // lower antenna, in km + double h_2__km = MAX(h_tx__meter, h_rx__meter) / 1000; // higher antenna, in km + + double a_e__km = a_0__km * 1 / (1 - 0.04665 * exp(0.005577 * N_s)); // effective earth radius, in km + + double theta__rad = d__km / a_e__km; + + double k = 2.0 * PI / (lambda__meter / 1000); // wave number, in rad/km + + double nu = pow(a_e__km * k / 2.0, THIRD); // Intermediate value nu + + // dielectric ground constant. See Eq (17) DeMinco 99-368 + complex eta = complex(epsilon, -sigma / (epsilon_0 * 2 * PI * f__hz)); + + // Find the surface impedance, DeMinco 99-368 Eqn (15) + complex delta = sqrt(eta - 1.0); + if (pol == POLARIZATION__VERTICAL) + delta /= eta; + + complex q = -nu * j * delta; // intermediate value q + + // Determine which smooth earth method is used; SG3 Groundwave Handbook, Eq 15 + double d_test__km = 80 * pow(f__mhz, -THIRD); + + double E_gw; + if (d__km < d_test__km) + { + E_gw = FlatEarthCurveCorrection(delta, q, h_1__km, h_2__km, d__km, k); + result->method = METHOD__FLAT_EARTH_CURVE; + } + else + { + E_gw = ResidueSeries(d__km, k, h_1__km, h_2__km, nu, theta__rad, q); + result->method = METHOD__RESIUDE_SERIES; + } + + // Antenna gains + double G_tx__dbi = 4.77; + double G_rx__dbi = 4.77; + + double G_tx = pow(10, G_tx__dbi / 10); + + // Un-normalize the electric field strength + double E_0 = sqrt(ETA * (P_tx__watt * G_tx) / (4.0 * PI)) / d__km; // V/km or mV/m + E_gw = E_gw * E_0; + + // Calculate the basic transmission loss using (derived using Friis Transmission Equation with Electric Field Strength) + // Pt Gt * Pt * ETA * 4*PI * f^2 + // L = ---- = --------------------------- and convert to dB + // Pr E^2 * c^2 + // with all values entered using base units: Watts, Hz, and V/m + // basic transmission loss is not a function of power/gain, but since electric field strength E_gw is a function of (Gt * Pt), + // and Lbtl is a function of 1/E_gw, we add in (Gt * Pt) to remove its effects + result->A_btl__db = 10 * log10(P_tx__watt * G_tx) + + 10 * log10(ETA * 4 * PI) + + 20 * log10(f__hz) + - 20 * log10(E_gw / 1000) + - 20 * log10(C); + + // the 60 constant comes from converting field strength from mV/m to dB(uV/m) thus 20*log10(1e3) + result->E_dBuVm = 60 + 20 * log10(E_gw); + + // Note power is a function of frequency. 42.8 comes from MHz to hz, power in dBm, and the remainder from + // the collection of constants in the derivation of the below equation. + result->P_rx__dbm = result->E_dBuVm + G_rx__dbi - 20.0*log10(f__hz) + 42.8; + + return SUCCESS; +} \ No newline at end of file diff --git a/src/ResidueSeries.cpp b/src/ResidueSeries.cpp new file mode 100644 index 0000000..eb9c113 --- /dev/null +++ b/src/ResidueSeries.cpp @@ -0,0 +1,83 @@ +#include "..\include\LFMF.h" + +/*============================================================================= + | + | Description: Calculates the groundwave field strength using the + | Residue Series method + | + | Input: d__km - Path distance, in km + | k - Wavenumber, in rad/km + | h_1__km - Height of the higher antenna, in km + | h_2__km - Height of the lower antenna, in km + | nu - Intermediate value, + | pow(a_e__km * k / 2.0, THIRD); + | theta__rad - Angular distance of path, in radians + | q - Intermediate value -j*nu*delta + | + | Returns: E_gw - Normalized field strength in mV/m + | + *===========================================================================*/ +double ResidueSeries(double d__km, double k, double h_1__km, double h_2__km, double nu, double theta__rad, complex q) +{ + complex DW2[200], W2[200]; // dummy variables + complex G; + + complex j = complex(0.0, 1.0); + + complex T[200], W1[200], W[200]; + double yHigh, yLow; + + complex GW = complex(0, 0); // initial ground wave + + // Associated argument for the height-gain function H_1[h_1] + yHigh = k * h_2__km / nu; + + // Associated argument for the height-gain function H_2[h_2] + yLow = k * h_1__km / nu; + + double x = nu * theta__rad; + + for (int i = 0; i < 200; i++) + { + T[i] = WiRoot(i + 1, &DW2[i], q, &W2[i], WONE, WAIT); // find the (i+1)th root of Airy function for given q + W1[i] = Airy(T[i], WONE, WAIT); // Airy function of (i)th root + + if (h_1__km > 0) + { + W[i] = Airy(T[i] - yLow, WONE, WAIT) / W1[i]; //height gain function H_1(h_1) eqn.(22) from NTIA report 99-368 + + if (h_2__km > 0) + W[i] *= Airy(T[i] - yHigh, WONE, WAIT) / W1[i]; //H_1(h_1)*H_1(h_2) + } + else if (h_2__km > 0) + W[i] = Airy(T[i] - yHigh, WONE, WAIT) / W1[i]; + else + W[i] = complex(1, 0); + + // W[i] is the coefficient of the distance factor for the i-th + W[i] /= (T[i] - (q*q)); // H_1(h_1)*H_1(h_2)/(t_i-q^2) eqn.26 from NTIA report 99-368 + G = W[i] * exp(-1.0*j*x*T[i]); // sum of exp(-j*x*t_i)*W[i] eqn.26 from NTIA report 99-368 + GW += G; // sum the series + + if (i != 0) + { + if (((abs((GW*GW).real())) + (abs((GW*GW).imag()))) == 0) // when the ground wave is too small, close to 0 + return 0; // end the loop and output E = 0 + else if (((abs((G / GW).real())) + (abs((G / GW).imag()))) < 0.0005) // when the new G is too small compared to its series sum + { + // when the new G is too small compared to its series sum, it's ok to stop the loop + // because adding small number to a significant big one doesn't affect their sum. + //J1 = i; + break; + } + } + } + + // field strength. complex(sqrt(PI/2)) = sqrt(pi)*e(-j*PI/4) + complex Ew = sqrt(x)*complex(sqrt(PI / 2), -sqrt(PI / 2))*GW; + + double E_gw = abs(Ew); // take the magnitude of the result + + return E_gw; + +}; \ No newline at end of file diff --git a/src/ValidateInputs.cpp b/src/ValidateInputs.cpp new file mode 100644 index 0000000..949751e --- /dev/null +++ b/src/ValidateInputs.cpp @@ -0,0 +1,56 @@ +#include "..\include\LFMF.h" + +/*============================================================================= + | + | Description: Perform input parameter validation + | + | Input: h_tx__meter - Height of the transmitter, in meter + | h_rx__meter - Height of the receiver, in meter + | f__mhz - Frequency, in MHz + | P_tx__watt - Transmitter power, in Watts + | N_s - Surface refractivity, in N-Units + | d__km - Path distance, in km + | epsilon - Relative permittivity + | sigma - Conductivity + | pol - Polarization + | + 0 : POLARIZATION__HORIZONTAL + | + 1 : POLARIZATION__VERTICAL + | + | Outputs: [None] + | + | Returns: error - Error code + | + *===========================================================================*/ +int ValidateInput(double h_tx__meter, double h_rx__meter, double f__mhz, double P_tx__watt, + double N_s, double d__km, double epsilon, double sigma, int pol) +{ + if (h_tx__meter < 0 || h_tx__meter > 50) + return ERROR__TX_TERMINAL_HEIGHT; + + if (h_rx__meter < 0 || h_rx__meter > 50) + return ERROR__RX_TERMINAL_HEIGHT; + + if (f__mhz < 0.01 || f__mhz > 30) + return ERROR__FREQUENCY; + + if (P_tx__watt <= 0) + return ERROR__TX_POWER; + + if (N_s < 250 || N_s > 400) + return ERROR__SURFACE_REFRACTIVITY; + + if (d__km < 0.001 || d__km > 10000) + return ERROR__PATH_DISTANCE; + + if (epsilon < 1) + return ERROR__EPSILON; + + if (sigma <= 0) + return ERROR__SIGMA; + + if (pol != POLARIZATION__HORIZONTAL && + pol != POLARIZATION__VERTICAL) + return ERROR__POLARIZATION; + + return SUCCESS; +} \ No newline at end of file diff --git a/src/WiRoot.cpp b/src/WiRoot.cpp new file mode 100644 index 0000000..21591ca --- /dev/null +++ b/src/WiRoot.cpp @@ -0,0 +1,241 @@ +#include "..\include\LFMF.h" + +/*============================================================================= + | + | Description: This routine finds the roots to the equation + | + | Wi'(ti) - q*Wi(ti) = 0 + | + | The parameter i selects the ith root of the equation. The + | function Wi(ti) is the "Airy function of the third kind" + | as defined by Hufford[1] and Wait. The root is found by + | iteration starting from a real root. + | + | Note: Although roots that are found for W1 (Wait) and + | Wi(2) (Hufford) will be equal, and the roots found for + | W2 (Wait) and Wi(1) (Hufford) will be equal, the return + | values for *Wi and *DWi will not be the same. The input + | parameters for kind and scale are used here as they + | are in Airy() for consistency. + | + | References: "Airy Functions of the third kind" are found in equation 38 + | of NTIA Report 87-219 "A General Theory of Radio + | Propagation through a Stratified Atmosphere", George + | Hufford, July 1987 + | + | Input: i - The ith complex root of + | Wi'(2)(ti) - q*Wi(2)(ti) starting with 1 + | DWi - Derivative of "Airy function of the + | third kind" Wi'(2)(ti) + | q - Intermediate value -j*nu*delta + | Wi - "Airy function of the third kind" Wi(2)(ti) + | kind - Kind of Airy function to use + | scaling - Type of scaling to use + | + | Outputs: DWi - Derivative of "Airy function of the + | third kind" Wi'(2)(ti) + | Wi - "Airy function of the third kind" Wi(2)(ti) + | + | Returns: tw - ith complex root of the "Airy function of + | the third kind" + | + *===========================================================================*/ +complex WiRoot(int i, complex *DWi, complex q, complex *Wi, int kind, int scaling) +{ + complex ph; // Airy root phase + complex ct; // Temp + complex ti; // the ith complex root of Wi'(2)(ti) - q*Wi(2)(ti) = 0 + + complex tw; // Return variable + complex T; // Temp + complex A; // Temp + + double t, tt; // Temp + double eps; // Temp + + int cnt; // Temp + int dkind; // Temp + + // From the NIST DLMF (Digital Library of Mathematical Functions) + // http://dlmf.nist.gov/ + // The first 10 roots of Ai(Z) and Ai'(Z) can be found in: Table 9.9.1: Zeros of Ai and Ai. + // Note: That ak is the root of Ai(ak) while akp (ak') is the root of Ai'(akp) + + // Root of the Airy function, Ai(ak) + // TZERO(I) in GWINT.FOR + double akp[] = { -1.0187929716, + -3.2481975822, + -4.8200992112, + -6.1633073556, + -7.3721772550, + -8.4884867340, + -9.5354490524, + -10.5276603970, + -11.4750666335, + -12.3847883718, + -13.2636395229 }; + + // Root of the derivative of Airy function, Ai'(akp) + // TINFIN(I) in GWINT.FOR + double ak[] = { -2.3381074105, + -4.0879494441, + -5.5205698281, + -6.7867080901, + -7.9441335871, + -9.0226508533, + -10.0401743416, + -11.0085243037, + -11.9360255632, + -12.8287867529, + -13.6914890352 }; + + // Verify that the input data is correct + // Make sure that the desired root is greater than or equal to one + if (i <= 0) + { + // There is an input parameter error; printf("WiRoot(): The root must be >= 0 (%d)\n", i); + tw = complex(-998.8, -998.8); + return tw; + }; + + if ((scaling != HUFFORD) && (scaling != WAIT)) + { + // There is an input parameter error; printf("WiRoot(): The scaling must be HUFFORD (%d) or WAIT (%d)\n", HUFFORD, WAIT); + tw = complex(-997.7, -997.7); + return tw; + }; + + if ((kind != WTWO) && (kind != WONE)) + { + // There is an input parameter error; printf("WiRoot(): The kind must be W1 (%d) or W2 (%d)\n", WONE, WTWO); + tw = complex(-996.6, -996.6); + return tw; + }; + // Input parameters verified + + // Initialize the Wi and Wi'(z)functions + *DWi = complex(0.0, 0.0); // Wi'(z) + *Wi = complex(0.0, 0.0); // Wi(z) + + // This routine starts with a real root of the Airy function to find the complex root + // The real root has to be turned into a complex number. + + // ph is a factor that is used to find the root of the Wi function + // Determine what scaling the user wants and which Wi function is used to set ph and + // the dkind flag. This will allow that the real root that starts this process can be + // scaled appropriately. + // This is the simmilar to the initial scaling that is done in Airy() + // Note that W1 Wait = Wi(2) Hufford and W2 Wait = Wi(1) Hufford + // So the following inequalities keep this all straight + if ((kind == WONE && scaling == HUFFORD) || (kind == WTWO && scaling == WAIT)) + { + // Wi(1)(Z) in Eqn 38 Hufford NTIA Report 87-219 or Wait W2 + ph = complex(cos(-2.0*PI / 3.0), sin(-2.0*PI / 3.0)); + // Set the dkind flag + if (scaling == WAIT) + { + dkind = DWTWO; + } + else + { + dkind = DWONE; + }; + } + else if ((kind == WTWO && scaling == HUFFORD) || (kind == WONE && scaling == WAIT)) + { + // Wi(2)(Z) in Eqn 38 Hufford NTIA Report 87-219 or Wait W1 + ph = complex(cos(2.0*PI / 3.0), sin(2.0*PI / 3.0)); + if (scaling == WAIT) + { + dkind = DWONE; + } + else + { + dkind = DWTWO; + }; + }; + + + // Note: The zeros of the Airy functions i[ak'] and Ak'[ak], ak' and ak, are on the negative real axis. + // This is why 4*i+3 and 4*i+1 are used here instead of 4*k-3 and 4*k-1 which are + // used in 9.9.8 and 9.9.6 in NIST DLMF. We are finding the ith negative root here. + if (pow(abs(q), 3.0) <= 4 * (i - 1) + 3) { + // Small Z, use ak' as the first guess (Ak(ak') = 0) + if (i <= 10) + { + // The desired root is less than 10 so it is in the array above + tt = akp[i - 1]; + } + else + { + // The desired root is a higher order than those given in the ak array above + // so we will approximate it from the first three terms of NIST DLMF 9.9.1.9 + // First find the argument (9.9.8) used in 9.9.1.9 for the ith negative root of Ai'(ak). + t = (3.0 / 8.0)*PI*(4.0*(i - 1) + 1); + tt = -1.0*pow(t, 2.0 / 3.0)*(1.0 - ((7.0 / 48.0)*pow(t, -2.0)) + ((35.0 / 288.0)*pow(t, -4.0))); + }; + // Make the real Airy root complex + ti = tt * ph; + // t is now the solution for q = 0, The next step is the first Newton iteration + ti = ti + q / ti; + } + else + { + // Large q, use ak as the first guess (Ai'(ak) = 0) + if (i <= 10) + { + // The desired root is less than 10 so it is in the array above + tt = ak[i - 1]; + } + else + { + // The desired root must be approximated from the first three terms of NIST DLMF 9.9.1.8 + // First find the argument (9.9.6) used in 9.9.1.8 for the ith negative root of Ai(ak). + t = (3.0 / 8.0)*PI*(4.0*(i - 1) + 3.0); + tt = -1.0*pow(t, 2.0 / 3.0)*(1.0 + ((5.0 / 48.0)*pow(t, -2.0)) - ((5.0 / 36.0)*pow(t, -4.0))); + }; + ti = tt * ph; + // t is now the solution for Z = infinity. Next step the first newton iteration + ti = ti + 1.0 / q; + }; + + cnt = 0; // Set the iteration counter + eps = 0.5e-6; // Set the error desired for the iteration + + // Now iterate by Newton's method + + ////////////////////////////////////////////////////////////////////// + // Note: We can use the following from + // Berry and Christman "The Path Integrals of LF/VLF Wave Hop Theory" + // Radio Science Vol. 69D, No. 11 , November 1965 + // Eqn (14) E(t) = W2'(t) - q W2(t) + // Eqn (39) E'(t) = t W2(t) - q W2'(t) + ////////////////////////////////////////////////////////////////////// + + do { + + // f(q) = Wi'(ti) - q*Wi(ti) + *Wi = Airy(ti, kind, scaling); + // f'(q) = tw*Wi(ti) - q*Wi'(ti); + *DWi = Airy(ti, dkind, scaling); + // The Newton correction factor for iteration f(q)/f'(q) + A = (*DWi - q * (*Wi)) / (ti*(*Wi) - q * (*DWi)); + ti = ti - A; // New root guess ti + cnt++; // Increment the counter + + } while ((cnt <= 25) && ((abs((A / ti).real()) + (abs((A / ti).imag())) > eps))); + + // Check to see if there if the loop converged on an answer + if (cnt == 26) // The cnt that fails is an arbitrary number most converge in ~5 tries + { + // No Convergence return 0 + j*0 as the root as TW() did + tw = complex(0.0, 0.0); + } + else + { + // Converged! + tw = ti; + }; + + return tw; +}; \ No newline at end of file diff --git a/src/werf.cpp b/src/werf.cpp new file mode 100644 index 0000000..73160a6 --- /dev/null +++ b/src/werf.cpp @@ -0,0 +1,83 @@ +#include "..\include\LFMF.h" + +/*============================================================================= + | + | Description: Calculates the complex error function, which is more + | commonly known as the Faddeeva function (a.k.a. w(z) in + | Abramowitz and Stegun Eqn 7.1.3). This routine makes the + | calculation in the first quadrant then transforms to the + | appropriate quadrant indicated by the input parameter, qi. + | For large argument the asymptotic form of w(z) is used + | while for small argument it is solved by using the + | Maclaurin series expansion of the Dawson's integral is used. + | + | Input: qi - Input argument + | + | Outputs: [None] + | + | Returns: w - Complex value of w(qi) ( exp(-qi^2)*erfc(-i*qi) ) + | + *===========================================================================*/ +complex werf(complex qi) +{ + complex za, a, w, jzn; + + complex j = complex(0.0, 1.0); + + // Force qi into the first quadrant + complex z = complex(abs(qi.real()), abs(qi.imag())); + + if (abs(qi) > 2.0) + { + complex zsq = pow(z, 2); + + // ref: + // - https://dlmf.nist.gov/7.9 + // - https://arxiv.org/ftp/arxiv/papers/1711/1711.05291.pdf + w = z * ((0.4613135279 * j) / (zsq - 0.1901635092) + (0.09999216168 * j) / (zsq - 1.7844927485) + (0.0028838938748 * j) / (zsq - 5.52534374379)); + } + else + { + // Determine the Faddeeva function, w(z), using the method given in + // "Computation of the Complex Error Function" J.A.C. Weideman + // SAIM Vol 31, No. 5, pp 1497-1518, October 1994 + // from Eqn (4) in the above: + // w(z) = ... = exp(-x^2) + ((i*2)/sqrt(PI))*daw(x) + complex term = z; + double num = 1.0; + double denom = 1.0; + double sign = 1.0; + complex Csum = z; + + // Determine the Maclaurin series of the Dawson's integral + for (int i = 1; i < WERF_TERMS; i++) { + term = term * z * z; // z^(2n+1) + num = num * 2; // 2^n + denom = denom * (2 * i + 1); // (2n+1)!! + sign = -sign; // (-1)^n + Csum = Csum + sign * (num / denom) * term; + } + + // Given Csum is series expansion of the Dawson's integral + // find the Faddeeva function from: + // w(z) = ... = exp(-x^2) + ((2*i)/sqrt(PI))*daw(x) + w = exp(-1.0 * z * z) + ((2.0 * j) / sqrt(PI)) * Csum; + }; + + // Now put the result in the correct quadrant given by qi + if (qi.imag() < 0.0) + { + // Use 7.1.11 w(-z) = 2*exp(-z^2) - w(z) + // If you are in the lower half plane + w = 2.0 * exp(-1.0 * z * z) - w; + if (qi.real() > 0.0) + w = conj(w); + } + else + { + if (qi.real() < 0.0) + w = conj(w); + }; + + return w; +}; \ No newline at end of file diff --git a/win32/LFMF.def b/win32/LFMF.def new file mode 100644 index 0000000..5c296c9 --- /dev/null +++ b/win32/LFMF.def @@ -0,0 +1,3 @@ +LIBRARY +EXPORTS + LFMF \ No newline at end of file diff --git a/win32/LFMF.rc b/win32/LFMF.rc new file mode 100644 index 0000000..9d7b895 --- /dev/null +++ b/win32/LFMF.rc @@ -0,0 +1,99 @@ +// Microsoft Visual C++ generated resource script. +// +#include "resource.h" + +#define APSTUDIO_READONLY_SYMBOLS +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 2 resource. +// +#include "winres.h" + +///////////////////////////////////////////////////////////////////////////// +#undef APSTUDIO_READONLY_SYMBOLS + +///////////////////////////////////////////////////////////////////////////// +// English (United States) resources + +#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU) +LANGUAGE LANG_ENGLISH, SUBLANG_ENGLISH_US +#pragma code_page(1252) + +#ifdef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// TEXTINCLUDE +// + +1 TEXTINCLUDE +BEGIN + "resource.h\0" +END + +2 TEXTINCLUDE +BEGIN + "#include ""winres.h""\r\n" + "\0" +END + +3 TEXTINCLUDE +BEGIN + "\r\n" + "\0" +END + +#endif // APSTUDIO_INVOKED + + +///////////////////////////////////////////////////////////////////////////// +// +// Version +// + +VS_VERSION_INFO VERSIONINFO + FILEVERSION 1,0,0,0 + PRODUCTVERSION 1,0,0,0 + FILEFLAGSMASK 0x3fL +#ifdef _DEBUG + FILEFLAGS 0x1L +#else + FILEFLAGS 0x0L +#endif + FILEOS 0x40004L + FILETYPE 0x2L + FILESUBTYPE 0x0L +BEGIN + BLOCK "StringFileInfo" + BEGIN + BLOCK "040904b0" + BEGIN + VALUE "CompanyName", "The Institute for Telecommunication Sciences" + VALUE "FileDescription", "Low Frequency / Medium Frequency (LF/MF) Propagation Model" + VALUE "FileVersion", "1.0.0.0" + VALUE "InternalName", "LFMF.dll" + VALUE "OriginalFilename", "LFMF.dll" + VALUE "ProductName", "Low Frequency / Medium Frequency (LF/MF) Propagation Model" + VALUE "ProductVersion", "1.0.0.0" + END + END + BLOCK "VarFileInfo" + BEGIN + VALUE "Translation", 0x409, 1200 + END +END + +#endif // English (United States) resources +///////////////////////////////////////////////////////////////////////////// + + + +#ifndef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 3 resource. +// + + +///////////////////////////////////////////////////////////////////////////// +#endif // not APSTUDIO_INVOKED + diff --git a/win32/LFMF.sln b/win32/LFMF.sln new file mode 100644 index 0000000..02f8307 --- /dev/null +++ b/win32/LFMF.sln @@ -0,0 +1,31 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +# Visual Studio Version 16 +VisualStudioVersion = 16.0.30320.27 +MinimumVisualStudioVersion = 10.0.40219.1 +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "LFMF", "LFMF.vcxproj", "{9A9F0547-8206-4C93-A0A6-545F693F6118}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|x64 = Debug|x64 + Debug|x86 = Debug|x86 + Release|x64 = Release|x64 + Release|x86 = Release|x86 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Debug|x64.ActiveCfg = Debug|x64 + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Debug|x64.Build.0 = Debug|x64 + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Debug|x86.ActiveCfg = Debug|Win32 + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Debug|x86.Build.0 = Debug|Win32 + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Release|x64.ActiveCfg = Release|x64 + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Release|x64.Build.0 = Release|x64 + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Release|x86.ActiveCfg = Release|Win32 + {9A9F0547-8206-4C93-A0A6-545F693F6118}.Release|x86.Build.0 = Release|Win32 + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection + GlobalSection(ExtensibilityGlobals) = postSolution + SolutionGuid = {D501E6CA-935A-42E0-816A-F56CD1DE92D5} + EndGlobalSection +EndGlobal diff --git a/win32/LFMF.vcxproj b/win32/LFMF.vcxproj new file mode 100644 index 0000000..c27eaa7 --- /dev/null +++ b/win32/LFMF.vcxproj @@ -0,0 +1,195 @@ + + + + + Debug + Win32 + + + Release + Win32 + + + Debug + x64 + + + Release + x64 + + + + + + + + + + + + + + + + + + + + + + + 16.0 + Win32Proj + {9a9f0547-8206-4c93-a0a6-545f693f6118} + LFMF + 10.0 + + + + DynamicLibrary + true + v142 + Unicode + + + DynamicLibrary + false + v142 + true + Unicode + + + DynamicLibrary + true + v142 + Unicode + + + DynamicLibrary + false + v142 + true + Unicode + + + + + + + + + + + + + + + + + + + + + true + $(SolutionDir)..\bin\x86\Debug + + + false + $(SolutionDir)..\bin\x86\Release + + + true + $(SolutionDir)..\bin\x64\Debug + + + false + $(SolutionDir)..\bin\x64\Release + + + + Level3 + true + WIN32;_DEBUG;LFMF_EXPORTS;_WINDOWS;_USRDLL;%(PreprocessorDefinitions) + true + NotUsing + + + MultiThreadedDebug + StdCall + + + Windows + true + false + LFMF.def + + + + + Level3 + true + true + true + WIN32;NDEBUG;LFMF_EXPORTS;_WINDOWS;_USRDLL;%(PreprocessorDefinitions) + true + NotUsing + + + MultiThreaded + StdCall + + + Windows + true + true + true + false + LFMF.def + + + + + Level3 + true + _DEBUG;LFMF_EXPORTS;_WINDOWS;_USRDLL;%(PreprocessorDefinitions) + true + NotUsing + + + MultiThreadedDebug + StdCall + + + Windows + true + false + LFMF.def + + + + + Level3 + true + true + true + NDEBUG;LFMF_EXPORTS;_WINDOWS;_USRDLL;%(PreprocessorDefinitions) + true + NotUsing + + + MultiThreaded + StdCall + + + Windows + true + true + true + false + LFMF.def + + + + + + \ No newline at end of file diff --git a/win32/resource.h b/win32/resource.h new file mode 100644 index 0000000..7b6a391 --- /dev/null +++ b/win32/resource.h @@ -0,0 +1,14 @@ +//{{NO_DEPENDENCIES}} +// Microsoft Visual C++ generated include file. +// Used by LFMF.rc + +// Next default values for new objects +// +#ifdef APSTUDIO_INVOKED +#ifndef APSTUDIO_READONLY_SYMBOLS +#define _APS_NEXT_RESOURCE_VALUE 101 +#define _APS_NEXT_COMMAND_VALUE 40001 +#define _APS_NEXT_CONTROL_VALUE 1001 +#define _APS_NEXT_SYMED_VALUE 101 +#endif +#endif