EMTElement.cpp

#include "EMTElement.h"
 
#include <ranges>
#include <wx/app.h>
 
#include "../../forms/EMTElementForm.h"
#include "../../editors/Workspace.h"
#include "../../utils/PropertiesData.h"
#include "../../extLibs/fftw/fftw3.h"
#include "../../elements/GCText.h"
#include "../../utils/Path.h"
 
#include <wx/dcgraph.h>
#include <wx/textfile.h>
#include <wx/dir.h>
#include <wx/process.h>
#include <wx/wfstream.h>
#include <wx/datstrm.h>
#include <wx/busyinfo.h>
 
EMTElement::EMTElement()
{
}
 
EMTElement::EMTElement(wxString name)
{
    m_data.name = name;
}
 
EMTElement::~EMTElement()
{
}
 
Element* EMTElement::GetCopy()
{
    EMTElement* copy = new EMTElement();
    *copy = *this;
    return copy;
}
 
bool EMTElement::AddParent(Element* parent, wxPoint2DDouble position)
{
    if (parent) {
        m_parentList.push_back(parent);
        parent->AddChild(this);
        wxPoint2DDouble parentPt =
            parent->RotateAtPosition(position, -parent->GetAngle());        // Rotate click to horizontal position.
        parentPt.m_y = parent->GetPosition().m_y;                           // Centralize on bus.
        parentPt = parent->RotateAtPosition(parentPt, parent->GetAngle());  // Rotate back.
 
        m_position = parentPt + wxPoint2DDouble(0.0, 100.0);  // Shifts the position to the down of the bus.
        m_width = 50;
        m_height = 50;
        m_rect = wxRect2DDouble(m_position.m_x - m_width / 2.0, m_position.m_y - m_height / 2.0, m_width, m_height);
 
        m_pointList.push_back(parentPt);
        m_pointList.push_back(GetSwitchPoint(parent, parentPt, m_position));
        m_pointList.push_back(m_position + wxPoint2DDouble(0.0, -m_height / 2.0 - 10.0));
        m_pointList.push_back(m_position + wxPoint2DDouble(0.0, -m_height / 2.0));
 
        m_inserted = true;
 
        wxRect2DDouble genRect(0, 0, 0, 0);
        m_switchRect.push_back(genRect);  // Push a general rectangle.
        UpdateSwitches();
 
        // Base data
        auto data = static_cast<Bus*>(parent)->GetElectricalData();
        m_data.puVoltage = data.voltage;
        m_data.baseVoltage = GetValueFromUnit(data.nominalVoltage, data.nominalVoltageUnit);
        m_data.frequency = data.stabFreq;
 
        return true;
    }
    return false;
}
 
void EMTElement::DrawDC(GUIColour* guiColour, wxPoint2DDouble translation, double scale, wxGraphicsContext* gc) const
{
    wxColour elementColour;
    if (m_online) {
        if (m_dynEvent)
            elementColour = guiColour->eventElement;
        else
            elementColour = guiColour->enabled;
    }
    else
        elementColour = guiColour->disable;
 
    if (m_inserted) {
 
        if (m_selected) {
            gc->SetPen(wxPen(guiColour->selection, 2.0 + m_borderSize * 2.0));
            gc->SetBrush(*wxTRANSPARENT_BRUSH);
 
            gc->StrokeLines(m_pointList.size(), &m_pointList[0]);
 
            // Push the current matrix on stack.
            gc->PushState();
            // Rotate the matrix around the object position.
            gc->Translate(m_position.m_x, m_position.m_y);
            gc->Rotate(wxDegToRad(m_angle));
            gc->Translate(-m_position.m_x, -m_position.m_y);
 
            //gc->DrawRectangle(m_position.m_x - m_width / 2.0, m_position.m_y - m_height / 2.0, m_width, m_height);
            gc->DrawRoundedRectangle(m_position.m_x - m_width / 2.0, m_position.m_y - m_height / 2.0, m_width, m_height, 10.0);
 
            gc->PopState();
 
            // Draw node selection.
            gc->SetPen(*wxTRANSPARENT_PEN);
            gc->SetBrush(wxBrush(guiColour->selection));
            DrawDCCircle(m_pointList[0], 5.0 + m_borderSize / scale, 10, gc);
        }
        // Draw EMTElement (layer 2).
        // Draw node.
        gc->SetPen(*wxTRANSPARENT_PEN);
        gc->SetBrush(wxBrush(elementColour));
        DrawDCCircle(m_pointList[0], 5.0, 10, gc);
 
        gc->SetPen(wxPen(elementColour, 2));
        gc->SetBrush(*wxTRANSPARENT_BRUSH);
        gc->StrokeLines(m_pointList.size(), &m_pointList[0]);
 
        DrawDCSwitches(guiColour, gc);
 
        // Push the current matrix on stack.
        gc->PushState();
        // Rotate the matrix around the object position.
        gc->Translate(m_position.m_x, m_position.m_y);
        gc->Rotate(wxDegToRad(m_angle));
        gc->Translate(-m_position.m_x, -m_position.m_y);
 
        gc->SetPen(wxPen(elementColour, 2));
        gc->SetBrush(wxBrush(guiColour->background));
        //gc->DrawRectangle(m_position.m_x - m_width / 2.0, m_position.m_y - m_height / 2.0, m_width, m_height);
        gc->DrawRoundedRectangle(m_position.m_x - m_width / 2.0, m_position.m_y - m_height / 2.0, m_width, m_height, 10.0);
 
        //gc->SetPen(*wxBLACK_PEN);
        //wxFont font(10, wxFONTFAMILY_ROMAN, wxFONTSTYLE_NORMAL, wxFONTWEIGHT_BOLD);
        wxFont font;
        font.SetFaceName(wxT("CMU Serif"));
        font.SetPointSize(10);
        font.MakeBold();
        gc->SetFont(font, m_online ? wxColour(255, 60, 0) : m_offlineElementColour);
        double textWidth, textHeight;
        gc->GetTextExtent(_("EMT"), &textWidth, &textHeight);
        gc->DrawText(_("EMT"), m_position.m_x - textWidth / 2.0, m_position.m_y - textHeight / 2.0 + 20.0);
 
        wxGCDC* gcdc = new wxGCDC(gc);
        gcdc->SetPen(wxPen(m_online ? wxColour(0, 60, 255) : m_offlineElementColour, 2));
        std::vector<wxPoint> ptList;
        for (double x = m_position.m_x - m_width / 2.0 + 2; x < (m_position.m_x + m_width / 2.0); x += (m_width - 4.0) / 6.0) {
            ptList.emplace_back(x, m_position.m_y + std::sin((x - (m_position.m_x - m_width / 2.0 + 2)) / m_width * 6 * M_PI) * 20.0);
        }
        gcdc->DrawSpline(ptList.size(), ptList.data());
 
        gc->PopState();
    }
}
 
