diff --git a/common_source/mathlib/math_lib_common.h b/common_source/mathlib/math_lib_common.h
index 8f088b635ca83f7a04e310ef697de3b076f43ef0..13228f5dd6053542b87d76da8736d1ff6a676d8a 100644
--- a/common_source/mathlib/math_lib_common.h
+++ b/common_source/mathlib/math_lib_common.h
@@ -109,6 +109,11 @@ inline std::complex<double> myConj(std::complex<double> x) {
return conj(x);
}
+inline std::complex<double> IMAGFROMREAL ( double x)
+{
+ return std::complex<double>(x,0);
+}
+
inline double myReal ( std::complex<double> x )
{
return real ( x );
diff --git a/program/source/interpol/interpol.cc b/program/source/interpol/interpol.cc
index 1c4a441db51a44866367654a65b36112b36b4d05..db21a5ddec8f8fffcd87f189919b918c65f9bc57 100644
--- a/program/source/interpol/interpol.cc
+++ b/program/source/interpol/interpol.cc
@@ -123,7 +123,7 @@ Interpolate_on_structured_grid::~Interpolate_on_structured_grid() {
delete[] lambda;
}
-void Interpolate_on_structured_grid::update_Interpolate_on_structured_grid(Blockgrid &blockgrid_)
+void Interpolate_on_structured_grid::update_Interpolate_on_structured_grid(Blockgrid &blockgrid_, bool onlyOnSurfaceZ)
{
int Nx, Ny, Nz;
int typ;
@@ -186,9 +186,15 @@ void Interpolate_on_structured_grid::update_Interpolate_on_structured_grid(Block
Ny = blockgrid->Give_Ny_hexahedron(id_hex);
Nz = blockgrid->Give_Nz_hexahedron(id_hex);
- for(int k=0;k<Nz;++k)
+ int limitForZ = Nz;
+ if (onlyOnSurfaceZ)
+ {limitForZ = 1;}
+
+ for(int k=0;k<limitForZ;++k)
for(int j=0;j<Ny;++j)
- for(int i=0;i<Nx;++i) {
+ for(int i=0;i<Nx;++i)
+ //if (k == 0 || k == (Nz-1) || j == 0 || j == (Ny-1) || i == 0 || i == (Nx-1))
+ {
// corner points of general hex-cell
cWSD = blockgrid->Give_coord_hexahedron(id_hex,i, j, k );
cESD = blockgrid->Give_coord_hexahedron(id_hex,i+1,j ,k );
@@ -3564,6 +3570,10 @@ void Interpolate_direct::find_surface_cell(D3vector v)
setPrevIndex();
return;
}
+ else
+ {
+ std::cout << "check";
+ }
}
@@ -4176,15 +4186,15 @@ int Interpolate_direct::checkBoxSurface(int id_Hex, int i, int j, int k, D3vecto
//std::cout << "new best : "; lamTemp.Print();
if( typ != -1 && new_lam_better(lambda,lamTemp))
{
- if (MIN(lamTemp) < 0.0 || MAX(lamTemp) >1.0 )
- {
- badLambdaFound = true;
- //lamTemp.Print();
- }
- else
- {
- badLambdaFound = false;
- }
+// if (MIN(lamTemp) < 0.0 || MAX(lamTemp) >1.0 )
+// {
+// badLambdaFound = true;
+// //lamTemp.Print();
+// }
+// else
+// {
+// badLambdaFound = false;
+// }
lambda = lamTemp;
typ_tet = typ;
diff --git a/program/source/interpol/interpol.h b/program/source/interpol/interpol.h
index f581b48579140cee3a7569ba56d2a12ae78857e4..266c683d494b03120c186cb9d1e8d85ef1fedc29 100644
--- a/program/source/interpol/interpol.h
+++ b/program/source/interpol/interpol.h
@@ -79,7 +79,7 @@ class Interpolate_on_structured_grid {
template <class DTyp>
void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
- void update_Interpolate_on_structured_grid(Blockgrid& blockgrid_);
+ void update_Interpolate_on_structured_grid(Blockgrid& blockgrid_, bool onlyOnSurfaceZ = false);
void setInterpolateDirect(Interpolate_direct *id_){id = id_;}
int getNx(){return nx;}
diff --git a/program2D/main.cc b/program2D/main.cc
index ac4a890b1e6d6f72f7923d581aaf8a21bbaedff5..658349b8d404dcda20a88d0ef0bfe764037ba02a 100644
--- a/program2D/main.cc
+++ b/program2D/main.cc
@@ -289,6 +289,8 @@ void spectrumPlaneWavePropagation(double distance,
temp = Expi( 1.0*MYREAL(temp));
varIn = varIn * temp;
+ temp = varIn;
+
@@ -321,6 +323,25 @@ void applyLens(VariableFFT& varIn,VariableFFT& varOPL, double wavenumber = 1)
double C_4;
double focalLength;
+void estimateBeamWaist(VariableFFT& varIn,
+ X_coordinate2d& X,
+ Y_coordinate2d& Y)
+{
+
+ Variable2D<double> unity(*(varIn.Give_blockgrid()));
+ Variable2D<double> Intensity(*(varIn.Give_blockgrid()));
+ unity = 1.0;
+ Function2d1<double, std::complex<double> > absolute(ABS);
+ Intensity = absolute(varIn) * absolute(varIn);
+
+
+ double power = product(Intensity,unity);
+ double medianX = product(Intensity * X * X,unity) / power;
+ 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;
+}
+
void estimateBeamQuality(VariableFFT& varIn,
VariableFFT& Kx,
VariableFFT& Ky,
@@ -329,8 +350,8 @@ void estimateBeamQuality(VariableFFT& varIn,
VariableFFT& temp,
double lambda)
{
- std::cout << "Beamquality not estimated! " << std::endl;
- return;
+// std::cout << "Beamquality not estimated! " << std::endl;
+// return;
Function2d1<double, std::complex<double> > absolute(ABS);
Function2d1<double, std::complex<double> > myRealFunc(myReal);
@@ -529,11 +550,10 @@ void estimateBeamQuality(VariableFFT& varIn,
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 << "M2diffX " << M2diffX << std::endl;
+ std::cout << "M2abbX " << M2abbX << 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;
@@ -552,7 +572,7 @@ void estimateBeamQuality(VariableFFT& varIn,
int main(int argc, char** argv) {
std::ofstream DATEI;
- int n = 5;
+ int n = 8;
double R1 = 100;
double R2 = 100;
@@ -571,9 +591,10 @@ int main(int argc, char** argv) {
distance = 150 ;
dz = distance / double(distanceIncrements);
double geometrySize = 5.0; //[mm]
- geometrySize = 8.0;
+ geometrySize = 2.0;
radiusGauss = geometrySize / 4.0; //[mm]
+ radiusGauss = 0.2;
rayleighrange = M_PI *radiusGauss*radiusGauss /lambda;
std::cout << "rayleighrange " << rayleighrange<< std::endl;
distanceFromWaist = -rayleighrange;
@@ -606,22 +627,33 @@ int main(int argc, char** argv) {
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();
+ 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);
- 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()));
+ UGFrom3DSlice slice(thermalRefractiveIndex3D->Give_Ug(),0.0);
+
+ Blockgrid2DFrom3D D2block(&slice,thermalRefractiveIndex3D->Give_blockgrid());
+ Variable2D<std::complex<double> > vardDeltaNComplex(D2block);
+ Variable2D<double > varDeltaN(D2block);
IteratorZDirection zIterator(thermalRefractiveIndex3D->Give_blockgrid());
- zIterator.gotoFront();
- thermalRefractiveIndex3D->interpolateSlizeZ(&varDoubleRefr, &zIterator);
+ zIterator.next();
+ varDeltaN.interpolateSlizeZ(thermalRefractiveIndex3D,0.0);
+
+ vardDeltaNComplex = varDeltaN;
VariableFFT varE(n,n,geo);
+
+ VariableFFT deltaN(varE);
+
+ deltaN.interpolate(vardDeltaNComplex,0.0);
+
VariableFFT varIn(varE);
VariableFFT varIntensity(varE);
VariableFFT varPhase(varE);
@@ -703,6 +735,7 @@ int main(int argc, char** argv) {
double alpha = 2.0 * M_PI / distance / 2.0 ;//200.0;
alpha = alpha / 2.0;
+ double n0 = 1.83;
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)));
@@ -747,67 +780,72 @@ int main(int argc, char** argv) {
std::cout << "undersampling " << std::endl;
}
//virtualLensCorrection(distance,wavenumberAverage,varIn,varWavenumber);
- applyLens(varIn, varPhaseLens,wavenumber);
- estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);
+// 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++)
+
+ for (int iter = 0 ; iter < zIterator.getSize();iter++)
{
- std::cout << "progress : " << double(iter) / double(distanceIncrements) * 100 << "%" << std::endl;
+
+ std::cout << "progress : " << double(iter) / double(zIterator.getSize()) * 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,
+ varDeltaN.interpolateSlizeZ(thermalRefractiveIndex3D,zIterator.get_z());
+ vardDeltaNComplex = varDeltaN;
+ deltaN.interpolate(vardDeltaNComplex,0.0);
+
+
+
+// deltaN = 0.0;
+// n0 = 1.0;
+
+ varWavenumber = (n0 + deltaN) * 2.0 * M_PI / lambda;
+ double wavenumberAverage_ = 0.1 * Minimum( (myRealFunc(varWavenumber))) + 0.9 * Maximum( (myRealFunc(varWavenumber)));
+
+
+ virtualLensCorrection(zIterator.get_hz(),wavenumberAverage_,varIn,varWavenumber, temp);
+ spectrumPlaneWavePropagation( zIterator.get_hz(),
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;
+ int plotevery = 10;
if (iter % plotevery == 0)
{
- estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);
+ //estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);
+ //estimateBeamWaist(varIn,X,Y);
+
+
std::ofstream DATEIA;
- DATEIA.open("/local/er96apow/FFT_results/varEAtPlane_with_aberration"+std::to_string(iter *dz)+".vtk");
+ DATEIA.open("/local/er96apow/FFT_results/fibercryst_lowres_"+std::to_string(iter)+".vtk");
varIntensity = absolute(varIn) * absolute(varIn);
varIntensity.Print_VTK(DATEIA);
DATEIA.close();
-// std::ofstream DATEIB;
+ DATEIA.open("/local/er96apow/FFT_results/fibercryst_fftspace_lowres_"+std::to_string(iter)+".vtk");
+ temp.Print_VTK(DATEIA);
+ DATEIA.close();
+
+// vardDeltaNComplex.interpolateSlizeZ(varIntensity,0.0);
-// DATEIB.open("varEFFT"+std::to_string(iter / plotevery)+".vtk");
-// varE.Print_VTK(DATEIB);
-// DATEIB.close();
- if (iter == 101)
- iter += 500;
+// DATEIA.open("/local/er96apow/FFT_results/interpolated_"+std::to_string(iter)+".vtk");
+// //varIntensity = absolute(varIn) * absolute(varIn);
+// vardDeltaNComplex.Print_VTK(DATEIA);
+// DATEIA.close();
+// if (iter == 101)
+// iter += 500;
}
+
+ zIterator.next();
}
@@ -836,24 +874,24 @@ int main(int argc, char** argv) {
// wavenumber,
// varIn,
// varE);
- spectrumPlaneWave(gauss,
- distance,
- wavenumber,
- varIn,
- varE,X,Y,Kx,Ky,temp);
+// spectrumPlaneWave(gauss,
+// distance,
+// wavenumber,
+// varIn,
+// varE,X,Y,Kx,Ky,temp);
- // varE = varE / L_infty( varE);
+// // varE = varE / L_infty( varE);
- std::ofstream DATEIA;
+// std::ofstream DATEIA;
- DATEIA.open("varIn.vtk");
- varIn.Print_VTK(DATEIA);
- DATEIA.close();
+// DATEIA.open("varIn.vtk");
+// varIn.Print_VTK(DATEIA);
+// DATEIA.close();
- DATEIA.open("varE.vtk");
- varE.Print_VTK(DATEIA);
- DATEIA.close();
+// DATEIA.open("varE.vtk");
+// varE.Print_VTK(DATEIA);
+// DATEIA.close();
// DATEIA.open("varEimag.vtk");
// varE.Print_VTK(DATEIA,imPart);
diff --git a/program2D/source/extemp/interpolateFrom3D.h b/program2D/source/extemp/interpolateFrom3D.h
index 0bddfd09ce0bb90ecc2d3458645ec52adc087b88..dfdcf447f301f617ddf6101dd4739b39acf40214 100644
--- a/program2D/source/extemp/interpolateFrom3D.h
+++ b/program2D/source/extemp/interpolateFrom3D.h
@@ -75,6 +75,7 @@ public:
int get_k_intern() { return k; } ///> nur nummer in einem Block
int get_nummerBlock() { return nummerBlock; }
double get_z() { return z_Wert; }
+ double get_hz() { return hz; }
TypeOfUG giveTypeUg() {return typeUg;}
private: