В процессе

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

int main()
{
    int M=100, N=100;
    long double V0, K0 = 0, S0 = 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);
    for (int i = 1; i < 6; i++)
    {
        R = 1 - pow(10, (-1)*i);
        cout << "R=" << R << endl;
        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;
            // Вычисление K0
        /*  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);
                    K0 = K0 + (h * H * etta * 1.0 / xminusy);
                }
            } */
            // Вычисление S0
            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;
                    y1pr = cos(un) * cos(vm);
                    y2pr = sin(un) * cos(vm);
                    y3pr = (-1)*sin(vm);
                    betta = sqrt(pow(y1pr, 2) + pow(y2pr, 2) + pow(y3pr, 2));

                    // Вычисление I(H)






                    tetta = 1.0 / betta;
                    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);
                    S0 = S0 + (etta * tetta);
                }
            }

            cout << "K0=" << K0 << endl;
            if (abs(x) < 1)
            {
                V0 = 4 * pi;
            }
            else if (abs(x) > 1)
            {
                V0 = 4.0 * pi / abs(x);
            }

            cout << "V0=" << V0 << endl;
            cout << "______________________" << endl;
    }

    return 0;
}