From ffafb640d74095fe83460ae599789994cf8ea1ea Mon Sep 17 00:00:00 2001 From: jmilo Date: Sun, 7 Jun 2026 18:00:48 +0300 Subject: [PATCH 1/9] Upgrading the gaussian code to adjust the wind speed according to the source height --- hera/simulations/gaussian/gasCloud.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py index 9a9efee8b..e0fe3ed7d 100644 --- a/hera/simulations/gaussian/gasCloud.py +++ b/hera/simulations/gaussian/gasCloud.py @@ -304,7 +304,9 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): def getDF_noQ_xarray(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*m, dz=1*m, dt=1*min): stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) + # u = tonumber(meteorology.u10, m/s) + u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) + u = round(u, 2) xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy, m)) ycoordRange = numpy.arange(tonumber(miny, m), tonumber(maxy, m), tonumber(dxdy, m)) @@ -333,7 +335,9 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz tcoordRange = numpy.arange(0,tonumber(timeSpan,min),tonumber(dt,min)) stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) + # u = tonumber(meteorology.u10, m/s) + u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) + u = round(u, 2) inversion = tonumber(meteorology.inversion, m) TX = self._getTXterm(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) @@ -350,7 +354,9 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz def getDosageFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) + # u = tonumber(meteorology.u10, m/s) + u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) + u = round(u, 2) inversion = tonumber(meteorology.inversion, m) xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy, m)) From c5265bdcb622fcf1dc387d526ebbd411389ebbf3 Mon Sep 17 00:00:00 2001 From: jmilo Date: Thu, 18 Jun 2026 16:29:41 +0300 Subject: [PATCH 2/9] fixing bug with DF --- hera/simulations/gaussian/gasCloud.py | 59 +++++++++++++++++++++------ 1 file changed, 46 insertions(+), 13 deletions(-) diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py index e0fe3ed7d..a62006fa4 100644 --- a/hera/simulations/gaussian/gasCloud.py +++ b/hera/simulations/gaussian/gasCloud.py @@ -280,8 +280,28 @@ def _getXZterm_ones(self, xcoordRange, zcoordRange): return XR_vertical - def _getDF(self, stability, u, xcoordRange, zcoordRange): + # def _getDF(self, stability, u, 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 + # """ + # + # v = tonumber(0.003*m/s, m/min) #Deposition velocity. By default we take this value to be 0.003 [m/s] + # v = tonumber(0.003 * m / s, m/s) + # u_ms = tonumber(u, m / s) + # X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') + # sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] + # H = tonumber(self.sourceHeight, m) + # dx = xcoordRange[1] - xcoordRange[0] + # DF = (numpy.e**(numpy.cumsum(1/(sigmaZ*numpy.e**(0.5*(H/sigmaZ)**2))*dx, axis=0)))**((-v/u_ms)*numpy.sqrt(2/numpy.pi)) + # + # return DF + def _getDF(self, stability, u, xcoordRange, zcoordRange): """ :param stability: The stability state :param u: Wind speed @@ -289,14 +309,30 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): :param zcoordRange: z-coordinates of the grid :return: The Depletion Factor """ - - v = tonumber(0.003*m/s, m/min) #Deposition velocity. By default we take this value to be 0.003 [m/s] + # Deposition velocity normalized to standard m/s units + v = tonumber(0.003 * m / s, m / s) #Deposition velocity. By default we take this value to be 0.003 [m/s] + u_ms = tonumber(u, m / s) + print(u_ms) X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] H = tonumber(self.sourceHeight, 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_ms) * 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 @@ -304,9 +340,8 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): def getDF_noQ_xarray(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*m, dz=1*m, dt=1*min): stability = meteorology.stability - # u = tonumber(meteorology.u10, m/s) - u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) - u = round(u, 2) + u = meteorology.u10 + # u = meteorology.getWindVelocity(height=self.sourceHeight) xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy, m)) ycoordRange = numpy.arange(tonumber(miny, m), tonumber(maxy, m), tonumber(dxdy, m)) @@ -335,9 +370,8 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz tcoordRange = numpy.arange(0,tonumber(timeSpan,min),tonumber(dt,min)) stability = meteorology.stability - # u = tonumber(meteorology.u10, m/s) - u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) - u = round(u, 2) + u = tonumber(meteorology.u10, m/min) + # u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) inversion = tonumber(meteorology.inversion, m) TX = self._getTXterm(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) @@ -354,9 +388,8 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz def getDosageFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): stability = meteorology.stability - # u = tonumber(meteorology.u10, m/s) - u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) - u = round(u, 2) + u = tonumber(meteorology.u10, m/min) + # u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) inversion = tonumber(meteorology.inversion, m) xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy, m)) From 99b9fef97825ad8a68bc822eba296c9bbee4b170 Mon Sep 17 00:00:00 2001 From: jmilo Date: Thu, 18 Jun 2026 16:31:45 +0300 Subject: [PATCH 3/9] fixing bug with DF --- hera/simulations/gaussian/gasCloud.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py index a62006fa4..e9cb236a3 100644 --- a/hera/simulations/gaussian/gasCloud.py +++ b/hera/simulations/gaussian/gasCloud.py @@ -163,7 +163,7 @@ def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dx Note - This vector should have the same length as tcoordRange (described within the function) :return: """ - + 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)) From 8ee3430adf7a3b0092217464188e1a0bda956a40 Mon Sep 17 00:00:00 2001 From: jmilo Date: Thu, 18 Jun 2026 17:08:51 +0300 Subject: [PATCH 4/9] fixing bug with DF --- hera/simulations/gaussian/gasCloud.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py index e9cb236a3..ba40e67fd 100644 --- a/hera/simulations/gaussian/gasCloud.py +++ b/hera/simulations/gaussian/gasCloud.py @@ -310,9 +310,8 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): :return: The Depletion Factor """ # Deposition velocity normalized to standard m/s units - v = tonumber(0.003 * m / s, m / s) #Deposition velocity. By default we take this value to be 0.003 [m/s] - u_ms = tonumber(u, m / s) - print(u_ms) + v = tonumber(0.003 * m / s, m / min) #Deposition velocity. By default we take this value to be 0.003 [m/s] + u = tonumber(u, m / min) X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] @@ -328,7 +327,7 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): exponent_sum = numpy.cumsum(inner_term, axis=0) # 3. Apply the constants and the negative multiplier in log-space - outer_multiplier = (-v / u_ms) * numpy.sqrt(2 / numpy.pi) + outer_multiplier = (-v / u) * numpy.sqrt(2 / numpy.pi) total_exponent = exponent_sum * outer_multiplier # 4. Bring back to standard space safely @@ -370,11 +369,11 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz tcoordRange = numpy.arange(0,tonumber(timeSpan,min),tonumber(dt,min)) stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) + u = meteorology.u10 # u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) inversion = tonumber(meteorology.inversion, m) - TX = self._getTXterm(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) + TX = self._getTXterm(stability=stability, u=tonumber(u, m/min), 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) @@ -388,7 +387,7 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz def getDosageFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) + u = meteorology.u10 # u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) inversion = tonumber(meteorology.inversion, m) @@ -397,7 +396,7 @@ def getDosageFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, zcoordRange = numpy.arange(tonumber(minz, m), tonumber(maxz, m), tonumber(dz, m)) tcoordRange = numpy.arange(0, tonumber(timeSpan, min), tonumber(dt, min)) - TX = self._getTXDosage(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) + TX = self._getTXDosage(stability=stability, u=tonumber(u, m/min), 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) From a82c1a0eb4823a451d15a019ca09a7e9cf4c4732 Mon Sep 17 00:00:00 2001 From: jmilo Date: Thu, 18 Jun 2026 17:12:19 +0300 Subject: [PATCH 5/9] fixing bug with DF --- hera/simulations/gaussian/gasCloud.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py index ba40e67fd..53f80c1e2 100644 --- a/hera/simulations/gaussian/gasCloud.py +++ b/hera/simulations/gaussian/gasCloud.py @@ -310,8 +310,8 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): :return: The Depletion Factor """ # Deposition velocity normalized to standard m/s units - v = tonumber(0.003 * m / s, m / min) #Deposition velocity. By default we take this value to be 0.003 [m/s] - u = tonumber(u, m / min) + v = tonumber(0.003 * m/s, m/min) #Deposition velocity. By default we take this value to be 0.003 [m/s] + u = tonumber(u, m/min) X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] From dcc02d848bcd4f7a62b230b036e34073e9f8ab9a Mon Sep 17 00:00:00 2001 From: jmilo Date: Sun, 21 Jun 2026 10:14:54 +0300 Subject: [PATCH 6/9] general fixes, some lost progress --- hera/simulations/gaussian/toolkit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hera/simulations/gaussian/toolkit.py b/hera/simulations/gaussian/toolkit.py index 50c41787d..c81911937 100644 --- a/hera/simulations/gaussian/toolkit.py +++ b/hera/simulations/gaussian/toolkit.py @@ -310,7 +310,7 @@ 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.ylabel(r"Dosage $\left[\frac{1}{m^3} \cdot min\right]$") plt.title(f"Dosage per distance. y={y}[m], z={z}[m], time={time}[min]") plt.grid() plt.show() @@ -335,7 +335,7 @@ 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"Dosage 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() plt.show() From 89addc4b91fc072c215388396500c830e224046f Mon Sep 17 00:00:00 2001 From: jmilo Date: Tue, 23 Jun 2026 11:34:04 +0300 Subject: [PATCH 7/9] general fixes, some lost progress adding limits to the plots in the presentations, adding the option to choose the wind profile type --- hera/simulations/gaussian/Meteorology.py | 13 +- hera/simulations/gaussian/gasCloud.py | 236 +++++++++++------------ hera/simulations/gaussian/toolkit.py | 112 ++++++++--- 3 files changed, 216 insertions(+), 145 deletions(-) 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 b17f46df6..e10580221 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, sigmaType): """ Parameters @@ -32,11 +32,18 @@ def __init__(self, sourceQ, sourceHeight, initialCloudSize, sigmaType): self.initialCloudSize = initialCloudSize self.sigmaType = sigmaType self.sourceQ = sourceQ + self.meteorology = meteorology + if wind_profile_type == 'HotSpot': + self.u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) + elif wind_profile_type == 'default': + self.u = self.meteorology.getWindVelocity(height=self.sourceHeight) + else: + raise "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,sigmaType): """ Return the type of the release based on the units of Q Parameters @@ -69,16 +76,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,sigmaType=sigmaType) - def _getTXterm(self, stability, u, xcoordRange, tcoordRange): + def _getTXterm(self, u, 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 +94,23 @@ def _getTXterm(self, stability, u, xcoordRange, tcoordRange): ------- The X component of the Gaussian concentration formula. """ + + u = tonumber(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 +119,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,14 +127,13 @@ def _getXYterm(self, stability, xcoordRange, ycoordRange): - def _getXZterm(self, stability, inversion, xcoordRange, zcoordRange, numOfReflections): + def _getXZterm(self, inversion, 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. @@ -138,7 +145,7 @@ def _getXZterm(self, stability, inversion, xcoordRange, zcoordRange, numOfReflec """ 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 @@ -167,7 +174,7 @@ def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dx 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)) + 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') @@ -180,13 +187,12 @@ def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dx - def _getTXDosage(self, stability, u, xcoordRange, tcoordRange): + def _getTXDosage(self, u, 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 +201,10 @@ def _getTXDosage(self, stability, u, xcoordRange, tcoordRange): ------- The X component of the Gaussian dosage formula. """ + + u = tonumber(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,30 +289,9 @@ def _getXZterm_ones(self, xcoordRange, zcoordRange): return XR_vertical - # def _getDF(self, stability, u, 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 - # """ - # - # v = tonumber(0.003*m/s, m/min) #Deposition velocity. By default we take this value to be 0.003 [m/s] - # v = tonumber(0.003 * m / s, m/s) - # u_ms = tonumber(u, m / s) - # X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') - # sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] - # H = tonumber(self.sourceHeight, m) - # dx = xcoordRange[1] - xcoordRange[0] - # DF = (numpy.e**(numpy.cumsum(1/(sigmaZ*numpy.e**(0.5*(H/sigmaZ)**2))*dx, axis=0)))**((-v/u_ms)*numpy.sqrt(2/numpy.pi)) - # - # return DF - - def _getDF(self, stability, u, xcoordRange, zcoordRange): - """ - :param stability: The stability state + + def _getDF(self, u, xcoordRange, zcoordRange): + """ :param u: Wind speed :param xcoordRange: x-coordinates of the grid :param zcoordRange: z-coordinates of the grid @@ -314,9 +301,8 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): 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(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 @@ -339,10 +325,13 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): - 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 = meteorology.u10.m_as(ureg.m/ureg.min) - # u = meteorology.getWindVelocity(height=self.sourceHeight) + def getDF_noQ_xarray(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*ureg.m, dz=1*ureg.m, + dt=1*ureg.min, wind_profile_type='HotSpot'): + + if wind_profile_type == 'HotSpot': + u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) + else: + u = self.meteorology.getWindVelocity(height=self.sourceHeight) 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)) @@ -352,7 +341,7 @@ def getDF_noQ_xarray(self, meteorology, minx, miny, minz, maxx, maxy, maxz, time 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(u=u, xcoordRange=xcoordRange, zcoordRange=zcoordRange) return TX*XY*(XZ*DF) @@ -362,22 +351,23 @@ 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, minx, miny, minz, maxx, maxy, maxz, timeSpan, + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, wind_profile_type='HotSpot'): 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)) - stability = meteorology.stability - u = meteorology.u10.m_as(ureg.m/ureg.min) - # u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) - inversion = meteorology.inversion.m_as(ureg.m) + if wind_profile_type == 'HotSpot': + u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) + else: + u = self.meteorology.getWindVelocity(height=self.sourceHeight) + inversion = self.meteorology.inversion.m_as(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, + TX = self._getTXterm(u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) + XY = self._getXYterm(xcoordRange=xcoordRange, ycoordRange=ycoordRange) + XZ = self._getXZterm(inversion=inversion, xcoordRange=xcoordRange, zcoordRange=zcoordRange, numOfReflections=numOfReflections) ret = TX*XY*XZ @@ -386,41 +376,43 @@ 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 = meteorology.u10.m_as(ureg.m/ureg.min) - # u = tonumber(meteorology.getWindVelocity(height=self.sourceHeight), m / min) - inversion = meteorology.inversion.m_as(ureg.m) + def getDosageFromMinMaxRange_inst_noQ(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + + if wind_profile_type == 'HotSpot': + u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) + else: + u = self.meteorology.getWindVelocity(height=self.sourceHeight) + inversion = self.meteorology.inversion.m_as(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)) - 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, + TX = self._getTXDosage(u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) + XY = self._getXYterm(xcoordRange=xcoordRange, ycoordRange=ycoordRange) + XZ = self._getXZterm(inversion=inversion, 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(u=u, 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, + def getDosageFromMinMaxRange_inst_NoERF_noQ(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, - dz=dz, dt=dt,numOfReflections=numOfReflections) + dz=dz, dt=dt,numOfReflections=numOfReflections,wind_profile_type=wind_profile_type) 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, + D_F = self.getDF_noQ_xarray(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt) ret = D_without_Q*D_F @@ -429,9 +421,9 @@ def getDosageFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, return ret - def getConcentrationFromMinMaxRange_inst(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, + def getConcentrationFromMinMaxRange_inst(self, 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, + C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, dz=dz, dt=dt,numOfReflections=numOfReflections) @@ -440,11 +432,12 @@ def getConcentrationFromMinMaxRange_inst(self, meteorology, minx, miny, minz, ma 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, minx, miny, minz, maxx, maxy, maxz, timeSpan, + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, + maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, + dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, + DF=DF, wind_profile_type=wind_profile_type) ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q ret.attrs['Q'] = 1*ureg.mg*ureg.min/ureg.m**3 @@ -452,12 +445,12 @@ 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, minx, miny, minz, maxx, maxy, maxz, timeSpan, + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): - D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, + D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy,maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) + dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q ret.attrs['Q'] = 1*ureg.mg*ureg.min/ureg.m**3 @@ -476,12 +469,13 @@ 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, 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, wind_profile_type='HotSpot'): - D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, + D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, - dt=dt, numOfReflections=numOfReflections, DF=DF) + dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) out_units = unumToPint(outputUnits) ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q @@ -491,12 +485,12 @@ 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, 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, wind_profile_type='HotSpot'): - D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, + D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, - dt=dt, numOfReflections=numOfReflections, DF=DF) + dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) out_units = unumToPint(outputUnits) currentUnites = unumToPint(D_without_Q.attrs['Q']) factor = currentUnites.m_as(out_units) @@ -506,12 +500,13 @@ 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=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, + outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False, wind_profile_type='HotSpot'): - D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, + D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) + dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) out_units = unumToPint(outputUnits) ret = tonumber(self.sourceQ, ureg.mg) * D_without_Q @@ -521,12 +516,13 @@ 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=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, + outputUnits=ureg.s/ureg.m**3, DF=False, wind_profile_type='HotSpot'): - D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, + D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) + dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) out_units = unumToPint(outputUnits) currentUnites = unumToPint(D_without_Q.attrs['Q']) @@ -538,7 +534,7 @@ 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, + def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, 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): """ @@ -558,7 +554,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, + D_F = self.getDF_noQ_xarray(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt) D_without_Q = D_without_Q * D_F @@ -570,7 +566,7 @@ 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, + def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, 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): """ @@ -582,7 +578,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, + D_F = self.getDF_noQ_xarray(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt) D_without_Q = D_without_Q*D_F @@ -609,8 +605,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: @@ -623,8 +619,8 @@ def get_TIAC_from_dist(self, data, y, z, dist_list): 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=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): """ Returns ------- @@ -632,15 +628,15 @@ 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, + C_without_Q = self.getDosageFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, - numOfReflections=numOfReflections, DF=DF) + numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) 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=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): """ Returns ------- @@ -648,28 +644,28 @@ 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, + C_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) + dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) 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=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): - C_without_Q = self.getConcentrationFromMinMaxRange_cont(meteorology=meteorology, minx=minx, miny=miny, minz=minz, + C_without_Q = self.getConcentrationFromMinMaxRange_cont(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) + dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) 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=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): - C_without_Q = self.getConcentrationFromMinMaxRange_cont_NoERF(meteorology=meteorology, minx=minx, miny=miny, minz=minz, + C_without_Q = self.getConcentrationFromMinMaxRange_cont_NoERF(minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) + dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) 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 9715b9d7f..7d4a46d57 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 @@ -24,6 +26,8 @@ def __init__(self, projectName: str, filesDirectory: str = None): self._sigmaDict = dict(briggsRural=BriggsRural) self._presentation = presentationLayer() + self.spaceTime = dict() + def getSigmaType(self,sigmaName): """ @@ -56,22 +60,69 @@ 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, stabilityType, windSpeed, maxx, dt, dxdy_multiplier, + dz=0.5, + minimal_maxy=0*ureg.m, + initialCloudSize=[0.1*ureg.m,0.1*ureg.m,0.1*ureg.m]): + + # 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=stabilityType, + 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 + + # when we only care about the concentration along the x-axis: Set the parameter minimal_maxy=0*m. + + 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, sigmaTypeName="briggsRural"): """ Parameters @@ -86,6 +137,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 +146,9 @@ 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,sigmaType=sigmaType) return gascloud @@ -227,7 +281,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 +310,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 +341,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 +367,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}{m^3} \cdot min\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 +397,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}{m^3} \cdot min\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() From 0d83aff84d4403bf9e02621b0e55de4f7b01961d Mon Sep 17 00:00:00 2001 From: jmilo Date: Tue, 23 Jun 2026 15:54:25 +0300 Subject: [PATCH 8/9] general fixes, some lost progress adding limits to the plots in the presentations, adding the option to choose the wind profile type --- hera/simulations/gaussian/gasCloud.py | 211 ++++++++++---------------- hera/simulations/gaussian/toolkit.py | 24 +-- 2 files changed, 93 insertions(+), 142 deletions(-) diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py index e10580221..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, meteorology, wind_profile_type, sigmaType): + def __init__(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind_profile_type, spaceTime, sigmaType): """ Parameters @@ -26,6 +26,8 @@ def __init__(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind_pr 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 @@ -33,17 +35,18 @@ def __init__(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind_pr self.sigmaType = sigmaType self.sourceQ = sourceQ self.meteorology = meteorology + self.spaceTime = spaceTime if wind_profile_type == 'HotSpot': - self.u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) + self.u = self.meteorology.getWindVelocity_hotSpot(height=sourceHeight) elif wind_profile_type == 'default': - self.u = self.meteorology.getWindVelocity(height=self.sourceHeight) + self.u = self.meteorology.getWindVelocity(height=sourceHeight) else: - raise "wind_profile_type must be either 'default' or 'HotSpot'" + raise ValueError("wind_profile_type must be either 'default' or 'HotSpot'") @staticmethod - def createGasCloud(sourceQ,sourceHeight,initialCloudSize,meteorology,wind_profile_type,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 @@ -77,10 +80,10 @@ def createGasCloud(sourceQ,sourceHeight,initialCloudSize,meteorology,wind_profil returnCls = instantaneousReleaseGasCloud if instantaneous else continuousReleaseGasCloud return returnCls(sourceQ=sourceQ,sourceHeight=sourceHeight,initialCloudSize=initialCloudSize, - meteorology=meteorology,wind_profile_type=wind_profile_type,sigmaType=sigmaType) + meteorology=meteorology,wind_profile_type=wind_profile_type, spaceTime=spaceTime, sigmaType=sigmaType) - def _getTXterm(self, u, xcoordRange, tcoordRange): + def _getTXterm(self, xcoordRange, tcoordRange): """ Parameters ---------- @@ -95,7 +98,7 @@ def _getTXterm(self, u, xcoordRange, tcoordRange): The X component of the Gaussian concentration formula. """ - u = tonumber(u, ureg.m / ureg.min) + u = tonumber(self.u, ureg.m / ureg.min) T, X = numpy.meshgrid(tcoordRange, xcoordRange, indexing='ij') 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)) @@ -127,14 +130,13 @@ def _getXYterm(self, xcoordRange, ycoordRange): - def _getXZterm(self, 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. - 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. @@ -143,6 +145,8 @@ def _getXZterm(self, inversion, xcoordRange, zcoordRange, numOfReflections): ------- 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=self.meteorology.stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] @@ -187,7 +191,7 @@ def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dx - def _getTXDosage(self, u, xcoordRange, tcoordRange): + def _getTXDosage(self, xcoordRange, tcoordRange): """ Parameters ---------- @@ -202,7 +206,7 @@ def _getTXDosage(self, u, xcoordRange, tcoordRange): The X component of the Gaussian dosage formula. """ - u = tonumber(u, ureg.m / ureg.min) + u = tonumber(self.u, ureg.m / ureg.min) T, X = numpy.meshgrid(tcoordRange, xcoordRange, indexing='ij') 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))) @@ -290,7 +294,7 @@ def _getXZterm_ones(self, xcoordRange, zcoordRange): - def _getDF(self, u, xcoordRange, zcoordRange): + def _getDF(self, xcoordRange, zcoordRange): """ :param u: Wind speed :param xcoordRange: x-coordinates of the grid @@ -299,7 +303,7 @@ def _getDF(self, u, xcoordRange, zcoordRange): """ # 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(u, ureg.m/ureg.min) + u = tonumber(self.u, ureg.m/ureg.min) X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') sigmaZ = self.sigmaType.getSigma(x=X, stability=self.meteorology.stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] @@ -325,23 +329,18 @@ def _getDF(self, u, xcoordRange, zcoordRange): - def getDF_noQ_xarray(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*ureg.m, dz=1*ureg.m, - dt=1*ureg.min, wind_profile_type='HotSpot'): + def getDF_noQ_xarray(self): - if wind_profile_type == 'HotSpot': - u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) - else: - u = self.meteorology.getWindVelocity(height=self.sourceHeight) - - 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(u=u, xcoordRange=xcoordRange, zcoordRange=zcoordRange) + DF = self._getDF(xcoordRange=xcoordRange, zcoordRange=zcoordRange) return TX*XY*(XZ*DF) @@ -351,24 +350,17 @@ def getDF_noQ_xarray(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10 class instantaneousReleaseGasCloud(abstractGasCloud): - def getConcentrationFromMinMaxRange_inst_noQ(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, wind_profile_type='HotSpot'): - - 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 getConcentrationFromMinMaxRange_inst_noQ(self, numOfReflections=3): - if wind_profile_type == 'HotSpot': - u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) - else: - u = self.meteorology.getWindVelocity(height=self.sourceHeight) - inversion = self.meteorology.inversion.m_as(ureg.m) + 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(u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) + TX = self._getTXterm(xcoordRange=xcoordRange, tcoordRange=tcoordRange) XY = self._getXYterm(xcoordRange=xcoordRange, ycoordRange=ycoordRange) - XZ = self._getXZterm(inversion=inversion, xcoordRange=xcoordRange, zcoordRange=zcoordRange, - numOfReflections=numOfReflections) + XZ = self._getXZterm(xcoordRange=xcoordRange, zcoordRange=zcoordRange, numOfReflections=numOfReflections) ret = TX*XY*XZ ret.attrs['Q'] = 1/ureg.m**3 @@ -376,44 +368,32 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, minx, miny, minz, maxx, maxy, return ret - def getDosageFromMinMaxRange_inst_noQ(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): - - if wind_profile_type == 'HotSpot': - u = self.meteorology.getWindVelocity_hotSpot(height=self.sourceHeight) - else: - u = self.meteorology.getWindVelocity(height=self.sourceHeight) - inversion = self.meteorology.inversion.m_as(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(u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) + TX = self._getTXDosage(xcoordRange=xcoordRange, tcoordRange=tcoordRange) XY = self._getXYterm(xcoordRange=xcoordRange, ycoordRange=ycoordRange) - XZ = self._getXZterm(inversion=inversion, xcoordRange=xcoordRange, zcoordRange=zcoordRange, - numOfReflections=numOfReflections) + XZ = self._getXZterm(xcoordRange=xcoordRange, zcoordRange=zcoordRange, numOfReflections=numOfReflections) D_F = 1 if DF: - D_F = self._getDF(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, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): - C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, - dz=dz, dt=dt,numOfReflections=numOfReflections,wind_profile_type=wind_profile_type) + 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(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 @@ -421,23 +401,16 @@ def getDosageFromMinMaxRange_inst_NoERF_noQ(self, minx, miny, minz, maxx, maxy, return ret - def getConcentrationFromMinMaxRange_inst(self, 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(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, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): - D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, - dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, - DF=DF, wind_profile_type=wind_profile_type) + 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 @@ -445,12 +418,9 @@ def getDosageFromMinMaxRange_inst(self, minx, miny, minz, maxx, maxy, maxz, time - def getDosageFromMinMaxRange_inst_NoERF(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + def getDosageFromMinMaxRange_inst_NoERF(self, numOfReflections=3, DF=False): - D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy,maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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 @@ -469,13 +439,9 @@ def concentrationConversion_mass_to_Bq(self, C, outputUnits, specificActivity): return C_Bq - def getTIACFromMinMaxRange_inst(self,specifitActivity, 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, wind_profile_type='HotSpot'): + def getTIACFromMinMaxRange_inst(self,specifitActivity, numOfReflections=3, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False): - D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, - maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, - dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(numOfReflections=numOfReflections, DF=DF) out_units = unumToPint(outputUnits) ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q @@ -485,12 +451,9 @@ def getTIACFromMinMaxRange_inst(self,specifitActivity, minx, miny, minz, maxx, m return ret - def getTIACFromMinMaxRange_inst_noQ(self, 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, wind_profile_type='HotSpot'): + def getTIACFromMinMaxRange_inst_noQ(self, numOfReflections=3, outputUnits=ureg.s/ureg.m**3, DF=False): - D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(minx=minx, miny=miny, minz=minz, maxx=maxx, - maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, - dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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) @@ -500,13 +463,10 @@ def getTIACFromMinMaxRange_inst_noQ(self, minx, miny, minz, maxx, maxy, maxz, ti - def getTIACFromMinMaxRange_inst_NoERF(self, specifitActivity, 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, wind_profile_type='HotSpot'): + 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(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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 @@ -516,13 +476,10 @@ def getTIACFromMinMaxRange_inst_NoERF(self, specifitActivity, minx, miny, minz, return ret - def getTIACFromMinMaxRange_inst_NoERF_noQ(self, 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, wind_profile_type='HotSpot'): + 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(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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']) @@ -534,8 +491,8 @@ def getTIACFromMinMaxRange_inst_NoERF_noQ(self, minx, miny, minz, maxx, maxy, ma - def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, 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]. @@ -554,8 +511,7 @@ def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, minx, miny, m D_F = 1 if DF: - D_F = self.getDF_noQ_xarray(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 @@ -566,8 +522,8 @@ def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, minx, miny, m return ret - def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, 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]. @@ -578,8 +534,7 @@ def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, minx, miny, minz, maxx, D_F = 1 if DF: - D_F = self.getDF_noQ_xarray(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] @@ -619,8 +574,8 @@ def get_TIAC_from_dist(self, data, y, z, dist_list): class continuousReleaseGasCloud(abstractGasCloud): - def getConcentrationFromMinMaxRange_cont(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + def getConcentrationFromMinMaxRange_cont(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt, + numOfReflections=3, DF=False): """ Returns ------- @@ -628,15 +583,13 @@ def getConcentrationFromMinMaxRange_cont(self, minx, miny, minz, maxx, maxy, max 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(minx=minx, miny=miny, minz=minz, maxx=maxx, - maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, - numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + def getConcentrationFromMinMaxRange_cont_NoERF(self, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy, dz, dt, + numOfReflections=3, DF=False): """ Returns ------- @@ -644,28 +597,22 @@ def getConcentrationFromMinMaxRange_cont_NoERF(self, minx, miny, minz, maxx, max 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(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + 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(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False, wind_profile_type='HotSpot'): + 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(minx=minx, miny=miny, minz=minz, - maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, - dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF, wind_profile_type=wind_profile_type) + 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 7d4a46d57..d98fba53e 100644 --- a/hera/simulations/gaussian/toolkit.py +++ b/hera/simulations/gaussian/toolkit.py @@ -26,7 +26,6 @@ def __init__(self, projectName: str, filesDirectory: str = None): self._sigmaDict = dict(briggsRural=BriggsRural) self._presentation = presentationLayer() - self.spaceTime = dict() def getSigmaType(self,sigmaName): @@ -76,10 +75,15 @@ def getMeteorologyFromURefHeight(self, u, refHeight, inversion, verticalProfileT - def getSpaceTime(self, stabilityType, windSpeed, maxx, dt, dxdy_multiplier, - dz=0.5, - minimal_maxy=0*ureg.m, - initialCloudSize=[0.1*ureg.m,0.1*ureg.m,0.1*ureg.m]): + 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). @@ -96,7 +100,7 @@ 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=stabilityType, + sigmaX = self.getSigmaType(sigmaName='briggsRural').getSigma(x=x, stability=meteorology.stability, sigma0=initialCloudSize, units=False)['sigmaX'][0] return x - 3 * sigmaX - maxx @@ -108,8 +112,6 @@ def get_function(x, maxx): timeSpan = x_timeSpan / (windSpeed * 60 * (ureg.s / ureg.min)) timeSpan = numpy.ceil(timeSpan.m_as(ureg.min)) * ureg.min - # when we only care about the concentration along the x-axis: Set the parameter minimal_maxy=0*m. - spaceTime = { 'minx': 0 * ureg.m, 'maxx': maxx, 'miny': -maxy, 'maxy': maxy + 1 * ureg.m, @@ -122,7 +124,8 @@ def get_function(x, maxx): - def getGasCloud(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind_profile_type, sigmaTypeName="briggsRural"): + def getGasCloud(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind_profile_type, + spaceTime, sigmaTypeName="briggsRural"): """ Parameters @@ -148,7 +151,8 @@ def getGasCloud(self, sourceQ, sourceHeight, initialCloudSize, meteorology, wind sigmaType = self.getSigmaType(sigmaTypeName) gascloud = abstractGasCloud.createGasCloud(sourceQ=sourceQ,sourceHeight=sourceHeight, initialCloudSize=initialCloudSize,meteorology=meteorology, - wind_profile_type=wind_profile_type,sigmaType=sigmaType) + wind_profile_type=wind_profile_type, spaceTime=spaceTime, + sigmaType=sigmaType) return gascloud From d9c0a371a352b2654f0a37c9f2a8371801264eb0 Mon Sep 17 00:00:00 2001 From: jmilo Date: Tue, 14 Jul 2026 11:27:18 +0300 Subject: [PATCH 9/9] reverting back to transposing the x and y axes. --- hera/simulations/LSM/template.py | 2 ++ 1 file changed, 2 insertions(+) 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: