diff --git a/hera/simulations/LSM/template.py b/hera/simulations/LSM/template.py
index 66a33eee1..98baa86aa 100644
--- a/hera/simulations/LSM/template.py
+++ b/hera/simulations/LSM/template.py
@@ -293,6 +293,7 @@ def run(self,topography=None, stations=None,canopy=None,params=dict(),deposition
                     finalxarray = xarray.concat(L, dim="datetime")
                     new_coords = dict(x=finalxarray.x - xshift, y=finalxarray.y - yshift)
                     finalxarray = finalxarray.assign_coords(coords=new_coords)
+                    finalxarray = finalxarray.transpose("datetime", "x","y","z")
                     finalxarray.to_netcdf(os.path.join(netcdf_output, "data%s.nc" % i))
                     L = []
                     i += 1
@@ -302,6 +303,7 @@ def run(self,topography=None, stations=None,canopy=None,params=dict(),deposition
 
             new_coords = dict(x=finalxarray.x-xshift,y=finalxarray.y-yshift)
             finalxarray= finalxarray.assign_coords(coords=new_coords)
+            finalxarray = finalxarray.transpose("datetime", "x", "y", "z")
 
             logger.info(f"saved xarray in {netcdf_output}")
             if not self.forceKeep:
diff --git a/hera/simulations/gaussian/Meteorology.py b/hera/simulations/gaussian/Meteorology.py
index 731909acb..8c9e9bd74 100644
--- a/hera/simulations/gaussian/Meteorology.py
+++ b/hera/simulations/gaussian/Meteorology.py
@@ -270,6 +270,13 @@ def getWindVelocity(self, height):
         return self.u_refHeight * (height / refHeight) ** self.wind_p
 
 
+    def getWindVelocity_hotSpot(self, height):
+        facfor_dict = dict(A=0.07, B=0.07, C=0.1, D=0.15, E=0.35, F=0.55)
+        u_H = self.u10*(height/(10*ureg.m))**facfor_dict[self.stability]  # The value 10*m in the denominator is the reference height for u10.
+        # u_H_float = round(u_H.m_as(ureg.m/ureg.s), 2)
+        # ret = u_H_float*ureg.m/ureg.s
+        return u_H
+
 class StandardMeteorolgyConstant_log(StandardMeteorolgyConstant_powerLaw):
 
     def getWindVelocity(self, height):
@@ -333,7 +340,7 @@ def __init__(self):
                                 log=StandardMeteorolgyConstant_log,
                                 uniformWind=StandardMeteorolgyConstant_uniformWind)
 
-    def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=20*ureg.degC, stability="D", z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC):
+    def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=ureg.Quantity(20, ureg.degC), stability="D", z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=ureg.Quantity(35, ureg.degC)):
         """
            Creating a meteorology object.
 
@@ -348,8 +355,8 @@ def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", tempe
         return self.meteorology[verticalProfileType](u10=u10, inversion=inversion, temperature=temperature, stability=stability,
                                                      z0=z0, ustar=ustar, skinSurfaceTemperature=skinSurfaceTemperature)
 
-    def getMeteorologyFromURefHeight(self, u, inversion, refHeight, verticalProfileType="log", temperature=20*ureg.degC, stability="D",
-                                     z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC):
+    def getMeteorologyFromURefHeight(self, u, inversion, refHeight, verticalProfileType="log", temperature=ureg.Quantity(20, ureg.degC), stability="D",
+                                     z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=ureg.Quantity(35, ureg.degC)):
         """
            Creating a meteorology object.
 
diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py
index 9cc1ff92d..37add92ec 100644
--- a/hera/simulations/gaussian/gasCloud.py
+++ b/hera/simulations/gaussian/gasCloud.py
@@ -11,7 +11,7 @@
 
 class abstractGasCloud:
 
-    def __init__(self, sourceQ, sourceHeight, initialCloudSize, sigmaType):
+    def __init__(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind_profile_type, spaceTime, sigmaType):
         """
 
         Parameters
@@ -26,17 +26,27 @@ def __init__(self, sourceQ, sourceHeight, initialCloudSize, sigmaType):
 
         sourceHeight : pint Quantity
         initialCloudSize : 3-touple pint Quantity, the sigmas in each axis.
+        wind_profile_type: str,
+        spaceTime: dictionary containing the keys:
         sigmaType : The sigma type, for example from Briggs, rural/urban.
         """
         self.sourceHeight = sourceHeight
         self.initialCloudSize = initialCloudSize
         self.sigmaType = sigmaType
         self.sourceQ = sourceQ
+        self.meteorology = meteorology
+        self.spaceTime = spaceTime
 
+        if wind_profile_type == 'HotSpot':
+            self.u = self.meteorology.getWindVelocity_hotSpot(height=sourceHeight)
+        elif wind_profile_type == 'default':
+            self.u = self.meteorology.getWindVelocity(height=sourceHeight)
+        else:
+            raise ValueError("wind_profile_type must be either 'default' or 'HotSpot'")
 
 
     @staticmethod
