#include "nearToFarField.h"
#include <vtkQuad.h>
#include <vtkTriangle.h>
#include <vtkPointDataToCellData.h>
#include <vtkDoubleArray.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkGeometryFilter.h>
#include <vtkPolyDataWriter.h>
#include <vtkProbeFilter.h>
#include <vtkIdList.h>
#include <vtkExtractUnstructuredGrid.h>
#include <vtkPlane.h>
#include <vtkCutter.h>
#include "vtkForCenos.hpp"
#include "nlohmann/json.hpp"
#include <complex>
#include <boost/filesystem.hpp>


using json = nlohmann::json;

nearToFarField::nearToFarField() {
    radiated_power = 0;
    x_0 = 0;
    y_0 = 0;
    z_0 = 0;
    Nphi = 60;
    Ntheta = 30;
}

vtkSmartPointer<vtkUnstructuredGrid> nearToFarField::createSphere(double x_0, double y_0, double z_0, double r_sphere, int Ntheta, int Nphi)
{
    vtkSmartPointer<vtkPoints> points =
        vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkUnstructuredGrid> ugrid =
        vtkSmartPointer<vtkUnstructuredGrid>::New();
    const double PI = std::acos(-1);

    const double thetaRange = PI; //from 0 to PI
    const double phiRange = 2 * PI; //from 0 to 2*PI
    double dtheta = thetaRange / Ntheta;
    double dphi = phiRange / Nphi;

    points->InsertNextPoint(x_0, y_0, z_0 + r_sphere);
    phiValues.push_back(0.);
    thetaValues.push_back(0.);
    for (int i = 1; i < Ntheta; i++)
    {
        for (int j = 0; j < Nphi; j++)
        {
            double x = x_0 + r_sphere * sin(i * dtheta) * cos(j * dphi);
            double y = y_0 + r_sphere * sin(i * dtheta) * sin(j * dphi);
            double z = z_0 + r_sphere * cos(i * dtheta);
            points->InsertNextPoint(x, y, z);
            phiValues.push_back(j * dphi);
            thetaValues.push_back(i * dtheta);
        }
    }

    points->InsertNextPoint(x_0, y_0, z_0 - r_sphere);
    phiValues.push_back(0.);
    thetaValues.push_back(PI);
    ugrid->SetPoints(points);

    //top layer
    for (int j = 0; j < Nphi - 1; j++)
    {
        vtkSmartPointer<vtkTriangle> tri =
            vtkSmartPointer<vtkTriangle>::New();
        tri->GetPointIds()->SetId(0, 0);
        tri->GetPointIds()->SetId(1, j + 1);
        tri->GetPointIds()->SetId(2, j + 2);
        ugrid->InsertNextCell(tri->GetCellType(), tri->GetPointIds());
    }

    {
        vtkSmartPointer<vtkTriangle> tri =
            vtkSmartPointer<vtkTriangle>::New();
        tri->GetPointIds()->SetId(0, 0);
        tri->GetPointIds()->SetId(1, Nphi);
        tri->GetPointIds()->SetId(2, 1);
        ugrid->InsertNextCell(tri->GetCellType(), tri->GetPointIds());
    }


    //main part
    for (int i = 1; i < Ntheta - 1; i++)
    {
        for (int j = 0; j < Nphi - 1; j++)
        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();
            quad->GetPointIds()->SetId(0, (i - 1) * Nphi + j + 1);
            quad->GetPointIds()->SetId(1, (i - 1) * Nphi + j + 1 + 1);
            quad->GetPointIds()->SetId(2, i * Nphi + j + 1 + 1);
            quad->GetPointIds()->SetId(3, i * Nphi + j + 1);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }

        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();
            quad->GetPointIds()->SetId(0, (i - 1) * Nphi + 1);
            quad->GetPointIds()->SetId(1, (i - 1) * Nphi + Nphi);
            quad->GetPointIds()->SetId(2, i * Nphi + Nphi);
            quad->GetPointIds()->SetId(3, i * Nphi + 1);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }
    }


    //bottom layer
    for (int j = 0; j < Nphi - 1; j++)
    {
        vtkSmartPointer<vtkTriangle> tri =
            vtkSmartPointer<vtkTriangle>::New();
        tri->GetPointIds()->SetId(0, (Ntheta - 1) * Nphi + 1);
        tri->GetPointIds()->SetId(1, (Ntheta - 2) * Nphi + j + 1);
        tri->GetPointIds()->SetId(2, (Ntheta - 2) * Nphi + j + 2);
        ugrid->InsertNextCell(tri->GetCellType(), tri->GetPointIds());
    }

    {
        vtkSmartPointer<vtkTriangle> tri =
            vtkSmartPointer<vtkTriangle>::New();
        tri->GetPointIds()->SetId(0, (Ntheta - 1) * Nphi + 1);
        tri->GetPointIds()->SetId(1, (Ntheta - 2) * Nphi + Nphi);
        tri->GetPointIds()->SetId(2, (Ntheta - 2) * Nphi + 1);
        ugrid->InsertNextCell(tri->GetCellType(), tri->GetPointIds());
    }
    return ugrid;
}

