// ------------------------------------------------------------
// main.cc
//
// ------------------------------------------------------------
#define _USE_MATH_DEFINES
#include <cmath>
#include <complex>
#include <iostream>
#include <string>
#include <vector>
#include "source/ugblock2D.h"
#include "ugblock.h"
#include "source/ugblock2D3D.h"


using std::complex;
using namespace ::_COLSAMM_;

#define L2NORM(vector) sqrt(product(vector,vector).real())
#define RESTART 0


double isZero(double x)
{
    return fabs(x) < 1e-10?1.0:0.0;
}

complex<double> I(0.0,1.0);

std::complex<double> expi(double x) {
   return exp(I*x);
}

std::complex<double> sinExp(double x) {
   return sin(x);
}

std::complex<double> cosExp(double x) {
   return cos(x);
}

double realPart(std::complex<double> x) {
   return x.real();
}

double imPart(std::complex<double> x) {
   return x.imag();
}

double complexAngle(std::complex<double> x) {
    return (std::arg(x)+2.0*M_PI);
 }

double radiusLoch;
std::complex<double> loch(double x, double y) {
   if(sqrt(x*x+y*y) < radiusLoch) return 1.0;
   return 0.0;
}

double radiusGauss;
double curvature;
double distanceFromWaist;
double lambda;
double rayleighrange;
std::complex<double> gauss(double x, double y) {
    double k = 2.0 * M_PI / lambda;
    double waist = radiusGauss * sqrt(1+pow(distanceFromWaist / rayleighrange,2));
    if (curvature == 0.0 || std::isnan(curvature))
    {
        return exp(-(x*x+y*y)/(waist*waist));
    }
    else
    {
        return exp(-(x*x+y*y)/(waist*waist))*expi(-1.0 * k * (x*x+y*y) / (2.0 * curvature));
    }
}
std::complex<double> spalt(double x, double y) {
   if(fabs(x) < radiusLoch) return 1.0;
   return 0.0;
}
double myRealSqrt(std::complex<double> z) {
    return sqrt(z.real());
}

std::complex<double> myLimitToOne(std::complex<double> in)
{
    if (in.real() < 1.0)
    {
        return std::complex<double>(1,0);
    }
    return in;
}

void CalcFresnel(std::complex<double> (*aperture) (double x, double y),
                 double distance,
                 double wavenumber,
                 VariableFFT& varIn,
                 VariableFFT& varFar) {
    Blockgrid2D* blockGridRec = varFar.Give_blockgrid();

    X_coordinate2d X(*blockGridRec);
    Y_coordinate2d Y(*blockGridRec);

    Function2d2<std::complex<double>,double> Aperture(aperture);
    Function2d1<std::complex<double>,double> Expi(expi);

    varFar = Aperture(X,Y) * Expi(wavenumber * (X*X+Y*Y) / (2.0 * distance));
    varIn  = varFar;
    varFar.FFT();
}


void spectrumPlaneWave(std::complex<double> (*aperture) (double x, double y),
                       double distance,
                       double wavenumber,
                       VariableFFT& varIn,
                       VariableFFT& varFar,
                       X_coordinate2d& X,
                       Y_coordinate2d& Y,
                       VariableFFT& Kx,
                       VariableFFT& Ky,
                       VariableFFT& temp) {
          Blockgrid2D* blockGridRec = varFar.Give_blockgrid();

//          X_coordinate2d X(*blockGridRec);
//          Y_coordinate2d Y(*blockGridRec);
//          VariableFFT Kx(varFar);
//          VariableFFT Ky(varFar);

//          double dx = varIn.getHx();
//          double dy = varIn.getHy();

//          double ukx = 2.0 * M_PI / varIn.getSizeX();
//          double uky = 2.0 * M_PI / varIn.getSizeY();

//          Kx = X / dx * ukx ;
//          Ky = Y / dy * uky ;

          Function2d2<std::complex<double>,double> Aperture(aperture);
          Function2d1<double, std::complex<double> > absolute(ABS);
          Function2d1<double, std::complex<double> > myRealFunc(myReal);
          Function2d1<std::complex<double>,double> Expi(expi);
          Function2d1<double, std::complex<double> > Sqrt(myRealSqrt);


          //VariableFFT temp(varFar);
          temp = Sqrt(wavenumber*wavenumber - Kx*Kx - Ky*Ky) * (distance) ;
          std::ofstream DATEIA;

          DATEIA.open("varWaveK.vtk");
          temp.Print_VTK(DATEIA);
          DATEIA.close();

          //varIn = varIn * Expi(myRealFunc(temp));

          varFar = Aperture(X,Y);
          varIn  = varFar;
          varFar.FFT();
          varFar = varFar * Expi( myRealFunc(temp));
          varFar.inversFFT();
      }