-    def createGasCloud(sourceQ,sourceHeight,initialCloudSize,sigmaType):
+    def createGasCloud(sourceQ,sourceHeight,initialCloudSize,meteorology,wind_profile_type,spaceTime,sigmaType):
         """
             Return the type of the release based on the units of Q
         Parameters
@@ -69,16 +79,16 @@ def createGasCloud(sourceQ,sourceHeight,initialCloudSize,sigmaType):
 
         returnCls = instantaneousReleaseGasCloud if instantaneous else continuousReleaseGasCloud
 
-        return returnCls(sourceQ=sourceQ,sourceHeight=sourceHeight,initialCloudSize=initialCloudSize,sigmaType=sigmaType)
+        return returnCls(sourceQ=sourceQ,sourceHeight=sourceHeight,initialCloudSize=initialCloudSize,
+                         meteorology=meteorology,wind_profile_type=wind_profile_type, spaceTime=spaceTime, sigmaType=sigmaType)
 
 
-    def _getTXterm(self, stability, u, xcoordRange, tcoordRange):
+    def _getTXterm(self, xcoordRange, tcoordRange):
         """
         Parameters
         ----------
         initialCloudSize : 3-tuple of float/pint Quantity (default m)
             The initial cloud size (standard deviation) in the x,y and z dimensions.
-        stability : the stability class
         u: the wind speed / pint Quantity (default m/s)
         xcoordRange : Tuple in numpy.arange format, unitless.
         tcoordRange : Tuple in numpy.arange format, unitless.
@@ -87,22 +97,23 @@ def _getTXterm(self, stability, u, xcoordRange, tcoordRange):
         -------
         The X component of the Gaussian concentration formula.
         """
+
+        u = tonumber(self.u, ureg.m / ureg.min)
         T, X = numpy.meshgrid(tcoordRange, xcoordRange, indexing='ij')
-        sigmaX = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaX']
+        sigmaX = self.sigmaType.getSigma(x=X, stability=self.meteorology.stability, sigma0=self.initialCloudSize, units=False)['sigmaX']
         downwind = (1 / (numpy.sqrt(2 * numpy.pi) * sigmaX)) * numpy.exp((-(X - u * T) ** 2) / (2 * sigmaX ** 2))
         XR_downwind = xarray.DataArray(downwind, dims=("time", "x"), coords={"time": tcoordRange, "x": xcoordRange})
         return XR_downwind
 
 
 
-    def _getXYterm(self, stability, xcoordRange, ycoordRange):
+    def _getXYterm(self, xcoordRange, ycoordRange):
         """
 
         Parameters
         ----------
         initialCloudSize : 3-tuple of float/pint Quantity (default m)
             The initial cloud size (standard deviation) in the x,y and z dimensions.
-        stability : the stability class
         xcoordRange : Tuple in numpy.arange format, unitless.
         ycoordRange : Tuple in numpy.arange format, unitless.
 
@@ -111,7 +122,7 @@ def _getXYterm(self, stability, xcoordRange, ycoordRange):
         The Y component of the Gaussian concentration formula.
         """
         X, Y = numpy.meshgrid(xcoordRange, ycoordRange, indexing='ij')
-        sigmaY = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaY']
+        sigmaY = self.sigmaType.getSigma(x=X, stability=self.meteorology.stability, sigma0=self.initialCloudSize, units=False)['sigmaY']
         crosswind = (1 / (numpy.sqrt(2 * numpy.pi) * sigmaY)) * numpy.exp((-(Y - 0) ** 2) / (2 * sigmaY ** 2))
         XR_crosswind = xarray.DataArray(crosswind, dims=("x", "y"), coords={"x": xcoordRange, "y": ycoordRange})
         return XR_crosswind
@@ -119,15 +130,13 @@ def _getXYterm(self, stability, xcoordRange, ycoordRange):
 
 
 
-    def _getXZterm(self, stability, inversion, xcoordRange, zcoordRange, numOfReflections):
+    def _getXZterm(self, xcoordRange, zcoordRange, numOfReflections):
         """
 
         Parameters
         ----------
         initialCloudSize : 3-tuple of float/pint Quantity (default m)
             The initial cloud size (standard deviation) in the x,y and z dimensions.
-        stability : the stability class
-        inversion: the inversion height / pint Quantity (default m)
         xcoordRange : Tuple in numpy.arange format, unitless.
         zcoordRange : Tuple in numpy.arange format, unitless.
         numOfReflections : The number of reflections of the summation of the Z component.
@@ -136,9 +145,11 @@ def _getXZterm(self, stability, inversion, xcoordRange, zcoordRange, numOfReflec
         -------
         The Z component of the Gaussian concentration formula.
         """
+        inversion = self.meteorology.inversion.m_as(ureg.m)
+
         sourceHeight = tonumber(self.sourceHeight, ureg.m)
         X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij')
-        sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ']
+        sigmaZ = self.sigmaType.getSigma(x=X, stability=self.meteorology.stability, sigma0=self.initialCloudSize, units=False)['sigmaZ']
 
         nSum = numpy.arange(-numOfReflections, numOfReflections + 1, 1)
         vertical = 0
@@ -165,10 +176,10 @@ def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dx
         :return:
         """
 
-        xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy,ureg.m))
-        ycoordRange = numpy.arange(tonumber(miny,ureg.m),tonumber(maxy,ureg.m),tonumber(dxdy,ureg.m))
-        zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz,ureg.m))
-        tcoordRange = numpy.arange(0,tonumber(timeSpan,ureg.min),tonumber(dt,ureg.min))
+        xcoordRange = numpy.arange(minx.m_as(ureg.m), maxx.m_as(ureg.m), dxdy.m_as(ureg.m))
+        ycoordRange = numpy.arange(miny.m_as(ureg.m),maxy.m_as(ureg.m),dxdy.m_as(ureg.m))
+        zcoordRange = numpy.arange(minz.m_as(ureg.m), maxz.m_as(ureg.m), dz.m_as(ureg.m))
+        tcoordRange = numpy.arange(0,timeSpan.m_as(ureg.min),dt.m_as(ureg.min))
 
         frac, X = numpy.meshgrid(fracVector, xcoordRange, indexing='ij')
         XR_downwind = xarray.DataArray(frac, dims=("time", "x"), coords={"time": tcoordRange, "x": xcoordRange})
