diff --git a/program2D/main.cc b/program2D/main.cc
index bd13fa835060bc3cee1d35414e200c288de8dbbb..ed42c5c70396be83cd452fb409325f4ca9d0faf6 100644
--- a/program2D/main.cc
+++ b/program2D/main.cc
@@ -38,22 +38,57 @@ std::complex<double> cosExp(double x) {
}
double realPart(std::complex<double> x) {
- return x.real();;
+ return x.real();
}
double imPart(std::complex<double> x) {
- return x.imag();;
+ 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)
+ {
+ 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(abs(x) < radiusLoch) return 1.0;
+ 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,
@@ -72,47 +107,442 @@ void CalcFresnel(std::complex<double> (*aperture) (double x, double y),
varFar.FFT();
}
-double myRealSqrt(std::complex<double> z) {
- return sqrt(z.real());
-}
void spectrumPlaneWave(std::complex<double> (*aperture) (double x, double y),
double distance,
double wavenumber,
VariableFFT& varIn,
- VariableFFT& varFar) {
+ 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);
+// X_coordinate2d X(*blockGridRec);
+// Y_coordinate2d Y(*blockGridRec);
+// VariableFFT Kx(varFar);
+// VariableFFT Ky(varFar);
- double dx = varIn.getHx();
- double dy = varIn.getHy();
+// double dx = varIn.getHx();
+// double dy = varIn.getHy();
- double ukx = 2.0 * M_PI / varIn.getSizeX();
- double uky = 2.0 * M_PI / varIn.getSizeY();
+// double ukx = 2.0 * M_PI / varIn.getSizeX();
+// double uky = 2.0 * M_PI / varIn.getSizeY();
- Kx = X / dx * ukx ;
- Ky = Y / dy * uky ;
+// 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( Sqrt(wavenumber*wavenumber - Kx*Kx - Ky*Ky) * (distance));
+ varFar = varFar * Expi( myRealFunc(temp));
varFar.inversFFT();
}
-void estimateBeamQuality(VariableFFT& varIn, double lambda)
+
+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)
{
- CalcContinuousArg;
+ 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;
}
@@ -122,10 +552,33 @@ int main(int argc, char** argv) {
int n = 10;
- double distance = 1 ; //[mm]
- double lambda = 1.0e-3; //[mm]
- radiusLoch = 0.1; //[mm]
- double geometrySize = 0.8; //[mm]
+ 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 wavenumber = 2.0 * M_PI / lambda;
cout << " Test CalcFresnel!!" << endl;
@@ -150,12 +603,195 @@ int main(int argc, char** argv) {
double radiusGeo = 0.5 * geometrySize;
Rechteck geo(-radiusGeo, -radiusGeo, radiusGeo, radiusGeo);
+
+
+
+
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);
+
+ 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;
+// std::ofstream DATEIG;
+// 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);
+
+// DATEIG.open("varE___before.vtk");
+// varIn.Print_VTK(DATEIG);
+// DATEIG.close();
+// 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("varEAtPlane_newshift"+std::to_string(iter / plotevery)+".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;
+ }
+ }
@@ -178,13 +814,18 @@ int main(int argc, char** argv) {
// varIn,
// varE);
- spectrumPlaneWave(loch,
+// spectrumPlaneWave(loch,
+// distance,
+// wavenumber,
+// varIn,
+// varE);
+ spectrumPlaneWave(gauss,
distance,
wavenumber,
varIn,
- varE);
+ varE,X,Y,Kx,Ky,temp);
- varE = varE / L_infty( varE);
+ // varE = varE / L_infty( varE);
std::ofstream DATEIA;