© 2021 Neveles. All rights reserved.

#include <iostream>
#include <math.h>
#include <iomanip>

#include "cmath.h"
#include "quanc8.c"
#include "zeroin.c"
#include "rkf45.c"

double const u = 1 / 82.3;

double const A = 1.2;
double const B = 0;
double const C = 0;
double D = 0;

double const T = 12;
double const H = 0.4; // лучше брать меньше

double const ERR_FOR_QUANC = 1.0e-10;
double const TOL_FOR_ZEROIN = 1.0e-10;
double const ERR_FOR_RKF = 1.0e-4; // Больше брать нельзя
double const MAXNFE = 30000;

double funForQuanc8(double z)
{
  return 1 / (exp(1.8 * z * z) + D);
}

double quanc(double x)
{
  double a = 0.0;
  double b = 1.0;
  double abserr = ERR_FOR_QUANC;
  double relerr = ERR_FOR_QUANC;
  double result = 0.0;
  double errest;
  int nofun;
  double posn;
  int flag;

  double funForQuanc8(double z);

  D = x;
  quanc8(funForQuanc8, a, b, abserr, relerr, &result, &errest, &nofun, &posn, &flag);

  return result - 0.3644632 * x;
}

double calculateDWithZeroin()
{
  double d1 = 0.0;
  double d2 = 3.0;
  double tol = TOL_FOR_ZEROIN;
  int flag;

  double quanc(double D);

  double resultD = zeroin(d1, d2, quanc, tol, &flag);

  std::cout.precision(8);
  std::cout << "\nTOL = " << TOL_FOR_ZEROIN << '\n';
  std::cout << "D   = " << std::setw(9) << D;
  std::cout << "\nERR for rkf = " << ERR_FOR_RKF;
  std::cout << "\nMAXNFE      = " << MAXNFE;
  std::cout << "\n\n";

  return resultD;
}

int funForRKF(int n, double t, double* x, double* dxdt)
{
  dxdt[0] = x[1]; // m1' = m2
  dxdt[2] = x[3]; // m3' = m4

  dxdt[1] = -(1 - u) * x[0] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
    (x[0] - 1) / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));

  dxdt[3] = -(1 - u) * x[2] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
    x[2] / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));

  return 0;
}

void RKF45(double* x0array, double* x1array, double* x2array, double* x3array)
{
  int n = 4;

  double x[4];
  double xp[4];

  double t = 0.0;
  double tout = 0.0;

  double relerr = ERR_FOR_RKF;
  double abserr = ERR_FOR_RKF;

  double h = 0.0;

  int nfe = 0;
  int maxnfe = 5000;

  int flag = 1;

  int step = 0;

  double D = calculateDWithZeroin();

  int  fail = 0;
  rkfinit(n, &fail);

  if (fail != 0)
  {
    return;
  }

  x[0] = A;
  x[1] = B;
  x[2] = C;
  x[3] = D;

  std::cout << "   t        x          x'         y          y'\n";
  std::cout << "------------------------------------------------------\n";
  for (step = 1; step <= T / H; ++step)
  {
    tout = H * step;
    t = tout - H;

    int funForRKF(int n, double t, double* x, double* dxdt);

    rkf45(funForRKF, n, x, xp, &t, tout, &relerr, abserr, &h, &nfe, maxnfe, &flag);

    printf("%5.1f %10.6f %10.6f %10.6f %10.6f \n", t, x[0], x[1], x[2], x[3]);

    if (flag != 2)
    {
      std::cout << "ERROR IN FLAG (flag = " << flag << ")\n";

      break;
    }
  }

  std::cout << "------------------------------------------------------\n";


  rkfend();

  std::cout << "nfe  : " << nfe << '\n';
  std::cout << "h    : " << h << '\n';
  std::cout << "flag : " << flag << '\n';

  std::cout << "------------------------------------------------------\n\n";
}

int main()
{
  double rkf45x0Array[30] = { 0 };
  double rkf45x1Array[30] = { 0 };
  double rkf45x2Array[30] = { 0 };
  double rkf45x3Array[30] = { 0 };

  RKF45(rkf45x0Array, rkf45x1Array, rkf45x2Array, rkf45x3Array);

  return 0;
}