diff --git a/src/lagrangian/dieselSpray/Make/files b/src/lagrangian/dieselSpray/Make/files
index 86dd21c0bcd473295a4737b9821555c70fd0af84..7fd1ca7f2a48980a44c05565ed59f2dd78d908cb 100644
--- a/src/lagrangian/dieselSpray/Make/files
+++ b/src/lagrangian/dieselSpray/Make/files
@@ -26,6 +26,7 @@ $(injector)/injector/injector.C
 $(injector)/injector/injectorIO.C
 $(injector)/injectorType/injectorType.C
 $(injector)/unitInjector/unitInjector.C
+$(injector)/multiHoleInjector/multiHoleInjector.C
 $(injector)/commonRailInjector/commonRailInjector.C
 $(injector)/swirlInjector/swirlInjector.C
 $(injector)/definedInjector/definedInjector.C
diff --git a/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.C b/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.C
index 9249c8ac65cba2198cc7abaee228c32399b377bc..0f00cb79eb9ddb199b687d45177e81d7580a0440 100644
--- a/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.C
+++ b/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.C
@@ -187,13 +187,14 @@ Foam::label Foam::commonRailInjector::nParcelsToInject
     return nParcels;
 }
 
-const Foam::vector Foam::commonRailInjector::position() const
+const Foam::vector Foam::commonRailInjector::position(const label n) const
 {
     return position_;
 }
 
 Foam::vector Foam::commonRailInjector::position
 (
+    const label n,
     const scalar time,
     const bool twoD,
     const scalar angleOfWedge,
@@ -236,12 +237,21 @@ Foam::vector Foam::commonRailInjector::position
     return position_;
 }
 
+Foam::label Foam::commonRailInjector::nHoles() const
+{
+    return 1;
+}
+
 Foam::scalar Foam::commonRailInjector::d() const
 {
     return d_;
 }
 
-const Foam::vector& Foam::commonRailInjector::direction() const
+const Foam::vector& Foam::commonRailInjector::direction
+(
+    const label i,
+    const scalar time
+) const
 {
     return direction_;
 }
@@ -358,7 +368,7 @@ void Foam::commonRailInjector::correctProfiles
 
     forAll(velocityProfile_, i)
     {
-        scalar Pinj = getTableValue(injectionPressureProfile_, massFlowRateProfile_[i][0]);
+        scalar Pinj = getTableValue(injectionPressureProfile_, velocityProfile_[i][0]);
         scalar Vinj = sqrt(2.0*(Pinj - referencePressure)/rho);
         scalar mfr = massFlowRateProfile_[i][1]/(rho*A);
         scalar Cd = mfr/Vinj;
@@ -367,4 +377,14 @@ void Foam::commonRailInjector::correctProfiles
     }
 }
 
+Foam::vector Foam::commonRailInjector::tan1(const label n) const
+{
+    return tangentialInjectionVector1_;
+}
+
+Foam::vector Foam::commonRailInjector::tan2(const label n) const
+{
+    return tangentialInjectionVector2_;
+}
+
 // ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.H b/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.H
index fcb222bc879b644c33992e0c2f6bfdaded3e6cf3..c829799029c84e7761308304d52113de71258ef6 100644
--- a/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.H
+++ b/src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.H
@@ -132,11 +132,12 @@ public:
         ) const;
 
         //- Return the injection position
-        const vector position() const;
+        const vector position(const label n) const;
 
         //- Return the injection position
         vector position
         (
+            const label n,
             const scalar time,
             const bool twoD,
             const scalar angleOfWedge,
@@ -146,11 +147,18 @@ public:
             Random& rndGen
         ) const;
     
+        //- Return the number of holes
+        label nHoles() const;
+
         //- Return the injector diameter
         scalar d() const;
 
         //- Return the injection direction
-        const vector& direction() const;
+        const vector& direction
+        (
+            const label i,
+            const scalar time
+        ) const;
 
         //- Return the mass of the injected particle
         scalar mass
@@ -206,6 +214,9 @@ public:
         List<pair> CdProfile() const;
         scalar Cd(const scalar time) const;
 
+        vector tan1(const label n) const;
+        vector tan2(const label n) const;
+
         void correctProfiles
         (
             const liquidMixture& fuel,
diff --git a/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.C b/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.C
index 89960f39782b03a15fd18188d71a1139a3544010..8d422973cfa8c640d692bcb1c22ef6ab63041490 100644
--- a/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.C
+++ b/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.C
@@ -178,13 +178,14 @@ Foam::label Foam::definedInjector::nParcelsToInject
     return nParcels;
 }
 
-const Foam::vector Foam::definedInjector::position() const
+const Foam::vector Foam::definedInjector::position(const label n) const
 {
     return position_;
 }
 
 Foam::vector Foam::definedInjector::position
 (
+    const label n,
     const scalar time,
     const bool twoD,
     const scalar angleOfWedge,
@@ -227,12 +228,21 @@ Foam::vector Foam::definedInjector::position
     return position_;
 }
 
+Foam::label Foam::definedInjector::nHoles() const
+{
+    return 1;
+}
+
 Foam::scalar Foam::definedInjector::d() const
 {
     return d_;
 }
 
-const Foam::vector& Foam::definedInjector::direction() const
+const Foam::vector& Foam::definedInjector::direction
+(
+    const label i,
+    const scalar time
+) const
 {
     return direction_;
 }
@@ -344,4 +354,14 @@ void Foam::definedInjector::correctProfiles
     }
 }
 
+Foam::vector Foam::definedInjector::tan1(const label n) const
+{
+    return tangentialInjectionVector1_;
+}
+
+Foam::vector Foam::definedInjector::tan2(const label n) const
+{
+    return tangentialInjectionVector2_;
+}
+
 // ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.H b/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.H
index 6389b94339b791e2e188d120b590879aeda8867f..66ca317b3ae549075a5a4c79af350abba1b39f5f 100644
--- a/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.H
+++ b/src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.H
@@ -132,11 +132,12 @@ public:
         ) const;
 
         //- Return the injection position
-        const vector position() const;
+        const vector position(const label n) const;
 
         //- Return the injection position
         vector position
         (
+            const label n,
             const scalar time,
             const bool twoD,
             const scalar angleOfWedge,
@@ -146,11 +147,18 @@ public:
             Random& rndGen
         ) const;
     
+        //- Return the number of holes
+        label nHoles() const;
+
         //- Return the injector diameter
         scalar d() const;
 
         //- Return the injection direction
-        const vector& direction() const;
+        const vector& direction
+        (
+            const label i,
+            const scalar time
+        ) const;
 
         //- Return the mass of the injected particle
         scalar mass
@@ -210,6 +218,9 @@ public:
 
         scalar Cd(const scalar time) const;
 
+        vector tan1(const label n) const;
+        vector tan2(const label n) const;
+
         void correctProfiles
         (
             const liquidMixture& fuel,
diff --git a/src/lagrangian/dieselSpray/injector/injectorType/injectorType.H b/src/lagrangian/dieselSpray/injector/injectorType/injectorType.H
index 878cc8897d3b16ac233f1b466cfbbc9c411b351a..ab789f5213f69af10a244e1a825ac2281352af06 100644
--- a/src/lagrangian/dieselSpray/injector/injectorType/injectorType.H
+++ b/src/lagrangian/dieselSpray/injector/injectorType/injectorType.H
@@ -114,11 +114,12 @@ public:
         ) const = 0;
 
         //- Return the injection position
-        virtual const vector position() const = 0;
+        virtual const vector position(const label n) const = 0;
 
         //- Return the injection position
         virtual vector position
         (
+            const label n,
             const scalar time,
             const bool twoD,
             const scalar angleOfWedge,
@@ -127,12 +128,19 @@ public:
             const vector& axisOfWedgeNormal,
             Random& rndGen
         ) const = 0;
+
+        //- Return the number of holes
+        virtual label nHoles() const = 0;
     
         //- Return the injector diameter
         virtual scalar d() const = 0;
 
-        //- Return the injection direction
-        virtual const vector& direction() const = 0;
+        //- Return the injection direction for hole i
+        virtual const vector& direction
+        (
+            const label i, 
+            const scalar time
+        ) const = 0;
 
         //- Return the mass of the injected liquid between times
         virtual scalar mass
@@ -189,6 +197,12 @@ public:
 
         virtual bool pressureIndependentVelocity() const = 0;
 