void virtualLensCorrection(double dz, double wavenumberAverage,
                           VariableFFT& varIn,
                           VariableFFT& varWavenumber,
                           VariableFFT& temp)
{





    Function2d1<std::complex<double>,double> Expi(expi);
    Function2d1<double, std::complex<double> > absolute(ABS);
    Function2d1<double, std::complex<double> > myRealFunc(myReal);

    temp = myRealFunc(varWavenumber - wavenumberAverage ) * dz;
//    std::ofstream DATEIC;
//    DATEIC.open("varTEMPAVERAGE.vtk");
//    temp.Print_VTK(DATEIC);
//    DATEIC.close();



//    temp = myRealFunc(varWavenumber - minimumWavenmumber ) * dz;
//    DATEIC.open("varTEMPMINIMA.vtk");
//    temp.Print_VTK(DATEIC);
//    DATEIC.close();
    //sampling
    double samplingTest =Maximum(myRealFunc(temp)) ;
    double samplingMax = M_PI / samplingTest;
    if (samplingMax < 1.0)
    {
        std::cout << "undersampling \n";
        std::cout << "dz = " << dz << "\n";
        std::cout << "dzMax = " << samplingMax * dz << std::endl;
    }

//    std::ofstream DATEIB;

//    DATEIB.open("varTEMPTEST.vtk");
//    temp.Print_VTK(DATEIB);
//    DATEIB.close();



    varIn = varIn * Expi(1.0*myRealFunc(temp));

//    std::ofstream DATEIA;

//    DATEIA.open("varInAfterCorrection.vtk");
//    varIn.Print_VTK(DATEIA);
//    DATEIA.close();


}

void powerTest(VariableFFT& varIn)
{
    Function2d1<double, std::complex<double> > absolute(ABS);
    double power = product(absolute(varIn),absolute(varIn));
    std::cout << "power  = " << power << "\n";

}
void spectrumPlaneWavePropagation(double distance,
                       double wavenumber,
                       VariableFFT& varIn,
                       VariableFFT& varFar,
                       VariableFFT& Kx,
                       VariableFFT& Ky,
                       VariableFFT& temp ) {
         // Blockgrid2D* blockGridRec = varFar.Give_blockgrid();

//          X_coordinate2d X(*blockGridRec);
//          Y_coordinate2d Y(*blockGridRec);
//          VariableFFT Kx(varFar);
//          VariableFFT Ky(varFar);

//          double dx = varIn.getHx();
//          double dy = varIn.getHy();

//          double ukx = 2.0 * M_PI / varIn.getSizeX();
//          double uky = 2.0 * M_PI / varIn.getSizeY();

//          Kx = X / dx * ukx ;
//          Ky = Y / dy * uky ;


                  Function2d1<double, std::complex<double> > MYREAL(myReal);

          double maxKx = Maximum(MYREAL(Kx));
          double maxKy = Maximum(MYREAL(Ky));

          Function2d1<std::complex<double>,double> Expi(expi);
          Function2d1<double, std::complex<double> > Sqrt(myRealSqrt);
          Function2d1<double, std::complex<double> > absolute(ABS);

         // VariableFFT temp(varIn);
          temp = wavenumber*wavenumber - Kx*Kx - Ky*Ky;
          temp = Sqrt(temp) * (distance);

//          std::ofstream DATEIA;
//          DATEIA.open("varTempA.vtk");
//          temp.Print_VTK(DATEIA);
//          DATEIA.close();


          if ((wavenumber*wavenumber - maxKx*maxKx - maxKy*maxKy) < 0)
          {
              std::cout << "warning : negative sqrt() ! decrease step size or increase resolution" << std::endl;
          }

         // varIn.FFTShift();

          varIn.FFT();
          //varIn.FFTShift();


//          DATEIB.open("varIn.FFTShift-FFT-FFTShift.vtk");
//          varIn.Print_VTK(DATEIB);
//          DATEIB.close();

//          std::ofstream DATEIC;
//          DATEIC.open("varTempWithOutExp.vtk");
//          temp.Print_VTK(DATEIC);
//          DATEIC.close();
          varFar = varIn;
          temp = Expi( 1.0*MYREAL(temp));
          varIn = varIn * temp;




          //varIn.FFTShift();
          varIn.inversFFT();
          //varIn.FFTShift();

//          std::ofstream DATEID;
//          DATEID.open("varTempD.vtk");
//          varIn.Print_VTK(DATEID);
//          DATEID.close();
//varFar = varFar;
         // varFar = varIn;
      }
