Fixed issues associated with Python2

Also decided to include a copy of PyMieScatt with the package.  This is
because it had problems with Python 2 support and also the default accuracy settings will have to be changed

BuildDocs/phonana.md

@@ -255,7 +255,8 @@ Phonana doesn't need any extra installation above what is needed to run PDielec,
 
     sudo pip install vtk
     sudo pip install pyvtk
-    sudo pip instal  vtkinterface
+    sudo pip install  vtkinterface
+    sudo pip install  mayavi
 
 The installation of these modules can be performed as root, if all users on the machine require access to them, or they can be installed in the user's file system using 'pip --user' instead of 'pip'
 

Python/Calculator.py

@@ -465,9 +465,10 @@ def mie_scattering(dielectric_medium, dielecv, shape, L, vf, size, size_mu, size
        vf is the volume fraction of filler
        Mie only works for spherical particles, so shape, and L parameters are ignored
        The routine returns the effective dielectric constant"""
-    import PyMieScatt as ps
+    # import Python.PyMieScatt as ps
     from scipy.integrate import trapz
     from scipy.stats import lognorm
+    from Python.PyMieScatt.Mie import MieS1S2
     # define i as a complex number
     i = complex(0,1)
     # We need to taken account of the change in wavelength and the change in size parameter due to the 
@@ -522,7 +523,7 @@ def mie_scattering(dielectric_medium, dielecv, shape, L, vf, size, size_mu, size
                 # The size parameter is 2pi r / lambda 
                 x = 2 * PI * r / lambda_vacuum_mu
                 # Calculate the S1 and S2 scattering factors, and store in a list
-                s1,s2 = ps.MieS1S2(refractive_index, x*refractive_index_medium, 1)
+                s1,s2 = MieS1S2(refractive_index, x*refractive_index_medium, 1)
                 s1_factors.append(s1)
             # Now integrate
             s1 = trapz(s1_factors*ndp,dp)
@@ -539,7 +540,7 @@ def mie_scattering(dielectric_medium, dielecv, shape, L, vf, size, size_mu, size
         else:
             # Calculate the scattering factors at 0 degrees 
             #jk print("refractive_index, size, refractive_index_medium", refractive_index, size, refractive_index_medium)
-            s1,s2 = ps.MieS1S2(refractive_index, size*refractive_index_medium, 1)
+            s1,s2 = MieS1S2(refractive_index, size*refractive_index_medium, 1)
         # qext,qsca,qabs,g,qpr,qback,qratio = ps.AutoMieQ(refractive_index, wavelength_nm, diameter_nm)
         #jk print('s1,s2',s1,s2)
         #jk print('qext,qsca,qabs',qext,qsca,qabs)

Python/PyMieScatt/CoreShell.py

@@ -0,0 +1,139 @@
+# -*- coding: utf-8 -*-
+# http://pymiescatt.readthedocs.io/en/latest/forwardCS.html
+import numpy as np
+from scipy.special import jv, yv
+from PyMieScatt.Mie import MieQ, MiePiTau
+
+def MieQCoreShell(mCore,mShell,wavelength,dCore,dShell):
+#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#MieQCoreShell
+  xCore = np.pi*dCore/wavelength
+  xShell = np.pi*dShell/wavelength
+  if xCore==xShell:
+    return MieQ(mCore,wavelength,dShell)
+  elif xCore==0:
+    return MieQ(mShell,wavelength,dShell)
+  elif mCore==mShell:
+    return MieQ(mCore,wavelength,dShell)
+  elif xCore>0:
+    nmax = np.round(2+xShell+4*(xShell**(1/3)))
+    n = np.arange(1,nmax+1)
+    n1 = 2*n+1
+    n2 = n*(n+2)/(n+1)
+    n3 = n1/(n*(n+1))
+    xShell2 = xShell**2
+
+    an, bn = CoreShell_ab(mCore,mShell,xCore,xShell)
+
+    qext = (2/xShell2)*np.sum(n1*(an.real+bn.real))
+    qsca = (2/xShell2)*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*xShell2))*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/xShell2)*(np.abs(np.sum(n1*((-1)**n)*(an-bn)))**2)
+    qratio = qback/qsca
+
+    return qext, qsca, qabs, g, qpr, qback, qratio
+
+def CoreShell_ab(mCore,mShell,xCore,xShell):
+#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShell_ab
+  m = mShell/mCore
+  u = mCore*xCore
+  v = mShell*xCore
+  w = mShell*xShell
+
+  mx = max(np.abs(mCore*xShell),np.abs(mShell*xShell))
+  nmax = np.round(2+xShell+4*(xShell**(1/3)))
+  nmx = np.round(max(nmax,mx)+16)
+  n = np.arange(1,nmax+1)
+  nu = n+0.5
+
+  sv = np.sqrt(0.5*np.pi*v)
+  sw = np.sqrt(0.5*np.pi*w)
+  sy = np.sqrt(0.5*np.pi*xShell)
+
+  pv = sv*jv(nu,v)
+  pw = sw*jv(nu,w)
+  py = sy*jv(nu,xShell)
+
+  chv = -sv*yv(nu,v)
+  chw = -sw*yv(nu,w)
+  chy = -sy*yv(nu,xShell)
+
+  p1y = np.append([np.sin(xShell)], [py[0:int(nmax)-1]])
+  ch1y = np.append([np.cos(xShell)], [chy[0:int(nmax)-1]])
+  gsy = py-(0+1.0j)*chy
+  gs1y = p1y-(0+1.0j)*ch1y
+
+  # B&H Equation 4.89
+  Dnu = np.zeros((int(nmx)),dtype=complex)
+  Dnv = np.zeros((int(nmx)),dtype=complex)
+  Dnw = np.zeros((int(nmx)),dtype=complex)
+  for i in range(int(nmx)-1,1,-1):
+    Dnu[i-1] = i/u-1/(Dnu[i]+i/u)
+    Dnv[i-1] = i/v-1/(Dnv[i]+i/v)
+    Dnw[i-1] = i/w-1/(Dnw[i]+i/w)
+
+  Du = Dnu[1:int(nmax)+1]
+  Dv = Dnv[1:int(nmax)+1]
+  Dw = Dnw[1:int(nmax)+1]
+
+  uu = m*Du-Dv
+  vv = Du/m-Dv
+  fv = pv/chv
+
+  dns = ((uu*fv/pw)/(uu*(pw-chw*fv)+(pw/pv)/chv))+Dw
+  gns = ((vv*fv/pw)/(vv*(pw-chw*fv)+(pw/pv)/chv))+Dw
+  a1 = dns/mShell+n/xShell
+  b1 = mShell*gns+n/xShell
+
+  an = (py*a1-p1y)/(gsy*a1-gs1y)
+  bn = (py*b1-p1y)/(gsy*b1-gs1y)
+
+  return an, bn
+
+def CoreShellScatteringFunction(mCore,mShell,wavelength,dCore,dShell,minAngle=0, maxAngle=180, angularResolution=0.5, normed=False):
+#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellScatteringFunction
+  xCore = np.pi*dCore/wavelength
+  xShell = np.pi*dShell/wavelength
+  theta = np.linspace(minAngle,maxAngle,int((maxAngle-minAngle)/angularResolution))*np.pi/180
+  thetaSteps = len(theta)
+  SL = np.zeros(thetaSteps)
+  SR = np.zeros(thetaSteps)
+  SU = np.zeros(thetaSteps)
+  for j in range(thetaSteps):
+    u = np.cos(theta[j])
+    S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,u)
+    SL[j] = (np.sum((np.conjugate(S1)*S1))).real
+    SR[j] = (np.sum((np.conjugate(S2)*S2))).real
+    SU[j] = (SR[j]+SL[j])/2
+  if normed:
+    SL /= np.max(SL)
+    SR /= np.max(SR)
+    SU /= np.max(SU)
+  return theta,SL,SR,SU
+
+def CoreShellS1S2(mCore,mShell,xCore,xShell,mu):
+#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellS1S2
+  nmax = np.round(2+xShell+4*(xShell**(1/3)))
+  an,bn = CoreShell_ab(mCore,mShell,xCore,xShell)
+  pin,taun = MiePiTau(mu,nmax)
+  n = np.arange(1,int(nmax)+1)
+  n2 = (2*n+1)/(n*(n+1))
+  pin *= n2
+  taun *= n2
+  S1=np.sum(an*np.conjugate(pin))+np.sum(bn*np.conjugate(taun))
+  S2=np.sum(an*np.conjugate(taun))+np.sum(bn*np.conjugate(pin))
+  return S1,S2
+
+def CoreShellMatrixElements(mCore,mShell,xCore,xShell,mu):
+#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellMatrixElements
+  S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,mu)
+  S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
+  S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
+  S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
+  S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
+  return S11, S12, S33, S34