+        //- Return a vector perpendicular to the injection direction and tan2 for hole n
+        virtual vector tan1(const label n) const = 0;
+
+        //- Return a vector perpendicular to the injection direction and tan1 for hole n
+        virtual vector tan2(const label n) const = 0;
+
         scalar getTableValue
         (
             const List<pair>& table,
diff --git a/src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.C b/src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.C
new file mode 100644
index 0000000000000000000000000000000000000000..54b4c4f4860bf68a6b6b4a0ea3a83fb05537f360
--- /dev/null
+++ b/src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.C
@@ -0,0 +1,420 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiHoleInjector.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Random.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+namespace Foam
+{
+
+defineTypeNameAndDebug(multiHoleInjector, 0);
+
+addToRunTimeSelectionTable
+(
+    injectorType,
+    multiHoleInjector,
+    dictionary
+);
+}
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::multiHoleInjector::multiHoleInjector
+(
+    const Foam::Time& t,
+    const Foam::dictionary& dict
+)
+:
+    injectorType(t, dict),
+    propsDict_(dict.subDict(typeName + "Props")),
+    centerPosition_(propsDict_.lookup("position")),
+    xyAngle_(readScalar(propsDict_.lookup("xyAngle"))),
+    zAngle_(readScalar(propsDict_.lookup("zAngle"))),
+    nHoles_(readLabel(propsDict_.lookup("nHoles"))),
+    umbrellaAngle_(readScalar(propsDict_.lookup("umbrellaAngle"))),
+    nozzleTipDiameter_(readScalar(propsDict_.lookup("nozzleTipDiameter"))),
+    angleSpacing_(propsDict_.lookup("angleSpacing")),
+    d_(readScalar(propsDict_.lookup("diameter"))),
+    Cd_(readScalar(propsDict_.lookup("Cd"))),
+    mass_(readScalar(propsDict_.lookup("mass"))),
+    nParcels_(readLabel(propsDict_.lookup("nParcels"))),
+    X_(propsDict_.lookup("X")),
+    massFlowRateProfile_(propsDict_.lookup("massFlowRateProfile")),
+    velocityProfile_(massFlowRateProfile_),
+    injectionPressureProfile_(massFlowRateProfile_),
+    CdProfile_(massFlowRateProfile_),
+    TProfile_(propsDict_.lookup("temperatureProfile")),
+    averageParcelMass_(nHoles_*mass_/nParcels_),
+    direction_(nHoles_),
+    position_(nHoles_),
+    pressureIndependentVelocity_(true),
+    tangentialInjectionVector1_(nHoles_),
+    tangentialInjectionVector2_(nHoles_)
+{
+
+
+    // check if time entries for soi and eoi match
+    if (mag(massFlowRateProfile_[0][0]-TProfile_[0][0]) > SMALL)
+    {
+        FatalError << "multiHoleInjector::multiHoleInjector(const time& t, const dictionary dict) " << endl
+            << " start-times do not match for TemperatureProfile and massFlowRateProfile."
+            << abort(FatalError);
+    }
+
+    if (mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]-TProfile_[TProfile_.size()-1][0]) > SMALL)
+    {
+        FatalError << "multiHoleInjector::multiHoleInjector(const time& t, const dictionary dict) " << endl
+            << " end-times do not match for TemperatureProfile and massFlowRateProfile."
+            << abort(FatalError);
+    }
+
+    // convert CA to real time
+    forAll(massFlowRateProfile_, i)
+    {
+        massFlowRateProfile_[i][0] = t.userTimeToTime(massFlowRateProfile_[i][0]);
+        velocityProfile_[i][0] = massFlowRateProfile_[i][0];
+        injectionPressureProfile_[i][0] = massFlowRateProfile_[i][0];
+    }
+    
+    forAll(TProfile_, i)
+    {
+        TProfile_[i][0] = t.userTimeToTime(TProfile_[i][0]);
+    }
+
+    scalar integratedMFR = integrateTable(massFlowRateProfile_);
+
+    forAll(massFlowRateProfile_, i)
+    {
+        // correct the massFlowRateProfile to match the injected mass
+        massFlowRateProfile_[i][1] *= mass_/integratedMFR;
+        
+        CdProfile_[i][0] = massFlowRateProfile_[i][0];
+        CdProfile_[i][1] = Cd_;
+    }
+
+    setTangentialVectors();
+
+    // check molar fractions
+    scalar Xsum = 0.0;
+    forAll(X_, i)
+    {
+        Xsum += X_[i];
+    }
+
+    if (mag(Xsum - 1.0) > SMALL)
+    {
+        Info << "Warning!!!\n multiHoleInjector::multiHoleInjector(const time& t, Istream& is)"
+            << "X does not add up to 1.0, correcting molar fractions."
+            << endl;
+        forAll(X_, i)
+        {
+            X_[i] /= Xsum;
+        }
+    }
+
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::multiHoleInjector::~multiHoleInjector()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::multiHoleInjector::setTangentialVectors()
+{
+    scalar pi180 = mathematicalConstant::pi/180.0;
+    scalar alpha = xyAngle_*pi180;
+    scalar phi = zAngle_*pi180;
+
+    vector xp(cos(alpha), sin(alpha), 0.0);
+    vector zp(cos(alpha)*sin(phi), sin(alpha)*sin(phi), cos(phi));
+    if (mag(zp-xp) < 1.0e-15)
+    {
+        xp = vector(0.0, 0.0, -1.0);
+        xp -= (xp & zp)*zp;
+        xp /= mag(xp);
+    }
+    vector yp = zp^xp;
+
+//    Info << "xp = " << xp << endl;
+//    Info << "yp = " << yp << endl;
+//    Info << "zp = " << zp << endl;
+
+    scalar angle = 0.0;
+    scalar u = umbrellaAngle_*pi180/2.0;
+    for(label i=0; i<nHoles_; i++)
+    {
+        angle += angleSpacing_[i];
+        scalar v = angle*pi180;
+        direction_[i] = cos(v)*sin(u)*xp + sin(v)*sin(u)*yp + cos(u)*zp;
+        vector dp = direction_[i] - (direction_[i] & zp)*direction_[i];
+        if (mag(dp) > SMALL)
+        {
+            dp /= mag(dp);
+        }
+        position_[i] = centerPosition_ + 0.5*nozzleTipDiameter_*dp;
+//        Info << "i = " << i << ", dir = " << direction_[i] << ", pos = " << position_[i] << endl;
+    }
+
+    Random rndGen(label(0));
+
+    for(label i=0; i<nHoles_; i++)
+    {
+        vector tangent(vector::zero);
+        scalar magV = 0;
+        while (magV < SMALL)
+        {
+            vector testThis = rndGen.vector01();
+            
+            tangent = testThis - (testThis & direction_[i])*direction_[i];
+            magV = mag(tangent);
+        }
+
+        tangentialInjectionVector1_[i] = tangent/magV;
+        tangentialInjectionVector2_[i] = direction_[i] ^ tangentialInjectionVector1_[i];
+
+    }   
+}
+
+
+Foam::label Foam::multiHoleInjector::nParcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+) const
+{
+
+    scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
+    label nParcels = label(mInj/averageParcelMass_ + 0.49);
+    
+    return nParcels;
+}
+
+const Foam::vector Foam::multiHoleInjector::position(const label n) const
+{
+    return position_[n];
+}
+
+Foam::vector Foam::multiHoleInjector::position
+(
+    const label n,
+    const scalar time,
+    const bool twoD,
+    const scalar angleOfWedge,
+    const vector& axisOfSymmetry,
+    const vector& axisOfWedge,
+    const vector& axisOfWedgeNormal,
+    Random& rndGen
+) const
+{
+    if (twoD)
+    {
+        scalar is = position_[n] & axisOfSymmetry;
+        scalar magInj = mag(position_[n] - is*axisOfSymmetry);
+
+        vector halfWedge =
+            axisOfWedge*cos(0.5*angleOfWedge)
+          + axisOfWedgeNormal*sin(0.5*angleOfWedge);
+        halfWedge /= mag(halfWedge);
+
+        return (is*axisOfSymmetry + magInj*halfWedge);
+    }
+    else
+    {
+        // otherwise, disc injection
+        scalar iRadius = d_*rndGen.scalar01();
+        scalar iAngle = 2.0*mathematicalConstant::pi*rndGen.scalar01();
+
+        return
+        ( 
+            position_[n]
+          + iRadius
+          * (
+              tangentialInjectionVector1_[n]*cos(iAngle)
+            + tangentialInjectionVector2_[n]*sin(iAngle)
+          )
+        );
+        
+    }
+    return position_[0];
+}
+
+Foam::label Foam::multiHoleInjector::nHoles() const
+{
+    return nHoles_;
+}
+
+Foam::scalar Foam::multiHoleInjector::d() const
+{
+    return d_;
+}
+
+const Foam::vector& Foam::multiHoleInjector::direction
+(
+    const label i,
+    const scalar time
+) const
+{
+    return direction_[i];
+}
+
+Foam::scalar Foam::multiHoleInjector::mass
+(
+    const scalar time0,
+    const scalar time1,
+    const bool twoD,
+    const scalar angleOfWedge
+) const
+{
+    scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
+
+    // correct mass if calculation is 2D 
+    if (twoD)
+    {
+        mInj *= 0.5*angleOfWedge/mathematicalConstant::pi;
+    }
+
+    return mInj;
+}
+
+Foam::scalar Foam::multiHoleInjector::mass() const
+{
+    return mass_;
+}
+
+const Foam::scalarField& Foam::multiHoleInjector::X() const
+{
+    return X_;
+}
+
+Foam::List<Foam::multiHoleInjector::pair> Foam::multiHoleInjector::T() const
+{
+    return TProfile_;
+}
+
+Foam::scalar Foam::multiHoleInjector::T(const scalar time) const
+{
+    return getTableValue(TProfile_, time);
+}
+
+Foam::scalar Foam::multiHoleInjector::tsoi() const
+{
+    return massFlowRateProfile_[0][0];
+}
+
+Foam::scalar Foam::multiHoleInjector::teoi() const
+{
+    return massFlowRateProfile_[massFlowRateProfile_.size()-1][0];
+}
+
+Foam::scalar Foam::multiHoleInjector::massFlowRate
+(
+    const scalar time
+) const
+{
+    return getTableValue(massFlowRateProfile_, time);
+}
+
+Foam::scalar Foam::multiHoleInjector::injectionPressure
+(
+    const scalar time
+) const
+{
+    return getTableValue(injectionPressureProfile_, time);
+}
+
+Foam::scalar Foam::multiHoleInjector::velocity
+(
+    const scalar time
+) const
+{
+    return getTableValue(velocityProfile_, time);
+}
+
+Foam::List<Foam::multiHoleInjector::pair> Foam::multiHoleInjector::CdProfile() const
+{
+    return CdProfile_;
+}
+
+Foam::scalar Foam::multiHoleInjector::Cd
+(
+    const scalar time
+) const
+{
+    return Cd_;
+}
+
+Foam::scalar Foam::multiHoleInjector::fractionOfInjection(const scalar time) const
+{
+    return integrateTable(massFlowRateProfile_, time)/mass_;
+}
+
+Foam::scalar Foam::multiHoleInjector::injectedMass
+(
+    const scalar t
+) const
+{
+    return mass_*fractionOfInjection(t);
+}
+
+
+void Foam::multiHoleInjector::correctProfiles
+(
+    const liquidMixture& fuel,
+    const scalar referencePressure
+)
+{
+
+    scalar A = nHoles_*0.25*mathematicalConstant::pi*pow(d_, 2.0);
+
+    forAll(velocityProfile_, i)
+    {
+        scalar time = velocityProfile_[i][0];
+        scalar rho = fuel.rho(referencePressure, T(time), X_);
+        scalar v = massFlowRateProfile_[i][1]/(Cd_*rho*A);
+        velocityProfile_[i][1] = v;
+        injectionPressureProfile_[i][1] = referencePressure + 0.5*rho*v*v;
+    }
+}
+
+Foam::vector Foam::multiHoleInjector::tan1(const label n) const
+{
+    return tangentialInjectionVector1_[n];
+}
+
+Foam::vector Foam::multiHoleInjector::tan2(const label n) const
+{
+    return tangentialInjectionVector2_[n];
+}
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.H b/src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.H
new file mode 100644
index 0000000000000000000000000000000000000000..30233d72a409d9c37aa31488f6e192419fb0978b
--- /dev/null
+++ b/src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.H
@@ -0,0 +1,248 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::multiHoleInjector
+
+Description
+    The unit injector
+
+SourceFiles
+    multiHoleInjectorI.H
+    multiHoleInjector.C
+    multiHoleInjectorIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiHoleInjector_H
+#define multiHoleInjector_H
+
+#include "injectorType.H"
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class multiHoleInjector Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiHoleInjector
+:
+    public injectorType
+{
+
+private:
+
+    typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
+
+    // Private data
+
+        dictionary propsDict_;
+
+        vector centerPosition_;
+        scalar xyAngle_;
+        scalar zAngle_;
+        label nHoles_;
+        scalar umbrellaAngle_;
+        scalar nozzleTipDiameter_;
+        List<scalar> angleSpacing_;
+        scalar d_;
+        scalar Cd_;
+        scalar mass_;
+        label nParcels_;
+        scalarField X_;
+        List<pair> massFlowRateProfile_;
+        List<pair> velocityProfile_;
+        List<pair> injectionPressureProfile_;
+        List<pair> CdProfile_;
+        List<pair> TProfile_;
+        scalar averageParcelMass_;
+        List<vector> direction_;
+        List<vector> position_;
+
+        bool pressureIndependentVelocity_;
+
+        //- two orthogonal vectors that are also orthogonal
+        //  to the injection direction
+        List<vector> tangentialInjectionVector1_, tangentialInjectionVector2_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        multiHoleInjector(const multiHoleInjector&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const multiHoleInjector&);
+
+        //- Create two vectors orthonoal to each other
+        //  and the injection vector
+        void setTangentialVectors();
+
+        //- Return the fraction of the total injected liquid
+        scalar fractionOfInjection(const scalar time) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiHoleInjector");
+
+
+    // Constructors
+
+        //- Construct from components
+        multiHoleInjector
+        (
+            const Time& t,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        ~multiHoleInjector();
+
+
+    // Member Functions
+
+        //- Return number of particles to inject
+        label nParcelsToInject
+        (
+            const scalar t0,
+            const scalar t1
+        ) const;
+
+        //- Return the injection position
+        const vector position(const label n) const;
+
+        //- Return the injection position
+        vector position
+        (
+            const label n,
+            const scalar time,
+            const bool twoD,
+            const scalar angleOfWedge,
+            const vector& axisOfSymmetry,
+            const vector& axisOfWedge,
+            const vector& axisOfWedgeNormal,
+            Random& rndGen
+        ) const;
+    
+        //- Return the number of holes
+        label nHoles() const;
+
+        //- Return the injector diameter
+        scalar d() const;
+
+        //- Return the injection direction
+        const vector& direction
+        (
+            const label i,
+            const scalar time
+        ) const;
+
+        //- Return the mass of the injected particle
+        scalar mass
+        (
+            const scalar t0,
+            const scalar t1,
+            const bool twoD,
+            const scalar angleOfWedge
+        ) const;
+
+        //- Return the mass injected by the injector
+        scalar mass() const;
+
+        //- Return the fuel mass fractions of the injected particle
+        const scalarField& X() const;
+
+        //- Return the temperature profile of the injected particle
+        List<pair> T() const;
+
+        //- Return the temperature of the injected particle
+        scalar T(const scalar time) const;
+
+        //- Return the start-of-injection time
+        scalar tsoi() const;
+
+        //- Return the end-of-injection time
+        scalar teoi() const;
+
+        //- Return the injected liquid mass
+        scalar injectedMass(const scalar t) const;
+
+        List<pair> massFlowRateProfile() const
+        {
+            return massFlowRateProfile_;
+        }
+
+        scalar massFlowRate(const scalar time) const;
+
+        List<pair> injectionPressureProfile() const
+        {
+            return injectionPressureProfile_;
+        }
+    
+        scalar injectionPressure(const scalar time) const;
+
+        List<pair> velocityProfile() const
+        {
+            return velocityProfile_;
+        }
+        
+        scalar velocity(const scalar time) const;
+
+        List<pair> CdProfile() const;
+        scalar Cd(const scalar time) const;
+
+        vector tan1(const label n) const;
+        vector tan2(const label n) const;
+
+        void correctProfiles
+        (
+            const liquidMixture& fuel,
+            const scalar referencePressure
+        );
+
+        bool pressureIndependentVelocity() const
+        {
+            return pressureIndependentVelocity_;
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.C b/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.C
index ef3196695d18242f594af429beeb16eed262e6ba..192319f8ced148c3004fceb550623053e638705b 100644
--- a/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.C
+++ b/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.C
@@ -187,13 +187,14 @@ Foam::label Foam::swirlInjector::nParcelsToInject
     return nParcels;
 }
 
-const Foam::vector Foam::swirlInjector::position() const
+const Foam::vector Foam::swirlInjector::position(const label n) const
 {
     return position_;
 }
 
 Foam::vector Foam::swirlInjector::position
 (
+    const label n,
     const scalar time,
     const bool twoD,
     const scalar angleOfWedge,
@@ -236,12 +237,21 @@ Foam::vector Foam::swirlInjector::position
     return position_;
 }
 
+Foam::label Foam::swirlInjector::nHoles() const
+{
+    return 1;
+}
+
 Foam::scalar Foam::swirlInjector::d() const
 {
     return d_;
 }
 
-const Foam::vector& Foam::swirlInjector::direction() const
+const Foam::vector& Foam::swirlInjector::direction
+(
+    const label i,
+    const scalar time
+) const
 {
     return direction_;
 }
@@ -370,4 +380,14 @@ void Foam::swirlInjector::correctProfiles
     }
 }
 
+Foam::vector Foam::swirlInjector::tan1(const label n) const
+{
+    return tangentialInjectionVector1_;
+}
+
+Foam::vector Foam::swirlInjector::tan2(const label n) const
+{
+    return tangentialInjectionVector2_;
+}
+
 // ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.H b/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.H
