PySDM_examples.Arabas_et_al_2023.plots

  1import matplotlib
  2import numpy as np
  3from matplotlib import pyplot
  4from PySDM_examples.Arabas_et_al_2023.curved_text import CurvedText
  5from PySDM_examples.Arabas_et_al_2023.frozen_fraction import FrozenFraction
  6
  7from PySDM.physics import si
  8
  9labels = {True: "singular/INAS", False: "time-dependent/ABIFM"}
 10colors = {True: "black", False: "teal"}
 11qi_unit = si.g / si.m**3
 12
 13
 14def make_temperature_plot(data):
 15    pyplot.xlabel("time [s]")
 16
 17    xy1 = pyplot.gca()
 18
 19    xy1.set_ylabel("temperature [K]")
 20    xy1.set_ylim(200, 300)
 21    datum = data[0]["products"]
 22    xy1.plot(datum["t"], datum["T"], marker=".", label="T", color="black")
 23
 24    xy2 = xy1.twinx()
 25    plotted = {singular: False for singular in (True, False)}
 26    for v in data:
 27        datum = v["products"]
 28        xy2.plot(
 29            datum["t"],
 30            np.asarray(datum["qi"]) / qi_unit,  # marker='.',
 31            label=(
 32                f"Monte-Carlo ({labels[v['singular']]})"
 33                if not plotted[v["singular"]]
 34                else ""
 35            ),
 36            color=colors[v["singular"]],
 37        )
 38        plotted[v["singular"]] = True
 39    xy2.set_ylabel("ice water content [g/m3]")
 40
 41    xy1.grid()
 42    xy1.legend()  # bbox_to_anchor=(.2, 1.15))
 43    xy2.legend()  # bbox_to_anchor=(.7, 1.15))
 44
 45
 46def make_freezing_spec_plot(
 47    data,
 48    formulae,
 49    volume,
 50    droplet_volume,
 51    total_particle_number,
 52    surf_spec,
 53    cooling_rate_K_min=None,
 54):
 55    pyplot.xlabel("temperature [K]")
 56    plotted = {singular: False for singular in (True, False)}
 57
 58    prim = pyplot.gca()
 59    for v in data:
 60        datum = v["products"]
 61        color = colors[v["singular"]]
 62        prim.plot(
 63            datum["T"],
 64            np.asarray(datum["qi"]) / qi_unit,
 65            marker=".",
 66            linewidth=0.333,
 67            label=f"{labels[v['singular']]}" if not plotted[v["singular"]] else "",
 68            color=color,
 69        )
 70        plotted[v["singular"]] = True
 71
 72    ff = FrozenFraction(
 73        volume=volume,
 74        droplet_volume=droplet_volume,
 75        total_particle_number=total_particle_number,
 76        rho_w=formulae.constants.rho_w,
 77    )
 78    twin = prim.secondary_yaxis(
 79        "right",
 80        functions=(lambda x: ff.qi2ff(x * qi_unit), lambda x: ff.ff2qi(x) / qi_unit),
 81    )
 82    twin.set_ylabel("frozen fraction")
 83
 84    T = np.linspace(max(datum["T"]), min(datum["T"]))
 85    for multiplier, color in {0.1: "yellow", 1: "brown", 10: "orange"}.items():
 86        qi = (
 87            ff.ff2qi(
 88                formulae.freezing_temperature_spectrum.cdf(
 89                    T, multiplier * surf_spec.median
 90                )
 91            )
 92            / qi_unit
 93        )
 94        prim.plot(
 95            T,
 96            qi,
 97            label="" if multiplier != 1 else "singular CDFs for median surface",
 98            linewidth=2.5,
 99            color=color,
100            linestyle="--",
101        )
102        if multiplier != 1:
103            _ = CurvedText(
104                x=T.squeeze(),
105                y=qi.squeeze(),
106                text=f"                      {multiplier}x median S",
107                va="bottom",
108                color="black",
109                axes=prim,
110            )
111    title = f"$σ_g$=exp({np.log(surf_spec.s_geom):.3g})"
112    if cooling_rate_K_min is not None:
113        title += f", c={cooling_rate_K_min} K/min"
114    prim.set_title(title)
115    # prim.set_ylabel('ice water content [$g/m^3$]')
116    prim.set_yticks([])
117    prim.set_xlim(T[0], T[-1])
118    prim.legend(bbox_to_anchor=(1.02, -0.2))
119    prim.grid()
120
121
122def make_pdf_plot(A_spec, Shima_T_fz, A_range, T_range):
123    N = 256
124    T_space = np.linspace(*T_range, N)
125    A_space = np.linspace(*A_range, N)
126    grid = np.meshgrid(A_space, T_space)
127    sampled_pdf = Shima_T_fz(grid[1], grid[0]) * A_spec.pdf(grid[0])
128
129    fig = pyplot.figure(
130        figsize=(4.5, 6),
131    )
132    ax = fig.add_subplot(111)
133    ax.set_ylabel("freezing temperature [K]")
134    ax.set_yticks(np.linspace(*T_range, num=5, endpoint=True))
135    ax.set_xlabel("insoluble surface [μm$^2$]")
136
137    data = sampled_pdf * si.um**2
138    data[data == 0] = np.nan
139    cnt = ax.contourf(
140        grid[0] / si.um**2,
141        grid[1],
142        data,
143        norm=matplotlib.colors.LogNorm(),
144        cmap="Blues",
145        levels=np.logspace(-3, 0, 7),
146    )
147    cbar = pyplot.colorbar(cnt, ticks=[0.001, 0.01, 0.1, 1.0], orientation="horizontal")
148    cbar.set_label("pdf [K$^{-1}$ μm$^{-2}$]")
149    ax.set_title(f"$σ_g$=exp({np.log(A_spec.s_geom):.3g})")
150
151    ax_histx = ax.inset_axes([0, 1.05, 1, 0.25], sharex=ax)
152    ax_histy = ax.inset_axes([1.05, 0, 0.25, 1], sharey=ax)
153    ax_histx.tick_params(axis="x", labelbottom=False, bottom=False)
154    ax_histx.tick_params(axis="y", labelleft=False, left=False)
155    ax_histy.tick_params(axis="y", labelleft=False, left=False)
156    ax_histy.tick_params(axis="x", labelbottom=False, bottom=False)
157    ax_histx.plot(A_space / si.um**2, np.sum(sampled_pdf, axis=0), color="teal")
158    ax_histy.plot(np.sum(sampled_pdf, axis=1), T_space, color="black")
159
160    pyplot.grid()