Обратная задача 2

#include <iostream>
#include <cassert>
#include <fstream>
#include <vector>
#include <cmath>
#include <random>


using namespace std;

#define M_PI 3.1415926
#define EPS  0.0000001

struct Point3d {
    double x;
    double y;
    double z;
};

struct Range {
    Point3d start;
    Point3d end;
};

Range source; // Источник
vector<Range>  receivers;  // Приёмники
vector<double> potentials; // Разность потенциалов измеренное на каждом приёмнике

double calculatedSigma = 0.01;
vector<double> calculatedPotential;

constexpr double goldenSigma = 0.1;
constexpr double goldenI = 1;




double distance(Point3d first, Point3d second) {
    return sqrt(pow(first.x - second.x, 2) +
                pow(first.y - second.y, 2) +
                pow(first.z - second.z, 2));
}

void calculateNextPotentials() {
    // Генерация синтетический данных
    for (size_t i = 0; i < receivers.size(); i++) {
        double U = goldenI / (2 * M_PI * calculatedSigma) * ((1 / distance(source.end, receivers[i].start) - 1 / distance(source.start, receivers[i].start)) -
            (1 / distance(source.end, receivers[i].end) - 1 / distance(source.start, receivers[i].end)));
        calculatedPotential[i] = U;
        //std::cout << "calculatedPotential[" << i << "] = " << U << std::endl;
    }
}

// Минимизируемый функционал
//TODOOOOOOOOOO
double getF() {
    double F = 0;

    for (size_t i = 0; i < receivers.size(); i++) {
        F += 1 / potentials[i] * pow(potentials[i] - calculatedPotential[i], 2);
    }

    return F;
}

int main() {
    std::random_device dev;
    std::mt19937 rng(dev());
    std::uniform_int_distribution<std::mt19937::result_type> dist(90, 110);

    ifstream input("example1.txt");
    assert(input.is_open() && "Error while opening file");

    size_t sourcesCount;
    size_t receiversCount;

    input >> sourcesCount;
    input >> receiversCount;

    calculatedPotential.resize(receiversCount);

    // Считываем источник
    input >> source.start.x >> source.start.y >> source.start.z;
    input >> source.end.x >> source.end.y >> source.end.z;

    // Считываем приёмники
    for (size_t i = 0; i < receiversCount; i++) {
        Range receiver;

        input >> receiver.start.x >> receiver.start.y >> receiver.start.z;
        input >> receiver.end.x >> receiver.end.y >> receiver.end.z;

        receivers.emplace_back(receiver);
    }


    // Считываем разность потенциалов на каждом приёмнике
    for (size_t i = 0; i < receiversCount; i++) {
        double potential;
        input >> potential;

        double noise = (dist(rng) / 100.0);
        potential *= noise;

        potentials.emplace_back(potential);
    }

    vector<double> diff(receivers.size(), 0);
    calculateNextPotentials();

    std::cout << calculatedSigma << " " << getF() << endl;
    // ---------------------- //
    double prevSigma = calculatedSigma;
    do {
        prevSigma = calculatedSigma;

        double elemA = 0; // Элемент СЛАУ
        double elemB = 0; // Правая часть уравнения
        for (size_t i = 0; i < receivers.size(); i++) {
            diff[i] = -goldenI / (2 * M_PI * pow(calculatedSigma,2)) * ((1 / distance(source.end, receivers[i].start) - 1 / distance(source.start, receivers[i].start)) -
                (1 / distance(source.end, receivers[i].end) - 1 / distance(source.start, receivers[i].end)));

            elemA += pow((1 / potentials[i]) * diff[i], 2);
            elemB -= pow(1 / potentials[i], 2) * diff[i] * (calculatedPotential[i] - potentials[i]);
        }
        double deltaSigma = elemB / elemA;
        calculatedSigma += deltaSigma;

        // Пересчитываем potenrial
        calculateNextPotentials();

        std::cout << calculatedSigma << " " << getF() << std::endl;
    } while (abs(prevSigma - calculatedSigma) > EPS);

    system("pause");
    return 0;
}