Merge pull request ComputationalRadiationPhysics#851 from PrometheusPi/fix-phasePlaneWave

avoid non zero E-field integral in plane wave

Author: ax3l

src/picongpu/include/fields/laserProfiles/laserPlaneWave.hpp

@@ -19,9 +19,7 @@
  */
 
 
-
-#ifndef LASERPLANEWAVE_HPP
-#define	LASERPLANEWAVE_HPP
+#pragma once
 
 #include "types.h"
 #include "simulation_defines.hpp"
@@ -30,13 +28,38 @@ namespace picongpu
 {
     /** plane wave (use periodic boundaries!)
      *
-     *  no phase shifts, no spacial envelope
+     *  no transverse spacial envelope
+     *  based on the electric potential
+     *  Phi = Phi_0 * exp(0.5 * (x-x_0)^2 / sigma^2) * cos(k*(x - x_0) - phi)
+     *  by applying -grad Phi = -d/dx Phi = E(x)
+     *  we get:
+     *  E = -Phi_0 * exp(0.5 * (x-x_0)^2 / sigma^2) * [k*sin(k*(x - x_0) - phi) + x/sigma^2 * cos(k*(x - x_0) - phi)]
+     *
+     *  This approach ensures that int_{-infinity}^{+infinity} E(x) = 0 for any phase
+     *  if we have no transverse profile as we have with this plane wave train
+     *
+     *  Since PIConGPU requires a temporally defined electric field, we use:
+     *  t = x/c and (x-x_0)/sigma = (t-t_0)/tau and k*(x-x_0) = omega*(t-t_0) with omega/k = c and tau * c = sigma
+     *  and get:
+     *  E = -Phi_0*omega/c * exp(0.5 * (t-t_0)^2 / tau^2) * [sin(omega*(t - t_0) - phi) + t/(omega*tau^2) * cos(omega*(t - t_0) - phi)]
+     *  and define:
+     *    E_0 = -Phi_0*omega/c
+     *    integrationCorrectionFactor = t/(omega*tau^2)
+     *
+     *  Please consider:
+     *   1) The above formulae does only apply to a Gaussian envelope. If the plateau length is
+     *      not zero, the integral over the volume will only vanish if the plateau length is
+     *      a multiple of the wavelength.
+     *   2) Since we define our envelope by a sigma of the laser intensity, 
+     *      tau = PULSE_LENGTH / sqrt(2)
      */
     namespace laserPlaneWave
     {
-
-        /** Compute the
+        /** calculates longitudinal field distribution
          *
+         * @param currentStep
+         * @param phase
+         * @return
          */
         HINLINE float3_X laserLongitudinal( uint32_t currentStep, float_X& phase )
         {
@@ -46,48 +69,50 @@ namespace picongpu
             double envelope = double(AMPLITUDE );
             float3_X elong(float3_X::create(0.0));
 
-            // a NON-symmetric (starting with phase=0) pulse will be initialized at position z=0 for
-            // a time of RAMP_INIT * PULSE_LENGTH + LASER_NOFOCUS_CONSTANT = INIT_TIME.
-            // we shift the complete pulse for the half of this time to start with
-            // the front of the laser pulse.
             const double mue = 0.5 * RAMP_INIT * PULSE_LENGTH;
 
             const double w = 2.0 * PI * f;
+            const double tau = PULSE_LENGTH / sqrt( 2.0 );
 
             const double endUpramp = mue;
             const double startDownramp = mue + LASER_NOFOCUS_CONSTANT;
 
-
+            double integrationCorrectionFactor = 0.0;
 
             if( runTime > startDownramp )
             {
                 // downramp = end
-                const double exponent =
-                    ( ( runTime - startDownramp )
-                      / PULSE_LENGTH / sqrt( 2.0 ) );
+                const double exponent = (runTime - startDownramp) / tau;
                 envelope *= exp( -0.5 * exponent * exponent );
+                integrationCorrectionFactor = ( runTime - startDownramp )/ (w*tau*tau);
             }
             else if ( runTime < endUpramp )
             {
                 // upramp = start
-                const double exponent = ( ( runTime - endUpramp ) / PULSE_LENGTH / sqrt( 2.0 ) );
+                const double exponent = (runTime - endUpramp) / tau;
                 envelope *= exp( -0.5 * exponent * exponent );
-
+                integrationCorrectionFactor = ( runTime - endUpramp )/ (w*tau*tau);
             }
 
             const double timeOszi = runTime - endUpramp;
+            const double t_and_phase = w * timeOszi + LASER_PHASE;
+            // to understand both components [sin(...) + t/tau^2 * cos(...)] see description above
             if( Polarisation == LINEAR_X )
             {
-                elong.x() = float_X( envelope * math::sin( w * timeOszi + LASER_PHASE));
+              elong.x() = float_X( envelope * (math::sin(t_and_phase)
+                          + math::cos(t_and_phase) * integrationCorrectionFactor));
             }
             else if( Polarisation == LINEAR_Z)
             {
-                elong.z() = float_X( envelope * math::sin( w * timeOszi + LASER_PHASE));
+              elong.z() = float_X( envelope * (math::sin(t_and_phase)
+                          + math::cos(t_and_phase) * integrationCorrectionFactor));
             }
             else if( Polarisation == CIRCULAR )
             {
-                elong.x() = float_X( envelope / sqrt(2.0)  * math::sin( w * timeOszi + LASER_PHASE));
-                elong.z() = float_X( envelope / sqrt(2.0)  * math::cos( w * timeOszi + LASER_PHASE));
+                elong.x() = float_X( envelope / sqrt(2.0) * (math::sin(t_and_phase)
+                            + math::cos(t_and_phase) * integrationCorrectionFactor));
+                elong.z() = float_X( envelope / sqrt(2.0) * (math::cos(t_and_phase)
+                            - math::sin(t_and_phase) * integrationCorrectionFactor));
             }
 
 
@@ -96,7 +121,7 @@ namespace picongpu
             return elong;
         }
 
-        /**
+        /** calculates transverse field distribution
          *
          * @param elong
          * @param phase
@@ -112,7 +137,3 @@ namespace picongpu
     }
 }
 
-#endif	/* LASERPLANEWAVE_HPP */
-
-
-