An example of code is as follows.

struct SingleTPM{P<:AbstractThermoParams, S<:HeatConductionSolver ,T_UPPER<:BoundaryCondition, T_LOWER<:BoundaryCondition} <: ThermoPhysicalModel
    shape          ::ShapeModel
    thermo_params  ::P

    flux           ::Matrix{Float64}  # (n_face, 3)
    temperature    ::Matrix{Float64}  # (n_depth, n_face)

    face_forces    ::Vector{SVector{3, Float64}}
    force          ::MVector{3, Float64}
    torque         ::MVector{3, Float64}

    SELF_SHADOWING ::Bool
    SELF_HEATING   ::Bool
    SOLVER         ::S
    BC_UPPER       ::T_UPPER
    BC_LOWER       ::T_LOWER
end

• Create a periodic ephemerides with a rotation period
• Computation of shadow is performed only once and stored.
• Save the simulation data at the final rotation

以下のようなものは早めに取り組めるかと思いました。

• 参考文献
• 図(コードの理解に役立つもの。図の説明まで言語化できると良いが、最初から完璧でなくてもOK)

Need to fix function update_thermal_force!. Total emittance from facet i should be E = A_B * F_sun + A_B * F_scat + A_TH * F_rad + εσT^4.

• loadobjのstaticは無意味なキーワード引数のように見えました。
  

  AsteroidThermoPhysicalModels.jl/src/obj.jl

  Line 12 in b179aee

  function loadobj(shapepath::String; scale=1, static=true, message=true)

• loadobjのmessageがBoolなのは変数名からわかりにくい(Stringを期待してしまう)ので、enable_messageのような名前の方が良さそうに思います。
  

  AsteroidThermoPhysicalModels.jl/src/obj.jl

  Line 12 in b179aee

  function loadobj(shapepath::String; scale=1, static=true, message=true)

Add [-] to dimensionless values.

This line

can be replaced with Unitful features.

julia> using SPICE, Unitful

julia> P = SPICE.convrt(7.63262, "hours", "seconds")
27477.432

julia> uconvert(u"s", 7.63262u"hr")
27477.432 s

julia> convert(Float64, (7.63262u"hr") / (1u"s"))
27477.432

配列の長さにnを使ってインデックスにiを使うことが多いです。なのでit, iz, nt, nzのように下付き文字を避けるのが良いかも知れないですね。
(下付き文字は制限が大きいので、特別な意味を込めたいときに使うことが多いです。)

#65ですが、A_B = (A_B isa Real ? A_B : A_B[j]) は型安定でないので避けた方が良さそうですね。後のPRをすべてmergeした後に修正するか悩ましいですが。。
f031538#diff-c540b482ce6af68518b0196d98161b7c55cc1fe6852e390ef64e29836f54f4a3R388