Add functions

CMakeLists.txt

@@ -0,0 +1,18 @@
+cmake_minimum_required(VERSION 3.5)
+
+project(
+    Elleptic_integrals
+    LANGUAGES CXX
+)
+
+set(CMAKE_CXX_STANDARD 20)
+
+# enable_testing()
+add_subdirectory(googletest)
+add_subdirectory(modules)
+add_subdirectory(bin)
+# add_executable(test_app main.cpp)
+
+enable_testing()
+add_subdirectory(TEST)
+

Test/CMakeLists.txt

@@ -0,0 +1,28 @@
+# include(FetchContent)
+
+# FetchContent_Declare(
+#   googletest
+#   GIT_REPOSITORY https://github.com/google/googletest.git
+#   GIT_TAG        release-1.12.1
+# )
+
+# set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
+# FetchContent_MakeAvailable(googletest)
+
+enable_testing()
+
+add_executable(
+  Elleptic_integralsTest
+  Elleptic_integralsTest.cpp
+)
+
+target_link_libraries(
+  Elleptic_integralsTest
+  GTest::gtest_main
+)
+
+target_include_directories(Elleptic_integralsTest PUBLIC ${PROJECT_SOURCE_DIR})
+
+include(GoogleTest)
+
+gtest_discover_tests(Elleptic_integralsTest)

Test/Elleptic_integralsTest.cpp

@@ -0,0 +1,76 @@
+#include <gtest/gtest.h>
+#include "modules/Elliptic_integrals.cpp"
+#include <cmath>
+
+class Elliptic_IntegralsTest : public testing::Test {
+    protected:
+
+    Elliptic_IntegralsTest() 
+    {
+
+    }
+
+    ~Elliptic_IntegralsTest() override 
+    {
+
+    }
+
+    void SetUp() override 
+    {
+    }
+
+    void TearDown() override 
+    {
+    }
+
+};
+
+TEST_F(Elliptic_IntegralsTest, Elliptic_0_0)
+{
+    Elliptic_Integrals ElInt;
+
+    std::vector<double> res = ElInt.Elliptic(0, 0);
+    EXPECT_DOUBLE_EQ(res[0], 0) << "F(0, 0)";
+    EXPECT_DOUBLE_EQ(res[1], 0) << "E(0, 0)";
+}
+
+TEST_F(Elliptic_IntegralsTest, Elliptic_90_0)
+{
+    Elliptic_Integrals ElInt;
+
+    std::vector<double> res = ElInt.Elliptic(90, 0);
+    EXPECT_DOUBLE_EQ(res[0], M_PI / 2) << "F(90, 0)";
+    EXPECT_DOUBLE_EQ(res[1], M_PI / 2)<< "E(90, 0)";
+}
+
+TEST_F(Elliptic_IntegralsTest, Elliptic_0_1)
+{
+    Elliptic_Integrals ElInt;
+
+    std::vector<double> res = ElInt.Elliptic(0, 1);
+    EXPECT_DOUBLE_EQ(res[0], 0) << "F(0, 1)";
+    EXPECT_DOUBLE_EQ(res[1], 0) << "E(0, 1)";
+}
+
+TEST_F(Elliptic_IntegralsTest, Elliptic_90_1)
+{
+    Elliptic_Integrals ElInt;
+
+    std::vector<double> res = ElInt.Elliptic(90, 1);
+    EXPECT_DOUBLE_EQ(res[0], std::numeric_limits<double>::infinity()) << "F(90, 1)";
+    EXPECT_DOUBLE_EQ(res[1], 1)<< "E(90, 1)";
+}
+
+TEST_F(Elliptic_IntegralsTest, Elliptic_Throw) 
+{
+    Elliptic_Integrals ellIntegrals;
+    std::vector<double> result = ellIntegrals.Elliptic(50, 2); 
+
+    EXPECT_EQ(result.size(), 0); // if k > 1 or k < - 1 vector of Elliptic(phi, k) should be empty
+}
+
+int main(int argc, char **argv) 
+{
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}

bin/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_executable(${PROJECT_NAME} main.cpp)
+
+
+target_link_libraries(${PROJECT_NAME} PRIVATE Calc)
+target_include_directories(${PROJECT_NAME} PUBLIC ${PROJECT_SOURCE_DIR})

bin/main.cpp

@@ -0,0 +1,94 @@
+#include "modules/Elliptic_integrals.h"
+#include <iomanip>
+#include <string>
+
+
+void PrintEF(int index)
+{
+  Elliptic_Integrals EL;
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;
+  std::cout << std::setw(2) << "alpha/phi";
+  for (int j=0; j <=6; j++)
+  {
+    std::cout << std::setw(11) << std::right << std::setprecision(6) << j * 5;
+  }
+  std::cout << std::endl;
+  // std::cout << "      ";
+  for (double alpha = 0; alpha <= 90; alpha = alpha + 2)
+  {
+    double k = sin(alpha * M_PI / 180);
+    if ((int) alpha % 10 == 0 && alpha != 0)
+    {
+      std::cout << std::endl;
+    }
+    std::cout  << std::setw(9) << int(alpha) << ' ';
+    for (double phi = 0; phi <= 30; phi = phi + 5)
+    {
+      std::vector<double> res = EL.Elliptic(phi, k);
+      std::cout << std::fixed << std::setprecision(8) << res[index] << " ";
+    }
+    std::cout << std::endl;
+  }
+
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;
+
+  std::cout << std::setw(2) << "alpha/phi";
+  for (int j=7; j <=12; j++)
+  {
+    std::cout << std::setw(11) << std::right << std::setprecision(6) << j * 5;
+  }
+  std::cout << std::endl;
+  // std::cout << "      ";
+  for (double alpha = 0; alpha <= 90; alpha = alpha + 2)
+  {
+    double k = sin(alpha * M_PI / 180);
+    if ((int) alpha % 10 == 0 && alpha != 0)
+    {
+      std::cout << std::endl;
+    }
+    std::cout  << std::setw(9) << int(alpha) << ' ';
+    for (double phi = 35; phi <= 60; phi = phi + 5)
+    {
+      std::vector<double> res = EL.Elliptic(phi, k);
+      std::cout << std::fixed << std::setprecision(8) << res[index] << " ";
+    }
+    std::cout << std::endl;
+  }
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;
+
+  std::cout << std::setw(2) << "alpha/phi";
+  for (int j=13; j <=18; j++)
+  {
+    std::cout << std::setw(11) << std::right << std::setprecision(6) << j * 5;
+  }
+  std::cout << std::endl;
+  // std::cout << "      ";
+  for (double alpha = 0; alpha <= 90; alpha = alpha + 2)
+  {
+    double k = sin(alpha * M_PI / 180);
+    if ((int) alpha % 10 == 0 && alpha != 0)
+    {
+      std::cout << std::endl;
+    }
+    std::cout  << std::setw(9) << int(alpha) << ' ';
+    for (double phi = 65; phi <= 90; phi = phi + 5)
+    {
+      std::vector<double> res = EL.Elliptic(phi, k);
+      std::cout << std::fixed << std::setprecision(8) << res[index] << " ";
+    }
+    std::cout << std::endl;
+  }
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;
+}
+
+int main()
+{
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;
+  std::cout << "Elliptic integral of the first kind F(phi, sin(alpha))" << std::endl;
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;  
+  PrintEF(0);
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;
+  std::cout << "Elliptic integral of the second kind E(phi, sin(alpha))" << std::endl;
+  std::cout << "------------------------------------------------------------------------------------------" << std::endl;  
+  PrintEF(1);
+}

googletest

@@ -0,0 +1 @@
+Subproject commit f8d7d77c06936315286eb55f8de22cd23c188571

modules/CMakeLists.txt

@@ -0,0 +1 @@
+add_library(Calc Elliptic_integrals.cpp Elliptic_integrals.h)

modules/Elliptic_integrals.cpp

@@ -0,0 +1,160 @@
+#include "Elliptic_integrals.h"
+
+// std::string Error()
+// {
+//     std::string mes = "Error! k should be from -1 to 1";
+//     std::exception(mes);
+// }
+
+/**
+ * Purpose: Constructor
+ * 
+ * @return None
+ * @throws None
+ */
+
+Elliptic_Integrals::Elliptic_Integrals()
+{
+
+}
+
+/**
+ * Purpose: Destructor
+ * 
+ * @return None
+ * @throws None
+ */
+
+Elliptic_Integrals::~Elliptic_Integrals()
+{
+    
+}
+
+/**
+ * Purpose: Compute complete elliptic integrals K(k) and E(k)
+ * 
+ * @param k - modulus k (-1 <= k <= 1)
+ * @return Vector with values of K(k) (the first place) and E(k) (the second place)
+ * @throws "Error! k should be from -1 to 1"
+ */
+
+std::vector<double> Elliptic_Integrals::CompleteElliptic(const double k)
+{
+    std::vector<double> result;
+
+    double PK = 1.0 - k * k;
+
+    if ((k > 1) || (k < -1))
+    {
+        std::cerr << "Parameter k must be in the range [-1, 1]" << std::endl;
+        return result;
+    }
+    
+    if (k == 1.0)
+    {
+        double CK = std::numeric_limits<double>::infinity();
+        double CE = 1.0;
+        result.push_back(CK);
+        result.push_back(CE);
+    }
+    else
+    {
+        double AK = (((0.01451196212 * PK + 0.03742563713) * PK + 0.03590092383) * PK + 0.09666344259) * PK + 1.38629436112;
+        double BK = (((0.00441787012 * PK + 0.03328355346) * PK + 0.06880248576) * PK + 0.12498593597) * PK + 0.5;
+        double CK = AK - BK * log(PK);
+        double AE = (((0.01736506451 * PK + 0.04757383546) * PK + 0.06260601220) * PK + 0.44325141463) * PK + 1.0;
+        double BE = (((0.00526449639 * PK + 0.04069697526) * PK + 0.09200180037) * PK + 0.24998368310) * PK;
+        double CE = AE - BE * log(PK);
+        result.push_back(CK);
+        result.push_back(CE);
+    }
+    return result;
+}
+
+
+/**
+ * Purpose: Compute elliptic integrals F(phi, k) and E(phi, k)
+ * 
+ * @param k - modulus k (-1 <= k <= 1)
+ * @param phi - argument (in degrees)
+ * @return Vector with values of F(phi, k)) (the first place) and E(phi, k) (the second place)
+ * @throws "Error! k should be from -1 to 1"
+ */
+
+std::vector<double> Elliptic_Integrals::Elliptic(const double phi, const double k)
+{
+    std::vector<double> result;
+
+    if ((k > 1) || (k < -1))
+    {
+        std::cerr << "Parameter k must be in the range [-1, 1]" << std::endl;
+        return result;
+    }
+
+    double G = 0.0;
+    double pi = M_PI; 
+    double A0 = 1.0;
+    double B0 = sqrt(1.0 - k * k);
+    double D0 = pi / 180.0 * phi;
+    double R = k * k;
+
+    if (k == 1 and phi == 90)
+    {
+        double FE = std::numeric_limits<double>::infinity();
+        double EE = 1.0;
+        result.push_back(FE);
+        result.push_back(EE);
+    }
+    else if (k == 1)
+    {
+        double FE = log((1.0 + sin(D0)) / cos(D0));
+        double EE = sin(D0);
+        result.push_back(FE);
+        result.push_back(EE);
+    }
+    else
+    {
+        double FAC = 1.0;
+        double A;
+        double B;
+        double C;
+        double D;
+        for (int k = 1; k <= 1000; k++)
+        {
+            A = (A0 + B0) / 2.0;
+            B = sqrt(A0 * B0);
+            C = (A0 - B0) / 2.0;
+            FAC = 2.0 * FAC;
+            R = R + FAC * C * C;
+            if (phi != 90)
+            {
+                D = D0 + atan((B0 / A0) * tan(D0));
+                G = G + C * sin(D);
+                D0 = D + pi * int(D / pi + 0.5);
+            }
+            A0 = A;
+            B0 = B;
+            if (C < 1e-7) goto stop;
+            continue;
+        }
+        stop:
+        double CK = pi / (2.0 * A);
+        double CE = pi * (2.0 - R) / (4.0 * A);
+        if (phi == 90)
+        {
+            double FE = CK;
+            double EE = CE;
+            result.push_back(FE);
+            result.push_back(EE);
+        }
+        else
+        {
+            double FE = D / (FAC * A);
+            double EE = FE * CE / CK + G;
+            result.push_back(FE);
+            result.push_back(EE);
+        }
+    }
+    return result;
+}
+

modules/Elliptic_integrals.h

@@ -0,0 +1,17 @@
+#define _USE_MATH_DEFINES 
+
+#pragma once
+#include <iostream> 
+#include <string>
+#include <vector>
+#include <cmath>
+
+class Elliptic_Integrals
+{
+public:
+    Elliptic_Integrals();
+    ~Elliptic_Integrals();
+    // Error(const char* err);
+    std::vector<double> CompleteElliptic(const double k);
+    std::vector<double> Elliptic(const double k, const double phi);
+};