diff --git a/.gitignore b/.gitignore index 9f410ca..e683b71 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,7 @@ target/ # Python generics .env +__pycache__/ python/__pycache__/ python/kalast/kalast.abi3.so python/kalast/__pycache__/__init__.cpython-311.pyc @@ -17,6 +18,7 @@ python/kalast/__pycache__/__init__.cpython-311.pyc assets/perso/ bin/releases/ runs/ +out/ kalast.log kalast.err -kalast.out \ No newline at end of file +kalast.out diff --git a/cfg/cfg.toml b/cfg/cfg.toml index 2e9beb2..19dfa93 100644 --- a/cfg/cfg.toml +++ b/cfg/cfg.toml @@ -1,23 +1,53 @@ [window] -ambient = [0.07, 0.07, 0.07] -export_frames = true -shadows = true +colormap = { name = "inferno", vmax = 400.0, scalar = "temperature" } +# export_frames = true +# high_dpi = true +# width = 1024 +# height = 768 -[scene] -sun.position.cartesian = [100, 0, 0] +[simulation] +routines = "thermal" +step = 198 +duration = 94780800 -[scene.camera] -name = "test" -position.cartesian = [10, 0, 0] -projection.perspective = 0.5236 -# projection = "orthographic" +[simulation.export] +step = 198 +duration = 172800 +period = 172800 +cooldown_start = 94608000 -[[bodies]] -name = "A" -state.cartesian.position = [0, 0, 0] -material.albedo = 0.1 +# 1728000 = 20 days +# 1900800 = 22 days +# 31536000 = 365 days +# 31708800 = 367 days +# 94608000 = 365*3 days +# 94780800 = 365*3+2 days + +[scene] +sun.position.cartesian = [149597870.7, 0, 0] +camera.position.cartesian = [5, 0, 0] [[bodies]] -name = "B" -state.cartesian.position = [3, 0, 0] -material.albedo = 0.1 +name = "Simple" +color = "data" +interior.grid1d.increasing = { skin = "one", m = 10, n = 5, b = 20 } +record = { columns = [0], mesh = true, depth = true } +spin = { period = 86400 } +temperature.scalar = 1 + +[bodies.mesh] +shape.shape = "triangle" +position = [0.0, -0.3333, 0.3333] + +# 1) Identity +# 2 to 3) Rotation 90° around X and Y axis +# orientation = [1, 0, 0, 0, 1, 0, 0, 0, 1] +# orientation = [1, 0, 0, 0, 0, 1, 0, -1, 0] +orientation = [0, 0, -1, 0, 1, 0, 1, 0, 0] + +[bodies.material] +albedo = 0.1 +emissivity = 0.9 +thermal_inertia = 500.0 +density = 2000.0 +heat_capacity = 600.0 diff --git a/examples/apophis-gravity/cfg/cfg.toml b/examples/apophis-gravity/cfg/cfg.toml new file mode 100644 index 0000000..3252c85 --- /dev/null +++ b/examples/apophis-gravity/cfg/cfg.toml @@ -0,0 +1,27 @@ +[window] +# colormap = { name = "cividis", vmax = 100, scalar = "file" } +colormap = { name = "cividis", vmax = 50, scalar = "file" } +export_frames = true +high_dpi = true +# width = 1024 +# height = 768 + +[simulation] +step = 60 +duration = 12000 + +[simulation.export] +step = 60 +duration = 12000 +period = 12000 +cooldown_start = 0 + +[scene] +camera.position.cartesian = [0, 1, 0] + +[[bodies]] +name = "Apophis" +mesh = { shape.path = "/Users/gregoireh/data/kalast/runs/Apophis flyby Hakan/input/apophis_new.obj", factor = [1e-3, 1e-3, 1e-3] } +color = "data" +# file_data = { path = "/Users/gregoireh/data/kalast/runs/Gravity Elisa/input/results_slope_change_percentage-2.csv" } +file_data = { path = "/Users/gregoireh/data/kalast/runs/Gravity Elisa/input/results_slope_change-2.csv" } \ No newline at end of file diff --git a/examples/apophis/cfg/cfg.toml b/examples/apophis/cfg/cfg.toml new file mode 100644 index 0000000..01e0460 --- /dev/null +++ b/examples/apophis/cfg/cfg.toml @@ -0,0 +1,47 @@ +[window] +colormap = { name = "inferno", vmax = 400.0, scalar = "temperature" } +export_frames = true +high_dpi = true +# width = 1024 +# height = 768 + +[simulation] +routines = "thermal" +step = 600 +duration = 25920000 +file = { path = "/Users/gregoireh/Downloads/hakan/forkalast-2.csv", behavior = "loop" } + +[simulation.export] +step = 60 +duration = 109800 +period = 1296000 +cooldown_start = 13399200 + +# 1296000 = 360h +# 12960000 = 3600h +# 25920000 = 3600h * 2 +# 109800 = 30.5h +# 439200 = 30.5h * 4 +# 549000 = 30.5h * 5 +# 1098000 = 30.5h * 10 +# 13399200 = 30.5h * 4 + 3600h * 1 +# 117079200 = 30.5h * 4 + 3600h * 9 + +[scene] +sun.position = "file" +camera.position.cartesian = [-1, 0, 0] + +[[bodies]] +name = "Apophis" +color = "data" +mesh = { shape.path = "/Users/gregoireh/Downloads/hakan/apophis_new.obj", factor = [1e-3, 1e-3, 1e-3] } +interior.grid1d.linear = { size = 40, a = 2e-2 } +state = "file" +record = { rows = [0], mesh = true, depth = true } + +[bodies.material] +albedo = 0.3 +emissivity = 0.9 +thermal_inertia = 500.0 +density = 2100.0 +heat_capacity = 600.0 diff --git a/examples/juventas/cfg/cfg.toml b/examples/juventas/cfg/cfg.toml index 6141622..dac1fe0 100644 --- a/examples/juventas/cfg/cfg.toml +++ b/examples/juventas/cfg/cfg.toml @@ -1,10 +1,8 @@ [window] ambient = [0.07, 0.07, 0.07] export_frames = true -# width = 2048 -# height = 1536 -width = 1024 -height = 768 +width = 2048 +height = 1536 shadows = true [spice] @@ -13,12 +11,20 @@ frame = "JUVENTAS_NAVCAM" origin = "DIDYMOS" [simulation] -start = "2027 MAR 23 16:02:00.0 TDB" +# start = "2027 MAR 23 16:02:00.0 TDB" +# start = "2027 MAR 26 16:16:10.0 TDB" +# start = "2027 MAR 30 16:16:10.0 TDB" +# start = "2027 APR 02 16:16:10.0" +# start = "2027 APR 05 16:16:10.0" +# start = "2027 MAY 03 16:16:10.0" +start = "2027 MAY 06 06:08:47.0" + + +duration = 2678400 step = 60 -duration = 259080 -# 86400 * 60 = 5184000 -# 86400 * 3 = 259200 -# ^-120 = 259080 +export = { step = 60, duration = 2678400 } + +# 86400 * 31 = 2678400 [scene] sun.position = "spice" @@ -29,8 +35,7 @@ name = "JUVENTAS" position = "spice" up = [1, 0, 0] direction.vector = [0, 0, 1] -# projection.perspective = 0.4974188 # 28.5° -projection.perspective = 0.09 +projection.perspective = 0.4974188 # 28.5° [[bodies]] name = "DIDYMOS" diff --git a/plot/__init__.py b/plot/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/plot/apophis.py b/plot/apophis.py new file mode 100644 index 0000000..555aa96 --- /dev/null +++ b/plot/apophis.py @@ -0,0 +1,72 @@ +from pathlib import Path + +import numpy +from natsort import natsorted +from pyarrow import csv + +import plot + +d = dict() + +path_data = Path("/Users/gregoireh/data/kalast") +path_runs = path_data / "runs" +path_draft = path_data / "graphs/draft" + +path_run = path_runs / "Apophis flyby Hakan" + +path_cfg = path_run / "cfg/cfg.toml" +path_simu = path_run / "simu" +path_rec = path_simu / "rec" + +d["nf"] = 3996 +d["body"] = "Apophis" +path_setup_body = path_simu / f"{d['body']}" +path_progress = path_simu / "progress.csv" + +path_sph = path_setup_body / "mesh.csv" +read_options = csv.ReadOptions( + column_names=[ + "x", + "y", + "z", + "lon", + "lat", + "rad", + ], + skip_rows=1, +) +tab = csv.read_csv(path_sph, read_options) +df = tab.to_pandas() +d["centers"] = numpy.array( + [df["x"].to_numpy(), df["y"].to_numpy(), df["z"].to_numpy()] +).T +d["sph"] = numpy.array( + [df["lon"].to_numpy(), df["lat"].to_numpy(), df["rad"].to_numpy()] +).T + +path_depth = path_setup_body / "depth.csv" +read_options = csv.ReadOptions( + column_names=[ + "depth", + ], + skip_rows=1, +) +tab = csv.read_csv(path_depth, read_options) +df = tab.to_pandas() +d["depth"] = df["depth"].to_numpy() + +list_path_date = [p for p in natsorted((path_rec).glob("*"), key=str) if p.is_dir()] +list_elapsed = [int(p.name) for p in list_path_date] +it_list_path_date = iter(list_path_date) + +path_date = next(it_list_path_date) + +p_csv = path_date / d["body"] / "temperatures/temperatures-all.csv" +read_options = csv.ReadOptions() +tab = csv.read_csv(p_csv, read_options) +df = tab.to_pandas() +d["tmp"] = df["tmp"].array.reshape((d["nf"], -1)).T +print(d["tmp"].shape) +print(d["tmp"][0].min(), d["tmp"][0].max(), d["tmp"][0].mean()) + +plot.surface.plot(d, save=True) \ No newline at end of file diff --git a/plot/main.mplstyle b/plot/main.mplstyle new file mode 100644 index 0000000..45a5752 --- /dev/null +++ b/plot/main.mplstyle @@ -0,0 +1,26 @@ +figure.figsize: 4.0, 4.0 +figure.autolayout: True + +svg.fonttype: none +font.family: Helvetica, Arial, sans-serif + +xtick.top: True +xtick.bottom: True +xtick.major.size: 4.0 +xtick.major.width: 1.0 +xtick.minor.size: 2.0 +xtick.minor.width: 1.0 +xtick.direction: in + +ytick.left: True +ytick.right: True +ytick.major.size: 4.0 +ytick.major.width: 1.0 +ytick.minor.size: 2.0 +ytick.minor.width: 1.0 +ytick.direction: in + +# grid.color: "#f8f8f8ff" +# grid.linewidth: 0.25 + +savefig.format: svg \ No newline at end of file diff --git a/plot/plot/__init__.py b/plot/plot/__init__.py new file mode 100644 index 0000000..f536672 --- /dev/null +++ b/plot/plot/__init__.py @@ -0,0 +1 @@ +from plot import daily, depth, smap diff --git a/plot/plot/daily.py b/plot/plot/daily.py new file mode 100644 index 0000000..5d4c228 --- /dev/null +++ b/plot/plot/daily.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python + +from pathlib import Path + +import matplotlib +import numpy +import scipy +import util +from matplotlib import pyplot, ticker +from mpl_toolkits.basemap import Basemap + +Data = dict[str, numpy.array] +Data_m = dict[str, str] + + +def plot(d: dict[str, Data], save: bool = False): + pyplot.style.use("plot/main.mplstyle") + fig, ax = pyplot.subplots(figsize=(15, 7.3)) + + ax.set_xlabel("Time elapsed (h)") + ax.set_ylabel("Temperature (K)") + + dt = 200 + xmax = d["tmp-cols"].shape[0] * dt + x = numpy.arange(0, xmax, dt) / 3600 + y = d["tmp-cols"][:, 0] + + ax.plot(x, y) + + ax.grid(True) + + ax.set_xlim(0, x.max()) + ax.set_ylim(200, 400) + + # loc = ticker.MultipleLocator(base=30.0) + # ax.xaxis.set_major_locator(loc) + + loc = ticker.MultipleLocator(base=50) + ax.yaxis.set_major_locator(loc) + + if save: + out = Path("out/surface") + if not out.exists(): + out.mkdir(parents=True) + + fig.savefig(out / "daily.png", bbox_inches="tight", dpi=300) + + pyplot.show() \ No newline at end of file diff --git a/plot/plot/depth.py b/plot/plot/depth.py new file mode 100644 index 0000000..1b31c9a --- /dev/null +++ b/plot/plot/depth.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python + +from pathlib import Path + +import matplotlib +import numpy +import scipy +import util +from matplotlib import pyplot, ticker +from mpl_toolkits.basemap import Basemap + +Data = dict[str, numpy.array] +Data_m = dict[str, str] + + +def plot(d: dict[str, Data], save: bool = False): + pyplot.style.use("plot/main.mplstyle") + fig, ax = pyplot.subplots(figsize=(15, 7.3)) + + ax.set_xlabel("Temperature (K)") + ax.set_ylabel("Depth (m)") + + x1 = d["tmp-cols"].min(axis=0) + x2 = d["tmp-cols"].mean(axis=0) + x3 = d["tmp-cols"].max(axis=0) + y = d["depth"] + + ax.plot(x1, y) + ax.plot(x2, y) + ax.plot(x3, y) + + ax.grid(True) + + ax.set_xlim(200, 400) + ax.set_ylim(1.5, y[1]) + + ax.set(yscale="log") + + # loc = ticker.MultipleLocator(base=30.0) + # ax.xaxis.set_major_locator(loc) + + loc = ticker.MultipleLocator(base=50) + ax.xaxis.set_major_locator(loc) + + if save: + out = Path("out/surface") + if not out.exists(): + out.mkdir(parents=True) + + fig.savefig(out / "depth.png", bbox_inches="tight", dpi=300) + + pyplot.show() \ No newline at end of file diff --git a/plot/plot/smap.py b/plot/plot/smap.py new file mode 100644 index 0000000..03be470 --- /dev/null +++ b/plot/plot/smap.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python + +from pathlib import Path + +import matplotlib +import numpy +import scipy +import util +from matplotlib import pyplot, ticker +from mpl_toolkits.basemap import Basemap + +Data = dict[str, numpy.array] +Data_m = dict[str, str] + +cbar_bottom = True +interp = True +proj = "cyl" + + +def plot(d: dict[str, Data], save: bool = False): + pyplot.style.use("plot/main.mplstyle") + fig, ax = pyplot.subplots(figsize=(15, 7.3)) + + if proj == "cyl": + map = Basemap( + ax=ax, + projection="cyl", + llcrnrlat=-90, + urcrnrlat=90, + llcrnrlon=-180, + urcrnrlon=180, + resolution="c", + ) + elif proj == "moll": + map = Basemap( + ax=ax, + projection="moll", + lon_0=0, + resolution="c", + ) + elif proj == "ortho": + map = Basemap(ax=ax, projection="ortho", lon_0=40, lat_0=20, resolution="l") + elif proj == "nsper": + map = Basemap( + ax=ax, + projection="nsper", + lon_0=40, + lat_0=40, + satellite_height=3000 * 1000.0, + resolution="l", + ) + + if proj is not None and proj != "cyl": + map.drawparallels(numpy.arange(-90.0, 120.0, 30.0)) + map.drawmeridians(numpy.arange(0.0, 420.0, 60.0)) + + if proj is None or proj == "cyl": + # ax.set_xlabel("Longitude (°)") + # ax.set_ylabel("Latitude (°)") + + ax.set_xlim(-180, 180) + ax.set_ylim(-90, 90) + + cmin = 0 + cmax = 400 + vmin = d["tmp"].min() + vmax = d["tmp"].max() + + norm = matplotlib.colors.Normalize(vmin=cmin, vmax=cmax) + cmap = matplotlib.cm.inferno + + lvl = 1 + lvls = numpy.arange(cmin, cmax + lvl / 10, lvl, dtype=int) + + p = d["sph"] * util.DPR + xy = p[:, :2] + z = d["tmp"][0] + + # n = z.size + nx = 50 + ny = nx + + xgrid = numpy.linspace(-numpy.pi, numpy.pi, nx) * util.DPR + ygrid = numpy.linspace(-numpy.pi / 2, numpy.pi / 2, ny) * util.DPR + grid = numpy.array(numpy.meshgrid(xgrid, ygrid)) + + fgrid = grid.reshape(2, -1).T + zflat = scipy.interpolate.RBFInterpolator(xy, z)(fgrid) + zgrid = zflat.reshape(nx, ny) + + # tmap = map.contourf(*grid, zgrid, levels=lvls, cmap=cmap, norm=norm, latlon=True) + + # pcolormesh + # contourf + tmap = map.contourf( + *grid, + zgrid, + cmap=cmap, + norm=norm, + latlon=True, + levels=lvls, + # shading="gouraud", + ) + + # p = ax.scatter(*xy.T, c=z, s=50, ec="k", vmin=0, vmax=400) + + if cbar_bottom: + orientation = "horizontal" + pad = 0.05 + shrink = 0.4 + aspect = 30 + else: + orientation = "vertical" + pad = 0.02 + shrink = 0.855 + aspect = 23 + + _ = fig.colorbar( + tmap, + label="Temperature (K)", + orientation=orientation, + shrink=shrink, + pad=pad, + aspect=aspect, + ) + + cax = fig.axes[-1] + + loc = ticker.MultipleLocator(base=50.0) + if not cbar_bottom: + cax.plot([0, 0.15], [vmin, vmin], color="k", linewidth=0.9) + cax.plot([0, 0.15], [vmax, vmax], color="k", linewidth=0.9) + cax.text(-0.08, vmin, f"{vmin:.0f}", ha="right", va="center_baseline") + cax.text(-0.08, vmax, f"{vmax:.0f}", ha="right", va="center_baseline") + cax.yaxis.set_major_locator(loc) + else: + cax.plot([vmin, vmin], [0.75, 1.0], color="k", linewidth=0.9) + cax.plot([vmax, vmax], [0.75, 1.0], color="k", linewidth=0.9) + cax.text(vmin, 1.1, f"{vmin:.0f}", ha="center", va="bottom") + cax.text(vmax, 1.1, f"{vmax:.0f}", ha="center", va="bottom") + cax.xaxis.set_major_locator(loc) + + if proj is None or proj == "cyl": + loc = ticker.MultipleLocator(base=30.0) + ax.xaxis.set_major_locator(loc) + + loc = ticker.MultipleLocator(base=30) + ax.yaxis.set_major_locator(loc) + + if save: + out = Path("out/surface") + if not out.exists(): + out.mkdir(parents=True) + + fig.savefig(out / "temperature.png", bbox_inches="tight", dpi=300) + + pyplot.show() \ No newline at end of file diff --git a/plot/simple.py b/plot/simple.py new file mode 100644 index 0000000..46c6212 --- /dev/null +++ b/plot/simple.py @@ -0,0 +1,72 @@ +from pathlib import Path + +import numpy +from natsort import natsorted +from pyarrow import csv + +import plot + +d = dict() + +path_data = Path("/Users/gregoireh/data/kalast") +path_desktop = Path("/Users/gregoireh/Desktop/kalast-runs") +path_runs = path_data / "runs" +path_draft = path_data / "graphs/draft" + +# path_run = path_runs / "Apophis flyby Hakan" +path_run = path_desktop / "run #29" + +path_cfg = path_run / "cfg/cfg.toml" +path_simu = path_run / "simu" +path_rec = path_simu / "rec" + +d["nf"] = 1 +d["body"] = "Simple" +path_setup_body = path_simu / f"{d['body']}" +path_progress = path_simu / "progress.csv" + +path_sph = path_setup_body / "mesh.csv" +read_options = csv.ReadOptions( + column_names=["x", "y", "z", "lon", "lat", "rad"], skip_rows=1 +) +tab = csv.read_csv(path_sph, read_options) +df = tab.to_pandas() +d["centers"] = numpy.array( + [df["x"].to_numpy(), df["y"].to_numpy(), df["z"].to_numpy()] +).T +d["sph"] = numpy.array( + [df["lon"].to_numpy(), df["lat"].to_numpy(), df["rad"].to_numpy()] +).T + +path_depth = path_setup_body / "depth.csv" +read_options = csv.ReadOptions(column_names=["depth"], skip_rows=1) +tab = csv.read_csv(path_depth, read_options) +df = tab.to_pandas() +d["depth"] = df["depth"].to_numpy() +d["nz"] = d["depth"].size + +list_path_date = [p for p in natsorted((path_rec).glob("*"), key=str) if p.is_dir()] +list_elapsed = [int(p.name) for p in list_path_date] +it_list_path_date = iter(list_path_date) + +path_date = next(it_list_path_date) + +p_csv = path_date / d["body"] / "temperatures/temperatures-all.csv" +read_options = csv.ReadOptions() +tab = csv.read_csv(p_csv, read_options) +df = tab.to_pandas() +d["tmp-all"] = df["tmp"].array.reshape((d["nf"], -1)).T +print(d["tmp-all"].shape) +print(d["tmp-all"][0].min(), d["tmp-all"][0].max(), d["tmp-all"][0].mean()) + +p_csv = path_date / d["body"] / "temperatures/temperatures-columns.csv" +read_options = csv.ReadOptions() +tab = csv.read_csv(p_csv, read_options) +df = tab.to_pandas() +# d["tmp-cols"] = df["0"].array.reshape((d["nz"], -1)).T +d["tmp-cols"] = df["0"].array.reshape((-1, d["nz"])) +print(d["tmp-cols"].shape) + +# plot.smap.plot(d, save=True) +plot.daily.plot(d, save=True) +plot.depth.plot(d, save=True) \ No newline at end of file diff --git a/plot/util.py b/plot/util.py new file mode 100644 index 0000000..c5d04e8 --- /dev/null +++ b/plot/util.py @@ -0,0 +1,96 @@ +import glm +import numpy + +RPD = numpy.pi / 180 +DPR = 1 / RPD + +axes = ["x", "y", "z"] +XYZ = numpy.eye(3) + +AU_km = 1.495978707e8 +solar_constant = 1367 +stephan_boltzmann = 5.670374419e-8 + +fmt_time_spice = "YYYY-MM-DD HR:MN:SC.### ::RND ::UTC" + +# LPO: Mars ~1% fovcov ~9h30mn before flyby (~7800 pixels) +# Mars ~2.5% fovcov ~6h0mn before flyby +# Mars ~5% fovcov ~4h0mn before flyby +# Mars ~10% fovcov ~3h0mn before flyby +# Mars ~25% fovcov ~2h0mn before flyby +# Mars ~50% fovcov ~1h0mn before flyby +dates = dict( + LPC_LAUNCH="2024-10-27 14:15:18.680", + LPC_MGA="2025-03-16 19:14:50.8", + LPO_LAUNCH="2024-10-07 16:27:53.417", + LPO_MGA="2025-03-12 23:35:50.8", + ARRIVAL="2026 DEC 02 01:59:59.999", + HERA_START="2027 FEB 09 16:07:00.0 TDB", + HERA_END="2027 JUL 24 20:06:00.0 TDB", + JUVENTAS_START="2027 MAR 24 12:00:00.0 TDB", + JUVENTAS_END="2027 MAR 26 16:01:00.0 TDB", +) + +scenario = "LPC" + +tiri_id = -91200 +tiri_pxx = 1024 +tiri_pxy = 768 +tiri_fovx = 13.0 * RPD +tiri_fovy = 10.0 * RPD + +afc_pxx = 1024 +afc_pxy = 1024 +afc_fovx = 5.47 * RPD +afc_fovy = 5.47 * RPD + +earth_orbit_period = 365.25 * 86400 +earth_diameter = 12742 + +moon_orbit_period = 29.5 * 86400 +moon_diameter = 3474.8 + +mars_id = 499 +mars_orbit_period = 687 * 86400 +mars_diameter = 6779 +mars_radii = numpy.array([6378.1366, 6378.1366, 6356.7519]) + +phobos_orbit_period = 7 * 3600 + 39 * 60 +phobos_diameter = 22.533 + +deimos_orbit_period = 30.312 * 3600 +deimos_diameter = 12.4 + + +# u1 & u2 are unit vectors +def angle_between(u1, u2): + return numpy.arccos(numpy.clip(u1 @ u2, -1.0, 1.0)) + + +def mat_angle_axis(angle: float, axis: numpy.array): + m = glm.mat3(glm.rotate(angle, glm.vec3(axis))) + return numpy.array(m.to_list()) + + +def sign(x): + return "-" if x < 0 else "" + + +def hava(x: float, y: float) -> tuple[str, str]: + if x < 0: + ha = "right" + if x > 0: + ha = "left" + if y < 0: + va = "top" + if y > 0: + va = "bottom" + return ha, va + + +def sha(ha: str) -> float: + return 1 if ha == "left" else -1 + + +def sva(ha: str) -> float: + return 1 if ha == "bottom" else -1 \ No newline at end of file diff --git a/preferences.toml b/preferences.toml index 82e7a89..cd5b21f 100644 --- a/preferences.toml +++ b/preferences.toml @@ -4,4 +4,5 @@ touchpad_controls = true [preferences.debug] config = true general = true - +window = false +thermal_stats = false diff --git a/src/body/material.rs b/src/body/material.rs index 5ac54d8..9c4f316 100644 --- a/src/body/material.rs +++ b/src/body/material.rs @@ -2,17 +2,20 @@ use crate::util::*; use serde::{Deserialize, Serialize}; -#[allow(unused)] -fn thermal_skin_depth(diffusivity: Float, period: Float) -> Float { +pub fn thermal_skin_depth_one(diffusivity: Float, period: Float) -> Float { (diffusivity * period / PI).sqrt() } +pub fn thermal_skin_depth_two_pi(diffusivity: Float, period: Float) -> Float { + (4.0 * PI * diffusivity * period).sqrt() +} + /// units: W.m^{-1}.K^{-1} -fn conductivity(thermal_inertia: Float, density: Float, heat_capacity: Float) -> Float { +pub fn conductivity(thermal_inertia: Float, density: Float, heat_capacity: Float) -> Float { thermal_inertia.powi(2) / (density * heat_capacity) } -fn diffusivity(conductivity: Float, density: Float, heat_capacity: Float) -> Float { +pub fn diffusivity(conductivity: Float, density: Float, heat_capacity: Float) -> Float { conductivity / (density * heat_capacity) } @@ -84,7 +87,11 @@ impl Material { diffusivity(self.conductivity(), self.density, self.heat_capacity) } - pub fn thermal_skin_depth(&self, period: Float) -> Float { - thermal_skin_depth(self.diffusivity(), period) + pub fn thermal_skin_depth_one(&self, period: Float) -> Float { + thermal_skin_depth_one(self.diffusivity(), period) + } + + pub fn thermal_skin_depth_two_pi(&self, period: Float) -> Float { + thermal_skin_depth_two_pi(self.diffusivity(), period) } } diff --git a/src/config/body.rs b/src/config/body.rs index e6ae520..9146c31 100644 --- a/src/config/body.rs +++ b/src/config/body.rs @@ -196,6 +196,9 @@ pub struct Body { #[serde(default)] pub record: Record, + #[serde(default)] + pub file_data: Option, + #[serde(flatten)] extra: HashMap, } @@ -264,11 +267,21 @@ pub struct Mesh { #[serde(default)] pub shape: MeshSource, - /// Resize factor to be applied to the mesh. + /// Resize correction to be applied to the mesh. /// Default is `[1, 1, 1]`. #[serde(default = "default_mesh_factor")] pub factor: Vec3, + /// Position correction to be applied to the mesh. + /// Default is `[0, 0, 0]`. + #[serde(default = "default_vec3_zeros")] + pub position: Vec3, + + /// Orientation correction to be applied to the mesh. + /// Default Mat3 identity. + #[serde(default = "default_mat3_identity")] + pub orientation: Mat3, + /// Wether to render vertex- (smooth) or facet-wise (flat). /// Default is flat, smooth is `false`. /// Smooth does not work for thermophysical model, only for viewer. @@ -281,15 +294,13 @@ impl Default for Mesh { Self { shape: MeshSource::default(), factor: default_mesh_factor(), + position: default_vec3_zeros(), + orientation: default_mat3_identity(), smooth: false, } } } -fn default_mesh_factor() -> Vec3 { - vec3(1.0, 1.0, 1.0) -} - #[derive(Clone, Debug, Serialize, Deserialize)] pub enum MeshSource { #[serde(rename = "shape")] @@ -420,11 +431,28 @@ pub enum InteriorGrid1D { #[serde(rename = "exp")] Exp { size: usize, a: Float }, + #[serde(rename = "increasing")] + Increasing { + skin: SkinDepth, + m: usize, + n: usize, + b: usize, + }, + /// Unimplemented. #[serde(rename = "file")] File { path: PathBuf }, } +#[derive(Clone, Debug, Serialize, Deserialize)] +pub enum SkinDepth { + #[serde(rename = "one")] + One, + + #[serde(rename = "two_pi")] + TwoPi, +} + /** # Configuration of the Spin of the Body. @@ -742,9 +770,6 @@ pub struct FileSetup { #[serde(default)] pub path: Option, - #[serde(default)] - pub row_multiplier: Option, - #[serde(default)] pub behavior: Option, } @@ -814,7 +839,7 @@ impl FileColumnsOut { } } -#[derive(Clone, Debug, Serialize, Deserialize)] +#[derive(Clone, Debug, Serialize, Deserialize, PartialEq)] pub enum FileBehavior { #[serde(rename = "stop")] Stop, @@ -922,4 +947,22 @@ pub struct Record { #[serde(default)] pub cells: Vec, + + #[serde(default)] + pub mesh: bool, + + #[serde(default)] + pub depth: bool, +} + +fn default_mesh_factor() -> Vec3 { + vec3(1.0, 1.0, 1.0) +} + +fn default_vec3_zeros() -> Vec3 { + Vec3::zeros() +} + +fn default_mat3_identity() -> Mat3 { + Mat3::identity() } diff --git a/src/config/preferences.rs b/src/config/preferences.rs index e70f494..152355f 100644 --- a/src/config/preferences.rs +++ b/src/config/preferences.rs @@ -20,6 +20,7 @@ pub struct Debug { pub general: Option, pub window: Option, pub simulation: Option, + pub thermal_stats: Option, } #[derive(Clone, Debug, Serialize, Deserialize, Default)] diff --git a/src/config/window.rs b/src/config/window.rs index b0610a7..a180ad3 100644 --- a/src/config/window.rs +++ b/src/config/window.rs @@ -146,6 +146,9 @@ pub enum CfgScalar { #[serde(rename = "temperature")] Temperature, + + #[serde(rename = "file")] + File, } impl Default for CfgScalar { diff --git a/src/simu/converge.rs b/src/simu/converge.rs index 600078e..c02c2a4 100644 --- a/src/simu/converge.rs +++ b/src/simu/converge.rs @@ -48,7 +48,7 @@ pub fn check>( let tsd = asteroid.surface.faces[c0_ii] .vertex .material - .thermal_skin_depth(cb.spin.period); + .thermal_skin_depth_one(cb.spin.period); let ii_tsd_next = zdepth.iter().position(|&x| x > tsd).unwrap(); let c0 = df.column(&c0_name).unwrap().f64().unwrap(); diff --git a/src/simu/export.rs b/src/simu/export.rs index 4f5cc39..7cc83a2 100644 --- a/src/simu/export.rs +++ b/src/simu/export.rs @@ -67,11 +67,17 @@ impl Export { self.exporting = true; self.exporting_started_elapsed = elapsed as _; self.remaining_duration_export = cfg.simulation.export.duration as _; + println!("Detected export time."); + println!("Simulation time step: {}", time.time_step); time.set_time_step(cfg.simulation.export.step); + println!("Export time step: {}", time.time_step); } else if self.cooldown_export - (dt as i64) < 0 { // So export does not really start here, but the time step is adapted to not miss the beginning of export // (in case export time step is smaller than simulation time step). + println!("Detected pre-export time."); + println!("Simulation time step: {}", time.time_step); time.set_time_step(cfg.simulation.export.step); + println!("Export time step: {}", time.time_step); } } @@ -121,7 +127,10 @@ impl Export { self.is_first_it_export = true; self.cooldown_export = (cfg.simulation.export.period - cfg.simulation.export.duration) as _; + println!("End of export."); + println!("Export time step: {}", time.time_step); time.set_time_step(cfg.simulation.step); + println!("Simulation time step: {}", time.time_step); // let _cvg = kalast::simu::converge::check_all(&mut bodies, &folder_tpm, &cfg.time.export); } @@ -150,13 +159,8 @@ impl Export { self.exporting_started_elapsed as _, &config.bodies[body].name, ); - let folder_tpm = folders.simu_rec_time_body_temperatures( - self.exporting_started_elapsed as _, - &config.bodies[body].name, - ); let folder_img = folders.simu_rec_time_frames(self.exporting_started_elapsed as _); fs::create_dir_all(&folder_state).unwrap(); - fs::create_dir_all(&folder_tpm).unwrap(); fs::create_dir_all(&folder_img).unwrap(); let mut df = df!( diff --git a/src/simu/folders.rs b/src/simu/folders.rs index ee29078..dc630f5 100644 --- a/src/simu/folders.rs +++ b/src/simu/folders.rs @@ -87,7 +87,7 @@ impl FoldersRun { return; } - fs::create_dir_all(&self.path).unwrap(); - fs::copy(path, &self.path).unwrap(); + // fs::create_dir_all(&self.path).unwrap(); + // fs::copy(path, &self.path).unwrap(); } } diff --git a/src/simu/routines/core.rs b/src/simu/routines/core.rs index 5b74214..c671561 100644 --- a/src/simu/routines/core.rs +++ b/src/simu/routines/core.rs @@ -1,12 +1,16 @@ +use std::{fs, path::Path}; + use crate::{ compute_cosine_emission_angle, compute_cosine_incidence_angle, compute_cosine_phase_angle, - config::Body, config::CfgCamera, config::CfgCameraDirection, config::CfgCameraPosition, - config::CfgScalar, config::CfgSun, config::CfgSunPosition, config::Config, - config::FileBehavior, config::FileColumns, config::FileColumnsOut, config::FileSetup, - config::FrameCenter, config::SpicePosition, config::SpiceState, config::State, - config::StateCartesian, config::DEFAULT_ABCORR, config::DEFAULT_FRAME, find_ref_orbit, - matrix_orientation_obliquity, matrix_spin, position_in_inertial_frame, update_colormap_scalar, - util::*, AirlessBody, BodyData, FoldersRun, MovementMode, Time, Window, WindowScene, + config::{ + Body, CfgCamera, CfgCameraDirection, CfgCameraPosition, CfgScalar, CfgSun, CfgSunPosition, + Config, FileBehavior, FileColumns, FileColumnsOut, FileSetup, FrameCenter, SpicePosition, + SpiceState, State, StateCartesian, DEFAULT_ABCORR, DEFAULT_FRAME, + }, + find_ref_orbit, matrix_orientation_obliquity, matrix_spin, position_in_inertial_frame, + update_colormap_scalar, + util::*, + AirlessBody, BodyData, FoldersRun, Interior, MovementMode, Time, Window, WindowScene, }; use downcast_rs::{impl_downcast, DowncastSync}; @@ -25,6 +29,104 @@ pub trait Routines: DowncastSync { ) { } + fn fn_update_file_index_core( + &self, + _config: &Config, + time: &mut Time, + file: &Path, + opt_loop: bool, + ) { + let df = match CsvReader::from_path(&file) { + Ok(reader) => reader.has_header(true).finish().unwrap(), + Err(e) => panic!("{}", e), + }; + + let elapsed = df.column("elapsed").unwrap(); + let mut time_elapsed = time.elapsed_seconds(); + let mut found = false; + let mut current_time_file = None; + let mut next_tried_time_file = None; + + if opt_loop { + let last = elapsed + .iter() + .last() + .unwrap() + .try_extract::() + .unwrap(); + time_elapsed = time_elapsed % last; + } + + if let Some(&(mut index)) = time.file_row.as_ref() { + current_time_file = Some(elapsed.get(index).unwrap().try_extract::().unwrap()); + + index += 1; + + if opt_loop { + index = index % df.shape().0; + } + + if let Ok(try_next) = elapsed.get(index) { + let try_next = try_next.try_extract::().unwrap(); + next_tried_time_file = Some(try_next); + + if try_next == time_elapsed { + time.file_row = Some(index); + found = true; + } + } + } + + // attention avec skip car si modulo applied alors on va tout skip sans trouver le bon au début + // let skip = time.file_row.unwrap_or_default(); + let skip = 0; + + if !found { + time.file_row = Some( + elapsed + .iter() + .skip(skip) + .map(|t| t.try_extract::().unwrap()) + .position(|t| time_elapsed == t) + .expect(&format!( + "Elapsed time is not synchronised with file. Current elapsed time is {}{}. \ + From reading the file, it was impossible to find a time matching simulation time. \ + For more info, current: {:?} and next tried: {:?}.", + time.elapsed_seconds(), + opt_loop.then(|| format!( + "(Looping on file, modulo applied to time elapsed: {})", + time_elapsed + )).unwrap_or("".to_string()), + current_time_file, + next_tried_time_file + )), + ); + } + + // println!( + // "{:?} {} {} {}", + // time.file_row, + // time.elapsed_seconds(), + // time_elapsed, + // time.time_step + // ); + } + + fn fn_update_file_index(&self, config: &Config, time: &mut Time) { + let mut files = vec![config.simulation.file.as_ref()]; + for body in &config.bodies { + files.push(body.file_data.as_ref()); + } + + for file in &files { + if let Some(file) = file { + let opt_loop = + file.behavior.clone().unwrap_or(FileBehavior::Stop) == FileBehavior::Loop; + self.fn_update_file_index_core(config, time, file.path.as_ref().unwrap(), opt_loop); + } + } + } + fn fn_update_scene_core(&self, config: &Config, time: &Time, scene: &mut WindowScene) { let elapsed_from_start = time.elapsed_seconds_from_start(); @@ -371,7 +473,7 @@ pub trait Routines: DowncastSync { fn fn_update_body_data( &mut self, - _cfg: &Config, + _config: &Config, _body: usize, _bodies: &mut [AirlessBody], _bodies_data: &mut [BodyData], @@ -380,15 +482,41 @@ pub trait Routines: DowncastSync { ) { } + fn read_file_data(&self, config: &Config, body: usize, time: &Time) -> Vec { + if let Some(file) = config.bodies[body].file_data.as_ref() { + let df = match CsvReader::from_path(&file.path.as_ref().unwrap()) { + Ok(reader) => reader.has_header(true).finish().unwrap(), + Err(e) => panic!("{}", e), + }; + + let row_index = time.file_row.unwrap(); + df.get_row(row_index) + .unwrap() + .0 + .into_iter() + .skip(1) + .map(|v| { + v.cast(&DataType::Float64) + .unwrap() + .try_extract::() + .unwrap() + }) + .collect_vec() + } else { + panic!("No file data found."); + } + } + fn fn_update_body_colormap( &self, + config: &Config, body: usize, bodies: &mut [AirlessBody], pre_computed_bodies: &[BodyData], - cfg: &Config, + time: &Time, win: &Window, ) { - let scalars = match &cfg.window.colormap.scalar { + let scalars = match &config.window.colormap.scalar { Some(CfgScalar::AngleIncidence) => compute_cosine_incidence_angle( &bodies[body], &pre_computed_bodies[body].normals, @@ -408,24 +536,114 @@ pub trait Routines: DowncastSync { ) .map(|a| a.acos() * DPR), None => return, + Some(CfgScalar::File) => { + let v = self.read_file_data(config, body, time); + DRVector::from_row_slice(&v) + } _ => unreachable!(), }; - update_colormap_scalar(win, cfg, scalars.as_slice(), &mut bodies[body], body); + update_colormap_scalar(win, config, scalars.as_slice(), &mut bodies[body], body); } fn fn_update_body( &mut self, - cfg: &Config, + config: &Config, body: usize, bodies: &mut [AirlessBody], bodies_data: &mut [BodyData], time: &Time, window: &Window, ) { - self.fn_update_body_matrix_model(cfg, body, bodies, bodies_data, time); - self.fn_update_body_data(cfg, body, bodies, bodies_data, time, &window.scene.borrow()); - self.fn_update_body_colormap(body, bodies, bodies_data, cfg, window); + self.fn_update_body_matrix_model(config, body, bodies, bodies_data, time); + self.fn_update_body_data( + config, + body, + bodies, + bodies_data, + time, + &window.scene.borrow(), + ); + self.fn_update_body_colormap(config, body, bodies, bodies_data, time, window); + } + + fn fn_export_body_once_core( + &self, + config: &Config, + body: usize, + bodies: &[AirlessBody], + _bodies_data: &[BodyData], + folders: &FoldersRun, + ) { + if config.bodies[body].record.mesh { + let faces = bodies[body].surface.faces.clone(); + let sph = faces.iter().map(|f| f.vertex.sph()).collect_vec(); + let mut df = df!( + "x" => faces.iter().map(|f| f.vertex.position.x).collect_vec(), + "y" => faces.iter().map(|f| f.vertex.position.y).collect_vec(), + "z" => faces.iter().map(|f| f.vertex.position.z).collect_vec(), + "lon" => sph.iter().map(|sph| sph[1]).collect_vec(), + "lat" => sph.iter().map(|sph| sph[2]).collect_vec(), + "rad" => sph.iter().map(|sph| sph[0]).collect_vec(), + ) + .unwrap(); + + let folder_simu = folders.simu_body(&config.bodies[body].name); + fs::create_dir_all(&folder_simu).unwrap(); + + let mut file = std::fs::File::options() + .append(true) + .create(true) + .open(folder_simu.join("mesh.csv")) + .unwrap(); + CsvWriter::new(&mut file) + .include_header(true) + .finish(&mut df) + .unwrap(); + } + + if config.bodies[body].record.depth { + let mut depth = None; + + if let Some(interior) = bodies[body].interior.as_ref() { + match interior { + Interior::Grid(grid) => { + depth = Some(grid.depth.clone()); + } + } + } + + if let Some(depth) = depth { + let mut df = df!( + "depth" => depth, + ) + .unwrap(); + + let folder_simu = folders.simu_body(&config.bodies[body].name); + fs::create_dir_all(&folder_simu).unwrap(); + + let mut file = std::fs::File::options() + .append(true) + .create(true) + .open(folder_simu.join("depth.csv")) + .unwrap(); + CsvWriter::new(&mut file) + .include_header(true) + .finish(&mut df) + .unwrap(); + } + } + } + + fn fn_export_body_once( + &self, + config: &Config, + body: usize, + bodies: &[AirlessBody], + bodies_data: &[BodyData], + folders: &FoldersRun, + ) { + self.fn_export_body_once_core(config, body, bodies, bodies_data, folders); } fn fn_export_iteration( @@ -465,14 +683,17 @@ pub trait Routines: DowncastSync { fn fn_render( &mut self, - _cfg: &Config, + config: &Config, bodies: &mut [AirlessBody], - _time: &Time, + time: &Time, win: &mut Window, ) { // self.win.update_vaos(self.bodies.iter_mut().map(|b| &mut b.asteroid_mut().surface)); win.render_asteroids(&bodies); win.swap_window(); + if let Some(true) = config.preferences.debug.window { + println!("New frame: {}.", time.elapsed_seconds()); + } } } @@ -560,17 +781,12 @@ pub fn read_state_file( time: &Time, columns: &[FileColumns], ) -> Vec { - let mut row_index = time.iteration() * file.row_multiplier.unwrap_or(1); - let df = match CsvReader::from_path(&file.path.as_ref().unwrap()) { - Ok(reader) => reader.has_header(false).finish().unwrap(), + Ok(reader) => reader.has_header(true).finish().unwrap(), Err(e) => panic!("{}", e), }; - if let Some(FileBehavior::Loop) = file.behavior { - row_index = row_index % df.shape().0; - }; - + let row_index = time.file_row.unwrap(); let row = df .get_row(row_index) .unwrap() diff --git a/src/simu/routines/thermal.rs b/src/simu/routines/thermal.rs index 31b7ace..b1d717c 100644 --- a/src/simu/routines/thermal.rs +++ b/src/simu/routines/thermal.rs @@ -75,7 +75,7 @@ impl ThermalBodyData { let dx2 = DRVector::from_row_slice( &depth_grid .iter() - .skip(1) + // .skip(1) .take(depth_size - 1) .tuple_windows() .map(|(a, b)| (b - a).powi(2)) @@ -170,7 +170,7 @@ pub trait RoutinesThermal: Routines { &thermal_data.emissivities, &thermal_data.conductivities, thermal_data.tmp.rows(1, 2).as_view(), - depth_grid[2], + depth_grid[1], ) } @@ -235,7 +235,7 @@ impl Routines for RoutinesThermalDefault { // be in meters for computation unless explicitely mentioned. fn fn_update_body_data( &mut self, - cfg: &Config, + config: &Config, body: usize, bodies: &mut [AirlessBody], bodies_data: &mut [BodyData], @@ -245,8 +245,13 @@ impl Routines for RoutinesThermalDefault { let sun_position = scene.light.position; if time.iteration == 0 { - self.fn_compute_initial_temperatures(&bodies[body], &cfg.bodies[body], &sun_position); - // return; + let tmp = self.fn_compute_initial_temperatures( + &bodies[body], + &config.bodies[body], + &sun_position, + ); + self.data[body].tmp = tmp; + return; } let dt = time.used_time_step(); @@ -294,16 +299,47 @@ impl Routines for RoutinesThermalDefault { .tmp .set_row(curr_size + 1, &temperatures_bottom); + if let Some(true) = config.preferences.debug.thermal_stats { + println!( + "Update: {:.0} SF: {:.1}±({:.1})/{:.1}/{:.1} | T: {:.1}±({:.1})/{:.1}/{:.1}", + time.elapsed_seconds(), + fluxes_solar.mean(), + fluxes_solar.variance().sqrt(), + fluxes_solar.max(), + fluxes_solar.min(), + temperatures_surface.mean(), + temperatures_surface.variance().sqrt(), + temperatures_surface.max(), + temperatures_surface.min() + ); + println!( + "more T0: {:.1}±({:.1})/{:.1}/{:.1} | T1: {:.1}±({:.1})/{:.1}/{:.1} | T2: {:.1}±({:.1})/{:.1}/{:.1}", + self.data[body].tmp.row(0).mean(), + self.data[body].tmp.row(0).variance().sqrt(), + self.data[body].tmp.row(0).max(), + self.data[body].tmp.row(0).min(), + self.data[body].tmp.row(1).mean(), + self.data[body].tmp.row(1).variance().sqrt(), + self.data[body].tmp.row(1).max(), + self.data[body].tmp.row(1).min(), + self.data[body].tmp.row(2).mean(), + self.data[body].tmp.row(2).variance().sqrt(), + self.data[body].tmp.row(2).max(), + self.data[body].tmp.row(2).min() + ); + } + self.data[body].fluxes = fluxes; self.data[body].fluxes_solar = fluxes_solar; } fn fn_update_body_colormap( &self, + cfg: &Config, body: usize, bodies: &mut [AirlessBody], pre_computed_bodies: &[BodyData], - cfg: &Config, + _time: &Time, win: &Window, ) { let cmap = &cfg.window.colormap; @@ -331,6 +367,7 @@ impl Routines for RoutinesThermalDefault { Some(CfgScalar::FluxEmitted) => unimplemented!(), Some(CfgScalar::FluxSelf) => unimplemented!(), Some(CfgScalar::FluxMutual) => unimplemented!(), + Some(CfgScalar::File) => unimplemented!(), None | Some(CfgScalar::Temperature) => self.data[body].tmp.row(0).into_owned(), }; diff --git a/src/simu/scenario.rs b/src/simu/scenario.rs index d4605b6..aedca1b 100644 --- a/src/simu/scenario.rs +++ b/src/simu/scenario.rs @@ -1,11 +1,15 @@ use crate::{ - check_if_latest_version, config, config::CfgRoutines, config::Config, config::InteriorGrid1D, - path_cfg_folder, read_surface_main, util::*, AirlessBody, BodyData, Export, FoldersRun, - FrameEvent, Result, Routines, RoutinesThermalDefault, RoutinesViewerDefault, Time, Window, - KEY_BACKWARD, KEY_FORWARD, KEY_LEFT, KEY_RIGHT, SENSITIVITY, + check_if_latest_version, + config::{self, CfgRoutines, Config, InteriorGrid1D}, + path_cfg_folder, read_surface_main, thermal_skin_depth_one, thermal_skin_depth_two_pi, + util::*, + AirlessBody, BodyData, Export, FoldersRun, FrameEvent, Result, Routines, + RoutinesThermalDefault, RoutinesViewerDefault, Time, Window, KEY_BACKWARD, KEY_FORWARD, + KEY_LEFT, KEY_RIGHT, SENSITIVITY, }; use chrono::Utc; +use config::SkinDepth; use itertools::Itertools; use sdl2::keyboard::Keycode; @@ -50,7 +54,9 @@ impl Scenario { } let mut folders = FoldersRun::new(&config); - folders.save_cfgs(&path_cfg); + // folders.save_cfgs(&path_cfg); + folders.save_cfgs("cfg"); + // folders.save_src(&path_mainrs); let bodies = vec![]; @@ -170,6 +176,26 @@ impl Scenario { InteriorGrid1D::File { path } => { asteroid.with_interior_grid_from_file(path) } + InteriorGrid1D::Increasing { skin, m, n, b } => { + let diffusivity = + asteroid.surface.faces[0].vertex.material.diffusivity(); + let period = cb.spin.period; + let zs = match skin { + SkinDepth::One => thermal_skin_depth_one(diffusivity, period), + SkinDepth::TwoPi => thermal_skin_depth_two_pi(diffusivity, period), + }; + let zmax = zs * *b as Float; + + let mut z = vec![0.0]; + let mut dz = zs / *m as Float; + let mut ii = 0; + while z[ii] < zmax { + dz = dz * (1.0 + 1.0 / *n as Float); + z.push(z[ii] + dz); + ii += 1; + } + asteroid.with_interior_grid_depth(z) + } }, }, }; @@ -200,6 +226,16 @@ impl Scenario { self.routines .init(&self.config, &mut self.bodies, &self.time, &mut self.win); + for body in 0..self.bodies.len() { + self.routines.fn_export_body_once( + &self.config, + body, + &mut self.bodies, + &self.pre_computed_bodies, + &self.folders, + ) + } + 'main_loop: loop { // Register keyboard and mouse interactions. let event = self.win.events(); @@ -223,6 +259,9 @@ impl Scenario { let elapsed = self.time.elapsed_seconds(); let jd = self.time.jd(); + self.routines + .fn_update_file_index(&self.config, &mut self.time); + self.routines.fn_update_scene( &self.config, &self.time, diff --git a/src/simu/time.rs b/src/simu/time.rs index d325bf9..8bbb5ec 100644 --- a/src/simu/time.rs +++ b/src/simu/time.rs @@ -9,6 +9,7 @@ pub struct Time { pub(crate) time_start: usize, pub(crate) real_time: NaiveTime, pub(crate) used_time_step: usize, + pub(crate) file_row: Option, } impl Time { @@ -20,6 +21,7 @@ impl Time { time_start: 0, real_time: Utc::now().time(), used_time_step: 0, + file_row: None, } } @@ -31,6 +33,7 @@ impl Time { time_start: self.time_start, real_time: self.real_time, used_time_step: self.used_time_step, + file_row: self.file_row, } } @@ -42,6 +45,7 @@ impl Time { time_start: start, real_time: self.real_time, used_time_step: self.used_time_step, + file_row: self.file_row, } } diff --git a/src/simu/util.rs b/src/simu/util.rs index 59f8cea..b449a44 100644 --- a/src/simu/util.rs +++ b/src/simu/util.rs @@ -22,9 +22,12 @@ pub fn read_surface(cb: &cg::Body, kind: cg::MeshKind) -> Result { let surface = builder .update_all(|v| { - v.position.x *= mesh.factor.x; - v.position.y *= mesh.factor.y; - v.position.z *= mesh.factor.z; + // v.position.x = mesh.orientation * mesh.factor.x * v.position.x + mesh.position.x; + // v.position.y *= mesh.factor.y; + // v.position.z *= mesh.factor.z; + + v.position = mesh.orientation * v.position.component_mul(&mesh.factor) + mesh.position; + v.material = Material { albedo: cb.material.albedo, emissivity: cb.material.emissivity, diff --git a/src/thermal/bound.rs b/src/thermal/bound.rs index 6fea76c..7069570 100644 --- a/src/thermal/bound.rs +++ b/src/thermal/bound.rs @@ -21,36 +21,55 @@ pub fn newton_method_temperature( depth: Float, ) -> DRVector { let mut index = 0; - let mut old_temperatures = temperatures.clone_owned(); + let mut temperatures = temperatures.clone_owned(); loop { - let new_temperatures = &old_temperatures - - newton_method_temperature_function( - &old_temperatures, - fluxes, - emissivities, - conductivities, - subsurface_temperatures, - depth, - ) - .component_div(&newton_method_temperature_derivative( - &old_temperatures, - emissivities, - conductivities, - depth, - )); + let f = newton_method_temperature_function( + &temperatures, + fluxes, + emissivities, + conductivities, + subsurface_temperatures, + depth, + ); - if index > NUMBER_ITERATION_FAIL - || (&new_temperatures - &old_temperatures) - .abs() - .iter() - .all(|&dt| dt < NEWTON_METHOD_THRESHOLD) - { + let df = newton_method_temperature_derivative( + &temperatures, + emissivities, + conductivities, + depth, + ); + + let delta = (-&f).component_div(&df); + temperatures += δ + + /* + println!( + "Newton f: {:.1}±({:.1})/{:.1}/{:.1} | df: {:.1}±({:.1})/{:.1}/{:.1} | dT: {:.1}±({:.1})/{:.1}/{:.1}", + f.mean(), + f.variance().sqrt(), + f.max(), + f.min(), + df.mean(), + df.variance().sqrt(), + df.max(), + df.min(), + delta.mean(), + delta.variance().sqrt(), + delta.max(), + delta.min(), + ); + */ + + if index > NUMBER_ITERATION_FAIL { + panic!("Newton method never converged."); + } + + if delta.abs().iter().all(|&dt| dt < NEWTON_METHOD_THRESHOLD) { // println!("index: {}", index); - return new_temperatures; + return temperatures; } - old_temperatures = new_temperatures; index += 1; } } @@ -67,7 +86,7 @@ pub fn newton_method_temperature_function( + conductivities.component_mul( &(-subsurface_temperatures.row(1) + 4.0 * subsurface_temperatures.row(0) - 3.0 * temperatures), - ) / depth + ) / (2.0 * depth) } pub fn newton_method_temperature_derivative( @@ -77,5 +96,5 @@ pub fn newton_method_temperature_derivative( depth: Float, ) -> DRVector { -4.0 * STEFAN_BOLTZMANN * emissivities.component_mul(&temperatures.map(|t| t.powi(3))) - - 3.0 / depth * conductivities + - 3.0 / (2.0 * depth) * conductivities } diff --git a/src/util.rs b/src/util.rs index 1f0eba0..6da356d 100644 --- a/src/util.rs +++ b/src/util.rs @@ -109,7 +109,8 @@ pub const RADIANS_PER_SECOND_DMS: Float = RADIANS_PER_MINUTE_DMS / 60.0; pub const NUMBER_ITERATION_FAIL: usize = 1e4 as usize; /// Threshold that defines the convergence condition of the numerical Newton method. -pub const NEWTON_METHOD_THRESHOLD: Float = 1e-4; +// pub const NEWTON_METHOD_THRESHOLD: Float = 1e-4; +pub const NEWTON_METHOD_THRESHOLD: Float = 0.1; pub const VERSION: &str = env!("CARGO_PKG_VERSION");