FFT example updated, 140W pump, fibercryst

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();