@@ -180,13 +191,12 @@ def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dx
 
 
 
-    def _getTXDosage(self, stability, u, xcoordRange, tcoordRange):
+    def _getTXDosage(self, xcoordRange, tcoordRange):
         """
         Parameters
         ----------
         initialCloudSize : 3-tuple of float/pint Quantity (default m)
             The initial cloud size (standard deviation) in the x,y and z dimensions.
-        stability : the stability class
         u: the wind speed / pint Quantity (default m/s)
         xcoordRange : Tuple in numpy.arange format, unitless.
         tcoordRange : Tuple in numpy.arange format, unitless.
@@ -195,8 +205,10 @@ def _getTXDosage(self, stability, u, xcoordRange, tcoordRange):
         -------
         The X component of the Gaussian dosage formula.
         """
+
+        u = tonumber(self.u, ureg.m / ureg.min)
         T, X = numpy.meshgrid(tcoordRange, xcoordRange, indexing='ij')
-        sigmaX = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaX']
+        sigmaX = self.sigmaType.getSigma(x=X, stability=self.meteorology.stability, sigma0=self.initialCloudSize, units=False)['sigmaX']
         downwind_erf = (1/(2*u))*(special.erf(X/(numpy.sqrt(2)*sigmaX))-special.erf((X-u*T)/(numpy.sqrt(2)*sigmaX)))
         XR_downwind_erf = xarray.DataArray(downwind_erf, dims=( "time", "x"), coords={"time": tcoordRange, "x": xcoordRange})
         return XR_downwind_erf
@@ -281,41 +293,54 @@ def _getXZterm_ones(self, xcoordRange, zcoordRange):
         return XR_vertical
 
 
-    def _getDF(self, stability, u, xcoordRange, zcoordRange):
 
+    def _getDF(self, xcoordRange, zcoordRange):
         """
-        :param stability: The stability state
         :param u: Wind speed
         :param xcoordRange: x-coordinates of the grid
         :param zcoordRange: z-coordinates of the grid
         :return: The Depletion Factor
         """
-
+        # Deposition velocity normalized to standard m/s units
         v = tonumber(0.003*ureg.m/ureg.s, ureg.m/ureg.min) #Deposition velocity. By default we take this value to be 0.003 [m/s]
+        u = tonumber(self.u, ureg.m/ureg.min)
 
         X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij')
-        sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ']
+        sigmaZ = self.sigmaType.getSigma(x=X, stability=self.meteorology.stability, sigma0=self.initialCloudSize, units=False)['sigmaZ']
         H = tonumber(self.sourceHeight, ureg.m)
+
+        # Grid spacing along the downwind x-axis
         dx = xcoordRange[1] - xcoordRange[0]
-        DF = (numpy.e**(numpy.cumsum(1/(sigmaZ*numpy.e**(0.5*(H/sigmaZ)**2))*dx, axis=0)))**((-v/u)*numpy.sqrt(2/numpy.pi))
+
+        # 1. Calculate the inner component of the integral
+        inner_term = 1 / (sigmaZ * numpy.exp(0.5 * (H / sigmaZ) ** 2)) * dx
+
+        # 2. Integrate along the downwind distance axis (x-axis corresponds to axis=0)
+        exponent_sum = numpy.cumsum(inner_term, axis=0)
+
+        # 3. Apply the constants and the negative multiplier in log-space
+        outer_multiplier = (-v / u) * numpy.sqrt(2 / numpy.pi)
+        total_exponent = exponent_sum * outer_multiplier
+
+        # 4. Bring back to standard space safely
+        DF = numpy.exp(total_exponent)
 
         return DF
 
 
 
-    def getDF_noQ_xarray(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min):
-        stability = meteorology.stability
-        u = tonumber(meteorology.u10, ureg.m/ureg.min)
+    def getDF_noQ_xarray(self):
 
-        xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy, ureg.m))
-        ycoordRange = numpy.arange(tonumber(miny, ureg.m), tonumber(maxy, ureg.m), tonumber(dxdy, ureg.m))
-        zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz, ureg.m))
-        tcoordRange = numpy.arange(0, tonumber(timeSpan, ureg.min), tonumber(dt, ureg.min))
+        ST = self.spaceTime
+        xcoordRange = numpy.arange(tonumber(ST['minx'], ureg.m), tonumber(ST['maxx'], ureg.m), tonumber(ST['dxdy'], ureg.m))
+        ycoordRange = numpy.arange(tonumber(ST['miny'], ureg.m), tonumber(ST['maxy'], ureg.m), tonumber(ST['dxdy'], ureg.m))
+        zcoordRange = numpy.arange(tonumber(ST['minz'], ureg.m), tonumber(ST['maxz'], ureg.m), tonumber(ST['dz'], ureg.m))
+        tcoordRange = numpy.arange(0, tonumber(ST['timeSpan'], ureg.min), tonumber(ST['dt'], ureg.min))
 
         TX = self._getTXterm_ones(xcoordRange=xcoordRange, tcoordRange=tcoordRange)
         XY = self._getXYterm_ones(xcoordRange=xcoordRange, ycoordRange=ycoordRange)
         XZ = self._getXZterm_ones(xcoordRange=xcoordRange, zcoordRange=zcoordRange)
-        DF = self._getDF(stability=stability, u=u, xcoordRange=xcoordRange, zcoordRange=zcoordRange)
+        DF = self._getDF(xcoordRange=xcoordRange, zcoordRange=zcoordRange)
 
         return TX*XY*(XZ*DF)
 
