PySDM_examples.Shima_et_al_2009.spectrum_plotter
1import matplotlib 2import numpy as np 3from matplotlib import pyplot 4from open_atmos_jupyter_utils import show_plot 5from packaging import version 6from PySDM_examples.Shima_et_al_2009.error_measure import error_measure 7 8from PySDM.physics.constants import si 9 10_matplotlib_version_3_3_3 = version.parse("3.3.0") 11_matplotlib_version_actual = version.parse(matplotlib.__version__) 12 13 14class SpectrumColors: 15 def __init__(self, begining="#2cbdfe", end="#b317b1"): 16 self.b = begining 17 self.e = end 18 19 def __call__(self, value: float): 20 bR, bG, bB = int(self.b[1:3], 16), int(self.b[3:5], 16), int(self.b[5:7], 16) 21 eR, eG, eB = int(self.e[1:3], 16), int(self.e[3:5], 16), int(self.e[5:7], 16) 22 R = bR + int((eR - bR) * value) 23 G = bG + int((eG - bG) * value) 24 B = bB + int((eB - bB) * value) 25 result = f"#{hex(R)[2:4]}{hex(G)[2:4]}{hex(B)[2:4]}" 26 return result 27 28 29class SpectrumPlotter: 30 def __init__(self, settings, title=None, grid=True, legend=True, log_base=10): 31 self.settings = settings 32 self.format = "pdf" 33 self.colors = SpectrumColors() 34 self.smooth = False 35 self.smooth_scope = 2 36 self.legend = legend 37 self.grid = grid 38 self.title = title 39 self.xlabel = "particle radius [µm]" 40 self.ylabel = "dm/dlnr [g/m^3/(unit dr/r)]" 41 self.log_base = log_base 42 self.ax = pyplot 43 self.fig = pyplot 44 self.finished = False 45 46 def finish(self): 47 if self.finished: 48 return 49 self.finished = True 50 if self.grid: 51 self.ax.grid() 52 53 base_arg = { 54 "base" 55 + ( 56 "x" if _matplotlib_version_actual < _matplotlib_version_3_3_3 else "" 57 ): self.log_base 58 } 59 if self.title is not None: 60 try: 61 self.ax.title(self.title) 62 except TypeError: 63 self.ax.set_title(self.title) 64 try: 65 self.ax.xscale("log", **base_arg) 66 self.ax.xlabel(self.xlabel) 67 self.ax.ylabel(self.ylabel) 68 except AttributeError: 69 self.ax.set_xscale("log", **base_arg) 70 self.ax.set_xlabel(self.xlabel) 71 self.ax.set_ylabel(self.ylabel) 72 if self.legend: 73 self.ax.legend() 74 75 def show(self): 76 self.finish() 77 pyplot.tight_layout() 78 show_plot() 79 80 def save(self, file): 81 self.finish() 82 pyplot.savefig(file, format=self.format) 83 84 def plot( 85 self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False 86 ): 87 error = self.plot_analytic_solution(self.settings, t, spectrum, title) 88 if label is not None and add_error_to_label: 89 label += f" error={error:.4g}" 90 self.plot_data(self.settings, t, spectrum, label, color) 91 return error 92 93 def plot_analytic_solution(self, settings, t, spectrum, title): 94 if t == 0: 95 analytic_solution = settings.spectrum.size_distribution 96 else: 97 98 def analytic_solution(x): 99 return settings.norm_factor * settings.kernel.analytic_solution( 100 x=x, t=t, x_0=settings.X0, N_0=settings.n_part 101 ) 102 103 volume_bins_edges = self.settings.formulae.trivia.volume( 104 settings.radius_bins_edges 105 ) 106 dm = np.diff(volume_bins_edges) 107 dr = np.diff(settings.radius_bins_edges) 108 109 pdf_m_x = volume_bins_edges[:-1] + dm / 2 110 pdf_m_y = analytic_solution(pdf_m_x) 111 112 pdf_r_x = settings.radius_bins_edges[:-1] + dr / 2 113 pdf_r_y = pdf_m_y * dm / dr * pdf_r_x 114 115 x = pdf_r_x * si.metres / si.micrometres 116 y_true = ( 117 pdf_r_y 118 * self.settings.formulae.trivia.volume(radius=pdf_r_x) 119 * settings.rho 120 / settings.dv 121 * si.kilograms 122 / si.grams 123 ) 124 125 self.ax.plot(x, y_true, color="black") 126 127 if spectrum is not None: 128 y = spectrum * si.kilograms / si.grams 129 error = error_measure(y, y_true, x) 130 self.title = ( 131 title or f"error measure: {error:.2f}" 132 ) # TODO #327 relative error 133 return error 134 return None 135 136 def plot_data(self, settings, t, spectrum, label, color): 137 if self.smooth: 138 scope = self.smooth_scope 139 if t != 0: 140 new = np.copy(spectrum) 141 for _ in range(2): 142 for i in range(scope, len(spectrum) - scope): 143 new[i] = np.mean(spectrum[i - scope : i + scope + 1]) 144 scope = 1 145 for i in range(scope, len(spectrum) - scope): 146 spectrum[i] = np.mean(new[i - scope : i + scope + 1]) 147 148 x = settings.radius_bins_edges[:-scope] 149 dx = np.diff(x) 150 self.ax.plot( 151 (x[:-1] + dx / 2) * si.metres / si.micrometres, 152 spectrum[:-scope] * si.kilograms / si.grams, 153 label=label or f"t = {t}s", 154 color=color 155 or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)), 156 ) 157 else: 158 self.ax.step( 159 settings.radius_bins_edges[:-1] * si.metres / si.micrometres, 160 spectrum * si.kilograms / si.grams, 161 where="post", 162 label=label or f"t = {t}s", 163 color=color 164 or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)), 165 )
class
SpectrumColors:
15class SpectrumColors: 16 def __init__(self, begining="#2cbdfe", end="#b317b1"): 17 self.b = begining 18 self.e = end 19 20 def __call__(self, value: float): 21 bR, bG, bB = int(self.b[1:3], 16), int(self.b[3:5], 16), int(self.b[5:7], 16) 22 eR, eG, eB = int(self.e[1:3], 16), int(self.e[3:5], 16), int(self.e[5:7], 16) 23 R = bR + int((eR - bR) * value) 24 G = bG + int((eG - bG) * value) 25 B = bB + int((eB - bB) * value) 26 result = f"#{hex(R)[2:4]}{hex(G)[2:4]}{hex(B)[2:4]}" 27 return result
class
SpectrumPlotter:
30class SpectrumPlotter: 31 def __init__(self, settings, title=None, grid=True, legend=True, log_base=10): 32 self.settings = settings 33 self.format = "pdf" 34 self.colors = SpectrumColors() 35 self.smooth = False 36 self.smooth_scope = 2 37 self.legend = legend 38 self.grid = grid 39 self.title = title 40 self.xlabel = "particle radius [µm]" 41 self.ylabel = "dm/dlnr [g/m^3/(unit dr/r)]" 42 self.log_base = log_base 43 self.ax = pyplot 44 self.fig = pyplot 45 self.finished = False 46 47 def finish(self): 48 if self.finished: 49 return 50 self.finished = True 51 if self.grid: 52 self.ax.grid() 53 54 base_arg = { 55 "base" 56 + ( 57 "x" if _matplotlib_version_actual < _matplotlib_version_3_3_3 else "" 58 ): self.log_base 59 } 60 if self.title is not None: 61 try: 62 self.ax.title(self.title) 63 except TypeError: 64 self.ax.set_title(self.title) 65 try: 66 self.ax.xscale("log", **base_arg) 67 self.ax.xlabel(self.xlabel) 68 self.ax.ylabel(self.ylabel) 69 except AttributeError: 70 self.ax.set_xscale("log", **base_arg) 71 self.ax.set_xlabel(self.xlabel) 72 self.ax.set_ylabel(self.ylabel) 73 if self.legend: 74 self.ax.legend() 75 76 def show(self): 77 self.finish() 78 pyplot.tight_layout() 79 show_plot() 80 81 def save(self, file): 82 self.finish() 83 pyplot.savefig(file, format=self.format) 84 85 def plot( 86 self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False 87 ): 88 error = self.plot_analytic_solution(self.settings, t, spectrum, title) 89 if label is not None and add_error_to_label: 90 label += f" error={error:.4g}" 91 self.plot_data(self.settings, t, spectrum, label, color) 92 return error 93 94 def plot_analytic_solution(self, settings, t, spectrum, title): 95 if t == 0: 96 analytic_solution = settings.spectrum.size_distribution 97 else: 98 99 def analytic_solution(x): 100 return settings.norm_factor * settings.kernel.analytic_solution( 101 x=x, t=t, x_0=settings.X0, N_0=settings.n_part 102 ) 103 104 volume_bins_edges = self.settings.formulae.trivia.volume( 105 settings.radius_bins_edges 106 ) 107 dm = np.diff(volume_bins_edges) 108 dr = np.diff(settings.radius_bins_edges) 109 110 pdf_m_x = volume_bins_edges[:-1] + dm / 2 111 pdf_m_y = analytic_solution(pdf_m_x) 112 113 pdf_r_x = settings.radius_bins_edges[:-1] + dr / 2 114 pdf_r_y = pdf_m_y * dm / dr * pdf_r_x 115 116 x = pdf_r_x * si.metres / si.micrometres 117 y_true = ( 118 pdf_r_y 119 * self.settings.formulae.trivia.volume(radius=pdf_r_x) 120 * settings.rho 121 / settings.dv 122 * si.kilograms 123 / si.grams 124 ) 125 126 self.ax.plot(x, y_true, color="black") 127 128 if spectrum is not None: 129 y = spectrum * si.kilograms / si.grams 130 error = error_measure(y, y_true, x) 131 self.title = ( 132 title or f"error measure: {error:.2f}" 133 ) # TODO #327 relative error 134 return error 135 return None 136 137 def plot_data(self, settings, t, spectrum, label, color): 138 if self.smooth: 139 scope = self.smooth_scope 140 if t != 0: 141 new = np.copy(spectrum) 142 for _ in range(2): 143 for i in range(scope, len(spectrum) - scope): 144 new[i] = np.mean(spectrum[i - scope : i + scope + 1]) 145 scope = 1 146 for i in range(scope, len(spectrum) - scope): 147 spectrum[i] = np.mean(new[i - scope : i + scope + 1]) 148 149 x = settings.radius_bins_edges[:-scope] 150 dx = np.diff(x) 151 self.ax.plot( 152 (x[:-1] + dx / 2) * si.metres / si.micrometres, 153 spectrum[:-scope] * si.kilograms / si.grams, 154 label=label or f"t = {t}s", 155 color=color 156 or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)), 157 ) 158 else: 159 self.ax.step( 160 settings.radius_bins_edges[:-1] * si.metres / si.micrometres, 161 spectrum * si.kilograms / si.grams, 162 where="post", 163 label=label or f"t = {t}s", 164 color=color 165 or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)), 166 )
SpectrumPlotter(settings, title=None, grid=True, legend=True, log_base=10)
31 def __init__(self, settings, title=None, grid=True, legend=True, log_base=10): 32 self.settings = settings 33 self.format = "pdf" 34 self.colors = SpectrumColors() 35 self.smooth = False 36 self.smooth_scope = 2 37 self.legend = legend 38 self.grid = grid 39 self.title = title 40 self.xlabel = "particle radius [µm]" 41 self.ylabel = "dm/dlnr [g/m^3/(unit dr/r)]" 42 self.log_base = log_base 43 self.ax = pyplot 44 self.fig = pyplot 45 self.finished = False
def
finish(self):
47 def finish(self): 48 if self.finished: 49 return 50 self.finished = True 51 if self.grid: 52 self.ax.grid() 53 54 base_arg = { 55 "base" 56 + ( 57 "x" if _matplotlib_version_actual < _matplotlib_version_3_3_3 else "" 58 ): self.log_base 59 } 60 if self.title is not None: 61 try: 62 self.ax.title(self.title) 63 except TypeError: 64 self.ax.set_title(self.title) 65 try: 66 self.ax.xscale("log", **base_arg) 67 self.ax.xlabel(self.xlabel) 68 self.ax.ylabel(self.ylabel) 69 except AttributeError: 70 self.ax.set_xscale("log", **base_arg) 71 self.ax.set_xlabel(self.xlabel) 72 self.ax.set_ylabel(self.ylabel) 73 if self.legend: 74 self.ax.legend()
def
plot( self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False):
85 def plot( 86 self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False 87 ): 88 error = self.plot_analytic_solution(self.settings, t, spectrum, title) 89 if label is not None and add_error_to_label: 90 label += f" error={error:.4g}" 91 self.plot_data(self.settings, t, spectrum, label, color) 92 return error
def
plot_analytic_solution(self, settings, t, spectrum, title):
94 def plot_analytic_solution(self, settings, t, spectrum, title): 95 if t == 0: 96 analytic_solution = settings.spectrum.size_distribution 97 else: 98 99 def analytic_solution(x): 100 return settings.norm_factor * settings.kernel.analytic_solution( 101 x=x, t=t, x_0=settings.X0, N_0=settings.n_part 102 ) 103 104 volume_bins_edges = self.settings.formulae.trivia.volume( 105 settings.radius_bins_edges 106 ) 107 dm = np.diff(volume_bins_edges) 108 dr = np.diff(settings.radius_bins_edges) 109 110 pdf_m_x = volume_bins_edges[:-1] + dm / 2 111 pdf_m_y = analytic_solution(pdf_m_x) 112 113 pdf_r_x = settings.radius_bins_edges[:-1] + dr / 2 114 pdf_r_y = pdf_m_y * dm / dr * pdf_r_x 115 116 x = pdf_r_x * si.metres / si.micrometres 117 y_true = ( 118 pdf_r_y 119 * self.settings.formulae.trivia.volume(radius=pdf_r_x) 120 * settings.rho 121 / settings.dv 122 * si.kilograms 123 / si.grams 124 ) 125 126 self.ax.plot(x, y_true, color="black") 127 128 if spectrum is not None: 129 y = spectrum * si.kilograms / si.grams 130 error = error_measure(y, y_true, x) 131 self.title = ( 132 title or f"error measure: {error:.2f}" 133 ) # TODO #327 relative error 134 return error 135 return None
def
plot_data(self, settings, t, spectrum, label, color):
137 def plot_data(self, settings, t, spectrum, label, color): 138 if self.smooth: 139 scope = self.smooth_scope 140 if t != 0: 141 new = np.copy(spectrum) 142 for _ in range(2): 143 for i in range(scope, len(spectrum) - scope): 144 new[i] = np.mean(spectrum[i - scope : i + scope + 1]) 145 scope = 1 146 for i in range(scope, len(spectrum) - scope): 147 spectrum[i] = np.mean(new[i - scope : i + scope + 1]) 148 149 x = settings.radius_bins_edges[:-scope] 150 dx = np.diff(x) 151 self.ax.plot( 152 (x[:-1] + dx / 2) * si.metres / si.micrometres, 153 spectrum[:-scope] * si.kilograms / si.grams, 154 label=label or f"t = {t}s", 155 color=color 156 or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)), 157 ) 158 else: 159 self.ax.step( 160 settings.radius_bins_edges[:-1] * si.metres / si.micrometres, 161 spectrum * si.kilograms / si.grams, 162 where="post", 163 label=label or f"t = {t}s", 164 color=color 165 or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)), 166 )