#include "modules/Elliptic_integrals.h"
#include <iomanip>
#include <string>


void PrintEF(int index)
{
  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;

  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 = 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;

  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 = 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;

  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 = Elliptic(phi, k);
      std::cout << std::fixed << std::setprecision(8) << res[index] << " ";
    }
    std::cout << std::endl;
  }
  std::cout << "------------------------------------------------------------------------------------------" << std::endl;
}

/**
 * Purpose: A table for comparison with Irene A. Stegun, Milton Abramowitz, "Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables", page 613-618
 * 
 * @return A table for comparison
 */

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