@@ -325,22 +350,17 @@ def getDF_noQ_xarray(self, meteorology, minx, miny, minz, maxx, maxy, maxz, time
 class instantaneousReleaseGasCloud(abstractGasCloud):
 
 
-    def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                        dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3):
+    def getConcentrationFromMinMaxRange_inst_noQ(self, numOfReflections=3):
 
-        xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy,ureg.m))
-        ycoordRange = numpy.arange(tonumber(miny,ureg.m),tonumber(maxy,ureg.m),tonumber(dxdy,ureg.m))
-        zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz,ureg.m))
-        tcoordRange = numpy.arange(0,tonumber(timeSpan,ureg.min),tonumber(dt,ureg.min))
+        ST = self.spaceTime
+        xcoordRange = numpy.arange(tonumber(ST['minx'], ureg.m), tonumber(ST['maxx'], ureg.m), tonumber(ST['dxdy'], ureg.m))
+        ycoordRange = numpy.arange(tonumber(ST['miny'], ureg.m), tonumber(ST['maxy'], ureg.m), tonumber(ST['dxdy'], ureg.m))
+        zcoordRange = numpy.arange(tonumber(ST['minz'], ureg.m), tonumber(ST['maxz'], ureg.m), tonumber(ST['dz'], ureg.m))
+        tcoordRange = numpy.arange(0, tonumber(ST['timeSpan'], ureg.min), tonumber(ST['dt'], ureg.min))
 
-        stability = meteorology.stability
-        u = tonumber(meteorology.u10, ureg.m/ureg.min)
-        inversion = tonumber(meteorology.inversion, ureg.m)
-
-        TX = self._getTXterm(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange)
-        XY = self._getXYterm(stability=stability, xcoordRange=xcoordRange, ycoordRange=ycoordRange)
-        XZ = self._getXZterm(stability=stability, inversion=inversion, xcoordRange=xcoordRange, zcoordRange=zcoordRange,
-                             numOfReflections=numOfReflections)
+        TX = self._getTXterm(xcoordRange=xcoordRange, tcoordRange=tcoordRange)
+        XY = self._getXYterm(xcoordRange=xcoordRange, ycoordRange=ycoordRange)
+        XZ = self._getXZterm(xcoordRange=xcoordRange, zcoordRange=zcoordRange, numOfReflections=numOfReflections)
 
         ret = TX*XY*XZ
         ret.attrs['Q'] = 1/ureg.m**3
@@ -348,41 +368,32 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz
         return ret
 
 
-    def getDosageFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                           dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
-        stability = meteorology.stability
-        u = tonumber(meteorology.u10, ureg.m/ureg.min)
-        inversion = tonumber(meteorology.inversion, ureg.m)
-
-        xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy, ureg.m))
-        ycoordRange = numpy.arange(tonumber(miny, ureg.m), tonumber(maxy, ureg.m), tonumber(dxdy, ureg.m))
-        zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz, ureg.m))
-        tcoordRange = numpy.arange(0, tonumber(timeSpan, ureg.min), tonumber(dt, ureg.min))
+    def getDosageFromMinMaxRange_inst_noQ(self, numOfReflections=3, DF=False):
+        ST = self.spaceTime
+        xcoordRange = numpy.arange(tonumber(ST['minx'], ureg.m), tonumber(ST['maxx'], ureg.m), tonumber(ST['dxdy'], ureg.m))
+        ycoordRange = numpy.arange(tonumber(ST['miny'], ureg.m), tonumber(ST['maxy'], ureg.m), tonumber(ST['dxdy'], ureg.m))
+        zcoordRange = numpy.arange(tonumber(ST['minz'], ureg.m), tonumber(ST['maxz'], ureg.m), tonumber(ST['dz'], ureg.m))
+        tcoordRange = numpy.arange(0, tonumber(ST['timeSpan'], ureg.min), tonumber(ST['dt'], ureg.min))
 
-        TX = self._getTXDosage(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange)
-        XY = self._getXYterm(stability=meteorology.stability, xcoordRange=xcoordRange, ycoordRange=ycoordRange)
-        XZ = self._getXZterm(stability=stability, inversion=inversion, xcoordRange=xcoordRange, zcoordRange=zcoordRange,
-                             numOfReflections=numOfReflections)
+        TX = self._getTXDosage(xcoordRange=xcoordRange, tcoordRange=tcoordRange)
+        XY = self._getXYterm(xcoordRange=xcoordRange, ycoordRange=ycoordRange)
+        XZ = self._getXZterm(xcoordRange=xcoordRange, zcoordRange=zcoordRange, numOfReflections=numOfReflections)
         D_F = 1
         if DF:
-            D_F = self._getDF(stability=stability, u=u, xcoordRange=xcoordRange, zcoordRange=zcoordRange)
+            D_F = self._getDF(xcoordRange=xcoordRange, zcoordRange=zcoordRange)
 
         ret = TX*XY*(XZ*D_F)
         ret.attrs['Q'] = 1*ureg.min/ureg.m**3
 
         return ret
 
-    def getDosageFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                           dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
-        C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                                    maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan,dxdy=dxdy,
-                                                                    dz=dz, dt=dt,numOfReflections=numOfReflections)
+    def getDosageFromMinMaxRange_inst_NoERF_noQ(self, numOfReflections=3, DF=False):
+        C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(numOfReflections=numOfReflections)
 
         D_without_Q = self.trapezoidal_integration(data=C_without_Q)
         D_F = 1
         if DF:
