#include "circulation.h"
#include <vtkLine.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 <vtkAppendFilter.h>
#include <vtkTransform.h>
#include <vtkTransformFilter.h>
#include <vtkTransformPolyDataFilter.h>

#include <complex>
#include <math.h>


circulation::circulation()
{
    this->sourceGrids =
        vtkSmartPointer<vtkMultiBlockDataSet>::New();
}


std::complex<double> circulation::execute()
{
    std::complex<double> circValue = 0;

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



    vtkSmartPointer<vtkGeometryFilter> geometryFilter1 =
        vtkSmartPointer<vtkGeometryFilter>::New();

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


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

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


    for (vtkIdType i = 0; i < circulationData->GetNumberOfCells(); i++)
    {
        vtkSmartPointer<vtkCell>  c = circulationData->GetCell(i);
        std::vector<std::complex<double>> h_field = std::vector<std::complex<double>>({
            std::complex<double>(h_re->GetTuple(i)[0], h_im->GetTuple(i)[0]),
            std::complex<double>(h_re->GetTuple(i)[1], h_im->GetTuple(i)[1]),
            std::complex<double>(h_re->GetTuple(i)[2], h_im->GetTuple(i)[2]) });

        vtkSmartPointer<vtkPoints> pts = c->GetPoints();
        double pt0[3], pt1[3];
        pts->GetPoint(0, pt0);
        pts->GetPoint(1, pt1);
        circValue += (pt1[0] - pt0[0]) * h_field[0] + (pt1[1] - pt0[1]) * h_field[1] + (pt1[2] - pt0[2]) * h_field[2];
    }
    return circValue;
}

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

void circulation::addSourceGrid(vtkSmartPointer<vtkUnstructuredGrid> sg)
{
    int i = sourceGrids->GetNumberOfBlocks();
    sourceGrids->SetBlock(i, sg);
}

void circulation::setCirculationGrid(vtkSmartPointer<vtkUnstructuredGrid> cg)
{
    circulationGrid = cg;
}

void circulation::setDirectionVector(double x, double y, double z)
{
    vec_x = x;
    vec_y = y;
    vec_z = z;
}



vtkSmartPointer<vtkDataSet> circulation::getCirculationPathData()
{
    std::vector<double> boundBox = getBoundingBox(circulationGrid);

    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 x_0 = (maxx + minx) / 2;
    double y_0 = (maxy + miny) / 2;
    double z_0 = (maxz + minz) / 2;

    std::vector<double> lengthVec = std::vector<double>({ delta_x * (1 - this->vec_x), delta_y * (1 - this->vec_y), delta_z * (1 - this->vec_z) });

    //results are sensitive to the chosen radius!!!!
    double radius = 1.5 * std::sqrt(pow(lengthVec[0], 2) + pow(lengthVec[1], 2) + pow(lengthVec[2], 2));

    vtkSmartPointer<vtkUnstructuredGrid> ugrid = createCirculationPath(x_0, y_0, z_0, radius, this->vec_x, this->vec_y, this->vec_z);

    //create single unstructured grid from multiblock
    vtkSmartPointer<vtkAppendFilter> append =
        vtkSmartPointer<vtkAppendFilter>::New();
    for (int i = 0; i < sourceGrids->GetNumberOfBlocks(); i++)
    {
        append->AddInputData(sourceGrids->GetBlock(i));
    }
    append->Update();

    // write the result
    vtkSmartPointer<vtkUnstructuredGrid> combinedGrid = append->GetOutput();

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

    vtkSmartPointer<vtkDataSet> newGrid =
        vtkSmartPointer<vtkDataSet>(probeFilter->GetOutput());


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

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

    return sourceGridwithCellsData;
}


vtkSmartPointer<vtkUnstructuredGrid> circulation::createCirculationPath(double x_0, double y_0, double z_0, double radius, double v_x, double v_y, double v_z)
{
    vtkSmartPointer<vtkPoints> points =
        vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkUnstructuredGrid> ugrid =
        vtkSmartPointer<vtkUnstructuredGrid>::New();
    int N = 50;
    const double PI = std::acos(-1);
    const double da = 2 * PI / N;

    //calculatate rotation axis and angle
    // axis calculated by vector product of (0,0,1) and (v_x,v_z,v_z)
    // angle - by dot product divided by length of same two vectors
    double rotAxis[3];
    double rotAngle;
    rotAxis[0] = -v_y;
    rotAxis[1] = v_x;
    if ((v_x == 0.0) && (v_y == 0.0))
        rotAxis[2] = 1;
    else
        rotAxis[2] = 0;
    rotAngle = std::acos(v_z / std::sqrt(std::pow(v_x, 2) + std::pow(v_y, 2) + std::pow(v_z, 2)));

    vtkSmartPointer<vtkTransform> transform =
        vtkSmartPointer<vtkTransform>::New();
    transform->RotateWXYZ(rotAngle * 180 / PI, rotAxis);

    // transform x_0, y_0, z_0
    double in[3] = { x_0, y_0, z_0 };
    double newZero[3];
    vtkSmartPointer<vtkTransform> zeroPointTransform =
        vtkSmartPointer<vtkTransform>::New();
    zeroPointTransform->RotateWXYZ(-rotAngle * 180 / PI, rotAxis);

    zeroPointTransform->TransformPoint(in, newZero);

    transform->Translate(newZero);
    
    for (int i = 0; i <= N; i++)
    {
        double pt[3];

        pt[0] = radius * cos(i * da) ;
        pt[1] = radius * sin(i * da);
        pt[2] = 0.0 ;

        points->InsertNextPoint(pt);
    }

    for (int i = 0; i < N; i++)
    {
        vtkSmartPointer<vtkLine> line =
            vtkSmartPointer<vtkLine>::New();
        line->GetPointIds()->SetId(0, i);
        line->GetPointIds()->SetId(1, i+1);
        ugrid->InsertNextCell(line->GetCellType(), line->GetPointIds());
    }

    ugrid->SetPoints(points);

    vtkSmartPointer<vtkTransformFilter > transformFilter =
        vtkSmartPointer<vtkTransformFilter >::New();
    transformFilter->SetTransform(transform);
    transformFilter->SetInputData(ugrid);
    transformFilter->Update();

    return transformFilter->GetUnstructuredGridOutput();
}



std::vector<double> circulation::getBoundingBox(vtkSmartPointer<vtkUnstructuredGrid> grid)
{
    double xmin = +1e200, xmax = -1e200, ymin = +1e200, ymax = -1e200, zmin = +1e200, zmax = -1e200;

    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 });
}