void EMTElement::DrawDC(GUIColour* guiColour, wxPoint2DDouble translation, double scale, wxDC& dc) const
{
    wxColour elementColour;
    if (m_online) {
        if (m_dynEvent)
            elementColour = guiColour->eventElement;
        else
            elementColour = guiColour->enabled;
    }
    else
        elementColour = guiColour->disable;
 
    std::vector<wxPoint> pointListInt;
    for (auto& pt : m_pointList) {
        pointListInt.emplace_back(static_cast<int>(pt.m_x), static_cast<int>(pt.m_y));
    }
 
    if (m_inserted) {
 
        if (m_selected) {
            dc.SetPen(wxPen(guiColour->selection, 2 + m_borderSize * 2.0));
            dc.SetBrush(*wxTRANSPARENT_BRUSH);
 
            dc.DrawLines(pointListInt.size(), &pointListInt[0]);
 
            DrawDCRoundedRectRotated(dc, m_position, m_width, m_height, 10.0, m_angle);
 
            // Draw node selection.
            dc.SetPen(*wxTRANSPARENT_PEN);
            dc.SetBrush(wxBrush(guiColour->selection));
            DrawDCCircle(m_pointList[0], 5.0 + m_borderSize / scale, dc);
        }
        // Draw EMTElement (layer 2).
        // Draw node.
        dc.SetPen(*wxTRANSPARENT_PEN);
        dc.SetBrush(wxBrush(elementColour));
        DrawDCCircle(pointListInt[0], 5.0, dc);
 
        dc.SetPen(wxPen(elementColour, 2));
        dc.SetBrush(*wxTRANSPARENT_BRUSH);
        dc.DrawLines(pointListInt.size(), &pointListInt[0]);
 
        DrawDCSwitches(guiColour, dc);
 
        dc.SetPen(wxPen(elementColour, 2));
        dc.SetBrush(wxBrush(guiColour->background));
        DrawDCRoundedRectRotated(dc, m_position, m_width, m_height, 10.0, m_angle);
 
        wxFont font;
        font.SetFaceName(wxT("CMU Serif"));
        font.SetPointSize(10);
        font.MakeBold();
        GCText gcText;
        gcText.SetFont(font);
        gcText.SetText(_("EMT"));
        wxPoint pt = RotateAround(m_position + wxPoint2DDouble(0.0, 20.0), m_position, m_angle);
        gcText.Draw(pt, gcText.GetWidth(), gcText.GetHeight(), dc, m_angle, m_online ? wxColour(255, 60, 0) : m_offlineElementColour);
 
        dc.SetPen(wxPen(m_online ? wxColour(0, 60, 255) : m_offlineElementColour, 2));
        std::vector<wxPoint> ptList;
        for (double x = m_position.m_x - m_width / 2.0 + 2; x < (m_position.m_x + m_width / 2.0); x += (m_width - 4.0) / 6.0) {
            wxPoint pt = RotateAround(wxPoint2DDouble(x, m_position.m_y + std::sin((x - (m_position.m_x - m_width / 2.0 + 2)) / m_width * 6 * M_PI) * 20.0), m_position, m_angle);
            ptList.emplace_back(pt);
        }
        dc.DrawSpline(ptList.size(), ptList.data());
    }
}
 
void EMTElement::Rotate(bool clockwise)
{
    double rotAngle = m_rotationAngle;
    if (!clockwise) rotAngle = -m_rotationAngle;
 
    m_angle += rotAngle;
    if (m_angle >= 360 || m_angle <= -360) m_angle = 0.0;
    m_pointList[2] = RotateAtPosition(m_pointList[2], rotAngle);
    m_pointList[3] = RotateAtPosition(m_pointList[3], rotAngle);
    UpdateSwitchesPosition();
}
 
bool EMTElement::GetContextMenu(wxMenu& menu)
{
    menu.Append(ID_EDIT_ELEMENT, _("Edit Electromagnetic Transient"));
    GeneralMenuItens(menu);
    return true;
}
 
wxString EMTElement::GetTipText() const
{
    wxString tipText = m_data.name;
    tipText += wxT("\n");
 
    tipText += wxString::Format(_("\nP = %.5f p.u."), m_data.power.real());
    tipText += wxString::Format(_("\nQ = %.5f p.u."), m_data.power.imag());
    if (auto itCurrrent = m_data.currHarmonics.find(1); itCurrrent != m_data.currHarmonics.end())
        tipText += wxString::Format(_("\nI = %.5f A"), std::abs(itCurrrent->second));
 
    wxString harmonicsInfo = _("\n\nHarmonics info:");
    bool hasHarmonics = false;
    for (auto& [order, current] : m_data.currHarmonics) {
        if (order != 1) {
            hasHarmonics = true;
            harmonicsInfo += wxString::Format(_("\nIh(%d): %.5f%s%.2f%s A"), order, std::abs(current), wxString(L'\u2220'), wxRadToDeg(std::arg(current)), wxString(L'\u00B0'));
        }
    }
 
    if (hasHarmonics) tipText += harmonicsInfo;
 
    return tipText;
}
 
bool EMTElement::ShowForm(wxWindow* parent, Element* element, wxWindow* workspace)
{
    Workspace* ws = dynamic_cast<Workspace*>(workspace);
    if (!ws) return false;
 
    EMTElementForm emtForm(parent, this, ws);
    emtForm.SetTitle(_("Electromagnetic Transient"));
    emtForm.CenterOnParent();
    if (emtForm.ShowModal() == wxID_OK) {
        return true;
    }
    return false;
}
 
