#include "MyElectricField.hh"
#include <fstream>
#include <sstream>
#include <cmath>
#include <iostream>

MyElectricField::MyElectricField(const G4String& filename)
{
    LoadFieldData(filename);
}

MyElectricField::~MyElectricField() {}

void MyElectricField::GetFieldValue(const G4double point[4], G4double* field) const
{
    InterpolateField(point, field);
}

void MyElectricField::LoadFieldData(const G4String& filename)
{
    std::ifstream infile(filename);
    if (!infile.is_open())
    {
        G4cerr << "Error opening file: " << filename << G4endl;
        return;
    }
    G4cout << "Successfully opened field data file: " << filename << G4endl;

    G4double x, y, z, Ex, Ey, Ez;
    std::string line;
    while (std::getline(infile, line))
    {
        if (line[0] == '%') continue; // Skip header lines
        std::istringstream iss(line);
        if (!(iss >> x >> y >> z >> Ex >> Ey >> Ez)) { break; }
        fieldData.push_back({x, y, z, Ex, Ey, Ez});
    }
    infile.close();
}

void MyElectricField::InterpolateField(const G4double point[4], G4double* field) const
{
    // nearest-neighbor interpolation for demonstration purposes
    G4double minDist = std::numeric_limits<G4double>::max();
    const FieldPoint* nearestPoint = nullptr;

    for (const auto& fp : fieldData)
    {
        G4double dx = point[0] - fp.x;
        G4double dy = point[1] - fp.y;
        G4double dz = point[2] - fp.z;
        G4double dist = std::sqrt(dx*dx + dy*dy + dz*dz);

        if (dist < minDist)
        {
            minDist = dist;
            nearestPoint = &fp;
        }
    }

    if (nearestPoint)
    {
        field[0] = nearestPoint->Ex;
        field[1] = nearestPoint->Ey;
        field[2] = nearestPoint->Ez;
    }
    else
    {
        field[0] = field[1] = field[2] = 0.0;
    }
}

// function to print field values at specified points
void MyElectricField::PrintFieldValues(const std::vector<G4ThreeVector>& points) const 
{
    for (const auto& point : points) 
    {
        G4double g4point[4] = { point.x(), point.y(), point.z(), 0.0 };
        G4double field[3] = { 0.0, 0.0, 0.0 };
        GetFieldValue(g4point, field);
        G4cout << "Field at (" << point.x()  << " mm, " << point.y()<< " mm, " << point.z() << " mm) : ("
               << field[0]  << " V/m, " << field[1] << " V/m, " 
               << field[2]  << " V/m)" << G4endl;
    }
}