Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Comparing changes

Choose two branches to see what's changed or to start a new pull request. If you need to, you can also compare across forks.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also compare across forks.
base fork: ScattPy/scikits.scattpy
...
head fork: ScattPy/scikits.scattpy
Checking mergeability… Don't worry, you can still create the pull request.
  • 6 commits
  • 10 files changed
  • 0 commit comments
  • 2 contributors
View
202 scikits/scattpy/spheroidal.py
@@ -1,94 +1,9 @@
+from scipy import *
import scipy
import scipy.integrate
-import scipy.special as special
-import scipy.linalg as linalg
-import f_spheroid
-
-
-from scipy import *
-
-# ---- Calculation of norm for spheroidal angular functions ----
-optimize_norm={}
-def get_norm_cv(m, n, c, cv, type):
- sum = 0
- k = 0
- key = str(m)+' '+str(n)+' '+str(c)+' '+str(cv)
- if key in optimize_norm:
- return optimize_norm[key]
-
- d = f_spheroid.sdmn(m, n, c, cv, type)
- isEven = (n - m) % 2
- if(isEven == 1):
- k = 1
-
- while (k < d.shape[0]):
- if (type == 1):
- value = d[k / 2] * d[k / 2] / (2.0 * k + 2.0 * m + 1.0) * get_norm_factorial(k, m)
- else:
- value = d[k / 2] * d[k / 2] / (2.0 * k + 2.0 * m + 1.0) * get_norm_factorial(k, m)
- sum += value
- k += 2
-
- sum = sqrt(2.0 * sum)
- optimize_norm[key] = sum
- return sum
-
-def get_norm_factorial(k, m):
- fact = 1l
- for i in range(1, 2 * m + 1):
- fact *= k + i
- return fact
-
-
-def get_pro_norm_cv(m, n, c, cv):
- return get_norm_cv(m, n, c, cv, 1)
-
-
-def get_obl_norm_cv(m, n, c, cv):
- return get_norm_cv(m, n, c, cv, -1)
-
+import scipy.linalg
-def get_pro_norm(m, n, c):
- return get_pro_norm_cv(m, n, c, special.pro_cv(m, n, c))
-
-
-def get_obl_norm(m, n, c):
- return get_obl_norm_cv(m, n, c, special.obl_cv(m, n, c))
-
-# ---- Shortcut for different types of spheroids
-
-def rad1_cv(m, n, c, cv, type, x):
- if type == 1:
- return special.pro_rad1_cv(m, n, c, cv, x)
- elif type == -1:
- return special.obl_rad1_cv(m, n, c, cv, x)
-
-
-def rad2_cv(m, n, c, cv, type, x):
- if type == 1:
- return special.pro_rad2_cv(m, n, c, cv, x)
- elif type == -1:
- return special.obl_rad2_cv(m, n, c, cv, x)
-
-
-def ang1_cv(m, n, c, cv, type, x):
- if type == 1:
- value = special.pro_ang1_cv(m, n, c, cv, x)
- elif type == -1:
- value = special.obl_ang1_cv(m, n, c, cv, x)
- return value / get_norm_cv(m, n, c, cv, type)
-
-
-def get_cv(m, n, c, type):
- if type == 1:
- return special.pro_cv(m, n, c)
- elif type == -1:
- return special.obl_cv(m, n, c)
-
-#according to 1.41 Komarov, Slavyanov "Spheroidal funtions"
-def rad3_cv(m, n, c, cv, type, x):
- return [rad1_cv(m, n, c, cv, type, x)[0] + 1j * rad2_cv(m, n, c, cv, type, x)[0],
- rad1_cv(m, n, c, cv, type, x)[1] + 1j * rad2_cv(m, n, c, cv, type, x)[1]]
+from spheroidal_svm import *
#quad integration for complex numbers
def quad(func, a, b, **kwargs):
@@ -100,61 +15,40 @@ def imag_func(x):
imag_integral = scipy.integrate.quad(imag_func, a, b, **kwargs)
return real_integral[0] + 1j*imag_integral[0]
-# ---- Generation of A matrices
-
-#according to (86)
-def get_A11(particle, c1, nmax):
- return get_A(particle, c1, c1, 3, nmax).transpose()
-
-
-def get_A12(particle, c1, c2, nmax):
- return -get_A(particle, c2, c1, 1, nmax).transpose()
-
-
-def get_A10(particle, c1, nmax):
- return -get_A(particle, c1, c1, 1, nmax).transpose()
-
-
-def get_A21(particle, c1, nmax):
- return get_B(particle, c1, c1, 3, nmax).transpose()
-
-
-def get_A22(particle, c1, c2, nmax):
- return -get_C(particle, c2, c1, nmax).transpose()
-
-
-def get_A20(particle, c1, nmax):
- return -get_B(particle, c1, c1, 1, nmax).transpose()
-
#according to (81)
-def get_a_functions(m, n, c, cv, type, rank, particle):
+def get_a_function(m, n, c, cv, rank, particle):
+ type = particle.type
if rank == 1:
return lambda nu: rad1_cv(m, n, c, cv, type, particle.function(nu))[0] * ang1_cv(m, n, c, cv, type, nu)[0]
elif rank == 3:
return lambda nu: rad3_cv(m, n, c, cv, type, particle.function(nu))[0] * ang1_cv(m, n, c, cv, type, nu)[0]
#according to (81)
-def get_b_functions(m, n, c, cv, type, rank, particle):
+def get_b_function(m, n, c, cv, rank, particle):
+ type = particle.type
if rank == 1:
return lambda nu: (metric_nu(nu, particle) / metric_psi(nu, particle)
* rad1_cv(m, n, c, cv, type, particle.function(nu))[1] * ang1_cv(m, n, c, cv, type, nu)[0]
- particle.derivative(nu) * metric_psi(nu, particle) / metric_nu(nu, particle)
- * rad1_cv(m, n, c, cv, type, particle.function(nu))[0] * ang1_cv(m, n, c, cv, type, nu)[1])
+ * rad1_cv(m, n, c, cv, type, particle.function(nu))[0] * ang1_cv(m, n, c, cv, type, nu)[1]) \
+ / get_integral_metric(particle)(nu)
+
elif rank == 3:
return lambda nu: (metric_nu(nu, particle) / metric_psi(nu, particle)
* rad3_cv(m, n, c, cv, type, particle.function(nu))[1] * ang1_cv(m, n, c, cv, type, nu)[0]
- particle.derivative(nu) * metric_psi(nu, particle) / metric_nu(nu, particle)
- * rad3_cv(m, n, c, cv, type, particle.function(nu))[0] * ang1_cv(m, n, c, cv, type, nu)[1])
+ * rad3_cv(m, n, c, cv, type, particle.function(nu))[0] * ang1_cv(m, n, c, cv, type, nu)[1])\
+ / get_integral_metric(particle)(nu)
#according to (82)
-def get_c_functions(m, n, c, cv, type, rank, particle):
+def get_c_function(m, n, c, cv, rank, particle):
eps = particle.eps
#according to (30)
+ #is there a minus?
delta = lambda nu: metric_phi(nu, particle)
#according to (28)
- return lambda nu: get_b_functions(m, n, c, cv, type, 1, particle)(nu) / eps -\
- (1.0 / eps - 1.0) * IzIt(nu,particle) / delta(nu) * get_a_functions(m, n, c, cv, type, 1, particle)(nu) \
- * get_integral_metric(particle)(nu)
+ return lambda nu: get_b_function(m, n, c, cv, 1, particle)(nu) / eps -\
+ (1.0 / eps - 1.0) * IzIt(nu,particle) / delta(nu) * get_a_function(m, n, c, cv, 1, particle)(nu)
#-------------------------Metric coefficients ------------------------------------------------
@@ -193,62 +87,8 @@ def RIn(nu,particle):
/ get_integral_metric(particle)(nu)
def RIt(nu,particle):
return (particle.d / 2.0)**2 * (particle.derivative(nu)*particle.function(nu)+particle.type*nu) / get_integral_metric(particle)(nu)
-#----------------------------------------------------------------------------------------------
-
-#according to (83)
-def get_A(particle, c2, c1, rank, nmax):
- A = zeros((nmax, nmax),dtype=complex)
- type = particle.type
- m = 1
- for i in range(nmax):
- for k in range(nmax):
- l = i + m
- n = k + m
- cv_l = get_cv(m, l, c2, type)
- cv_n = get_cv(m, n, c1, type)
- func = lambda nu: get_a_functions(m, l, c2, cv_l, type, rank, particle)(nu) * ang1_cv(m, n, c1, cv_n, particle.type, nu)[0]
- A[i][k] = quad(func, -1, 1)
- return A
-
-#according to (84)
-def get_Z(get_z_functions, particle, c2, c1, rank, nmax):
- Z = zeros((nmax, nmax),dtype=complex)
- type = particle.type
- m = 1
- for i in range(nmax):
- for k in range(nmax):
- l = i + m
- n = k + m
- cv_l = get_cv(m, l, c2, type)
- cv_n = get_cv(m, n, c1, type)
- func = lambda nu: get_z_functions(m, l, c2, cv_l, type, rank, particle)(nu) *\
- ang1_cv(m, n, c1, cv_n, particle.type, nu)[0] * metric_phi(nu, particle)
- Z[i][k] = quad(func, -1, 1)
- return Z
-
-
-def get_C(particle, c2, c1, nmax):
- return get_Z(get_c_functions, particle, c2, c1, 1, nmax)
-
-
-def get_B(particle, c2, c1, rank, nmax):
- return get_Z(get_b_functions, particle, c2, c1, rank, nmax)
-
-#-------Solve equation and find solution of scattering by SVM
-
-#according to (85)
-def get_fullA(particle, c1, c2, nmax):
- A11 = get_A11(particle, c1, nmax)
- A12 = get_A12(particle, c1, c2, nmax)
- A21 = get_A21(particle, c1, nmax)
- A22 = get_A22(particle, c1, c2, nmax)
-
- return bmat([[A11, A12], [A21, A22]])
-
-
-def get_fullB(particle, c, nmax):
- return bmat([[get_A10(particle, c, nmax)], [get_A20(particle, c, nmax)]])
+#-------Solve equation and find solution of scattering
def get_Bin(inputWave, particle, c, nmax):
b = zeros((nmax, 1),dtype=complex)
@@ -258,8 +98,8 @@ def get_Bin(inputWave, particle, c, nmax):
return b
#Return b_sca and b_int. b_sca = result[0] and b_int = result[1]
-def getSolution(particle, inputWave, c1, c2, nmax):
- A = get_fullA(particle, c1, c2, nmax)
- B = get_fullB(particle, c1, nmax) * get_Bin(inputWave, particle, c1, nmax)
- x = linalg.solve(A, B)
+def getSolution(method,particle, inputWave, c1, c2, nmax):
+ A = method.get_fullA()
+ B = method.get_fullB() * get_Bin(inputWave, particle, c1, nmax)
+ x = scipy.linalg.solve(A, B)
return (x[0:nmax + 1], x[nmax + 1:])
View
83 scikits/scattpy/spheroidal_ebcm.py
@@ -0,0 +1,83 @@
+from numpy import zeros, asarray, mat
+import scipy
+
+from spheroidal_functions import *
+
+import spheroidal
+from spheroidal_svm import SpheroidalSVM
+
+
+class SpheroidalEBCM:
+ def __init__(self, particle, c2, c1, nmax):
+ self.particle = particle
+ self.c2 = c2
+ self.c1 = c1
+ self.nmax = nmax
+
+ #according to (88)
+ def get_BSi(self):
+ Bs = zeros((self.nmax, self.nmax), dtype=complex)
+ type = self.particle.type
+ m = 1
+ for i in range(self.nmax):
+ for k in range(self.nmax):
+ n = i + m
+ l = k + m
+ cv_l = get_cv(m, l, self.c2, type)
+ cv_n = get_cv(m, n, self.c1, type)
+ func = lambda nu: (spheroidal.get_a_function(m, l, self.c2, cv_l, 1, self.particle)(nu) *\
+ spheroidal.get_b_function(m, n, self.c1, cv_n, 3, self.particle)(nu) -\
+ spheroidal.get_c_function(m, l, self.c2, cv_l, 1, self.particle)(nu) *\
+ spheroidal.get_a_function(m, n, self.c1, cv_n, 3, self.particle)(nu)) * spheroidal.metric_phi(nu, self.particle) *\
+ spheroidal.get_integral_metric(self.particle)(nu)
+
+ Bs[i][k] = spheroidal.quad(func, -1, 1)
+ return 1j * mat(Bs)
+
+ #according to (88)
+ def get_BRi(self):
+ Br = zeros((self.nmax, self.nmax), dtype=complex)
+ type = self.particle.type
+ m = 1
+ for i in range(self.nmax):
+ for k in range(self.nmax):
+ n = i + m
+ l = k + m
+ cv_l = get_cv(m, l, self.c2, type)
+ cv_n = get_cv(m, n, self.c1, type)
+ func = lambda nu: (spheroidal.get_a_function(m, l, self.c2, cv_l, 1, self.particle)(nu) * \
+ spheroidal.get_b_function(m, n, self.c1, cv_n, 1, self.particle)(nu) - \
+ spheroidal.get_c_function(m, l, self.c2, cv_l, 1, self.particle)(nu) * \
+ spheroidal.get_a_function(m, n, self.c1, cv_n, 1, self.particle)(nu)) * spheroidal.metric_phi(nu, self.particle) *\
+ spheroidal.get_integral_metric(self.particle)(nu)
+
+ Br[i][k] = spheroidal.quad(func, -1, 1)
+ return -1j * mat(Br)
+
+ #according to (99)
+ def get_BSm(self):
+ svm = SpheroidalSVM(self.particle, self.c2, self.c1, self.nmax)
+ Bs = 1j * ( svm.get_B(3) * svm.get_A(self.c2, self.c1, 1).T -
+ svm.get_A(self.c1, self.c1, 3) * svm.get_C().T)
+ return Bs
+
+ #according to (99)
+ def get_BRm(self):
+ svm = SpheroidalSVM(self.particle, self.c2, self.c1, self.nmax)
+ Br = -1j * ( svm.get_B(1) * svm.get_A(self.c2, self.c1, 1).T -
+ svm.get_A(self.c1, self.c1, 1) * svm.get_C().T)
+ return Br
+
+
+ #Return b_sca and b_int. b_sca = result[0] and b_int = result[1]
+ def getSolution(self, inputWave):
+ b_in = spheroidal.get_Bin(inputWave, self.particle, self.c1, self.nmax)
+ b_int = scipy.linalg.solve(self.get_BSi(), -b_in)
+ b_sca = asarray(self.get_BRi() * b_int)
+ return (b_sca, b_int)
+
+ def getMatrixSolution(self, inputWave):
+ b_in = spheroidal.get_Bin(inputWave, self.particle, self.c1, self.nmax)
+ b_int = scipy.linalg.solve(self.get_BSm(), -b_in)
+ b_sca = asarray(self.get_BRm() * b_int)
+ return (b_sca, b_int)
View
88 scikits/scattpy/spheroidal_functions.py
@@ -0,0 +1,88 @@
+import scipy.special as special
+
+from numpy.lib.scimath import sqrt
+
+import f_spheroid
+
+# ---- Calculation of norm for spheroidal angular functions ----
+
+optimize_norm={}
+def get_norm_cv(m, n, c, cv, type):
+ sum = 0
+ k = 0
+ key = str(m)+' '+str(n)+' '+str(c)+' '+str(cv)
+ if key in optimize_norm:
+ return optimize_norm[key]
+
+ d = f_spheroid.sdmn(m, n, c, cv, type)
+ isEven = (n - m) % 2
+ if(isEven == 1):
+ k = 1
+
+ while (k < d.shape[0]):
+ if (type == 1):
+ value = d[k / 2] * d[k / 2] / (2.0 * k + 2.0 * m + 1.0) * get_norm_factorial(k, m)
+ else:
+ value = d[k / 2] * d[k / 2] / (2.0 * k + 2.0 * m + 1.0) * get_norm_factorial(k, m)
+ sum += value
+ k += 2
+
+ sum = sqrt(2.0 * sum)
+ optimize_norm[key] = sum
+ return sum
+
+def get_norm_factorial(k, m):
+ fact = 1l
+ for i in range(1, 2 * m + 1):
+ fact *= k + i
+ return fact
+
+def get_pro_norm_cv(m, n, c, cv):
+ return get_norm_cv(m, n, c, cv, 1)
+
+
+def get_obl_norm_cv(m, n, c, cv):
+ return get_norm_cv(m, n, c, cv, -1)
+
+
+def get_pro_norm(m, n, c):
+ return get_pro_norm_cv(m, n, c, special.pro_cv(m, n, c))
+
+
+def get_obl_norm(m, n, c):
+ return get_obl_norm_cv(m, n, c, special.obl_cv(m, n, c))
+
+# ---- Shortcut for different types of spheroids
+
+def rad1_cv(m, n, c, cv, type, x):
+ if type == 1:
+ return special.pro_rad1_cv(m, n, c, cv, x)
+ elif type == -1:
+ return special.obl_rad1_cv(m, n, c, cv, x)
+
+
+def rad2_cv(m, n, c, cv, type, x):
+ if type == 1:
+ return special.pro_rad2_cv(m, n, c, cv, x)
+ elif type == -1:
+ return special.obl_rad2_cv(m, n, c, cv, x)
+
+
+def ang1_cv(m, n, c, cv, type, x):
+ if type == 1:
+ value = special.pro_ang1_cv(m, n, c, cv, x)
+ elif type == -1:
+ value = special.obl_ang1_cv(m, n, c, cv, x)
+ return value / get_norm_cv(m, n, c, cv, type)
+
+
+def get_cv(m, n, c, type):
+ if type == 1:
+ return special.pro_cv(m, n, c)
+ elif type == -1:
+ return special.obl_cv(m, n, c)
+
+#according to 1.41 Komarov, Slavyanov "Spheroidal funtions"
+def rad3_cv(m, n, c, cv, type, x):
+ return [rad1_cv(m, n, c, cv, type, x)[0] + 1j * rad2_cv(m, n, c, cv, type, x)[0],
+ rad1_cv(m, n, c, cv, type, x)[1] + 1j * rad2_cv(m, n, c, cv, type, x)[1]]
View
24 scikits/scattpy/spheroidal_particles.py
@@ -2,9 +2,9 @@
import spheroidal
-class ProlateSpheroid:
- type = 1
- d = 4.
+class Spheroid:
+
+ d = 2.
c = 5.
def __init__(self,psi=2.,eps=1.,c=1,derivative=0):
@@ -19,22 +19,12 @@ def function(self, nu):
def derivative(self, nu):
return self.deriv
+class ProlateSpheroid(Spheroid):
+ type = 1
-class OblateSpheroid:
- type = -1
- d = 1
- eps = 1
- c = 1
-
- def __init__(self,psi=2,derivative=0):
- self.psi = psi
- self.deriv = derivative
-
- def function(self, nu):
- return self.psi
- def derivative(self, nu):
- return self.deriv
+class OblateSpheroid(Spheroid):
+ type = -1
#Incedent waves
class TEInputWave:
View
98 scikits/scattpy/spheroidal_svm.py
@@ -0,0 +1,98 @@
+from numpy import zeros, bmat, mat
+import scipy
+
+from spheroidal_functions import *
+
+import spheroidal
+
+
+class SpheroidalSVM:
+
+ def __init__(self,particle,c2,c1,nmax):
+ self.particle=particle
+ self.c2 = c2
+ self.c1 = c1
+ self.nmax = nmax
+
+ #according to (83)
+ def get_A(self, c2,c1,rank):
+ A = zeros((self.nmax, self.nmax),dtype=complex)
+ type = self.particle.type
+ m = 1
+ for i in range(self.nmax):
+ for k in range(self.nmax):
+ l = i + m
+ n = k + m
+ cv_l = get_cv(m, l, c2, type)
+ cv_n = get_cv(m, n, c1, type)
+ func = lambda nu: spheroidal.get_a_function(m, l, c2, cv_l, rank, self.particle)(nu) * ang1_cv(m, n, c1, cv_n, type, nu)[0]
+ A[i][k] = spheroidal.quad(func, -1, 1)
+ return mat(A)
+
+ #according to (84)
+ def __get_Z(self,z_function,c2,c1, rank):
+ Z = zeros((self.nmax, self.nmax),dtype=complex)
+ type = self.particle.type
+ m = 1
+ for i in range(self.nmax):
+ for k in range(self.nmax):
+ l = i + m
+ n = k + m
+ cv_l = get_cv(m, l, c2, type)
+ cv_n = get_cv(m, n, c1, type)
+ func = lambda nu: z_function(m, l, c2, cv_l, rank, self.particle)(nu) *\
+ ang1_cv(m, n, c1, cv_n, type, nu)[0] * spheroidal.metric_phi(nu, self.particle) * \
+ spheroidal.get_integral_metric(self.particle)(nu)
+ Z[i][k] = spheroidal.quad(func, -1, 1)
+ return mat(Z)
+
+
+ def get_C(self):
+ return self.__get_Z(spheroidal.get_c_function, self.c2, self.c1, 1)
+
+
+ def get_B(self, rank):
+ return self.__get_Z(spheroidal.get_b_function, self.c1, self.c1, rank)
+
+ # ---- Generation of A matrices
+
+ #according to (86)
+ def get_A11(self):
+ return self.get_A(self.c1, self.c1, 3).T
+
+
+ def get_A12(self):
+ return -self.get_A(self.c2, self.c1, 1).T
+
+
+ def get_A10(self):
+ return -self.get_A(self.c1, self.c1, 1).T
+
+
+ def get_A21(self):
+ return self.get_B(3).T
+
+
+ def get_A22(self):
+ return -self.get_C().T
+
+
+ def get_A20(self):
+ return -self.get_B(1).T
+
+ #according to (85)
+ def get_fullB(self):
+ return bmat([[self.get_A10()], [self.get_A20()]])
+
+ def get_fullA(self):
+ A11 = self.get_A11()
+ A12 = self.get_A12()
+ A21 = self.get_A21()
+ A22 = self.get_A22()
+ return bmat([[A11, A12], [A21, A22]])
+
+ def getSolution(self,inputWave):
+ A = self.get_fullA()
+ B = self.get_fullB() * spheroidal.get_Bin(inputWave, self.particle, self.c1, self.nmax)
+ x = -scipy.linalg.solve(A, B)
+ return (x[0:self.nmax + 1], x[self.nmax + 1:])
View
19 test/draw_spheroidal.py
@@ -6,9 +6,9 @@
import matplotlib.pyplot as plt
import numpy as np
-from scikits.scattpy.spheroidal import *
-from scikits.scattpy.spheroidal_particles import *
-from scikits.scattpy.properties import *
+from spheroidal import *
+from spheroidal_particles import *
+from properties import *
def plot_graphics():
y=[]
@@ -18,16 +18,18 @@ def plot_graphics():
alpha = pi / 4
k = 1
particle = ProlateSpheroid(psi=2,c=c1,derivative=0,eps=1)
- for i in range(6,32,2):
+ for i in range(6,52,2):
nmax = i
+ print nmax
start = time.time()
- b_sca = getSolution(particle, TMInputWave(alpha), c1, c2, nmax)[0]
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ b_sca = getSolution(svm,particle, TMInputWave(alpha), c1, c2, nmax)[0]
C_ext = -getCext(particle, alpha, k, b_sca, nmax)[0]
C_sca = getCsca(k, b_sca, nmax)[0]
delta=(C_ext-C_sca)/(C_ext+C_sca)
execution_time.append(time.time() - start)
y.append(delta)
- x = np.arange(6,12,2)
+ x = np.arange(6,52,2)
y= np.fabs(y)
plt.plot(x,np.log10(y))
plt.savefig('error')
@@ -35,4 +37,7 @@ def plot_graphics():
plt.plot(x,np.log10(execution_time))
plt.savefig('time')
-plot_graphics()
+if __name__ == '__main__':
+ plot_graphics()
+
+
View
57 test/ebcm_test.py
@@ -0,0 +1,57 @@
+import unittest
+
+from numpy.linalg import det
+import numpy
+
+from spheroidal import *
+from spheroidal_ebcm import SpheroidalEBCM
+from spheroidal_particles import *
+from properties import *
+
+#main integration tests for non-absorbing particle
+class testEBCM(unittest.TestCase):
+
+ def test_IntegrationCase1(self):
+ nmax = 6
+ c1 = 1
+ c2 = 2
+ alpha = pi / 3
+ k = 1
+ particle = ProlateSpheroid(psi=10,c=c1,derivative=0,eps=10)
+ ebcm = SpheroidalEBCM(particle,c2,c1,nmax)
+ b_sca = ebcm.getSolution(TMInputWave(alpha))[0]
+ C_ext = getCext(particle, alpha, k, b_sca, nmax)[0]
+ C_sca = getCsca(k, b_sca, nmax)[0]
+ print C_ext, C_sca
+ delta = (C_ext-C_sca)/(C_ext+C_sca)
+ print delta
+ self.assertAlmostEqual(delta,0,5)
+ self.assertAlmostEqual(C_ext,C_sca,5)
+ self.assertTrue(C_ext>0.)
+ self.assertTrue(C_sca>0.)
+
+ def test_IntegrationCase2(self):
+ nmax = 10
+ c1 = 1
+ c2 = 6
+ alpha = pi / 3
+ k = 1
+ particle = ProlateSpheroid(psi=2,c=c1,derivative=0,eps=2)
+ ebcm = SpheroidalEBCM(particle,c2,c1,nmax)
+ b_sca = ebcm.getMatrixSolution(TMInputWave(alpha))[0]
+ C_ext = getCext(particle, alpha, k, b_sca, nmax)[0]
+ C_sca = getCsca(k, b_sca, nmax)[0]
+ delta = (C_ext-C_sca)/(C_ext+C_sca)
+ print delta
+ self.assertAlmostEqual(delta,0,5)
+ self.assertAlmostEqual(C_ext,C_sca,5)
+ self.assertTrue(C_ext>0.)
+ self.assertTrue(C_sca>0.)
+
+if __name__ == '__main__':
+ #import cProfile
+ #cProfile.run('unittest.main()','profiler_results')
+ unittest.main()
+
+
+
View
366 test/spheroidal_matrix_test.py
@@ -4,6 +4,7 @@
#tests without using derivatives
#prolate spheroid
from spheroidal_particles import TMInputWave
+from spheroidal_svm import *
class Particle:
type = 1
@@ -21,244 +22,132 @@ def derivative(self, nu):
class testSpheroidalMatrixA(unittest.TestCase):
- def test_get_a_functionR1_1(self):
- n = 1;
- m = 0;
- c = 2.5;
- z = 0.5;
- type = 1;
- particle = Particle(2.5)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 1, particle)(z), -0.1280857374177537)
-
- def test_get_a_functionR1_2(self):
- n = 3;
- m = 2;
- c = 2.5;
- z = -0.5;
- type = 1;
- particle = Particle(1.5)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 1, particle)(z), -0.1547625029193905)
-
- def test_get_a_functionR1_3(self):
- n = 5;
- m = 0;
- c = 0.5;
- z = -0.9;
- type = 1;
- particle = Particle(5)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 1, particle)(z), 0.000642762007655166)
-
- def test_get_a_functionR1_4(self):
- n = 4;
- m = 2;
- c = 4.5;
- z = -0.6;
- type = 1;
- psi = 6
- particle = Particle(psi)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 1, particle)(z), 0.03417296826676850, 5)
-
- def test_get_a_functionR3_1(self):
- n = 1;
- m = 0;
- c = 2.5;
- z = 0.5;
- type = 1;
- particle = Particle(2.5)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 3, particle)(z),
- -0.1280857374177537 + 0.01321180926927677j)
-
- def test_get_a_functionR3_2(self):
- n = 3;
- m = 2;
- c = 2.5;
- z = -0.5;
- type = 1;
- particle = Particle(1.5)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 3, particle)(z), -0.154763 + 0.461291j, 5)
-
- def test_get_a_functionR3_3(self):
- n = 5;
- m = 0;
- c = 0.5;
- z = -0.9;
- type = 1;
- particle = Particle(5)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 3, particle)(z),
- 0.000642762007655166 - 0.539400871910197j)
-
- def test_get_a_functionR3_4(self):
- n = 4;
- m = 2;
- c = 4.5;
- z = -0.6;
- type = 1;
- psi = 6
- particle = Particle(psi)
- cv = get_cv(m, n, c, type)
- self.assertAlmostEqual(get_a_functions(m, n, c, cv, type, 3, particle)(z),
- 0.03417296826676850 + 0.01608088386790631j, 5)
-
- def test_get_A1(self):
- particle = Particle(5)
- c1 = 0.3;
- c2 = 0.4;
- nmax = 1
- self.assertAlmostEqual(get_A(particle, c2, c1, 1, nmax)[0][0], 0.432861, 5)
-
- def test_get_A2(self):
- particle = Particle(3)
- c1 = 2.5;
- c2 = 1.5;
- nmax = 2
- A = get_A(particle, c2, c1, 1, nmax)
- self.assertAlmostEqual(A[0][0], 0.0450505, 5)
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[1][1], 0.238367, 5)
-
- def test_get_A3(self):
- particle = Particle(1.5)
- c1 = 0.5;
- c2 = 3.5;
- nmax = 3
- A = get_A(particle, c2, c1, 1, nmax)
- self.assertAlmostEqual(A[0][0], 0.0800241, 5)
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[1][1], 0.229468, 5)
- self.assertAlmostEqual(A[2][2], 0.248043, 5)
- self.assertAlmostEqual(A[0][2], -0.0150865)
- self.assertAlmostEqual(A[2][0], 0.0474487, 5)
- self.assertAlmostEqual(A[2][1], 0)
- self.assertAlmostEqual(A[1][2], 0)
def test_getA11Simplest(self):
psi = 1.5
c = 2
particle = Particle(psi,d=1,eps=1)
nmax = 3
- A = get_A11(particle,c,nmax)
+ svm = SpheroidalSVM(particle,c,c,nmax)
+ A = svm.get_A11()
cv = get_cv(1,1,c,particle.type)
- self.assertAlmostEqual(A[0][0],rad3_cv(1,1,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[0,0],rad3_cv(1,1,c,cv,particle.type,psi)[0])
cv = get_cv(1,2,c,particle.type)
- self.assertAlmostEqual(A[1][1],rad3_cv(1,2,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[1,1],rad3_cv(1,2,c,cv,particle.type,psi)[0])
cv = get_cv(1,3,c,particle.type)
- self.assertAlmostEqual(A[2][2],rad3_cv(1,3,c,cv,particle.type,psi)[0])
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[0][2], 0)
- self.assertAlmostEqual(A[2][0], 0)
- self.assertAlmostEqual(A[2][1], 0)
- self.assertAlmostEqual(A[1][2], 0)
+ self.assertAlmostEqual(A[2,2],rad3_cv(1,3,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[1,0], 0)
+ self.assertAlmostEqual(A[0,1], 0)
+ self.assertAlmostEqual(A[0,2], 0)
+ self.assertAlmostEqual(A[2,0], 0)
+ self.assertAlmostEqual(A[2,1], 0)
+ self.assertAlmostEqual(A[1,2], 0)
def test_getA12Simplest(self):
psi = 3
c = 2.5
particle = Particle(psi,d=2.0,eps=1)
nmax = 3
- A = get_A12(particle,c,c,nmax)
+ svm = SpheroidalSVM(particle,c,c,nmax)
+ A = svm.get_A12()
cv = get_cv(1,1,c,particle.type)
- self.assertAlmostEqual(A[0][0],-rad1_cv(1,1,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[0,0],-rad1_cv(1,1,c,cv,particle.type,psi)[0])
cv = get_cv(1,2,c,particle.type)
- self.assertAlmostEqual(A[1][1],-rad1_cv(1,2,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[1,1],-rad1_cv(1,2,c,cv,particle.type,psi)[0])
cv = get_cv(1,3,c,particle.type)
- self.assertAlmostEqual(A[2][2],-rad1_cv(1,3,c,cv,particle.type,psi)[0])
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[0][2], 0)
- self.assertAlmostEqual(A[2][0], 0)
- self.assertAlmostEqual(A[2][1], 0)
- self.assertAlmostEqual(A[1][2], 0)
+ self.assertAlmostEqual(A[2,2],-rad1_cv(1,3,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[1,0], 0)
+ self.assertAlmostEqual(A[0,1], 0)
+ self.assertAlmostEqual(A[0,2], 0)
+ self.assertAlmostEqual(A[2,0], 0)
+ self.assertAlmostEqual(A[2,1], 0)
+ self.assertAlmostEqual(A[1,2], 0)
def test_getA21Simplest(self):
psi = 3
c = 2.5
particle = Particle(psi,d=2.0,eps=1)
nmax = 3
- A = get_A21(particle,c,nmax)
+ svm = SpheroidalSVM(particle,c,c,nmax)
+ A = svm.get_A21()
cv = get_cv(1,1,c,particle.type)
coeff = lambda nu: metric_phi(nu,particle) * metric_nu(nu,particle) / metric_psi(nu,particle) * ang1_cv(1, 1, c, cv, particle.type, nu)[0] * ang1_cv(1, 1, c, cv, particle.type, nu)[0]
coeff = quad(coeff,-1,1)
- self.assertAlmostEqual(A[0][0],rad3_cv(1,1,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[0,0],rad3_cv(1,1,c,cv,particle.type,psi)[1] * coeff)
cv = get_cv(1,2,c,particle.type)
- self.assertAlmostEqual(A[1][1],rad3_cv(1,2,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[1,1],rad3_cv(1,2,c,cv,particle.type,psi)[1] * coeff)
cv = get_cv(1,3,c,particle.type)
- self.assertAlmostEqual(A[2][2],rad3_cv(1,3,c,cv,particle.type,psi)[1] * coeff)
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[0][2], 0)
- self.assertAlmostEqual(A[2][0], 0)
- self.assertAlmostEqual(A[2][1], 0)
- self.assertAlmostEqual(A[1][2], 0)
+ self.assertAlmostEqual(A[2,2],rad3_cv(1,3,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[1,0], 0)
+ self.assertAlmostEqual(A[0,1], 0)
+ self.assertAlmostEqual(A[0,2], 0)
+ self.assertAlmostEqual(A[2,0], 0)
+ self.assertAlmostEqual(A[2,1], 0)
+ self.assertAlmostEqual(A[1,2], 0)
def test_getA22Simplest(self):
psi = 3
c = 2.5
particle = Particle(psi,d=2.0,eps=1)
nmax = 3
- A = get_A22(particle,c,c,nmax)
+ svm = SpheroidalSVM(particle,c,c,nmax)
+ A = svm.get_A22()
cv = get_cv(1,1,c,particle.type)
coeff = lambda nu: metric_phi(nu,particle) * metric_nu(nu,particle) / metric_psi(nu,particle) * ang1_cv(1, 1, c, cv, particle.type, nu)[0] * ang1_cv(1, 1, c, cv, particle.type, nu)[0]
coeff = - quad(coeff,-1,1)
- self.assertAlmostEqual(A[0][0],rad1_cv(1,1,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[0,0],rad1_cv(1,1,c,cv,particle.type,psi)[1] * coeff)
cv = get_cv(1,2,c,particle.type)
- self.assertAlmostEqual(A[1][1],rad1_cv(1,2,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[1,1],rad1_cv(1,2,c,cv,particle.type,psi)[1] * coeff)
cv = get_cv(1,3,c,particle.type)
- self.assertAlmostEqual(A[2][2],rad1_cv(1,3,c,cv,particle.type,psi)[1] * coeff)
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[0][2], 0)
- self.assertAlmostEqual(A[2][0], 0)
- self.assertAlmostEqual(A[2][1], 0)
- self.assertAlmostEqual(A[1][2], 0)
+ self.assertAlmostEqual(A[2,2],rad1_cv(1,3,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[1,0], 0)
+ self.assertAlmostEqual(A[0,1], 0)
+ self.assertAlmostEqual(A[0,2], 0)
+ self.assertAlmostEqual(A[2,0], 0)
+ self.assertAlmostEqual(A[2,1], 0)
+ self.assertAlmostEqual(A[1,2], 0)
def test_getA10Simplest(self):
psi = 3
c = 2.5
particle = Particle(psi,d=2.0,eps=1)
nmax = 3
- A = get_A10(particle,c,nmax)
+ svm = SpheroidalSVM(particle,c,c,nmax)
+ A = svm.get_A10()
cv = get_cv(1,1,c,particle.type)
- self.assertAlmostEqual(A[0][0],-rad1_cv(1,1,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[0,0],-rad1_cv(1,1,c,cv,particle.type,psi)[0])
cv = get_cv(1,2,c,particle.type)
- self.assertAlmostEqual(A[1][1],-rad1_cv(1,2,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[1,1],-rad1_cv(1,2,c,cv,particle.type,psi)[0])
cv = get_cv(1,3,c,particle.type)
- self.assertAlmostEqual(A[2][2],-rad1_cv(1,3,c,cv,particle.type,psi)[0])
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[0][2], 0)
- self.assertAlmostEqual(A[2][0], 0)
- self.assertAlmostEqual(A[2][1], 0)
- self.assertAlmostEqual(A[1][2], 0)
+ self.assertAlmostEqual(A[2,2],-rad1_cv(1,3,c,cv,particle.type,psi)[0])
+ self.assertAlmostEqual(A[1,0], 0)
+ self.assertAlmostEqual(A[0,1], 0)
+ self.assertAlmostEqual(A[0,2], 0)
+ self.assertAlmostEqual(A[2,0], 0)
+ self.assertAlmostEqual(A[2,1], 0)
+ self.assertAlmostEqual(A[1,2], 0)
def test_getA20Simplest(self):
psi = 3
c = 2.5
particle = Particle(psi,d=2.0,eps=1)
nmax = 3
- A = get_A20(particle,c,nmax)
+ svm = SpheroidalSVM(particle,c,c,nmax)
+ A = svm.get_A20()
cv = get_cv(1,1,c,particle.type)
coeff = lambda nu: metric_phi(nu,particle) * metric_nu(nu,particle) / metric_psi(nu,particle) * ang1_cv(1, 1, c, cv, particle.type, nu)[0] * ang1_cv(1, 1, c, cv, particle.type, nu)[0]
coeff = - quad(coeff,-1,1)
- self.assertAlmostEqual(A[0][0],rad1_cv(1,1,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[0,0],rad1_cv(1,1,c,cv,particle.type,psi)[1] * coeff)
cv = get_cv(1,2,c,particle.type)
- self.assertAlmostEqual(A[1][1],rad1_cv(1,2,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[1,1],rad1_cv(1,2,c,cv,particle.type,psi)[1] * coeff)
cv = get_cv(1,3,c,particle.type)
- self.assertAlmostEqual(A[2][2],rad1_cv(1,3,c,cv,particle.type,psi)[1] * coeff)
- self.assertAlmostEqual(A[1][0], 0)
- self.assertAlmostEqual(A[0][1], 0)
- self.assertAlmostEqual(A[0][2], 0)
- self.assertAlmostEqual(A[2][0], 0)
- self.assertAlmostEqual(A[2][1], 0)
- self.assertAlmostEqual(A[1][2], 0)
+ self.assertAlmostEqual(A[2,2],rad1_cv(1,3,c,cv,particle.type,psi)[1] * coeff)
+ self.assertAlmostEqual(A[1,0], 0)
+ self.assertAlmostEqual(A[0,1], 0)
+ self.assertAlmostEqual(A[0,2], 0)
+ self.assertAlmostEqual(A[2,0], 0)
+ self.assertAlmostEqual(A[2,1], 0)
+ self.assertAlmostEqual(A[1,2], 0)
def test_get_fullA1(self):
c1 = 0.5;
@@ -267,11 +156,12 @@ def test_get_fullA1(self):
d = 1.5;
psi = 1.5
particle = Particle(psi, d)
- A = get_fullA(particle, c1, c2, nmax)
- self.assertAlmostEqual(A[0, 0], get_A(particle, c1, c1, 3, nmax))
- self.assertAlmostEqual(A[0, 1], -get_A(particle, c2, c1, 1, nmax))
- self.assertAlmostEqual(A[1, 0], get_B(particle, c1, c1, 3, nmax))
- self.assertAlmostEqual(A[1, 1], -get_C(particle, c2, c1, nmax))
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ A = svm.get_fullA()
+ self.assertAlmostEqual(A[0, 0], svm.get_A11())
+ self.assertAlmostEqual(A[0, 1], svm.get_A12())
+ self.assertAlmostEqual(A[1, 0], svm.get_A21())
+ self.assertAlmostEqual(A[1, 1], svm.get_A22())
def test_get_fullA2(self):
c1 = 0.5;
@@ -280,11 +170,12 @@ def test_get_fullA2(self):
d = 1.5;
psi = 1.5
particle = Particle(psi, d)
- A = get_fullA(particle, c1, c2, nmax)
- A11 = get_A(particle, c1, c1, 3, nmax)
- A12 = -get_A(particle, c2, c1, 1, nmax)
- A21 = get_B(particle, c1, c1, 3, nmax)
- A22 = -get_C(particle, c2, c1, nmax)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ A = svm.get_fullA()
+ A11 = svm.get_A11()
+ A12 = svm.get_A12()
+ A21 = svm.get_A21()
+ A22 = svm.get_A22()
self.assertAlmostEqual(A[0, 0], A11[0, 0])
self.assertAlmostEqual(A[0, 1], A11[0, 1])
self.assertAlmostEqual(A[1, 0], A11[1, 0])
@@ -305,55 +196,16 @@ def test_get_fullA2(self):
class testSpheroidalMatrixZ(unittest.TestCase):
- def test_get_B1(self):
- d = 8;
- psi = 6
- c1 = 0.3;
- c2 = 0.4;
- nmax = 1
- particle = Particle(psi, d)
- self.assertAlmostEqual(get_B(particle, c2, c1, 1, nmax)[0][0], -3.69601, 5)
-
- def test_get_B2(self):
- d = 2;
- psi = 3
- c1 = 4;
- c2 = 5;
- nmax = 2
- particle = Particle(psi, d)
- B = get_B(particle, c2, c1, 1, nmax)
- self.assertAlmostEqual(B[0][0], 2.79893, 5)
- self.assertAlmostEqual(B[1][0], 0, 5)
- self.assertAlmostEqual(B[0][1], 0, 5)
- self.assertAlmostEqual(B[1][1], 1.55881, 5)
-
- def test_get_B3(self):
- d = 1;
- psi = 10
- c1 = 0.5;
- c2 = 5;
- nmax = 3
- particle = Particle(psi, d)
- B = get_B(particle, c2, c1, 1, nmax)
- self.assertAlmostEqual(B[0][0], -2.040523, 4)
- self.assertAlmostEqual(B[1][0], 0, 5)
- self.assertAlmostEqual(B[0][1], 0, 5)
- self.assertAlmostEqual(B[1][1], -4.41099, 5)
- self.assertAlmostEqual(B[1][2], 0, 5)
- self.assertAlmostEqual(B[2][1], 0, 5)
- self.assertAlmostEqual(B[0][2], 0.618216, 5)
- self.assertAlmostEqual(B[2][0], 0.373256, 5)
- self.assertAlmostEqual(B[2][2], 1.16732, 5)
-
def test_getFullB1(self):
c1 = 4.5;
nmax = 1;
d = 2.5;
psi = 3.5
particle = Particle(psi, d)
- B = get_fullB(particle, c1, nmax)
- self.assertAlmostEqual(B[0], -get_A(particle, c1, c1, 1, nmax))
- self.assertAlmostEqual(B[1], -get_B(particle, c1, c1, 1, nmax))
+ svm = SpheroidalSVM(particle,c1,c1,nmax)
+ B = svm.get_fullB()
+ self.assertAlmostEqual(B[0], svm.get_A10())
+ self.assertAlmostEqual(B[1], svm.get_A20())
def test_getFullB2(self):
c1 = 4.5;
@@ -361,9 +213,10 @@ def test_getFullB2(self):
d = 2.5;
psi = 3.5
particle = Particle(psi, d)
- B = get_fullB(particle, c1, nmax)
- B1 = -get_A(particle, c1, c1, 1, nmax)
- B2 = -get_B(particle, c1, c1, 1, nmax)
+ svm = SpheroidalSVM(particle,c1,c1,nmax)
+ B = svm.get_fullB()
+ B1 = svm.get_A10()
+ B2 = svm.get_A20()
self.assertAlmostEqual(B[0, 0], B1[0, 0])
self.assertAlmostEqual(B[0, 1], B1[0, 1])
self.assertAlmostEqual(B[1, 0], B1[1, 0])
@@ -371,47 +224,4 @@ def test_getFullB2(self):
self.assertAlmostEqual(B[2, 0], B2[0, 0])
self.assertAlmostEqual(B[2, 1], B2[0, 1])
self.assertAlmostEqual(B[3, 0], B2[1, 0])
- self.assertAlmostEqual(B[3, 1], B2[1, 1])
-
- def test_get_C1(self):
- d = 8;
- eps = 0.7;
- psi = 6
- c1 = 0.3;
- c2 = 0.4;
- nmax = 1
- particle = Particle(psi, d, eps)
- self.assertAlmostEqual(get_C(particle, c2, c1, nmax)[0][0], -9.65321, 5)
-
- def test_get_C2(self):
- d = 2;
- psi = 3;
- eps = 5
- c1 = 4;
- c2 = 5;
- nmax = 2
- particle = Particle(psi, d, eps)
- C = get_C(particle, c2, c1, nmax)
- self.assertAlmostEqual(C[0][0], 0.590814, 5)
- self.assertAlmostEqual(C[1][0], 0, 5)
- self.assertAlmostEqual(C[0][1], 0, 5)
- self.assertAlmostEqual(C[1][1], 0.166018, 5)
-
- def test_get_C3(self):
- d = 1;
- eps = 0.3;
- psi = 10
- c1 = 0.5;
- c2 = 5;
- nmax = 3
- particle = Particle(psi, d, eps)
- C = get_C(particle, c2, c1, nmax)
- self.assertAlmostEqual(C[0][0], -6.60149, 4)
- self.assertAlmostEqual(C[1][0], 0, 5)
- self.assertAlmostEqual(C[0][1], 0, 5)
- self.assertAlmostEqual(C[1][1], -14.78301, 5)
- self.assertAlmostEqual(C[1][2], 0, 5)
- self.assertAlmostEqual(C[2][1], 0, 5)
- self.assertAlmostEqual(C[0][2], 2.00006, 5)
- self.assertAlmostEqual(C[2][0], 1.17881, 5)
- self.assertAlmostEqual(C[2][2], 3.68663, 5)
+ self.assertAlmostEqual(B[3, 1], B2[1, 1])
View
1  test/spheroidal_test.py
@@ -1,5 +1,6 @@
import numpy
import unittest
+import scipy.linalg
from scipy.special import *
View
118 test/svm_test.py
@@ -1,5 +1,7 @@
import unittest
+from numpy import transpose
+
from spheroidal import *
from spheroidal_particles import *
from properties import *
@@ -7,17 +9,121 @@
#main integration tests for non-absorbing particle
class testSVM(unittest.TestCase):
- def test_IntegrationCase(self):
- nmax = 10
+ #according to (100)
+ #holds if k=1 so d should be 2 * c1
+ def test_relation1(self):
+ nmax = 6
+ c1 = 1
+ c2 = 2
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=1)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ relation = svm.get_B(3) * svm.get_A(c1,c1,1).T - svm.get_A(c1,c1,3) * svm.get_B(1).T
+ right_part = 1j*mat(eye(nmax))
+ for i in range(0,nmax):
+ for k in range(0,nmax):
+ self.assertAlmostEqual(relation[i,k],right_part[i,k])
+
+ def test_relation2(self):
+ nmax = 6
+ c1 = 1.3
+ c2 = 2
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=1)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ relation = svm.get_B(3) * svm.get_A(c1,c1,3).T - svm.get_A(c1,c1,3) * svm.get_B(3).T
+ for i in range(0,nmax):
+ for k in range(0,nmax):
+ self.assertAlmostEqual(relation[i,k],0)
+
+ def test_relation3(self):
+ nmax = 6
+ c1 = 1.3
+ c2 = 2
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=1)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ relation = svm.get_B(1) * svm.get_A(c1,c1,1).T - svm.get_A(c1,c1,1) * svm.get_B(1).T
+ for i in range(0,nmax):
+ for k in range(0,nmax):
+ self.assertAlmostEqual(relation[i,k],0)
+
+ def test_relation4(self):
+ nmax = 6
+ c1 = 1
+ c2 = 2
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=1)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ relation = svm.get_A(c1,c1,1).T * svm.get_B(3) - svm.get_A(c1,c1,3).T * svm.get_B(1)
+ right_part = 1j*mat(eye(nmax))
+ for i in range(0,nmax):
+ for k in range(0,nmax):
+ self.assertAlmostEqual(relation[i,k],right_part[i,k])
+
+ def test_relation5(self):
+ nmax = 6
+ c1 = 1.3
+ c2 = 2
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=1)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ relation = svm.get_A(c1,c1,3).T * svm.get_A(c1,c1,1) - svm.get_A(c1,c1,1).T * svm.get_A(c1,c1,3)
+ for i in range(0,nmax):
+ for k in range(0,nmax):
+ self.assertAlmostEqual(relation[i,k], 0)
+
+ def test_relation6(self):
+ nmax = 6
c1 = 1.3
c2 = 2
- alpha = pi / 4
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=1)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ relation = svm.get_B(3).T * svm.get_B(1) - svm.get_B(1).T * svm.get_B(3)
+ for i in range(0,nmax):
+ for k in range(0,nmax):
+ self.assertAlmostEqual(relation[i,k], 0)
+
+ def test_relation7(self):
+ nmax = 6
+ c1 = 1
+ c2 = 2
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=1)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ relation = svm.get_A(c1,c1,1) * svm.get_B(3).T - svm.get_B(1).T * svm.get_A(c1,c1,3)
+ right_part = 1j*mat(eye(nmax))
+ for i in range(0,nmax):
+ for k in range(0,nmax):
+ self.assertAlmostEqual(relation[i,k], right_part[i,k])
+
+
+ def test_IntegrationCase1(self):
+ nmax = 6
+ c1 = 1
+ c2 = 2
+ alpha = pi / 3
+ k = 1
+ particle = ProlateSpheroid(psi=2,c=c1,derivative=0,eps=2)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ b_sca = svm.getSolution(TMInputWave(alpha))[0]
+ C_ext = getCext(particle, alpha, k, b_sca, nmax)[0]
+ C_sca = getCsca(k, b_sca, nmax)[0]
+ print C_ext, C_sca
+ delta = (C_ext-C_sca)/(C_ext+C_sca)
+ print delta
+ self.assertAlmostEqual(delta,0,5)
+ self.assertAlmostEqual(C_ext,C_sca,5)
+ self.assertTrue(C_ext>0.)
+ self.assertTrue(C_sca>0.)
+
+ def test_IntegrationCase2(self):
+ nmax = 10
+ c1 = 1
+ c2 = 6
+ alpha = pi / 3
k = 1
- particle = ProlateSpheroid(psi=2,c=c1,derivative=0,eps=1)
- b_sca = getSolution(particle, TMInputWave(alpha), c1, c2, nmax)[0]
- C_ext = -getCext(particle, alpha, k, b_sca, nmax)[0]
+ particle = OblateSpheroid(psi=2,c=c1,derivative=0,eps=2)
+ svm = SpheroidalSVM(particle,c2,c1,nmax)
+ b_sca = svm.getSolution(TMInputWave(alpha))[0]
+ C_ext = getCext(particle, alpha, k, b_sca, nmax)[0]
C_sca = getCsca(k, b_sca, nmax)[0]
delta = (C_ext-C_sca)/(C_ext+C_sca)
+ print delta
self.assertAlmostEqual(delta,0,5)
self.assertAlmostEqual(C_ext,C_sca,5)
self.assertTrue(C_ext>0.)

No commit comments for this range

Something went wrong with that request. Please try again.