rapidxml::xml_node<>* EMTElement::SaveElement(rapidxml::xml_document<>& doc, rapidxml::xml_node<>* elementListNode)
{
    auto elementNode = XMLParser::AppendNode(doc, elementListNode, "EMTElement");
    XMLParser::SetNodeAttribute(doc, elementNode, "ID", m_elementID);
 
    SaveCADProperties(doc, elementNode);
 
    auto electricalProp = XMLParser::AppendNode(doc, elementNode, "ElectricalProperties");
    auto isOnline = XMLParser::AppendNode(doc, electricalProp, "IsOnline");
    XMLParser::SetNodeValue(doc, isOnline, m_online);
    auto name = XMLParser::AppendNode(doc, electricalProp, "Name");
    XMLParser::SetNodeValue(doc, name, m_data.name);
    auto atpFilePath = XMLParser::AppendNode(doc, electricalProp, "ATPFilePath");
    XMLParser::SetNodeValue(doc, atpFilePath, m_data.atpFile.GetFullPath());
    auto atpNodeName = XMLParser::AppendNode(doc, electricalProp, "ATPNodeName");
    XMLParser::SetNodeValue(doc, atpNodeName, m_data.atpNodeName);
    auto stepSize = XMLParser::AppendNode(doc, electricalProp, "StepSize");
    XMLParser::SetNodeValue(doc, stepSize, m_data.stepSize);
    auto cyclesToSS = XMLParser::AppendNode(doc, electricalProp, "CyclesToSS");
    XMLParser::SetNodeValue(doc, cyclesToSS, m_data.cyclesToSS);
    auto recordFrequency = XMLParser::AppendNode(doc, electricalProp, "RecordFrequency");
    XMLParser::SetNodeValue(doc, recordFrequency, m_data.recordFrequency);
    auto useMedianFilter = XMLParser::AppendNode(doc, electricalProp, "UseMedianFilter");
    XMLParser::SetNodeValue(doc, useMedianFilter, m_data.useMedianFilter);
    auto numMaxHarmonics = XMLParser::AppendNode(doc, electricalProp, "NumMaxHarmonics");
    XMLParser::SetNodeValue(doc, numMaxHarmonics, m_data.numMaxHarmonics);
 
    return elementNode;
}
 
bool EMTElement::OpenElement(rapidxml::xml_node<>* elementNode, std::vector<Element*> parentList)
{
    if (!OpenCADProperties(elementNode, parentList)) return false;
 
    auto electricalProp = elementNode->first_node("ElectricalProperties");
    if (!electricalProp) return false;
 
    SetOnline(XMLParser::GetNodeValueInt(electricalProp, "IsOnline"));
    m_data.name = electricalProp->first_node("Name")->value();
    m_data.atpFile = wxFileName(electricalProp->first_node("ATPFilePath")->value());
    m_data.atpNodeName = electricalProp->first_node("ATPNodeName")->value();
    m_data.stepSize = XMLParser::GetNodeValueDouble(electricalProp, "StepSize");
    m_data.cyclesToSS = XMLParser::GetNodeValueInt(electricalProp, "CyclesToSS");
    m_data.recordFrequency = XMLParser::GetNodeValueInt(electricalProp, "RecordFrequency");
    m_data.useMedianFilter = XMLParser::GetNodeValueInt(electricalProp, "UseMedianFilter");
    m_data.numMaxHarmonics = XMLParser::GetNodeValueInt(electricalProp, "NumMaxHarmonics");
 
    m_inserted = true;
    return true;
}
 
wxArrayString EMTElement::GetATPNodes(wxArrayString atpFile)
{
    wxArrayString nodeList;
    uint32_t mode = 0;
    std::vector<wxString> mode0Cards{ "/OUTPUT", "BLANK", "TACS", "INCLUDE", "TRANSFORMER" };
    std::vector<wxString> mode1Cards{ "/SOURCE" };
    std::vector<wxString> mode2Cards{ "/BRANCH", "/SWITCH", };
    wxString ground = "000000";
    for (wxString line : atpFile) {
        if (line.IsEmpty()) continue;
        if (tolower(line[0]) == 'c') continue; // Skip comments
 
        bool cardHeader = false;
        for (wxString card : mode0Cards) {
            if (line.Find(card) != wxNOT_FOUND) {
                mode = 0;
                cardHeader = true;
                break;
            }
        }
        for (wxString card : mode1Cards) {
            if (line.Find(card) != wxNOT_FOUND) {
                mode = 1;
                cardHeader = true;
                break;
            }
        }
        for (wxString card : mode2Cards) {
            if (line.Find(card) != wxNOT_FOUND) {
                mode = 2;
                cardHeader = true;
                break;
            }
        }
        if (cardHeader) continue;
 
 
        wxString node = "";
        switch (mode)
        {
        case 1: {
            node = line(2, 6).Trim();
            node.Replace(" ", "0");
            if (node != ground && node[5] >= 'A') { // Ignore ground and monophasic nodes
                if (nodeList.Index(node(0, 5)) == wxNOT_FOUND)
                    nodeList.Add(node(0, 5));
            }
        } break;
        case 2: {
            node = line(2, 6);
            if (node != ground && node[5] >= 'A') { // Ignore ground and monophasic nodes
                if (nodeList.Index(node(0, 5)) == wxNOT_FOUND)
                    nodeList.Add(node(0, 5));
            }
            node = line(8, 6);
            if (node != ground && node[5] >= 'A') { // Ignore ground and monophasic nodes
                if (nodeList.Index(node(0, 5)) == wxNOT_FOUND)
                    nodeList.Add(node(0, 5));
            }
        } break;
        }
    }
    return nodeList;
}
 
bool EMTElement::SetATPParameter(wxTextFile& atpFile, const wxString& card, const int& line, const int& initPos, const int& size, const wxString& value)
{
    bool foundCard = false;
    int currLine = -1;
    wxString lineStr = "";
    lineStr = atpFile.GetFirstLine();
    while (!atpFile.Eof())
    {
        if (tolower(lineStr[0]) == 'c') {
            lineStr = atpFile.GetNextLine();
            continue; // Skip comments
        }
        if (lineStr.Find(card) != wxNOT_FOUND) {
            foundCard = true;
            currLine = -1;
        }
        if (foundCard) currLine++;
        if (foundCard && currLine == line) {
            for (int i = initPos; i < initPos + size; i++) {
                lineStr[i] = value[i - initPos];
            }
            atpFile.RemoveLine(atpFile.GetCurrentLine());
            atpFile.InsertLine(lineStr, atpFile.GetCurrentLine());
            return true;
        }
 
        lineStr = atpFile.GetNextLine();
    }
    return false;
}
 