vtkSmartPointer<vtkUnstructuredGrid> nearToFarField::createBox(double minx, double maxx, double miny, double maxy, double minz, double maxz, int Nx, int Ny, int Nz)
{
    double dx = (maxx - minx) / Nx;
    double dy = (maxy - miny) / Ny;
    double dz = (maxz - minz) / Nz;

    vtkSmartPointer<vtkDoubleArray> normalData =
        vtkSmartPointer<vtkDoubleArray>::New();
    normalData->SetNumberOfComponents(3);
    normalData->SetName("Normals");


    vtkSmartPointer<vtkPoints> points =
        vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkUnstructuredGrid> ugrid =
        vtkSmartPointer<vtkUnstructuredGrid>::New();

    int pointCounter = 0;
    std::vector<std::vector<int> > minx_plane;
    std::vector<std::vector<int> > miny_plane;
    std::vector<std::vector<int> > minz_plane;
    std::vector<std::vector<int> > maxx_plane;
    std::vector<std::vector<int> > maxy_plane;
    std::vector<std::vector<int> > maxz_plane;
    minx_plane.resize(Ny + 1);
    miny_plane.resize(Nx + 1);
    minz_plane.resize(Nx + 1);
    maxx_plane.resize(Ny + 1);
    maxy_plane.resize(Nx + 1);
    maxz_plane.resize(Nx + 1);

    for (int i = 0; i <= Ny; i++)
    {
        minx_plane[i].resize(Nz + 1);
        maxx_plane[i].resize(Nz + 1);
    }
    for (int i = 0; i <= Nx; i++)
    {
        miny_plane[i].resize(Nz + 1);
        maxy_plane[i].resize(Nz + 1);
    }
    for (int i = 0; i <= Nx; i++)
    {
        minz_plane[i].resize(Ny + 1);
        maxz_plane[i].resize(Ny + 1);
    }



    //minz plane
    for (int i = 0; i <= Nx; i++)
    {
        for (int j = 0; j <= Ny; j++)
        {
            points->InsertNextPoint(minx + i * dx, miny + j * dy, minz);
            minz_plane[i][j] = pointCounter;
            pointCounter++;

        }
    }

    // miny plane
    for (int i = 0; i <= Nx; i++)
    {
        for (int k = 0; k <= Nz; k++)
        {
            points->InsertNextPoint(minx + i * dx, miny, minz + k * dz);
            miny_plane[i][k] = pointCounter;
            pointCounter++;

        }
    }

    //minx plane
    for (int j = 0; j <= Ny; j++)
    {
        for (int k = 0; k <= Nz; k++)
        {
            points->InsertNextPoint(minx, miny + j * dy, minz + k * dz);
            minx_plane[j][k] = pointCounter;
            pointCounter++;
        }
    }


    //maxz plane
    for (int i = 0; i <= Nx; i++)
    {
        for (int j = 0; j <= Ny; j++)
        {
            points->InsertNextPoint(minx + i * dx, miny + j * dy, maxz);
            maxz_plane[i][j] = pointCounter;
            pointCounter++;
        }
    }

    // maxy plane
    for (int i = 0; i <= Nx; i++)
    {
        for (int k = 0; k <= Nz; k++)
        {
            points->InsertNextPoint(minx + i * dx, maxy, minz + k * dz);
            maxy_plane[i][k] = pointCounter;
            pointCounter++;
        }
    }

    //maxx plane
    for (int j = 0; j <= Ny; j++)
    {
        for (int k = 0; k <= Nz; k++)
        {
            points->InsertNextPoint(maxx, miny + j * dy, minz + k * dz);
            maxx_plane[j][k] = pointCounter;
            pointCounter++;
        }
    }

    ugrid->SetPoints(points);


    //minz plane
    for (int i = 0; i < Nx; i++)
    {
        for (int j = 0; j < Ny; j++)
        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();
            quad->GetPointIds()->SetId(3, minz_plane[i][j]);
            quad->GetPointIds()->SetId(2, minz_plane[i + 1][j]);
            quad->GetPointIds()->SetId(1, minz_plane[i + 1][j + 1]);
            quad->GetPointIds()->SetId(0, minz_plane[i][j + 1]);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }
    }

    // miny plane
    for (int i = 0; i < Nx; i++)
    {
        for (int k = 0; k < Nz; k++)
        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();
            quad->GetPointIds()->SetId(0, miny_plane[i][k]);
            quad->GetPointIds()->SetId(1, miny_plane[i + 1][k]);
            quad->GetPointIds()->SetId(2, miny_plane[i + 1][k + 1]);
            quad->GetPointIds()->SetId(3, miny_plane[i][k + 1]);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }
    }

    //minx plane
    for (int j = 0; j < Ny; j++)
    {
        for (int k = 0; k < Nz; k++)
        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();
            quad->GetPointIds()->SetId(3, minx_plane[j][k]);
            quad->GetPointIds()->SetId(2, minx_plane[j + 1][k]);
            quad->GetPointIds()->SetId(1, minx_plane[j + 1][k + 1]);
            quad->GetPointIds()->SetId(0, minx_plane[j][k + 1]);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }
    }

    //maxz plane
    for (int i = 0; i < Nx; i++)
    {
        for (int j = 0; j < Ny; j++)
        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();
            quad->GetPointIds()->SetId(0, maxz_plane[i][j]);
            quad->GetPointIds()->SetId(1, maxz_plane[i + 1][j]);
            quad->GetPointIds()->SetId(2, maxz_plane[i + 1][j + 1]);
            quad->GetPointIds()->SetId(3, maxz_plane[i][j + 1]);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }
    }

    // maxy plane
    for (int i = 0; i < Nx; i++)
    {
        for (int k = 0; k < Nz; k++)
        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();
            quad->GetPointIds()->SetId(3, maxy_plane[i][k]);
            quad->GetPointIds()->SetId(2, maxy_plane[i + 1][k]);
            quad->GetPointIds()->SetId(1, maxy_plane[i + 1][k + 1]);
            quad->GetPointIds()->SetId(0, maxy_plane[i][k + 1]);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }
    }

    //maxx plane
    for (int j = 0; j < Ny; j++)
    {
        for (int k = 0; k < Nz; k++)
        {
            vtkSmartPointer<vtkQuad> quad =
                vtkSmartPointer<vtkQuad>::New();

            quad->GetPointIds()->SetId(0, maxx_plane[j][k]);
            quad->GetPointIds()->SetId(1, maxx_plane[j + 1][k]);
            quad->GetPointIds()->SetId(2, maxx_plane[j + 1][k + 1]);
            quad->GetPointIds()->SetId(3, maxx_plane[j][k + 1]);
            ugrid->InsertNextCell(quad->GetCellType(), quad->GetPointIds());
        }
    }
    return ugrid;
}

void nearToFarField::readConfig()
{
    char const* tmp = getenv("CENOS_CONFIG");
    if (!tmp) {
        std::cout << "variable CENOS_CONFIG not set " << "\n";
    }

    std::string ntffConfigJson = std::string(tmp) + "\\jsons\\ntff.json";
    json j;

    if (boost::filesystem::exists(ntffConfigJson))
    {
        std::ifstream infile(ntffConfigJson, std::ifstream::in);
        if (infile.is_open()) {
            std::stringstream strStream;
            strStream << infile.rdbuf();
            std::string str = strStream.str();
            std::istringstream iss(str);
            iss >> j;
            Ntheta = j["Ntheta"].get<int>();
            Nphi = j["Nphi"].get<int>();
        }
        infile.close();
    }

}

std::vector<vtkSmartPointer<vtkUnstructuredGrid>> nearToFarField::execute()
{
    const double PI = 3.14159;
    readConfig();

    vtkSmartPointer<vtkDataSet> enclosedData = this->getNearFieldData(enclosingGrids);
    std::vector<double> bbox = getBoundingBox(enclosingGrids);

    //TODO find x_0, y_0, z_0, r_sphere
    x_0 = (bbox[1] + bbox[0]) / 2;
    y_0 = (bbox[3] + bbox[2]) / 2;
    z_0 = (bbox[5] + bbox[4]) / 2;
    double r_sphere = 1;
    vtkSmartPointer<vtkUnstructuredGrid> ugrid = createSphere(x_0, y_0, z_0, r_sphere, Nphi, Ntheta);


    vtkSmartPointer<vtkPointDataToCellData> pointToCell =
        vtkSmartPointer<vtkPointDataToCellData>::New();
    pointToCell->SetInputData(enclosedData);
    pointToCell->PassPointDataOn();
    pointToCell->Update();

    vtkSmartPointer<vtkUnstructuredGrid> sourceGridwithCellsData = pointToCell->GetUnstructuredGridOutput();

    vtkSmartPointer<vtkDoubleArray> normals = vtkSmartPointer<vtkDoubleArray>::New();
    normals->SetNumberOfComponents(3);
    normals->SetName("Normals");

    vtkSmartPointer<vtkDoubleArray> areas = vtkSmartPointer<vtkDoubleArray>::New();
    areas->SetNumberOfComponents(1);
    areas->SetName("Areas");

    const vtkIdType numSourceCells = sourceGridwithCellsData->GetNumberOfCells();
    for (vtkIdType i = 0; i < numSourceCells; i++)
    {

        vtkSmartPointer<vtkCell>  c = sourceGridwithCellsData->GetCell(i);
        double area = 0;
        if (c->IsA("vtkTriangle"))
        {
            vtkSmartPointer<vtkPoints> pts = c->GetPoints();
            double normal[3] = { 0, 0, 0 };
            double pt0[3], pt1[3], pt2[3];
            pts->GetPoint(0, pt0);
            pts->GetPoint(1, pt1);
            pts->GetPoint(2, pt2);
            std::vector<double> v1 = { pt1[0] - pt0[0], pt1[1] - pt0[1], pt1[2] - pt0[2] };
            std::vector<double> v2 = { pt2[0] - pt0[0], pt2[1] - pt0[1], pt2[2] - pt0[2] };
            area = 0.5 * sqrt(pow((v1[1] * v2[2] - v1[2] * v2[1]), 2)
                + pow((v1[2] * v2[0] - v1[0] * v2[2]), 2)
                + pow((v1[0] * v2[1] - v1[1] * v2[0]), 2));

            vtkTriangle::ComputeNormal(pt0, pt1, pt2, normal);
            normals->InsertTuple3(i, normal[0], normal[1], normal[2]);
            areas->InsertTuple1(i, area);
        }
        else if (c->IsA("vtkQuad"))
        {
            vtkSmartPointer<vtkPoints> pts = c->GetPoints();
            double normal[3] = { 0, 0, 0 };
            double pt0[3], pt1[3], pt2[3], pt3[3];
            pts->GetPoint(0, pt0);
            pts->GetPoint(1, pt1);
            pts->GetPoint(2, pt2);
            pts->GetPoint(3, pt3);
            std::vector<double> v1 = { pt1[0] - pt0[0], pt1[1] - pt0[1], pt1[2] - pt0[2] };
            std::vector<double> v2 = { pt2[0] - pt0[0], pt2[1] - pt0[1], pt2[2] - pt0[2] };
            std::vector<double> v3 = { pt3[0] - pt0[0], pt3[1] - pt0[1], pt3[2] - pt0[2] };

            double area1 = 0.5 * sqrt(pow((v1[1] * v2[2] - v1[2] * v2[1]), 2)
                + pow((v1[2] * v2[0] - v1[0] * v2[2]), 2)
                + pow((v1[0] * v2[1] - v1[1] * v2[0]), 2));
            double area2 = 0.5 * sqrt(pow((v3[1] * v2[2] - v3[2] * v2[1]), 2)
                + pow((v3[2] * v2[0] - v3[0] * v2[2]), 2)
                + pow((v3[0] * v2[1] - v3[1] * v2[0]), 2));
            area = area1 + area2;

            vtkTriangle::ComputeNormal(pt0, pt1, pt2, normal);
            normals->InsertTuple3(i, normal[0], normal[1], normal[2]);
            areas->InsertTuple1(i, area);
        }
    }

    sourceGridwithCellsData->GetCellData()->AddArray(normals);
    sourceGridwithCellsData->GetCellData()->AddArray(areas);

    int nCELLS = sourceGridwithCellsData->GetNumberOfCells();

    std::vector<std::vector<std::complex<double>> > magneticCurrent;
    std::vector<std::vector<std::complex<double>> > electricCurrent;
    std::vector<std::vector<std::complex<double>> > magneticField;
    std::vector<std::vector<std::complex<double>> > electricField;
    magneticCurrent.resize(nCELLS);
    electricCurrent.resize(nCELLS);    
    magneticField.resize(nCELLS);
    electricField.resize(nCELLS);
    for (vtkIdType el_i = 0; el_i < nCELLS; el_i++)
    {
        magneticCurrent[el_i].resize(3);
        electricCurrent[el_i].resize(3);
        magneticField[el_i].resize(3);
        electricField[el_i].resize(3);
    }

    /// calculate magnetic and electric currents

    vtkSmartPointer<vtkDataArray> h_re = sourceGridwithCellsData->GetCellData()->GetArray("Magnetic_Field_re,_[A/m]");
    vtkSmartPointer<vtkDataArray> h_im = sourceGridwithCellsData->GetCellData()->GetArray("Magnetic_Field_im,_[A/m]");

    vtkSmartPointer<vtkDataArray> e_re = sourceGridwithCellsData->GetCellData()->GetArray("Electric_Field_re,_[V/m]");
    vtkSmartPointer<vtkDataArray> e_im = sourceGridwithCellsData->GetCellData()->GetArray("Electric_Field_im,_[V/m]");
    vtkSmartPointer<vtkDataArray> normal = sourceGridwithCellsData->GetCellData()->GetArray("Normals");


    vtkSmartPointer<vtkDoubleArray> elcurr = vtkSmartPointer<vtkDoubleArray>::New();
    elcurr->SetNumberOfComponents(3);
    elcurr->SetName("El_current");

    vtkSmartPointer<vtkDoubleArray> magcurr = vtkSmartPointer<vtkDoubleArray>::New();
    magcurr->SetNumberOfComponents(3);
    magcurr->SetName("Mag_current");

    for (vtkIdType el_i = 0; el_i < nCELLS; el_i++)
    {
        vtkSmartPointer<vtkCell>  c = sourceGridwithCellsData->GetCell(el_i);
        double normDir[3];
        normal->GetTuple(el_i, normDir);
        std::vector<std::complex<double>> h_field = std::vector<std::complex<double>>({
            std::complex<double>(h_re->GetTuple(el_i)[0], h_im->GetTuple(el_i)[0]),
            std::complex<double>(h_re->GetTuple(el_i)[1], h_im->GetTuple(el_i)[1]),
            std::complex<double>(h_re->GetTuple(el_i)[2], h_im->GetTuple(el_i)[2]) });

        //magnetic current = n x h
        electricCurrent[el_i][0] = normDir[1] * h_field[2] - normDir[2] * h_field[1];
        electricCurrent[el_i][1] = normDir[2] * h_field[0] - normDir[0] * h_field[2];
        electricCurrent[el_i][2] = normDir[0] * h_field[1] - normDir[1] * h_field[0];
        magneticField[el_i] = { h_field[0], h_field[1], h_field[2] };

        elcurr->InsertNextTuple3(imag(electricCurrent[el_i][0]), imag(electricCurrent[el_i][1]), imag(electricCurrent[el_i][2]));

        std::vector<std::complex<double>> e_field = std::vector<std::complex<double>>({
            std::complex<double>(e_re->GetTuple(el_i)[0], e_im->GetTuple(el_i)[0]),
            std::complex<double>(e_re->GetTuple(el_i)[1], e_im->GetTuple(el_i)[1]),
            std::complex<double>(e_re->GetTuple(el_i)[2], e_im->GetTuple(el_i)[2]) });

        //electric current = -n x e
        magneticCurrent[el_i][0] = normDir[2] * e_field[1] - normDir[1] * e_field[2];
        magneticCurrent[el_i][1] = normDir[0] * e_field[2] - normDir[2] * e_field[0];
        magneticCurrent[el_i][2] = normDir[1] * e_field[0] - normDir[0] * e_field[1];
        electricField[el_i] = { e_field[0], e_field[1], e_field[2] };
        magcurr->InsertNextTuple3( real(magneticCurrent[el_i][0]), real(magneticCurrent[el_i][1]), real(magneticCurrent[el_i][2]) );

    }

    sourceGridwithCellsData->GetCellData()->AddArray(elcurr);
    sourceGridwithCellsData->GetCellData()->AddArray(magcurr);

 /*
    //uncomment to write out test_block.vtk
    //It contains data interpolated on the bounding box
    vtkSmartPointer<vtkGeometryFilter> geometryFilter1 =
        vtkSmartPointer<vtkGeometryFilter>::New();

    vtkSmartPointer<vtkPolyDataWriter> polywriter1 =
        vtkSmartPointer<vtkPolyDataWriter>::New();

    geometryFilter1->SetInputData(sourceGridwithCellsData);
    geometryFilter1->Update();

    std::string filename1 = path + "/test_block.vtk";
    polywriter1->SetFileName(filename1.c_str());
    polywriter1->SetInputData(geometryFilter1->GetOutput());
    polywriter1->Update();
    polywriter1->Write();
    */

    /// finished calculate currents


    vtkSmartPointer<vtkDoubleArray> farField = vtkSmartPointer<vtkDoubleArray>::New();
    farField->SetNumberOfComponents(1);
    farField->SetName("Far_Field_rE_[V]");

    vtkSmartPointer<vtkDoubleArray> radIntensity = vtkSmartPointer<vtkDoubleArray>::New();
    radIntensity->SetNumberOfComponents(1);
    radIntensity->SetName("Radiation_Intensity_U_[W/srd]");

    vtkSmartPointer<vtkDoubleArray> directivity = vtkSmartPointer<vtkDoubleArray>::New();
    directivity->SetNumberOfComponents(1);
    directivity->SetName("Directivity_D_[dB]");

    //insert values
    radiated_power = 0;
    for (int i = 0; i < phiValues.size(); i++)
    {
        std::vector<double> data = getFarField(magneticCurrent, electricCurrent, k0, std::vector<double>({ thetaValues[i], phiValues[i] }), r_sphere, sourceGridwithCellsData);
        farField->InsertNextValue(data[0]);
        radIntensity->InsertNextValue(data[1]);
        directivity->InsertNextValue(data[2]);
        radiated_power += data[1] * sin(thetaValues[i]) * PI / Ntheta * 2 * PI / Nphi;
    }

    double tuple[1];
    for (vtkIdType t = 0; t < phiValues.size(); t++)
    {
        directivity->GetTuple(t, tuple);
        directivity->SetTuple1(t, 10*log10(4 * PI * tuple[0] / radiated_power));
    }



    vtkSmartPointer<vtkPoints> pts = ugrid->GetPoints();
    
         
    ugrid->GetPointData()->AddArray(farField);
    ugrid->GetPointData()->AddArray(radIntensity);
    ugrid->GetPointData()->AddArray(directivity);

    writePlaneFields(ugrid, "XY");
    writePlaneFields(ugrid, "YZ");
    writePlaneFields(ugrid, "ZX");

    vtkSmartPointer<vtkUnstructuredGrid> ugrid_rad = vtkSmartPointer<vtkUnstructuredGrid>::New();
    ugrid_rad->DeepCopy(ugrid);
    vtkSmartPointer<vtkUnstructuredGrid> ugrid_dir = vtkSmartPointer<vtkUnstructuredGrid>::New();
    ugrid_dir->DeepCopy(ugrid);

    vtkSmartPointer<vtkDoubleArray> blockId1 = vtkSmartPointer<vtkDoubleArray>::New();
    blockId1->SetNumberOfComponents(1);
    blockId1->SetName("BlockId");
    vtkSmartPointer<vtkDoubleArray> blockId2 = vtkSmartPointer<vtkDoubleArray>::New();
    blockId2->SetNumberOfComponents(1);
    blockId2->SetName("BlockId");
    vtkSmartPointer<vtkDoubleArray> blockId3 = vtkSmartPointer<vtkDoubleArray>::New();
    blockId3->SetNumberOfComponents(1);
    blockId3->SetName("BlockId");

    for (int i = 0; i < ugrid->GetNumberOfCells(); i++)
    {
        blockId1->InsertNextValue(1.0);
        blockId2->InsertNextValue(2.0);
        blockId3->InsertNextValue(3.0);
    }


    double pt[3];
    //get min, max of field
    double minMax[2];
    farField->GetRange(minMax);
    double dmm = minMax[1] - minMax[0];
    for (vtkIdType t = 0; t < phiValues.size(); t++)
    {
        farField->GetTuple(t, tuple);
        double scaleFact = (tuple[0] - minMax[0]) / dmm;
        pt[0] = x_0 + scaleFact * r_sphere * sin(thetaValues[t]) * cos(phiValues[t]);
        pt[1] = y_0 + scaleFact * r_sphere * sin(thetaValues[t]) * sin(phiValues[t]);
        pt[2] = z_0 + scaleFact * r_sphere * cos(thetaValues[t]);
        ugrid->GetPoints()->SetPoint(t, pt);
    }
    ugrid->GetCellData()->AddArray(blockId1);

    //radiant intensity
    radIntensity->GetRange(minMax);
    dmm = minMax[1] - minMax[0];
    for (vtkIdType t = 0; t < phiValues.size(); t++)
    {
        radIntensity->GetTuple(t, tuple);
        double scaleFact = (tuple[0] - minMax[0]) / dmm;
        pt[0] = x_0 + scaleFact * r_sphere * sin(thetaValues[t]) * cos(phiValues[t]);
        pt[1] = y_0 + scaleFact * r_sphere * sin(thetaValues[t]) * sin(phiValues[t]);
        pt[2] = z_0 + scaleFact * r_sphere * cos(thetaValues[t]);
        ugrid_rad->GetPoints()->SetPoint(t, pt);
    }
    ugrid_rad->GetCellData()->AddArray(blockId2);

    //directivity
    directivity->GetRange(minMax);
    dmm = minMax[1] - minMax[0];
    for (vtkIdType t = 0; t < phiValues.size(); t++)
    {
        directivity->GetTuple(t, tuple);
        double scaleFact = (tuple[0] - minMax[0]) / dmm;
        pt[0] = x_0 + scaleFact * r_sphere * sin(thetaValues[t]) * cos(phiValues[t]);
        pt[1] = y_0 + scaleFact * r_sphere * sin(thetaValues[t]) * sin(phiValues[t]);
        pt[2] = z_0 + scaleFact * r_sphere * cos(thetaValues[t]);
        ugrid_dir->GetPoints()->SetPoint(t, pt);
    }
    ugrid_dir->GetCellData()->AddArray(blockId3);



    
    //uncomment to write out farField.vtk
    vtkSmartPointer<vtkGeometryFilter> geometryFilter =
        vtkSmartPointer<vtkGeometryFilter>::New();

    vtkSmartPointer<vtkPolyDataWriter> polywriter =
        vtkSmartPointer<vtkPolyDataWriter>::New();

    geometryFilter->SetInputData(ugrid_dir);
    geometryFilter->Update();


    std::string filename = path + "/farField.vtk";
    polywriter->SetFileName(filename.c_str());
    polywriter->SetInputData(geometryFilter->GetOutput());
    polywriter->Update();
    polywriter->Write();
    
    std::vector<vtkSmartPointer<vtkUnstructuredGrid>> returnVec;
    returnVec.push_back(ugrid);
    returnVec.push_back(ugrid_rad);
    returnVec.push_back(ugrid_dir);
    return returnVec;
}


