From 43365a75e48a6dcb43f880587e3c9212d0de7f6d Mon Sep 17 00:00:00 2001 From: Vishal Bagadia Date: Wed, 8 Nov 2023 00:33:22 -0500 Subject: [PATCH 1/3] Added Integration Method Flag for Mie_SD and Mie_LogNormal; Includes original trapezoid integration and added support for simpson integration --- PyMieScatt/Mie.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/PyMieScatt/Mie.py b/PyMieScatt/Mie.py index 958245a..5e1118d 100644 --- a/PyMieScatt/Mie.py +++ b/PyMieScatt/Mie.py @@ -3,6 +3,7 @@ import numpy as np from scipy.special import jv, yv from scipy.integrate import trapz +from scipy.integrate import simpson import warnings def coerceDType(d): @@ -252,7 +253,7 @@ def AutoMie_ab(m,x): else: return Mie_ab(m,x) -def Mie_SD(m, wavelength, dp, ndp, nMedium=1.0, SMPS=True, interpolate=False, asDict=False): +def Mie_SD(m, wavelength, dp, ndp, integrate_method:str, nMedium=1.0, SMPS=True, interpolate=False, asDict=False): # http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD nMedium = nMedium.real m /= nMedium @@ -283,7 +284,7 @@ def Mie_SD(m, wavelength, dp, ndp, nMedium=1.0, SMPS=True, interpolate=False, as Bratio = np.sum(Q_ratio*aSDn) bigG = np.sum(g*Q_sca*aSDn)/np.sum(Q_sca*aSDn) Bpr = Bext - bigG*Bsca - else: + elif integrate_method=='trapz': Bext = trapz(Q_ext*aSDn,dp) Bsca = trapz(Q_sca*aSDn,dp) Babs = Bext-Bsca @@ -291,6 +292,14 @@ def Mie_SD(m, wavelength, dp, ndp, nMedium=1.0, SMPS=True, interpolate=False, as Bratio = trapz(Q_ratio*aSDn,dp) bigG = trapz(g*Q_sca*aSDn,dp)/trapz(Q_sca*aSDn,dp) Bpr = Bext - bigG*Bsca + elif integrate_method=='simpson': + Bext = simpson(Q_ext*aSDn,dp) + Bsca = simpson(Q_sca*aSDn,dp) + Babs = Bext-Bsca + Bback = simpson(Q_back*aSDn,dp) + Bratio = simpson(Q_ratio*aSDn,dp) + bigG = simpson(g*Q_sca*aSDn,dp)/simpson(Q_sca*aSDn,dp) + Bpr = Bext - bigG*Bsca if asDict: return dict(Bext=Bext, Bsca=Bsca, Babs=Babs, G=bigG, Bpr=Bpr, Bback=Bback, Bratio=Bratio) @@ -489,7 +498,7 @@ def MieQ_withSizeParameterRange(m, nMedium=1.0, xRange=(1,10), nx=1000, logX=Fal qratio = np.array([q[6] for q in _qD]) return xValues, qext, qsca, qabs, g, qpr, qback, qratio -def Mie_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,nMedium=1.0, numberOfBins=10000,lower=1,upper=1000,gamma=[1],returnDistribution=False,decomposeMultimodal=False,asDict=False): +def Mie_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_method:str,nMedium=1.0, numberOfBins=10000,lower=1,upper=1000,gamma=[1],returnDistribution=False,decomposeMultimodal=False,asDict=False): # http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_Lognormal nMedium = nMedium.real m /= nMedium @@ -519,7 +528,7 @@ def Mie_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,nMedium=1.0, ndp = numberOfParticles*ithPart(1,dp,geoMean,geoStdDev) if ndp[-1]>np.max(ndp)/100 or ndp[0]>np.max(ndp)/100: warnings.warn("Warning: distribution may not be compact on the specified interval. Consider using a higher upper bound.") - Bext, Bsca, Babs, bigG, Bpr, Bback, Bratio = Mie_SD(m,wavelength,dp,ndp,SMPS=False) + Bext, Bsca, Babs, bigG, Bpr, Bback, Bratio = Mie_SD(m,wavelength,dp,ndp,integrate_method,SMPS=False) if returnDistribution: if decomposeMultimodal: if asDict==True: From 61463499d53ae3fb2e693ae9ac8ade784cb5555e Mon Sep 17 00:00:00 2001 From: Vishal Bagadia Date: Wed, 8 Nov 2023 00:39:20 -0500 Subject: [PATCH 2/3] Additional Size Distribution descriptive parameters added, incl. Surface Area and Volume descriptors for multimodal size distributions and lognormal size distributions. --- PyMieScatt/Mie.py | 161 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 161 insertions(+) diff --git a/PyMieScatt/Mie.py b/PyMieScatt/Mie.py index 5e1118d..57d38e1 100644 --- a/PyMieScatt/Mie.py +++ b/PyMieScatt/Mie.py @@ -545,3 +545,164 @@ def Mie_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_met return dict(Bext=Bext, Bsca=Bsca, Babs=Babs, bigG=bigG, Bpr=Bpr, Bback=Bback, Bratio=Bratio) else: return Bext, Bsca, Babs, bigG, Bpr, Bback, Bratio + +def SurfaceArea_SD(m, wavelength, dp, ndp,ndpi,integrate_method:str, nMedium=1.0, SMPS=True, interpolate=False, asDict=False): + '''Returns Surface Area of particle for given Size Distribution in units of um2/cm3. Returns total Surface Area (SArea), Surface Area Distribution (dSD/ddp), and Surface Area Distribution per mode (dSD/ddpi). Choice of integration method''' +# http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD + nMedium = nMedium.real + m /= nMedium + wavelength /= nMedium + dp = coerceDType(dp) + ndp = coerceDType(ndp) + _length = np.size(dp) + + # scaling of 1e-6 to cast in units of um2/cm3. + aSDn = 4*np.pi*((dp/2)**2)*ndp*(1e-6) + # ith_SD = lambda ndpi, dp : 4*np.pi*((dp/2)**2)*ndpi*(1E-6) + SDdp = aSDn + SDdpi = [] + for mode in ndpi: + aSDni = 4*np.pi*((dp/2)**2)*mode*(1e-6) + SDdpi.append(aSDni) +# _logdp = np.log10(dp) + if SMPS: + SArea = np.sum(aSDn) + elif integrate_method=='trapz': + SArea = trapz(aSDn,dp) + elif integrate_method=='simpson': + SArea = simpson(aSDn,dp) + + if asDict: + return dict(SArea=SArea,SDdp=SDdp,SDdpi=SDdpi) + else: + return SArea,SDdp,SDdpi + +def SurfaceArea_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_method:str,nMedium=1.0, numberOfBins=10000,lower=1,upper=1000,gamma=[1],returnDistribution=False,decomposeMultimodal=False,asDict=False): + '''Returns Surface Area and Surface Area Size Distribution Parameters, and input distribution particle size and number if specified. Returns surface area in units of um2/cm3.''' +# http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_Lognormal + nMedium = nMedium.real + m /= nMedium + wavelength /= nMedium + ithPart = lambda gammai, dp, dpgi, sigmagi: (gammai/(np.sqrt(2*np.pi)*np.log(sigmagi)*dp))*np.exp(-(np.log(dp)-np.log(dpgi))**2/(2*np.log(sigmagi)**2)) + dp = np.logspace(np.log10(lower),np.log10(upper),numberOfBins) + if all([type(x) in [list, tuple, np.ndarray] for x in [geoStdDev, geoMean]]): + # multimodal + if len(gamma)==1 and (len(geoStdDev)==len(geoMean)>1): + # gamma is distributed equally among modes + gamma = [1 for x in geoStdDev] + gamma = [float(x/np.sum(gamma)) for x in gamma] + ndpi = [numberOfParticles*ithPart(g,dp,dpg,sg) for g,dpg,sg in zip(gamma,geoMean,geoStdDev)] + ndp = np.sum(ndpi,axis=0) + elif len(gamma)==len(geoStdDev)==len(geoMean): + # gamma is fully specified for each mode + gamma = [float(x/np.sum(gamma)) for x in gamma] + ndpi = [numberOfParticles*ithPart(g,dp,dpg,sg) for g,dpg,sg in zip(gamma,geoMean,geoStdDev)] + ndp = np.sum(ndpi,axis=0) + else: + # user fucked up + warnings.warn("Not enough parameters to fully specify each mode.") + return None + else: + # unimodal + decomposeMultimodal = False + ndp = numberOfParticles*ithPart(1,dp,geoMean,geoStdDev) + if ndp[-1]>np.max(ndp)/100 or ndp[0]>np.max(ndp)/100: + warnings.warn("Warning: distribution may not be compact on the specified interval. Consider using a higher upper bound.") + + SArea, SDdp, SDdpi = SurfaceArea_SD(m,wavelength,dp,ndp,ndpi,integrate_method,SMPS=False) + if returnDistribution: + if decomposeMultimodal: + if asDict==True: + return dict(SArea,SDdp,SDdpi), dp, ndp, ndpi + else: + return SArea,SDdp,SDdpi,dp,ndp,ndpi + else: + if asDict==True: + return dict(SArea,SDdp), dp, ndp + else: + return SArea,SDdp,dp,ndp + else: + if asDict==True: + return dict(SArea) + else: + return SArea + +def Volume_SD(m, wavelength, dp, ndp, ndpi,integrate_method:str, nMedium=1.0, SMPS=True, interpolate=False, asDict=False): + '''Returns Volume of particle for given Size Distribution in units of um3/cm3. Returns total Volume (V), Volume Distribution (dV/ddp), and Volume Distribution per mode (dV/ddpi). Choice of integration method''' +# http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD + nMedium = nMedium.real + m /= nMedium + wavelength /= nMedium + dp = coerceDType(dp) + ndp = coerceDType(ndp) + _length = np.size(dp) + + # scaling of 1e-9 to cast in units of um3/cm3 + aVn = (4/3)*np.pi*((dp/2)**3)*ndp*(1e-9) + dVdp = aVn + dVdpi = [] + for mode in ndpi: + aVni = (4/3)*np.pi*((dp/2)**3)*mode*(1e-9) + dVdpi.append(aVni) +# _logdp = np.log10(dp) + if SMPS: + Volume = np.sum(aVn) + elif integrate_method=='trapz': + Volume = trapz(aVn,dp) + elif integrate_method=='simpson': + Volume = simpson(aVn,dp) + + if asDict: + return dict(Volume=Volume,dVdp=dVdp,dVdpi=dVdpi) + else: + return Volume,dVdp,dVdpi + +def Volume_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_method:str,nMedium=1.0, numberOfBins=10000,lower=1,upper=1000,gamma=[1],returnDistribution=False,decomposeMultimodal=False,asDict=False): + '''Returns Volume and Volume Size Distribution parameters, and input distribution particle size and number if specified. Returns Volume in units of um3/cm3.''' +# http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_Lognormal + nMedium = nMedium.real + m /= nMedium + wavelength /= nMedium + ithPart = lambda gammai, dp, dpgi, sigmagi: (gammai/(np.sqrt(2*np.pi)*np.log(sigmagi)*dp))*np.exp(-(np.log(dp)-np.log(dpgi))**2/(2*np.log(sigmagi)**2)) + dp = np.logspace(np.log10(lower),np.log10(upper),numberOfBins) + if all([type(x) in [list, tuple, np.ndarray] for x in [geoStdDev, geoMean]]): + # multimodal + if len(gamma)==1 and (len(geoStdDev)==len(geoMean)>1): + # gamma is distributed equally among modes + gamma = [1 for x in geoStdDev] + gamma = [float(x/np.sum(gamma)) for x in gamma] + ndpi = [numberOfParticles*ithPart(g,dp,dpg,sg) for g,dpg,sg in zip(gamma,geoMean,geoStdDev)] + ndp = np.sum(ndpi,axis=0) + elif len(gamma)==len(geoStdDev)==len(geoMean): + # gamma is fully specified for each mode + gamma = [float(x/np.sum(gamma)) for x in gamma] + ndpi = [numberOfParticles*ithPart(g,dp,dpg,sg) for g,dpg,sg in zip(gamma,geoMean,geoStdDev)] + ndp = np.sum(ndpi,axis=0) + else: + # user fucked up + warnings.warn("Not enough parameters to fully specify each mode.") + return None + else: + # unimodal + decomposeMultimodal = False + ndp = numberOfParticles*ithPart(1,dp,geoMean,geoStdDev) + if ndp[-1]>np.max(ndp)/100 or ndp[0]>np.max(ndp)/100: + warnings.warn("Warning: distribution may not be compact on the specified interval. Consider using a higher upper bound.") + + Volume,dVdp,dVdpi = Volume_SD(m,wavelength,dp,ndp,ndpi,integrate_method,SMPS=False) + if returnDistribution: + if decomposeMultimodal: + if asDict==True: + return dict(Volume=Volume,dVdp=dVdp,dVdpi=dVdpi), dp, ndp, ndpi + else: + return Volume, dVdp,dVdpi,dp,ndp,ndpi + else: + if asDict==True: + return dict(Volume=Volume,dVdp=dVdp), dp, ndp + else: + return Volume,dVdp,dp,ndp + else: + if asDict==True: + return dict(Volume) + else: + return Volume \ No newline at end of file From c36b270c5d6d1a0689a5ff4501f9394282f4abd4 Mon Sep 17 00:00:00 2001 From: Vishal Bagadia Date: Wed, 8 Nov 2023 10:25:37 -0500 Subject: [PATCH 3/3] Fixed SDist functions for SA and Vol --- PyMieScatt/Mie.py | 50 ++++++++++++++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/PyMieScatt/Mie.py b/PyMieScatt/Mie.py index 57d38e1..25e59bc 100644 --- a/PyMieScatt/Mie.py +++ b/PyMieScatt/Mie.py @@ -546,7 +546,7 @@ def Mie_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_met else: return Bext, Bsca, Babs, bigG, Bpr, Bback, Bratio -def SurfaceArea_SD(m, wavelength, dp, ndp,ndpi,integrate_method:str, nMedium=1.0, SMPS=True, interpolate=False, asDict=False): +def SurfaceArea_SD(m, wavelength, dp, ndp,integrate_method:str, nMedium=1.0, ndpi=False, SMPS=True, interpolate=False, asDict=False): '''Returns Surface Area of particle for given Size Distribution in units of um2/cm3. Returns total Surface Area (SArea), Surface Area Distribution (dSD/ddp), and Surface Area Distribution per mode (dSD/ddpi). Choice of integration method''' # http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD nMedium = nMedium.real @@ -560,10 +560,11 @@ def SurfaceArea_SD(m, wavelength, dp, ndp,ndpi,integrate_method:str, nMedium=1.0 aSDn = 4*np.pi*((dp/2)**2)*ndp*(1e-6) # ith_SD = lambda ndpi, dp : 4*np.pi*((dp/2)**2)*ndpi*(1E-6) SDdp = aSDn - SDdpi = [] - for mode in ndpi: - aSDni = 4*np.pi*((dp/2)**2)*mode*(1e-6) - SDdpi.append(aSDni) + if ndpi: + SDdpi = [] + for mode in ndpi: + aSDni = 4*np.pi*((dp/2)**2)*mode*(1e-6) + SDdpi.append(aSDni) # _logdp = np.log10(dp) if SMPS: SArea = np.sum(aSDn) @@ -572,10 +573,16 @@ def SurfaceArea_SD(m, wavelength, dp, ndp,ndpi,integrate_method:str, nMedium=1.0 elif integrate_method=='simpson': SArea = simpson(aSDn,dp) - if asDict: - return dict(SArea=SArea,SDdp=SDdp,SDdpi=SDdpi) + if ndpi: + if asDict: + return dict(SArea=SArea,SDdp=SDdp,SDdpi=SDdpi) + else: + return SArea,SDdp,SDdpi else: - return SArea,SDdp,SDdpi + if asDict: + return dict(SArea=SArea,SDdp=SDdp) + else: + return SArea,SDdp def SurfaceArea_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_method:str,nMedium=1.0, numberOfBins=10000,lower=1,upper=1000,gamma=[1],returnDistribution=False,decomposeMultimodal=False,asDict=False): '''Returns Surface Area and Surface Area Size Distribution Parameters, and input distribution particle size and number if specified. Returns surface area in units of um2/cm3.''' @@ -609,7 +616,7 @@ def SurfaceArea_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integ if ndp[-1]>np.max(ndp)/100 or ndp[0]>np.max(ndp)/100: warnings.warn("Warning: distribution may not be compact on the specified interval. Consider using a higher upper bound.") - SArea, SDdp, SDdpi = SurfaceArea_SD(m,wavelength,dp,ndp,ndpi,integrate_method,SMPS=False) + SArea, SDdp, SDdpi = SurfaceArea_SD(m,wavelength,dp,ndp,integrate_method,ndpi=ndpi,SMPS=False) if returnDistribution: if decomposeMultimodal: if asDict==True: @@ -627,7 +634,7 @@ def SurfaceArea_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integ else: return SArea -def Volume_SD(m, wavelength, dp, ndp, ndpi,integrate_method:str, nMedium=1.0, SMPS=True, interpolate=False, asDict=False): +def Volume_SD(m, wavelength, dp, ndp,integrate_method:str, nMedium=1.0, ndpi=False,SMPS=True, interpolate=False, asDict=False): '''Returns Volume of particle for given Size Distribution in units of um3/cm3. Returns total Volume (V), Volume Distribution (dV/ddp), and Volume Distribution per mode (dV/ddpi). Choice of integration method''' # http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD nMedium = nMedium.real @@ -640,10 +647,11 @@ def Volume_SD(m, wavelength, dp, ndp, ndpi,integrate_method:str, nMedium=1.0, S # scaling of 1e-9 to cast in units of um3/cm3 aVn = (4/3)*np.pi*((dp/2)**3)*ndp*(1e-9) dVdp = aVn - dVdpi = [] - for mode in ndpi: - aVni = (4/3)*np.pi*((dp/2)**3)*mode*(1e-9) - dVdpi.append(aVni) + if ndpi: + dVdpi = [] + for mode in ndpi: + aVni = (4/3)*np.pi*((dp/2)**3)*mode*(1e-9) + dVdpi.append(aVni) # _logdp = np.log10(dp) if SMPS: Volume = np.sum(aVn) @@ -652,10 +660,16 @@ def Volume_SD(m, wavelength, dp, ndp, ndpi,integrate_method:str, nMedium=1.0, S elif integrate_method=='simpson': Volume = simpson(aVn,dp) - if asDict: - return dict(Volume=Volume,dVdp=dVdp,dVdpi=dVdpi) + if ndpi: + if asDict: + return dict(Volume=Volume,dVdp=dVdp,dVdpi=dVdpi) + else: + return Volume,dVdp,dVdpi else: - return Volume,dVdp,dVdpi + if asDict: + return dict(Volume=Volume,dVdp=dVdp) + else: + return Volume,dVdp def Volume_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_method:str,nMedium=1.0, numberOfBins=10000,lower=1,upper=1000,gamma=[1],returnDistribution=False,decomposeMultimodal=False,asDict=False): '''Returns Volume and Volume Size Distribution parameters, and input distribution particle size and number if specified. Returns Volume in units of um3/cm3.''' @@ -689,7 +703,7 @@ def Volume_Lognormal(m,wavelength,geoStdDev,geoMean,numberOfParticles,integrate_ if ndp[-1]>np.max(ndp)/100 or ndp[0]>np.max(ndp)/100: warnings.warn("Warning: distribution may not be compact on the specified interval. Consider using a higher upper bound.") - Volume,dVdp,dVdpi = Volume_SD(m,wavelength,dp,ndp,ndpi,integrate_method,SMPS=False) + Volume,dVdp,dVdpi = Volume_SD(m,wavelength,dp,ndp,integrate_method,ndpi=ndpi,SMPS=False) if returnDistribution: if decomposeMultimodal: if asDict==True: