From 2b73ac6e5a6a30d990bb7a84560a304fd67b2219 Mon Sep 17 00:00:00 2001 From: John Kendrick Date: Sun, 14 Nov 2021 11:24:31 +0000 Subject: [PATCH] Resolving a problem that seems to be associated with the introduction of AutoMie_ab into the MieS1S2 routine. Looking at the code for AutoMie_ab there is a hard crossover and this seems to have introduced the problem. There is also a check for a crossover in both MieQ and AutoMieQ which could cause a conflict. I have introduced a module wide variable called crossover and given it a value of 0.01. I have removed the crossover parameter from the AutoMieQ function call. I now use the module wide parameter to test for crossover in AutoMieQ and AutoMie_ab I have removed the crossover check in MieQ and left the full Mie expansion. The user is now able to change the crossover value him or herself. import PyMieScatt.Mie as Mie Mie.crossover = 0.001 This allows the use some control over the crossover if they are unhappy with the default setting. --- PyMieScatt/Mie.py | 91 +++++++++++++++++++++++------------------------ 1 file changed, 45 insertions(+), 46 deletions(-) diff --git a/PyMieScatt/Mie.py b/PyMieScatt/Mie.py index 958245a..7ddb99b 100644 --- a/PyMieScatt/Mie.py +++ b/PyMieScatt/Mie.py @@ -5,6 +5,10 @@ from scipy.integrate import trapz import warnings +# Parameter to determine wavelength crossover between methods for treating +# long and short wavelengths +crossover = 0.01 + def coerceDType(d): if type(d) is not np.ndarray: return np.array(d) @@ -17,51 +21,46 @@ def MieQ(m, wavelength, diameter, nMedium=1.0, asDict=False, asCrossSection=Fals m /= nMedium wavelength /= nMedium x = np.pi*diameter/wavelength - if x==0: - return 0, 0, 0, 1.5, 0, 0, 0 - elif x<=0.05: - return RayleighMieQ(m, wavelength, diameter, nMedium, asDict) - elif x>0.05: - nmax = np.round(2+x+4*(x**(1/3))) - n = np.arange(1,nmax+1) - n1 = 2*n+1 - n2 = n*(n+2)/(n+1) - n3 = n1/(n*(n+1)) - x2 = x**2 - - an,bn = Mie_ab(m,x) - - qext = (2/x2)*np.sum(n1*(an.real+bn.real)) - qsca = (2/x2)*np.sum(n1*(an.real**2+an.imag**2+bn.real**2+bn.imag**2)) - qabs = qext-qsca - - g1 = [an.real[1:int(nmax)], - an.imag[1:int(nmax)], - bn.real[1:int(nmax)], - bn.imag[1:int(nmax)]] - g1 = [np.append(x, 0.0) for x in g1] - g = (4/(qsca*x2))*np.sum((n2*(an.real*g1[0]+an.imag*g1[1]+bn.real*g1[2]+bn.imag*g1[3]))+(n3*(an.real*bn.real+an.imag*bn.imag))) - - qpr = qext-qsca*g - qback = (1/x2)*(np.abs(np.sum(n1*((-1)**n)*(an-bn)))**2) - qratio = qback/qsca - if asCrossSection: - css = np.pi*(diameter/2)**2 - cext = css*qext - csca = css*qsca - cabs = css*qabs - cpr = css*qpr - cback = css*qback - cratio = css*qratio - if asDict: - return dict(Cext=cext,Csca=csca,Cabs=cabs,g=g,Cpr=cpr,Cback=cback,Cratio=cratio) - else: - return cext, csca, cabs, g, cpr, cback, cratio + nmax = np.round(2+x+4*(x**(1/3))) + n = np.arange(1,nmax+1) + n1 = 2*n+1 + n2 = n*(n+2)/(n+1) + n3 = n1/(n*(n+1)) + x2 = x**2 + + an,bn = Mie_ab(m,x) + + qext = (2/x2)*np.sum(n1*(an.real+bn.real)) + qsca = (2/x2)*np.sum(n1*(an.real**2+an.imag**2+bn.real**2+bn.imag**2)) + qabs = qext-qsca + + g1 = [an.real[1:int(nmax)], + an.imag[1:int(nmax)], + bn.real[1:int(nmax)], + bn.imag[1:int(nmax)]] + g1 = [np.append(x, 0.0) for x in g1] + g = (4/(qsca*x2))*np.sum((n2*(an.real*g1[0]+an.imag*g1[1]+bn.real*g1[2]+bn.imag*g1[3]))+(n3*(an.real*bn.real+an.imag*bn.imag))) + + qpr = qext-qsca*g + qback = (1/x2)*(np.abs(np.sum(n1*((-1)**n)*(an-bn)))**2) + qratio = qback/qsca + if asCrossSection: + css = np.pi*(diameter/2)**2 + cext = css*qext + csca = css*qsca + cabs = css*qabs + cpr = css*qpr + cback = css*qback + cratio = css*qratio + if asDict: + return dict(Cext=cext,Csca=csca,Cabs=cabs,g=g,Cpr=cpr,Cback=cback,Cratio=cratio) else: - if asDict: - return dict(Qext=qext,Qsca=qsca,Qabs=qabs,g=g,Qpr=qpr,Qback=qback,Qratio=qratio) - else: - return qext, qsca, qabs, g, qpr, qback, qratio + return cext, csca, cabs, g, cpr, cback, cratio + else: + if asDict: + return dict(Qext=qext,Qsca=qsca,Qabs=qabs,g=g,Qpr=qpr,Qback=qback,Qratio=qratio) + else: + return qext, qsca, qabs, g, qpr, qback, qratio def Mie_ab(m,x): # http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_ab @@ -169,7 +168,7 @@ def RayleighMieQ(m, wavelength, diameter, nMedium=1.0, asDict=False, asCrossSect else: return qext, qsca, qabs, g, qpr, qback, qratio -def AutoMieQ(m, wavelength, diameter, nMedium=1.0, crossover=0.01, asDict=False, asCrossSection=False): +def AutoMieQ(m, wavelength, diameter, nMedium=1.0, asDict=False, asCrossSection=False): # http://pymiescatt.readthedocs.io/en/latest/forward.html#AutoMieQ nMedium = nMedium.real m_eff = m / nMedium @@ -247,7 +246,7 @@ def LowFrequencyMie_ab(m,x): return an,bn def AutoMie_ab(m,x): - if x<0.5: + if x