void nearToFarField::writePlaneFields(vtkSmartPointer<vtkUnstructuredGrid> ugrid, std::string planeName)
{
    vtkSmartPointer<vtkPlane> plane =
        vtkSmartPointer<vtkPlane>::New();
    std::string angleName;
    plane->SetOrigin(x_0, y_0, z_0);
    if ((planeName == "XY") || (planeName == "YX"))
    {
        plane->SetNormal(0, 0, 1);
        angleName = "phi";
    }
    else if ((planeName == "XZ") || (planeName == "ZX"))
    {
        plane->SetNormal(0, 1, 0);
    }
    else if ((planeName == "YZ") || (planeName == "ZY")){
        plane->SetNormal(1, 0, 0);
        angleName = "theta";}
    else{
        return;}

    // Create cutter
    vtkSmartPointer<vtkCutter> cutter = vtkSmartPointer<vtkCutter>::New();
    cutter->SetCutFunction(plane);
    cutter->SetInputData(ugrid);
    cutter->Update();
    vtkSmartPointer<vtkPolyData> cutterPolyData = cutter->GetOutput();
    
    for (int j = 0; j < cutterPolyData->GetPointData()->GetNumberOfArrays(); j++) // Write a sperate csv file for each paramter
    {
        vtkSmartPointer<vtkDataArray> dataArray = cutterPolyData->GetPointData()->GetArray(j);
        std::string csvName = path + "/plane" + planeName + "_" + dataArray->GetName() + ".csv";
        std::ofstream planeCsv;
        double PI = acos(-1);

        // get theta data for specific plane
        int nrofpoints = cutterPolyData->GetNumberOfPoints();
        std::vector<std::pair<double, int>> thetaAngles;
        for (int i = 0; i < nrofpoints; i++)
        {
            if ((planeName == "XY") || (planeName == "YX"))
                thetaAngles.push_back(std::make_pair(atan2(cutterPolyData->GetPoint(i)[1] - y_0, cutterPolyData->GetPoint(i)[0] - x_0) * 180 / PI, i));
            else if ((planeName == "XZ") || (planeName == "ZX"))
                thetaAngles.push_back(std::make_pair(atan2(cutterPolyData->GetPoint(i)[0] - x_0, cutterPolyData->GetPoint(i)[2] - z_0) * 180 / PI, i));
            else if ((planeName == "YZ") || (planeName == "ZY"))
                thetaAngles.push_back(std::make_pair(atan2(cutterPolyData->GetPoint(i)[1] - y_0, cutterPolyData->GetPoint(i)[2] - z_0) * 180 / PI, i));
        }
        std::sort(thetaAngles.begin(), thetaAngles.end());


        std::string outString;

        if (isFirstFrequency == true) { //create new file and write "time" and all theta values in degrees
            planeCsv.open(csvName);
            outString.append("time");
            for (int i = 0; i < nrofpoints; i++){ //degree values
                outString.append(", " + std::to_string(thetaAngles[i].first));
            }
            outString.append("\n");

            outString.append(std::to_string(f0));
            for (int i = 0; i < dataArray->GetNumberOfValues(); i++) {
                outString.append(", " + std::to_string(dataArray->GetTuple1(thetaAngles[i].second)));
            }
            outString.append("\n");
        }
        else {
            planeCsv.open(csvName, std::ios_base::app); //add line with frequency and parameter values
            outString.append(std::to_string(f0));
            for (int i = 0; i < dataArray->GetNumberOfValues(); i++){
                outString.append(", " + std::to_string(dataArray->GetTuple1(thetaAngles[i].second)));
            }
            outString.append("\n");
        }
        planeCsv << outString;
        planeCsv.close();
    }
}

