Skip to content

GitLab

Commit 110859af authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

added exmaple : reads refractive index (deltaN) from 3v var and propagates it along the crystal

parent 22d6b44f
Hide whitespace changes
Inline Side-by-side
Showing with 133 additions and 79 deletions
common_source/mathlib/math_lib_common.h
View file @ 110859af
......@@ -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 );
......
program/source/interpol/interpol.cc
View file @ 110859af
......@@ -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;
......
program/source/interpol/interpol.h
View file @ 110859af
......@@ -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;}
......
program2D/main.cc
View file @ 110859af
......@@ -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);
......
program2D/source/extemp/interpolateFrom3D.h
View file @ 110859af
......@@ -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:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment