Спасите!!!

#include "pch.h"
#include <iostream>
#include <conio.h>
#include <string>
#include <cstring>
#include <math.h>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <algorithm>
using namespace std;

double epsilonJ (double H, double Q, double betta)
{
    double epsilon;
    epsilon = (H / 2.0) + (Q * 1.0 / (betta * betta));
    return epsilon;
}

double betta0J (double delta, double H, double P, double betta)
{
    double betta0;
    betta0 = ((delta * H) + (2 * P)) / (2.0 * betta * betta);
    //cout << '\a';
    return betta0;
}

double hi0J (double H, double Q, double r, double betta)
{
    double hi0;
    hi0 = (H * H / 4.0) + (((H * Q) + (r * r)) / (betta * betta));
    return hi0;
}

double I3(double epsilon1, double deltaplus, double deltaminus, double zplus, double zminus)
{
    double I3;
    if (abs(deltaplus) > pow(10, -10))
    {
        if (((epsilon1 * epsilon1) - (deltaplus * deltaminus)) > pow(10, -10)) //+
        {
            I3 = (log(abs(((deltaplus * zplus) + epsilon1 - sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus))) / (((deltaplus * zplus) + epsilon1 + sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus)))))) / (2 * sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus)))) - (log(abs(((deltaplus * zminus) + epsilon1 - sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus))) / (((deltaplus * zminus) + epsilon1 + sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus)))))) / (2 * sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus))));
        }

        else if (((epsilon1 * epsilon1) - (deltaplus * deltaminus)) < pow(10, -10))
        {
            I3 = (atan(((deltaplus * zplus) + epsilon1) / (sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1)))) / sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1))) - (atan(((deltaplus * zminus) + epsilon1) / (sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1)))) / sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1)));
        }

        else if (abs((epsilon1 * epsilon1) - (deltaplus * deltaminus)) <= pow(10, -10))
        {
            I3 = ((-1.0) / ((deltaplus * zplus) + epsilon1)) + (1.0 / (deltaplus * zminus + epsilon1));
        }
    }
    else if (abs(deltaplus) <= pow(10, -10))
    {
        if (abs(epsilon1) > pow(10, -10))
        {
            I3 = (log(abs((2 * epsilon1 * zplus) + deltaminus)) / (2.0 * epsilon1)) - (log(abs((2 * epsilon1 * zminus) + deltaminus)) / (2 * epsilon1));
        }
        else
        {
            I3 = (zplus / deltaminus) - (zminus / deltaminus);
        }
    }
    return I3;
}

double I2(double deltaplus, double zplus, double zminus, double epsilon1, double deltaminus)
{
    double I;
    I = I3(epsilon1, deltaplus, deltaminus, zplus, zminus);
    double I2;
    //if (deltaplus != 0)
    //{
        I2 = ((zplus - (epsilon1 * log(abs((deltaplus * zplus * zplus) + (2 * epsilon1 * zplus) + deltaminus)) / deltaplus) - (zminus - (epsilon1 * log(abs((deltaplus * zminus * zminus) + (2 * epsilon1 * zminus) + deltaminus)) / deltaplus))) / deltaplus) + (2 * epsilon1 * epsilon1 - (deltaplus * deltaminus)) * I / (deltaplus * deltaplus);
    /*}
    else
    {
        I2 = 0;
    } */
    return I2;
}

double I4(double deltaminus, double zplus, double zminus, double epsilon1, double deltaplus)
{
    double I4;
    double I;
    I = I3(epsilon1, deltaplus, deltaminus, zplus, zminus);
    if (abs(deltaplus) > pow(10,-10))
    {
        //if (deltaminus != 0)
        //{
            I4 = ((((-1.0) / zplus) - (epsilon1 * log(zplus * zplus / abs((deltaplus * zplus * zplus) + (2 * epsilon1 * zplus) + deltaminus)) / deltaminus) - ((-1.0) / zminus) - (epsilon1 * log(zminus * zminus / abs((deltaplus * zminus * zminus) + (2 * epsilon1 * zminus) + deltaminus)) / deltaminus)) / deltaminus) + (((2 * epsilon1 * epsilon1) - (deltaplus * deltaminus)) * I / (deltaminus * deltaminus));
        /*}
        else
        {
            I4 = 0;
        } */
    }
    else
    {
        if (abs(epsilon1) > pow(10,-10))
        {
            I4 = (((-1) / (deltaminus * zplus)) + ((2 * epsilon1 * log(abs(((2 * epsilon1 * zplus) + deltaminus) / zplus))) / (deltaminus * deltaminus))) - (((-1) / (deltaminus * zminus)) + (2 * epsilon1 * log(abs((2 * epsilon1 * zminus + deltaminus) / zminus))) / (deltaminus * deltaminus));
        }

        else
        {
            I4 = ((-1.0) / (deltaminus * zplus)) + (1.0 / (deltaminus * zminus));
        }
    }

    return I4;
}