index 8c4853b9b7b257dd8c37c2b21a5ad5e46aa5efcd..4cb6cf1da56fd1a1f56ce256b77a1c29b4d0606c 100644
--- a/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.H
+++ b/src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.H
@@ -146,11 +146,12 @@ public:
         ) const;
 
         //- Return the injection position
-        const vector position() const;
+        const vector position(const label n) const;
 
         //- Return the injection position
         vector position
         (
+            const label n,
             const scalar time,
             const bool twoD,
             const scalar angleOfWedge,
@@ -160,11 +161,18 @@ public:
             Random& rndGen
         ) const;
     
+        //- Return the number of holes
+        label nHoles() const;
+
         //- Return the injector diameter
         scalar d() const;
 
         //- Return the injection direction
-        const vector& direction() const;
+        const vector& direction
+        (
+            const label i,
+            const scalar time
+        ) const;
 
         //- Return the mass of the injected particle
         scalar mass
@@ -220,6 +228,9 @@ public:
         //- Return the injected liquid mass
         scalar injectedMass(const scalar t) const;
 
+        vector tan1(const label n) const;
+        vector tan2(const label n) const;
+
         void correctProfiles
         (
             const liquidMixture& fuel,
diff --git a/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C b/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C
index 0429537b2f94e4064638b406163c7d5ca4fca92f..9363379cabaed45dec6768991e50686a6b5dee4f 100644
--- a/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C
+++ b/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C
@@ -60,17 +60,32 @@ Foam::unitInjector::unitInjector
     d_(readScalar(propsDict_.lookup("diameter"))),
     Cd_(readScalar(propsDict_.lookup("Cd"))),
     mass_(readScalar(propsDict_.lookup("mass"))),
-    T_(readScalar(propsDict_.lookup("temperature"))),
     nParcels_(readLabel(propsDict_.lookup("nParcels"))),
     X_(propsDict_.lookup("X")),
     massFlowRateProfile_(propsDict_.lookup("massFlowRateProfile")),
     velocityProfile_(massFlowRateProfile_),
     injectionPressureProfile_(massFlowRateProfile_),
     CdProfile_(massFlowRateProfile_),
-    TProfile_(massFlowRateProfile_),
+    TProfile_(propsDict_.lookup("temperatureProfile")),
     averageParcelMass_(mass_/nParcels_),
     pressureIndependentVelocity_(true)
 {
+
+    // check if time entries for soi and eoi match
+    if (mag(massFlowRateProfile_[0][0]-TProfile_[0][0]) > SMALL)
+    {
+        FatalError << "unitInjector::unitInjector(const time& t, const dictionary dict) " << endl
+            << " start-times do not match for TemperatureProfile and massFlowRateProfile."
+            << abort(FatalError);
+    }
+
+    if (mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]-TProfile_[TProfile_.size()-1][0]) > SMALL)
+    {
+        FatalError << "unitInjector::unitInjector(const time& t, const dictionary dict) " << endl
+            << " end-times do not match for TemperatureProfile and massFlowRateProfile."
+            << abort(FatalError);
+    }
+
     // convert CA to real time
     forAll(massFlowRateProfile_, i)
     {
@@ -79,6 +94,11 @@ Foam::unitInjector::unitInjector
         injectionPressureProfile_[i][0] = massFlowRateProfile_[i][0];
     }
     
+    forAll(TProfile_, i)
+    {
+        TProfile_[i][0] = t.userTimeToTime(TProfile_[i][0]);
+    }
+
     scalar integratedMFR = integrateTable(massFlowRateProfile_);
 
     forAll(massFlowRateProfile_, i)
@@ -86,9 +106,6 @@ Foam::unitInjector::unitInjector
         // correct the massFlowRateProfile to match the injected mass
         massFlowRateProfile_[i][1] *= mass_/integratedMFR;
         
-        TProfile_[i][0] = massFlowRateProfile_[i][0];
-        TProfile_[i][1] = T_;
-
         CdProfile_[i][0] = massFlowRateProfile_[i][0];
         CdProfile_[i][1] = Cd_;
     }
@@ -159,13 +176,14 @@ Foam::label Foam::unitInjector::nParcelsToInject
     return nParcels;
 }
 
-const Foam::vector Foam::unitInjector::position() const
+const Foam::vector Foam::unitInjector::position(const label n) const
 {
     return position_;
 }
 
 Foam::vector Foam::unitInjector::position
 (
+    const label n,
     const scalar time,
     const bool twoD,
     const scalar angleOfWedge,
@@ -208,12 +226,21 @@ Foam::vector Foam::unitInjector::position
     return position_;
 }
 
+Foam::label Foam::unitInjector::nHoles() const
+{
+    return 1;
+}
+
 Foam::scalar Foam::unitInjector::d() const
 {
     return d_;
 }
 
-const Foam::vector& Foam::unitInjector::direction() const
+const Foam::vector& Foam::unitInjector::direction
+(
+    const label i,
+    const scalar time
+) const
 {
     return direction_;
 }
@@ -254,7 +281,7 @@ Foam::List<Foam::unitInjector::pair> Foam::unitInjector::T() const
 
 Foam::scalar Foam::unitInjector::T(const scalar time) const
 {
-    return T_;
+    return getTableValue(TProfile_, time);
 }
 
 Foam::scalar Foam::unitInjector::tsoi() const
@@ -328,14 +355,24 @@ void Foam::unitInjector::correctProfiles
     scalar A = 0.25*mathematicalConstant::pi*pow(d_, 2.0);
     scalar pDummy = 1.0e+5;
 
-    scalar rho = fuel.rho(pDummy, T_, X_);
-
     forAll(velocityProfile_, i)
     {
+        scalar time = velocityProfile_[i][0];
+        scalar rho = fuel.rho(pDummy, T(time), X_);
         scalar v = massFlowRateProfile_[i][1]/(Cd_*rho*A);
         velocityProfile_[i][1] = v;
         injectionPressureProfile_[i][1] = referencePressure + 0.5*rho*v*v;
     }
 }
 
+Foam::vector Foam::unitInjector::tan1(const label n) const
+{
+    return tangentialInjectionVector1_;
+}
+
+Foam::vector Foam::unitInjector::tan2(const label n) const
+{
+    return tangentialInjectionVector2_;
+}
+
 // ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.H b/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.H
index 057700aad783026e317c9869c9b2d001c4f393c0..c7592417437a0b4c1d16e17b49d2f638081397a7 100644
--- a/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.H
+++ b/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.H
@@ -68,7 +68,6 @@ private:
         scalar d_;
         scalar Cd_;
         scalar mass_;
-        scalar T_;
         label nParcels_;
         scalarField X_;
         List<pair> massFlowRateProfile_;
@@ -132,11 +131,12 @@ public:
         ) const;
 
         //- Return the injection position
-        const vector position() const;
+        const vector position(const label n) const;
 
         //- Return the injection position
         vector position
         (
+            const label n,
             const scalar time,
             const bool twoD,
             const scalar angleOfWedge,
@@ -146,11 +146,18 @@ public:
             Random& rndGen
         ) const;
     
+        //- Return the number of holes
+        label nHoles() const;
+
         //- Return the injector diameter
         scalar d() const;
 
         //- Return the injection direction
-        const vector& direction() const;
+        const vector& direction
+        (
+            const label i,
+            const scalar time
+        ) const;
 
         //- Return the mass of the injected particle
         scalar mass
@@ -206,6 +213,9 @@ public:
         List<pair> CdProfile() const;
         scalar Cd(const scalar time) const;
 
+        vector tan1(const label n) const;
+        vector tan2(const label n) const;
+
         void correctProfiles
         (
             const liquidMixture& fuel,
diff --git a/src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C b/src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C
index 1737b9c20db7878ea64f0a027a50aa1d1c3c9019..07c6362d16f40e2484eb2f5db7c268ed9fc2728c 100644
--- a/src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C
+++ b/src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C
@@ -86,6 +86,7 @@ void parcel::setRelaxationTimes
     W = 1.0/W;
 
     scalarField Xf(Nf, 0.0);
+    scalarField Yf(Nf, 0.0);
     scalarField psat(Nf, 0.0);
     scalarField msat(Nf, 0.0);
 
@@ -94,6 +95,7 @@ void parcel::setRelaxationTimes
         label j = sDB.liquidToGasIndex()[i];
         scalar Y = sDB.composition().Y()[j][celli];        
         scalar Wi = sDB.gasProperties()[j].W();
+        Yf[i] = Y;
         Xf[i] = Y*W/Wi;
         psat[i] = fuels.properties()[i].pv(pressure, temperature);
         msat[i] = min(1.0, psat[i]/pressure)*Wi/W;
@@ -116,6 +118,22 @@ void parcel::setRelaxationTimes
     scalar rhoFuelVap = pressureAtSurface*fuels.W(X())/(specie::RR*Tf);
 
     scalarField Xs(sDB.fuels().Xs(pressure, temperature, T(), Xf, X()));
+    scalarField Ys(Nf, 0.0);
+    scalar Wliq = 0.0;
+
+    for(label i=0; i<Nf; i++)
+    {
+        label j = sDB.liquidToGasIndex()[i];
+        scalar Wi = sDB.gasProperties()[j].W();
+        Wliq += Xs[i]*Wi;
+    }
+
+    for(label i=0; i<Nf; i++)
+    {
+        label j = sDB.liquidToGasIndex()[i];
+        scalar Wi = sDB.gasProperties()[j].W();
+        Ys[i] = Xs[i]*Wi/Wliq;
+    }
 
     scalar Reynolds = Re(Up, nuf);
     scalar Prandtl = Pr(cpMixture, muf, kMixture);
@@ -238,7 +256,6 @@ void parcel::setRelaxationTimes
                     dp0 = dp;
                 }
                 
-                label j = sDB.liquidToGasIndex()[i];
                 scalar vapourSurfaceEnthalpy = 0.0;
                 scalar vapourFarEnthalpy = 0.0;
                 
@@ -266,8 +283,6 @@ void parcel::setRelaxationTimes
                     vapourSurfaceEnthalpy,
                     vapourFarEnthalpy,
                     cpMixture,
-                    Xs[i],
-                    Xf[j],
                     temperature,
                     kLiquid
                 );