bool EMTElement::AddConnectionToNode(wxTextFile& atpFile, const wxString& node)
{
    // Get connnected bus data
    bool hasBusData = false;
    BusElectricalData busData;
    if (!m_parentList.empty()) {
        if (m_parentList[0] != nullptr) {
            Bus* bus = static_cast<Bus*>(m_parentList[0]);
            busData = bus->GetElectricalData();
            hasBusData = true;
        }
    }
 
    std::map<wxString, wxString> masks;
 
    // Load masks from T94.atp file
    wxTextFile t94File(Paths::GetDataPath() + "/atp/T94.atp");
    if (t94File.Open()) {
        wxString current;
 
        for (size_t i = 0; i < t94File.GetLineCount(); i++) {
            wxString line = t94File.GetLine(i);
 
            if (line.StartsWith("[[") && line.EndsWith("]]")) {
                current = line.Mid(2, line.Length() - 4);
                masks[current].Clear();
                continue;
            }
 
            if (!current.IsEmpty())
                masks[current] += line + "\n";
        }
    }
    t94File.Close();
 
    for (const auto& [key, value] : masks) {
        masks[key].RemoveLast(); // remove last newline character
    }
 
    // Load the FOREIGN MODELS from PSPMOD.mod file
    wxTextFile pspmodFile(Paths::GetDataPath() + "/atp/PSPMOD.mod");
    wxString pspmodString = "";
    if (pspmodFile.Open()) {
        pspmodString = pspmodFile.GetFirstLine();
        while (!pspmodFile.Eof())
        {
            pspmodString += pspmodFile.GetNextLine() + "\n";
        }
    }
    pspmodFile.Close();
 
    masks["BRANCH"].Replace("T94L_", "PSPNC");
    masks["BRANCH"].Replace("<STEP>", ATPField(m_data.pspStepSize, 6));
 
    masks["SWITCH"].Replace("T94L_", "PSPNC");
 
    //wxString switchMask = "  PSPNC%c%s%c                                        MEASURING                %d";
    //wxString sourceMask = "14PSPNC%c  %s%s%s                           -1.      100.";
 
    wxString lineStr = "";
    int branchCardPos = -1, switchCardPos = -1, sourceCardPos = -1, modelsCardPos = -1, outputCardPos = -1;
    int blankBranchCardPos = -1, blankModelsCardPos = -1;
    lineStr = atpFile.GetFirstLine();
    while (!atpFile.Eof())
    {
        if (lineStr.IsEmpty()) {
            lineStr = atpFile.GetNextLine();
            continue;
        }
        if (tolower(lineStr[0]) == 'c') {
            lineStr = atpFile.GetNextLine();
            continue; // Skip comments
        }
 
        if (lineStr.Find("/SWITCH") != wxNOT_FOUND) {
            switchCardPos = atpFile.GetCurrentLine() + 1;
        }
        if (lineStr.Find("/SOURCE") != wxNOT_FOUND) {
            sourceCardPos = atpFile.GetCurrentLine() + 1;
        }
        if (lineStr.Find("/BRANCH") != wxNOT_FOUND) {
            branchCardPos = atpFile.GetCurrentLine() + 1;
        }
        if (lineStr.Find("/MODELS") != wxNOT_FOUND) {
            modelsCardPos = atpFile.GetCurrentLine() + 1;
        }
        if (lineStr.Find("/OUTPUT") != wxNOT_FOUND) {
            outputCardPos = atpFile.GetCurrentLine() + 1;
        }
        if (lineStr.Find("BLANK BRANCH") != wxNOT_FOUND) {
            blankBranchCardPos = atpFile.GetCurrentLine() + 1;
        }
        if (lineStr.Find("BLANK MODELS") != wxNOT_FOUND) {
            blankModelsCardPos = atpFile.GetCurrentLine() + 1;
        }
 
        lineStr = atpFile.GetNextLine();
    }
 
    if (branchCardPos < 0) return false;
    if (outputCardPos < 0) return false;
 
    if (modelsCardPos < 0 && branchCardPos > 0) {
        // Models card must be the first card
        atpFile.InsertLine("/MODELS", branchCardPos - 1);
        atpFile.InsertLine("MODELS", branchCardPos);
        atpFile.InsertLine("ENDMODELS", branchCardPos + 1);
        modelsCardPos = branchCardPos - 1;
        outputCardPos += 3;
        blankBranchCardPos += 3;
        if (blankModelsCardPos < 0)
        {
            atpFile.InsertLine("BLANK MODELS", blankBranchCardPos - 1);
            blankModelsCardPos = blankBranchCardPos;
            blankBranchCardPos++;
 
        }
    }
 
    if (switchCardPos < 0 && outputCardPos > 0) {
        atpFile.InsertLine("/SWITCH", outputCardPos - 1);
        switchCardPos = outputCardPos;
        outputCardPos++;
        blankBranchCardPos++;
    }
 
    if (sourceCardPos < 0 && outputCardPos > 0) {
        atpFile.InsertLine("/SOURCE", outputCardPos - 1);
        sourceCardPos = outputCardPos;
        outputCardPos++;
        blankBranchCardPos++;
    }
 
    //for (char i = 'A'; i <= 'C'; ++i) {
    //  lineStr = wxString::Format(switchMask, i, node, i, i == 'A' ? 1 : 0);
    //  atpFile.InsertLine(lineStr, switchCardPos + (i - 'A'));
    //}
    //sourceCardPos += 3;
 
    double voltage = std::abs(m_data.puVoltage) * m_data.baseVoltage * (std::sqrt(2) / std::sqrt(3)); // phase-ground, peak value
    //wxString ampl = wxString::FromCDouble(voltage, 3); // Amplitude in pu
    //wxString freq = wxString::FromCDouble(m_data.frequency, 5); // Frequency in Hz
    // Insert spaces before to complete 10 characters
    //while (freq.Length() < 10) freq = " " + freq;
    //while (ampl.Length() < 10) ampl = " " + ampl;
    for (char i = 'A'; i <= 'C'; ++i) {
        wxString tag = wxString::Format("<AMP.%c   >", i);
        masks["SOURCE"].Replace(tag, ATPField(voltage, 10));
 
        tag = wxString::Format("<FREQ.%c  >", i);
        masks["SOURCE"].Replace(tag, ATPField(m_data.frequency, 10));
 
        double angle = std::arg(m_data.puVoltage) * 180.0 / M_PI - 120.0 * (i - 'A'); // Angle in degrees
        tag = wxString::Format("<PHASE.%c >", i);
        masks["SOURCE"].Replace(tag, ATPField(angle, 10));
    }
 
    lineStr = atpFile.GetFirstLine();
    while (!atpFile.Eof())
    {
        if (lineStr.IsEmpty()) {
            lineStr = atpFile.GetNextLine();
            continue;
        }
        if (tolower(lineStr[0]) == 'c') {
            lineStr = atpFile.GetNextLine();
            continue; // Skip comments
        }
 
        if (lineStr.Find(node) != wxNOT_FOUND) {
            lineStr.Replace(node, "PSPNC");
            atpFile.RemoveLine(atpFile.GetCurrentLine());
            atpFile.InsertLine(lineStr, atpFile.GetCurrentLine());
        }
 
        if (lineStr.Find("/MODELS") != wxNOT_FOUND)
            atpFile.InsertLine(pspmodString, atpFile.GetCurrentLine() + 2);
 
        if (lineStr.Find("/SWITCH") != wxNOT_FOUND)
            atpFile.InsertLine(masks["SWITCH"], atpFile.GetCurrentLine() + 1);
 
        if (lineStr.Find("/SOURCE") != wxNOT_FOUND)
            atpFile.InsertLine(masks["SOURCE"], atpFile.GetCurrentLine() + 1);
 
        if (lineStr.Find("/BRANCH") != wxNOT_FOUND)
            atpFile.InsertLine(masks["BRANCH"], atpFile.GetCurrentLine() + 1);
 
        lineStr = atpFile.GetNextLine();
    }
 
    lineStr = atpFile.GetFirstLine();
    while (!atpFile.Eof())
    {
        lineStr += atpFile.GetNextLine() + "\n";
 
    }
 
 
    //int cardPos = sourceCardPos;
    //for (char i = 'A'; i <= 'C'; ++i) {
    //  wxString angle = wxString::FromCDouble(std::arg(m_data.puVoltage) * 180.0 / M_PI - 120.0 * (i - 'A'), 5); // Angle in degrees
    //
    //  // Insert spaces before to complete 10 characters
    //  while (angle.Length() < 10) angle = " " + angle;
    //
    //  //lineStr = wxString::Format(sourceMask, i, ampl, freq, angle);
    //  atpFile.InsertLine(lineStr, cardPos);
    //  cardPos++;
    //
    //  if (hasBusData) {
    //      // Insert harmonic voltage
    //      for (size_t j = 0; j < busData.harmonicOrder.size(); ++j) {
    //          int order = busData.harmonicOrder[j];
    //          if (order == 1) continue; // Skip fundamental
    //
    //          std::complex<double> harmVoltage = busData.harmonicVoltage[j];
    //
    //          voltage = std::abs(busData.harmonicVoltage[j]) * m_data.baseVoltage * (std::sqrt(2) / std::sqrt(3));
    //          wxString amplH = wxString::FromCDouble(voltage, 3);
    //          while (amplH.Length() < 10) amplH = " " + amplH;
    //
    //          wxString freqH = wxString::FromCDouble(m_data.frequency * static_cast<double>(order), 5);
    //          while (freqH.Length() < 10) freqH = " " + freqH;
    //
    //          double angleValue = std::arg(harmVoltage) * 180.0 / M_PI;
    //          if (order % 3 == 1) { // Positive sequence
    //              angleValue += -120.0 * (i - 'A');
    //          }
    //          else if (order % 3 == 2) { // Negative sequence
    //              angleValue += 120.0 * (i - 'A');
    //          }
    //          // Zero sequence doesn't change the angle because it's the same for all phases
    //
    //          wxString angleH = wxString::FromCDouble(angleValue, 5);
    //          while (angleH.Length() < 10) angleH = " " + angleH;
    //
    //          //lineStr = wxString::Format(sourceMask, i, amplH, freqH, angleH);
    //          atpFile.InsertLine(lineStr, cardPos);
    //          cardPos++;
    //      }
    //  }
    //}
 
    return true;
}
 
