Functions to calculate themal radiation flux (#145)

* Add `radiation.jl`

* Add test for `radiation.jl`

* Add a fractal topography model for test

* Update radiation.jl

* Update radiation.jl

* Update radiation.jl

* Update blackbody_radiation.jl

Add test for local TPM and thermal radiation calculation

* Update src/radiation.jl

Co-authored-by: Yuto Horikawa <hyrodium@gmail.com>

* Update blackbody_radiation.jl

Add a comment

* Update src/radiation.jl

Co-Authored-By: Yuto Horikawa <hyrodium@gmail.com>

* Update blackbody_radiation.jl

Just add a comment.

* Update radiation.jl

`L` denotes radiation [W/m²].

* Update radiation.jl

Just correct some comments.

* Update blackbody_radiation.jl

Added the value tests of the thermal radiation.

* Rename `test/blackbody_radiation.jl` to `test/thermal_radiation.jl`

* Change function names

- blackbody_radiation -> blackbody_radiance
- thermal_radiation -> thermal_radiance

* Change file names

src/radiation.jl -> src/thermal_radiation.jl

* Update thermal_radiation.jl

Add a test of blackbody radiation

* Update AsteroidThermoPhysicalModels.jl

* Update thermal_radiation.jl

* Update thermal_radiation.jl

---------

Co-authored-by: Yuto Horikawa <hyrodium@gmail.com>

.gitignore

@@ -15,8 +15,7 @@
 /test/TPM_Didymos/output
 
 # Keep these shape models
+!fractal_v2572_f5000.obj
 !icosahedron.obj
-!single_face.obj
 !ryugu_test.obj
-
-
+!single_face.obj

src/AsteroidThermoPhysicalModels.jl

@@ -18,6 +18,8 @@ import GeometryBasics
 const SOLAR_CONST = 1366.0   # Solar constant, Φ☉ [W/m^2]
 const c₀ = 299792458.0       # Speed of light [m/s]
 const σ_SB = 5.670374419e-8  # Stefan–Boltzmann constant [W/m^2/K^4]
+const h = 6.62607015e-34     # Planck constant [J⋅s]
+const k_B = 1.380649e-23     # Boltzmann constant [J/K]
 
 include("obj.jl")
 include("shape.jl")
@@ -29,6 +31,7 @@ include("TPM.jl")
 include("heat_conduction.jl")
 include("energy_flux.jl")
 include("non_grav.jl")
+include("thermal_radiation.jl")
 export thermal_skin_depth, thermal_inertia, init_temperature!, run_TPM!
 
 include("roughness.jl")

src/thermal_radiation.jl

@@ -0,0 +1,75 @@
+"""
+    blackbody_radiance(λ, T) -> L_λ
+
+Accroding to Planck's law, calculate the spectral intensity of blackbody radiation
+at wavelength `λ` and temperature `T`.
+
+# Arguments
+- `λ` : Wavelength [m]
+- `T` : Temperature [K]
+
+# Return
+- `L_λ` : Spectral radiance [W/m²/m/steradian]
+
+cf. https://github.com/JuliaAstro/Planck.jl/blob/main/src/Planck.jl
+"""
+blackbody_radiance(λ, T) = 2 * h * c₀^2 / λ^5 / expm1(h * c₀ / (λ * k_B * T))
+
+"""
+    blackbody_radiance(T) -> L
+
+According to Stefan-Boltzmann law, calculate the total radiance of blackbody radiation
+at temperature `T`, integrated over all wavelength.
+
+# Arguments
+- `T` : Temperature [K]
+
+# Return
+- `L` : Radiance [W/m²/steradian]
+"""
+blackbody_radiance(T) = σ_SB * T^4  ## TODO: テストを追加する
+
+
+"""
+    thermal_radiance(shape, emissivities, temperatures, obs) -> L
+
+Calculate the radiance from the temperature distribution based on a shape model.
+
+# Arguments
+- `shape`        : Shape model of an asteroid
+- `emissivities` : Emissivity of each facet of the shape model [-]
+- `temperatures` : Temperature of each facet of the shape model [K]
+- `obs`          : Position vector of the observer in the same coordinate system as `shape` [m]
+
+# Return
+- `L` : Radiance [W/m²]
+"""
+function thermal_radiance(shape::ShapeModel, emissivities::AbstractVector{<:Real}, temperatures::AbstractVector{<:Real}, obs::StaticVector{3, <:Real})
+    if length(temperatures) != length(shape.faces)
+        throw(ArgumentError("Length of `temperatures` must be equal to the number of faces of the shape model."))
+    end
+
+    L = 0.0  # Radiance [W/m²]
+    for i in eachindex(shape.faces)
+        c = shape.face_centers[i]
+        n̂ = shape.face_normals[i]
+        a = shape.face_areas[i]
+
+        ε = emissivities[i]
+        T = temperatures[i]
+
+        ## Direction and distance from the facet center to the observer
+        d̂ = normalize(obs - c)
+        d = norm(obs - c)
+
+        cosθ = n̂ ⋅ d̂
+        cosθ < 0 && continue  # Ignore the facet if the observer is below the horizon.
+
+        ## TODO: Ray-trace for each facet
+
+        ## TODO: Consider the solid angle of the facet.
+        ## The coefficient may change from π to 2π or 4π.
+        L += ε * σ_SB * T^4 * a * cosθ / (π * d^2)
+    end
+    return L
+end

test/runtests.jl

@@ -12,9 +12,11 @@ using DataFrames
 
 Aqua.test_all(AsteroidThermoPhysicalModels, ambiguities=false)
 
-include("find_visiblefacets.jl")
 include("TPM_Ryugu/TPM_Ryugu.jl")
 include("TPM_non-uniform_thermoparams/TPM_non-uniform_thermoparams.jl")
 include("TPM_zero-conductivity/TPM_zero-conductivity.jl")
 include("TPM_Didymos/TPM_Didymos.jl")
+
+include("thermal_radiation.jl")
+include("find_visiblefacets.jl")
 include("heat_conduction_1D.jl")