-
Notifications
You must be signed in to change notification settings - Fork 5
/
TPM_Ryugu.jl
107 lines (90 loc) · 4.02 KB
/
TPM_Ryugu.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
# See https://github.com/Astroshaper/Astroshaper-examples/tree/main/TPM_Ryugu for more information.
@testset "TPM_Ryugu" begin
msg = """\n
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
| Test: TPM_Ryugu |
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
"""
println(msg)
##= Download Files =##
paths_kernel = [
"lsk/naif0012.tls",
"pck/hyb2_ryugu_shape_v20190328.tpc",
"fk/hyb2_ryugu_v01.tf",
"spk/2162173_Ryugu.bsp",
]
paths_shape = [
"SHAPE_SFM_49k_v20180804.obj",
]
for path_kernel in paths_kernel
url_kernel = "https://data.darts.isas.jaxa.jp/pub/hayabusa2/spice_bundle/spice_kernels/$(path_kernel)"
filepath = joinpath("kernel", path_kernel)
mkpath(dirname(filepath))
isfile(filepath) || Downloads.download(url_kernel, filepath)
end
for path_shape in paths_shape
url_shape = "https://data.darts.isas.jaxa.jp/pub/hayabusa2/paper/Watanabe_2019/$(path_shape)"
filepath = joinpath("shape", path_shape)
mkpath(dirname(filepath))
isfile(filepath) || Downloads.download(url_shape, filepath)
end
##= Load data with SPICE =##
for path_kernel in paths_kernel
filepath = joinpath("kernel", path_kernel)
SPICE.furnsh(filepath)
end
##= Ephemerides =##
P = SPICE.convrt(7.63262, "hours", "seconds") # Rotation period of Ryugu
et_begin = SPICE.utc2et("2018-07-01T00:00:00") # Start time of TPM
et_end = et_begin + 2P # End time of TPM
step = P / 360 # Time step of TPM, corresponding to 1 deg rotation
et_range = et_begin : step : et_end
"""
- `time` : Ephemeris times
- `sun` : Sun's position in the RYUGU_FIXED frame
"""
ephem = (
time = collect(et_range),
sun = [SVector{3}(SPICE.spkpos("SUN", et, "RYUGU_FIXED", "None", "RYUGU")[1]) * 1000 for et in et_range],
)
SPICE.kclear()
##= Load obj file =##
path_obj = joinpath("shape", "ryugu_test.obj") # Small model for test
# path_obj = joinpath("shape", "SHAPE_SFM_49k_v20180804.obj")
shape = AsteroidThermoPhysicalModels.load_shape_obj(path_obj; scale=1000, find_visible_facets=true)
##= Thermal properties =##
k = 0.1
ρ = 1270.0
Cₚ = 600.0
l = AsteroidThermoPhysicalModels.thermal_skin_depth(P, k, ρ, Cₚ)
Γ = AsteroidThermoPhysicalModels.thermal_inertia(k, ρ, Cₚ)
thermo_params = AsteroidThermoPhysicalModels.thermoparams(
P = P,
l = l,
Γ = Γ,
A_B = 0.04, # Bolometric Bond albedo
A_TH = 0.0,
ε = 1.0,
z_max = 0.6,
Nz = 41,
)
println(thermo_params)
##= Setting of TPM =##
stpm = AsteroidThermoPhysicalModels.SingleTPM(shape, thermo_params;
SELF_SHADOWING = true,
SELF_HEATING = true,
SOLVER = AsteroidThermoPhysicalModels.ForwardEulerSolver(thermo_params),
BC_UPPER = AsteroidThermoPhysicalModels.RadiationBoundaryCondition(),
BC_LOWER = AsteroidThermoPhysicalModels.InsulationBoundaryCondition(),
)
AsteroidThermoPhysicalModels.init_temperature!(stpm, 200)
##= Run TPM =##
time_begin = ephem.time[end] - P # Time to start storing temperature
time_end = ephem.time[end] # Time to end storing temperature
face_ID = [1, 2, 3, 4, 10] # Face indices at which you want to save underground temperature
result = AsteroidThermoPhysicalModels.run_TPM!(stpm, ephem, time_begin, time_end, face_ID)
##= Save TPM result =##
savedir = "TPM_Ryugu"
mkpath(savedir)
AsteroidThermoPhysicalModels.export_TPM_results(savedir, result, stpm, ephem)
end