lens creation function aufgerÃaeumt, beispiel fur fft angepasst (gerade fehler...

lens creation function aufgerÃaeumt, beispiel fur fft angepasst (gerade fehler drin, verstaerkung geht nicht

program2D/main.cc

@@ -31,6 +31,11 @@ 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);
 }
@@ -160,6 +165,40 @@ void spectrumPlaneWave(std::complex<double> (*aperture) (double x, double y),
           varFar.inversFFT();
       }
 
+void amplification(double dz,
+                   VariableFFT& pumppower,
+                   VariableFFT& varIn,
+                   VariableFFT& varWavenumber,
+                   VariableFFT& temp)
+{
+    double c = 3e10;// mm / s
+    double planck = 6.626e-34;
+    double emissionCrosssection = 7.7* 1e-18;//mm²
+
+    VariableFFT photonDensity(temp);
+    Function2d1<double, std::complex<double> > absolute(ABS);
+    photonDensity = absolute(varIn) * absolute(varIn)/ (planck * c / lambda) / dz;
+    temp = pumppower / (planck * c / lambda); // N /1 / mm³
+    // falsch : alle pumpphotonen werden zur verstärkung genutzt -> anpassen
+
+    std::cout << "pumppower "<< L_infty(pumppower) << std::endl;
+    std::cout << "inversion "<< L_infty(temp) << std::endl;
+    std::cout << "photonDensity "<< L_infty(photonDensity) << std::endl;
+    Function2d1<std::complex<double>,std::complex<double>> Exp(expComplex);
+    temp= emissionCrosssection * temp * photonDensity * dz   ;
+    std::cout << "arg exp temp "<< L_infty(temp) << std::endl;
+    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();
+
+    varIn = varIn * temp;
+
+
+}
 
 void virtualLensCorrection(double dz, double wavenumberAverage,
                            VariableFFT& varIn,
@@ -338,6 +377,7 @@ void estimateBeamWaist(VariableFFT& varIn,
     double power = product(Intensity,unity);
     double medianX = product(Intensity * X * X,unity) / power;
     double medianY = product(Intensity * Y * Y,unity) / power;
+    std::cout << "power = " << power <<std::endl;
     std::cout << "beamwaistX = " << 2.0 * sqrt(medianX) <<std::endl;
     std::cout << "beamwaistY = " << 2.0 * sqrt(medianY) <<std::endl;
 }
@@ -591,7 +631,7 @@ int main(int argc, char** argv) {
     distance = 150 ;
     dz = distance / double(distanceIncrements);
     double geometrySize  =  5.0; //[mm]
-    geometrySize =  2.0;
+    geometrySize =  1.0;
 
     radiusGauss = geometrySize / 4.0;   //[mm]
     radiusGauss = 0.2;
@@ -627,8 +667,10 @@ 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"));
+    VTK_Reader reader(         QString("/local/er96apow/FAUbox/Promotion/Vectorial_BPM/fibercryst_exmaple/RefractionIndexTherm_last.vtk"));
+    VTK_Reader readerPumppower(QString("/local/er96apow/FAUbox/Promotion/Vectorial_BPM/fibercryst_exmaple/abspower.vtk"));
     Variable<double> * thermalRefractiveIndex3D = reader.give_first_variable();
+    Variable<double> * pumpPowerRaytracing = readerPumppower.give_first_variable();
 
 
 //    Unstructured_grid *ug = new Cylinder(2,1,20);
@@ -638,10 +680,15 @@ int main(int argc, char** argv) {
 
 
     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);
@@ -651,6 +698,7 @@ int main(int argc, char** argv) {
     VariableFFT varE(n,n,geo);
 
     VariableFFT deltaN(varE);
+    VariableFFT pumppower(varE);
 
     deltaN.interpolate(vardDeltaNComplex,0.0);
 
@@ -795,8 +843,13 @@ int main(int argc, char** argv) {
         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);
+
 
 
 
@@ -808,6 +861,7 @@ int main(int argc, char** argv) {
 
 
         virtualLensCorrection(zIterator.get_hz(),wavenumberAverage_,varIn,varWavenumber, temp);
+        amplification(zIterator.get_hz(),pumppower,varIn,varWavenumber, temp);
         spectrumPlaneWavePropagation( zIterator.get_hz(),
                      wavenumber,
                      varIn,
@@ -820,7 +874,7 @@ int main(int argc, char** argv) {
         {
 
             //estimateBeamQuality(varIn,Kx,Ky,X,Y,temp,lambda);
-            //estimateBeamWaist(varIn,X,Y);
+            estimateBeamWaist(varIn,X,Y);
 
 
 
@@ -831,9 +885,9 @@ int main(int argc, char** argv) {
             varIntensity.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();
+//            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);