From f8c65b4e99f82883d4cfaab348d45f8f4df5fabc Mon Sep 17 00:00:00 2001
From: iwia06 <christoph.pflaum@fau.de>
Date: Fri, 30 Jul 2021 12:32:08 +0200
Subject: [PATCH] im wesentlichen bei FFT variablen scale eingefuegt
---
program2D/main.cc | 1143 ++-------------------
program2D/source/extemp/variableFFT.cc | 19 +
program2D/source/extemp/variableFFT.h | 13 +-
program2D/source/grid/blockgrid2D.h | 3 +
program2D/source/interpol/interpolTwoD.cc | 1 +
5 files changed, 103 insertions(+), 1076 deletions(-)
diff --git a/program2D/main.cc b/program2D/main.cc
index 7cf6766..1253275 100644
--- a/program2D/main.cc
+++ b/program2D/main.cc
@@ -9,70 +9,15 @@
#include <string>
#include <vector>
#include <memory>
-#include "ugblock.h"
+//#include "ugblock.h"
#include "source/ugblock2D.h"
-#include "source/ugblock2D3D.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> expComplex(std::complex<double> x) {
- return exp(x.real());
-}
-
-
-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 lochDouble(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) {
//gauss electric fielöd
double k = 2.0 * M_PI / lambda;
@@ -83,1037 +28,95 @@ std::complex<double> gauss(double x, double y) {
}
else
{
- return exp(-(x*x+y*y)/(waist*waist))*expi(-1.0 * k * (x*x+y*y) / (2.0 * curvature));
+// 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 amplification(double dz,
- VariableFFT& pumppower,
- VariableFFT& varIn,
- VariableFFT& varWavenumber,
- VariableFFT& temp)
-{
-
- // varIn : Intensity = |varIn|^2
- // also das E-Feld so skaliert, dass es quadriert direkt die Intensität ergibt.
-
- double c = 3e11;// mm / s
- //c = 3e8;// m / s
- double planck = 6.626e-34;
- double emissionCrosssection = 2.5* 1e-17;//mm²
- double upperLevelLifetime = 230e-6;
- double Ntot = 1.39e+17; //( 1 / mm^3)
- double lambdaPump = 808e-6;
- // lambda = 1064e-6;
- //doping 0.1&
- //Ntot = 1.3e+26; //( 1 / m^3)
-
- //nd:yag 1%
- upperLevelLifetime = 225e-6;
- Ntot = 1.39e+17;
- //pumppower : W / mm³
-
- VariableFFT photonDensity(temp);
- VariableFFT inversion(temp);
- VariableFFT pumpPhotons(temp);
- Function2d1<double, std::complex<double> > absolute(ABS);
- // photonDensity : Intensity / energie_photon , unit: 1 / (s mm^2)
- photonDensity = absolute(varIn) * absolute(varIn)/ (planck * c / lambda);
- // pumpPhotons : pumppower / energie_photon , unit: 1 / (s mm^3)
- pumpPhotons = pumppower / (planck * c / lambdaPump); // N / s mm³
-
- // exp complex: normale exp funktion, die aber den realteil des arguments in den exponenten packe, also exp(real(arg)) und einen complexen wert (imag = 0) zurückgibt
- //noetig, weil multipliziert mit komplexer variable
- Function2d1<std::complex<double>,std::complex<double>> Exp(expComplex);
-
- //inversion : stimmt das so?
- inversion = pumpPhotons / (emissionCrosssection * photonDensity + 1.0 / upperLevelLifetime + pumpPhotons / Ntot);
- //inversion = pumpPhotons / (emissionCrosssection * photonDensity * 2.0 + 1.0 / upperLevelLifetime + pumpPhotons / Ntot);
-
- //temp : argument, welches in den exponenten zur verstärkung kommt
- temp= emissionCrosssection * inversion * dz ;
-
-
- Function2d1<double, std::complex<double> > Sqrt(myRealSqrt);
-
- //verstärkungsschritt. hier muss mit der wurzel der verstärkung multipliziert werden, da varIn die wurzel der Intensität ist
- temp = Sqrt(Exp(temp));
- //temp = Sqrt(1.0 + emissionCrosssection * inversion * dz);
-
- varIn = varIn * temp;
-
-
-}
-
-void virtualLensCorrection(double dz, double wavenumberAverage,
- VariableFFT& varIn,
- VariableFFT& varWavenumber,
- VariableFFT& temp)
-{
- std::cout << "virtualLensCorrection skpped\n";
- return;
-
-
-
-
-
- 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();
+*/
+double fRealPart(std::complex< double > x) {
+ return x.real();
}
-void powerTest(VariableFFT& varIn)
-{
- Function2d1<double, std::complex<double> > absolute(ABS);
- double power = product(absolute(varIn),absolute(varIn));
- std::cout << "power = " << power << "\n";
-
+double fImagPart(std::complex< double > x) {
+ return x.imag();
}
-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;
-
- temp = varIn;
-
-
-
-
- //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);
-
-
+std::complex< double > fexpi(double x) {
+ std::complex< double > I(0.0,1.0);
+ return exp(I*x);
}
-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 powerHelper = product(absolute(Intensity),absolute(unity));
-
- double medianX = product(Intensity * X * X,unity) / powerHelper;
- double medianY = product(Intensity * Y * Y,unity) / powerHelper;
- powerHelper = powerHelper *varIn.getHx() * varIn.getHy();
- std::cout << "power = " << powerHelper <<std::endl;
- std::cout << "beamwaistX = " << 2.0 * sqrt(medianX) <<std::endl;
- std::cout << "beamwaistY = " << 2.0 * sqrt(medianY) <<std::endl;
+double fexp(double x) {
+ return exp(x);
}
-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);
-
- double incREAL = varIn.getHx();
- double incFREQ = 2.0 * M_PI / varIn.getSizeY();
- incFREQ = 0.000248;
-
- incFREQ = varIn.getHx() / varIn.getNx();
- double powerT = product(Intensity,unity);// * incREAL * incREAL;
- double powerFFTT = product(IntensityFFT,unity);// * incFREQ * incFREQ;
- Intensity = Intensity / powerT;
- IntensityFFT = IntensityFFT / powerFFTT;
- 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;
-
- f = myRealFunc(IntensityFFT * X * Kx);
- double medianXphiX_FFT = product(f,unity)/ waveNumber/ waveNumber;
- f = myRealFunc(absolute(varIn) * absolute(temp) * Y * Ky);// fehler hier!
- f = myRealFunc(IntensityFFT * Intensity * Y);// fehler hier!
- double medianYphiY_FFT = product(f,unity)/ 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_FFT = 2.0 * waveNumber * sqrt(fabs(medianX * phiX_FFT - medianXphiX_FFT*medianXphiX_FFT));
- double M2Y_FFT = 2.0 * waveNumber * sqrt(fabs(medianY * phiY_FFT - medianYphiY_FFT*medianYphiY_FFT));
- 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 << "M2X_FFT " << M2X_FFT << std::endl;
- std::cout << "M2Y_FFT " << M2Y_FFT << 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;
-// std::cout << "M2TotalY " << sqrt(M2abbY*M2abbY+M2diffY*M2diffY) << std::endl;
-
-// std::ofstream DATEIR;
-// DATEIR.open("diffDX.vtk");
-// ux.Print_VTK(DATEIR);
-// DATEIR.close();
-
- //CalcContinuousArg;
-
-
+double SG_extern;
+double stretch_extern;
+double rSG(double x, double y) {
+ double r = sqrt(x*x+y*y);
+ r = r / stretch_extern;
+ return pow(r,SG_extern);
}
int main(int argc, char** argv) {
std::ofstream DATEI;
- int n = 8;
-
- 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.1; //[mm]
-
-
- distance = fabs(distanceFromWaist * 2.0);
-
- distance = 150 ;
- dz = distance / double(distanceIncrements);
- double geometrySize = 5.0; //[mm]
- geometrySize = 1.0;
-
- radiusGauss = geometrySize / 4.0; //[mm]
- radiusGauss = 0.2;
- 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;
+ double lambda = 1064e-6; //[mm]
+
+
+ Function2d1< double, std::complex< double > > realPart(fRealPart);
+ Function2d1< double, std::complex< double > > imagPart(fImagPart);
+
+ Function2d1< std::complex< double >, double > Expi(fexpi);
+ Function2d1< double, double > Exp(fexp);
+ Function2d2< double, double > powSG(rSG);
+
+ double radiusGeo = 0.1;
+// double radiusGeo = 1.0;
Rechteck geo(-radiusGeo, -radiusGeo, radiusGeo, radiusGeo);
+ double w = 0.0705;
+ stretch_extern = 1.41421;
+// double w = 0.5;
+ double SG = 15.6;
+// double SG = 2.0;
+ SG_extern = SG;
+ double R = 2.03;
+
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;
-
-
- //wenn daten gelesen werden
-// VTK_Reader reader( QString("C:/Users/rall/FAUbox/Promotion/Vectorial_BPM/fibercryst_exmaple/RefractionIndexTherm_last.vtk"));
-// VTK_Reader readerPumppower(QString("C:/Users/rall/FAUbox/Promotion/Vectorial_BPM/fibercryst_exmaple/abspower.vtk"));
-// Variable<double> * thermalRefractiveIndex3D = reader.give_first_variable();
-// Variable<double> * pumpPowerRaytracing = readerPumppower.give_first_variable();
- //wenn analytisch beschrieben
- Cylinder *cyl = new Cylinder(0.5,0.25,20.0);
- Blockgrid *bg = new Blockgrid(cyl,8,8,128);
- Variable<double> *thermalRefractiveIndex3D = new Variable<double>(*bg);
- Variable<double> *pumpPowerRaytracing = new Variable<double>(*bg);
-
- Function2<double> flatTopPump(lochDouble);
- *pumpPowerRaytracing = *pumpPowerRaytracing * 140.0;
-
- std::cout << "manually setting pumppower:" << std::endl;
-
- //homogeneous
- //*pumpPowerRaytracing = 140.0 / (0.5*0.5*M_PI*20.0);
- X_coordinate X3D(*(pumpPowerRaytracing->Give_blockgrid()));
- Y_coordinate Y3D(*(pumpPowerRaytracing->Give_blockgrid()));
- *pumpPowerRaytracing =flatTopPump(X3D,Y3D)* 140.0 / (radiusLoch*radiusLoch*M_PI*20.0);
-
-
-
-// Unstructured_grid *ug = new Cylinder(2,1,20);
-// Blockgrid *bg = new Blockgrid(ug,10,10,20);
-// Variable<double> * thermalRefractiveIndex3D = new Variable<double>(*bg);
-
-
-
- UGFrom3DSlice slice(thermalRefractiveIndex3D->Give_Ug(),0.0);
- UGFrom3DSlice slicePumppower(pumpPowerRaytracing->Give_Ug(),0.0);
-
- Blockgrid2DFrom3D D2block(&slice,thermalRefractiveIndex3D->Give_blockgrid());
- Blockgrid2DFrom3D D2blockPumppower(&slicePumppower,pumpPowerRaytracing->Give_blockgrid());
- Variable2D<std::complex<double> > vardDeltaNComplex(D2block);
- Variable2D<std::complex<double> > varPumppowerComplex(D2blockPumppower);
-
- Variable2D<double > varDeltaN(D2block);
- Variable2D<double > varPumppower(D2blockPumppower);
-
- IteratorZDirection zIterator(thermalRefractiveIndex3D->Give_blockgrid());
- zIterator.next();
- varDeltaN.interpolateSlizeZ(thermalRefractiveIndex3D,0.0);
-
-
- vardDeltaNComplex = varDeltaN;
- VariableFFT varE(n,n,geo);
-
- VariableFFT deltaN(varE);
- VariableFFT pumppower(varE);
-
- deltaN.interpolate(vardDeltaNComplex,0.0);
-
- VariableFFT varIn(varE);
- VariableFFT unit(varE);
- unit = 1.0;
- 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 varPumpslice(varE);
-
-
- Function2d2<std::complex<double>,double> flatTopSlice(loch);
- varPumpslice = 140.0 * flatTopSlice(X,Y) / (radiusLoch*radiusLoch*M_PI*20.0);
-
- 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 ;
-
-
- Function2d1<double, std::complex<double> > absolute(ABS);
- double incFREQ = (fabs(Maximum(absolute(Kx) * 2.0))) /double(Kx.getSizeX()) ;
-
- 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);
-
- double power = 1.0; // unit : Watt
-
- double amp = sqrt(2.0 * power / M_PI / radiusGauss / radiusGauss); // unit = Watt / mm^2
- varE = varE * amp;
-
-// 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;
-
- double n0 = 1.8230;
- 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);
-
- //varRefr = 1.83+deltaN ;
-
- 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;
-
-
-
- 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);
-
-
- double totalPumpPower = 0;
- Function2d1<double, std::complex<double> > absolute2(ABS);
-
- temp = absolute2(varE)*absolute2(varE);
- double powerHelper = product(absolute2(temp),absolute2(unit));
- powerHelper = powerHelper *varE.getHx() * varE.getHy();
- std::cout << "power before : "<<powerHelper << std::endl;
- estimateBeamWaist(varIn,X,Y);
-// spectrumPlaneWave(gauss,
-// 118.105,
-// //wavenumber,
-// 2.0 * M_PI / lambda,
-// varIn,
-// varE,X,Y,Kx,Ky,temp);
-// spectrumPlaneWavePropagation( 118.105,
-// wavenumber * 2,
-// varIn,
-// varE,Kx,Ky,temp);
-// CalcFresnel(gauss,
-// 118.105,
-// wavenumber,
-// varIn,
-// varE);
-// estimateBeamWaist(varIn,X,Y);
-// return 0;
-
-
- double distanceTotal{0};
- int iterCounter{0};
- for (int iter = 0 ; iter < zIterator.getSize()-1;iter++)
- {
- std::cout << "iter " << iter<< "\n";
-
- std::cout << "progress : " << double(iter) / double(zIterator.getSize()) * 100 << "%" << std::endl;
- std::ofstream DATEIQ;
- varDeltaN.interpolateSlizeZ(thermalRefractiveIndex3D,zIterator.get_z());
-
- varPumppower.interpolateSlizeZ(pumpPowerRaytracing, zIterator.get_z());
- vardDeltaNComplex = varDeltaN;
- varPumppowerComplex = varPumppower;
- deltaN.interpolate(vardDeltaNComplex,0.0);
- pumppower.interpolate(varPumppowerComplex,0.0);
-
- //varPumpslice
- std::cout << "pumpslice direct on fft domain calculated\n";
- pumppower =varPumpslice;
- totalPumpPower += product(absolute2(pumppower),absolute(unit)) * varE.getHx() * varE.getHy() * zIterator.get_hz();
-
-
-// deltaN = 0.0;
-// n0 = 1.0;
-
- std::cout << "refrindex gradeint neglected" << std::endl;
- //varWavenumber = (n0 + deltaN) * 2.0 * M_PI / lambda;
- varWavenumber = (n0) * 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);
-
-
- double fac = 1.0;
- iter--;
- for (int iterFactor = 0 ; iterFactor < 1.0 / fac;iterFactor++ )
- {
- iter++;
- }
- for (int iterDiscret = 0 ; iterDiscret < fac; iterDiscret++)
- {
- iterCounter++;
- distanceTotal += zIterator.get_hz() / fac;
- amplification(zIterator.get_hz() / fac,pumppower,varIn,varWavenumber, temp);
- spectrumPlaneWavePropagation( zIterator.get_hz() / fac,
- wavenumberAverage_,
- varIn,
- varE,Kx,Ky,temp);
- }
-
-
-
-
- int plotevery = 100;
- if (iter % plotevery == 0)
- {
-
- //estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);
- estimateBeamWaist(varIn,X,Y);
-
-
-
- std::ofstream DATEIA;
-
- DATEIA.open("FFT_results/fibercryst_lowres_"+std::to_string(iter)+".vtk");
- varIntensity = absolute(varIn) * absolute(varIn);
- varIntensity.Print_VTK(DATEIA);
- DATEIA.close();
- DATEIA.open("FFT_results/fibercryst_lowres_Field_"+std::to_string(iter)+".vtk");
- varIntensity = absolute(varIn);
- varIntensity.Print_VTK(DATEIA);
- DATEIA.close();
- DATEIA.open("FFT_results/pump_power"+std::to_string(iter)+".vtk");
- pumppower.Print_VTK(DATEIA);
- DATEIA.close();
-
-// 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);
-
-// 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();
- }
- estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);
- temp = absolute2(varIn)*absolute2(varIn);
- powerHelper = product(absolute2(temp),absolute2(unit));
- powerHelper = powerHelper *varE.getHx() * varE.getHy();
- std::cout << "power after : "<<powerHelper << std::endl;
-
- std::cout << "total pump power : " << totalPumpPower << std::endl;
- std::cout << "distanceTotal " << distanceTotal << std::endl;
- std::cout << "iterCounter " << iterCounter << std::endl;
-// 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;
-
+ int depth = 8;
+
+ VariableFFT testU (depth, depth, geo);
+ Blockgrid2D* bg = testU.Give_blockgrid();
+
+ Variable2D<double> RealPart(*bg);
+ Variable2D<double> ImagPart(*bg);
+
+ X_coordinate2d X(*bg);
+ Y_coordinate2d Y(*bg);
+
+ testU = Exp(-powSG(X/w,Y/w));
+// testU = Exp(-powSG(X,Y)/(w*w));
+
+ testU = testU * Expi(-wavenumber*(X*X+Y*Y)/(2.0*R) );
+
+// testU = Y;
+
+ DATEI.open("testU.vtk");
+ testU.Print_VTK(DATEI);
+ DATEI.close();
+
+ RealPart = realPart(testU);
+ ImagPart = imagPart(testU);
+
+ DATEI.open("realU.vtk");
+ RealPart.Print_VTK(DATEI);
+ DATEI.close();
+
+ DATEI.open("imagU.vtk");
+ ImagPart.Print_VTK(DATEI);
+ DATEI.close();
+
}
diff --git a/program2D/source/extemp/variableFFT.cc b/program2D/source/extemp/variableFFT.cc
index 49c380a..a17e592 100644
--- a/program2D/source/extemp/variableFFT.cc
+++ b/program2D/source/extemp/variableFFT.cc
@@ -125,6 +125,25 @@ void VariableFFT::init() {
assignmentBlockgrid = NULL;
}
+void VariableFFT::setDomainStrech(double strechX, double strechY) {
+ Blockgrid2D* blockGridFFT = Give_blockgrid();
+ blockGridFFT->SetDomainStrech(strechX,strechY);
+}
+
+double VariableFFT::getHx() const {
+ return Hx * Give_blockgrid()->getStrechX();
+}
+double VariableFFT::getHy() const {
+ return Hy * Give_blockgrid()->getStrechY();
+}
+double VariableFFT::getNx() { return Nx; }
+double VariableFFT::getNy() { return Ny; }
+double VariableFFT::getSizeX() {
+ return Hx*Nx * Give_blockgrid()->getStrechX();
+}
+double VariableFFT::getSizeY() {
+ return Hy*Ny * Give_blockgrid()->getStrechY();
+}
VariableFFT::~VariableFFT() {
if(ownBlockGrid) blockgrid;
diff --git a/program2D/source/extemp/variableFFT.h b/program2D/source/extemp/variableFFT.h
index cb4dd48..8477929 100644
--- a/program2D/source/extemp/variableFFT.h
+++ b/program2D/source/extemp/variableFFT.h
@@ -50,6 +50,7 @@ public:
VariableFFT(VariableFFT* other);
VariableFFT(const VariableFFT& other);
~VariableFFT();
+ void setDomainStrech(double strechX, double strechY);
void FFT();
void inversFFT();
@@ -67,12 +68,12 @@ public:
std::complex<double> defaultInterpolation);
- double getHx() const { return Hx; }
- double getHy() const { return Hy; }
- double getNx() const { return Nx; }
- double getNy() const { return Ny; }
- double getSizeX() const { return Hx*Nx; }
- double getSizeY() const { return Hy*Ny; }
+ double getHx() const;
+ double getHy() const;
+ double getNx();
+ double getNy();
+ double getSizeX();
+ double getSizeY();
/*
diff --git a/program2D/source/grid/blockgrid2D.h b/program2D/source/grid/blockgrid2D.h
index adaf24e..6249cca 100644
--- a/program2D/source/grid/blockgrid2D.h
+++ b/program2D/source/grid/blockgrid2D.h
@@ -61,6 +61,9 @@ class Blockgrid2D {
* @param strechDomainY streched die Geometrie um Faktor strechDomainY in y-Richtung
*/
void SetDomainStrech(double strechDomainX_,double strechDomainY_) { strechDomainX = strechDomainX_; strechDomainY = strechDomainY_; }
+ double getStrechX() const { return strechDomainX; }
+ double getStrechY() const { return strechDomainY; }
+
inline bool inStrechedDomain(const D2vector& point_coord);
double MinimumStrechedX() { return ug->MinimumX()*strechDomainX; };
double MinimumStrechedY() { return ug->MinimumY()*strechDomainY; };
diff --git a/program2D/source/interpol/interpolTwoD.cc b/program2D/source/interpol/interpolTwoD.cc
index 9b89887..83f4bca 100644
--- a/program2D/source/interpol/interpolTwoD.cc
+++ b/program2D/source/interpol/interpolTwoD.cc
@@ -906,6 +906,7 @@ void VectorOfArrowsCreatorComplex::createVectorList(QVector<QPointF>& startP, QV
double signX = 1.0;
double signY = 1.0;
+//test NOW
if(valueUx.real() < 0.0) signX = -1.0;
if(valueUy.real() < 0.0) signY = -1.0;
--
GitLab