fieldplot.py 5.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168
  1. #!/usr/bin/env python
  2. # -*- coding: UTF-8 -*-
  3. #
  4. # Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
  5. #
  6. # This file is part of python-scattnlay
  7. #
  8. # This program is free software: you can redistribute it and/or modify
  9. # it under the terms of the GNU General Public License as published by
  10. # the Free Software Foundation, either version 3 of the License, or
  11. # (at your option) any later version.
  12. #
  13. # This program is distributed in the hope that it will be useful,
  14. # but WITHOUT ANY WARRANTY; without even the implied warranty of
  15. # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16. # GNU General Public License for more details.
  17. #
  18. # The only additional remark is that we expect that all publications
  19. # describing work using this software, or all commercial products
  20. # using it, cite the following reference:
  21. # [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
  22. # a multilayered sphere," Computer Physics Communications,
  23. # vol. 180, Nov. 2009, pp. 2348-2354.
  24. #
  25. # You should have received a copy of the GNU General Public License
  26. # along with this program. If not, see <http://www.gnu.org/licenses/>.
  27. # Several functions to plot field and streamlines (power flow lines).
  28. import scattnlay
  29. from scattnlay import fieldnlay
  30. from scattnlay import scattnlay
  31. import numpy as np
  32. import cmath
  33. def unit_vector(vector):
  34. """ Returns the unit vector of the vector. """
  35. return vector / np.linalg.norm(vector)
  36. def angle_between(v1, v2):
  37. """ Returns the angle in radians between vectors 'v1' and 'v2'::
  38. >>> angle_between((1, 0, 0), (0, 1, 0))
  39. 1.5707963267948966
  40. >>> angle_between((1, 0, 0), (1, 0, 0))
  41. 0.0
  42. >>> angle_between((1, 0, 0), (-1, 0, 0))
  43. 3.141592653589793
  44. """
  45. v1_u = unit_vector(v1)
  46. v2_u = unit_vector(v2)
  47. angle = np.arccos(np.dot(v1_u, v2_u))
  48. if np.isnan(angle):
  49. if (v1_u == v2_u).all():
  50. return 0.0
  51. else:
  52. return np.pi
  53. return angle
  54. ###############################################################################
  55. def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m):
  56. # Initial position
  57. flow_x = [x0]
  58. flow_y = [y0]
  59. flow_z = [z0]
  60. max_step = x[-1]/3
  61. min_step = x[0]/2000
  62. # max_step = min_step
  63. step = min_step*2.0
  64. if max_step < min_step:
  65. max_step = min_step
  66. coord = np.vstack(([flow_x[-1]], [flow_y[-1]], [flow_z[-1]])).transpose()
  67. terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord)
  68. Ec, Hc = E[0, 0, :], H[0, 0, :]
  69. S = np.cross(Ec, Hc.conjugate()).real
  70. Snorm_prev = S/np.linalg.norm(S)
  71. length = 0
  72. dpos = step
  73. while length < max_length:
  74. if step<max_step:
  75. step = step*2.0
  76. r = np.sqrt(flow_x[-1]**2 + flow_y[-1]**2 + flow_z[-1]**2)
  77. while step > min_step:
  78. #Evaluate displacement from previous poynting vector
  79. dpos = step
  80. dx = dpos*Snorm_prev[0];
  81. dy = dpos*Snorm_prev[1];
  82. dz = dpos*Snorm_prev[2];
  83. #Test the next position not to turn more than max_angle
  84. coord = np.vstack(([flow_x[-1]+dx], [flow_y[-1]+dy], [flow_z[-1]+dz])).transpose()
  85. terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord)
  86. Ec, Hc = E[0, 0, :], H[0, 0, :]
  87. Eth = max(np.absolute(Ec))/1e10
  88. Hth = max(np.absolute(Hc))/1e10
  89. for i in xrange(0,len(Ec)):
  90. if abs(Ec[i]) < Eth:
  91. Ec[i] = 0+0j
  92. if abs(Hc[i]) < Hth:
  93. Hc[i] = 0+0j
  94. S = np.cross(Ec, Hc.conjugate()).real
  95. Snorm = S/np.linalg.norm(S)
  96. #Snorm = Snorm.real
  97. #Snorm = np.absolute(Snorm)
  98. diff = Snorm-Snorm_prev
  99. if np.linalg.norm(diff)<0.05:
  100. # angle = angle_between(Snorm, Snorm_prev)
  101. # if abs(angle) < max_angle:
  102. break
  103. step = step/2.0
  104. #3. Save result
  105. Snorm_prev = Snorm
  106. dx = dpos*Snorm_prev[0];
  107. dy = dpos*Snorm_prev[1];
  108. dz = dpos*Snorm_prev[2];
  109. length = length + step
  110. flow_x.append(flow_x[-1] + dx)
  111. flow_y.append(flow_y[-1] + dy)
  112. flow_z.append(flow_z[-1] + dz)
  113. return np.array(flow_x), np.array(flow_y), np.array(flow_z)
  114. ###############################################################################
  115. def GetField(crossplane, npts, factor, x, m):
  116. """
  117. crossplane: XZ, YZ, XY
  118. npts: number of point in each direction
  119. factor: ratio of plotting size to outer size of the particle
  120. x: size parameters for particle layers
  121. m: relative index values for particle layers
  122. """
  123. scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
  124. zero = np.zeros(npts*npts, dtype = np.float64)
  125. if crossplane=='XZ':
  126. coordX, coordZ = np.meshgrid(scan, scan)
  127. coordX.resize(npts*npts)
  128. coordZ.resize(npts*npts)
  129. coordY = zero
  130. coordPlot1 = coordX
  131. coordPlot2 = coordZ
  132. elif crossplane=='YZ':
  133. coordY, coordZ = np.meshgrid(scan, scan)
  134. coordY.resize(npts*npts)
  135. coordZ.resize(npts*npts)
  136. coordX = zero
  137. coordPlot1 = coordY
  138. coordPlot2 = coordZ
  139. elif crossplane=='XY':
  140. coordX, coordY = np.meshgrid(scan, scan)
  141. coordX.resize(npts*npts)
  142. coordY.resize(npts*npts)
  143. coordZ = zero
  144. coordPlot1 = coordY
  145. coordPlot2 = coordX
  146. coord = np.vstack((coordX, coordY, coordZ)).transpose()
  147. terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord)
  148. Ec = E[0, :, :]
  149. Hc = H[0, :, :]
  150. P=[]
  151. P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real, range(0, len(E[0]))))
  152. # for n in range(0, len(E[0])):
  153. # P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
  154. return Ec, Hc, P, coordPlot1, coordPlot2