void EMTElement::MedianFilter(std::vector<double>& data)
{
    // Based on: http://www.librow.com/articles/article-1
 
    const size_t n = data.size();
 
    if (n < 2)
        return;
 
    std::vector<double> extension(n + 4);
 
    std::memcpy(extension.data() + 2, data.data(), n * sizeof(double));
 
    extension[0] = data[1];
    extension[1] = data[0];
    extension[n + 2] = data[n - 1];
    extension[n + 3] = data[n - 2];
 
    DoMedianFilter(extension.data(), data, static_cast<int>(n + 4));
 
    //size_t n = data.size();
    //std::vector<double> result(data);
    //
    //if (data.empty() || n < 1)
    //  return result;
    //if (n == 1)
    //{
    //  result[0] = data[0];
    //  return result;
    //}
    //double* extension = new double[n + 4];
    //memcpy(extension + 2, data.data(), n * sizeof(double));
    //for (int i = 0; i < 2; ++i)
    //{
    //  extension[i] = data[1 - i];
    //  extension[n + 2 + i] = data[n - 1 - i];
    //}
    // DoMedianFilter(extension, result, n + 4);
    //delete[] extension;
    //
    //return result;
}
 
bool EMTElement::CalculateCurrent(wxString& errorMsg, const Mode& mode, const bool& saveRawData, const bool& saveFFTData)
{
    wxFileName fileName(m_data.atpFile);
    if (!fileName.IsOk()) {
        errorMsg = wxString::Format(_("Invalid ATP file path for the electromagnetic element \"%s\"."), m_data.name);
        return false;
    }
 
    //double fundFreq = properties.GetSimulationPropertiesData().stabilityFrequency;
    double fundFreq = m_data.frequency;
    double timeSim = (1.0 / fundFreq) * static_cast<double>(m_data.cyclesToSS);
 
    //wxString atpFolder = properties.GetGeneralPropertiesData().atpPath.GetPath();
    wxString atpFolder = m_data.atpPath.GetPath();
    // Set fundamental frequency to EMT element
    //m_data.frequency = fundFreq;
 
    wxExecuteEnv env;
    env.cwd = atpFolder;
 
    // Save the ATP file in ATP work folder
    wxTextFile origFile(m_data.atpFile.GetFullPath());
    //fileName.SetPath(properties.GetGeneralPropertiesData().atpWorkFolder);
    fileName.SetPath(m_data.atpWorkFolder);
    wxString name = fileName.GetFullPath();
    wxTextFile copyFile(fileName.GetFullPath());
    if (origFile.Open()) {
        copyFile.Create();
        for (size_t i = 0; i < origFile.GetLineCount(); i++) {
            copyFile.AddLine(origFile.GetLine(i));
        }
        origFile.Close();
 
        wxString stepSizeStr = wxString::Format("%.2E", m_data.stepSize);
        stepSizeStr.Replace(wxT(","), wxT("."));
        wxString timeStr = wxString::Format("%.6f", timeSim);
        timeStr.Replace(wxT(","), wxT("."));
        SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 1, 0, 8, stepSizeStr);
        SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 1, 8, 8, timeStr);
        SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 2, 0, 8, wxT("99999999"));
        SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 2, 8, 8, wxString::Format("%8d", m_data.recordFrequency));
 
        // Add source and switch to selected node
        AddConnectionToNode(copyFile, m_data.atpNodeName);
 
        copyFile.Write();
    }
    else {
        errorMsg = wxString::Format(_("Fail to open ATP file of the electromagnetic element \"%s\"."), m_data.name);
        return false;
    }
 
 
    ATPPSPBridge bridge;
 
    //wxProcess* proc = new wxProcess();
 
    //wxString cmd = properties.GetGeneralPropertiesData().atpPath.GetFullPath() + wxT(" both ") + fileName.GetFullPath() + wxT(" s -R");
    wxString cmd = m_data.atpPath.GetFullPath() + wxT(" both ") + fileName.GetFullPath() + wxT(" s -R");
    long pid = wxExecute(cmd, wxEXEC_ASYNC | wxEXEC_HIDE_CONSOLE | wxEXEC_NODISABLE, nullptr, &env);
    //long pid = wxExecute(cmd, wxEXEC_SYNC | wxEXEC_HIDE_CONSOLE | wxEXEC_NODISABLE, nullptr, &env);
    //long pid = wxExecute(cmd, wxEXEC_ASYNC | wxEXEC_NODISABLE, nullptr, &env);
    //long pid = wxExecute(cmd, wxEXEC_ASYNC | wxEXEC_NODISABLE, proc, &env);
 
    if (!bridge.Connect(5000, 100))
    {
        if (!CheckLISFile(fileName, errorMsg))
        {
            errorMsg = _("Unable to connect to ATP.\n"
                "Please make sure you are using an ATP version with PSP support.");
        }
 
        return false;
    }
    //wxEndBusyCursor();
    bridge.SetProcessId(pid);
 
 
    auto& data = bridge.GetSharedData();
    data.vrms = std::abs(m_data.puVoltage) * m_data.baseVoltage / std::sqrt(3);
    data.phase = std::arg(m_data.puVoltage) * 180.0 / M_PI;
    data.mode = mode;
    //m_data.atpData.clear();
    m_data.atpSampleData.clear();
 
    //bridge.ReleaseATP();
    while (bridge.WaitATP())
    {
        if (!saveRawData || !saveFFTData) {
            bridge.ReleaseATP();
            continue;
        }
 
        data = bridge.GetSharedData();
 
 
        // Process ATP data
        if (data.stepCount > 0)
        {
            const double dt = data.atpStepsize;
            const double t0 = data.t - (data.stepCount - 1) * dt;
 
            //m_data.atpData.reserve(m_data.atpData.size() + data.stepCount);
            m_data.atpSampleData.reserve(m_data.atpSampleData.size() + data.stepCount);
 
            std::vector<double> ia, ib, ic;
 
            if (m_data.useMedianFilter)
            {
                ia.reserve(data.stepCount);
                ib.reserve(data.stepCount);
                ic.reserve(data.stepCount);
 
                for (uint32_t i = 0; i < data.stepCount; ++i)
                {
                    ia.push_back(data.samples[i].current[0]);
                    ib.push_back(data.samples[i].current[1]);
                    ic.push_back(data.samples[i].current[2]);
                }
 
                MedianFilter(ia);
                MedianFilter(ib);
                MedianFilter(ic);
            }
 
            for (uint32_t i = 0; i < data.stepCount; ++i)
            {
                const double t = t0 + i * dt;
 
                const double currentA = m_data.useMedianFilter ? -ia[i] : -data.samples[i].current[0];
                const double currentB = m_data.useMedianFilter ? -ib[i] : -data.samples[i].current[1];
                const double currentC = m_data.useMedianFilter ? -ic[i] : -data.samples[i].current[2];
 
                const double iaMag = std::hypot(-data.phasor.Id[0], -data.phasor.Iq[0]);
                const double ibMag = std::hypot(-data.phasor.Id[1], -data.phasor.Iq[1]);
                const double icMag = std::hypot(-data.phasor.Id[2], -data.phasor.Iq[2]);
                const double iaAng = std::atan2(-data.phasor.Iq[0], -data.phasor.Id[0]);
                const double ibAng = std::atan2(-data.phasor.Iq[1], -data.phasor.Id[1]);
                const double icAng = std::atan2(-data.phasor.Iq[2], -data.phasor.Id[2]);
 
                //m_data.atpData.emplace_back(t, currentA);
                m_data.atpSampleData.push_back({
                    .t = t,
                    .current = { currentA, currentB, currentC },
                    .mag = {iaMag, ibMag, icMag},
                    .angle = {iaAng, ibAng, icAng}
                    });
            }
        }
 
        bridge.ReleaseATP();
    }
 
    if (data.terminate < 0)
    {
        CheckLISFile(fileName, errorMsg);
        return false;
    }
    constexpr size_t phase = 0; // 0 = A, 1 = B, 2 = C
 
    // Save current
    m_data.current = std::complex<double>(-data.phasor.Id[phase], -data.phasor.Iq[phase]);
 
    // Delete all .tmp files
    wxDir dir(atpFolder);
    if (dir.IsOpened()) {
        wxString file;
        bool cont = dir.GetFirst(&file, wxT("*.tmp"), wxDIR_FILES);
        while (cont) {
            wxRemoveFile(atpFolder + wxT("/") + file);
            cont = dir.GetNext(&file);
        }
    }
 
    if (saveFFTData)
    {
        // FFT
        double dataStepSize = m_data.atpSampleData[1].t - m_data.atpSampleData[0].t;
        bool useRemainder = false;
        size_t sampleCount = m_data.atpSampleData.size();
        size_t n = static_cast<size_t>(std::ceil(1.0 / (dataStepSize * fundFreq)));
 
        // Due to ATP logic, the number of samples can be greater than the number
        // of samples stored. In this case, follow the remainder approach.
        if (n > sampleCount)
        {
            n = sampleCount;
            useRemainder = true;
        }
 
        double fs = 1.0 / dataStepSize;
        double df = fs / static_cast<double>(n);
 
        if (useRemainder) {
            double rmder = std::abs(std::remainder(fundFreq, df));
            size_t minRmderN = n;
 
            while (rmder > 1e-3)
            {
                --n;
                df = fs / static_cast<double>(n);
                if (std::abs(std::remainder(fundFreq, df)) < rmder)
                    minRmderN = n;
 
                rmder = std::abs(std::remainder(fundFreq, df));
                if (n == 0) {
                    n = minRmderN;
                    df = fs / static_cast<double>(n);
                    break;
                }
            }
        }
 
        double dtWindow = m_data.atpSampleData.back().t - m_data.atpSampleData[sampleCount - n].t;
 
        fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (n / 2 + 1));
 
        double* in = (double*)fftw_malloc(sizeof(double) * n);
 
        // Populate FFT input
        const size_t first = sampleCount - n;
        for (size_t i = 0; i < n; ++i) {
            in[i] = m_data.atpSampleData[first + i].current[phase];
        }
 
        fftw_plan p = fftw_plan_dft_r2c_1d(n, in, out, FFTW_ESTIMATE);
 
        fftw_execute(p);
        fftw_destroy_plan(p);
 
        double ampCorrection = 2.0 / static_cast<double>(n);
 
        // Identify fundamental frequency
        int fundIndex = -1;
        double fundMagnitude = 0.0;
 
        for (size_t i = 0; i < (n / 2 + 1); ++i) {
            double freq = static_cast<double>(i) * df;
 
            if (std::remainder(freq, fundFreq)) {
                fundIndex = static_cast<int>(i);
 
                fundMagnitude =
                    std::hypot(out[i][0], out[i][1]) * ampCorrection;
 
                break;
            }
        }
 
        if (fundIndex < 0) {
            errorMsg = wxString::Format(_("Fail to identify the fundamental frequency of the electromagnetic element \"%s\"."), m_data.name);
 
            fftw_free(in);
            fftw_free(out);
 
            return false;
        }
 
        // Fundamental and harmonics
        m_data.currHarmonics.clear();
 
        for (size_t i = 0; i < (n / 2 + 1); ++i) {
            double freq = static_cast<double>(i) * df;
            int order = static_cast<int>(std::round(freq / fundFreq));
 
            if (order > m_data.numMaxHarmonics)
                break;
 
            double magnitude = std::hypot(out[i][0], out[i][1]) * ampCorrection;
 
            if ((magnitude / fundMagnitude) > (m_data.harmonicsThreshold / 100.0)) {
                m_data.currHarmonics[order] =
                    std::complex<double>(out[i][0], out[i][1]) *
                    ampCorrection / std::sqrt(2.0);
            }
        }
 
        m_data.inFFTData.clear();
        m_data.outFFTData.clear();
 
        for (size_t i = 0; i < n; ++i)
        {
            m_data.inFFTData.emplace_back(m_data.atpSampleData[first + i].t, in[i]);
        }
 
        for (size_t i = 0; i < (n / 2 + 1); ++i)
        {
            std::complex<double> value(out[i][0] * ampCorrection, out[i][1] * ampCorrection);
            double freq = static_cast<double>(i) * df;
            m_data.outFFTData.emplace_back(freq, value);
        }
 
        fftw_free(in);
        fftw_free(out);
    }
 
    return true;
 
    errorMsg = "Not implemented yet.\n EMTElement.cpp, line 772ish";
    return false;
 
    fileName.SetFullName(fileName.GetName() + wxT(".pl4"));
    wxFileInputStream pl4File(fileName.GetFullPath());
    // https://github.com/ldemattos/readPL4/wiki/PISA's-PL4-Format---Time-domain
    if (pl4File.IsOk() && pl4File.GetSize() > 0) {
        wxDataInputStream store(pl4File);
        store.UseBasicPrecisions();
        float* buffer = new float[2];
        std::vector<double> timeVec;
        std::vector<double> valueVec;
        double oldTime = 0.0;
        int read16Bytes = 0;
        while (!pl4File.Eof()) {
            store.ReadFloat(buffer, 2);
            if (read16Bytes > 11 && (buffer[0] >= oldTime)) {
                oldTime = buffer[0];
                timeVec.emplace_back(buffer[0]);
                valueVec.emplace_back(buffer[1]);
            }
            read16Bytes++;
        }
        delete[] buffer;
 
        if (m_data.useMedianFilter)
            MedianFilter(valueVec);
 
        // FFT
        double dataStepSize = timeVec[1] - timeVec[0]; // Real sampling time from data
        bool useRemainder = false;
        size_t n = ceil(1.0 / (dataStepSize * fundFreq)); // Number of samples
        // Due to ATP logic, the number of samples can be greater than the number of samples in the PL4 file
        // In this case, I choose to follow the remainder approach
        if (n > valueVec.size()) {
            n = valueVec.size();
            useRemainder = true;
        }
 
        double fs = 1.0 / dataStepSize; // Sampling frequency
        double df = fs / static_cast<double>(n); // Frequency resolution
 
        if (useRemainder) {
            double rmder = abs(remainder(fundFreq, df));
            int minRmderN = n;
            while (rmder > 1e-3) {
                n--;
                df = fs / static_cast<double>(n);
 
                if (abs(remainder(fundFreq, df)) < rmder) minRmderN = n;
 
                rmder = abs(remainder(fundFreq, df));
 
                if (n <= 0) {
                    n = minRmderN;
                    df = fs / static_cast<double>(n);
                    break;
                }
            }
        }
 
        double dtWindow = timeVec[timeVec.size() - 1] - timeVec[timeVec.size() - n];
 
        fftw_complex* out;
        double* in;
        fftw_plan p;
        out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (n / 2 + 1));
        in = (double*)fftw_malloc(sizeof(double) * n);
        // Populate input array
        for (size_t i = (valueVec.size() - n); i < valueVec.size(); ++i) {
            size_t index = i - (valueVec.size() - n);
            in[index] = valueVec[i];
        }
 
        p = fftw_plan_dft_r2c_1d(n, in, out, FFTW_ESTIMATE);
        fftw_execute(p); /* repeat as needed */
        fftw_destroy_plan(p);
 
        double ampCorrection = 2.0 / static_cast<double>(n);
 
        // Identify fundamental frequency
        int fundIndex = -1;
        double fundMagnitude = 0.0;
        for (size_t i = 0; i < (n / 2 + 1); i++) {
            double freq = static_cast<double>(i) * df;
            if (remainder(freq, fundFreq)) {
                fundIndex = i;
                fundMagnitude = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]) * ampCorrection;
                break;
            }
        }
        if (fundIndex < 0) {
            errorMsg = wxString::Format(_("Fail to identify the fundamental frequency of the electromagnetic element \"%s\"."), m_data.name);
            return false;
        }
 
        // Fundamental and harmonics
        m_data.currHarmonics.clear();
        //m_data.currHarmonicsOrder.clear();
        for (size_t i = 0; i < (n / 2 + 1); i++) {
            double freq = static_cast<double>(i) * df;
            int order = static_cast<int>(round(freq / fundFreq));
 
            if (order > m_data.numMaxHarmonics) break; // Stop if the order is greater than the maximum number of harmonics
 
            double magnitude = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]) * ampCorrection;
            if ((magnitude / fundMagnitude) > (m_data.harmonicsThreshold / 100.0)) {
                //m_data.currHarmonics.emplace_back(std::complex<double>(out[i][0], out[i][1]) * ampCorrection);
                //m_data.currHarmonicsOrder.emplace_back(order);
                m_data.currHarmonics[order] = std::complex<double>(out[i][0], out[i][1]) * ampCorrection / sqrt(2.0); // RMS value
            }
        }
 
        //if (saveFFTData) {
        //  m_data.atpData.clear();
        //  m_data.inFFTData.clear();
        //  m_data.outFFTData.clear();
        //  for (size_t i = 0; i < valueVec.size(); i++) {
        //      m_data.atpData.emplace_back(std::make_pair(timeVec[i], valueVec[i]));
        //  }
        //  for (size_t i = 0; i < n; i++) {
        //      m_data.inFFTData.emplace_back(std::make_pair(timeVec[i + timeVec.size() - n], in[i]));
        //  }
        //  for (size_t i = 0; i < (n / 2 + 1); i++) {
        //      std::complex<double> value(out[i][0] * ampCorrection, out[i][1] * ampCorrection);
        //      double freq = static_cast<double>(i) * df;
        //      m_data.outFFTData.emplace_back(std::make_pair(freq, value));
        //  }
        //}
        //
        fftw_free(in);
        fftw_free(out);
 
    }
    else {
        wxString atpErrorMsg = _("No error found.");
 
        fileName.SetFullName(fileName.GetName() + wxT(".lis"));
        wxTextFile lisFile(fileName.GetFullPath());
 
        if (!lisFile.Exists()) {
            errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP program did not return any error messages."), m_data.name);
            return false;
        }
        if (lisFile.Open()) {
            wxString line = lisFile.GetFirstLine() + "\n";
            bool foundError = false;
            int lineCount = -1;
            atpErrorMsg.Clear();
            while (!lisFile.Eof())
            {
                if (line.Find(wxT("KILL ")) != wxNOT_FOUND) {
                    foundError = true;
                    lineCount++;
                }
                if (line.Find(wxT("----------")) != wxNOT_FOUND) {
                    foundError = false;
                }
                if (foundError && lineCount == 1)
                    atpErrorMsg += line + " ";
                line = lisFile.GetNextLine();
            }
            lisFile.Close();
        }
        else
        {
            errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP program did not return any error messages."), m_data.name);
            return false;
        }
 
        errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP returned the following error message:\n\"%s\""), m_data.name, atpErrorMsg);
        return false;
    }
    return true;
}
 
