Hi John,
I am printing field values at random points to make sure the interpolation is fine. It does print out what the values I expect, 1000 V/m for Ex and 0 for the rest. I am doing a simple nearest-neighbor interpolation in this case. Here’s my definition of MyElectricField:
#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; } // parses the line into variables, if fail, loop breaks
fieldData.push_back({x, y, z, Ex, Ey, Ez}); // adds parsed values to the fieldData vector
}
infile.close();
}
void MyElectricField::InterpolateField(const G4double point[4], G4double* field) const
{
// simple nearest-neighbor interpolation
G4double minDist = std::numeric_limits<G4double>::max(); //minDist is set to a large value
const FieldPoint* nearestPoint = nullptr; // initialize pointer
for (const auto& fp : fieldData) // loop over all field data points
{
// calculate the distance between the query point and the current field data point
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); // compute the euclidean distance to the input point
if (dist < minDist)
{
minDist = dist;
nearestPoint = &fp;
}
}
if (nearestPoint) //if nearest point is found, the field array is set to values of the nearestPoint
{
field[0] = nearestPoint->Ex;
field[1] = nearestPoint->Ey;
field[2] = nearestPoint->Ez;
}
else
{
G4cerr << "Warning: no nearest point found." << G4endl; // if no nearestPoint found, values remains unchanged
}
}
// debugging 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;
}
}
and for completeness my corresponding .hh file:
#ifndef MYELECTRICFIELD_HH
#define MYELECTRICFIELD_HH
#include "G4ElectricField.hh"
#include "globals.hh"
#include <vector>
class MyElectricField : public G4ElectricField
{
public:
MyElectricField(const G4String& filename);
virtual ~MyElectricField();
virtual void GetFieldValue(const G4double point[4], G4double* field) const override;
// method to print field values at specified points
void PrintFieldValues(const std::vector<G4ThreeVector>& points) const;
private:
struct FieldPoint {
G4double x, y, z;
G4double Ex, Ey, Ez;
};
std::vector<FieldPoint> fieldData;
void LoadFieldData(const G4String& filename);
void InterpolateField(const G4double point[4], G4double* field) const;
};
#endif