Dear CS10-Gitlab-users, on Thursday, Feb 3 there will be maintenance. That will lead to a downtime of the CS10-Gitlab-service including Subversion and Mattermost chat from 09:30. This might take the whole day since we don't know how long it is going to take. We are sorry for the inconvenience! Best regards, CS10-Admin-Team

Commit f8c65b4e authored by Christoph Pflaum's avatar Christoph Pflaum
Browse files

im wesentlichen bei FFT variablen scale eingefuegt

parent b020896d
......@@ -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));