PySDM_examples.Berry_1967.example_fig_6

  1import numpy as np
  2from matplotlib import pyplot as plt
  3from scipy.optimize import minimize
  4
  5from PySDM.backends import CPU
  6from PySDM.physics import constants as const
  7
  8backend = CPU()
  9um = const.si.um
 10
 11
 12def full_params(params):
 13    return (1, 1, *params)
 14
 15
 16def print_collection_efficiency_portrait(params):
 17    size = 2 * 5.236
 18    plt.figure(num=1, figsize=(size / 2, size))
 19    points = 200
 20    Y_c = np.zeros(points)
 21    x_values = np.linspace(0, 1, points)
 22    pair = np.array([0.0, 0.0])
 23    is_first_in_pair = np.array([True, False])
 24    idx = np.arange(len(is_first_in_pair))
 25    for r in (
 26        radius * const.si.um for radius in (8, 10, 14, 16, 20, 30, 40, 60, 80, 136)
 27    ):
 28        pair[0] = r
 29        for i, _ in enumerate(x_values):
 30            pair[1] = x_values[i] * r
 31            backend._linear_collection_efficiency_body(
 32                params=full_params(params),
 33                output=Y_c[i : i + 1],
 34                radii=pair,
 35                is_first_in_pair=is_first_in_pair,
 36                idx=idx,
 37                length=2,
 38                unit=1 * um,
 39            )
 40        plt.plot(x_values, Y_c, label=f"{r / const.si.um}um")
 41    xticks = np.arange(11) / 10
 42    yticks = np.arange(21) / 10
 43    plt.xticks(xticks, xticks)
 44    plt.yticks(yticks, yticks)
 45    plt.grid()
 46    plt.savefig("berry_fig_6.pdf", format="pdf")
 47
 48
 49p = np.array((0.2, 0.4, 0.6, 0.8, 0.9, 1.0))
 50expected = [
 51    [0, 0, 8.3, 0, 0, 0],  # 8
 52    [0, 36, 55.6, 24, 0, 0],  # 10
 53    [27, 85, 105.3, 96.2, 49.6, 0],  # 14
 54    [46.6, 94.7, 116.5, 107.5, 84.2, 0],  # 16
 55    [70, 108.2, 130.8, 136, 112.8, 0],  # 20
 56    [117, 138, 158, 179, 189, 0],  # 136
 57]
 58expected = np.array(expected)
 59expected /= 100
 60radii = np.array((8 * um, 10 * um, 14 * um, 16 * um, 20 * um, 136 * um))
 61
 62
 63def Y_c_portrait(
 64    params,
 65):
 66    Y_c = np.zeros(expected.shape)
 67    is_first_in_pair = np.array([True, False])
 68    idx = np.arange(len(is_first_in_pair))
 69    pair = np.array([0.0, 0.0])
 70    for i, _ in enumerate(radii):
 71        pair[0] = radii[i]
 72        for j, __ in enumerate(p):
 73            pair[1] = p[j] * radii[i]
 74            backend._linear_collection_efficiency_body(
 75                params=full_params(params),
 76                output=Y_c[i : i + 1, j],
 77                radii=pair,
 78                is_first_in_pair=is_first_in_pair,
 79                idx=idx,
 80                length=2,
 81                unit=1 * um,
 82            )
 83    return Y_c
 84
 85
 86def error(params):
 87    Y_c = Y_c_portrait(params)
 88    return np.sum((Y_c - expected) ** 2) + 0 * np.abs(Y_c[0, 2] - expected[0, 2])
 89
 90
 91def error2(params):
 92    Y_c = Y_c_portrait(params)
 93    return np.sum((Y_c[:-1] - expected[:-1]) ** 2) + 1000 * np.abs(
 94        Y_c[0, 2] - expected[0, 2]
 95    )
 96
 97
 98if __name__ == "__main__":
 99    x0 = np.array([-27, 1.65, -58, 1.9, 1, 1.13, 1, 1, 0.004, 4, 8])
100    x = np.array([-7, 1.78, -20.5, 1.73, 0.26, 1.47, 1, 0.82, -0.003, 4.4, 8])
101    print_collection_efficiency_portrait(x)
102    res = minimize(error, x)
103    print_collection_efficiency_portrait(res.x)
104    res = minimize(error2, res.x)
105    print_collection_efficiency_portrait(res.x)