void EMTElement::UpdateData(const PropertiesData* properties, bool updateVoltageBase)
{
    if (properties != nullptr) {
        m_data.frequency = properties->GetSimulationPropertiesData().stabilityFrequency;
        m_data.atpPath = properties->GetGeneralPropertiesData().atpPath;
        m_data.atpWorkFolder = properties->GetGeneralPropertiesData().atpWorkFolder;
        m_data.pspStepSize = properties->GetSimulationPropertiesData().timeStep;
    }
    if (!m_parentList.empty()) {
        Bus* bus = static_cast<Bus*>(m_parentList[0]);
        if (bus != nullptr) {
            auto busData = bus->GetElectricalData();
            std::complex<double> voltage = std::complex<double>(1.0, 0.0);
            if (abs(busData.voltage) > 1e-3)
                m_data.puVoltage = busData.voltage;
            if (updateVoltageBase) {
                m_data.baseVoltage = GetValueFromUnit(busData.nominalVoltage, busData.nominalVoltageUnit);
            }
        }
    }
 
}
 
void EMTElement::SetNominalVoltage(std::vector<double> nominalVoltage, std::vector<ElectricalUnit> nominalVoltageUnit)
{
    if (!nominalVoltage.empty()) {
        m_data.baseVoltage = GetValueFromUnit(nominalVoltage[0], nominalVoltageUnit[0]);
    }
}
 
void EMTElement::DoMedianFilter(double* extension, std::vector<double>& result, const int& n)
{
    //   Move window through all elements of the signal
    for (int i = 2; i < n - 2; ++i)
    {
        //   Pick up window elements
        double window[5];
        for (int j = 0; j < 5; ++j)
            window[j] = extension[i - 2 + j];
        //   Order elements (only half of them)
        for (int j = 0; j < 3; ++j)
        {
            //   Find position of minimum element
            int min = j;
            for (int k = j + 1; k < 5; ++k)
                if (window[k] < window[min])
                    min = k;
            //   Put found minimum element in its place
            std::swap(window[j], window[min]);
        }
        //   Get result - the middle element
        result[i - 2] = window[2];
    }
    //return result;
}
 
wxString EMTElement::ATPField(double value, size_t width)
{
    wxString best;
 
    // Try fixed-point notation
    for (int prec = 15; prec >= 0; --prec)
    {
        wxString s = wxString::Format("%.*f", prec, value);
        s.Replace(",", ".");
 
        while (s.Contains(".") && s.EndsWith("0"))
            s.RemoveLast();
        if (s.EndsWith("."))
            s.RemoveLast();
 
        if (s.Length() <= width)
        {
            best = s;
            break;
        }
    }
 
    // Try scientific notation
    for (int prec = 15; prec >= 0; --prec)
    {
        wxString s = wxString::Format("%.*E", prec, value);
        s.Replace(",", ".");
 
        if (prec == 0)
            s.Replace("E", ".E");
 
        if (s.Length() <= width)
        {
            // Keep the representation with the highest useful precision
            if (best.empty() || s.Length() > best.Length())
                best = s;
 
            break;
        }
    }
 
    if (best.empty())
        return wxString('*', width);
 
    best.Prepend(wxString(' ', width - best.Length()));
    return best;
}
 
bool EMTElement::CheckLISFile(wxFileName fileName, wxString& errorMsg) const
{
    wxString atpErrorMsg = _("No error found.");
 
    fileName.SetFullName(fileName.GetName() + wxT(".lis"));
    wxTextFile lisFile(fileName.GetFullPath());
 
    if (!lisFile.Exists()) {
        errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP program did not return any error messages."), m_data.name);
        return false;
    }
    if (lisFile.Open()) {
        wxString line = lisFile.GetFirstLine() + "\n";
        bool foundError = false;
        int lineCount = -1;
        atpErrorMsg.Clear();
        while (!lisFile.Eof())
        {
            if (line.Find(wxT("KILL ")) != wxNOT_FOUND) {
                foundError = true;
                lineCount++;
            }
            if (line.Find(wxT("----------")) != wxNOT_FOUND) {
                foundError = false;
            }
            if (foundError && lineCount == 1)
                atpErrorMsg += line + " ";
            line = lisFile.GetNextLine();
        }
        lisFile.Close();
    }
    else
    {
        errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP program did not return any error messages."), m_data.name);
        return true;
    }
 
    errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP returned the following error message:\n\"%s\""), m_data.name, atpErrorMsg);
    return true;
}