4 files changed, 75 insertions, 0 deletions
diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package
index b825924c6..edcf402c9 100644
--- a/Source/Initialization/Make.package
+++ b/Source/Initialization/Make.package
@@ -1,4 +1,5 @@
 CEXE_sources += CustomDensityProb.cpp
+CEXE_sources += PlasmaProfiles.cpp
 CEXE_sources += WarpXInitData.cpp
 CEXE_sources += CustomMomentumProb.cpp
 CEXE_sources += PlasmaInjector.cpp
diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H
index f6b76e8b6..cd5802a55 100644
--- a/Source/Initialization/PlasmaInjector.H
+++ b/Source/Initialization/PlasmaInjector.H
@@ -10,6 +10,8 @@
 #include "AMReX_ParmParse.H"
 #include "AMReX_Utility.H"
 
+enum class predefined_profile_flag { Null, parabolic_channel };
+
 ///
 /// PlasmaDensityProfile describes how the charge density
 /// is set in particle initialization. Subclasses must define a
@@ -59,6 +61,24 @@ private:
 };
 
 ///
+/// This describes predefined density distributions.
+///
+class PredefinedDensityProfile : public PlasmaDensityProfile
+{
+public:
+    PredefinedDensityProfile(const std::string& species_name);
+    virtual amrex::Real getDensity(amrex::Real x,
+                                   amrex::Real y,
+                                   amrex::Real z) const override;
+    amrex::Real ParabolicChannel(amrex::Real x,
+                                 amrex::Real y,
+                                 amrex::Real z) const;
+private:
+    predefined_profile_flag which_profile = predefined_profile_flag::Null;
+    amrex::Vector<amrex::Real> params;
+};
+
+///
 /// This describes a density function parsed in the input file. 
 ///
 class ParseDensityProfile : public PlasmaDensityProfile
diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp
index 0da9318de..52b5d8b61 100644
--- a/Source/Initialization/PlasmaInjector.cpp
+++ b/Source/Initialization/PlasmaInjector.cpp
@@ -70,6 +70,17 @@ CustomDensityProfile::CustomDensityProfile(const std::string& species_name)
     pp.getarr("custom_profile_params", params);
 }
 
+PredefinedDensityProfile::PredefinedDensityProfile(const std::string& species_name)
+{
+    ParmParse pp(species_name);
+    std::string which_profile_s;
+    pp.getarr("predefined_profile_params", params);
+    pp.query("predefined_profile_name", which_profile_s);
+    if (which_profile_s == "parabolic_channel"){
+        which_profile = predefined_profile_flag::parabolic_channel;
+    }
+}
+
 ParseDensityProfile::ParseDensityProfile(std::string parse_density_function)
     : _parse_density_function(parse_density_function)
 {
@@ -333,6 +344,8 @@ void PlasmaInjector::parseDensity(ParmParse pp){
         rho_prof.reset(new ConstantDensityProfile(density));
     } else if (rho_prof_s == "custom") {
         rho_prof.reset(new CustomDensityProfile(species_name));
+    } else if (rho_prof_s == "predefined") {
+        rho_prof.reset(new PredefinedDensityProfile(species_name));
     } else if (rho_prof_s == "parse_density_function") {
         pp.get("density_function(x,y,z)", str_density_function);
         rho_prof.reset(new ParseDensityProfile(str_density_function));
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;
+}