void applyLens(VariableFFT& varIn,VariableFFT& varOPL, double wavenumber = 1)
{

        Function2d1<std::complex<double>,double> Expi(expi);
        Function2d1<double, std::complex<double> > absolute(ABS);
        Function2d1<double, std::complex<double> > myRealPart(realPart);



        varIn = varIn * Expi(-1.0*myRealPart(varOPL) * wavenumber);



}

double C_4;
double focalLength;

void estimateBeamQuality(VariableFFT& varIn,
                         VariableFFT& Kx,
                         VariableFFT& Ky,
                         X_coordinate2d& X,
                         Y_coordinate2d& Y,
                         VariableFFT& temp,
                         double lambda)
{
    std::cout << "Beamquality not estimated! " << std::endl;
    return;

    Function2d1<double, std::complex<double> > absolute(ABS);
    Function2d1<double, std::complex<double> > myRealFunc(myReal);
    Function2d1<double, std::complex<double> > winkel(complexAngle);
    Local_stiffness_matrix2D<double> dx(*(varIn.Give_blockgrid()));
    dx.Calculate(d_dx(v_())*w_());
    Local_stiffness_matrix2D<double> dy(*(varIn.Give_blockgrid()));
    dy.Calculate(d_dy(v_())*w_());

    CalcContinuousArg cont(*(varIn.Give_blockgrid()));
    int gaussIterations = 100;
    Variable2D<double> f(*(varIn.Give_blockgrid()));
    Variable2D<double> TEMP(*(varIn.Give_blockgrid()));
     Variable2D<double>  ux(*(varIn.Give_blockgrid()));
     Variable2D<double>  uy(*(varIn.Give_blockgrid()));
     Variable2D<double>  px(*(varIn.Give_blockgrid()));
     Variable2D<double>  py(*(varIn.Give_blockgrid()));
     Variable2D<double>  unity(*(varIn.Give_blockgrid()));
     Variable2D<double>  Intensity(*(varIn.Give_blockgrid()));
     Variable2D<double>  IntensityFFT(*(varIn.Give_blockgrid()));
     Variable2D<double>  phase(*(varIn.Give_blockgrid()));
     Variable2D<double>  REALPART(*(varIn.Give_blockgrid()));

//     X_coordinate2d X(*(varIn.Give_blockgrid()));
//     Y_coordinate2d Y(*(varIn.Give_blockgrid()));
//    VariableFFT temp(varIn);


    unity = 1.0;
    Intensity = absolute(varIn) * absolute(varIn);

    temp = varIn;
    temp.FFT();
    IntensityFFT = absolute(temp) * absolute(temp);
    std::ofstream DATEIX;



//    phase = winkel(varIn) - 2.0*M_PI;
//    DATEIX.open("diff_phase_noncontinuous.vtk");
//    phase.Print_VTK(DATEIX);
//    DATEIX.close();

    cont.calcArg(phase,varIn);
    phase = phase * (-1.0);

//    DATEIX.open("diff_phase_from_calcArg.vtk");
//    phase.Print_VTK(DATEIX);
//    DATEIX.close();
//    DATEIX.open("diff_field.vtk");
//    varIn.Print_VTK(DATEIX);
//    DATEIX.close();


    double waveNumber = 2.0 * M_PI / lambda;
    phase = phase / waveNumber;
    //phase = phase;

     Local_stiffness_matrix2D<double> helm(*(varIn.Give_blockgrid())); // also in FEA_Def and solverFEA
    helm.Calculate(v_()*w_());


     (ux) = 0;
    Function2d1<double, std::complex<double> > myRealPart(realPart);
    REALPART = (myRealPart(varIn));
     f = (dy)(Intensity);

     for(int i=0;i<gaussIterations;++i)
     {
         TEMP = ux;
         (ux) = (ux) - ((helm)(ux) -f) / (helm).diag();
         if (L_infty(TEMP-ux) < 1e-13)
             i = gaussIterations;
     }
     f = (dx)(Intensity);
     for(int i=0;i<gaussIterations;++i)
     {
         TEMP = uy;
         (uy) = (uy) - ((helm)(uy) -f) / (helm).diag();
         if (L_infty(TEMP-uy) < 1e-13)
             i = gaussIterations;
     }
     f = (dx)(phase);
     for(int i=0;i<gaussIterations;++i)
     {
         TEMP = px;
         (px) = (px) - ((helm)(px) -f) / (helm).diag();
         if (L_infty(TEMP-px) < 1e-13)
             i = gaussIterations;
     }
     f = (dy)(phase);
     for(int i=0;i<gaussIterations;++i)
     {
         TEMP = py;
         (py) = (py) - ((helm)(py) -f) / (helm).diag();
         if (L_infty(TEMP-py) < 1e-13)
             i = gaussIterations;
     }


//         DATEIX.open("diff_phase_continuous.vtk");
//         px.Print_VTK(DATEIX);
//         DATEIX.close();
//     px = 2.0*X / focalLength / 2.0 + (X*X*X * 4.0 + 4.0 * X*Y*Y) * C_4;
//     //px = px * waveNumber;
//     DATEIX.open("diff_phase_continuous_analytic.vtk");
//     px.Print_VTK(DATEIX);
//     DATEIX.close();
//     py = 2.0*Y / focalLength / 2.0 + (Y*Y*Y * 4.0 + 4.0 * X*X*Y) * C_4;



//     VariableFFT Kx(varIn);
//     VariableFFT Ky(varIn);

//     double delx = varIn.getHx();
//     double dely = varIn.getHy();

//     double ukx = 2.0 * M_PI / varIn.getSizeX();
//     double uky = 2.0 * M_PI / varIn.getSizeY();

//     Kx = (X-0.5*delx) / delx * ukx ;
//     Ky = (Y-0.5*dely) / dely * uky ;
//     DATEIX.open("diff_Ky.vtk");
//     Ky.Print_VTK(DATEIX);
//     DATEIX.close();
//     DATEIX.open("diff_I_FFT.vtk");
//     IntensityFFT.Print_VTK(DATEIX);
//     DATEIX.close();


     double power = product(Intensity,unity);
     double powerFFT = product(IntensityFFT,unity);

 //    std::cout << "ratio of dx / ux " << delx / ukx << std::endl;
     f = myRealFunc(IntensityFFT * Kx * Kx);
     double phiX_FFT = product(f,unity) / powerFFT / waveNumber/ waveNumber;
     f = myRealFunc(IntensityFFT * Ky * Ky);
     double phiY_FFT = product(f,unity) / powerFFT / waveNumber/ waveNumber;
     double medianX = product(Intensity * X * X,unity) / power;
     double medianX4 = product(Intensity * X * X * X * X,unity) / power;
     double medianX6 = product(Intensity * X * X * X * X* X * X,unity) / power;
     double betaX = sqrt((medianX * medianX6 - medianX4 * medianX4)/(medianX4 * medianX4));
     double C4abb = 16.0 * M_PI / lambda * 0.816 * C_4 * medianX4 ;
     C4abb = 16.0 * M_PI / lambda * C_4 * medianX4 ;
     std::cout << "aberration due to C4 " << C4abb << std::endl;
     double medianY = product(Intensity * Y * Y,unity) / power;
     std::cout << "beamwaistX = " << 2.0 * sqrt(medianX) <<std::endl;
     std::cout << "beamwaistY = " << 2.0 * sqrt(medianY) <<std::endl;

     f = 1.0 / 4.0 / waveNumber / waveNumber / power * ux * ux / Intensity;
     f = f + 1.0 / power * Intensity * px * px;
     double phiX = product(f,unity);
     f = 1.0 / 4.0 / waveNumber / waveNumber / power * uy * uy / Intensity;
     f = f + 1.0 / power * Intensity * py * py;
     double phiY = product(f,unity);
     f = 1.0 / power * Intensity * X * px;
     double medianXphiX = product(f,unity);
     f = 1.0 / power * Intensity * Y * py;
     double medianYphiY = product(f,unity);
     double M2X = 2.0 * waveNumber * sqrt(fabs(medianX * phiX_FFT - medianXphiX * medianXphiX));
     double M2Y = 2.0 * waveNumber * sqrt(fabs(medianY * phiY_FFT - medianYphiY * medianYphiY));
     double M2X_direct = 2.0 * waveNumber * sqrt(fabs(medianX * phiX - medianXphiX * medianXphiX));
     double M2Y_direct = 2.0 * waveNumber * sqrt(fabs(medianY * phiY - medianYphiY * medianYphiY));

     double REffX = medianX / medianXphiX;
     double REffY = medianY / medianYphiY;

     f = medianX / power * ux * ux / Intensity;
     double M2diffX = sqrt(product(f,unity));
     f = Intensity * (px - X / REffX)*(px - X / REffX);
     double M2abbX = sqrt(product(f,unity) / power) * 2.0 * waveNumber * medianY;
     f = medianY / power * uy * uy / Intensity;
     double M2diffY = sqrt(product(f,unity));
     f = Intensity * (py - Y / REffY)*(py - Y / REffY);
     double M2abbY = sqrt(product(f,unity) / power) * 2.0 * waveNumber * medianY;

//     std::cout << "medianX = " << medianX <<std::endl;
//     std::cout << "medianX4 = " << medianX4 <<std::endl;
//     std::cout << "medianX6 = " << medianX6 <<std::endl;
//     std::cout << "medianY = " << medianY <<std::endl;
//     std::cout << "phiX = " << phiX <<std::endl;
//     std::cout << "phiY = " << phiY <<std::endl;
//     std::cout << "phiX_FFT = " << phiX_FFT <<std::endl;
//     std::cout << "phiY_FFT = " << phiY_FFT <<std::endl;


//     std::cout << "REffX = " << REffX <<std::endl;
//     std::cout << "REffY = " << REffY <<std::endl;
//     std::cout << "medianXphiX = " << medianXphiX <<std::endl;
//     std::cout << "medianYphiY = " << medianYphiY <<std::endl;


     std::cout << "M2X_with_fft " << M2X << std::endl;
     std::cout << "M2Y_with_fft " << M2Y << std::endl;
     std::cout << "M2X_direct " << M2X_direct << std::endl;
     std::cout << "M2Y_direct " << M2Y_direct << std::endl;

//     std::cout << "M2diffX " << M2diffX << std::endl;
//     std::cout << "M2abbX " << M2abbX << std::endl;
//     std::cout << "M2TotalX " << sqrt(M2abbX*M2abbX+M2diffX*M2diffX) << std::endl;
//     std::cout << "M2diffY " << M2diffY << std::endl;
//     std::cout << "M2abbY " << M2abbY << std::endl;
//     std::cout << "M2TotalY " << sqrt(M2abbY*M2abbY+M2diffY*M2diffY) << std::endl;
//     std::cout << "C4 / M2Abb " << C4abb / M2abbY << std::endl;
//     std::cout << "C4 / M2Abb " << C4abb / M2abbX << std::endl;
//     std::cout << "M2TotalY " << sqrt(M2abbY*M2abbY+M2diffY*M2diffY) << std::endl;

//         std::ofstream DATEIR;
//         DATEIR.open("diffDX.vtk");
//         ux.Print_VTK(DATEIR);
//         DATEIR.close();

    //CalcContinuousArg;


}