diff --git a/src/lagrangian/dieselSpray/spray/findInjectorCell.H b/src/lagrangian/dieselSpray/spray/findInjectorCell.H
index c68f7deaf88ddb116dd264f5195ecba2fbd60d29..f1a8d7fc152c89feac4149d66366a0145ec9802e 100644
--- a/src/lagrangian/dieselSpray/spray/findInjectorCell.H
+++ b/src/lagrangian/dieselSpray/spray/findInjectorCell.H
@@ -22,7 +22,7 @@ reduce(foundCell, orOp<bool>());
 
 if (!foundCell)
 {
-    injectionPosition = it.position();
+    injectionPosition = it->position(n);
     injectorCell = mesh_.findCell(injectionPosition);
     
     if (injectorCell >= 0)
diff --git a/src/lagrangian/dieselSpray/spray/spray.C b/src/lagrangian/dieselSpray/spray/spray.C
index 4aa5c775efe57d18b891eec44d5259747f836a5d..fb9e16bc9b7ddc9f8c1588973946828cbb0fb2c5 100644
--- a/src/lagrangian/dieselSpray/spray/spray.C
+++ b/src/lagrangian/dieselSpray/spray/spray.C
@@ -52,6 +52,7 @@ defineTemplateTypeNameAndDebug(IOPtrList<injector>, 0);
 // Construct from components
 Foam::spray::spray
 (
+    const volPointInterpolation& vpi,
     const volVectorField& U,
     const volScalarField& rho,
     const volScalarField& p,
@@ -66,6 +67,7 @@ Foam::spray::spray
     runTime_(U.time()),
     time0_(runTime_.value()),
     mesh_(U.mesh()),
+    volPointInterpolation_(vpi),
     rndGen_(label(0)),
 
     U_(U),
@@ -261,7 +263,7 @@ Foam::spray::spray
         {
             FatalErrorIn
             (
-                "spray::spray(const volVectorField& U, "
+                "spray::spray(const pointMesh& pMesh, const volVectorField& U, "
                 "const volScalarField& rho, const volScalarField& p, "
                 "const volScalarField& T, const combustionMixture& composition,"
                 "const PtrList<specieProperties>& gaseousFuelProperties, "
diff --git a/src/lagrangian/dieselSpray/spray/spray.H b/src/lagrangian/dieselSpray/spray/spray.H
index 603d0c8b68a4abaffe1bb8089bce418698d2bb30..83c93c16f0d65a2a84c075da5fa6110787a905dd 100644
--- a/src/lagrangian/dieselSpray/spray/spray.H
+++ b/src/lagrangian/dieselSpray/spray/spray.H
@@ -36,6 +36,7 @@ Description
 #include "parcel.H"
 #include "injector.H"
 #include "IOPtrList.H"
+#include "volPointInterpolation.H"
 #include "interpolation.H"
 #include "liquid.H"
 #include "sprayThermoTypes.H"
@@ -75,6 +76,7 @@ class spray
             const Time& runTime_;
             scalar time0_;
             const fvMesh& mesh_;
+            const volPointInterpolation& volPointInterpolation_;
 
             //- Random number generator
             Random rndGen_;
@@ -187,6 +189,7 @@ public:
         //- Construct from components
         spray
         (
+            const volPointInterpolation& vpi,
             const volVectorField& U,
             const volScalarField& rho,
             const volScalarField& p,
diff --git a/src/lagrangian/dieselSpray/spray/sprayFunctions.C b/src/lagrangian/dieselSpray/spray/sprayFunctions.C
index 2a747d52f8cb26df973edaaae9972ad32cecd9b0..5159a8179c579d21be54bc9da07cbc5352808ed8 100644
--- a/src/lagrangian/dieselSpray/spray/sprayFunctions.C
+++ b/src/lagrangian/dieselSpray/spray/sprayFunctions.C
@@ -246,7 +246,23 @@ scalar spray::liquidPenetration
     const scalar prc
 ) const
 {
-    vector ip = injectors_[nozzlei].properties()->position();
+
+    label nHoles = injectors_[nozzlei].properties()->nHoles();
+    vector ip(vector::zero);
+    if (nHoles > 1)
+    {
+        for(label i=0;i<nHoles;i++)
+        {
+            ip += injectors_[nozzlei].properties()->position(i);
+        }
+        ip /= nHoles;
+    }
+    else
+    {
+        ip = injectors_[nozzlei].properties()->position(0);
+    }
+
+//    vector ip = injectors_[nozzlei].properties()->position();
     scalar d = 0.0;
     scalar mTot = 0.0;
 
diff --git a/src/lagrangian/dieselSpray/spray/sprayInject.C b/src/lagrangian/dieselSpray/spray/sprayInject.C
index 7b3a5d94c5d58c9b5adb188c15312d368dd58f13..a42acebf58e49297ca457755a180a27d78576e97 100644
--- a/src/lagrangian/dieselSpray/spray/sprayInject.C
+++ b/src/lagrangian/dieselSpray/spray/sprayInject.C
@@ -45,24 +45,31 @@ void spray::inject()
     // Inject the parcels for each injector sequentially
     forAll(injectors_, i)
     {
-        const injectorType& it = injectors_[i].properties();
+        autoPtr<injectorType>& it = injectors()[i].properties();
+        if (!it->pressureIndependentVelocity())
+        {
+            scalar referencePressure = p().average().value();
+            it->correctProfiles(fuels(), referencePressure);
+        }
+
+        const label nHoles = it->nHoles();
 
         // parcels have the same mass during a timestep
-        scalar mass = it.mass(time0, time, twoD_, angleOfWedge_);
+        scalar mass = it->mass(time0, time, twoD_, angleOfWedge_);
 
-        label Np = it.nParcelsToInject(time0, time);
+        label Np = it->nParcelsToInject(time0, time);
 
         if (mass > 0)
         {
             Np = max(1, Np);
-            scalar mp = mass/Np;
+            scalar mp = mass/Np/nHoles;
         
             // constT is only larger than zero for the first 
             // part of the injection
             scalar constT = max
             (
                 0.0,
-                it.tsoi() - time0
+                it->tsoi() - time0
             );
 
             // deltaT is the duration of injection during this timestep
@@ -71,8 +78,8 @@ void spray::inject()
                 runTime_.deltaT().value(),
                 min
                 (
-                    time - it.tsoi(),
-                    it.teoi() - time0
+                    time - it->tsoi(),
+                    it->teoi() - time0
                 )
             );
 
@@ -81,87 +88,92 @@ void spray::inject()
                 // calculate the time of injection for parcel 'j'
                 scalar toi = time0 + constT + deltaT*j/scalar(Np);
 
-                // calculate the velocity of the injected parcel
-                vector injectionPosition = it.position
-                (
-                    toi,
-                    twoD_,
-                    angleOfWedge_,
-                    axisOfSymmetry_,
-                    axisOfWedge_,
-                    axisOfWedgeNormal_,
-                    rndGen_
-                );
-                
-                scalar diameter = injection().d0(i, toi);
-                vector direction = injection().direction(i, toi, diameter);
-                vector U = injection().velocity(i, toi)*direction;
-
-                scalar symComponent = direction & axisOfSymmetry_;
-                vector normal = direction - symComponent*axisOfSymmetry_;
-                normal /= mag(normal);
-
-                // should be set from dict or model
-                scalar deviation = breakup().y0();
-                scalar ddev = breakup().yDot0();
-
-                label injectorCell = mesh_.findCell(injectionPosition);
-                
-#               include "findInjectorCell.H"
-
-                if (injectorCell >= 0)
+                for(label n=0; n<nHoles; n++)
                 {
-                    scalar liquidCore = 1.0;
-                
-                    // construct the parcel that is to be injected
 
-                    parcel* pPtr = new parcel
+                    // calculate the velocity of the injected parcel
+                    vector injectionPosition = it->position
                     (
-                        *this,
-                        injectionPosition,
-                        injectorCell,
-                        normal,
-                        diameter,
-                        it.T(toi),
-                        mp,
-                        deviation,
-                        ddev,
-                        0.0,
-                        0.0,
-                        0.0,
-                        liquidCore,
-                        scalar(i),
-                        U,
-                        vector::zero,
-                        it.X(),
-                        fuels_->components()
+                        n,
+                        toi,
+                        twoD_,
+                        angleOfWedge_,
+                        axisOfSymmetry_,
+                        axisOfWedge_,
+                        axisOfWedgeNormal_,
+                        rndGen_
                     );
+                
+                    scalar diameter = injection().d0(i, toi);
+                    vector direction = injection().direction(i, n, toi, diameter);
+                    vector U = injection().velocity(i, toi)*direction;
 
-                    injectedLiquidKE_ += 0.5*pPtr->m()*pow(mag(U), 2.0);
-                    
-                    scalar dt = time - toi;
-
-                    pPtr->stepFraction() =
-                        (runTime_.deltaT().value() - dt)
-                       /runTime_.deltaT().value();
+                    scalar symComponent = direction & axisOfSymmetry_;
+                    vector normal = direction - symComponent*axisOfSymmetry_;
+                    normal /= mag(normal);
 
-                    bool keepParcel = pPtr->move
-                    (
-                        *this
-                    );
+                    // should be set from dict or model
+                    scalar deviation = breakup().y0();
+                    scalar ddev = breakup().yDot0();
 
-                    if (keepParcel)
-                    {
-                        addParticle(pPtr);
-                    }
-                    else
+                    label injectorCell = mesh_.findCell(injectionPosition);
+                
+#                   include "findInjectorCell.H"
+   
+                    if (injectorCell >= 0)
                     {
-                        delete pPtr;
-                    }
-                }
-            }
-        }
-    }
+                        scalar liquidCore = 1.0;
+                
+                        // construct the parcel that is to be injected
+
+                        parcel* pPtr = new parcel
+                        (
+                            *this,
+                            injectionPosition,
+                            injectorCell,
+                            normal,
+                            diameter,
+                            it->T(toi),
+                            mp,
+                            deviation,
+                            ddev,
+                            0.0,
+                            0.0,
+                            0.0,
+                            liquidCore,
+                            scalar(i),
+                            U,
+                            vector::zero,
+                            it->X(),
+                            fuels_->components()
+                        );
+
+                        injectedLiquidKE_ += 0.5*pPtr->m()*pow(mag(U), 2.0);
+                    
+                        scalar dt = time - toi;
+
+                        pPtr->stepFraction() =
+                            (runTime_.deltaT().value() - dt)
+                           /runTime_.deltaT().value();
+
+                        bool keepParcel = pPtr->move
+                        (
+                            *this
+                        );
+  
+                        if (keepParcel)
+                        {
+                            addParticle(pPtr);
+                        }
+                        else
+                        {
+                            delete pPtr;
+                        }
+                    } // if (injectorCell....
+                } // for(label n=0...
+            } // for(label j=0....
+        } // if (mass>0)...
+    } // forAll(injectors)...
 
     time0_ = time;
 }
diff --git a/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/LISA/LISA.C b/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/LISA/LISA.C
index ffdeaf4e9ca6f1219446ed1de1ff9addaf1fe92a..47f15462ac6942be43562f5d12e1674c52a16cbd 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/LISA/LISA.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/LISA/LISA.C
@@ -140,12 +140,18 @@ void LISA::atomizeParcel
     const injectorType& it = 
         spray_.injectors()[label(p.injector())].properties();
 
-    const vector itPosition = it.position();
+    if (it.nHoles() > 1)
+    {
+        Info << "Warning: This atomization model is not suitable for multihole injector." << endl
+             << "Only the first hole will be used." << endl;
+    }
+
+    const vector itPosition = it.position(0);
     scalar pWalk = mag(p.position() - itPosition);
 
 //  Updating liquid sheet tickness... that is the droplet diameter
-
-    const vector direction = it.direction();
+ 
+    const vector direction = it.direction(0, spray_.runTime().value());
     
     scalar h = (p.position() - itPosition) & direction;
 
diff --git a/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/blobsSheetAtomization/blobsSheetAtomization.C b/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/blobsSheetAtomization/blobsSheetAtomization.C
index 1cb8db82ef48c8ab02696eb8a4256f9c0880d4d0..dc7f2b51f731f4f1a61140931a051f73edc68d2c 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/blobsSheetAtomization/blobsSheetAtomization.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/atomizationModel/blobsSheetAtomization/blobsSheetAtomization.C
@@ -112,7 +112,22 @@ void blobsSheetAtomization::atomizeParcel
     const injectorType& it =
         spray_.injectors()[label(p.injector())].properties();
 
-    const vector itPosition = it.position();
+    vector itPosition(vector::zero);
+    label nHoles = it.nHoles();
+    if (nHoles > 1)
+    {
+        for(label i=0; i<nHoles;i++)
+        {
+            itPosition += it.position(i);
+        }
+        itPosition /= nHoles;
+    }
+    else
+    {
+        itPosition = it.position(0);
+    }
+//    const vector itPosition = it.position();
+
 
     scalar lBU = B_ * sqrt
     (
diff --git a/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C b/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C
index 32e0ee4f49dbadf70bd9bfd69724b260d6b18504..8b4375234c03777f8bfd786ff9fe7102984e5620 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C
@@ -169,7 +169,7 @@ void reitzKHRT::breakupParcel
     // check if we have RT breakup
     if ((p.ct() > tauRT) && (lambdaRT < p.d()))
     {
-        // the RT breakup creates diameter/lmbdaRT new droplets
+        // the RT breakup creates diameter/lambdaRT new droplets
         p.ct() = -GREAT;
         scalar multiplier = p.d()/lambdaRT;
         scalar nDrops = multiplier*Np;
@@ -200,8 +200,6 @@ void reitzKHRT::breakupParcel
 
             scalar averageParcelMass = spray_.injectors()[injector].properties()->mass()/nParcels;
 
-            // NN. Since the parcel doesn't know from which injector
-            // it comes we use the first one to obtain a 'reference' mass
             if
             (
                 (p.ms()/averageParcelMass > msLimit_)
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.C b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.C
index fa84ee812410ab16fc34019595d5d9378da27755..4643c14c8edfac5403f0ae249e7ebefd5b3eb50c 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.C
@@ -196,8 +196,6 @@ scalar RutlandFlashBoil::boilingTime
     const scalar vapourSurfaceEnthalpy,
     const scalar vapourFarEnthalpy,
     const scalar cpGas,
-    const scalar Xs,
-    const scalar Xf,    
     const scalar temperature,
     const scalar kLiq
 ) const
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.H b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.H
index 010899df2e557e3c466ac2091fb3d84cccdce29d..6b1589ce5a554322694a85bc4b8f43e23d9ab4a5 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/RutlandFlashBoil/RutlandFlashBoil.H
@@ -135,8 +135,6 @@ public:
             const scalar vapourSurfaceEnthalpy,
             const scalar vapourFarEnthalpy,
             const scalar cpGas,
-            const scalar Xs,
-            const scalar Xf,
             const scalar temperature,
             const scalar kLiquid
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/evaporationModel/evaporationModel.H b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/evaporationModel/evaporationModel.H
index a2493609d4819fc417745fe6facd9abe24b26382..e336bc574b945b3b5b0e569ac4f499654b0ac503 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/evaporationModel/evaporationModel.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/evaporationModel/evaporationModel.H
@@ -143,8 +143,6 @@ public:
             const scalar vapourSurfaceEnthalpy,
             const scalar vapourFarEnthalpy,
             const scalar cpGas,
-            const scalar Xs,
-            const scalar Xf,
             const scalar temperature,
             const scalar kLiq
         ) const = 0;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.C b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.C
index ce506315ce1277448622aa217a549dd2c7e94e48..558ddda1594167316783008e90ee637f320d5f58 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.C
@@ -114,8 +114,6 @@ scalar noEvaporation::boilingTime
     const scalar,
     const scalar,
     const scalar,
-    const scalar,
-    const scalar,
     const scalar
 ) const
 {
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.H b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.H
index a17dd809632cfd360819392ecb281b3bbf0ce9fd..b4bdc272840072142cd4cba70c675847dfef9d9b 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/noEvaporation/noEvaporation.H
@@ -111,8 +111,6 @@ public:
             const scalar vapourSurfaceEnthalpy,
             const scalar vapourFarEnthalpy,
             const scalar cpGas,
-            const scalar Xs,
-            const scalar Xf,
             const scalar temperature,
             const scalar kLiq
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.C b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.C
index 944dd4df5ab8468f68b9e6aa52edf296d2fec55d..2c27286ca6d222dc71284106c83193c35980bd46 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.C
@@ -117,8 +117,6 @@ scalar saturateEvaporationModel::boilingTime
     const scalar, 
     const scalar, 
     const scalar, 
-    const scalar, 
-    const scalar, 
     const scalar 
 ) const
 {
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.H b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.H
index 1c02d95750a682844e0e6360e411b06ff292f880..51dfdf98acee248b654a5e8456ff0fbe281f6a17 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/saturateEvaporationModel/saturateEvaporationModel.H
@@ -123,8 +123,6 @@ public:
             const scalar, 
             const scalar, 
             const scalar, 
-            const scalar, 
-            const scalar, 
             const scalar 
         ) const;
 
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.C b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.C
index e1d39fceb6a3005ee8dcb19064bc6781efd9da7f..832dedc94b254293dff0f65979fd2f1544294206 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.C
@@ -181,8 +181,6 @@ scalar standardEvaporationModel::boilingTime
     const scalar, 
     const scalar, 
     const scalar, 
-    const scalar, 
-    const scalar, 
     const scalar 
 ) const
 {
diff --git a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.H b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.H
index 8527c5607d0f1a09b39da78f624662bc4986d505..a7c573b49811dddd29ef9745b0a2ab49b3339899 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/evaporationModel/standardEvaporationModel/standardEvaporationModel.H
@@ -127,8 +127,6 @@ public:
             const scalar, 
             const scalar, 
             const scalar, 
-            const scalar, 
-            const scalar, 
             const scalar 
         ) const;
 
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.C b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.C
index 71112ffa76fcd838fc38fb6f5fed2bf78cd9b190..f6740ec4fb7eaccc532625f6d49db9634cff3955 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.C
@@ -64,9 +64,7 @@ ChomiakInjector::ChomiakInjector
             sm.rndGen()
         )
     ),
-    maxSprayAngle_(ChomiakDict_.lookup("maxSprayConeAngle")),
-    tan1_(maxSprayAngle_.size()),
-    tan2_(maxSprayAngle_.size())
+    maxSprayAngle_(ChomiakDict_.lookup("maxSprayConeAngle"))
 {
 
     if (sm.injectors().size() != maxSprayAngle_.size())
@@ -77,24 +75,6 @@ ChomiakInjector::ChomiakInjector
             << abort(FatalError);
     }
 
-    forAll(sm.injectors(), i)
-    {
-        Random rndGen(label(0));
-        vector dir = sm.injectors()[i].properties()->direction();
-        scalar magV = 0.0;
-        vector tangent;
-        
-        while (magV < SMALL)
-        {
-            vector testThis = rndGen.vector01();
-            
-            tangent = testThis - (testThis & dir)*dir;
-            magV = mag(tangent);
-        }
-        
-        tan1_[i] = tangent/magV;
-        tan2_[i] = dir ^ tan1_[i];
-    }
     scalar referencePressure = sm.p().average().value();
 
     // correct velocityProfile