void nearToFarField::setFirstFreq(bool b)
{
    isFirstFrequency = b;
}

std::vector<double> nearToFarField::getFarField(std::vector<std::vector<std::complex<double>> > magneticCurrent, std::vector<std::vector<std::complex<double>> > electricCurrent,
    double k0, std::vector<double> angles, double radius, vtkSmartPointer<vtkUnstructuredGrid> sourceGridwithCellsData)
{
    double Z0 = 377.0; //free space impedance

    const double pi_ = std::acos(-1);
    const std::complex<double> i_(0, 1);

    double theta = angles[0];
    double phi = angles[1];
    std::vector<double> coords_s = { sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta) };

    std::complex<double>  N_theta = { 0,0 };
    std::complex<double> N_phi = { 0,0 };
    std::complex<double> L_theta = { 0,0 };
    std::complex<double> L_phi = { 0,0 };

    vtkSmartPointer<vtkCellArray> cellArray = sourceGridwithCellsData->GetCells();
    vtkSmartPointer<vtkDataArray> areas = sourceGridwithCellsData->GetCellData()->GetArray("Areas");

    int nCELLS = cellArray->GetNumberOfCells();

    for (vtkIdType el_i = 0; el_i < nCELLS; el_i++)
    {
        vtkSmartPointer<vtkCell>  c = sourceGridwithCellsData->GetCell(el_i);
        
        double area = areas->GetTuple1(el_i);

        //compute values of N and L
        size_t numNodes = c->GetNumberOfPoints();
        vtkSmartPointer<vtkPoints> pts = c->GetPoints();

        //get cell centroid coordinates
        double x = 0, y = 0, z = 0;
        for (vtkIdType nd_i = 0; nd_i < c->GetNumberOfPoints(); nd_i++)
        {
            x += pts->GetPoint(nd_i)[0] / c->GetNumberOfPoints();
            y += pts->GetPoint(nd_i)[1] / c->GetNumberOfPoints();
            z += pts->GetPoint(nd_i)[2] / c->GetNumberOfPoints();
        }


        double r_dot_r = x * coords_s[0] + y * coords_s[1] + z * coords_s[2];
        std::complex<double> exp_ikrr = exp(i_ * k0 * r_dot_r);

        N_theta += (electricCurrent[el_i][0] * cos(theta) * cos(phi) + electricCurrent[el_i][1] * cos(theta) * sin(phi)
            - electricCurrent[el_i][2] * sin(theta)) * exp_ikrr * area;

        N_phi += ( -electricCurrent[el_i][0] * sin(phi) + electricCurrent[el_i][1] * cos(phi)) * exp_ikrr * area;

        L_theta += (magneticCurrent[el_i][0] * cos(theta) * cos(phi) + magneticCurrent[el_i][1] * cos(theta) * sin(phi)
            - magneticCurrent[el_i][2] * sin(theta)) * exp_ikrr * area;

        L_phi += ( -magneticCurrent[el_i][0] * sin(phi) + magneticCurrent[el_i][1] * cos(phi)) * exp_ikrr * area;

    }


    std::complex<double> E_theta;
    std::complex<double> E_phi;
    double k0_over_4pir = k0 / (4 * pi_ * radius);
    std::complex<double> exp_ikr = exp(-i_ * k0 * radius);

    E_theta = -i_ * k0_over_4pir * (L_phi + Z0 * N_theta) * exp_ikr;
    E_phi = i_ * k0_over_4pir * (L_theta - Z0 * N_phi) * exp_ikr;

    double farField = sqrt ( std::abs(E_theta * conj(E_theta)) + std::abs(E_phi * std::conj(E_phi)) );
    double radIntensity = radius * radius / (2 * Z0) * (std::abs(E_theta * conj(E_theta)) + std::abs(E_phi * std::conj(E_phi)));
    double directivity = radius * radius / (2 * Z0) * (std::abs(E_theta * conj(E_theta)) + std::abs(E_phi * std::conj(E_phi)));
    std::vector<double> returnValues;
    returnValues.push_back(farField);
    returnValues.push_back(radIntensity);
    returnValues.push_back(directivity);
    return returnValues;
}


void nearToFarField::setWaveNumber(double k)
{
    k0 = k;
}

void nearToFarField::setFrequency(double f)
{
    f0 = f;
}

void nearToFarField::setPath(std::string p)
{
    path = p;
}

void nearToFarField::setSourceGrid(vtkSmartPointer<vtkUnstructuredGrid> sg)
{
    sourceGrid = sg;
}

void nearToFarField::addEnclosingGrid(vtkSmartPointer<vtkUnstructuredGrid> grid)
{
    enclosingGrids.push_back(grid);
}

double nearToFarField::getRadiatedPower()
{
    return radiated_power;
}

std::vector<double> nearToFarField::getBoundingBox(std::vector<vtkSmartPointer<vtkUnstructuredGrid>> gridList)
{
    double xmin = +1e200, xmax = -1e200, ymin = +1e200, ymax = -1e200, zmin = +1e200, zmax = -1e200;
    if (gridList.empty())
        return std::vector<double>({ 0, 0, 0, 0, 0, 0 });
    for (auto grid : gridList)
    {
        vtkSmartPointer<vtkPoints> pts = grid->GetPoints();
        pts->ComputeBounds();
        auto bounds = pts->GetBounds();
        if (xmin > bounds[0])
            xmin = bounds[0];
        if (xmax < bounds[1])
            xmax = bounds[1];
        if (ymin > bounds[2])
            ymin = bounds[2];
        if (ymax < bounds[3])
            ymax = bounds[3];
        if (zmin > bounds[4])
            zmin = bounds[4];
        if (zmax < bounds[5])
            zmax = bounds[5];
    }

    return std::vector<double>({ xmin, xmax, ymin, ymax, zmin, zmax });
}


