diff --git a/src/functionObjects/field/particleDistribution/particleDistribution.C b/src/functionObjects/field/particleDistribution/particleDistribution.C
index 0cd7827c2557fc3aae4b13d07d4cad70f85c24c1..71919b0933bb20a092f08de3ce27f533e476f12c 100644
--- a/src/functionObjects/field/particleDistribution/particleDistribution.C
+++ b/src/functionObjects/field/particleDistribution/particleDistribution.C
@@ -47,28 +47,6 @@ namespace functionObjects
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::functionObjects::particleDistribution::writeFileHeader
-(
-    Ostream& os
-) const
-{
-    writeHeader(os, "Particle distribution");
-    writeHeaderValue(os, "Cloud", cloudName_);
-    forAll(nameVsBinWidth_, i)
-    {
-        writeHeaderValue
-        (
-            os,
-            "Field:" + nameVsBinWidth_[i].first(),
-            nameVsBinWidth_[i].second()
-        );
-    }
-    os << endl;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::functionObjects::particleDistribution::particleDistribution
@@ -79,14 +57,14 @@ Foam::functionObjects::particleDistribution::particleDistribution
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
-    writeFile(runTime, name, typeName, dict),
+    writeFile(runTime, name),
     cloudName_("unknown-cloudName"),
     nameVsBinWidth_(),
     tagFieldName_("none"),
-    rndGen_(1234, -1)
+    rndGen_(1234, -1),
+    writerPtr_(nullptr)
 {
     read(dict);
-    writeFileHeader(file());
 }
 
 
@@ -105,6 +83,8 @@ bool Foam::functionObjects::particleDistribution::read(const dictionary& dict)
         dict.lookup("cloud") >> cloudName_;
         dict.lookup("nameVsBinWidth") >> nameVsBinWidth_;
         dict.readIfPresent("tagField", tagFieldName_);
+        word format(dict.lookup("setFormat"));
+        writerPtr_ = writer<scalar>::New(format);
 
         Info<< type() << " " << name() << " output:" << nl
             << "    Processing cloud : " << cloudName_ << nl
@@ -177,7 +157,6 @@ bool Foam::functionObjects::particleDistribution::write()
         }
     }
 
-    file() << "# Time: " << mesh_.time().timeName() << nl;
 
     bool ok = false;
     forAll(nameVsBinWidth_, i)
@@ -198,11 +177,6 @@ bool Foam::functionObjects::particleDistribution::write()
         }
     }
 
-    if (ok)
-    {
-        file() << nl;
-    }
-
     return true;
 }
 
@@ -220,11 +194,10 @@ void Foam::functionObjects::particleDistribution::generateDistribution
         return;
     }
 
-    Ostream& os = file();
-
+    word fName(fieldName);
     if (tag != -1)
     {
-        os  << tag << token::TAB;
+        fName = fName + '_' + Foam::name(tag);
     }
 
     distributionModels::general distribution
@@ -234,7 +207,31 @@ void Foam::functionObjects::particleDistribution::generateDistribution
         rndGen_
     );
 
-    os  << fieldName << distribution.writeDict(mesh_.time().timeName());
+    const Field<scalar> distX(distribution.x());
+    const Field<scalar> distY(distribution.y());
+
+    pointField xBin(distX.size(), Zero);
+    xBin.replace(0, distX);
+    const coordSet coords
+    (
+        fName,
+        "x",
+        xBin,
+        distX
+    );
+
+    const wordList fieldNames(1, fName);
+
+    fileName outputPath(baseTimeDir());
+    mkDir(outputPath);
+    OFstream graphFile(outputPath/writerPtr_->getFileName(coords, fieldNames));
+
+    Log << "    Writing distribution of " << fieldName
+        << " to " << graphFile.name() << endl;
+
+    List<const scalarField*> yPtrs(1);
+    yPtrs[0] = &distY;
+    writerPtr_->write(coords, fieldNames, yPtrs, graphFile);
 }
 
 
diff --git a/src/functionObjects/field/particleDistribution/particleDistribution.H b/src/functionObjects/field/particleDistribution/particleDistribution.H
index 37ece3aeb98928ce2979054620c35b2fd7e8d232..fec579f1938dfe5c6f8af601f72153c3d910accc 100644
--- a/src/functionObjects/field/particleDistribution/particleDistribution.H
+++ b/src/functionObjects/field/particleDistribution/particleDistribution.H
@@ -28,7 +28,7 @@ Group
     grpFieldFunctionObjects
 
 Description
