diff --git a/applications/solvers/additiveFoam/additiveFoam.C b/applications/solvers/additiveFoam/additiveFoam.C
index 5000591..b1afc92 100644
--- a/applications/solvers/additiveFoam/additiveFoam.C
+++ b/applications/solvers/additiveFoam/additiveFoam.C
@@ -42,6 +42,8 @@ Description
 #include "movingHeatSourceModel.H"
 #include "foamToExaCA/foamToExaCA.H"
 
+#include "interface/interfacePoints.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -67,12 +69,16 @@ int main(int argc, char *argv[])
 
     movingHeatSourceModel sources(mesh);
 
+    DynamicList<List<scalar>> interfaceData;
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
+        volScalarField R("R", mag(fvc::ddt(T)));
+
         #include "updateProperties.H"
 
         #include "readTimeControls.H"
@@ -94,6 +100,8 @@ int main(int argc, char *argv[])
         }
 
         #include "thermo/TEqn.H"
+
+        #include "interface/interface.H"
         
         ExaCA.update();
 
@@ -104,6 +112,8 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
+    #include "interface/interfaceWrite.H"
+
     ExaCA.write();
 
     return 0;
diff --git a/applications/solvers/additiveFoam/interface/interface.H b/applications/solvers/additiveFoam/interface/interface.H
new file mode 100644
index 0000000..f30460d
--- /dev/null
+++ b/applications/solvers/additiveFoam/interface/interface.H
@@ -0,0 +1,28 @@
+// write out interface dimensions
+Info<< "interface span " << "(" << Tliq.value() << "): "
+    << interfacePoints(T, Tliq.value()) << endl;
+
+Info<< "interface span " << "(" << Tsol.value() << "): "
+    << interfacePoints(T, Tsol.value()) << endl;
+
+// write out the maximum temperature
+Info<< "maximum temperature: " << gMax(T.internalField()) << endl;
+
+volScalarField G("G", mag(fvc::grad(T)));
+
+forAll(mesh.cells(), celli)
+{
+    if ( (T[celli] <= Tliq.value()) && (T.oldTime()[celli] > Tliq.value()) )
+    {
+        const point pt = mesh.C()[celli];
+
+        interfaceData.append(
+            {
+                pt[0],
+                pt[1],
+                pt[2],
+                R[celli] / G[celli],
+                G[celli]
+            });
+    }
+}
diff --git a/applications/solvers/additiveFoam/interface/interfacePoints.H b/applications/solvers/additiveFoam/interface/interfacePoints.H
new file mode 100644
index 0000000..7bb42da
--- /dev/null
+++ b/applications/solvers/additiveFoam/interface/interfacePoints.H
@@ -0,0 +1,81 @@
+vector interfacePoints(const volScalarField& field, const scalar& iso)
+{
+    const fvMesh& mesh = field.mesh();
+
+    // set local reference to mesh data
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+
+    const volVectorField& cc = mesh.C();
+
+    DynamicList<vector> positions(mesh.nFaces());
+
+    // check internal faces
+    for(label facei=0; facei < mesh.nInternalFaces(); facei++)
+    {        
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
+
+        scalar minFace = min(field[own], field[nei]);
+        scalar maxFace = max(field[own], field[nei]);
+
+        if ((minFace < iso) && (maxFace >= iso))
+        {
+            vector d = cc[nei] - cc[own];
+            vector p = cc[own] + d*(iso - field[own])/(field[nei] - field[own]);
+            positions.append(p);
+        }
+    }
+
+    // check boundary faces
+    const volScalarField::Boundary& fieldBf = field.boundaryField();
+
+    forAll(fieldBf, patchi)
+    {
+        const fvPatchScalarField& fieldPf = fieldBf[patchi];
+
+        const labelUList& faceCells = fieldPf.patch().faceCells();
+
+        if (fieldPf.coupled())
+        {
+            // processor boundary : interpolate across face
+            const vectorField ccn(cc.boundaryField()[patchi].patchNeighbourField());
+            const scalarField fn(fieldPf.patchNeighbourField());
+
+            forAll(faceCells, facei)
+            {
+                label own = faceCells[facei];
+
+                scalar minFace = min(field[own], fn[facei]);
+                scalar maxFace = max(field[own], fn[facei]);
+
+                if ((minFace < iso) && (maxFace >= iso))
+                {
+                    vector d = ccn[facei] -  cc[own];
+                    vector p = cc[own] + d*(iso - field[own])/(fn[facei] - field[own]);
+                    positions.append(p);
+                }
+            }
+        }
+        else
+        {
+            // physical boundary : take face centre in liquid
+            const vectorField fc(mesh.boundaryMesh()[patchi].faceCentres());
+            const scalarField pif(fieldPf.patchInternalField());
+
+            forAll(fc, i)
+            {
+                if (pif[i] >= iso)
+                {
+                    positions.append(fc[i]);
+                }
+            }
+        }
+    }
+
+    positions.shrink();
+
+    boundBox isoBb(positions);
+
+    return returnReduce(isoBb.span(), maxOp<vector>());
+}
diff --git a/applications/solvers/additiveFoam/interface/interfaceWrite.H b/applications/solvers/additiveFoam/interface/interfaceWrite.H
new file mode 100644
index 0000000..021d095
--- /dev/null
+++ b/applications/solvers/additiveFoam/interface/interfaceWrite.H
@@ -0,0 +1,22 @@
+const fileName interfacePath
+(
+    runTime.rootPath()/runTime.globalCaseName()/"solidificationData"
+);
+
+mkDir(interfacePath);
+
+OFstream os
+(
+    interfacePath + "/" + "data_" + Foam::name(Pstream::myProcNo()) + ".csv"
+);
+
+for(int i=0; i < interfaceData.size(); i++)
+{
+    int n = interfaceData[i].size()-1;
+
+    for(int j=0; j < n; j++)
+    {
+        os << interfaceData[i][j] << ",";
+    }
+    os << interfaceData[i][n] << endl;
+}
diff --git a/applications/solvers/additiveFoam/movingHeatSource/Make/files b/applications/solvers/additiveFoam/movingHeatSource/Make/files
index 044edf5..fda4be0 100644
--- a/applications/solvers/additiveFoam/movingHeatSource/Make/files
+++ b/applications/solvers/additiveFoam/movingHeatSource/Make/files
@@ -11,6 +11,7 @@ heatSourceModels/heatSourceModel/heatSourceModel.C
 heatSourceModels/heatSourceModel/heatSourceModelNew.C
 heatSourceModels/superGaussian/superGaussian.C
 heatSourceModels/modifiedSuperGaussian/modifiedSuperGaussian.C