int main(int argc, char** argv) {
    std::ofstream DATEI;

    int n = 5;

    double R1 = 100;
    double R2 = 100;

    focalLength = 100.0;

    int distanceIncrements = 15;
    double distance = 50 ;      //[mm]
    double dz = distance / double(distanceIncrements);
    lambda   = 1064e-6;  //[mm]
    radiusLoch = 0.2;   //[mm]


    distance = fabs(distanceFromWaist * 2.0);

    distance = 150 ;
    dz = distance / double(distanceIncrements);
    double geometrySize  =  5.0; //[mm]
    geometrySize =  8.0;

    radiusGauss = geometrySize / 4.0;   //[mm]
    rayleighrange = M_PI *radiusGauss*radiusGauss /lambda;
    std::cout << "rayleighrange " << rayleighrange<< std::endl;
    distanceFromWaist = -rayleighrange;
    distanceFromWaist = 0;
    curvature = distanceFromWaist * (1.0 + pow(rayleighrange / distanceFromWaist,2) );

    double radiusGeo = 0.5 * geometrySize;
    Rechteck geo(-radiusGeo, -radiusGeo, radiusGeo, radiusGeo);

    double wavenumber = 2.0 * M_PI / lambda;


    cout << " Test CalcFresnel!!" << endl;

    std::cout << "distance = " << distance << std::endl;
    std::cout << "lambda = " << lambda << std::endl;
    std::cout << "aperture = " << radiusLoch << std::endl;
    std::cout << "geometrySize = " << geometrySize << std::endl;

    std::cout << "Fresnelnumber = " << radiusLoch*radiusLoch / distance / lambda << std::endl;
    std::cout << "Fresnelnumber >> 1 : near field "      << std::endl;
    std::cout << "Fresnelnumber ~  1 : fresnel zone"     << std::endl;
    std::cout << "Fresnelnumber << 1 : fraunhofer zone"  << std::endl;

    double deltaX = geometrySize / pow(2,n-1);
    double deltaXFourierPlane = lambda * distance / pow(2,n-1) / deltaX;
    std::cout << "deltaX initial plane (z = 0) = " << deltaX << std::endl;
    std::cout << "deltaX fourier plane (z = " <<distance<<")  = " << deltaXFourierPlane << std::endl;
    std::cout << "dxInitial / dxFourier   = " << deltaX / deltaXFourierPlane << std::endl;
    std::cout << "curvature " << curvature << std::endl;


    //VTK_Reader reader(QString("/local/er96apow/FAUbox/Promotion/Vectorial_BPM/fibercryst_exmaple/RefractionIndexTherm_last.vtk"));
    //Variable<double> * thermalRefractiveIndex3D = reader.give_first_variable();


    Unstructured_grid *ug = new Cylinder(2,1,20);
    Blockgrid *bg = new Blockgrid(ug,10,10,20);
    Variable<double> * thermalRefractiveIndex3D = new Variable<double>(*bg);


    VariableFFT varSlice(n,n,geo);
    Variable2D<double> varDoubleRefr(*(varSlice.Give_blockgrid()));
    IteratorZDirection zIterator(thermalRefractiveIndex3D->Give_blockgrid());
    zIterator.gotoFront();
    thermalRefractiveIndex3D->interpolateSlizeZ(&varDoubleRefr, &zIterator);

    VariableFFT varE(n,n,geo);
    VariableFFT varIn(varE);
    VariableFFT varIntensity(varE);
    VariableFFT varPhase(varE);
    VariableFFT varTempX(varE);
    VariableFFT varTempY(varE);
    VariableFFT varRefr(varE);
    VariableFFT varPhaseLens(varE);
    VariableFFT temp(varE);
    VariableFFT varWavenumber(varE);
    Blockgrid2D* blockGridRec = varE.Give_blockgrid();
    X_coordinate2d X(*blockGridRec);
    Y_coordinate2d Y(*blockGridRec);
      VariableFFT Kx(varE);
      VariableFFT Ky(varE);


    double dx = varIn.getHx();
    double dy = varIn.getHy();

    double ukx = 2.0 * M_PI / varIn.getSizeX();
    double uky = 2.0 * M_PI / varIn.getSizeY();

    varTempX =(X  - 0.5*dx);
    varTempY =(Y  - 0.5*dy);
    Kx = varTempX / dx * ukx ;
    Ky = varTempY / dy * uky ;

    Function2d2<std::complex<double>,double> Aperture(gauss);
//    Function2d1<std::complex<double>,double> Expi(expi);
//    Function2d1<double, std::complex<double> > Sqrt(myRealSqrt);

    //sampling criterion


    varE = Aperture(X,Y);

    std::ofstream DATEIG;
    DATEIG.open("/local/er96apow/FFT_results/varE___before.vtk");
    varE.Print_VTK(DATEIG);
    DATEIG.close();

    C_4 = 0.00015;
    //C_4 = 0.0;
//    std::ofstream DATEIX;

//    DATEIX.open("varE_bef_fft.vtk");
//    varE.Print_VTK(DATEIX);
//    DATEIX.close();

    //varE.FFT();
//    std::ofstream DATEIR;
//    DATEIR.open("varE_fft.vtk");
//    varE.Print_VTK(DATEIR);
//    DATEIR.close();

//    varE.inversFFT();
//    std::ofstream DATEIH;
//    DATEIH.open("varE_fft_ifft.vtk");
//    varE.Print_VTK(DATEIH);
//    DATEIH.close();


    varIn = varE;

    Function2d1<double, double > wurzel(sqrt);
    //varPhaseLens =(-1.83+1.0) * (X*X+Y*Y)/2.0 * (1.0 / R1 + 1.0 / R2) *wavenumber;

    varPhaseLens = (X*X + Y*Y) / 2.0 / focalLength + C_4 * ( X * X + Y * Y ) * ( X * X + Y * Y );


    //varPhaseLens = (X*X + Y*Y) / 2.0 / focalLength;// + C_4 * ( X * X + Y * Y ) * ( X * X + Y * Y );
    //    *convergenceVariable = n_0 * wurzel(absolut( 1 - alpha * alpha * ( X * X + Y * Y ) + alpha * alpha * 4.0 * ( X * X + Y * Y ) * ( X * X + Y * Y ) )); //with C4 aberration


    //varPhaseLens = wurzel(X*X + Y*Y) / focalLength;


    Function2d1<std::complex<double>, std::complex<double> > limitToOne(myLimitToOne);
    double alpha = 2.0 * M_PI /  distance / 2.0 ;//200.0;
    alpha = alpha / 2.0;

    alpha = 1.0/1000.0 * 2.0 * M_PI ;
    varRefr = wurzel( 1.5 * 1.5 * (1.0 - alpha * alpha * (X*X+Y*Y)));
   // varRefr = wurzel( 1.5 * 1.5 * (1.0 - alpha * alpha * (X*X)));
    //varRefr = wurzel(1.8 * 1.8 * (1 - alpha * (X*X)));

    varRefr = limitToOne(varRefr);
    varWavenumber = 2.0 * M_PI * varRefr / lambda;

    Function2d1<double, std::complex<double> > myRealFunc(myReal);
    std::cout << "Minimum( (myRealFunc(varRefr))) " << Minimum( (myRealFunc(varRefr))) << std::endl;
    std::cout << "Minimum( X) " << Minimum( (myRealFunc(X))) << std::endl;
    std::cout << "Maximum( X) " << Maximum( (myRealFunc(X))) << std::endl;
    std::cout << " Maximum( (myRealFunc(varRefr))" << Maximum( (myRealFunc(varRefr))) << std::endl;

    Function2d1<double, std::complex<double> > absolute(ABS);


    double refrAverage = 0.1 *  Minimum( (myRealFunc(varRefr))) + 0.9 * Maximum( (myRealFunc(varRefr)));
    double wavenumberAverage = 2.0 * M_PI * refrAverage / lambda;

    //double focalLength = 1.0 / refrAverage / alpha / sin(alpha * distance);
    double peroid = 2.0 * M_PI / alpha;
    std::cout << "focal  length approx: " << focalLength<< std::endl;
    std::cout << "peroid length approx: " << peroid<< std::endl;
    std::cout << "z-dir      increment: " << dz<< std::endl;

//    DATEIG.open("varWavenumber.vtk");
//    varWavenumber.Print_VTK(DATEIG);
//    DATEIG.close();
//    DATEIG.open("varvarPhaseLens.vtk");
//    varPhaseLens.Print_VTK(DATEIG);
//    DATEIG.close();
//    DATEIG.open("varE___before.vtk");
//    varIn.Print_VTK(DATEIG);
//    DATEIG.close();

    double dzMax = M_PI / ( sqrt(wavenumberAverage*wavenumberAverage-pow((pow(2,n)/2-1) * ukx,2))
                           -sqrt(wavenumberAverage*wavenumberAverage-pow((pow(2,n)/2) * ukx,2))) ;

    if (dzMax < dz)
    {
        std::cout << "undersampling " << std::endl;
    }
    //virtualLensCorrection(distance,wavenumberAverage,varIn,varWavenumber);
    applyLens(varIn, varPhaseLens,wavenumber);
    estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);