-    Generates a particle distrbution for lagrangian data at a given time.
+    Generates a particle distribution for lagrangian data at a given time.
 
 Usage
     \verbatim
@@ -43,6 +43,7 @@ Usage
             (d 0.1)
             (U 10)
         );
+        setFormat   raw;
     }
     \endverbatim
 
@@ -53,6 +54,7 @@ Usage
         cloud       | Name of cloud to process | Yes        |
         nameVsBinWidth | List of cloud field vs bin width | Yes |
         tagField    | Name of cloud field to use to group particles | no | none
+        setFormat   | Output format            | yes         |
     \endtable
 
 See also
@@ -73,6 +75,7 @@ SourceFiles
 #include "scalarField.H"
 #include "cachedRandom.H"
 #include "Tuple2.H"
+#include "writer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -106,6 +109,9 @@ protected:
         //- Random number generator - used by distribution models
         cachedRandom rndGen_;
 
+        //- Writer
+        autoPtr<writer<scalar>> writerPtr_;
+
 
     // Protected Member Functions
 
@@ -133,9 +139,6 @@ protected:
         //- Disallow default bitwise assignment
         void operator=(const particleDistribution&) = delete;
 
-        //- Output file header information
-        virtual void writeFileHeader(Ostream& os) const;
-
 
 public:
 
diff --git a/src/lagrangian/distributionModels/general/general.C b/src/lagrangian/distributionModels/general/general.C
index 052bfbd1fb99e5169d019409dbc5b467057d24a6..87339642bf309e78dccdf806e06b4f820acedd26 100644
--- a/src/lagrangian/distributionModels/general/general.C
+++ b/src/lagrangian/distributionModels/general/general.C
@@ -247,17 +247,8 @@ Foam::dictionary Foam::distributionModels::general::writeDict
 {
     // dictionary dict = distributionModel::writeDict(dictName);
     dictionary dict(dictName);
-    List<scalar> data(xy_.size());
-    forAll(data, i)
-    {
-        data[i] = xy_[i][0];
-    }
-    dict.add("x", data);
-    forAll(data, i)
-    {
-        data[i] = xy_[i][1];
-    }
-    dict.add("y", data);
+    dict.add("x", x());
+    dict.add("y", y());
 
     return dict;
 }
@@ -280,6 +271,34 @@ void Foam::distributionModels::general::readDict(const dictionary& dict)
 }
 
 
+Foam::tmp<Foam::Field<Foam::scalar>>
+Foam::distributionModels::general::x() const
+{
+    tmp<Field<scalar>> tx(new Field<scalar>(xy_.size()));
+    scalarField& xi = tx.ref();
+    forAll(xy_, i)
+    {
+        xi[i] = xy_[i][0];
+    }
+
+    return tx;
+}
+
+
+Foam::tmp<Foam::Field<Foam::scalar>>
+Foam::distributionModels::general::y() const
+{
+    tmp<Field<scalar>> ty(new Field<scalar>(xy_.size()));
+    scalarField& yi = ty.ref();
+    forAll(xy_, i)
+    {
+        yi[i] = xy_[i][1];
+    }
+
+    return ty;
+}
+
+
 Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
diff --git a/src/lagrangian/distributionModels/general/general.H b/src/lagrangian/distributionModels/general/general.H
index 63d5d74d8911ddf614e35ded6781b23727a73d6e..db8a18194a411783329590075d82aea4eac8a105 100644
--- a/src/lagrangian/distributionModels/general/general.H
+++ b/src/lagrangian/distributionModels/general/general.H
@@ -38,6 +38,7 @@ SourceFiles
 #include "distributionModel.H"
 #include "Vector.H"
 #include "VectorSpace.H"
+#include "Field.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -121,6 +122,12 @@ public:
 
     // Member Functions
 
+        //- Bin boundaries
+        virtual tmp<Field<scalar>> x() const;
+
+        //- Probabilities
+        virtual tmp<Field<scalar>> y() const;
+
         //- Sample the distributionModel
         virtual scalar sample() const;
 
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/system/controlDict b/tutorials/lagrangian/sprayFoam/aachenBomb/system/controlDict
index b5bb5358985acc3098345e0c26ee78b35212fd74..d9b620638ab37f052228a816de511bcbad48248b 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/system/controlDict
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/system/controlDict
@@ -62,7 +62,7 @@ functions
             (d 1e-5)
             (U 10)
         );
-        distributionBinWidth 1e-5;
+        setFormat       raw;
     }
 }