vtkSmartPointer<vtkDataSet> nearToFarField::getNearFieldData(std::vector<vtkSmartPointer<vtkUnstructuredGrid>> gridList)
{
    std::vector<double> boundBox = getBoundingBox(gridList);

    int Nmax = 30;
    double minx = boundBox[0], miny = boundBox[2], minz = boundBox[4];
    double maxx = boundBox[1], maxy = boundBox[3], maxz = boundBox[5];
    double delta_x = maxx - minx, delta_y = maxy - miny, delta_z = maxz - minz;
    double delta_max = std::max(delta_x, std::max(delta_y, delta_z));

    int Nx = std::max(int(ceil(Nmax * delta_x / delta_max)), 3);
    int Ny = std::max(int(ceil(Nmax * delta_y / delta_max)), 3);
    int Nz = std::max(int(ceil(Nmax * delta_z / delta_max)), 3);

    std::vector<double> airB_boundbox = getBoundingBox(gridList);
    double minx_air = boundBox[0], miny_air = boundBox[2], minz_air = boundBox[4];
    double maxx_air = boundBox[1], maxy_air = boundBox[3], maxz_air = boundBox[5];

    vtkSmartPointer<vtkUnstructuredGrid> ugrid = createBox((minx + minx_air)/2, (maxx + maxx_air)/2, 
                                                            (miny + miny_air)/2, (maxy + maxy_air)/2, 
                                                            (minz + minz_air)/2, (maxz + maxz_air)/2, Nx, Ny, Nz);

    // Perform the interpolation
    vtkSmartPointer<vtkProbeFilter> probeFilter =
        vtkSmartPointer<vtkProbeFilter>::New();
    probeFilter->SetSourceData(sourceGrid);
    probeFilter->SetInputData(ugrid);
    probeFilter->Update();

    return probeFilter->GetOutput();
}