Skip to content

Commit

Permalink
finish tests for Mie simulation (close #1)
Browse files Browse the repository at this point in the history
  • Loading branch information
paulmueller committed Dec 20, 2017
1 parent c9ba15c commit 25efd67
Show file tree
Hide file tree
Showing 5 changed files with 124 additions and 48 deletions.
69 changes: 34 additions & 35 deletions qpsphere/models/_bhfield/wrap.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
"""BHFIELD sphere wrapper"""
import numbers
import os
import pathlib
import subprocess as sp
Expand All @@ -16,7 +15,7 @@ class BHFIELDExecutionError(BaseException):
pass


def clear_temp(wdir=".", rmdir=True):
def clear_temp(wdir, rmdir=True):
"""Remove all files in wdir"""
wdir = pathlib.Path(wdir)
extensions = ["*.log", "*.dat", "*.txt"]
Expand All @@ -29,9 +28,20 @@ def clear_temp(wdir=".", rmdir=True):
wdir.rmdir()


def load_field(wdir=".", size_grid=None):
def load_field(wdir, shape_grid):
"""Extract electric field from simulation file "V_0Ereim.dat"
Parameters
----------
wdir: str or pathlib.Path
path to the working directory
shape_grid: tuple of ints
the shape of the simulation data grid
Notes
-----
These are the files present in the working directory
bhfield.log:
log file containing scattering coeffs etc
(as well as definitions of the output fields)
Expand Down Expand Up @@ -95,33 +105,22 @@ def load_field(wdir=".", size_grid=None):

a = np.loadtxt(str(field_file))

if size_grid is None:
x = a[:, 0]
# y = a[:,1]
# z = a[:,2]
size_grid = int(np.sqrt(len(x)))

if isinstance(size_grid, numbers.Real):
size_grid = np.array((size_grid, size_grid), dtype=int)

assert size_grid[0] == int(
size_grid[0]), "resulting x-size is not an integer"
assert size_grid[1] == int(
size_grid[1]), "resulting y-size is not an integer"

size_grid = np.array(size_grid, dtype=int)
assert shape_grid[0] == int(
shape_grid[0]), "resulting x-size is not an integer"
assert shape_grid[1] == int(
shape_grid[1]), "resulting y-size is not an integer"

Ex = a[:, 3] + 1j * a[:, 6]
# Ey = a[:,4] + 1j*a[:,7]
# Ez = a[:,5] + 1j*a[:,8]

Exend = Ex.reshape((size_grid[1], size_grid[0])).transpose()
Exend = Ex.reshape((shape_grid[1], shape_grid[0])).transpose()
return Exend


def simulate_sphere(radius_sphere_um=2.5,
size_simulation_um=(7, 7),
size_grid=(50, 50),
size_simulation_um=[7, 7],
shape_grid=(50, 50),
refractive_index_medium=1.0,
refractive_index_sphere=1.01,
measurement_position_um=2.5,
Expand All @@ -134,11 +133,11 @@ def simulate_sphere(radius_sphere_um=2.5,
----------
radius_sphere_um : float
radius of sphere in um
size_simulation_um : tuple of floats
size_simulation_um : list of floats
Size of simulation volume in lateral dimension in um.
If a float is given, then a square simulation size is assumed. If
a tuple is given, then a rectangular shape is assumed.
size_grid : tuple of ints
shape_grid : tuple of ints
grid points in each lateral dimension.
If a float is given, then a square simulation size is assumed. If
a tuple is given, then a rectangular shape is assumed.
Expand All @@ -162,15 +161,15 @@ def simulate_sphere(radius_sphere_um=2.5,

# The size of the simulation must be zero
# if there is only one grid point.
if size_grid[0] == 1:
if shape_grid[0] == 1:
sizeum[0] = 0
if size_grid[1] == 1:
if shape_grid[1] == 1:
sizeum[1] = 0

assert np.allclose(size_grid[0], int(
size_grid[0])), "resulting x-size is not an integer"
assert np.allclose(size_grid[1], int(
size_grid[1])), "resulting y-size is not an integer"
assert np.allclose(shape_grid[0], int(
shape_grid[0])), "resulting x-size is not an integer"
assert np.allclose(shape_grid[1], int(
shape_grid[1])), "resulting y-size is not an integer"

while True:
# create temp dir
Expand All @@ -182,10 +181,10 @@ def simulate_sphere(radius_sphere_um=2.5,
wl=wavelength_um,
r_core=radius_sphere_um,
r_coat=radius_sphere_um,
n_grid_x=int(size_grid[0]),
n_grid_x=int(shape_grid[0]),
xspan_min=-sizeum[0] / 2 - offset_x_um,
xspan_max=sizeum[0] / 2 - offset_x_um,
n_grid_y=int(size_grid[1]),
n_grid_y=int(shape_grid[1]),
yspan_min=-sizeum[1] / 2 - offset_y_um,
yspan_max=sizeum[1] / 2 - offset_y_um,
n_grid_z=1,
Expand All @@ -202,21 +201,21 @@ def simulate_sphere(radius_sphere_um=2.5,
if arp:
raise
else:
msg = "bhfield: Standard precision failed. " +\
"Trying with arbitrary precision."
msg = "bhfield: Standard precision failed. " \
+ "Retrying with arbitrary precision."
warnings.warn(msg)
arp = True
clear_temp(wdir=wdir)
else:
break

result = load_field(wdir=wdir, size_grid=size_grid)
result = load_field(wdir=wdir, shape_grid=shape_grid)
clear_temp(wdir=wdir)

return result


def run_simulation(wdir=".", arp=True, **kwargs):
def run_simulation(wdir, arp=True, **kwargs):
"""
Example
-------
Expand Down
2 changes: 1 addition & 1 deletion qpsphere/models/mod_mie.py
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ def mie(radius=5e-6, sphere_index=1.339, medium_index=1.333,
"measurement_position_um": propd_um,
"wavelength_nm": wave_nm,
"size_simulation_um": size_um,
"size_grid": grid_size,
"shape_grid": grid_size,
"offset_x_um": offset_um[0],
"offset_y_um": offset_um[1]}

Expand Down
22 changes: 22 additions & 0 deletions tests/test_mie.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
import numpy as np

from qpsphere.models import mie


def test_basic():
data = np.array([0.61718511, 0.63344636, 0.39500855,
0.62834482, 0.5612372, 0.24405929,
0.39808633, 0.24597916, -0.00504312])

qpi = mie(grid_size=(3, 3),
center=(0, 0),
pixel_size=2e-6)
assert np.allclose(data, qpi.pha.flatten())


if __name__ == "__main__":
# Run all tests
loc = locals()
for key in list(loc.keys()):
if key.startswith("test_") and hasattr(loc[key], "__call__"):
loc[key]()
57 changes: 45 additions & 12 deletions tests/test_mie_bhfield.py
Original file line number Diff line number Diff line change
@@ -1,27 +1,60 @@
import warnings

import numpy as np

from qpsphere.models import _bhfield as bh


def test_default_simulation():
data = np.array([0.99915 + 0.0059473j, 1.01950 - 0.0067643j,
1.01950 - 0.0067643j, 0.99915 + 0.0059473j,
1.02340 - 0.0037607j, 0.97630 + 0.4906j,
0.97630 + 0.4906j, 1.02340 - 0.0037607j,
1.02340 - 0.0037607j, 0.97630 + 0.4906j,
0.97630 + 0.4906j, 1.02340 - 0.0037607j,
0.99915 + 0.0059473j, 1.01950 - 0.0067643j,
1.01950 - 0.0067643j, 0.99915 + 0.0059473j])

field = bh.simulate_sphere(shape_grid=(4, 4))
assert np.allclose(field.flatten(), data)


def test_force_arp_warning():
with warnings.catch_warnings(record=True) as rw:
bh.simulate_sphere(radius_sphere_um=7,
refractive_index_sphere=1.5,
size_simulation_um=(20, 20),
shape_grid=(4, 4),
arp=False)
msg = str(rw[0].message).lower()
assert msg.count("bhfield")
assert msg.count("standard precision failed")
assert msg.count("retrying with arbitrary precision")


def test_get_binary():
path1 = bh.fetch.get_binary(arp=False)
path2 = bh.fetch.get_binary(arp=True)
assert path1.exists()
assert path2.exists()


def test_default_simulation():
data = np.array([0.99915+0.0059473j, 1.01950-0.0067643j,
1.01950-0.0067643j, 0.99915+0.0059473j,
1.02340-0.0037607j, 0.97630+0.4906j,
0.97630+0.4906j, 1.02340-0.0037607j,
1.02340-0.0037607j, 0.97630+0.4906j,
0.97630+0.4906j, 1.02340-0.0037607j,
0.99915+0.0059473j, 1.01950-0.0067643j,
1.01950-0.0067643j, 0.99915+0.0059473j])

field = bh.simulate_sphere(size_grid=(4, 4))
assert np.allclose(field.flatten(), data)
def test_known_fail():
try:
bh.simulate_sphere(radius_sphere_um=50,
size_simulation_um=(70, 70),
shape_grid=(4, 4))
except bh.wrap.BHFIELDExecutionError:
pass
else:
assert False, "This simulation should not work with BHFIELD"


def test_shape():
f1 = bh.simulate_sphere(shape_grid=(1, 4))
f2 = bh.simulate_sphere(shape_grid=(4, 1))
assert f1.shape == (1, 4)
assert f2.shape == (4, 1)


if __name__ == "__main__":
Expand Down
22 changes: 22 additions & 0 deletions tests/test_proj.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
import numpy as np

from qpsphere.models import projection


def test_basic():
data = np.array([0.6854384, 0.62821467, 0.41126304,
0.62821467, 0.56522698, 0.30653737,
0.41126304, 0.30653737, 0.])

qpi = projection(grid_size=(3, 3),
center=(0, 0),
pixel_size=2e-6)
assert np.allclose(data, qpi.pha.flatten())


if __name__ == "__main__":
# Run all tests
loc = locals()
for key in list(loc.keys()):
if key.startswith("test_") and hasattr(loc[key], "__call__"):
loc[key]()

0 comments on commit 25efd67

Please sign in to comment.