Тест 3 "ур-е Гельмгольца" (внеш. сфера) (старая формула)

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

int main()
{
    int M = 50, N = 50;
    complex <double> V1;
    complex <double> K1;
    complex <double> p(0, 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, xminusy, R;
    double A = N * h, B = M * H;
    //double tetta, betta, y1pr, y2pr, y3pr;
    // double uq = pi / (1.0 * N), vl = pi / (2.0 * M);
    complex <double> z(0, 1.0); // мнимая единица
    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 * N; 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;
                //cout << "x1=" << x1 << endl;
                //cout << "x2=" << x2 << endl;
                //cout << "x3=" << x3 << endl;
                x = pow(pow(x1, 2) + pow(x2, 2) + pow(x3, 2), 0.5);
                //cout << "x=" << x << endl;
                // Вычисление K1
                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);
                        etta = sin(vm);
                        xminusy = pow(pow(x1 - y1, 2) + pow(x2 - y2, 2) + pow(x3 - y3, 2), 0.5);
                        K1 = K1 + ((h * H * etta * 1.0 * exp(z * xminusy) / xminusy) / (4 * pi));
                    }
                }
                if (abs(x) < 1)
                {
                    V1 = exp(z) * sin(abs(x)) / abs(x);
                }
                else if (abs(x) > 1)
                {
                    V1 = exp(z * abs(x)) * sin(1) / abs(x);
                }
                complex <double> maximum = max(abs(p), abs(K1 - V1) * 1.0 / abs(V1));
                K1 = 0;
                p = maximum;
            }
        }

        cout << "Относительная погрешность=" << abs(p) << endl;
        //cout << "V1=" << abs(V1) << endl;
        cout << "______________________" << endl;
    }

    return 0;
}