//        spectrumPlaneWavePropagation(100,
//                     wavenumber,
//                     varIn,
//                     varE,Kx,Ky,X,Y,temp);
//        DATEIG.open("varE___FFT.vtk");
//        varE.Print_VTK(DATEIG);
//        DATEIG.close();
//        DATEIG.open("varE___after.vtk");
//        varIn.Print_VTK(DATEIG);
//        DATEIG.close();
    for (int iter = 0 ; iter < distanceIncrements;iter++)
    {
        std::cout << "progress : " << double(iter) / double(distanceIncrements) * 100 << "%" << std::endl;
        std::ofstream DATEIQ;
//        DATEIQ.open("varIn____befCorrection.vtk");
//        varIn.Print_VTK(DATEIQ);
//        DATEIQ.close();

//        virtualLensCorrection(dz,wavenumberAverage,varIn,varWavenumber);
//        DATEIQ.open("varIn____AftCorrection.vtk");
//        varIn.Print_VTK(DATEIQ);
//        DATEIQ.close();
        spectrumPlaneWavePropagation(dz,
                     wavenumber,
                     varIn,
                     varE,Kx,Ky,temp);
//        DATEIQ.open("varIn____AftPropagation.vtk");
//        varIn.Print_VTK(DATEIQ);
//        DATEIQ.close();
//        powerTest(varIn);


        //virtualLensCorrection(dz,wavenumberAverage,varIn,varWavenumber);



        int plotevery = 1;
        if (iter % plotevery == 0)
        {

            estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);

            std::ofstream DATEIA;

            DATEIA.open("/local/er96apow/FFT_results/varEAtPlane_with_aberration"+std::to_string(iter *dz)+".vtk");
            varIntensity = absolute(varIn) * absolute(varIn);
            varIntensity.Print_VTK(DATEIA);
            DATEIA.close();

//            std::ofstream DATEIB;

//            DATEIB.open("varEFFT"+std::to_string(iter / plotevery)+".vtk");
//            varE.Print_VTK(DATEIB);
//            DATEIB.close();
            if (iter == 101)
                iter += 500;
        }
    }



//    VariableFFT varB(varA);
/*
    Blockgrid2D* blockGridRec = varE.Give_blockgrid();

    X_coordinate2d X(*blockGridRec);
    Y_coordinate2d Y(*blockGridRec);
*/






//    CalcFresnel(loch,
//                 distance,
//                 wavenumber,
//                 varIn,
//                 varE);

//    spectrumPlaneWave(loch,
//                 distance,
//                 wavenumber,
//                 varIn,
//                 varE);
    spectrumPlaneWave(gauss,
                 distance,
                 wavenumber,
                 varIn,
                 varE,X,Y,Kx,Ky,temp);

   // varE = varE / L_infty(  varE);


    std::ofstream DATEIA;

    DATEIA.open("varIn.vtk");
    varIn.Print_VTK(DATEIA);
    DATEIA.close();

    DATEIA.open("varE.vtk");
    varE.Print_VTK(DATEIA);
    DATEIA.close();

//    DATEIA.open("varEimag.vtk");
//	varE.Print_VTK(DATEIA,imPart);
//	DATEIA.close();

//    DATEIA.open("varEreal.vtk");
//	varE.Print_VTK(DATEIA,realPart);
//	DATEIA.close();

    cout << " Test CalcFresnel finished!!" << endl;

}