@@ -127,7 +107,8 @@ scalar ChomiakInjector::d0
 vector ChomiakInjector::direction
 (
     const label n,
-    const scalar,
+    const label hole,
+    const scalar time,
     const scalar d
 ) const
 {
@@ -162,13 +143,13 @@ vector ChomiakInjector::direction
     {
         normal = alpha*
         (
-            tan1_[n]*cos(beta) +
-            tan2_[n]*sin(beta)
+            injectors_[n].properties()->tan1(hole)*cos(beta) +
+            injectors_[n].properties()->tan2(hole)*sin(beta)
         );
     }
     
     // set the direction of injection by adding the normal vector
-    vector dir = dcorr*injectors_[n].properties()->direction() + normal;
+    vector dir = dcorr*injectors_[n].properties()->direction(hole, time) + normal;
     dir /= mag(dir);
 
     return dir;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.H
index 58194d6d977759e500ae9d1784ee98bcbd0170b1..fc754e63bf5b9919ac3f3824ce39770b2baaee9e 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/Chomiak/Chomiak.H
@@ -67,10 +67,6 @@ private:
         autoPtr<pdf> dropletPDF_;
         scalarList maxSprayAngle_;
 
-    // two perpendicular vectors, perpendicular to injection direction
-        List<vector> tan1_;
-        List<vector> tan2_;
-
 public:
 
     //- Runtime type information
@@ -100,7 +96,8 @@ public:
         //- Return the spray angle of the injector
         vector direction
         (
-            const label injector, 
+            const label injector,
+            const label hole, 
             const scalar time,
             const scalar d
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.C b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.C
index a15d1336916f6f47d24a3c615651fd9fba9a17f8..e6925c68c13bf48aa1a3cd09d8df3e038c0e1e44 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.C
@@ -67,9 +67,7 @@ blobsSwirlInjector::blobsSwirlInjector
     angle_(0.0),
     u_(0.0),
     x_(0.0),
-    h_(0.0),
-    tan1_(coneAngle_.size()),
-    tan2_(coneAngle_.size())
+    h_(0.0)
 {
 
     if (sm.injectors().size() != coneAngle_.size())
@@ -80,26 +78,6 @@ blobsSwirlInjector::blobsSwirlInjector
             << abort(FatalError);
     }
 
-
-    forAll(sm.injectors(), i)
-    {
-        Random rndGen(label(0));
-        vector dir = sm.injectors()[i].properties()->direction();
-        scalar magV = 0.0;
-        vector tangent;
-        
-        while (magV < SMALL)
-        {
-            vector testThis = rndGen.vector01();
-            
-            tangent = testThis - (testThis & dir)*dir;
-            magV = mag(tangent);
-        }
-        
-        tan1_[i] = tangent/magV;
-        tan2_[i] = dir ^ tan1_[i];
-    }
-
     scalar referencePressure = sm.p().average().value();
 
     // correct velocityProfile
@@ -156,8 +134,9 @@ scalar blobsSwirlInjector::d0
 vector blobsSwirlInjector::direction
 (
     const label n,
-    const scalar,
-    const scalar
+    const label hole,
+    const scalar time,
+    const scalar d
 ) const
 {
 
@@ -186,13 +165,13 @@ vector blobsSwirlInjector::direction
     {
         normal = alpha*
         (
-            tan1_[n]*cos(beta) +
-            tan2_[n]*sin(beta)
+            injectors_[n].properties()->tan1(hole)*cos(beta) +
+            injectors_[n].properties()->tan2(hole)*sin(beta)
         );
     }
     
     // set the direction of injection by adding the normal vector
-    vector dir = dcorr*injectors_[n].properties()->direction() + normal;
+    vector dir = dcorr*injectors_[n].properties()->direction(hole, time) + normal;
     dir /= mag(dir);
 
     return dir;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.H
index de0e2f278f719f44f55a39a31e1e0e8a0c0f615f..200a6cb1d42b1224f5285a5ab3ee3a04b817579b 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/blobsSwirl/blobsSwirlInjector.H
@@ -90,10 +90,6 @@ private:
 
         mutable scalar h_;
 
-    // two perpendicular vectors, perpendicular to injection direction
-        List<vector> tan1_;
-        List<vector> tan2_;
-
     // private member functions
 
         scalar kv
@@ -142,6 +138,7 @@ public:
         vector direction
         (
             const label injector,
+            const label hole,
             const scalar time,
             const scalar d
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.C b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.C
index b5bad617283ecc7c1a9b33af1de7b393550b9aac..52ca12ee6a73b071fa294dbf93eefb96c3c74a21 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.C
@@ -57,9 +57,7 @@ constInjector::constInjector
     injectorModel(dict, sm),
     specDict_(dict.subDict(typeName + "Coeffs")),
     dropletNozzleDiameterRatio_(specDict_.lookup("dropletNozzleDiameterRatio")),
-    sprayAngle_(specDict_.lookup("sprayAngle")),
-    tan1_(sprayAngle_.size()),
-    tan2_(sprayAngle_.size())
+    sprayAngle_(specDict_.lookup("sprayAngle"))
 {
     if (sm.injectors().size() != dropletNozzleDiameterRatio_.size())
     {
@@ -77,26 +75,6 @@ constInjector::constInjector
             << abort(FatalError);
     }
 
-
-    forAll(sm.injectors(), i)
-    {
-        Random rndGen(label(0));
-        vector dir = sm.injectors()[i].properties()->direction();
-        scalar magV = 0.0;
-        vector tangent;
-        
-        while (magV < SMALL)
-        {
-            vector testThis = rndGen.vector01();
-            
-            tangent = testThis - (testThis & dir)*dir;
-            magV = mag(tangent);
-        }
-        
-        tan1_[i] = tangent/magV;
-        tan2_[i] = dir ^ tan1_[i];
-    }
-
     scalar referencePressure = sm.p().average().value();
 
     // correct velocity and pressure profiles
@@ -129,12 +107,11 @@ scalar constInjector::d0
 vector constInjector::direction
 (
     const label n,
-    const scalar,
-    const scalar
+    const label hole,
+    const scalar time,
+    const scalar d
 ) const
 {
-    //    return sprayAngle_[n];
-
 
     /*
         randomly distribute parcels in a solid cone
@@ -179,13 +156,13 @@ vector constInjector::direction
     {
         normal = alpha*
         (
-            tan1_[n]*cos(beta) +
-            tan2_[n]*sin(beta)
+            injectors_[n].properties()->tan1(hole)*cos(beta) +
+            injectors_[n].properties()->tan2(hole)*sin(beta)
         );
     }
     
     // set the direction of injection by adding the normal vector
-    vector dir = dcorr*injectors_[n].properties()->direction() + normal;
+    vector dir = dcorr*injectors_[n].properties()->direction(n, time) + normal;
     dir /= mag(dir);
 
     return dir;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.H
index 67fe880c0f3c023e506a8081443634c547ca11d6..4a41808c4b6153fb4f6b63fce7199e3b29edd64d 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/constant/constInjector.H
@@ -59,10 +59,6 @@ private:
         scalarList dropletNozzleDiameterRatio_;
         scalarList sprayAngle_;
 
-    // two perpendicular vectors, perpendicular to injection direction
-        List<vector> tan1_;
-        List<vector> tan2_;
-
 public:
 
     //- Runtime type information
@@ -93,6 +89,7 @@ public:
         vector direction
         (
             const label injector, 
+            const label hole,
             const scalar time,
             const scalar d
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.C b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.C
index 6f06bb17528c7806a5fcfb30eb079244eb0e1711..740108fdf5ccd5db75d9e8b5e9bac463dd23a72a 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.C
@@ -65,9 +65,7 @@ definedHollowConeInjector::definedHollowConeInjector
         )
     ),
     innerConeAngle_(definedHollowConeDict_.lookup("innerConeAngle")),
-    outerConeAngle_(definedHollowConeDict_.lookup("outerConeAngle")),
-    tan1_(sm.injectors().size()),
-    tan2_(sm.injectors().size())
+    outerConeAngle_(definedHollowConeDict_.lookup("outerConeAngle"))
 {
 
     // convert CA to real time - inner cone angle
@@ -109,26 +107,6 @@ definedHollowConeInjector::definedHollowConeInjector
              << abort(FatalError);
     }
 
-    // initialise injectors
-    forAll(sm.injectors(), i)
-    {
-        Random rndGen(label(0));
-        vector dir = sm.injectors()[i].properties()->direction();
-        scalar magV = 0.0;
-        vector tangent;
-        
-        while (magV < SMALL)
-        {
-            vector testThis = rndGen.vector01();
-            
-            tangent = testThis - (testThis & dir)*dir;
-            magV = mag(tangent);
-        }
-        
-        tan1_[i] = tangent/magV;
-        tan2_[i] = dir ^ tan1_[i];
-    }
-
     scalar referencePressure = sm.p().average().value();
     // correct pressureProfile
     forAll(sm.injectors(), i)
@@ -162,6 +140,7 @@ scalar definedHollowConeInjector::d0
 vector definedHollowConeInjector::direction
 (
     const label n,
+    const label hole,
     const scalar t,
     const scalar d
 ) const
@@ -201,13 +180,13 @@ vector definedHollowConeInjector::direction
     {
         normal = alpha*
         (
-            tan1_[n]*cos(beta) +
-            tan2_[n]*sin(beta)
+            injectors_[n].properties()->tan1(hole)*cos(beta) +
+            injectors_[n].properties()->tan2(hole)*sin(beta)
         );
     }
     
     // set the direction of injection by adding the normal vector
-    vector dir = dcorr*injectors_[n].properties()->direction() + normal;
+    vector dir = dcorr*injectors_[n].properties()->direction(hole, t) + normal;
     // normailse direction vector
     dir /= mag(dir);
 
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.H
index 81023e49a9df4e1b8c961d7f8b326a51d87b80c4..83d9613667ca05ef8c26faf1f79eec711d2b6480 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedHollowCone/definedHollowCone.H
@@ -68,11 +68,6 @@ private:
         List<pair> innerConeAngle_;
         List<pair> outerConeAngle_;
 
-        // two perpendicular vectors, perpendicular to injection direction
-        List<vector> tan1_;
-        List<vector> tan2_;
-
-
 public:
 
     //- Runtime type information
@@ -102,7 +97,8 @@ public:
         //- Return the spray angle of the injector
         vector direction
         (
-            const label injector, 
+            const label injector,
+            const label hole,
             const scalar time,
             const scalar d
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.C b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.C
index 730eee8240c8f0c98780260d4d90e5c954a59f56..c3f56b07f985e5ab3c3ea7b8a54f51c3fff0a0db 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.C
@@ -61,32 +61,9 @@ definedPressureSwirlInjector::definedPressureSwirlInjector
     coneInterval_(definedPressureSwirlInjectorDict_.lookup("ConeInterval")),
     maxKv_(definedPressureSwirlInjectorDict_.lookup("maxKv")),
 
-    angle_(0.0),
-    tan1_(coneAngle_.size()),
-    tan2_(coneAngle_.size())
+    angle_(0.0)
 {
 
-
-
-    forAll(sm.injectors(), i)
-    {
-        Random rndGen(label(0));
-        vector dir = sm.injectors()[i].properties()->direction();
-        scalar magV = 0.0;
-        vector tangent;
-        
-        while (magV < SMALL)
-        {
-            vector testThis = rndGen.vector01();
-            
-            tangent = testThis - (testThis & dir)*dir;
-            magV = mag(tangent);
-        }
-        
-        tan1_[i] = tangent/magV;
-        tan2_[i] = dir ^ tan1_[i];
-    }
-
     scalar referencePressure = sm.p().average().value();
 
     // correct velocityProfile
@@ -219,8 +196,9 @@ scalar definedPressureSwirlInjector::d0
 vector definedPressureSwirlInjector::direction
 (
     const label n,
-    const scalar,
-    const scalar
+    const label hole,
+    const scalar time,
+    const scalar d
 ) const
 {
 
@@ -249,13 +227,13 @@ vector definedPressureSwirlInjector::direction
     {
         normal = alpha*
         (
-            tan1_[n]*cos(beta) +
-            tan2_[n]*sin(beta)
+            injectors_[n].properties()->tan1(hole)*cos(beta) +
+            injectors_[n].properties()->tan2(hole)*sin(beta)
         );
     }
     
     // set the direction of injection by adding the normal vector
-    vector dir = dcorr*injectors_[n].properties()->direction() + normal;
+    vector dir = dcorr*injectors_[n].properties()->direction(hole, time) + normal;
     dir /= mag(dir);
 
     return dir;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.H
index 14415d789b343fbd88fbd71f6bc6611d5a5d4e85..e84dd8f43ef618eb766c0263b60e55e2f5cd13f4 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/definedPressureSwirl/definedPressureSwirl.H
@@ -71,10 +71,6 @@ private:
     // The initial velocity for the parcels
         mutable scalar u_;
 
-    // two perpendicular vectors, perpendicular to injection direction
-        List<vector> tan1_;
-        List<vector> tan2_;
-
     // private member functions
 
         scalar kv
@@ -118,6 +114,7 @@ public:
         vector direction
         (
             const label injector,
+            const label hole,
             const scalar time,
             const scalar d
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.C b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.C
index 2121a803f11ee6c3ecd9859805069b37e284052f..001102e4f8591470fce95d690f34bb23ddcc2cba 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.C
@@ -65,9 +65,7 @@ hollowConeInjector::hollowConeInjector
         )
     ),
     innerAngle_(hollowConeDict_.lookup("innerConeAngle")),
-    outerAngle_(hollowConeDict_.lookup("outerConeAngle")),
-    tan1_(outerAngle_.size()),
-    tan2_(outerAngle_.size())
+    outerAngle_(hollowConeDict_.lookup("outerConeAngle"))
 {
 
     if (sm.injectors().size() != innerAngle_.size())
@@ -86,26 +84,8 @@ hollowConeInjector::hollowConeInjector
             << abort(FatalError);
     }
 
-    forAll(sm.injectors(), i)
-    {
-        Random rndGen(label(0));
-        vector dir = sm.injectors()[i].properties()->direction();
-        scalar magV = 0.0;
-        vector tangent;
-        
-        while (magV < SMALL)
-        {
-            vector testThis = rndGen.vector01();
-            
-            tangent = testThis - (testThis & dir)*dir;
-            magV = mag(tangent);
-        }
-        
-        tan1_[i] = tangent/magV;
-        tan2_[i] = dir ^ tan1_[i];
-    }
-
     scalar referencePressure = sm.ambientPressure();
+
     // correct velocityProfile
     forAll(sm.injectors(), i)
     {
@@ -136,8 +116,9 @@ scalar hollowConeInjector::d0
 vector hollowConeInjector::direction
 (
     const label n,
-    const scalar,
-    const scalar
+    const label hole,
+    const scalar time,
+    const scalar d
 ) const
 {
     scalar angle = innerAngle_[n] + rndGen_.scalar01()*(outerAngle_[n]-innerAngle_[n]);
@@ -166,13 +147,13 @@ vector hollowConeInjector::direction
     {
         normal = alpha*
         (
-            tan1_[n]*cos(beta) +
-            tan2_[n]*sin(beta)
+            injectors_[n].properties()->tan1(hole)*cos(beta) +
+            injectors_[n].properties()->tan2(hole)*sin(beta)
         );
     }
     
     // set the direction of injection by adding the normal vector
-    vector dir = dcorr*injectors_[n].properties()->direction() + normal;
+    vector dir = dcorr*injectors_[n].properties()->direction(hole, time) + normal;
     dir /= mag(dir);
 
     return dir;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.H
index e53b6afcbb9b9135cc88addbf79d319ac1a836ad..c0469b997c73ff4eb9d384a76430256cf5d1a0bc 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/hollowCone/hollowCone.H
@@ -63,10 +63,6 @@ private:
         scalarList innerAngle_;
         scalarList outerAngle_;
 
-    // two perpendicular vectors, perpendicular to injection direction
-        List<vector> tan1_;
-        List<vector> tan2_;
-
 public:
 
     //- Runtime type information
@@ -96,7 +92,8 @@ public:
         //- Return the spray angle of the injector
         vector direction
         (
-            const label injector, 
+            const label injector,
+            const label hole,
             const scalar time,
             const scalar d
         ) const;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/injectorModel/injectorModel.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/injectorModel/injectorModel.H
index b099208287ce7c304bac71a97752c2c92893241d..22bd419508a32dea6d20ffc66332b48af5d8ff74 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/injectorModel/injectorModel.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/injectorModel/injectorModel.H
@@ -118,6 +118,7 @@ public:
         virtual vector direction
         (
             const label injector, 
+            const label hole,
             const scalar time,
             const scalar d
         ) const = 0;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.C b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.C
index 83dba60703e261c8652b154e9ac5ee52349b233d..86d419c8d8d798b84bf55bd2993c109563640187 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.C
@@ -61,9 +61,7 @@ pressureSwirlInjector::pressureSwirlInjector
     coneInterval_(pressureSwirlInjectorDict_.lookup("ConeInterval")),
     maxKv_(pressureSwirlInjectorDict_.lookup("maxKv")),
 
-    angle_(0.0),
-    tan1_(coneAngle_.size()),
-    tan2_(coneAngle_.size())
+    angle_(0.0)
 {
 
     if (sm.injectors().size() != coneAngle_.size())
@@ -74,26 +72,6 @@ pressureSwirlInjector::pressureSwirlInjector
             << abort(FatalError);
     }
 
-
-    forAll(sm.injectors(), i)
-    {
-        Random rndGen(label(0));
-        vector dir = sm.injectors()[i].properties()->direction();
-        scalar magV = 0.0;
-        vector tangent;
-        
-        while (magV < SMALL)
-        {
-            vector testThis = rndGen.vector01();
-            
-            tangent = testThis - (testThis & dir)*dir;
-            magV = mag(tangent);
-        }
-        
-        tan1_[i] = tangent/magV;
-        tan2_[i] = dir ^ tan1_[i];
-    }
-
     scalar referencePressure = sm.p().average().value();
 
     // correct velocityProfile
@@ -147,8 +125,9 @@ scalar pressureSwirlInjector::d0
 vector pressureSwirlInjector::direction
 (
     const label n,
-    const scalar,
-    const scalar
+    const label hole,
+    const scalar time,
+    const scalar d
 ) const
 {
 
@@ -177,13 +156,13 @@ vector pressureSwirlInjector::direction
     {
         normal = alpha*
         (
-            tan1_[n]*cos(beta) +
-            tan2_[n]*sin(beta)
+            injectors_[n].properties()->tan1(hole)*cos(beta) +
+            injectors_[n].properties()->tan2(hole)*sin(beta)
         );
     }
     
     // set the direction of injection by adding the normal vector
-    vector dir = dcorr*injectors_[n].properties()->direction() + normal;
+    vector dir = dcorr*injectors_[n].properties()->direction(hole, time) + normal;
     dir /= mag(dir);
 
     return dir;
diff --git a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.H b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.H
index d57312871377c9bc0aa57dcc0a60c2034185a1a4..91e5bc9f1a118ec1fc30be5f4a926dd9006ea891 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.H
+++ b/src/lagrangian/dieselSpray/spraySubModels/injectorModel/pressureSwirl/pressureSwirlInjector.H
@@ -68,10 +68,6 @@ private:
     // The initial velocity for the parcels           
         mutable scalar u_;
 
-    // two perpendicular vectors, perpendicular to injection direction
-        List<vector> tan1_;
-        List<vector> tan2_;
-
     // private member functions
         
         scalar kv
@@ -113,7 +109,8 @@ public:
         //- Return the spray angle of the injector
         vector direction
         (
-            const label injector, 
+            const label injector,
+            const label hole, 
             const scalar time,
             const scalar d
         ) const;
diff --git a/src/thermophysicalModels/liquids/CH4N2O/CH4N2O.C b/src/thermophysicalModels/liquids/CH4N2O/CH4N2O.C
new file mode 100644
index 0000000000000000000000000000000000000000..5b7e4299c0424da370f0254735c77a4d26f1445e
--- /dev/null
+++ b/src/thermophysicalModels/liquids/CH4N2O/CH4N2O.C
@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+
+-------------------------------------------------------------------------------
+*/
+
+#include "CH4N2O.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(CH4N2O, 0);
+addToRunTimeSelectionTable(liquid, CH4N2O,);
+addToRunTimeSelectionTable(liquid, CH4N2O, Istream);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/liquids/CH4N2O/CH4N2O.H b/src/thermophysicalModels/liquids/CH4N2O/CH4N2O.H
new file mode 100644
index 0000000000000000000000000000000000000000..25c336f32875542a6132e41320fa773bdd05b9d3
--- /dev/null
+++ b/src/thermophysicalModels/liquids/CH4N2O/CH4N2O.H
@@ -0,0 +1,285 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::CH4N2O
+
+Description
+    urea, note that some of the properties are unavailable in the literature and have been copied from water.
+
+SourceFiles
+    CH4N2O.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CH4N2O_H
+#define CH4N2O_H
+
+#include "liquid.H"
+#include "NSRDSfunc0.H"
+#include "NSRDSfunc1.H"
+#include "NSRDSfunc2.H"
+#include "NSRDSfunc3.H"
+#include "NSRDSfunc4.H"
+#include "NSRDSfunc5.H"
+#include "NSRDSfunc6.H"
+#include "NSRDSfunc7.H"
+#include "NSRDSfunc14.H"
+#include "APIdiffCoefFunc.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class CH4N2O Declaration
+\*---------------------------------------------------------------------------*/
+
+class CH4N2O
+:
+    public liquid
+{
+    // Private data
+
+        NSRDSfunc0 rho_;
+        NSRDSfunc1 pv_;
+        NSRDSfunc6 hl_;
+        NSRDSfunc0 cp_;
+        NSRDSfunc0 h_;
+        NSRDSfunc7 cpg_;
+        NSRDSfunc4 B_;
+        NSRDSfunc1 mu_;
+        NSRDSfunc2 mug_;
+        NSRDSfunc0 K_;
+        NSRDSfunc2 Kg_;
+        NSRDSfunc6 sigma_;
+        APIdiffCoefFunc D_;
+
+public:
+
+    //- Runtime type information
+    TypeName("CH4N2O");
+
+
+    // Constructors
+
+        //- Construct null
+        CH4N2O()
+        :
+            liquid(60.056, 705.0, 9.050e+6, 0.218, 0.337, 405.85, 9.3131e+1, 465.0, 1.52e-29, 0.3449, 4.7813e+4),
+            rho_(1230.006936, 0, 0, 0, 0, 0),
+            pv_(12.06, -3992.0, 0, 0, 0),
+//            hl_(1463034.50113228, 0, 0, 0, 0, 0), 
+// NN.      we cant use constant heat of vapourisation, the below value is linear (sqrt) interpolation to critical temp
+            hl_(705.0, 2534249.0, 0.5, 0.0, 0.0, 0.0),
+            cp_(2006.46063673904, 0, 0, 0, 0, 0),
+            // NN: enthalpy, h_, is not used in the sprayModel.
+            // For consistency, the enthalpy is derived from hlat and hl.
+            // It is, however, convenient to have it available.
+            h_(-6154107.41641135, 2006.46063673904, 0, 0, 0, 0),
+
+            cpg_(811.875582789397, 2099.04089516451, 1627.3, 1603.63660583455, 724.41),
+            B_(-0.000383641934194752, 0.447249234048222, -469062.208605302, 5.5628080458239e+18, -2.3040162514986e+21),
+            mu_(-51.964, 3670.6, 5.7331, -5.3495e-29, 10),
+            mug_(2.6986e-06, 0.498, 1257.7, -19570),
+            K_(-0.4267, 0.0056903, -8.0065e-06, 1.815e-09, 0, 0),
+            Kg_(6.977e-05, 1.1243, 844.9, -148850),
+            sigma_(705.0, 1.0, 0.0, 0.0, 0.0, 0), // set to constant
+            D_(147.18, 20.1, 60.056, 28) // NN: Same as nHeptane
+        {}
+        CH4N2O
+        (
+            const liquid& l,
+            const NSRDSfunc0& density,
+            const NSRDSfunc1& vapourPressure,
+            const NSRDSfunc6& heatOfVapourisation,
+            const NSRDSfunc0& heatCapacity,
+            const NSRDSfunc0& enthalpy,
+            const NSRDSfunc7& idealGasHeatCapacity,
+            const NSRDSfunc4& secondVirialCoeff,
+            const NSRDSfunc1& dynamicViscosity,
+            const NSRDSfunc2& vapourDynamicViscosity,
+            const NSRDSfunc0& thermalConductivity,
+            const NSRDSfunc2& vapourThermalConductivity,
+            const NSRDSfunc6& surfaceTension,
+            const APIdiffCoefFunc& vapourDiffussivity
+        )
+        :
+            liquid(l),
+            rho_(density),
+            pv_(vapourPressure),
+            hl_(heatOfVapourisation),
+            cp_(heatCapacity),
+            h_(enthalpy),
+            cpg_(idealGasHeatCapacity),
+            B_(secondVirialCoeff),
+            mu_(dynamicViscosity),
+            mug_(vapourDynamicViscosity),
+            K_(thermalConductivity),
+            Kg_(vapourThermalConductivity),
+            sigma_(surfaceTension),
+            D_(vapourDiffussivity)
+        {}
+
+        //- Construct from Istream
+        CH4N2O(Istream& is)
+        :
+            liquid(is),
+            rho_(is),
+            pv_(is),
+            hl_(is),
+            cp_(is),
+            h_(is),
+            cpg_(is),
+            B_(is),
+            mu_(is),
+            mug_(is),
+            K_(is),
+            Kg_(is),
+            sigma_(is),
+            D_(is)
+        {}
+
+
+    // Member Functions
+
+        //- Liquid density [kg/m^3]
+        scalar rho(scalar p, scalar T) const
+        {
+            return rho_.f(p, T);
+        }
+
+        //- Vapour pressure [Pa]
+        scalar pv(scalar p, scalar T) const
+        {
+            return pv_.f(p, T);
+        }
+
+        //- Heat of vapourisation [J/kg]
+        scalar hl(scalar p, scalar T) const
+        {
+            return hl_.f(p, T);
+        }
+
+        //- Liquid heat capacity [J/(kg K)]
+        scalar cp(scalar p, scalar T) const
+        {
+            return cp_.f(p, T);
+        }
+
+        //- Liquid Enthalpy [J/(kg)]
+        scalar h(scalar p, scalar T) const
+        {
+            return h_.f(p, T);
+        }
+
+        //- Ideal gas heat capacity [J/(kg K)]
+        scalar cpg(scalar p, scalar T) const
+        {
+            return cpg_.f(p, T);
+        }
+
+        //- Second Virial Coefficient [m^3/kg]
+        scalar B(scalar p, scalar T) const
+        {
+            return B_.f(p, T);
+        }
+
+        //- Liquid viscosity [Pa s]
+        scalar mu(scalar p, scalar T) const
+        {
+            return mu_.f(p, T);
+        }
+
+        //- Vapour viscosity [Pa s]
+        scalar mug(scalar p, scalar T) const
+        {
+            return mug_.f(p, T);
+        }
+
+        //- Liquid thermal conductivity  [W/(m K)]
+        scalar K(scalar p, scalar T) const
+        {
+            return K_.f(p, T);
+        }
+
+        //- Vapour thermal conductivity  [W/(m K)]
+        scalar Kg(scalar p, scalar T) const
+        {
+            return Kg_.f(p, T);
+        }
+
+        //- Surface tension [N/m]
+        scalar sigma(scalar p, scalar T) const
+        {
+            return sigma_.f(p, T);
+        }
+
+        //- Vapour diffussivity [m2/s]
+        scalar D(scalar p, scalar T) const
+        {
+            return D_.f(p, T);
+        }
+
+
+        //- Write the function coefficients
+        void writeData(Ostream& os) const
+        {
+            liquid::writeData(os); os << nl;
+            rho_.writeData(os); os << nl;
+            pv_.writeData(os); os << nl;
+            hl_.writeData(os); os << nl;
+            cp_.writeData(os); os << nl;
+            cpg_.writeData(os); os << nl;
+            B_.writeData(os); os << nl;
+            mu_.writeData(os); os << nl;
+            mug_.writeData(os); os << nl;
+            K_.writeData(os); os << nl;
+            Kg_.writeData(os); os << nl;
+            sigma_.writeData(os); os << nl;
+            D_.writeData(os); os << endl;
+        }
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream& os, const CH4N2O& l)
+        {
+            l.writeData(os);
+            return os;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/thermophysicalModels/liquids/Make/files b/src/thermophysicalModels/liquids/Make/files
index 7d55ee40b759d37ad4064dc6cb66c51a5b0c5463..e2f27128d520e81ed9e7d4e8b9dee691724f1483 100644
--- a/src/thermophysicalModels/liquids/Make/files
+++ b/src/thermophysicalModels/liquids/Make/files
@@ -26,5 +26,8 @@ C2H5OH/C2H5OH.C
 Ar/Ar.C
 N2/N2.C
 MB/MB.C
+CH4N2O/CH4N2O.C
+nC3H8O/nC3H8O.C
+iC3H8O/iC3H8O.C
 
 LIB = $(FOAM_LIBBIN)/libliquids
diff --git a/src/thermophysicalModels/liquids/iC3H8O/iC3H8O.C b/src/thermophysicalModels/liquids/iC3H8O/iC3H8O.C
new file mode 100644
index 0000000000000000000000000000000000000000..2deb4001fc8e99632ed4053714be9159ca6e5c32
--- /dev/null
+++ b/src/thermophysicalModels/liquids/iC3H8O/iC3H8O.C
@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+
+-------------------------------------------------------------------------------
+*/
+
+#include "iC3H8O.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(iC3H8O, 0);
+addToRunTimeSelectionTable(liquid, iC3H8O,);
+addToRunTimeSelectionTable(liquid, iC3H8O, Istream);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/liquids/iC3H8O/iC3H8O.H b/src/thermophysicalModels/liquids/iC3H8O/iC3H8O.H
new file mode 100644
index 0000000000000000000000000000000000000000..4c286518879884a819d08de4445356fdcdfd3291
--- /dev/null
+++ b/src/thermophysicalModels/liquids/iC3H8O/iC3H8O.H
@@ -0,0 +1,269 @@
+// The FOAM Project // File: iC3H8O/iC3H8O.H
+/*
+-------------------------------------------------------------------------------
+ =========         | Class Interface
+ \\      /         |
+  \\    /          | Name:   iC3H8O
+   \\  /           | Family: iC3H8O
+    \\/            |
+    F ield         | FOAM version: 2.2
+    O peration     |
+    A and          | Copyright (C) 1991-2000 Nabla Ltd.
+    M anipulation  |          All Rights Reserved.
+-------------------------------------------------------------------------------
+CLASS
+    iC3H8O
+
+DESCRIPTION
+    iso-propanol
+
+
+*/
+// ------------------------------------------------------------------------- //
+
+#ifndef iC3H8O_H
+#define iC3H8O_H
+
+#include "liquid.H"
+#include "NSRDSfunc0.H"
+#include "NSRDSfunc1.H"
+#include "NSRDSfunc2.H"
+#include "NSRDSfunc3.H"
+#include "NSRDSfunc4.H"
+#include "NSRDSfunc5.H"
+#include "NSRDSfunc6.H"
+#include "NSRDSfunc7.H"
+#include "NSRDSfunc14.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class iC3H8O Declaration
+\*---------------------------------------------------------------------------*/
+
+class iC3H8O
+:
+    public liquid
+{
+    // Private data
+
+        NSRDSfunc5 rho_;
+        NSRDSfunc1 pv_;
+        NSRDSfunc6 hl_;
+        NSRDSfunc0 cp_;
+        NSRDSfunc0 h_;
+        NSRDSfunc7 cpg_;
+        NSRDSfunc4 B_;
+        NSRDSfunc1 mu_;
+        NSRDSfunc2 mug_;
+        NSRDSfunc0 K_;
+        NSRDSfunc2 Kg_;
+        NSRDSfunc0 sigma_;
+        NSRDSfunc1 D_;
+
+public:
+
+    //- Runtime type information
+    TypeName("iC3H8O");
+
+
+    // Constructors
+
+        //- Construct null
+        iC3H8O()
+        :
+            liquid(60.096, 508.31, 4.7643e+6, 0.22013, 0.248, 185.28, 3.20e-2, 355.41, 5.5372e-30, 0.6689, 2.3575e+4),
+            rho_(70.91328, 0.26475, 508.31, 0.243),
+            pv_(92.935, -8177.1, -10.031, 3.9988e-06, 2),
+            hl_(508.31, 948149.627263046, 0.087, 0.3007, 0, 0),
+            cp_(7760.91586794462, -68.3672790202343, 0.241380457933972, -0.000235057241746539, 0, 0),
+            // NN: enthalpy, h_, is not used in the sprayModel.
+            // For consistency, the enthalpy is derived from hlat and hl.
+            // It is, however, convenient to have it available.
+            h_(-6227786.27583977, 7760.91586794462, -34.1836395101172, 0.0804601526446574, -5.87643104366347e-05, 0),
+            cpg_(789.73642172524, 3219.8482428115, 1124, 1560.83599574015, 460),
+            B_(0.000502529286474973, -0.104665867944622, -717185.83599574, 3.3047124600639e+18, -1.43270766773163e+21),
+            mu_(-8.23, 2282.2, -0.98495, 0, 0),
+            mug_(1.993e-07, 0.7233, 178, 0),
+            K_(0.2029, -0.0002278, 0, 0, 0, 0),
+            Kg_(-80.642, -1.4549, -604.42, 0),
+            sigma_(0.03818, -3.818e-05, -6.51e-08, 0, 0, 0),
+            D_(4.75e-10, 1.75, 0.0, 0.0, 0.0) // NN. same as iC3H8O
+        {}
+        iC3H8O
+        (
+            const liquid& l,
+            const NSRDSfunc5& density,
+            const NSRDSfunc1& vapourPressure,
+            const NSRDSfunc6& heatOfVapourisation,
+            const NSRDSfunc0& heatCapacity,
+            const NSRDSfunc0& enthalpy,
+            const NSRDSfunc7& idealGasHeatCapacity,
+            const NSRDSfunc4& secondVirialCoeff,
+            const NSRDSfunc1& dynamicViscosity,
+            const NSRDSfunc2& vapourDynamicViscosity,
+            const NSRDSfunc0& thermalConductivity,
+            const NSRDSfunc2& vapourThermalConductivity,
+            const NSRDSfunc0& surfaceTension,
+            const NSRDSfunc1& vapourDiffussivity
+        )
+        :
+            liquid(l),
+            rho_(density),
+            pv_(vapourPressure),
+            hl_(heatOfVapourisation),
+            cp_(heatCapacity),
+            h_(enthalpy),
+            cpg_(idealGasHeatCapacity),
+            B_(secondVirialCoeff),
+            mu_(dynamicViscosity),
+            mug_(vapourDynamicViscosity),
+            K_(thermalConductivity),
+            Kg_(vapourThermalConductivity),
+            sigma_(surfaceTension),
+            D_(vapourDiffussivity)
+        {}
+
+        //- Construct from Istream
+        iC3H8O(Istream& is)
+        :
+            liquid(is),
+            rho_(is),
+            pv_(is),
+            hl_(is),
+            cp_(is),
+            h_(is),
+            cpg_(is),
+            B_(is),
+            mu_(is),
+            mug_(is),
+            K_(is),
+            Kg_(is),
+            sigma_(is),
+            D_(is)
+        {}
+
+
+    // Member Functions
+
+        //- Liquid density [kg/m^3]
+        scalar rho(scalar p, scalar T) const
+        {
+            return rho_.f(p, T);
+        }
+
+        //- Vapour pressure [Pa]
+        scalar pv(scalar p, scalar T) const
+        {
+            return pv_.f(p, T);
+        }
+
+        //- Heat of vapourisation [J/kg]
+        scalar hl(scalar p, scalar T) const
+        {
+            return hl_.f(p, T);
+        }
+
+        //- Liquid heat capacity [J/(kg K)]
+        scalar cp(scalar p, scalar T) const
+        {
+            return cp_.f(p, T);
+        }
+
+        //- Liquid Enthalpy [J/(kg)]
+        scalar h(scalar p, scalar T) const
+        {
+            return h_.f(p, T);
+        }
+
+        //- Ideal gas heat capacity [J/(kg K)]
+        scalar cpg(scalar p, scalar T) const
+        {
+            return cpg_.f(p, T);
+        }
+
+        //- Second Virial Coefficient [m^3/kg]
+        scalar B(scalar p, scalar T) const
+        {
+            return B_.f(p, T);
+        }
+
+        //- Liquid viscosity [Pa s]
+        scalar mu(scalar p, scalar T) const
+        {
+            return mu_.f(p, T);
+        }
+
+        //- Vapour viscosity [Pa s]
+        scalar mug(scalar p, scalar T) const
+        {
+            return mug_.f(p, T);
+        }
+
+        //- Liquid thermal conductivity  [W/(m K)]
+        scalar K(scalar p, scalar T) const
+        {
+            return K_.f(p, T);
+        }
+
+        //- Vapour thermal conductivity  [W/(m K)]
+        scalar Kg(scalar p, scalar T) const
+        {
+            return Kg_.f(p, T);
+        }
+
+        //- Surface tension [N/m]
+        scalar sigma(scalar p, scalar T) const
+        {
+            return sigma_.f(p, T);
+        }
+
+        //- Vapour diffussivity [m2/s]
+        scalar D(scalar p, scalar T) const
+        {
+            return D_.f(p, T);
+        }
+
+
+        //- Write the function coefficients
+        void writeData(Ostream& os) const
+        {
+            liquid::writeData(os); os << nl;
+            rho_.writeData(os); os << nl;
+            pv_.writeData(os); os << nl;
+            hl_.writeData(os); os << nl;
+            cp_.writeData(os); os << nl;
+            cpg_.writeData(os); os << nl;
+            B_.writeData(os); os << nl;
+            mu_.writeData(os); os << nl;
+            mug_.writeData(os); os << nl;
+            K_.writeData(os); os << nl;
+            Kg_.writeData(os); os << nl;
+            sigma_.writeData(os); os << nl;
+            D_.writeData(os); os << endl;
+        }
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream& os, const iC3H8O& l)
+        {
+            l.writeData(os);
+            return os;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/thermophysicalModels/liquids/nC3H8O/nC3H8O.C b/src/thermophysicalModels/liquids/nC3H8O/nC3H8O.C
new file mode 100644
index 0000000000000000000000000000000000000000..dd6e3ea75492164af7c58ea9b829cf8c0b76b32d
--- /dev/null
+++ b/src/thermophysicalModels/liquids/nC3H8O/nC3H8O.C
@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+
+-------------------------------------------------------------------------------
+*/
+
+#include "nC3H8O.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(nC3H8O, 0);
+addToRunTimeSelectionTable(liquid, nC3H8O,);
+addToRunTimeSelectionTable(liquid, nC3H8O, Istream);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/liquids/nC3H8O/nC3H8O.H b/src/thermophysicalModels/liquids/nC3H8O/nC3H8O.H
new file mode 100644
index 0000000000000000000000000000000000000000..2278e5dda24bf82804824398c81c903fc7eff92c
--- /dev/null
+++ b/src/thermophysicalModels/liquids/nC3H8O/nC3H8O.H
@@ -0,0 +1,267 @@
+/*
+-------------------------------------------------------------------------------
+ =========         | Class Interface
+ \\      /         |
+  \\    /          | Name:   nC3H8O
+   \\  /           | Family: nC3H8O
+    \\/            |
+    F ield         | FOAM version: 2.2
+    O peration     |
+    A and          | Copyright (C) 1991-2000 Nabla Ltd.
+    M anipulation  |          All Rights Reserved.
+-------------------------------------------------------------------------------
+CLASS
+    nC3H8O
+
+DESCRIPTION
+    propanol 
+
+*/
+// ------------------------------------------------------------------------- //
+
+#ifndef nC3H8O_H
+#define nC3H8O_H
+
+#include "liquid.H"
+#include "NSRDSfunc0.H"
+#include "NSRDSfunc1.H"
+#include "NSRDSfunc2.H"
+#include "NSRDSfunc3.H"
+#include "NSRDSfunc4.H"
+#include "NSRDSfunc5.H"
+#include "NSRDSfunc6.H"
+#include "NSRDSfunc7.H"
+#include "NSRDSfunc14.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class nC3H8O Declaration
+\*---------------------------------------------------------------------------*/
+
+class nC3H8O
+:
+    public liquid
+{
+    // Private data
+
+        NSRDSfunc5 rho_;
+        NSRDSfunc1 pv_;
+        NSRDSfunc6 hl_;
+        NSRDSfunc0 cp_;
+        NSRDSfunc0 h_;
+        NSRDSfunc7 cpg_;
+        NSRDSfunc4 B_;
+        NSRDSfunc1 mu_;
+        NSRDSfunc2 mug_;
+        NSRDSfunc0 K_;
+        NSRDSfunc2 Kg_;
+        NSRDSfunc0 sigma_;
+        NSRDSfunc1 D_;
+
+public:
+
+    //- Runtime type information
+    TypeName("nC3H8O");
+
+
+    // Constructors
+
+        //- Construct null
+        nC3H8O()
+        :
+            liquid(60.096, 536.71, 5.1696e+6, 0.21853, 0.253, 146.95, 6.5112e-7, 370.35, 5.6039e-30, 0.6279, 2.4557e+4),
+            rho_(75.300288, 0.272, 536.71, 0.2494),
+            pv_(77.46, -7960, -7.5235, 3e-07, 2),
+            hl_(536.71, 1098242.8115016, 0.647, -0.783, 0.613, 0),
+            cp_(216.320553780618, 18.5203674121406, -0.0751797124600639, 0.000126464323748669, 0, 0),
+            // NN: enthalpy, h_, is not used in the sprayModel.
+            // For consistency, the enthalpy is derived from hlat and hl.
+            // It is, however, convenient to have it available.
+            h_(-5533091.96851587, 216.320553780618, 9.26018370607029, -0.0250599041533546, 3.16160809371672e-05, 0),
+            cpg_(961.794462193823, 3467.78487752929, 1542, 2046.72523961661, 649),
+            B_(0.000933506389776358, -1.09325079872204, -531649.361022364, -2.32627795527157e+17, -3.81888977635783e+20),
+            mu_(0.571, 1521, -2.0894, 0, 0),
+            mug_(7.942e-07, 0.5491, 415.8, 0),
+            K_(0.204, -0.000169, 0, 0, 0, 0),
+            Kg_(-613.84, 0.7927, -1157400000, 0),
+            sigma_(0.04533, -6.88e-05, -1.6e-08, 0, 0, 0),
+            D_(4.75e-10, 1.75, 0.0, 0.0, 0.0) // NN. same as iC3H8O
+        {}
+        nC3H8O
+        (
+            const liquid& l,
+            const NSRDSfunc5& density,
+            const NSRDSfunc1& vapourPressure,
+            const NSRDSfunc6& heatOfVapourisation,
+            const NSRDSfunc0& heatCapacity,
+            const NSRDSfunc0& enthalpy,
+            const NSRDSfunc7& idealGasHeatCapacity,
+            const NSRDSfunc4& secondVirialCoeff,
+            const NSRDSfunc1& dynamicViscosity,
+            const NSRDSfunc2& vapourDynamicViscosity,
+            const NSRDSfunc0& thermalConductivity,
+            const NSRDSfunc2& vapourThermalConductivity,
+            const NSRDSfunc0& surfaceTension,
+            const NSRDSfunc1& vapourDiffussivity
+        )
+        :
+            liquid(l),
+            rho_(density),
+            pv_(vapourPressure),
+            hl_(heatOfVapourisation),
+            cp_(heatCapacity),
+            h_(enthalpy),
+            cpg_(idealGasHeatCapacity),
+            B_(secondVirialCoeff),
+            mu_(dynamicViscosity),
+            mug_(vapourDynamicViscosity),
+            K_(thermalConductivity),
+            Kg_(vapourThermalConductivity),
+            sigma_(surfaceTension),
+            D_(vapourDiffussivity)
+        {}
+
+        //- Construct from Istream
+        nC3H8O(Istream& is)
+        :
+            liquid(is),
+            rho_(is),
+            pv_(is),
+            hl_(is),
+            cp_(is),
+            h_(is),
+            cpg_(is),
+            B_(is),
+            mu_(is),
+            mug_(is),
+            K_(is),
+            Kg_(is),
+            sigma_(is),
+            D_(is)
+        {}
+
+
+    // Member Functions
+
+        //- Liquid density [kg/m^3]
+        scalar rho(scalar p, scalar T) const
+        {
+            return rho_.f(p, T);
+        }
+
+        //- Vapour pressure [Pa]
+        scalar pv(scalar p, scalar T) const
+        {
+            return pv_.f(p, T);
+        }
+
+        //- Heat of vapourisation [J/kg]
+        scalar hl(scalar p, scalar T) const
+        {
+            return hl_.f(p, T);
+        }
+
+        //- Liquid heat capacity [J/(kg K)]
+        scalar cp(scalar p, scalar T) const
+        {
+            return cp_.f(p, T);
+        }
+
+        //- Liquid Enthalpy [J/(kg)]
+        scalar h(scalar p, scalar T) const
+        {
+            return h_.f(p, T);
+        }
+
+        //- Ideal gas heat capacity [J/(kg K)]
+        scalar cpg(scalar p, scalar T) const
+        {
+            return cpg_.f(p, T);
+        }
+
+        //- Second Virial Coefficient [m^3/kg]
+        scalar B(scalar p, scalar T) const
+        {
+            return B_.f(p, T);
+        }
+
+        //- Liquid viscosity [Pa s]
+        scalar mu(scalar p, scalar T) const
+        {
+            return mu_.f(p, T);
+        }
+
+        //- Vapour viscosity [Pa s]
+        scalar mug(scalar p, scalar T) const
+        {
+            return mug_.f(p, T);
+        }
+
+        //- Liquid thermal conductivity  [W/(m K)]
+        scalar K(scalar p, scalar T) const
+        {
+            return K_.f(p, T);
+        }
+
+        //- Vapour thermal conductivity  [W/(m K)]
+        scalar Kg(scalar p, scalar T) const
+        {
+            return Kg_.f(p, T);
+        }
+
+        //- Surface tension [N/m]
+        scalar sigma(scalar p, scalar T) const
+        {
+            return sigma_.f(p, T);
+        }
+
+        //- Vapour diffussivity [m2/s]
+        scalar D(scalar p, scalar T) const
+        {
+            return D_.f(p, T);
+        }
+
+
+        //- Write the function coefficients
+        void writeData(Ostream& os) const
+        {
+            liquid::writeData(os); os << nl;
+            rho_.writeData(os); os << nl;
+            pv_.writeData(os); os << nl;
+            hl_.writeData(os); os << nl;
+            cp_.writeData(os); os << nl;
+            cpg_.writeData(os); os << nl;
+            B_.writeData(os); os << nl;
+            mu_.writeData(os); os << nl;
+            mug_.writeData(os); os << nl;
+            K_.writeData(os); os << nl;
+            Kg_.writeData(os); os << nl;
+            sigma_.writeData(os); os << nl;
+            D_.writeData(os); os << endl;
+        }
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream& os, const nC3H8O& l)
+        {
+            l.writeData(os);
+            return os;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/tutorials/dieselFoam/aachenBomb/constant/injectorProperties b/tutorials/dieselFoam/aachenBomb/constant/injectorProperties
index 66e47bddebc7818798a167b1153d09d63a31acd0..81b8b07bc6d0135f84a4c63b2b0b78dbc917dd7d 100644
--- a/tutorials/dieselFoam/aachenBomb/constant/injectorProperties
+++ b/tutorials/dieselFoam/aachenBomb/constant/injectorProperties
@@ -25,7 +25,6 @@ FoamFile
             diameter        0.00019;
             Cd              0.9;
             mass            6e-06;
-            temperature     320;
             nParcels        5000;
 
             X
@@ -67,6 +66,14 @@ FoamFile
                 (0.00120833 5.1737)
                 (0.00125 3.9213)
             );
+
+            temperatureProfile
+            (
+                (0.0      320.0)
+                (0.00125  320.0)
+            );
+
+
         }
 
         commonRailInjectorProps
diff --git a/tutorials/dieselFoam/aachenBomb/constant/polyMesh/boundary b/tutorials/dieselFoam/aachenBomb/constant/polyMesh/boundary
deleted file mode 100644
index 776ce2125360feefef39ba39f7b92e300f44f131..0000000000000000000000000000000000000000
--- a/tutorials/dieselFoam/aachenBomb/constant/polyMesh/boundary
+++ /dev/null
@@ -1,27 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      binary;
-    class       polyBoundaryMesh;
-    object      boundary;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-1
-(
-    walls
-    {
-        type            wall;
-        nFaces          19762;
-        startFace       494419;
-    }
-)
-
-// ************************************************************************* //