src/Algorithms/PartData.h

@@ -32,6 +32,8 @@
 //  by the compiler perform the correct operation.  For speed reasons
 //  they are not implemented.
 
+#include <Kokkos_Core.hpp>
+
 class PartData {
 
 public:
@@ -48,10 +50,10 @@ public:
     PartData();
 
     /// The constant charge per particle.
-    double getQ() const;
+    KOKKOS_INLINE_FUNCTION double getQ() const;
 
     /// The constant mass per particle.
-    double getM() const;
+    KOKKOS_INLINE_FUNCTION double getM() const;
 
     /// The constant reference momentum per particle.
     double getP() const;
@@ -101,12 +103,12 @@ protected:
 // Inline functions.
 // ------------------------------------------------------------------------
 
-inline double PartData::getQ() const {
+KOKKOS_INLINE_FUNCTION double PartData::getQ() const {
     return charge;
 }
 
 
-inline double PartData::getM() const {
+KOKKOS_INLINE_FUNCTION double PartData::getM() const {
     return mass;
 }
 

src/Algorithms/Tracker.cpp

@@ -63,7 +63,7 @@
 #include "Fields/BMultipoleField.h"
 
 // FIXME Remove headers and dynamic_cast in readOneBunchFromFile
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 
 #include <cfloat>
 #include <cmath>

src/Algorithms/Tracker.h

@@ -63,7 +63,7 @@
 #define CLASSIC_Tracker_HH
 
 #include "Algorithms/AbstractTracker.h"
-#include "PartBunch/PartBunch.hpp"
+#include "PartBunch/PartBunch.h"
 #include "Algorithms/PartData.h"
 
 #include "Utilities/ClassicField.h"

src/BasicActions/Option.h

@@ -20,6 +20,29 @@
 #ifndef OPAL_Option_HH
 #define OPAL_Option_HH
 
+#ifdef __CUDACC__
+#pragma push_macro("__cpp_consteval")
+#pragma push_macro("_NODISCARD")
+#pragma push_macro("__builtin_LINE")
+
+#define __cpp_consteval 201811L
+
+#ifdef _NODISCARD
+    #undef _NODISCARD
+    #define _NODISCARD
+#endif
+
+#define consteval constexpr
+
+#include <source_location>
+
+#undef consteval
+#pragma pop_macro("__cpp_consteval")
+#pragma pop_macro("_NODISCARD")
+#else
+#include <source_location>
+#endif
+
 #include "AbstractObjects/Action.h"
 #include "Utilities/Options.h"
 

src/BeamlineGeometry/Euclid3DGeometry.cpp

@@ -70,6 +70,10 @@ Euclid3D Euclid3DGeometry::getTransform(double /*fromS*/, double /*toS*/) const
     throw GeneralClassicException("Euclid3DGeometry::getTransform", "Not implemented");
 }
 
+Euclid3D Euclid3DGeometry::getTransform(double /*fromS*/) const {
+    throw GeneralClassicException("Euclid3DGeometry::getTransform", "Not implemented");
+}
+
 Euclid3D Euclid3DGeometry::getTotalTransform() const {
     return transformation_m;
 }
\ No newline at end of file

src/BeamlineGeometry/Euclid3DGeometry.h

@@ -70,6 +70,8 @@ class Euclid3DGeometry : public BGeometryBase {
     //  position [b]fromS[/b] to the position [b]toS[/b].
     virtual Euclid3D getTransform(double fromS, double toS) const;
 
+    virtual Euclid3D getTransform(double fromS) const;
+
     /// Get total transform from beginning to end
     //  Corresponds to the Euclid3D
     virtual Euclid3D getTotalTransform() const;

src/BeamlineGeometry/VarRadiusGeometry.cpp

@@ -45,3 +45,10 @@ Euclid3D VarRadiusGeometry::getTransform(double fromS, double toS) const {
     v.setZ(ref_to[1] - ref_from[1]);
     return v;
 }
+
+Euclid3D VarRadiusGeometry::getTransform(double fromS) const {
+    throw GeneralClassicException("Euclid3DGeometry::getTransform", "Not implemented");
+}
+
+
+

src/BeamlineGeometry/VarRadiusGeometry.h

@@ -116,6 +116,13 @@ public:
      *  to the entrance of the element
      *  Equivalent to getTransform(0.0, getEntrance())
      */
+    
+    virtual Euclid3D getTransform(double fromS) const;
+    /** Transform of the local coordinate system from the origin
+     *  to the entrance of the element
+     *  Equivalent to getTransform(0.0, getEntrance())
+     */
+    
     virtual Euclid3D getEntranceFrame() const;
     /** Transform of the local coordinate system from the origin
      *  to the exit of the element

src/CMakeLists.txt

@@ -89,14 +89,18 @@ add_library (libOPAL ${OPAL_SRCS})
 set_target_properties (libOPAL PROPERTIES OUTPUT_NAME OPAL)
 
 target_link_libraries (libOPAL
+    PUBLIC
     ${MPI_CXX_LIBRARIES}
-    )
+    ippl
+    Heffte
+)
 
 message (STATUS "OPAL_LIBS is: " ${CMAKE_BINARY_DIR})
 
 add_executable (opalx Main.cpp)
 target_link_libraries( opalx
-    ${CMAKE_BINARY_DIR}/src/libOPAL.a
+    PUBLIC
+    libOPAL
     ${KOKKOS_LIBRARY1}
     ${KOKKOS_LIBRARY2}
     ${KOKKOS_LIBRARY3}
@@ -113,6 +117,8 @@ target_link_libraries( opalx
     z
     ${CMAKE_DL_LIBS}
     ${MPI_CXX_LIBRARIES}
+#    ubsan
+#    asan
     )
 
 message (STATUS "src/CmakeLists IPPL_LIBRARY : ${IPPL_LIBRARY}")
@@ -133,4 +139,4 @@ set (HDRS
     changes.h
     )
 
-install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/src")
\ No newline at end of file
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/src")

src/Distribution/CMakeLists.txt

 set (_SRCS
     Distribution.cpp
+    FlatTop.cpp
+    Gaussian.cpp
     LaserProfile.cpp
+    MultiVariateGaussian.cpp
     )
 
 include_directories (
@@ -11,8 +14,11 @@ add_opal_sources(${_SRCS})
 
 set (HDRS
     Distribution.h
+    FlatTop.h
+    Gaussian.h
     LaserProfile.h
     MapGenerator.h
+    MultiVariateGaussian.h
     matrix_vector_operation.h
     )
 

src/Distribution/Distribution.cpp

@@ -59,10 +59,10 @@ using Base = ippl::ParticleBase<ippl::ParticleSpatialLayout<T, Dim>>;
 
 using view_type = typename ippl::detail::ViewType<ippl::Vector<double, Dim>, 1>::view_type;
 
-constexpr double SMALLESTCUTOFF = 1e-12;
-
 namespace DISTRIBUTION {
-    enum { TYPE, FNAME, SIGMAX, SIGMAY, SIGMAZ, SIGMAPX, SIGMAPY, SIGMAPZ, SIZE, CUTOFFPX, CUTOFFPY, CUTOFFPZ, CUTOFFX, CUTOFFY, CUTOFFLONG };
+    enum { TYPE, FNAME, SIGMAX, SIGMAY, SIGMAZ, SIGMAPX, SIGMAPY, SIGMAPZ, CORR,
+           CUTOFFPX, CUTOFFPY, CUTOFFPZ, CUTOFFX, CUTOFFY, CUTOFFLONG, CORRX, CORRY,
+           CORRZ, CORRT, SIGMAT, TPULSEFWHM, TRISE, TFALL, FTOSCAMPLITUDE, FTOSCPERIODS, EMITTED, SIZE };
 }
 
 /*
@@ -83,7 +83,7 @@ Distribution::Distribution()
         "The DISTRIBUTION statement defines data for the 6D particle distribution."),
       distrTypeT_m(DistributionType::NODIST) {
     itsAttr[DISTRIBUTION::TYPE] =
-        Attributes::makePredefinedString("TYPE", "Distribution type.", {"GAUSS", "FROMFILE"});
+        Attributes::makePredefinedString("TYPE", "Distribution type.", {"GAUSS", "MULTIVARIATEGAUSS", "FLATTOP", "FROMFILE"});
 
     itsAttr[DISTRIBUTION::FNAME] =
         Attributes::makeString("FNAME", "File for reading in 6D particle coordinates.", "");
@@ -96,6 +96,41 @@ Distribution::Distribution()
     itsAttr[DISTRIBUTION::SIGMAPY] = Attributes::makeReal("SIGMAPY", "SIGMApy", 0.0);
     itsAttr[DISTRIBUTION::SIGMAPZ] = Attributes::makeReal("SIGMAPZ", "SIGMApz", 0.0);
 
+    itsAttr[DISTRIBUTION::CORR] = Attributes::makeRealArray("CORR", "r correlation");
+
+    itsAttr[DISTRIBUTION::CUTOFFPX] = Attributes::makeReal("CUTOFFPX", "Distribution cutoff px dimension in units of sigma.", 3.0);
+    itsAttr[DISTRIBUTION::CUTOFFPY] = Attributes::makeReal("CUTOFFPY", "Distribution cutoff py dimension in units of sigma.", 3.0);
+    itsAttr[DISTRIBUTION::CUTOFFPZ] = Attributes::makeReal("CUTOFFPZ", "Distribution cutoff pz dimension in units of sigma.", 3.0);
+
+    itsAttr[DISTRIBUTION::CUTOFFX] = Attributes::makeReal("CUTOFFX", "Distribution cutoff x direction in units of sigma.", 3.0);
+    itsAttr[DISTRIBUTION::CUTOFFY] = Attributes::makeReal("CUTOFFY", "Distribution cutoff r direction in units of sigma.", 3.0);
+    itsAttr[DISTRIBUTION::CUTOFFLONG] = Attributes::makeReal("CUTOFFLONG", "Distribution cutoff z or t direction in units of sigma.", 3.0);
+
+    itsAttr[DISTRIBUTION::CORRX] = Attributes::makeReal("CORRX", "x/px correlation, (R12 in transport notation).", 0.0);
+    itsAttr[DISTRIBUTION::CORRY] = Attributes::makeReal("CORRY", "y/py correlation, (R34 in transport notation).", 0.0);
+    itsAttr[DISTRIBUTION::CORRZ] = Attributes::makeReal("CORRZ", "z/pz correlation, (R56 in transport notation).", 0.0);
+    itsAttr[DISTRIBUTION::CORRT] = Attributes::makeReal("CORRT", "t/pt correlation, (R56 in transport notation).", 0.0);
+
+    itsAttr[DISTRIBUTION::SIGMAT] = Attributes::makeReal("SIGMAT", "SIGMAt (m)", 0.0);
+    itsAttr[DISTRIBUTION::TPULSEFWHM] = Attributes::makeReal("TPULSEFWHM", "Pulse FWHM for emitted distribution.", 0.0);
+    itsAttr[DISTRIBUTION::TRISE] = Attributes::makeReal("TRISE", "Rise time for emitted distribution.", 0.0);
+    itsAttr[DISTRIBUTION::TFALL] = Attributes::makeReal("TFALL", "Fall time for emitted distribution.", 0.0);
+
+    itsAttr[DISTRIBUTION::FTOSCAMPLITUDE]
+        = Attributes::makeReal("FTOSCAMPLITUDE", "Amplitude of oscillations superimposed "
+                               "on flat top portion of emitted GAUSS "
+                               "distribtuion (in percent of flat top "
+                               "amplitude)",0.0);
+
+    itsAttr[DISTRIBUTION::FTOSCPERIODS]
+        = Attributes::makeReal("FTOSCPERIODS", "Number of oscillations superimposed on "
+                               "flat top portion of emitted GAUSS "
+                               "distribution", 0.0);
+
+    itsAttr[DISTRIBUTION::EMITTED]
+        = Attributes::makeBool("EMITTED", "Emitted beam, from cathode, as opposed to "
+                               "an injected beam.", false);
+
     registerOwnership(AttributeHandler::STATEMENT);
 }
 
@@ -165,7 +200,19 @@ Inform& Distribution::printInfo(Inform& os) const {
     if (OpalData::getInstance()->inRestartRun()) {
         os << "* In restart. Distribution read in from .h5 file." << endl;
     } else {
-        printDistGauss(os);
+        switch (distrTypeT_m) {
+            case DistributionType::GAUSS:
+                printDistGauss(os);
+                break;
+            case DistributionType::MULTIVARIATEGAUSS:
+                printDistMultiVariateGauss(os);
+                break;
+            case DistributionType::FLATTOP:
+                printDistFlatTop(os);
+                break;
+            default:
+                throw OpalException("Distribution Param", "Unknown \"TYPE\" of \"DISTRIBUTION\"");
+         }
         os << "* " << endl;
         os << "* Distribution is injected." << endl;
     }
@@ -180,6 +227,14 @@ void Distribution::setAvrgPz(double avrgpz){
     avrgpz_m = avrgpz;
 }
 
+void Distribution::setTEmission(double tEmission) {
+        tEmission_m = tEmission;
+}
+
+double Distribution::getTEmission() const {
+        return tEmission_m;
+}
+
 void Distribution::setDistParametersGauss() {
     /*
      * Set distribution parameters. Do all the necessary checks depending
@@ -207,8 +262,191 @@ void Distribution::setDistParametersGauss() {
 
     setSigmaR_m();
     setSigmaP_m();
+
     avrgpz_m = 0.0;
 }
+
+void Distribution::setDistParametersMultiVariateGauss() {
+
+    cutoffR_m = 3.;
+    cutoffP_m = 3.;
+
+    // initialize the covariance matrix to identity
+    for (unsigned int i = 0; i < 6; ++ i) {
+        for (unsigned int j = 0; j < 6; ++ j) {
+            if (i==j)
+               correlationMatrix_m[i][j] = 1.0;
+            else
+               correlationMatrix_m[i][j] = 0.0;
+        }
+    }
+
+    // set diagonal elements first
+    setSigmaR_m();
+    setSigmaP_m();
+
+    for (unsigned int i = 0; i < 3; ++ i){
+        correlationMatrix_m[2*i  ][2*i  ] = sigmaR_m[i]*sigmaR_m[i];
+        correlationMatrix_m[2*i+1][2*i+1] = sigmaP_m[i]*sigmaP_m[i];
+    }
+
+    std::vector<double> cr = Attributes::getRealArray(itsAttr[DISTRIBUTION::CORR]);
+
+    if (!cr.empty()) {
+            // read off-diagonal correlation matrix from input file
+            if (cr.size() == 15) {
+                *gmsg << "* Use r to specify correlations" << endl;
+                unsigned int k = 0;
+                for (unsigned int i = 0; i < 5; ++ i) {
+                    for (unsigned int j = i + 1; j < 6; ++ j, ++ k) {
+                        correlationMatrix_m[j][i] = cr.at(k)*cr.at(k);
+                        correlationMatrix_m[i][j] = correlationMatrix_m[j][i]; // impose symmetry
+                    }
+                }
+            }
+            else {
+                throw OpalException("Distribution::SetDistParametersGauss",
+                                    "Inconsistent set of correlations specified, check manual");
+            }
+    }
+
+    avrgpz_m = 0.0;
+}
+
+void Distribution::setDistParametersFlatTop() {
+
+    cutoffR_m = 3.;
+    cutoffP_m = 3.;
+
+    // set diagonal elements first
+    setSigmaR_m();
+    setSigmaP_m();
+
+    // initialize the covariance matrix to identity
+    for (unsigned int i = 0; i < 6; ++ i) {
+        for (unsigned int j = 0; j < 6; ++ j) {
+            if (i==j)
+               correlationMatrix_m[i][j] = 1.0;
+            else
+               correlationMatrix_m[i][j] = 0.0;
+        }
+    }
+
+    cutoffR_m = ippl::Vector<double, 3>(
+        std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::CUTOFFX])),
+        std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::CUTOFFY])),
+        std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::CUTOFFLONG])));
+
+    correlationMatrix_m[1][0] = Attributes::getReal(itsAttr[DISTRIBUTION::CORRX]);
+    correlationMatrix_m[3][2] = Attributes::getReal(itsAttr[DISTRIBUTION::CORRY]);
+    correlationMatrix_m[5][4] = Attributes::getReal(itsAttr[DISTRIBUTION::CORRT]);
+
+    if (Attributes::getReal(itsAttr[DISTRIBUTION::CORRZ]) != 0.0)
+        correlationMatrix_m[5][4] = Attributes::getReal(itsAttr[DISTRIBUTION::CORRZ]);
+
+    emitting_m = Attributes::getBool(itsAttr[DISTRIBUTION::EMITTED]);
+
+    if (emitting_m) {
+        sigmaR_m[2] = 0.0;
+
+        sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::SIGMAT]));
+        sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::SIGMAT]));
+        tPulseLengthFWHM_m = std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::TPULSEFWHM]));
+
+        FTOSCAmplitude_m = std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::FTOSCAMPLITUDE]));
+        FTOSCPeriods_m = std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::FTOSCPERIODS]));
+
+        // If TRISE and TFALL are defined > 0.0 then these attributes
+        // override SIGMAT.
+        //
+        if (std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::TRISE])) > 0.0
+            || std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::TFALL])) > 0.0) {
+
+            double timeRatio = std::sqrt(2.0 * std::log(10.0)) - std::sqrt(2.0 * std::log(10.0 / 9.0));
+
+            if (std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::TRISE])) > 0.0)
+                sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::TRISE]))
+                    / timeRatio;
+
+            if (std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::TFALL])) > 0.0)
+                sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::TFALL]))
+                    / timeRatio;
+        }
+
+        // For an emitted beam, the longitudinal cutoff >= 0.
+        cutoffR_m[2] = std::abs(cutoffR_m[2]);
+
+    }
+    /*
+    cutoffR_m = Vector_t(Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFX]),
+                         Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFY]),
+                         Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFLONG]));
+
+    correlationMatrix_m(1, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRX]);
+    correlationMatrix_m(3, 2) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRY]);
+    correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRT]);
+
+    // CORRZ overrides CORRT.
+    if (Attributes::getReal(itsAttr[Attrib::Distribution::CORRZ]) != 0.0)
+        correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRZ]);
+
+    setSigmaR_m();
+    if (emitting_m) {
+        sigmaR_m[2] = 0.0;
+
+        sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT]));
+        sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT]));
+
+        tPulseLengthFWHM_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TPULSEFWHM]));
+
+        //
+         // If TRISE and TFALL are defined > 0.0 then these attributes
+         // override SIGMAT.
+         //
+        if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE])) > 0.0
+            || std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL])) > 0.0) {
+
+            double timeRatio = std::sqrt(2.0 * std::log(10.0)) - std::sqrt(2.0 * std::log(10.0 / 9.0));
+
+            if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE])) > 0.0)
+                sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE]))
+                    / timeRatio;
+
+            if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL])) > 0.0)
+                sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL]))
+                    / timeRatio;
+
+        }
+
+        // For an emitted beam, the longitudinal cutoff >= 0.
+        cutoffR_m[2] = std::abs(cutoffR_m[2]);
+    }
+
+    // Set laser profile/
+    laserProfileFileName_m = Attributes::getString(itsAttr[Attrib::Distribution::LASERPROFFN]);
+    if (!(laserProfileFileName_m == std::string(""))) {
+        laserImageName_m = Attributes::getString(itsAttr[Attrib::Distribution::IMAGENAME]);
+        laserIntensityCut_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::INTENSITYCUT]));
+        short flags = 0;
+        if (Attributes::getBool(itsAttr[Attrib::Distribution::FLIPX])) flags |= LaserProfile::FLIPX;
+        if (Attributes::getBool(itsAttr[Attrib::Distribution::FLIPY])) flags |= LaserProfile::FLIPY;
+        if (Attributes::getBool(itsAttr[Attrib::Distribution::ROTATE90])) flags |= LaserProfile::ROTATE90;
+        if (Attributes::getBool(itsAttr[Attrib::Distribution::ROTATE180])) flags |= LaserProfile::ROTATE180;
+        if (Attributes::getBool(itsAttr[Attrib::Distribution::ROTATE270])) flags |= LaserProfile::ROTATE270;
+
+        laserProfile_m = new LaserProfile(laserProfileFileName_m,
+                                          laserImageName_m,
+                                          laserIntensityCut_m,
+                                          flags);
+    }
+
+    // Legacy for ASTRAFLATTOPTH.
+    if (distrTypeT_m == DistributionType::ASTRAFLATTOPTH)
+        tRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE]));
+*/
+
+}
+
 void Distribution::createDistributionGauss(size_t numberOfParticles, double massIneV, ippl::ParticleAttrib<ippl::Vector<double, 3>>& R, ippl::ParticleAttrib<ippl::Vector<double, 3>>& P, std::shared_ptr<ParticleContainer_t> &pc, std::shared_ptr<FieldContainer_t> &fc, Vector_t<double, 3> nr) {
     // moved to Gaussian.hpp
 }
@@ -227,9 +465,54 @@ void Distribution::printDistGauss(Inform& os) const {
     os << "* SIGMAPZ    = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
 }
 
+void Distribution::printDistMultiVariateGauss(Inform& os)  const {
+    os << "* Distribution type: MULTIVARIATEGAUSS" << endl;
+    os << "* " << endl;
+    os << "* SIGMAX     = " << sigmaR_m[0] << " [m]" << endl;
+    os << "* SIGMAY     = " << sigmaR_m[1] << " [m]" << endl;
+    os << "* SIGMAZ     = " << sigmaR_m[2] << " [m]" << endl;
+    os << "* SIGMAPX    = " << sigmaP_m[0] << " [Beta Gamma]" << endl;
+    os << "* SIGMAPY    = " << sigmaP_m[1] << " [Beta Gamma]" << endl;
+    os << "* SIGMAPZ    = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
+
+    os << "* input cov matrix = ";
+    for (unsigned int i = 0; i < 6; ++ i) {
+        for (unsigned int j = 0; j < 6; ++ j) {
+            os << correlationMatrix_m[i][j] << " ";
+        }
+        os << endl << "                     ";
+    }
+}
+
+void Distribution::printDistFlatTop(Inform& os)  const {
+    os << "* Distribution type: FLATTOP" << endl;
+    os << "* " << endl;
+    os << "* SIGMAX     = " << sigmaR_m[0] << " [m]" << endl;
+    os << "* SIGMAY     = " << sigmaR_m[1] << " [m]" << endl;
+
+    if (emitting_m) {
+            os << "* Sigma Time Rise               = " << sigmaTRise_m
+               << " [sec]" << endl;
+            os << "* TPULSEFWHM                    = " << tPulseLengthFWHM_m
+               << " [sec]" << endl;
+            os << "* Sigma Time Fall               = " << sigmaTFall_m
+               << " [sec]" << endl;
+            os << "* Longitudinal cutoff           = " << cutoffR_m[2]
+               << " [units of Sigma Time]" << endl;
+            //os << "* Flat top modulation amplitude = "
+            //   << Attributes::getReal(itsAttr[DISTRIBUTION::FTOSCAMPLITUDE])
+            //   << " [Percent of distribution amplitude]" << endl;
+            //os << "* Flat top modulation periods   = "
+            //   << std::abs(Attributes::getReal(itsAttr[DISTRIBUTION::FTOSCPERIODS]))
+            //   << endl;
+    }
+    else{
+        os << "* SIGMAZ                        = " << sigmaR_m[2] << " [m]" << endl;
+    }
+}
+
 void Distribution::setAttributes() {
-    setSigmaR_m();
-    setSigmaP_m();
+    setDist();
 }
 
 void Distribution::setDist() {
@@ -240,6 +523,12 @@ void Distribution::setDist() {
         case DistributionType::GAUSS:
             setDistParametersGauss();
             break;
+        case DistributionType::MULTIVARIATEGAUSS:
+            setDistParametersMultiVariateGauss();
+            break;
+        case DistributionType::FLATTOP:
+            setDistParametersFlatTop();
+            break;
         default:
             throw OpalException("Distribution Param", "Unknown \"TYPE\" of \"DISTRIBUTION\"");
     }
@@ -247,7 +536,11 @@ void Distribution::setDist() {
 
 void Distribution::setDistType() {
     static const std::map<std::string, DistributionType> typeStringToDistType_s = {
-        {"NODIST", DistributionType::NODIST}, {"GAUSS", DistributionType::GAUSS}};
+        {"NODIST", DistributionType::NODIST},
+        {"GAUSS", DistributionType::GAUSS},
+        {"MULTIVARIATEGAUSS", DistributionType::MULTIVARIATEGAUSS},
+        {"FLATTOP", DistributionType::FLATTOP}
+    };
 
     distT_m = Attributes::getString(itsAttr[DISTRIBUTION::TYPE]);
 

src/Distribution/Distribution.h

@@ -49,7 +49,7 @@ class Beam;
 class Beamline;
 class H5PartWrapper;
 
-enum class DistributionType : short { NODIST = -1, GAUSS };
+enum class DistributionType : short { NODIST = -1, GAUSS, MULTIVARIATEGAUSS, FLATTOP, FROMFILE };
 
 using ParticleContainer_t = ParticleContainer<double, 3>;
 using FieldContainer_t = FieldContainer<double, 3>;
@@ -60,6 +60,8 @@ public:
 
     virtual ~Distribution();
 
+    using Matrix_t = ippl::Vector< ippl::Vector<double, 6>, 6>;
+
     virtual bool canReplaceBy(Object* object);
     virtual Distribution* clone(const std::string& name);
     virtual void execute();
@@ -88,6 +90,13 @@ public:
     ippl::Vector<double, 3> getSigmaR() const;
     ippl::Vector<double, 3> getSigmaP() const;
 
+    double getSigmaTRise() const;
+    double getSigmaTFall() const;
+    double getTPulseLengthFWHM() const;
+
+    double getFTOSCAmplitude() const;
+    double getFTOSCPeriods() const;
+
     void setDistType();
 
     void setDist();
@@ -99,6 +108,14 @@ public:
     ippl::Vector<double, 3> getCutoffR() const;
     ippl::Vector<double, 3> getCutoffP() const;
 
+    Matrix_t correlationMatrix_m;
+
+    bool emitting_m;                     /// Distribution is an emitted, and is currently
+                                         /// emitting, rather than an injected, beam.
+
+    double getTEmission() const;
+    void setTEmission(double tEmission);
+
 private:
     enum class EmissionModel : unsigned short { NONE, ASTRA, NONEQUIL };
 
@@ -166,11 +183,17 @@ private:
     // void initializeBeam(PartBunch_t* beam);
     void printDist(Inform& os, size_t numberOfParticles) const;
     void printDistGauss(Inform& os) const;
+    void printDistMultiVariateGauss(Inform& os) const;
+    void printDistFlatTop(Inform& os) const;
 
     void setAttributes();
 
     void setDistParametersGauss();
 
+    void setDistParametersMultiVariateGauss();
+
+    void setDistParametersFlatTop();
+
     void setSigmaR_m();
 
     void setSigmaP_m();
@@ -188,9 +211,18 @@ private:
 
     ippl::Vector<double, 3> pmean_m, xmean_m, sigmaR_m, sigmaP_m, cutoffR_m, cutoffP_m;
 
+    double sigmaTRise_m;
+    double sigmaTFall_m;
+    double tPulseLengthFWHM_m;
+
     DistributionType distrTypeT_m;
 
     double avrgpz_m;
+
+    double FTOSCAmplitude_m;
+    double FTOSCPeriods_m;
+
+    double tEmission_m;
 };
 
 inline Inform& operator<<(Inform& os, const Distribution& d) {
@@ -221,6 +253,26 @@ inline ippl::Vector<double, 3> Distribution::getCutoffP() const {
     return cutoffP_m;
 }
 
+inline double Distribution::getSigmaTRise() const {
+    return sigmaTRise_m;
+}
+
+inline double Distribution::getSigmaTFall() const {
+    return sigmaTFall_m;
+}
+
+inline double Distribution::getTPulseLengthFWHM() const {
+    return tPulseLengthFWHM_m;
+}
+
+inline double Distribution::getFTOSCAmplitude() const {
+    return FTOSCAmplitude_m;
+}
+
+inline double Distribution::getFTOSCPeriods() const {
+    return FTOSCPeriods_m;
+}
+
 inline DistributionType Distribution::getType() const {
     return distrTypeT_m;
 }

src/Distribution/SamplingBase.hpp

@@ -9,18 +9,32 @@ using FieldContainer_t = FieldContainer<double, 3>;
 using Distribution_t = Distribution;
 
 class SamplingBase{
-    protected:
-        std::shared_ptr<ParticleContainer_t> pc_m;
-        std::shared_ptr<FieldContainer_t> fc_m;
-        std::shared_ptr<Distribution_t> opalDist_m;
-        std::string samplingMethod_m;
-    public:
-        SamplingBase(std::shared_ptr<ParticleContainer_t> &pc, std::shared_ptr<FieldContainer_t> &fc, std::shared_ptr<Distribution_t> &dist)
-           : pc_m(pc), fc_m(fc), opalDist_m(dist) {}
-
-       virtual ~SamplingBase() {}
-
-       virtual void generateParticles(size_t& numberOfParticles, Vector_t<double, 3> nr) {}
+protected:
+    std::shared_ptr<ParticleContainer_t> pc_m;
+    std::shared_ptr<FieldContainer_t> fc_m;
+    std::shared_ptr<Distribution_t> opalDist_m;
+    std::string samplingMethod_m;
+public:
+    
+    SamplingBase(std::shared_ptr<ParticleContainer_t> &pc, std::shared_ptr<FieldContainer_t> &fc, std::shared_ptr<Distribution_t> &dist)
+        : pc_m(pc), fc_m(fc), opalDist_m(dist) {
+    }
+    
+    virtual ~SamplingBase() {}
+
+    virtual void generateParticles(size_t& numberOfParticles, Vector_t<double, 3> nr) {}
+
+    virtual void emitParticles(double t, double dt) {}
+
+    // testNumEmitParticles is purely made for testing and should be removed
+    virtual void testNumEmitParticles(size_t nsteps, double dt) {}
+
+    // testEmitParticles is purely made for testing and should be removed
+    virtual void testEmitParticles(size_t nsteps, double dt) {}
+
+    virtual void initDomainDecomp(double BoxIncr) {}
+
+    virtual void setWithDomainDecomp(bool withDomainDecomp) {}
 };
 #endif