-            D_F = self.getDF_noQ_xarray(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy,
-                                        maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt)
+            D_F = self.getDF_noQ_xarray()
 
         ret = D_without_Q*D_F
         ret.attrs['Q'] = 1*ureg.min/ureg.m**3
@@ -390,48 +401,16 @@ def getDosageFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz,
         return ret
 
 
-class instantaneousReleaseGasCloud(abstractGasCloud):
-
-
-    def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                        dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3):
-
-        xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy,m))
-        ycoordRange = numpy.arange(tonumber(miny,m),tonumber(maxy,m),tonumber(dxdy,m))
-        zcoordRange = numpy.arange(tonumber(minz, m), tonumber(maxz, m), tonumber(dz,m))
-        tcoordRange = numpy.arange(0,tonumber(timeSpan,min),tonumber(dt,min))
-
-        stability = meteorology.stability
-        u = tonumber(meteorology.u10, m/min)
-        inversion = tonumber(meteorology.inversion, m)
-
-        TX = self._getTXterm(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange)
-        XY = self._getXYterm(stability=stability, xcoordRange=xcoordRange, ycoordRange=ycoordRange)
-        XZ = self._getXZterm(stability=stability, inversion=inversion, xcoordRange=xcoordRange, zcoordRange=zcoordRange,
-                             numOfReflections=numOfReflections)
-
-        ret = TX*XY*XZ
-        ret.attrs['Q'] = 1/m**3
-
-        return ret
-
-
-    def getConcentrationFromMinMaxRange_inst(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                        dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
-        C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                                    maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan,dxdy=dxdy,
-                                                                    dz=dz, dt=dt,numOfReflections=numOfReflections)
+    def getConcentrationFromMinMaxRange_inst(self, numOfReflections=3):
+        C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(numOfReflections=numOfReflections)
 
         ret = tonumber(self.sourceQ, ureg.mg)*C_without_Q
         ret.attrs['Q'] = 1*ureg.mg/ureg.m**3
         return ret
 
 
-    def getDosageFromMinMaxRange_inst(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                           dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
-        D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx,
-                                                             maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz,
-                                                             dt=dt, numOfReflections=numOfReflections, DF=DF)
+    def getDosageFromMinMaxRange_inst(self, numOfReflections=3, DF=False):
+        D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(numOfReflections=numOfReflections, DF=DF)
 
         ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q
         ret.attrs['Q'] = 1*ureg.mg*ureg.min/ureg.m**3
@@ -439,12 +418,9 @@ def getDosageFromMinMaxRange_inst(self, meteorology, minx, miny, minz, maxx, max
 
 
 
-    def getDosageFromMinMaxRange_inst_NoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                           dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
+    def getDosageFromMinMaxRange_inst_NoERF(self, numOfReflections=3, DF=False):
 
-        D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                                   maxx=maxx, maxy=maxy,maxz=maxz, timeSpan=timeSpan,dxdy=dxdy,
-                                                                   dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF)
+        D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(**self.spaceTime, numOfReflections=numOfReflections, DF=DF)
 
         ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q
         ret.attrs['Q'] = 1*ureg.mg*ureg.min/ureg.m**3
@@ -463,12 +439,9 @@ def concentrationConversion_mass_to_Bq(self, C, outputUnits, specificActivity):
         return C_Bq
 
 
-    def getTIACFromMinMaxRange_inst(self,specifitActivity, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                           dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False):
+    def getTIACFromMinMaxRange_inst(self,specifitActivity, numOfReflections=3, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False):
 
-        D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx,
-                                                             maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz,
-                                                             dt=dt, numOfReflections=numOfReflections, DF=DF)
+        D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(numOfReflections=numOfReflections, DF=DF)
 
         out_units = unumToPint(outputUnits)
         ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q
@@ -478,12 +451,9 @@ def getTIACFromMinMaxRange_inst(self,specifitActivity, meteorology, minx, miny,
         return ret
 
 
-    def getTIACFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                           dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.s/ureg.m**3, DF=False):
+    def getTIACFromMinMaxRange_inst_noQ(self, numOfReflections=3, outputUnits=ureg.s/ureg.m**3, DF=False):
 
-        D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx,
-                                                             maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz,
-                                                             dt=dt, numOfReflections=numOfReflections, DF=DF)
+        D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(numOfReflections=numOfReflections, DF=DF)
         out_units = unumToPint(outputUnits)
         currentUnites = unumToPint(D_without_Q.attrs['Q'])
         factor = currentUnites.m_as(out_units)
@@ -493,12 +463,10 @@ def getTIACFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, m
 
 
 
-    def getTIACFromMinMaxRange_inst_NoERF(self, specifitActivity, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                    dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False):
+    def getTIACFromMinMaxRange_inst_NoERF(self, specifitActivity, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                          numOfReflections=3, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False):
 
-        D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                             maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy,
-                                                             dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF)
+        D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(**self.spaceTime, numOfReflections=numOfReflections, DF=DF)
 
         out_units = unumToPint(outputUnits)
         ret = tonumber(self.sourceQ, ureg.mg) * D_without_Q
@@ -508,12 +476,10 @@ def getTIACFromMinMaxRange_inst_NoERF(self, specifitActivity, meteorology, minx,
         return ret
 
 
-    def getTIACFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                    dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.s/ureg.m**3, DF=False):
+    def getTIACFromMinMaxRange_inst_NoERF_noQ(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                              numOfReflections=3, outputUnits=ureg.s/ureg.m**3, DF=False):
 
-        D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                             maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy,
-                                                             dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF)
+        D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(**self.spaceTime, numOfReflections=numOfReflections, DF=DF)
 
         out_units = unumToPint(outputUnits)
         currentUnites = unumToPint(D_without_Q.attrs['Q'])
@@ -525,8 +491,8 @@ def getTIACFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, m
 
 
 
-    def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                    dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False):
+    def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                            outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False):
         """
 
         :param C: xarray of concentrations in units of [mass/volume].
@@ -545,8 +511,7 @@ def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, meteorology,
 
         D_F = 1
         if DF:
-            D_F = self.getDF_noQ_xarray(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy,
-                                        maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt)
+            D_F = self.getDF_noQ_xarray(**self.spaceTime)
 
         D_without_Q = D_without_Q * D_F
 
@@ -557,8 +522,8 @@ def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, meteorology,
         return ret
 
 
-    def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                    dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, outputUnits=ureg.s/ureg.m**3, DF=False):
+    def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                                outputUnits=ureg.s/ureg.m**3, DF=False):
         """
 
         :param C_noQ: xarray of concentrations in units of [1/volume].
@@ -569,8 +534,7 @@ def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, meteorology, minx, miny
 
         D_F = 1
         if DF:
-            D_F = self.getDF_noQ_xarray(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy,
-                                        maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt)
+            D_F = self.getDF_noQ_xarray(**self.spaceTime)
 
         D_without_Q = D_without_Q*D_F
         D_without_Q.attrs['Q'] = 1*ureg.min/ureg.m**3 #need to verify that C_noQ was generated with time steps in [min], not [s]
@@ -596,8 +560,8 @@ def get_TIAC_from_dist(self, data, y, z, dist_list):
 
         """
         distances = numpy.array(data.squeeze().x)
-        TIAC_lastTime_DF = numpy.array(data.sel(y=y, z=z, time=data.time[-1], method='nearest'))
-        tuples = list(tuple(zip(distances, TIAC_lastTime_DF)))
+        TIAC_lastTime = numpy.array(data.sel(y=y, z=z, time=data.time[-1], method='nearest'))
+        tuples = list(tuple(zip(distances, TIAC_lastTime)))
 
         for x in dist_list:
             if x not in distances:
@@ -608,10 +572,10 @@ def get_TIAC_from_dist(self, data, y, z, dist_list):
 
 
 
-class continuousReleaseGasCloud(instantaneousReleaseGasCloud):
+class continuousReleaseGasCloud(abstractGasCloud):
 
-    def getConcentrationFromMinMaxRange_cont(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                        dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
+    def getConcentrationFromMinMaxRange_cont(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                             numOfReflections=3, DF=False):
         """
         Returns
         -------
@@ -619,15 +583,13 @@ def getConcentrationFromMinMaxRange_cont(self, meteorology, minx, miny, minz, ma
         since we assume the release rate is constant.
         Here we take the concentration xarray that was claculated using the error function (erf).
         """
-        C_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx,
-                                                        maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt,
-                                                        numOfReflections=numOfReflections, DF=DF)
+        C_without_Q = self.getDosageFromMinMaxRange_inst_noQ(**self.spaceTime,numOfReflections=numOfReflections, DF=DF)
 
         return tonumber(self.sourceQ, ureg.mg/ureg.s)*C_without_Q
 
 
-    def getConcentrationFromMinMaxRange_cont_NoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                        dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
+    def getConcentrationFromMinMaxRange_cont_NoERF(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                                   numOfReflections=3, DF=False):
         """
         Returns
         -------
@@ -635,28 +597,22 @@ def getConcentrationFromMinMaxRange_cont_NoERF(self, meteorology, minx, miny, mi
         since we assume the release rate is constant.
         Here we take the concentration xarray that was claculated without the error function (erf).
         """
-        C_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                          maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy,
-                                                          dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF)
+        C_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(**self.spaceTime, numOfReflections=numOfReflections, DF=DF)
         return tonumber(self.sourceQ, ureg.mg/ureg.s)*C_without_Q
 
 
-    def getDosageFromMinMaxRange_cont_NoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                        dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
+    def getDosageFromMinMaxRange_cont_NoERF(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                            numOfReflections=3, DF=False):
 
-        C_without_Q = self.getConcentrationFromMinMaxRange_cont(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                           maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy,
-                                                           dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF)
+        C_without_Q = self.getConcentrationFromMinMaxRange_cont(**self.spaceTime, numOfReflections=numOfReflections, DF=DF)
         D_without_Q = self.trapezoidal_integration(data=C_without_Q)
 
         return tonumber(self.sourceQ, ureg.mg/ureg.min) * D_without_Q
 
-    def getDosageFromMinMaxRange_cont_doubleNoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan,
-                                        dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False):
+    def getDosageFromMinMaxRange_cont_doubleNoERF(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt,
+                                                  numOfReflections=3, DF=False):
 
-        C_without_Q = self.getConcentrationFromMinMaxRange_cont_NoERF(meteorology=meteorology, minx=minx, miny=miny, minz=minz,
-                                                           maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy,
-                                                           dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF)
+        C_without_Q = self.getConcentrationFromMinMaxRange_cont_NoERF(**self.spaceTime, numOfReflections=numOfReflections, DF=DF)
         D_without_Q = self.trapezoidal_integration(data=C_without_Q)
 
         return tonumber(self.sourceQ, ureg.mg/ureg.min) * D_without_Q
diff --git a/hera/simulations/gaussian/toolkit.py b/hera/simulations/gaussian/toolkit.py
index bdfb6e5a5..d98fba53e 100644
--- a/hera/simulations/gaussian/toolkit.py
+++ b/hera/simulations/gaussian/toolkit.py
@@ -4,6 +4,8 @@
 from hera.simulations.gaussian.Meteorology import MeteorologyFactory
 from hera.utils.unitHandler import ureg, unumToPint, celsius, K
 import matplotlib.pyplot as plt
+import scipy
+import numpy
 
 
 
@@ -25,6 +27,7 @@ def __init__(self, projectName: str, filesDirectory: str = None):
 
         self._presentation = presentationLayer()
 
+
     def getSigmaType(self,sigmaName):
         """
 
@@ -56,22 +59,73 @@ def listSigmaTypes(self):
 
 
 
-
-
-    def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=20*ureg.degC, stability="D",
-                              z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC):
+    def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=ureg.Quantity(20, ureg.degC), stability="D",
+                              z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=ureg.Quantity(35, ureg.degC)):
         return MeteorologyFactory().getMeteorologyFromU10(u10=u10, inversion=inversion, verticalProfileType=verticalProfileType,
                     temperature=temperature, stability=stability, z0=z0, ustar=ustar, skinSurfaceTemperature=skinSurfaceTemperature)
 
 
-    def getMeteorologyFromURefHeight(self, u, refHeight, inversion, verticalProfileType="log", temperature=20*ureg.degC, stability="D",
-                              z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC):
+    def getMeteorologyFromURefHeight(self, u, refHeight, inversion, verticalProfileType="log", temperature=ureg.Quantity(20, ureg.degC), stability="D",
+                              z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=ureg.Quantity(35, ureg.degC)):
         return MeteorologyFactory().getMeteorologyFromURefHeight(u=u, refHeight=refHeight,  inversion=inversion,
                     verticalProfileType=verticalProfileType, temperature=temperature, stability=stability, z0=z0,
                     ustar=ustar,skinSurfaceTemperature=skinSurfaceTemperature)
 
 
-    def getGasCloud(self, sourceQ, sourceHeight, initialCloudSize, sigmaTypeName="briggsRural"):
+
+
+
+    def getSpaceTime(self, meteorology, sourceHeight, wind_profile_type, maxx, dt, dz, dxdy_multiplier, minimal_maxy, initialCloudSize):
+
+        # setting the wind speed at the source height according to the Hot Spot's manual, or Hera's built-in functions.
+        if wind_profile_type == 'HotSpot':
+            windSpeed = meteorology.getWindVelocity_hotSpot(height=sourceHeight)
+        elif wind_profile_type == 'default':
+            windSpeed = meteorology.getWindVelocity(height=sourceHeight)
+        else:
+            raise ValueError("wind_profile_type must be either 'default' or 'HotSpot'")
+
+        # it is preferable to take dxdy as a multiple of dt*windSpeed, so that the cloud center is directly above a grid-point.
+        # the dxdy_multiplier makes the cloud above every dxdy_multiplier-th grid-point (on the datetime axis).
+        dxdy = dt.m_as(ureg.s) * windSpeed.m_as(ureg.m / ureg.s) * dxdy_multiplier * ureg.m
+
+        # the Y-span should be the smallest multiple of dxdy that is greater/equal to the initially given maxy.
+        # the purpose is to limit the grid along the y-axis, to reduce runtime.
+        maxy = numpy.ceil(minimal_maxy.m_as(ureg.m) / dxdy.m_as(ureg.m)) * dxdy
+
+        # the time span should be such that the cloud has passed far enough beyond the maximum X (maxx),
+        # yet not too far as to calculate for times that don't affect the concentraion up to maxx.
+        # therefore we take a time span for which maxx is 3*sigmas away from the cloud center:
+        def find_x_for_timeSpan(maxx):
+            maxx = maxx.m_as(ureg.m)
+
+            def get_function(x, maxx):
+                sigmaX = self.getSigmaType(sigmaName='briggsRural').getSigma(x=x, stability=meteorology.stability,
+                                                                           sigma0=initialCloudSize,
+                                                                           units=False)['sigmaX'][0]
+                return x - 3 * sigmaX - maxx
+
+            root = scipy.optimize.root_scalar(get_function, bracket=[0, 10 ** 6], args=maxx).root
+            return root
+
+        x_timeSpan = find_x_for_timeSpan(maxx) * ureg.m
+        timeSpan = x_timeSpan / (windSpeed * 60 * (ureg.s / ureg.min))
+        timeSpan = numpy.ceil(timeSpan.m_as(ureg.min)) * ureg.min
+
+        spaceTime = {
+            'minx': 0 * ureg.m, 'maxx': maxx,
+            'miny': -maxy, 'maxy': maxy + 1 * ureg.m,
+            'minz': 0 * ureg.m, 'maxz': 2 * ureg.m,
+            'dxdy': dxdy, 'dz': dz,
+            'timeSpan': timeSpan, 'dt': dt}
+
+        return spaceTime
+
+
+
+
+    def getGasCloud(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind_profile_type,
+                    spaceTime, sigmaTypeName="briggsRural"):
         """
 
         Parameters
@@ -86,6 +140,7 @@ def getGasCloud(self, sourceQ, sourceHeight, initialCloudSize, sigmaTypeName="br
 
         sourceHeight : unum
         initialCloudSize : 3-touple unum, the sigmas in each axis.
+        wind_profile_type : str, gets either "HotSpot" or "default".
         sigmaTypeName : Name of the sigma type, for example from Briggs, rural/urban.
 
         Returns
@@ -94,7 +149,10 @@ def getGasCloud(self, sourceQ, sourceHeight, initialCloudSize, sigmaTypeName="br
 
         """
         sigmaType = self.getSigmaType(sigmaTypeName)
-        gascloud = abstractGasCloud.createGasCloud(sourceQ=sourceQ,sourceHeight=sourceHeight,initialCloudSize=initialCloudSize,sigmaType=sigmaType)
+        gascloud = abstractGasCloud.createGasCloud(sourceQ=sourceQ,sourceHeight=sourceHeight,
+                                                   initialCloudSize=initialCloudSize,meteorology=meteorology,
+                                                   wind_profile_type=wind_profile_type, spaceTime=spaceTime,
+                                                   sigmaType=sigmaType)
         return gascloud
 
 
@@ -227,7 +285,7 @@ def plotFixedPointDosageOverTime(self, D, x, y, z):
 
         plt.plot(time_array, dos_inst_t )
         plt.xlabel("Time from release $[min]$")
-        plt.ylabel(f"Dosage over time {units}")
+        plt.ylabel(f"TIAC over time {units}")
         plt.title(f"Receptor at x={x}[m], y={y}[m], z={z}[m].")
         plt.grid()
         plt.show()
@@ -256,10 +314,12 @@ def plotMaxConcentrationOverTime_noQ(self,C, y, z, x_min=None,x_max=None):
         plt.ylabel(r"Concentration $\left[\frac{1}{m^3}\right]$")
         plt.title(f"Maximum concentration over time. y={y}[m], z={z}[m]")
         plt.grid()
-        if x_min is not None and x_max is not None:
-            x_min = unumToPint(x_min).m_as(ureg.m)
-            x_max = unumToPint(x_max).m_as(ureg.m)
-            plt.xlim(x_min, x_max)
+
+        x_min = (C.x[0].values)*ureg.meter if x_min is None else x_min
+        x_max = (C.x[-1].values)*ureg.meter if x_max is None else x_max
+        x_min = unumToPint(x_min).m_as(ureg.m)
+        x_max = unumToPint(x_max).m_as(ureg.m)
+        plt.xlim(x_min, x_max)
         plt.show()
 
 
@@ -285,15 +345,16 @@ def plotFixedPointConcentrationOverTime_noQ(self, C, x, y, z, t_min=None,t_max=N
         plt.ylabel(r"Concentration $\left[\frac{1}{m^3}\right]$")
         plt.title(f"Receptor at x={x}[m], y={y}[m], z={z}[m].")
         plt.grid()
-        if t_min is not None and t_max is not None:
-            t_min = unumToPint(t_min).m_as(ureg.min)
-            t_max = unumToPint(t_max).m_as(ureg.min)
-            plt.xlim(t_min, t_max)
+        t_min = (C.time[0].values)*ureg.min if t_min is None else t_min
+        t_max = (C.time[-1].values)*ureg.min if t_max is None else t_max
+        t_min = unumToPint(t_min).m_as(ureg.min)
+        t_max = unumToPint(t_max).m_as(ureg.min)
+        plt.xlim(t_min, t_max)
         plt.show()
 
 
 
-    def plotTIACPerDistance_noQ(self,TIAC, y, z, time):
+    def plotTIACPerDistance_noQ(self,TIAC, y, z, time, x_min=None,x_max=None):
         """
         Given a line along the X-axis, this function plots the dosage over distance.
         :param TIAC: xarray of the TIAC in units of [1*time/volume]
@@ -310,13 +371,18 @@ def plotTIACPerDistance_noQ(self,TIAC, y, z, time):
         plt.plot(x_array, dos_x_inst )
         plt.xlabel("Distance from source $[m]$")
 
-        plt.ylabel(r"Dosage $\left[\frac{1*min}{m^3}\right]$")
-        plt.title(f"Dosage per distance. y={y}[m], z={z}[m], time={time}[min]")
+        plt.ylabel(r"TIAC $\left[\frac{1}{m^3} \cdot min\right]$")
+        plt.title(f"TIAC per distance. y={y}[m], z={z}[m], time={time}[min]")
         plt.grid()
+        x_min = (TIAC.x[0].values)*ureg.meter if x_min is None else x_min
+        x_max = (TIAC.x[-1].values)*ureg.meter if x_max is None else x_max
+        x_min = unumToPint(x_min).m_as(ureg.m)
+        x_max = unumToPint(x_max).m_as(ureg.m)
+        plt.xlim(x_min, x_max)
         plt.show()
 
 
-    def plotFixedPointTIACOverTime_noQ(self, TIAC, x, y, z):
+    def plotFixedPointTIACOverTime_noQ(self, TIAC, x, y, z, t_min=None,t_max=None):
         """
         Given a point in space, this function plots the dosage as a function of time.
         :param TIAC: xarray of the TIAC in units of [1*time/volume]
@@ -335,9 +401,15 @@ def plotFixedPointTIACOverTime_noQ(self, TIAC, x, y, z):
 
         plt.plot(time_array, dos_inst_t )
         plt.xlabel("Time from release $[min]$")
-        plt.ylabel(r"Dosage over time $\left[\frac{1*min}{m^3}\right]$")
+        plt.ylabel(r"TIAC over time $\left[\frac{1}{m^3} \cdot min\right]$")
         plt.title(f"Receptor at x={x}[m], y={y}[m], z={z}[m].")
         plt.grid()
+
+        t_min = (TIAC.time[0].values)*ureg.min if t_min is None else t_min
+        t_max = (TIAC.time[-1].values)*ureg.min if t_max is None else t_max
+        t_min = unumToPint(t_min).m_as(ureg.min)
+        t_max = unumToPint(t_max).m_as(ureg.min)
+        plt.xlim(t_min, t_max)
         plt.show()