double I1(double deltaplus, double hi1, double deltaminus, double zplus, double zminus, double epsilon1)
{
    double I1;
    double II, III, IIII;
    II = I2(deltaplus, zplus, zminus, epsilon1, deltaminus);
    III = I3(epsilon1, deltaplus, deltaminus, zplus, zminus);
    IIII = I4(deltaminus, zplus, zminus, epsilon1, deltaplus);
    I1 = (deltaplus * II / 2.0) - (hi1 * III) + (deltaminus * IIII / 2.0);
    return I1;
}

double I5(double U1, double U2, double k1, double k2)
{
    double I5;
    if (abs(k1) > pow(10,-10))
    {
        I5 = (-U2 + (U2 + (k2 / k1)) * log(abs(k1 * U2 + k2))) - (-U1 + (U1 + (k2 / k1)) * log(abs(k1 * U1 + k2)));
    }
    else
    {
        I5 = (U2 - U1) * log(abs(k2));
    }
    return I5;
}

int main()
{
    int M = 25, N = 25;
    long double V0, S0 = 0;
    long double p = 0;
    const double pi = 3.14159265358979;
    double x1, x2, x3, x, y1, y2, y3, y11, y21, y31, h = 2.0 * pi / N, H = pi / (1.0 * M), etta, R;
    double A = N * h, B = M * H;
    double tetta, betta, y1prv, y2prv, y3prv, alfa, y1pru, y2pru, y3pru, delta;
    double P, Q, r1, r2, r3, r;
    double delta0, alfa0, epsilon, betta0, hi0, hi1;
    double deltaplus, deltaminus, delta1, epsilon1;
    double IH, zminus, zplus, IminusH;
    setlocale(LC_ALL, "Rus");
    for (int i = 1; i < 6; i++)
    {
        R = 1 - pow(10, (-1)*i);
        cout << "R=" << R << endl;
        for (int q = 0; q <= 2; q++)
        {
            double uq = pi * q / (1.0 * N);

            for (int l = 0; l <= 2 * M; l++)
            {
                double vl = pi * l / (2.0 * M);
                y11 = sin(vl)*cos(uq); // !
                y21 = sin(vl)*sin(uq); // !
                y31 = cos(vl); // !
                x1 = y11 * R;
                x2 = y21 * R;
                x3 = y31 * R;
                x = sqrt((x1 * x1) + (x2 * x2) + (x3 * x3));

                for (int m = 0; m < M; m++)
                {
                    double vm = (m + 0.5)*H;
                    for (int n = 0; n < N; n++)
                    {
                        double un = (n + 0.5)*h;
                        y1 = sin(vm)*cos(un);
                        y2 = sin(vm)*sin(un);
                        y3 = cos(vm);
                        y1prv = cos(un) * cos(vm);
                        y2prv = sin(un) * cos(vm);
                        y3prv = (-1)*sin(vm);
                        y1pru = (-1) * sin(un) * sin(vm);
                        y2pru = sin(vm) * cos(un);
                        y3pru = 0;
                        r1 = y1 - x1;
                        r2 = y2 - x2;
                        r3 = y3 - x3;
                        r = sqrt(r1 * r1 + r2 * r2 + r3 * r3);
                        //cout << "r=" << r << endl;
                        betta = sqrt((y1prv * y1prv) + (y2prv * y2prv) + (y3prv * y3prv));
                        etta = sin(vm);
                        alfa = sqrt((y1pru * y1pru) + (y2pru * y2pru) + (y3pru * y3pru));
                        delta = (y1pru * y1prv) + (y2pru * y2prv) + (y3pru * y3prv);
                        P = (r1 * y1pru) + (r2 * y2pru) * (r3 * y3pru);
                        Q = (r1 * y1prv) + (r2 * y2prv) + (r3 * y3prv);
                        delta0 = delta * 1.0 / (betta * betta);
                        alfa0 = abs(alfa / betta); //!
                        epsilon = epsilonJ(H, Q, betta);
                        betta0 = betta0J(delta, H, P, betta);
                        hi0 = hi0J(H, Q, r, betta);
                        hi1 = sqrt(hi0 - (betta0 * betta0 / (alfa0 * alfa0)));
                        epsilon1 = epsilon - (delta0 * betta0 / (alfa0 * alfa0));
                        delta1 = delta0 * hi1 / alfa0;
                        deltaplus = hi1 + delta1;
                        deltaminus = hi1 - delta1;
                        if (abs(hi1) > pow(10,-10))
                        {
                            //IH = 0;
                            double dzettaplus = ((alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
                            double dzettaminus = (((-1) * alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
                            zplus = dzettaplus + sqrt((dzettaplus * dzettaplus) + 1);
                            zminus = dzettaminus + sqrt((dzettaminus * dzettaminus) + 1);
                            IH = hi1 * (((dzettaplus * log(abs(epsilon1 + (delta1 * dzettaplus) + hi1 * sqrt((dzettaplus * dzettaplus) + 1))) - (dzettaminus * log(abs(epsilon1 + (delta1 * dzettaminus) + hi1 * sqrt((dzettaminus * dzettaminus) + 1))))))  - I1(deltaplus, hi1, deltaminus, zplus, zminus, epsilon1) )/ alfa0;
                        }
                        else if (abs(hi1) <= pow(10,-10))
                        {
                            //IH = 0;
                            if (-h / 2 < (-betta0 / (alfa0 * alfa0)) && (-betta0 / (alfa0 * alfa0)) < h / 2)
                            {
                                IH = I5(-h / 2, (-betta0 / (alfa0 * alfa0)), delta0 - alfa0, epsilon - (betta0 / alfa0)) + I5((-betta0 / (alfa0 * alfa0)), h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
                            }
                            else if ((-betta0 / (alfa0 * alfa0)) >= h / 2)
                            {
                                IH = I5(-h / 2, h / 2, delta0 - alfa0, epsilon - (betta0 / alfa0));
                            }
                            else if ((-betta0 / (alfa0 * alfa0)) <= h / 2)
                            {
                                IH = I5(-h / 2, h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
                            }
                        }
                        epsilon = epsilonJ((-1) * H, Q, betta);
                        betta0 = betta0J(delta, (-1) * H, P, betta);
                        hi0 = hi0J((-1) * H, Q, r, betta);
                        hi1 = sqrt(hi0 - (betta0 * betta0 / (alfa0 * alfa0)));
                        epsilon1 = epsilon - (delta0 * betta0 / (alfa0 * alfa0));
                        delta1 = delta0 * hi1 / alfa0;
                        deltaplus = hi1 + delta1;
                        deltaminus = hi1 - delta1;
                        if (abs(hi1) > pow(10, -10))
                        {
                            //IminusH = 1;
                            double dzettaplus = ((alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
                            double dzettaminus = (((-1) * alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
                            zplus = dzettaplus + sqrt(dzettaplus * dzettaplus + 1);
                            zminus = dzettaminus + sqrt(dzettaminus * dzettaminus + 1);
                            IminusH = hi1 * ((dzettaplus * (log(abs(epsilon1 + (delta1 * dzettaplus) + hi1 * sqrt((dzettaplus * dzettaplus) + 1)))) - (dzettaminus * log(abs(epsilon1 + (delta1 * dzettaminus) + hi1 * sqrt((dzettaminus * dzettaminus) + 1))))) - I1(deltaplus, hi1, deltaminus, zplus, zminus, epsilon1)) / alfa0;
                        }
                        else if (abs(hi1) <= pow(10, -10))
                        {
                            //IminusH = 0;
                            if (-h / 2 < (-betta0 / (alfa0 * alfa0)) && (-betta0 / (alfa0 * alfa0)) < h / 2)
                            {
                                IminusH = I5(-h / 2, (-betta0 / (alfa0 * alfa0)), delta0 - alfa0, epsilon - (betta0 / alfa0)) + I5((-betta0 / (alfa0 * alfa0)), h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
                            }
                            else if ((-betta0 / (alfa0 * alfa0)) >= h / 2)
                            {
                                IminusH = I5(-h / 2, h / 2, delta0 - alfa0, epsilon - (betta0 / alfa0));
                            }
                            else if ((-betta0 / (alfa0 * alfa0)) <= h / 2)
                            {
                                IminusH = I5(-h / 2, h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
                            }
                            //cout << '\a';
                        }
                        //xminusy = pow(pow(x1 - y1, 2) + pow(x2 - y2, 2) + pow(x3 - y3, 2), 0.5); //r
                        tetta = 1.0 * (IH - IminusH)/ betta;
                        S0 = S0 + (etta * tetta);
                    }
                }
                if (abs(x) < 1)
                {
                    V0 = 4 * pi;
                }
                else if (abs(x) > 1)
                {
                    V0 = 4.0 * pi / abs(x);
                }
                //cout << S0 << endl;
                //cout << V0 << endl;
                long double maximum = max(p, abs(S0 - V0) * 1.0 / V0);
                S0 = 0;
                p = maximum;
            }
        }
        cout << "Относительная погрешность=" << p << endl;
        cout << "______________________" << endl;
    }

    return 0;
}