1 files changed, 41 insertions, 0 deletions
diff --git a/Source/Initialization/PlasmaProfiles.cpp b/Source/Initialization/PlasmaProfiles.cpp
new file mode 100644
index 000000000..d9d207f7e
--- /dev/null
+++ b/Source/Initialization/PlasmaProfiles.cpp
@@ -0,0 +1,41 @@
+#include <PlasmaInjector.H>
+#include <cmath>
+#include <iostream>
+#include <WarpXConst.H>
+
+using namespace amrex;
+
+Real PredefinedDensityProfile::getDensity(Real x, Real y, Real z) const {
+    Real n;
+    if ( which_profile == predefined_profile_flag::parabolic_channel ) {
+        n = ParabolicChannel(x,y,z);
+    }
+    return n;
+}
+
+///
+/// plateau between linear upramp and downramp, and parab transverse profile
+///
+Real PredefinedDensityProfile::ParabolicChannel(Real x, Real y, Real z) const {
+    //  params = [z_start   ramp_up   plateau   ramp_down   rc       n0]
+    Real z_start   = params[0];
+    Real ramp_up   = params[1];
+    Real plateau   = params[2];
+    Real ramp_down = params[3];
+    Real rc        = params[4];
+    Real n0        = params[5];
+    Real n;
+    Real kp = PhysConst::q_e/PhysConst::c*sqrt( n0/(PhysConst::m_e*PhysConst::ep0) );
+
+    if        ((z-z_start)>=0               and (z-z_start)<ramp_up                  ) {
+        n = (z-z_start)/ramp_up;
+    } else if ((z-z_start)>=ramp_up         and (z-z_start)<ramp_up+plateau          ) {
+        n = 1;
+    } else if ((z-z_start)>=ramp_up+plateau and (z-z_start)<ramp_up+plateau+ramp_down) {
+        n = 1-((z-z_start)-ramp_up-plateau)/ramp_down;
+    } else {
+        n = 0;
+    }
+    n *= n0*(1+4*(x*x+y*y)/(kp*kp*std::pow(rc,4)));
+    return n;
+}