diff --git a/program2D/main.cc b/program2D/main.cc
index 82a558c64198d61670d6bc0a130433b89bdfcadc..7cf67666a03457c8e24a80d5ef236ff9545d9781 100644
--- a/program2D/main.cc
+++ b/program2D/main.cc
@@ -8,6 +8,7 @@
#include <iostream>
#include <string>
#include <vector>
+#include <memory>
#include "ugblock.h"
#include "source/ugblock2D.h"
#include "source/ugblock2D3D.h"
@@ -62,12 +63,18 @@ std::complex<double> loch(double x, double y) {
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;
double waist = radiusGauss * sqrt(1+pow(distanceFromWaist / rayleighrange,2));
if (curvature == 0.0 || std::isnan(curvature))
@@ -171,6 +178,10 @@ void amplification(double dz,
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;
@@ -182,43 +193,38 @@ void amplification(double dz,
//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), dz richtig hier?
+ // 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³
-// std::cout << "pumppower "<< L_infty(pumppower) << std::endl;
-// std::cout << "inversion "<< L_infty(temp) << std::endl;
-// std::cout << "photonDensity "<< L_infty(photonDensity) << std::endl;
// 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 + 1.0 / upperLevelLifetime + pumpPhotons / Ntot);
+ //inversion = pumpPhotons / (emissionCrosssection * photonDensity * 2.0 + 1.0 / upperLevelLifetime + pumpPhotons / Ntot);
-// std::cout << "pumppower "<< L_infty(pumppower) << std::endl;
-// std::cout << "photonDensity "<< L_infty(photonDensity) << std::endl;
-// std::cout << "inversion "<< L_infty(inversion) << std::endl;
-// std::cout << "pumpPhotons "<< L_infty(pumpPhotons) << std::endl;
//temp : argument, welches in den exponenten zur verstärkung kommt
temp= emissionCrosssection * inversion * dz ;
- //std::cout << "arg exp temp "<< L_infty(temp) << std::endl;
- // e^(verstärkung)
- temp = Exp(temp);
-// std::cout << "exp temp "<< L_infty(temp) << std::endl;
-
-// std::ofstream DATEIC;
-// DATEIC.open("varGAIN.vtk");
-// temp.Print_VTK(DATEIC);
-// DATEIC.close();
- //verstärkungsschritt. hier eventuell sqrt(temp), da es sich hier nicht um die intensität, sondern sqrt(intensität) handelt?
+
+
+ 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;
@@ -229,6 +235,8 @@ void virtualLensCorrection(double dz, double wavenumberAverage,
VariableFFT& varWavenumber,
VariableFFT& temp)
{
+ std::cout << "virtualLensCorrection skpped\n";
+ return;
@@ -451,6 +459,16 @@ void estimateBeamQuality(VariableFFT& 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;
@@ -556,6 +574,14 @@ void estimateBeamQuality(VariableFFT& varIn,
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;
@@ -579,6 +605,8 @@ void estimateBeamQuality(VariableFFT& varIn,
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));
@@ -614,6 +642,9 @@ void estimateBeamQuality(VariableFFT& varIn,
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;
@@ -637,7 +668,7 @@ void estimateBeamQuality(VariableFFT& varIn,
int main(int argc, char** argv) {
std::ofstream DATEI;
- int n = 6;
+ int n = 8;
double R1 = 100;
double R2 = 100;
@@ -648,7 +679,7 @@ int main(int argc, char** argv) {
double distance = 50 ; //[mm]
double dz = distance / double(distanceIncrements);
lambda = 1064e-6; //[mm]
- radiusLoch = 0.2; //[mm]
+ radiusLoch = 0.1; //[mm]
distance = fabs(distanceFromWaist * 2.0);
@@ -692,10 +723,28 @@ int main(int argc, char** argv) {
std::cout << "curvature " << curvature << std::endl;
- 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 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);
@@ -714,6 +763,7 @@ int main(int argc, char** argv) {
Variable2D<double > varDeltaN(D2block);
Variable2D<double > varPumppower(D2blockPumppower);
+
IteratorZDirection zIterator(thermalRefractiveIndex3D->Give_blockgrid());
zIterator.next();
varDeltaN.interpolateSlizeZ(thermalRefractiveIndex3D,0.0);
@@ -738,9 +788,17 @@ int main(int argc, char** argv) {
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);
@@ -756,6 +814,10 @@ int main(int argc, char** argv) {
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);
@@ -770,10 +832,10 @@ int main(int argc, char** argv) {
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();
+// 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;
@@ -815,13 +877,16 @@ 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;
+ 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);
@@ -830,7 +895,6 @@ int main(int argc, char** argv) {
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)));
@@ -877,8 +941,31 @@ int main(int argc, char** argv) {
double powerHelper = product(absolute2(temp),absolute2(unit));
powerHelper = powerHelper *varE.getHx() * varE.getHy();
std::cout << "power before : "<<powerHelper << std::endl;
- for (int iter = 0 ; iter < zIterator.getSize();iter++)
+ 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;
@@ -890,27 +977,44 @@ int main(int argc, char** argv) {
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;
- varWavenumber = (n0 + deltaN) * 2.0 * M_PI / lambda;
+ 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);
- amplification(zIterator.get_hz(),pumppower,varIn,varWavenumber, temp);
- spectrumPlaneWavePropagation( zIterator.get_hz(),
- wavenumber,
- varIn,
- varE,Kx,Ky,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 = 10;
+
+
+ int plotevery = 100;
if (iter % plotevery == 0)
{
@@ -925,6 +1029,13 @@ int main(int argc, char** argv) {
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);
@@ -942,13 +1053,15 @@ int main(int argc, char** argv) {
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();