Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add _getindex and remove some ternary operators #144

Merged
merged 1 commit into from
Jul 6, 2024
Merged

add _getindex and remove some ternary operators #144

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions src/AsteroidThermoPhysicalModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,10 @@ 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]

# utils for NonUniformThermoParams and UniformThermoParams
_getindex(a::AbstractArray, i::Integer) = a[i]
_getindex(a::Real, ::Integer) = a

include("obj.jl")
include("shape.jl")
include("facet.jl")
Expand Down
24 changes: 12 additions & 12 deletions src/energy_flux.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@ Input energy per second on the whole surface [W]
function energy_in(stpm::SingleTPM)
E_in = 0.
for i in eachindex(stpm.shape.faces)
R_vis = (stpm.thermo_params.reflectance_vis isa Real ? stpm.thermo_params.reflectance_vis : stpm.thermo_params.reflectance_vis[i])
R_ir = (stpm.thermo_params.reflectance_ir isa Real ? stpm.thermo_params.reflectance_ir : stpm.thermo_params.reflectance_ir[i])
R_vis = _getindex(stpm.thermo_params.reflectance_vis, i)
R_ir = _getindex(stpm.thermo_params.reflectance_ir, i)
F_sun = stpm.flux[i, 1]
F_scat = stpm.flux[i, 2]
F_rad = stpm.flux[i, 3]
Expand All @@ -51,7 +51,7 @@ Output enegey per second from the whole surface [W]
function energy_out(stpm::SingleTPM)
E_out = 0.
for i in eachindex(stpm.shape.faces)
ε = (stpm.thermo_params.emissivity isa Real ? stpm.thermo_params.emissivity : stpm.thermo_params.emissivity[i])
ε = _getindex(stpm.thermo_params.emissivity, i)
T = stpm.temperature[begin, i] # Surface temperature
a = stpm.shape.face_areas[i]

Expand Down Expand Up @@ -151,7 +151,7 @@ function update_flux_scat_single!(stpm::SingleTPM)
for visiblefacet in stpm.shape.visiblefacets[i_face]
j = visiblefacet.id
fᵢⱼ = visiblefacet.f
R_vis = (stpm.thermo_params.reflectance_vis isa Real ? stpm.thermo_params.reflectance_vis : stpm.thermo_params.reflectance_vis[j])
R_vis = _getindex(stpm.thermo_params.reflectance_vis, j)

stpm.flux[i_face, 2] += fᵢⱼ * R_vis * stpm.flux[j, 1]
end
Expand Down Expand Up @@ -204,8 +204,8 @@ function update_flux_rad_single!(stpm::SingleTPM)
for visiblefacet in stpm.shape.visiblefacets[i]
j = visiblefacet.id
fᵢⱼ = visiblefacet.f
ε = (stpm.thermo_params.emissivity isa Real ? stpm.thermo_params.emissivity : stpm.thermo_params.emissivity[j])
R_ir = (stpm.thermo_params.reflectance_ir isa Real ? stpm.thermo_params.reflectance_ir : stpm.thermo_params.reflectance_ir[j])
ε = _getindex(stpm.thermo_params.emissivity, j)
R_ir = _getindex(stpm.thermo_params.reflectance_ir, j)
Tⱼ = stpm.temperature[begin, j]

stpm.flux[i, 3] += ε * σ_SB * (1 - R_ir) * fᵢⱼ * Tⱼ^4
Expand Down Expand Up @@ -421,12 +421,12 @@ function mutual_heating!(btpm::BinaryTPM, rₛ, R₂₁)
T₁ = btpm.pri.temperature[begin, i]
T₂ = btpm.sec.temperature[begin, j]

ε₁ = (thermo_params1.emissivity isa Real ? thermo_params1.emissivity : thermo_params1.emissivity[i])
ε₂ = (thermo_params2.emissivity isa Real ? thermo_params2.emissivity : thermo_params2.emissivity[j])
R_vis₁ = (thermo_params1.reflectance_vis isa Real ? thermo_params1.reflectance_vis : thermo_params1.reflectance_vis[i])
R_vis₂ = (thermo_params2.reflectance_vis isa Real ? thermo_params2.reflectance_vis : thermo_params2.reflectance_vis[j])
R_ir₁ = (thermo_params1.reflectance_ir isa Real ? thermo_params1.reflectance_ir : thermo_params1.reflectance_ir[i])
R_ir₂ = (thermo_params2.reflectance_ir isa Real ? thermo_params2.reflectance_ir : thermo_params2.reflectance_ir[j])
ε₁ = _getindex(thermo_params1.emissivity, i)
ε₂ = _getindex(thermo_params2.emissivity, j)
R_vis₁ = _getindex(thermo_params1.reflectance_vis, i)
R_vis₂ = _getindex(thermo_params2.reflectance_vis, j)
R_ir₁ = _getindex(thermo_params1.reflectance_ir, i)
R_ir₂ = _getindex(thermo_params2.reflectance_ir, j)

## Mutual heating by scattered light
btpm.pri.flux[i, 2] += f₁₂ * R_vis₂ * btpm.sec.flux[j, 1]
Expand Down
24 changes: 12 additions & 12 deletions src/heat_conduction.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,9 +67,9 @@ function forward_euler!(stpm::SingleTPM, Δt)
## Zero-conductivity (thermal inertia) case
if iszero(stpm.thermo_params.inertia)
for i_face in 1:n_face
R_vis = (stpm.thermo_params.reflectance_vis isa Real ? stpm.thermo_params.reflectance_vis : stpm.thermo_params.reflectance_vis[i_face] )
R_ir = (stpm.thermo_params.reflectance_ir isa Real ? stpm.thermo_params.reflectance_ir : stpm.thermo_params.reflectance_ir[i_face])
ε = (stpm.thermo_params.emissivity isa Real ? stpm.thermo_params.emissivity : stpm.thermo_params.emissivity[i_face])
R_vis = _getindex(stpm.thermo_params.reflectance_vis, i_face)
R_ir = _getindex(stpm.thermo_params.reflectance_ir, i_face)
ε = _getindex(stpm.thermo_params.emissivity, i_face)
εσ = ε * σ_SB

F_sun, F_scat, F_rad = stpm.flux[i_face, :]
Expand All @@ -82,7 +82,7 @@ function forward_euler!(stpm::SingleTPM, Δt)
for i_face in 1:n_face
P = stpm.thermo_params.period
Δz = stpm.thermo_params.Δz
l = (stpm.thermo_params.skindepth isa Real ? stpm.thermo_params.skindepth : stpm.thermo_params.skindepth[i_face])
l = _getindex(stpm.thermo_params.skindepth, i_face)

λ = (Δt/P) / (Δz/l)^2 / 4π
λ ≥ 0.5 && error("The forward Euler method is unstable because λ = $λ. This should be less than 0.5.")
Expand Down Expand Up @@ -116,7 +116,7 @@ function backward_euler!(stpm::SingleTPM, Δt)
# n_face = size(T, 2)

# for i_face in 1:n_face
# λ = (stpm.thermo_params.λ isa Real ? stpm.thermo_params.λ : stpm.thermo_params.λ[i_face])
# λ = _getindex(stpm.thermo_params.λ, i_face)

# stpm.SOLVER.a .= -λ
# stpm.SOLVER.a[begin] = 0
Expand Down Expand Up @@ -242,13 +242,13 @@ function update_upper_temperature!(stpm::SingleTPM, i::Integer)

#### Radiation boundary condition ####
if stpm.BC_UPPER isa RadiationBoundaryCondition
P = stpm.thermo_params.period
l = (stpm.thermo_params.skindepth isa Real ? stpm.thermo_params.skindepth : stpm.thermo_params.skindepth[i] )
Γ = (stpm.thermo_params.inertia isa Real ? stpm.thermo_params.inertia : stpm.thermo_params.inertia[i])
R_vis = (stpm.thermo_params.reflectance_vis isa Real ? stpm.thermo_params.reflectance_vis : stpm.thermo_params.reflectance_vis[i] )
R_ir = (stpm.thermo_params.reflectance_ir isa Real ? stpm.thermo_params.reflectance_ir : stpm.thermo_params.reflectance_ir[i])
ε = (stpm.thermo_params.emissivity isa Real ? stpm.thermo_params.emissivity : stpm.thermo_params.emissivity[i])
Δz = stpm.thermo_params.Δz
P = stpm.thermo_params.period
l = _getindex(stpm.thermo_params.skindepth, i)
Γ = _getindex(stpm.thermo_params.inertia, i)
R_vis = _getindex(stpm.thermo_params.reflectance_vis, i)
R_ir = _getindex(stpm.thermo_params.reflectance_ir, i)
ε = _getindex(stpm.thermo_params.emissivity, i)
Δz = stpm.thermo_params.Δz

F_sun, F_scat, F_rad = stpm.flux[i, :]
F_total = flux_total(R_vis, R_ir, F_sun, F_scat, F_rad)
Expand Down
8 changes: 4 additions & 4 deletions src/non_grav.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,10 @@ function update_thermal_force!(stpm::SingleTPM)

## Total emittance from face i , Eᵢ [W/m²].
## Note that both scattered light and thermal radiation are assumed to be isotropic.
R_vis = (stpm.thermo_params.reflectance_vis isa Real ? stpm.thermo_params.reflectance_vis : stpm.thermo_params.reflectance_vis[i])
R_ir = (stpm.thermo_params.reflectance_ir isa Real ? stpm.thermo_params.reflectance_ir : stpm.thermo_params.reflectance_ir[i])
ε = (stpm.thermo_params.emissivity isa Real ? stpm.thermo_params.emissivity : stpm.thermo_params.emissivity[i])
Eᵢ = R_vis * F_sun + R_vis * F_scat + R_ir * F_rad + ε * σ_SB * Tᵢ^4
R_vis = _getindex(stpm.thermo_params.reflectance_vis, i)
R_ir = _getindex(stpm.thermo_params.reflectance_ir, i)
ε = _getindex(stpm.thermo_params.emissivity, i)
Eᵢ = R_vis * F_sun + R_vis * F_scat + R_ir * F_rad + ε * σ_SB * Tᵢ^4

## Thermal force on each face
stpm.face_forces[i] = - 2/3 * Eᵢ * aᵢ / c₀ * n̂ᵢ # The first term normal to the face
Expand Down
Loading