+heatSourceModels/nLight/nLight.C
 
 movingHeatSourceModel/movingHeatSourceModel.C
 
diff --git a/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/heatSourceModel/heatSourceModel.C b/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/heatSourceModel/heatSourceModel.C
index 02a7d04..a9dad5d 100644
--- a/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/heatSourceModel/heatSourceModel.C
+++ b/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/heatSourceModel/heatSourceModel.C
@@ -239,7 +239,7 @@ Foam::heatSourceModel::qDot()
     );
     volScalarField& qDot_ = tqDot.ref();
 
-    // sample gaussian distribution at desired resolution
+    // sample distribution at desired resolution
     const scalar power_ = movingBeam_->power();
 
     if (power_ > small)
@@ -247,7 +247,7 @@ Foam::heatSourceModel::qDot()
         const vector position_ = movingBeam_->position();
 
         // udpate the absorbed power and heat source normalization term
-        const scalar aspectRatio = 
+        const scalar aspectRatio =
             dimensions_.z() / min(dimensions_.x(), dimensions_.y());
 
         dimensionedScalar absorbedPower
@@ -275,7 +275,7 @@ Foam::heatSourceModel::qDot()
         );
 
         const pointField& points = mesh_.points();
-        
+
         treeBoundBox beamBb
         (
             position_ - 1.5*dimensions_,
diff --git a/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/nLight/nLight.C b/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/nLight/nLight.C
new file mode 100644
index 0000000..686c0d7
--- /dev/null
+++ b/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/nLight/nLight.C
@@ -0,0 +1,247 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                Copyright (C) 2023 Oak Ridge National Laboratory                
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "nLight.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace heatSourceModels
+{
+    defineTypeNameAndDebug(nLight, 0);
+    addToRunTimeSelectionTable(heatSourceModel, nLight, dictionary);
+}
+}
+
+using Foam::constant::mathematical::pi;
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::heatSourceModels::nLight::nLight
+(
+    const word& sourceName,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    heatSourceModel(typeName, sourceName, dict, mesh),
+    mesh_(mesh),
+    spotCoeffs_(heatSourceModelCoeffs_.optionalSubDict("spotCoeffs")),
+    ringCoeffs_(heatSourceModelCoeffs_.optionalSubDict("ringCoeffs"))
+{
+    // Select parameters from specified mode
+    mode_ = heatSourceModelCoeffs_.lookup<label>("mode");
+
+    List<scalar> alphas_ =
+        {0.99, 0.83, 0.7, 0.54, 0.33, 0.21, 0.29};
+
+    List<scalar> spot4Sigma_ =
+        {156.99e-6, 158.23e-6, 165.31e-6, 182.28e-6, 181.15e-6, 214.63e-6, 402.04e-6};
+
+    List<scalar> ring4Sigma_ =
+        {48.49e-6, 112.15e-6, 116.07e-6, 135.32e-6, 116.74e-6, 112.46e-6, 116.62e-6};
+
+    dimensions_ = heatSourceModelCoeffs_.lookup<vector>("dimensions");
+
+    alpha_ = alphas_[mode_];
+
+    Info<< "Parameters for mode: " << mode_ << endl
+        << "\talpha: " << alpha_ << endl
+        << "\tD4s (spot): " << spot4Sigma_[mode_] << endl
+        << "\tD4s (ring): " << ring4Sigma_[mode_] << endl;
+    
+    //- Spot parameters
+    ks_ = spotCoeffs_.lookup<scalar>("k");
+    ms_ = spotCoeffs_.lookup<scalar>("m");
+    ds_ = vector(0.5*spot4Sigma_[mode_], 0.5*spot4Sigma_[mode_], dimensions_[2]);
+
+    //- Ring parameters
+    kr_ = ringCoeffs_.lookup<scalar>("k");
+    mr_ = ringCoeffs_.lookup<scalar>("m");
+    R_ = 152e-6;
+    r_ = 0.5 * ring4Sigma_[mode_];
+
+    // Calculated using second moment method from fitted profiles
+    List<scalar> D4sigma_ =
+        {162.12e-6, 235.39e-6, 283.97e-6, 283.97e-6,
+         340.0e-6, 385.18e-6, 413.13e-6, 438.85e-6};
+    
+    scalar R2sigma_ = 0.5*D4sigma_[mode_];    
+
+    //- Overwrite the dimensions used by the base class
+    dimensions_ = vector(R2sigma_, R2sigma_, ds_.z());
+    staticDimensions_ = dimensions_;
+}
+
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+inline Foam::scalar
+Foam::heatSourceModels::nLight::weight(const vector& d)
+{
+    // Spot parameters
+    const scalar as = Foam::pow(2.0, 1.0/ks_);
+    vector ss =
+        cmptDivide
+        (
+            vector(ds_.x(), ds_.y(), dimensions_.z()),
+            vector(as, as, 1.0)
+        );
+    
+    // Ring parameters
+    const scalar ar = Foam::pow(2.0, 1.0/kr_);
+    vector sr =
+        cmptDivide
+        (
+            vector(r_, r_, dimensions_.z()),
+            vector(ar, ar, 1.0)
+        );
+
+    // Check if cell is within heat source depth
+    if (d.z() < max(ss.z(), sr.z()))
+    {
+        // x-y distance from beam center
+        vector di = vector(d.x(), d.y(), 0.0);
+        
+        // Spot calcs
+        ss *= Foam::pow(1.0 - Foam::pow(d.z() / ss.z(), ms_), 1.0/ms_);
+
+        scalar xs = Foam::pow(magSqr(cmptDivide(di, ss)), ks_/2.0);
+        
+        // Ring calcs
+        const scalar rsz =
+            sr.x() * Foam::pow
+                     (
+                         1.0 - Foam::pow(d.z() / sr.z(), mr_),
+                         1.0 / mr_
+                     );
+        
+        scalar xr =
+            Foam::pow
+            (
+                Foam::sqr((mag(di) - R_) / rsz), 
+                kr_ / 2.0
+            );
+
+        // Return weights with power split and volume factored in
+        return alpha_ / As().value() / dimensions_.z() * Foam::exp(-xs) 
+            + (1.0 - alpha_) / Ar().value() / dimensions_.z() * Foam::exp(-xr);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+inline Foam::dimensionedScalar
+Foam::heatSourceModels::nLight::V0()
+{
+    //- Because of linear combination of shapes in nLight,
+    //  volume is integrated into the weight calculation
+    return dimensionedScalar("V0", dimVolume, 1.0);
+}
+
+inline Foam::dimensionedScalar
+Foam::heatSourceModels::nLight::As()
+{
+    const scalar a = Foam::pow(2.0, 1.0 / ks_);
+
+    const vector s =
+        cmptDivide
+        (
+            vector(ds_.x(), ds_.y(), dimensions_.z()),
+            vector(a, a, 1.0)
+        );
+
+    const dimensionedScalar As
+    (
+        "As",
+        dimArea,
+        s.x() * s.y() * pi * Foam::tgamma(1.0 + 2.0 / ks_)
+      * Foam::tgamma(1.0 + 1.0 / ms_) * Foam::tgamma(1.0 + 2.0 / ms_)
+      / Foam::tgamma(1.0 + 3.0 / ms_)
+    );
+
+    return As;
+}
+
+inline Foam::dimensionedScalar
+Foam::heatSourceModels::nLight::Ar()
+{
+    const scalar a = Foam::pow(2.0, 1.0 / kr_);
+    
+    const vector s =
+        cmptDivide
+        (
+            vector(r_, r_, dimensions_.z()),
+            vector(a, a, 1.0)
+        );
+    
+    const dimensionedScalar Ar
+    (
+        "Ar",
+        dimArea,
+        4.0 * pi * R_ * s.x() * Foam::tgamma(1.0 + 1.0 / kr_)
+      * Foam::tgamma(1.0 + 1.0 / mr_) * Foam::tgamma(1.0 + 1.0 / mr_)
+      / Foam::tgamma(1.0 + 2.0 / mr_)
+    );
+    
+    return Ar;
+}
+
+bool Foam::heatSourceModels::nLight::read()
+{
+    if (heatSourceModel::read())
+    {
+        heatSourceModelCoeffs_ = optionalSubDict(type() + "Coeffs");
+
+        //- General parameters
+        heatSourceModelCoeffs_.lookup("alpha") >> alpha_;
+        
+        //- Spot parameters
+        spotCoeffs_.lookup("k") >> ks_;
+        spotCoeffs_.lookup("m") >> ms_;
+        
+        //- Ring parameters
+        ringCoeffs_.lookup("k") >> kr_;
+        ringCoeffs_.lookup("m") >> kr_;
+        ringCoeffs_.lookup("R") >> R_;
+        ringCoeffs_.lookup("r") >> r_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/nLight/nLight.H b/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/nLight/nLight.H
new file mode 100644
index 0000000..6b1340b
--- /dev/null
+++ b/applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/nLight/nLight.H
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                Copyright (C) 2023 Oak Ridge National Laboratory                
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::heatSourceModels::nLight
+
+Description
+    Beam profile based on nLight system described in the listed reference.
+    The article describes a 2D profile which was extended to a generalized 3D
+    distribution by Jamie Stump and John Coleman for use with additiveFoam.
+    The exact profiles given in the nLight article may be approximated by
+    fitting the model coefficients in this class.
+
+    Reference:
+        Vanmunster, L., Goossens, L., Kinds, Y., van Hooreweder, B.,
+        & Vrancken, B. (2023). Effect of a variable laser beam profile 
+        system on productivity and surface quality of 316L stainless steel 
+        parts produced by Laser Powder Bed Fusion. CIRP Annals.
+        https://doi.org/10.1016/j.cirp.2023.04.035
+
+
+SourceFiles
+    nLight.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef nLight_H
+#define nLight_H
+
+#include "heatSourceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace heatSourceModels
+{
+
+/*---------------------------------------------------------------------------*\
+                                  Class nLight
+\*---------------------------------------------------------------------------*/
+
+class nLight
+:
+    public heatSourceModel
+{
+    // Private Data
+    
+        //- Pointer to mesh information
+        const fvMesh& mesh_;
+        
+        //- Spot/beam weight
+        scalar alpha_;
+        
+        //- Dictionary for spot shape parameters
+        const dictionary& spotCoeffs_;
+    
+        //- Laser mode [0-6]
+        scalar mode_;
+    
+        //- Spot distribution shape coefficient
+        scalar ks_;
+
+        //- Spot geometric shape coefficient
+        scalar ms_;
+        
+        //- Spot dimensions
+        vector ds_;
+                
+        //- Dictionary for ring shape parameters
+        const dictionary& ringCoeffs_;
+        
+        //- Ring distribution shape coefficient
+        scalar kr_;
+        
+        //- Ring geometric shape coefficient
+        scalar mr_;
+        
+        //- Ring center radius
+        scalar R_;
+        
+        //- Ring half width
+        scalar r_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("nLight");
+
+
+    // Constructors
+
+        //- Construct from components
+        nLight
+        (
+            const word& sourceName,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~nLight()
+    {}
+
+
+    // Member Functions
+
+        inline virtual scalar weight(const vector& d);
+
+        inline virtual dimensionedScalar V0();
+        
+        // Normalization factors for spot and ring integration
+        // NOTE: THESE DO NOT INCLUDE THE BEAM DEPTH!
+        inline virtual dimensionedScalar As();
+        inline virtual dimensionedScalar Ar();
+
+        //- Read the heatSourceProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace heatSourceModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/additiveFoam/movingHeatSource/movingHeatSourceModel/movingHeatSourceModel.C b/applications/solvers/additiveFoam/movingHeatSource/movingHeatSourceModel/movingHeatSourceModel.C
index ba998bc..6b7d804 100644
--- a/applications/solvers/additiveFoam/movingHeatSource/movingHeatSourceModel/movingHeatSourceModel.C
+++ b/applications/solvers/additiveFoam/movingHeatSource/movingHeatSourceModel/movingHeatSourceModel.C
@@ -57,7 +57,7 @@ Foam::movingHeatSourceModel::movingHeatSourceModel
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh_,
         dimensionedScalar(dimPower/dimVolume, 0.0)
diff --git a/tutorials/AMB2018-02-B/constant/heatSourceDict b/tutorials/AMB2018-02-B/constant/heatSourceDict
index 0122d92..1213cea 100644
--- a/tutorials/AMB2018-02-B/constant/heatSourceDict
+++ b/tutorials/AMB2018-02-B/constant/heatSourceDict
@@ -29,13 +29,40 @@ beam
     }
     
     heatSourceModel     superGaussian;
+    //heatSourceModel     nLight;
     
+    // Standard super-Gaussian spot profile
     superGaussianCoeffs
     {
         k               2.0;
         dimensions      (85.0e-6 85.0e-6 30e-6);
         nPoints         (10 10 10);
     }
+    
+    heatSourceModel     nLight;
+
+    // Profiles for nLight AFX-1000 lasers, supported modes 0-6.
+    nLightCoeffs
+    {
+        mode            5;
+        nPoints         (10 10 10);
+        dimensions      (78.5e-6 78.5e-6 20e-6);
+        
+        spotCoeffs
+        {
+            k           2.0;
+            m           3.0;
+        
+        }
+        ringCoeffs
+        {
+            k           2.0;
+            m           3.0;
+        }
+
+        transient       true;
+        isoValue        1726;
+    }
 }
 
 // ************************************************************************* //