From 86a665bfed8b253320e3ce2547f20cfd4225239e Mon Sep 17 00:00:00 2001 From: sunitisanghavi Date: Mon, 9 Oct 2023 21:36:01 -0700 Subject: [PATCH] linearize manually in prog. --- Project.toml | 1 + myfile.txt | 1 - src/CoreRT/CoreKernel/lin_doubling.jl | 209 +++++ src/CoreRT/CoreKernel/lin_elemental.jl | 417 ++++++++++ src/CoreRT/CoreKernel/lin_interaction.jl | 262 ++++++ src/CoreRT/CoreKernel/lin_rt_kernel.jl | 413 +++++++++ .../compEffectiveLayerProperties.jl | 15 +- .../lin_compEffectiveLayerProperties.jl | 213 +++++ src/CoreRT/Surfaces/lin_lambertian_surface.jl | 198 +++++ src/CoreRT/lin_model_from_parameters.jl | 2 + src/CoreRT/lin_types.jl | 786 ++++++++++++++++++ src/CoreRT/rt_helper_functions.jl | 4 +- src/CoreRT/rt_run.jl | 202 ++--- src/Scattering/lin_compute_NAI2.jl | 6 +- src/Scattering/lin_truncate_phase.jl | 189 ++--- src/Scattering/lin_types.jl | 92 ++ src/Scattering/types.jl | 2 + test/benchmarks/prototype_inelastic_OCO2.jl | 12 +- test/benchmarks/testAerosol.jl | 22 +- 19 files changed, 2812 insertions(+), 234 deletions(-) delete mode 100644 myfile.txt create mode 100644 src/CoreRT/CoreKernel/lin_doubling.jl create mode 100644 src/CoreRT/CoreKernel/lin_elemental.jl create mode 100644 src/CoreRT/CoreKernel/lin_interaction.jl create mode 100644 src/CoreRT/CoreKernel/lin_rt_kernel.jl create mode 100644 src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl create mode 100644 src/CoreRT/Surfaces/lin_lambertian_surface.jl create mode 100644 src/CoreRT/lin_types.jl create mode 100644 src/Scattering/lin_types.jl diff --git a/Project.toml b/Project.toml index 7ff35fb1..df14844b 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ForwardDiffChainRules = "c9556dd2-1aed-4cfe-8560-1557cf593001" +ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" InstrumentOperator = "9e589c1b-9e01-4e00-831a-aa39ce86e3ef" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" diff --git a/myfile.txt b/myfile.txt deleted file mode 100644 index 0ef9b9b1..00000000 --- a/myfile.txt +++ /dev/null @@ -1 +0,0 @@ -0.523599 1.308997 0.000000 0.816923 -1.0000000.523599 1.308997 0.000000 0.816923 -1.000000 \ No newline at end of file diff --git a/src/CoreRT/CoreKernel/lin_doubling.jl b/src/CoreRT/CoreKernel/lin_doubling.jl new file mode 100644 index 00000000..6ad510ea --- /dev/null +++ b/src/CoreRT/CoreKernel/lin_doubling.jl @@ -0,0 +1,209 @@ +#= + +This file contains RT doubling-related functions + +=# + +""" + $(FUNCTIONNAME)(pol_type, SFI, expk, ndoubl::Int, added_layer::AddedLayer, I_static::AbstractArray{FT}, + architecture) where {FT} + +Compute homogenous layer matrices from its elemental layer using Doubling +""" +function doubling_helper!(pol_type, + SFI, + expk, + dexpk_fctr, + ndoubl::Int, + added_layer::M, + lin_added_layer::M, + I_static::AbstractArray{FT}, + architecture) where {FT,M} + + # Unpack the added layer + @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, j₀⁺, j₀⁻ = added_layer + @unpack dr⁺⁻, dr⁻⁺, dt⁻⁻, dt⁺⁺, dj₀⁺, dj₀⁻ = lin_added_layer + #@show typeof(expk), typeof(I_static) + # Device architecture + dev = devi(architecture) + + # Note: short-circuit evaluation => return nothing evaluated iff ndoubl == 0 + ndoubl == 0 && return nothing + + # Geometric progression of reflections (1-RR)⁻¹ + gp_refl = similar(t⁺⁺) + tt⁺⁺_gp_refl = similar(t⁺⁺) + + dgp_refl = similar(dt⁺⁺) + dtt⁺⁺_gp_refl = similar(dt⁺⁺) + + # Dummy for source + j₁⁺ = similar(j₀⁺) + dj₁⁺ = similar(dj₀⁺) + # Dummy for J + j₁⁻ = similar(j₀⁻) + dj₁⁻ = similar(dj₀⁻) + + # Loop over number of doublings + for n = 1:ndoubl + # T⁺⁺(λ)[I - R⁺⁻(λ)R⁻⁺(λ)]⁻¹, for doubling R⁺⁻,R⁻⁺ and T⁺⁺,T⁻⁻ is identical + #@show typeof(gp_refl), typeof(I_static), typeof(I_static .- r⁻⁺), typeof(j₁⁺) + @timeit "Batch Inv Doubling" batch_inv!(gp_refl, I_static .- r⁻⁺ ⊠ r⁻⁺) + tt⁺⁺_gp_refl .= t⁺⁺ ⊠ gp_refl + for ctr=1:3 + tmp = -(dr⁻⁺[ctr,:,:,:] ⊠ r⁻⁺ + r⁻⁺ ⊠ dr⁻⁺[ctr,:,:,:]) + dgp_refl[ctr,:,:,:] .= tmp ⊠ gp_refl ⊠ tmp + dtt⁺⁺_gp_refl[ctr,:,:,:] .= dt⁺⁺[ctr,:,:,:] ⊠ gp_refl + + t⁺⁺ ⊠ dgp_refl[ctr,:,:,:] + end + + + # J⁺₂₁(λ) = J⁺₁₀(λ).exp(-τ(λ)/μ₀) + @inbounds @views j₁⁺[:,1,:] .= j₀⁺[:,1,:] .* expk' + # J⁻₁₂(λ) = J⁻₀₁(λ).exp(-τ(λ)/μ₀) + @inbounds @views j₁⁻[:,1,:] .= j₀⁻[:,1,:] .* expk' + for ctr=1:4 + dj₁⁺[ctr,:,1,:] .= (dj₀⁺[ctr,:,1,:] + j₀⁺[:,1,:]*dexpk_fctr) .* expk' + # J⁻₁₂(λ) = J⁻₀₁(λ).exp(-τ(λ)/μ₀) + dj₁⁻[ctr,:,1,:] .= (dj₀⁻[ctr,:,1,:] + j₀⁻[:,1,:]*dexpk_fctr) .* expk' + + # J⁻₀₂(λ) = J⁻₀₁(λ) + T⁻⁻₀₁(λ)[I - R⁻⁺₂₁(λ)R⁺⁻₀₁(λ)]⁻¹[J⁻₁₂(λ) + R⁻⁺₂₁(λ)J⁺₁₀(λ)] (see Eqs.8 in Raman paper draft) + dj₀⁻[ctr,:,1,:] .= dj₀⁻[ctr,:,1,:] + + (dtt⁺⁺_gp_refl[ctr,:,:,:] ⊠ (j₁⁻ + r⁻⁺ ⊠ j₀⁺)) + + (tt⁺⁺_gp_refl ⊠ (dj₁⁻[ctr,:,1,:] + dr⁻⁺[ctr,:,:,:] ⊠ j₀⁺ + r⁻⁺ ⊠ dj₀⁺[ctr,:,1,:])) + # J⁺₂₀(λ) = J⁺₂₁(λ) + T⁺⁺₂₁(λ)[I - R⁺⁻₀₁(λ)R⁻⁺₂₁(λ)]⁻¹[J⁺₁₀(λ) + R⁺⁻₀₁(λ)J⁻₁₂(λ)] (see Eqs.8 in Raman paper draft) + dj₀⁺[ctr,:,1,:] .= dj₁⁺[ctr,:,1,:] + + (dtt⁺⁺_gp_refl[ctr,:,:,:] ⊠ (j₀⁺ + r⁻⁺ ⊠ j₁⁻)) + + (tt⁺⁺_gp_refl ⊠ (dj₀⁺[ctr,:,1,:] + dr⁻⁺[ctr,:,:,:] ⊠ j₁⁻ + r⁻⁺ ⊠ dj₁⁻[ctr,:,1,:])) + end + + # J⁻₀₂(λ) = J⁻₀₁(λ) + T⁻⁻₀₁(λ)[I - R⁻⁺₂₁(λ)R⁺⁻₀₁(λ)]⁻¹[J⁻₁₂(λ) + R⁻⁺₂₁(λ)J⁺₁₀(λ)] (see Eqs.8 in Raman paper draft) + j₀⁻ .= j₀⁻ + (tt⁺⁺_gp_refl ⊠ (j₁⁻ + r⁻⁺ ⊠ j₀⁺)) + + # J⁺₂₀(λ) = J⁺₂₁(λ) + T⁺⁺₂₁(λ)[I - R⁺⁻₀₁(λ)R⁻⁺₂₁(λ)]⁻¹[J⁺₁₀(λ) + R⁺⁻₀₁(λ)J⁻₁₂(λ)] (see Eqs.8 in Raman paper draft) + j₀⁺ .= j₁⁺ + (tt⁺⁺_gp_refl ⊠ (j₀⁺ + r⁻⁺ ⊠ j₁⁻)) + expk .= expk.^2 + + # R⁻⁺₂₀(λ) = R⁻⁺₁₀(λ) + T⁻⁻₀₁(λ)[I - R⁻⁺₂₁(λ)R⁺⁻₀₁(λ)]⁻¹R⁻⁺₂₁(λ)T⁺⁺₁₀(λ) (see Eqs.8 in Raman paper draft) + r⁻⁺ .= r⁻⁺ + (tt⁺⁺_gp_refl ⊠ r⁻⁺ ⊠ t⁺⁺) + + # T⁺⁺₂₀(λ) = T⁺⁺₂₁(λ)[I - R⁺⁻₀₁(λ)R⁻⁺₂₁(λ)]⁻¹T⁺⁺₁₀(λ) (see Eqs.8 in Raman paper draft) + t⁺⁺ .= tt⁺⁺_gp_refl ⊠ t⁺⁺ + for ctr=1:3 + # R⁻⁺₂₀(λ) = R⁻⁺₁₀(λ) + T⁻⁻₀₁(λ)[I - R⁻⁺₂₁(λ)R⁺⁻₀₁(λ)]⁻¹R⁻⁺₂₁(λ)T⁺⁺₁₀(λ) (see Eqs.8 in Raman paper draft) + dr⁻⁺[ctr,:,:,:] .= dr⁻⁺[ctr,:,:,:] + + (dtt⁺⁺_gp_refl[ctr,:,:,:] ⊠ r⁻⁺ ⊠ t⁺⁺) + + (tt⁺⁺_gp_refl ⊠ dr⁻⁺[ctr,:,:,:] ⊠ t⁺⁺) + + (tt⁺⁺_gp_refl ⊠ r⁻⁺ ⊠ dt⁺⁺[ctr,:,:,:]) + + # T⁺⁺₂₀(λ) = T⁺⁺₂₁(λ)[I - R⁺⁻₀₁(λ)R⁻⁺₂₁(λ)]⁻¹T⁺⁺₁₀(λ) (see Eqs.8 in Raman paper draft) + dt⁺⁺[ctr,:,:,:] .= dtt⁺⁺_gp_refl[ctr,:,:,:] ⊠ t⁺⁺ + + tt⁺⁺_gp_refl ⊠ dt⁺⁺[ctr,:,:,:] + end + end + synchronize_if_gpu() + + # After doubling, revert D(DR)->R, where D = Diagonal{1,1,-1,-1} + apply_D_matrix!(pol_type.n, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, dr⁻⁺, dt⁺⁺, dr⁺⁻, dt⁻⁻) + + # For SFI, after doubling, revert D(DJ₀⁻)->J₀⁻ + apply_D_matrix_SFI!(pol_type.n, j₀⁻, dj₀⁻) + + return nothing +end + +function doubling!(pol_type, SFI, + expk, dexpk_fctr, + ndoubl::Int, + added_layer::AddedLayer{M},#{FT}, + lin_added_layer::linAddedLayer{M},#{FT}, + I_static::AbstractArray{FT}, + architecture) where {FT,M} + + doubling_helper!(pol_type, SFI, + expk, dexpk_fctr, + ndoubl, + added_layer, + lin_added_layer, + I_static, architecture) + synchronize_if_gpu() +end + +@kernel function apply_D!(n_stokes::Int, + r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, + dr⁻⁺, dt⁺⁺, dr⁺⁻, dt⁻⁻) + iμ, jμ, n = @index(Global, NTuple) + i = mod(iμ, n_stokes) + j = mod(jμ, n_stokes) + + if (i > 2) + r⁻⁺[iμ,jμ,n] = - r⁻⁺[iμ, jμ,n] + end + + if ((i <= 2) & (j <= 2)) | ((i > 2) & (j > 2)) + r⁺⁻[iμ,jμ,n] = r⁻⁺[iμ,jμ,n] + t⁻⁻[iμ,jμ,n] = t⁺⁺[iμ,jμ,n] + else + r⁺⁻[iμ,jμ,n] = - r⁻⁺[iμ,jμ,n] + t⁻⁻[iμ,jμ,n] = - t⁺⁺[iμ,jμ,n] + end + + for ctr=1:3 + if (i > 2) + dr⁻⁺[ctr,iμ,jμ,n] = - dr⁻⁺[ctr,iμ, jμ,n] + end + + if ((i <= 2) & (j <= 2)) | ((i > 2) & (j > 2)) + dr⁺⁻[ctr,iμ,jμ,n] = dr⁻⁺[ctr,iμ,jμ,n] + dt⁻⁻[ctr,iμ,jμ,n] = dt⁺⁺[ctr,iμ,jμ,n] + else + dr⁺⁻[ctr,iμ,jμ,n] = - dr⁻⁺[ctr,iμ,jμ,n] + dt⁻⁻[ctr,iμ,jμ,n] = - dt⁺⁺[ctr,iμ,jμ,n] + end + end + +end + +@kernel function apply_D_SFI!(n_stokes::Int, J₀⁻, dJ₀⁻) + iμ, _, n = @index(Global, NTuple) + i = mod(iμ, n_stokes) + if (i > 2) + J₀⁻[iμ, 1, n] = - J₀⁻[iμ, 1, n] + dJ₀⁻[1:4, iμ, 1, n] .= - dJ₀⁻[1:4, iμ, 1, n] + end +end + +function apply_D_matrix!(n_stokes::Int, + r⁻⁺::AbstractArray{FT,3}, t⁺⁺::AbstractArray{FT,3}, r⁺⁻::AbstractArray{FT,3}, t⁻⁻::AbstractArray{FT,3}, + dr⁻⁺::AbstractArray{FT,4}, dt⁺⁺::AbstractArray{FT,4}, dr⁺⁻::AbstractArray{FT,4}, dt⁻⁻::AbstractArray{FT,4}) where {FT} + if n_stokes == 1 + r⁺⁻ .= r⁻⁺ + t⁻⁻ .= t⁺⁺ + dr⁺⁻ .= dr⁻⁺ + dt⁻⁻ .= dt⁺⁺ + return nothing + else + device = devi(architecture(r⁻⁺)) + applyD_kernel! = apply_D!(device) + applyD_kernel!(n_stokes, + r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, + dr⁻⁺, dt⁺⁺, dr⁺⁻, dt⁻⁻, + ndrange=size(r⁻⁺)); +# wait(device, event); + synchronize_if_gpu(); + return nothing + end +end + + +function apply_D_matrix_SFI!(n_stokes::Int, + J₀⁻::AbstractArray{FT,3}, + dJ₀⁻::AbstractArray{FT,4}) where {FT} + n_stokes == 1 && return nothing + device = devi(architecture(J₀⁻)) + applyD_kernel! = apply_D_SFI!(device) + applyD_kernel!(n_stokes, J₀⁻, dJ₀⁻, ndrange=size(J₀⁻)); + # wait(device, event); + synchronize_if_gpu(); + nothing +end diff --git a/src/CoreRT/CoreKernel/lin_elemental.jl b/src/CoreRT/CoreKernel/lin_elemental.jl new file mode 100644 index 00000000..f5bb9d7d --- /dev/null +++ b/src/CoreRT/CoreKernel/lin_elemental.jl @@ -0,0 +1,417 @@ +#= + +This file contains RT elemental-related functions + +=# +#= +"Elemental single-scattering layer" +function elemental!(pol_type, SFI::Bool, + τ_sum::AbstractArray, #{FT2,1}, #Suniti + dτ_λ::AbstractArray{FT,1}, # dτ_λ: total optical depth of elemental layer (per λ) + dτ::FT, # dτ: scattering optical depth of elemental layer (scalar) + ϖ_λ::AbstractArray{FT,1}, # ϖ_λ: single scattering albedo of elemental layer (per λ, absorptions by gases included) + ϖ::FT, # ϖ: single scattering albedo of elemental layer (no trace gas absorption included) + Z⁺⁺::AbstractArray{FT,2}, # Z matrix + Z⁻⁺::AbstractArray{FT,2}, # Z matrix + m::Int, # m: fourier moment + ndoubl::Int, # ndoubl: number of doubling computations needed + scatter::Bool, # scatter: flag indicating scattering + quad_points::QuadPoints{FT2}, # struct with quadrature points, weights, + added_layer::Union{AddedLayer{FT},AddedLayerRS{FT}}, + I_static, + architecture) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2} + + @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, J₀⁺, J₀⁻ = added_layer + @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart, iμ₀ = quad_points + #@unpack ϖ_Cabannes = RS_type + arr_type = array_type(architecture) + # Need to check with paper nomenclature. This is basically eqs. 19-20 in vSmartMOM + # @show Array(τ_sum)[1], Array(dτ_λ)[1], Array(ϖ_λ)[1], Array(Z⁺⁺)[1,1] + # Later on, we can have Zs also vary with index, pretty easy here: + # Z⁺⁺_ = repeat(Z⁺⁺, 1, 1, 1) + Z⁺⁺_ = reshape(Z⁺⁺, (size(Z⁺⁺,1), size(Z⁺⁺,2),1)) + # Z⁻⁺_ = repeat(Z⁻⁺, 1, 1, 1) + Z⁻⁺_ = reshape(Z⁻⁺, (size(Z⁺⁺,1), size(Z⁺⁺,2),1)) + + D = Diagonal(arr_type(repeat(pol_type.D, size(qp_μ,1)))) + I₀_NquadN = arr_type(zeros(FT,size(qp_μN,1))); #incident irradiation + i_end = pol_type.n*iμ₀ + I₀_NquadN[iμ₀Nstart:i_end] = pol_type.I₀ + + device = devi(architecture) + + # If in scattering mode: + if scatter + + NquadN = length(qp_μN) + + # Needs explanation still, different weights: + # for m==0, ₀∫²ᵖⁱ cos²(mϕ)dϕ/4π = 0.5, while + # for m>0, ₀∫²ᵖⁱ cos²(mϕ)dϕ/4π = 0.25 + wct0 = m == 0 ? FT(0.50) * ϖ * dτ : FT(0.25) * ϖ * dτ + wct02 = m == 0 ? FT(0.50) : FT(0.25) + wct = m == 0 ? FT(0.50) * ϖ * wt_μN : FT(0.25) * ϖ * wt_μN + wct2 = m == 0 ? wt_μN/2 : wt_μN/4 + + # Get the diagonal matrices first + d_qp = Diagonal(1 ./ qp_μN) + d_wct = Diagonal(wct) + + # Calculate r⁻⁺ and t⁺⁺ + + # Version 1: no absorption in batch mode (initiation of a single scattering layer with no or low absorption) + if false #maximum(dτ_λ) < 0.0001 + # R⁻⁺₀₁(λ) = M⁻¹(0.5ϖₑ(λ)Z⁻⁺C)δ (See Eqs.7 in Raman paper draft) + r⁻⁺[:,:,:] .= d_qp * Z⁻⁺ * (d_wct * dτ) + # T⁺⁺₀₁(λ) = {I-M⁻¹[I - 0.5*ϖₑ(λ)Z⁺⁺C]}δ (See Eqs.7 in Raman paper draft) + t⁺⁺[:,:,:] .= I_static - (d_qp * ((I_static - Z⁺⁺ * d_wct) * dτ)) + if SFI + # Reminder: Add equation here what it does + expk = exp.(-τ_sum/qp_μ[iμ₀]) #exp(-τ(z)/μ₀) + # J₀⁺ = 0.5[1+δ(m,0)]M⁻¹ϖₑ(λ)Z⁺⁺τI₀exp(-τ(z)/μ₀) + J₀⁺[:,1,:] .= ((d_qp * Z⁺⁺ * I₀_NquadN) * wct0) .* expk' + # J₀⁻ = 0.5[1+δ(m,0)]M⁻¹ϖₑ(λ)Z⁻⁺τI₀exp(-τ(z)/μ₀) + J₀⁻[:,1,:] .= ((d_qp * Z⁻⁺ * I₀_NquadN) * wct0) .* expk' + + end + else + # Version 2: More computationally intensive definition of a single scattering layer with variable (0-∞) absorption + # Version 2: with absorption in batch mode, low tau_scatt but higher tau_total, needs different equations + kernel! = get_elem_rt!(device) + event = kernel!(r⁻⁺, t⁺⁺, ϖ_λ, dτ_λ, Z⁻⁺, Z⁺⁺, + qp_μN, wct2, ndrange=size(r⁻⁺)); + wait(device, event) + synchronize_if_gpu() + + if SFI + kernel! = get_elem_rt_SFI!(device) + event = kernel!(J₀⁺, J₀⁻, ϖ_λ, dτ_λ, τ_sum, Z⁻⁺, Z⁺⁺, qp_μN, ndoubl, wct02, pol_type.n, arr_type(pol_type.I₀), iμ₀, D, ndrange=size(J₀⁺)) + wait(device, event) + synchronize_if_gpu() + end + end + + # Apply D Matrix + apply_D_matrix_elemental!(ndoubl, pol_type.n, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻) + + if SFI + apply_D_matrix_elemental_SFI!(ndoubl, pol_type.n, J₀⁻) + end + else + # Note: τ is not defined here + t⁺⁺[:] = Diagonal{exp(-τ ./ qp_μN)} + t⁻⁻[:] = Diagonal{exp(-τ ./ qp_μN)} + end + #@pack! added_layer = r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, J₀⁺, J₀⁻ +end +=# +"Elemental single-scattering layer" +function elemental!(pol_type, SFI::Bool, + τ_sum::AbstractArray,#{FT2,1}, #Suniti + lin_τ_sum::AbstractArray, + dτ::AbstractArray, + lin_dτ::AbstractArray, + computed_layer_properties::CoreScatteringOpticalProperties, + lin_computed_layer_properties::linCoreScatteringOpticalProperties, + m::Int, # m: fourier moment + ndoubl::Int, # ndoubl: number of doubling computations needed + scatter::Bool, # scatter: flag indicating scattering + quad_points::QuadPoints{FT2}, # struct with quadrature points, weights, + added_layer::Union{AddedLayer{FT},AddedLayerRS{FT}}, + lin_added_layer::Union{linAddedLayer{FT},linAddedLayerRS{FT}} + architecture) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2} + + @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, j₀⁺, j₀⁻ = added_layer + # the following contain core RT derivatives with respect to τ, ϖ, and Z specific to the current layer and band. + # Suniti: carry out linearization with respect to the state vector only outside rt_kernel + @unpack dr⁺⁻, dr⁻⁺, dt⁻⁻, dt⁺⁺, dj₀⁺, dj₀⁻ = lin_added_layer + @unpack qp_μ, iμ₀, wt_μN, qp_μN = quad_points + @unpack τ, ϖ, Z⁺⁺, Z⁻⁺ = computed_layer_properties + #@unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺ = computed_layer_properties + + + #@show M + arr_type = array_type(architecture) + + # Need to check with paper nomenclature. This is basically eqs. 19-20 in vSmartMOM + I₀ = arr_type(pol_type.I₀) + D = Diagonal(arr_type(repeat(pol_type.D, size(qp_μ,1)))) + + device = devi(architecture) + + # If in scattering mode: + if scatter + # for m==0, ₀∫²ᵖⁱ cos²(mϕ)dϕ/4π = 0.5, while + # for m>0, ₀∫²ᵖⁱ cos²(mϕ)dϕ/4π = 0.25 + wct02 = m == 0 ? FT(0.50) : FT(0.25) + wct2 = m == 0 ? wt_μN/2 : wt_μN/4 + + # More computationally intensive definition of a single scattering layer with variable (0-∞) absorption + # with absorption in batch mode, low tau_scatt but higher tau_total, needs exact equations + kernel! = get_elem_rt!(device) + event = kernel!(r⁻⁺, t⁺⁺, + dr⁻⁺, dt⁺⁺, + ϖ, dτ, Z⁻⁺, Z⁺⁺, + #lin_ϖ, lin_dτ, lin_Z⁻⁺, lin_Z⁺⁺, + qp_μN, wct2, ndrange=size(r⁻⁺)); +# wait(device, event) + synchronize_if_gpu() + + # SFI part + kernel! = get_elem_rt_SFI!(device) + event = kernel!(j₀⁺, j₀⁻, + dj₀⁺, dj₀⁻, + ϖ, dτ, arr_type(τ_sum), Z⁻⁺, Z⁺⁺, + qp_μN, ndoubl, wct02, + pol_type.n, I₀, iμ₀, D, ndrange=size(j₀⁺)) + # wait(device, event) + synchronize_if_gpu() + + # Apply D Matrix + apply_D_matrix_elemental!(ndoubl, pol_type.n, + r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, + dr⁻⁺, dt⁺⁺, dr⁺⁻, dt⁻⁻) + + # apply D matrix for SFI + apply_D_matrix_elemental_SFI!(ndoubl, pol_type.n, j₀⁻, , dj₀⁻) + else + # Note: τ is not defined here + t⁺⁺ .= Diagonal{exp(-τ ./ qp_μN)} + t⁻⁻ .= t⁺⁺ #Diagonal{exp(-τ ./ qp_μN)} + # Derivativve wrt τ + dt⁺⁺[1,:,:] .= Diagonal{-exp(-τ ./ qp_μN)./qp_μN} + dt⁻⁻[1,:,:] .= dt⁺⁺ + dt⁺⁺[2:3,:,:] .= 0 + dt⁻⁻[2:3,:,:] .= 0 + end +end + +@kernel function get_elem_rt!(r⁻⁺, t⁺⁺, dr⁻⁺, dt⁺⁺, ϖ_λ, dτ_λ, Z⁻⁺, Z⁺⁺, μ, wct) + n2 = 1 + i, j, n = @index(Global, NTuple) + if size(Z⁻⁺,3)>1 + n2 = n + end + if (wct[j]>1.e-8) + tmpM = exp(-dτ_λ[n] * ((1 / μ[i]) + (1 / μ[j]))) + + # 𝐑⁻⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁻⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ+μⱼ)) ̇(1 - exp{-τ ̇(1/μᵢ + 1/μⱼ)}) ̇𝑤ⱼ + r⁻⁺[i,j,n] = + ϖ_λ[n] * Z⁻⁺[i,j,n2] * + (μ[j] / (μ[i] + μ[j])) * wct[j] * + (1 - tmpM) + + dr⁻⁺[1,i,j,n] = + ϖ_λ[n] * Z⁻⁺[i,j,n2] * + (1 / μ[i]) * wct[j] * + tmpM + + dr⁻⁺[2, i,j,n] = + Z⁻⁺[i,j,n2] * + (μ[j] / (μ[i] + μ[j])) * wct[j] * + (1 - tmpM) + + dr⁻⁺[3, i,j,n] = + ϖ_λ[n] * + (μ[j] / (μ[i] + μ[j])) * wct[j] * + (1 - tmpM) + + if (μ[i] == μ[j]) + # 𝐓⁺⁺(μᵢ, μᵢ) = (exp{-τ/μᵢ} + ϖ ̇𝐙⁺⁺(μᵢ, μᵢ) ̇(τ/μᵢ) ̇exp{-τ/μᵢ}) ̇𝑤ᵢ + if i == j + tmpM = exp(-dτ_λ[n] / μ[i]) + t⁺⁺[i,j,n] = + tmpM * + (1 + ϖ_λ[n] * Z⁺⁺[i,i,n2] * (dτ_λ[n] / μ[i]) * wct[i]) + dt⁺⁺[1,i,j,n] = ( - t⁺⁺[i,j,n] + tmpM * ϖ_λ[n] * Z⁺⁺[i,i,n2] * wct[i]) / μ[i] + dt⁺⁺[2,i,j,n] = (tmpM * Z⁺⁺[i,i,n2] * dτ_λ[n] * wct[i]) / μ[i] + dt⁺⁺[3,i,j,n] = (tmpM * ϖ_λ[n] * dτ_λ[n] * wct[i]) / μ[i] + else + t⁺⁺[i,j,n] = 0.0 + dt⁺⁺[1:3, i,j,n] = 0.0 + end + else + + # 𝐓⁺⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁺⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ-μⱼ)) ̇(exp{-τ/μᵢ} - exp{-τ/μⱼ}) ̇𝑤ⱼ + # (𝑖 ≠ 𝑗) + tmpM = (exp(-dτ_λ[n] / μ[i]) - exp(-dτ_λ[n] / μ[j])) + t⁺⁺[i,j,n] = + ϖ_λ[n] * Z⁺⁺[i,j,n2] * + #Z⁺⁺[i,j] * + (μ[j] / (μ[i] - μ[j])) * wct[j] * + tmpM + dt⁺⁺[1, i,j,n] = + -(t⁺⁺[i,j,n]/tmpM) * + (exp(-dτ_λ[n] / μ[i])/ μ[i] - exp(-dτ_λ[n] / μ[j])/ μ[j]) + dt⁺⁺[2, i,j,n] = + t⁺⁺[i,j,n]/ϖ_λ[n] + dt⁺⁺[3, i,j,n] = + t⁺⁺[i,j,n]/Z⁺⁺[i,j,n2] + end + else + r⁻⁺[i,j,n] = 0.0 + dr⁻⁺[1:3,i,j,n] = 0.0 + if i==j + t⁺⁺[i,j,n] = exp(-dτ_λ[n] / μ[i]) #Suniti + dt⁺⁺[1,i,j,n] = -t⁺⁺[i,j,n] / μ[i] + dt⁺⁺[2:3,i,j,n] = 0.0 + else + t⁺⁺[i,j,n] = 0.0 + dt⁺⁺[1:3,i,j,n] = 0.0 + end + end + nothing +end + +@kernel function get_elem_rt_SFI!(J₀⁺, J₀⁻, dJ₀⁺, dJ₀⁻, ϖ_λ, dτ_λ, τ_sum, Z⁻⁺, Z⁺⁺, μ, ndoubl, wct02, nStokes ,I₀, iμ0, D) + i_start = nStokes*(iμ0-1) + 1 + i_end = nStokes*iμ0 + + i, _, n = @index(Global, NTuple) ##Suniti: What are Global and Ntuple? + FT = eltype(I₀) + J₀⁺[i, 1, n]=0 + J₀⁻[i, 1, n]=0 + dJ₀⁺[1:3, i, 1, n]=0 + dJ₀⁻[1:3, i, 1, n]=0 + n2=1 + if size(Z⁻⁺,3)>1 + n2 = n + end + + Z⁺⁺_I₀ = FT(0.0); + Z⁻⁺_I₀ = FT(0.0); + dZ⁺⁺_I₀ = FT(0.0); + dZ⁻⁺_I₀ = FT(0.0); + + for ii = i_start:i_end + Z⁺⁺_I₀ += Z⁺⁺[i,ii,n2] * I₀[ii-i_start+1] + Z⁻⁺_I₀ += Z⁻⁺[i,ii,n2] * I₀[ii-i_start+1] + if(ii==i) + dZ⁺⁺_I₀ += I₀[ii-i_start+1] + dZ⁻⁺_I₀ += I₀[ii-i_start+1] + end + end + + if (i>=i_start) && (i<=i_end) + ctr = i-i_start+1 + # See Eq. 1.54 in Fell + # J₀⁺ = 0.25*(1+δ(m,0)) * ϖ(λ) * Z⁺⁺ * I₀ * (dτ(λ)/μ₀) * exp(-dτ(λ)/μ₀) + J₀⁺[i, 1, n] = wct02 * ϖ_λ[n] * Z⁺⁺_I₀ + * (dτ_λ[n] / μ[i]) * exp(-dτ_λ[n] / μ[i]) + + dJ₀⁺[1, i, 1, n] = (J₀⁺[i, 1, n] / dτ_λ[n]) * (1 - dτ_λ[n] / μ[i]) + dJ₀⁺[2, i, 1, n] = (J₀⁺[i, 1, n] / ϖ_λ[n]) + dJ₀⁺[3, i, 1, n] = (J₀⁺[i, 1, n] / Z⁺⁺_I₀) * dZ⁺⁺_I₀ + else + # J₀⁺ = 0.25*(1+δ(m,0)) * ϖ(λ) * Z⁺⁺ * I₀ * [μ₀ / (μᵢ - μ₀)] * [exp(-dτ(λ)/μᵢ) - exp(-dτ(λ)/μ₀)] + # See Eq. 1.53 in Fell + tmpM = (exp(-dτ_λ[n] / μ[i]) - exp(-dτ_λ[n] / μ[i_start])) + J₀⁺[i, 1, n] = + wct02 * ϖ_λ[n] * Z⁺⁺_I₀ * + (μ[i_start] / (μ[i] - μ[i_start])) * + tmpM + dJ₀⁺[1, i, 1, n] = - (J₀⁺[i, 1, n] / tmpM) * + (exp(-dτ_λ[n] / μ[i])/μ[i] - exp(-dτ_λ[n] / μ[i_start])/ μ[i_start]) + dJ₀⁺[2, i, 1, n] = (J₀⁺[i, 1, n] / ϖ_λ[n]) + dJ₀⁺[3, i, 1, n] = (J₀⁺[i, 1, n] / Z⁺⁺_I₀) * dZ⁺⁺_I₀ + end + #J₀⁻ = 0.25*(1+δ(m,0)) * ϖ(λ) * Z⁻⁺ * I₀ * [μ₀ / (μᵢ + μ₀)] * [1 - exp{-dτ(λ)(1/μᵢ + 1/μ₀)}] + # See Eq. 1.52 in Fell + tmpM = (1 - exp(-dτ_λ[n] * ((1 / μ[i]) + (1 / μ[i_start])))) + J₀⁻[i, 1, n] = wct02 * ϖ_λ[n] * Z⁻⁺_I₀ * + (μ[i_start] / (μ[i] + μ[i_start])) * tmpM + dJ₀⁻[1, i, 1, n] = (J₀⁻[i, 1, n] / tmpM) * + exp(-dτ_λ[n] * ((1 / μ[i]) + (1 / μ[i_start]))) + * ((1 / μ[i]) + (1 / μ[i_start])) + dJ₀⁻[2, i, 1, n] = J₀⁻[i, 1, n] / ϖ_λ[n] + dJ₀⁻[3, i, 1, n] = (J₀⁻[i, 1, n] / Z⁻⁺_I₀) * dZ⁻⁺_I₀ + + J₀⁺[i, 1, n] *= exp(-τ_sum[n]/μ[i_start]) # how to do this?! Add a fourth derivative to RT kernel elements (only for J terms) + J₀⁻[i, 1, n] *= exp(-τ_sum[n]/μ[i_start]) # 1: wrt τ, 2: wrt ϖ, 3: wrt Z, 4: wrt τ_sum + J₀⁺[4, i, 1, n] = - J₀⁺[i, 1, n]/μ[i_start] + J₀⁻[4, i, 1, n] = - J₀⁻[i, 1, n]/μ[i_start] + if ndoubl >= 1 + J₀⁻[i, 1, n] = D[i,i]*J₀⁻[i, 1, n] #D = Diagonal{1,1,-1,-1,...Nquad times} + dJ₀⁻[1, i, 1, n] = D[i,i]*dJ₀⁻[1, i, 1, n] + dJ₀⁻[2, i, 1, n] = D[i,i]*dJ₀⁻[2, i, 1, n] + dJ₀⁻[3, i, 1, n] = D[i,i]*dJ₀⁻[3, i, 1, n] + dJ₀⁻[4, i, 1, n] = D[i,i]*dJ₀⁻[4, i, 1, n] + end + nothing +end + +@kernel function apply_D_elemental!(ndoubl, pol_n, + r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, + dr⁻⁺, dt⁺⁺, dr⁺⁻, dt⁻⁻) + i, j, n = @index(Global, NTuple) + + if ndoubl < 1 + ii = mod(i, pol_n) + jj = mod(j, pol_n) + if ((ii <= 2) & (jj <= 2)) | ((ii > 2) & (jj > 2)) + r⁺⁻[i,j,n] = r⁻⁺[i,j,n] + t⁻⁻[i,j,n] = t⁺⁺[i,j,n] + dr⁺⁻[1:4,i,j,n] .= dr⁻⁺[1:4,i,j,n] + dt⁻⁻[1:4,i,j,n] .= dt⁺⁺[1:4,i,j,n] + else + r⁺⁻[i,j,n] = -r⁻⁺[i,j,n] + t⁻⁻[i,j,n] = -t⁺⁺[i,j,n] + dr⁺⁻[1:4,i,j,n] = -dr⁻⁺[1:4,i,j,n] + dt⁻⁻[1:4,i,j,n] = -dt⁺⁺[1:4,i,j,n] + end + else + if mod(i, pol_n) > 2 + r⁻⁺[i,j,n] = - r⁻⁺[i,j,n] + dr⁻⁺[1:4,i,j,n] .= - dr⁻⁺[1:4,i,j,n] + end + end + nothing +end + +@kernel function apply_D_elemental_SFI!(ndoubl, pol_n, J₀⁻, dJ₀⁻) + i, _, n = @index(Global, NTuple) + + if ndoubl>1 + if mod(i, pol_n) > 2 + J₀⁻[i, 1, n] = - J₀⁻[i, 1, n] + dJ₀⁻[1:4, i, 1, n] .= - dJ₀⁻[1:4, i, 1, n] + end + end + nothing +end + +function apply_D_matrix_elemental!(ndoubl::Int, n_stokes::Int, + r⁻⁺::AbstractArray{FT,3}, + t⁺⁺::AbstractArray{FT,3}, + r⁺⁻::AbstractArray{FT,3}, + t⁻⁻::AbstractArray{FT,3}, + dr⁻⁺::AbstractArray{FT,4}, + dt⁺⁺::AbstractArray{FT,4}, + dr⁺⁻::AbstractArray{FT,4}, + dt⁻⁻::AbstractArray{FT,4} + ) where {FT} + device = devi(architecture(r⁻⁺)) + applyD_kernel! = apply_D_elemental!(device) + applyD_kernel!(ndoubl,n_stokes, + r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, + dr⁻⁺, dt⁺⁺, dr⁺⁻, dt⁻⁻, + ndrange=size(r⁻⁺)); +# wait(device, event); + synchronize_if_gpu(); + return nothing +end + +function apply_D_matrix_elemental_SFI!(ndoubl::Int, n_stokes::Int, + J₀⁻::AbstractArray{FT,3}, + dJ₀⁻::AbstractArray{FT,4}) where {FT} + if ndoubl > 1 + return nothing + else + device = devi(architecture(J₀⁻)) + applyD_kernel! = apply_D_elemental_SFI!(device) + applyD_kernel!(ndoubl,n_stokes, J₀⁻, dJ₀⁻, ndrange=size(J₀⁻)); + # wait(device, event); + synchronize_if_gpu(); + return nothing + end +end \ No newline at end of file diff --git a/src/CoreRT/CoreKernel/lin_interaction.jl b/src/CoreRT/CoreKernel/lin_interaction.jl new file mode 100644 index 00000000..615eaac4 --- /dev/null +++ b/src/CoreRT/CoreKernel/lin_interaction.jl @@ -0,0 +1,262 @@ +#= + +This file contains RT interaction-related functions + +=# + +# No scattering in either the added layer or the composite layer +function interaction_helper!(::ScatteringInterface_00, SFI, + composite_layer::CompositeLayer{FT}, + lin_composite_layer::linCompositeLayer{FT}, + added_layer::AddedLayer{FT}, + lin_added_layer::linAddedLayer{FT}, + nparams::Int, + I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2} + @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, j₀⁺, j₀⁻ = added_layer + @unpack R⁻⁺, R⁺⁻, T⁺⁺, T⁻⁻, J₀⁺, J₀⁻ = composite_layer + @unpack dxr⁺⁻, dxr⁻⁺, dxt⁻⁻, dxt⁺⁺, dxj₀⁺, dxj₀⁻ = lin_added_layer + @unpack dR⁻⁺, dR⁺⁻, dT⁺⁺, dT⁻⁻, dJ₀⁺, dJ₀⁻ = lin_composite_layer + + for ctr=1:nparams + # Source Function + dJ₀⁺[ctr,:,1,:] .= dxj₀⁺[ctr,:,1,:] .+ + t⁺⁺[:,:,:] ⊠ dJ₀⁺[ctr,:,1,:] .+ + dxt⁺⁺[ctr,:,:,:] ⊠ J₀⁺[:,1,:] + dJ₀⁻[ctr,:,1,:] .= dJ₀⁻[ctr,:,1,:] .+ + T⁻⁻[:,:,:] ⊠ dxj₀⁻[ctr,:,1,:] .+ + dT⁻⁻[ctr,:,:,:] ⊠ j₀⁻[:,1,:] + + # Batched multiplication between added and composite + dT⁻⁻[ctr,:,:,:] .= + dxt⁻⁻[ctr,:,:,:] ⊠ T⁻⁻[:,:,:] .+ + t⁻⁻[:,:,:] ⊠ dT⁻⁻[ctr,:,:,:] + dT⁺⁺[ctr,:,:,:] .= + dxt⁺⁺[ctr,:,:,:] ⊠ T⁺⁺[:,:,:] .+ + t⁺⁺[:,:,:] ⊠ dT⁺⁺[ctr,:,:,:] + end + # Source Function + J₀⁺ .= j₀⁺ .+ t⁺⁺ ⊠ J₀⁺ + J₀⁻ .= J₀⁻ .+ T⁻⁻ ⊠ j₀⁻ + + # Batched multiplication between added and composite + T⁻⁻ .= t⁻⁻ ⊠ T⁻⁻ + T⁺⁺ .= t⁺⁺ ⊠ T⁺⁺ +end + +# No scattering in inhomogeneous composite layer. +# Scattering in homogeneous layer, added to bottom of the composite layer. +# Produces a new, scattering composite layer. +function interaction_helper!(::ScatteringInterface_01, SFI, + composite_layer::CompositeLayer{FT}, + lin_composite_layer::linCompositeLayer{FT}, + added_layer::AddedLayer{FT}, + lin_added_layer::linAddedLayer{FT}, + nparams::Int + I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2} + @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, j₀⁺, j₀⁻ = added_layer + @unpack R⁻⁺, R⁺⁻, T⁺⁺, T⁻⁻, J₀⁺, J₀⁻ = composite_layer + @unpack dxr⁺⁻, dxr⁻⁺, dxt⁻⁻, dxt⁺⁺, dxj₀⁺, dxj₀⁻ = lin_added_layer + @unpack dR⁻⁺, dR⁺⁻, dT⁺⁺, dT⁻⁻, dJ₀⁺, dJ₀⁻ = lin_composite_layer + + for ctr=1:nparams + # Source Function + dJ₀⁻[ctr,:,1,:] .= dJ₀⁻[ctr,:,1,:] .+ + dT⁻⁻[ctr,:,:,:] ⊠ (r⁻⁺ ⊠ J₀⁺ .+ j₀⁻) .+ + T⁻⁻[ctr,:,:,:] ⊠ + (dxr⁻⁺[ctr,:,:,:] ⊠ J₀⁺ .+ + r⁻⁺ ⊠ dJ₀⁺[ctr,:,1,:] .+ + dxj₀⁻[ctr,:,1,:]) + + dJ₀⁺[ctr,:,1,:] .= dxj₀⁺[ctr,:,1,:] .+ + dxt⁺⁺[ctr,:,:,:] ⊠ J₀⁺ .+ + t⁺⁺ ⊠ dJ₀⁺[ctr,:,1,:] + + # Batched multiplication between added and composite + dR⁻⁺[ctr,:,:,:] .= dT⁻⁻[ctr,:,:,:] ⊠ r⁻⁺ ⊠ T⁺⁺ .+ + T⁻⁻ ⊠ dxr⁻⁺[ctr,:,:,:] ⊠ T⁺⁺ .+ + T⁻⁻ ⊠ r⁻⁺ ⊠ dT⁺⁺[ctr,:,:,:] + dR⁺⁻[ctr,:,:,:] .= dxr⁺⁻[ctr,:,:,:] + dT⁺⁺[ctr,:,:,:] .= dxt⁺⁺[ctr,:,:,:] ⊠ T⁺⁺ .+ + t⁺⁺ ⊠ dT⁺⁺[ctr,:,:,:] + dT⁻⁻[ctr,:,:,:] .= dT⁻⁻[ctr,:,:,:] ⊠ t⁻⁻ .+ + T⁻⁻ ⊠ dxt⁻⁻[ctr,:,:,:] + end + # Source Function + J₀⁻ .= J₀⁻ .+ T⁻⁻ ⊠ (r⁻⁺ ⊠ J₀⁺ .+ j₀⁻) + J₀⁺ .= j₀⁺ .+ t⁺⁺ ⊠ J₀⁺ + + # Batched multiplication between added and composite + R⁻⁺ .= T⁻⁻ ⊠ r⁻⁺ ⊠ T⁺⁺ + R⁺⁻ .= r⁺⁻ + T⁺⁺ .= t⁺⁺ ⊠ T⁺⁺ + T⁻⁻ .= T⁻⁻ ⊠ t⁻⁻ +end + +# Scattering in inhomogeneous composite layer. +# no scattering in homogeneous layer which is +# added to the bottom of the composite layer. +# Produces a new, scattering composite layer. +function interaction_helper!(::ScatteringInterface_10, SFI, + composite_layer::CompositeLayer{FT}, + lin_composite_layer::linCompositeLayer{FT}, + added_layer::AddedLayer{FT}, + lin_added_layer::linAddedLayer{FT}, + nparams::Int, + I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2} + @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, j₀⁺, j₀⁻ = added_layer + @unpack R⁻⁺, R⁺⁻, T⁺⁺, T⁻⁻, J₀⁺, J₀⁻ = composite_layer + @unpack dxr⁺⁻, dxr⁻⁺, dxt⁻⁻, dxt⁺⁺, dxj₀⁺, dxj₀⁻ = lin_added_layer + @unpack dR⁻⁺, dR⁺⁻, dT⁺⁺, dT⁻⁻, dJ₀⁺, dJ₀⁻ = lin_composite_layer + + for ctr=1:nparams + # Source Function + dJ₀⁺[ctr,:,1,:] .= dxj₀⁺[ctr,:,1,:] .+ + dxt⁺⁺[ctr,:,:,:] ⊠ (J₀⁺ .+ R⁺⁻ ⊠ j₀⁻) .+ + t⁺⁺ ⊠ (dJ₀⁺[ctr,:,1,:] .+ dR⁺⁻[ctr,:,:,:] ⊠ j₀⁻ .+ R⁺⁻ ⊠ dxj₀⁻[ctr,:,1,:]) + dJ₀⁻[ctr,:,1,:] .= dJ₀⁻[ctr,:,1,:] .+ + dT⁻⁻[ctr,:,:,:] ⊠ j₀⁻ .+ T⁻⁻ ⊠ dxj₀⁻[ctr,:,1,:] + + # Batched multiplication between added and composite + dT⁺⁺[ctr,:,:,:] .= dxt⁺⁺[ctr,:,:,:] ⊠ T⁺⁺ .+ t⁺⁺ ⊠ dT⁺⁺[ctr,:,:,:] + dT⁻⁻[ctr,:,:,:] .= dT⁻⁻[ctr,:,:,:] ⊠ t⁻⁻ .+ T⁻⁻ ⊠ dxt⁻⁻[ctr,:,:,:] + dR⁺⁻[ctr,:,:,:] .= dxt⁺⁺[ctr,:,:,:] ⊠ R⁺⁻ ⊠ t⁻⁻ .+ + t⁺⁺ ⊠ (dR⁺⁻[ctr,:,:,:] ⊠ t⁻⁻ .+ R⁺⁻ ⊠ dxt⁻⁻[ctr,:,:,:]) + end + # Source Function + J₀⁺ .= j₀⁺ .+ t⁺⁺ ⊠ (J₀⁺ .+ R⁺⁻ ⊠ j₀⁻) + J₀⁻ .= J₀⁻ .+ T⁻⁻ ⊠ j₀⁻ + + # Batched multiplication between added and composite + T⁺⁺ .= t⁺⁺ ⊠ T⁺⁺ + T⁻⁻ .= T⁻⁻ ⊠ t⁻⁻ + R⁺⁻ .= t⁺⁺ ⊠ R⁺⁻ ⊠ t⁻⁻ +end + +# Scattering in inhomogeneous composite layer. +# Scattering in homogeneous layer which is added to the bottom of the composite layer. +# Produces a new, scattering composite layer. +function interaction_helper!(::ScatteringInterface_11, SFI, + composite_layer::CompositeLayer{FT}, + lin_composite_layer::linCompositeLayer{FT}, + added_layer::AddedLayer{FT}, + lin_added_layer::linAddedLayer{FT}, + nparams::Int, + I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2} + + @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, j₀⁺, j₀⁻ = added_layer #these are aliases to the respective struct elements + @unpack R⁻⁺, R⁺⁻, T⁺⁺, T⁻⁻, J₀⁺, J₀⁻ = composite_layer #these are aliases to the respective struct elements + @unpack dxr⁺⁻, dxr⁻⁺, dxt⁻⁻, dxt⁺⁺, dxj₀⁺, dxj₀⁻ = lin_added_layer + @unpack dR⁻⁺, dR⁺⁻, dT⁺⁺, dT⁻⁻, dJ₀⁺, dJ₀⁻ = lin_composite_layer + + # X₂₁ refers to added layer, X₁₀ to composite layer! + + # Used to store `(I - R⁺⁻ * r⁻⁺)⁻¹` + tmp_inv = similar(t⁺⁺) + + # Compute and store `(I - R⁺⁻ * r⁻⁺)⁻¹` + @timeit "interaction inv1" batch_inv!(tmp_inv, I_static .- r⁻⁺ ⊠ R⁺⁻) + # Temporary arrays: + + # T₁₂(I-R₀₁R₂₁)⁻¹ + T01_inv = T⁻⁻ ⊠ tmp_inv; + + dtmp_inv = zeros(nparams,size(t⁺⁺,1), ,size(t⁺⁺,2), ,size(t⁺⁺,3)) + dT01_inv = similar(dtmp_inv) + for ctr=1:nparams + dtmp_inv[ctr,:,:,:] = tmp_inv ⊠ + (dxr⁻⁺[ctr,:,:,:] ⊠ R⁺⁻ .+ + r⁻⁺ ⊠ dR⁺⁻[ctr,:,:,:]) ⊠ tmp_inv + dT01_inv[ctr,:,:,:] = dT⁻⁻[ctr,:,:,:] ⊠ tmp_inv .+ + T⁻⁻ ⊠ dtmp_inv[ctr,:,:,:] + + # J₀₂⁻ = J₀₁⁻ + T₀₁(1-R₂₁R₀₁)⁻¹(R₂₁J₁₀⁺+J₁₂⁻) + dJ₀⁻[ctr,:,1,:] .= dJ₀⁻[ctr,:,1,:] .+ + dT01_inv[ctr,:,:,:] ⊠ (r⁻⁺ ⊠ J₀⁺ .+ j₀⁻) .+ + T01_inv ⊠ (dxr⁻⁺[ctr,:,:,:] ⊠ J₀⁺ .+ + r⁻⁺ ⊠ dJ₀⁺[ctr,:,1,:] .+ + dxj₀⁻[ctr,:,1,:]) + + # R₂₀ = R₁₀ + T₀₁(I-R₂₁R₀₁)⁻¹ R₂₁T₁₀ + dR⁻⁺[ctr,:,:,:] .= dR⁻⁺[ctr,:,:,:] .+ + dT01_inv[ctr,:,:,:] ⊠ r⁻⁺ ⊠ T⁺⁺ .+ + T01_inv ⊠ (dxr⁻⁺[ctr,:,:,:] ⊠ T⁺⁺ .+ + r⁻⁺ ⊠ dT⁺⁺[ctr,:,:,:]) + + # T₀₂ = T₀₁(1-R₂₁R₀₁)⁻¹T₁₂ + dT⁻⁻[ctr,:,:,:] .= dT01_inv[ctr,:,:,:] ⊠ t⁻⁻ .+ + T01_inv ⊠ dxt⁻⁻[ctr,:,:,:] + end + + # J₀₂⁻ = J₀₁⁻ + T₀₁(1-R₂₁R₀₁)⁻¹(R₂₁J₁₀⁺+J₁₂⁻) + J₀⁻ .= J₀⁻ .+ T01_inv ⊠ (r⁻⁺ ⊠ J₀⁺ .+ j₀⁻) + + # R₂₀ = R₁₀ + T₀₁(I-R₂₁R₀₁)⁻¹ R₂₁T₁₀ + R⁻⁺ .= R⁻⁺ .+ T01_inv ⊠ r⁻⁺ ⊠ T⁺⁺ + + # T₀₂ = T₀₁(1-R₂₁R₀₁)⁻¹T₁₂ + T⁻⁻ .= T01_inv ⊠ t⁻⁻ + + # Repeating for mirror-reflected directions + + # Compute and store `(I - r⁻⁺ * R⁺⁻)⁻¹` + @timeit "interaction inv2" batch_inv!(tmp_inv, I_static .- R⁺⁻ ⊠ r⁻⁺) + # T₂₁(I-R₀₁R₂₁)⁻¹ + T21_inv = t⁺⁺ ⊠ tmp_inv + + dtmp_inv .= 0.0 + dT21_inv .= similar(dtmp_inv) + + for ctr=1:nparams + dtmp_inv[ctr,:,:,:] = tmp_inv ⊠ + (dR⁺⁻[ctr,:,:,:] ⊠ r⁻⁺ .+ + R⁺⁻ ⊠ dxr⁻⁺[ctr,:,:,:]) ⊠ tmp_inv + dT01_inv[ctr,:,:,:] = dxt⁺⁺[ctr,:,:,:] ⊠ tmp_inv .+ + t⁺⁺ ⊠ dtmp_inv[ctr,:,:,:] + + # J₂₀⁺ = J₂₁⁺ + T₂₁(I-R₀₁R₂₁)⁻¹(J₁₀ + R₀₁J₁₂⁻ ) + dJ₀⁺[ctr,:,1,:] .= dxj₀⁺[ctr,:,1,:] .+ + dT21_inv[ctr,:,:,:] ⊠ (R⁺⁻ ⊠ j₀⁻ .+ J₀⁺) .+ + T21_inv ⊠ (dR⁺⁻[ctr,:,:,:] ⊠ j₀⁻ .+ + R⁺⁻ ⊠ dxj₀⁻[ctr,:,1,:] .+ + dJ₀⁺[ctr,:,1,:]) + + # R₂₀ = R₁₀ + T₀₁(I-R₂₁R₀₁)⁻¹ R₂₁T₁₀ + dR⁻⁺[ctr,:,:,:] .= dxr⁺⁻[ctr,:,:,:] .+ + dT21_inv[ctr,:,:,:] ⊠ R⁺⁻ ⊠ t⁻⁻ .+ + T21_inv ⊠ (dR⁺⁻[ctr,:,:,:] ⊠ t⁻⁻ .+ + R⁺⁻ ⊠ dxt⁻⁻[ctr,:,:,:]) + + # T₀₂ = T₀₁(1-R₂₁R₀₁)⁻¹T₁₂ + dT⁺⁺[ctr,:,:,:] .= dT21_inv[ctr,:,:,:] ⊠ T⁺⁺ .+ + T21_inv ⊠ dT⁺⁺[ctr,:,:,:] + end + + + # J₂₀⁺ = J₂₁⁺ + T₂₁(I-R₀₁R₂₁)⁻¹(J₁₀ + R₀₁J₁₂⁻ ) + J₀⁺ .= j₀⁺ .+ T21_inv ⊠ (J₀⁺ .+ R⁺⁻ ⊠ j₀⁻) + + # T₂₀ = T₂₁(I-R₀₁R₂₁)⁻¹T₁₀ + T⁺⁺ .= T21_inv ⊠ T⁺⁺ + + # R₀₂ = R₁₂ + T₂₁(1-R₀₁R₂₁)⁻¹R₀₁T₁₂ + R⁺⁻ .= r⁺⁻ .+ T21_inv ⊠ R⁺⁻ ⊠ t⁻⁻ +end + +"Compute interaction between composite and added layers" +function interaction!(scattering_interface::AbstractScatteringInterface, SFI, +composite_layer::CompositeLayer{FT}, +lin_composite_layer::linCompositeLayer{FT}, +added_layer::AddedLayer{FT}, +lin_added_layer::linAddedLayer{FT}, +nparams::Int, +I_static::AbstractArray{FT2}) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2} + + interaction_helper!(scattering_interface, SFI, + composite_layer, + lin_composite_layer, + added_layer, + lin_added_layer, + nparams, + I_static) + synchronize_if_gpu() +end diff --git a/src/CoreRT/CoreKernel/lin_rt_kernel.jl b/src/CoreRT/CoreKernel/lin_rt_kernel.jl new file mode 100644 index 00000000..a59a7bd8 --- /dev/null +++ b/src/CoreRT/CoreKernel/lin_rt_kernel.jl @@ -0,0 +1,413 @@ +#= + +This file implements rt_kernel!, which performs the core RT routines (elemental, doubling, interaction) + +=# +#No Raman (default) +# Perform the Core RT routines (elemental, doubling, interaction) + +### +function rt_kernel!(RS_type::noRS{FT}, + pol_type, SFI, + added_layer, lin_added_layer, + composite_layer, lin_composite_layer, + computed_layer_properties::M, + lin_computed_layer_properties, + scattering_interface, + τ_sum, + lin_τ_sum, + m, quad_points, + I_static, + architecture, + qp_μN, iz) where {FT,M} + #@show array_type(architecture) + + @unpack qp_μ, μ₀ = quad_points + # Just unpack core optical properties from + @unpack τ, ϖ, Z⁺⁺, Z⁻⁺ = computed_layer_properties + @unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺ = lin_computed_layer_properties + + nparams = size(lin_τ,1) + # @show ndoubl + scatter = true # edit later + + dτ, lin_dτ, ndoubl, expk, dexpk_fctr = init_layer( + computed_layer_properties, + lin_computed_layer_properties, + quad_points, pol_type, architecture) + + # If there is scattering, perform the elemental and doubling steps + if scatter + #@show typeof(computed_layer_properties) + @timeit "elemental" elemental!(pol_type, SFI, + τ_sum, #lin_τ_sum, + dτ, #lin_dτ, + computed_layer_properties, + #lin_computed_layer_properties, + m, ndoubl, scatter, quad_points, + added_layer, + lin_added_layer, + architecture) + #println("Elemental done...") + @timeit "doubling" doubling!(pol_type, SFI, + expk, #lin_expk, + dexpk_fctr, + ndoubl, + added_layer, lin_added_layer, + I_static, architecture) + #@show added_layer.r⁻⁺[1:2,1,1], added_layer.r⁺⁻[1:2,1,1],added_layer.t⁺⁺[1:2,1,1], added_layer.t⁻⁻[1:2,1,1] + #@show dτ, ndoubl, expk + #println("Doubling done...") + else # This might not work yet on GPU! + # If not, there is no reflectance. Assign r/t appropriately + added_layer.r⁻⁺[:] .= 0; + added_layer.r⁺⁻[:] .= 0; + added_layer.j₀⁻[:] .= 0; + lin_added_layer.dr⁻⁺[:] .= 0; + lin_added_layer.dr⁺⁻[:] .= 0; + lin_added_layer.dj₀⁻[:] .= 0; + temp = Array(exp.(-τ_λ./qp_μN')) + dtemp = Array(-exp.(-τ_λ./qp_μN')./qp_μN') + #added_layer.t⁺⁺, added_layer.t⁻⁻ = (Diagonal(exp(-τ_λ / qp_μN)), Diagonal(exp(-τ_λ / qp_μN))) + for iλ = 1:length(τ_λ) + added_layer.t⁺⁺[:,:,iλ] = Diagonal(temp[iλ,:]); + added_layer.t⁻⁻[:,:,iλ] = Diagonal(temp[iλ,:]); + lin_added_layer.dt⁺⁺[1,:,:,iλ] = Diagonal(dtemp[iλ,:]); + lin_added_layer.dt⁻⁻[1,:,:,iλ] = Diagonal(dtemp[iλ,:]); + end + end + + # @assert !any(isnan.(added_layer.t⁺⁺)) + + # If this TOA, just copy the added layer into the composite layer + if (iz == 1) + + composite_layer.T⁺⁺[:], composite_layer.T⁻⁻[:] = (added_layer.t⁺⁺, added_layer.t⁻⁻) + composite_layer.R⁻⁺[:], composite_layer.R⁺⁻[:] = (added_layer.r⁻⁺, added_layer.r⁺⁻) + composite_layer.J₀⁺[:], composite_layer.J₀⁻[:] = (added_layer.j₀⁺, added_layer.j₀⁻) + + for ctr=1:nparams + lin_composite_layer.dT⁺⁺[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dt⁺⁺[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dt⁺⁺[2,:,:,:] + + lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dt⁺⁺[3,:,:,:] + + + lin_composite_layer.dT⁻⁻[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dt⁻⁻[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dt⁻⁻[2,:,:,:] + + lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dt⁻⁻[3,:,:,:] + + + lin_composite_layer.dR⁻⁺[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dr⁻⁺[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dr⁻⁺[2,:,:,:] + + lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dr⁻⁺[3,:,:,:] + + lin_composite_layer.dR⁺⁻[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dr⁺⁻[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dr⁺⁻[2,:,:,:] + + lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dr⁺⁻[3,:,:,:] + + lin_composite_layer.dJ₀⁺[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁺[1,:,1,:] + + lin_ϖ[ctr]*lin_added_layer.dj₀⁺[2,:,1,:] + + lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dj₀⁺[3,:,1,:] + + lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁺[4,:,1,:] + + lin_composite_layer.dJ₀⁻[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁻[1,:,1,:] + + lin_ϖ[ctr]*lin_added_layer.dj₀⁻[2,:,1,:] + + lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dj₀⁻[3,:,1,:] + + lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁻[4,:,1,:] + end + # If this is not the TOA, perform the interaction step + else + for ctr=1:nparams + lin_added_layer.dxt⁺⁺[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dt⁺⁺[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dt⁺⁺[2,:,:,:] + + lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dt⁺⁺[3,:,:,:] + + + lin_added_layer.dxt⁻⁻[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dt⁻⁻[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dt⁻⁻[2,:,:,:] + + lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dt⁻⁻[3,:,:,:] + + + lin_added_layer.dxr⁻⁺[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dr⁻⁺[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dr⁻⁺[2,:,:,:] + + lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dr⁻⁺[3,:,:,:] + + lin_added_layer.dxr⁺⁻[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dr⁺⁻[1,:,:,:] + + lin_ϖ[ctr]*lin_added_layer.dr⁺⁻[2,:,:,:] + + lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dr⁺⁻[3,:,:,:] + + lin_added_layer.dxj₀⁺[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁺[1,:,1,:] + + lin_ϖ[ctr]*lin_added_layer.dj₀⁺[2,:,1,:] + + lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dj₀⁺[3,:,1,:] + + lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁺[4,:,1,:] + + lin_added_layer.dxj₀⁻[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁻[1,:,1,:] + + lin_ϖ[ctr]*lin_added_layer.dj₀⁻[2,:,1,:] + + lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dj₀⁻[3,:,1,:] + + lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁻[4,:,1,:] + end + @timeit "interaction" interaction!(RS_type, scattering_interface, SFI, composite_layer, lin_composite_layer, added_layer, lin_added_layer, I_static) + end +end + +#= +function get_dtau_ndoubl(computed_layer_properties::CoreScatteringOpticalProperties, + lin_computed_layer_properties::linCoreScatteringOpticalProperties, + quad_points::QuadPoints{FT}) where {FT} + @unpack qp_μ = quad_points + @unpack τ, ϖ = computed_layer_properties + @unpack lin_τ = lin_computed_layer_properties + dτ_max = minimum([maximum(τ .* ϖ), FT(0.001) * minimum(qp_μ)]) + _, ndoubl = doubling_number(dτ_max, maximum(τ .* ϖ)) + # Compute dτ vector + dτ = τ ./ 2^ndoubl + lin_dτ = lin_τ / 2^ndoubl + + return dτ, lin_dτ, ndoubl +end + +function get_dtau_ndoubl(computed_layer_properties::CoreDirectionalScatteringOpticalProperties, + lin_computed_layer_properties::linCoreDirectionalScatteringOpticalProperties, + quad_points::QuadPoints{FT}) where {FT} + @unpack qp_μ,iμ₀ = quad_points + @unpack τ, ϖ, G = computed_layer_properties + @unpack lin_τ = lin_computed_layer_properties + gfct = Array(G)[iμ₀] + dτ_max = minimum([maximum(gfct * τ .* ϖ), FT(0.001) * minimum(qp_μ)]) + _, ndoubl = doubling_number(dτ_max, maximum(τ .* ϖ)) + # Compute dτ vector + dτ = τ ./ 2^ndoubl + lin_dτ = lin_τ / 2^ndoubl + + return dτ, lin_dτ, ndoubl +end + +function init_layer(computed_layer_properties::CoreDirectionalScatteringOpticalProperties, + lin_computed_layer_properties::linCoreDirectionalScatteringOpticalProperties, + quad_points, pol_type, architecture) + arr_type = array_type(architecture) + @unpack μ₀, iμ₀ = quad_points + @unpack G = computed_layer_properties + dτ, lin_dτ, ndoubl = get_dtau_ndoubl( + computed_layer_properties, + lin_computed_layer_properties, + quad_points) + gfct = Array(G)[iμ₀] + expk = exp.(-dτ*gfct/μ₀) + #lin_expk = -expk.*lin_dτ*gfct/μ₀ + dexpk_fctr = -gfct/μ₀ + + return dτ, lin_dτ, ndoubl, arr_type(expk), FT(dexpk_fctr) #arr_type(lin_expk) +end + +function init_layer(computed_layer_properties::CoreScatteringOpticalProperties, + lin_computed_layer_properties::linCoreScatteringOpticalProperties, + quad_points, pol_type, architecture) + arr_type = array_type(architecture) + @unpack μ₀ = quad_points + dτ, lin_dτ, ndoubl = get_dtau_ndoubl( + computed_layer_properties, + lin_computed_layer_properties, + quad_points) + expk = exp.(-dτ/μ₀) + #lin_expk = -expk.*lin_dτ/μ₀ + dexpk_fctr = -1/μ₀ + + return dτ, lin_dτ, ndoubl, arr_type(expk), FT(dexpk_fctr)#arr_type(lin_expk) +end + +=# +function rt_kernel!(RS_type::Union{RRS{FT}, VS_0to1{FT}, VS_1to0{FT}}, + pol_type, SFI, + added_layer, lin_added_layer, + composite_layer, lin_composite_layer, + computed_layer_properties::CoreScatteringOpticalProperties, + lin_computed_layer_properties::linCoreScatteringOpticalProperties, + scattering_interface, + τ_sum, + lin_τ_sum, + m, quad_points, + I_static, architecture, qp_μN, iz) where {FT} + @unpack qp_μ, μ₀ = quad_points + # Just unpack core optical properties from + @unpack τ, ϖ, Z⁺⁺, Z⁻⁺ = computed_layer_properties + @unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺ = lin_computed_layer_properties + + # SUNITI, check? Also, better to write function here + dτ_max = minimum([maximum(τ .* ϖ), FT(0.001) * minimum(qp_μ)]) + _, ndoubl = doubling_number(dτ_max, maximum(τ .* ϖ)) + scatter = true # edit later + arr_type = array_type(architecture) + # Compute dτ vector + dτ = τ ./ 2^ndoubl + expk = arr_type(exp.(-dτ /μ₀)) + + @unpack Z⁺⁺_λ₁λ₀, Z⁻⁺_λ₁λ₀ = RS_type + # If there is scattering, perform the elemental and doubling steps + if scatter + #@show τ, ϖ, RS_type.fscattRayl + + @timeit "elemental_inelastic" elemental_inelastic!(RS_type, + pol_type, SFI, + τ_sum, lin_τ_sum, + dτ, ϖ, + lin_dτ, lin_ϖ, + Z⁺⁺_λ₁λ₀, Z⁻⁺_λ₁λ₀, + m, ndoubl, scatter, + quad_points, + added_layer, lin_added_layer, + I_static, architecture) + #println("Elemental inelastic done...") + @timeit "elemental" elemental!(pol_type, SFI, + τ_sum, lin_τ_sum, + dτ, lin_dτ, + computed_layer_properties, + lin_computed_layer_properties, + m, ndoubl, scatter, quad_points, + added_layer, lin_added_layer, + architecture) + #println("Elemental done...") + @timeit "doubling_inelastic" doubling_inelastic!(RS_type, + pol_type, SFI, + expk, lin_expk, + ndoubl, + added_layer, lin_added_layer, + I_static, architecture) + #println("Doubling done...") + #@timeit "doubling" doubling!(pol_type, SFI, expk, ndoubl, added_layer, I_static, architecture) + else # This might not work yet on GPU! + # If not, there is no reflectance. Assign r/t appropriately + added_layer.r⁻⁺[:] .= 0; + added_layer.r⁺⁻[:] .= 0; + added_layer.j₀⁻[:] .= 0; + added_layer.ier⁻⁺[:] .= 0; + added_layer.ier⁺⁻[:] .= 0; + added_layer.ieJ₀⁻[:] .= 0; + added_layer.iet⁻⁻[:] .= 0; + added_layer.iet⁺⁺[:] .= 0; + added_layer.ieJ₀⁺[:] .= 0; + temp = Array(exp.(-τ_λ./qp_μN')) + #added_layer.t⁺⁺, added_layer.t⁻⁻ = (Diagonal(exp(-τ_λ / qp_μN)), Diagonal(exp(-τ_λ / qp_μN))) + for iλ = 1:length(τ_λ) + added_layer.t⁺⁺[:,:,iλ] = Diagonal(temp[iλ,:]); + added_layer.t⁻⁻[:,:,iλ] = Diagonal(temp[iλ,:]); + end + end + + # @assert !any(isnan.(added_layer.t⁺⁺)) + + # If this TOA, just copy the added layer into the composite layer + if (iz == 1) + composite_layer.T⁺⁺[:], composite_layer.T⁻⁻[:] = (added_layer.t⁺⁺, added_layer.t⁻⁻) + composite_layer.R⁻⁺[:], composite_layer.R⁺⁻[:] = (added_layer.r⁻⁺, added_layer.r⁺⁻) + composite_layer.J₀⁺[:], composite_layer.J₀⁻[:] = (added_layer.j₀⁺, added_layer.j₀⁻ ) + composite_layer.ieT⁺⁺[:], composite_layer.ieT⁻⁻[:] = (added_layer.iet⁺⁺, added_layer.iet⁻⁻) + composite_layer.ieR⁻⁺[:], composite_layer.ieR⁺⁻[:] = (added_layer.ier⁻⁺, added_layer.ier⁺⁻) + composite_layer.ieJ₀⁺[:], composite_layer.ieJ₀⁻[:] = (added_layer.ieJ₀⁺, added_layer.ieJ₀⁻ ) + + # If this is not the TOA, perform the interaction step + else + @timeit "interaction" interaction!(RS_type, scattering_interface, SFI, + composite_layer, lin_composite_layer, + added_layer, lin_added_layer, + I_static) + end +end + + + +function rt_kernel!( + RS_type::Union{RRS_plus{FT}, VS_0to1_plus{FT}, VS_1to0_plus{FT}}, + pol_type, SFI, + added_layer, + composite_layer, + computed_layer_properties::CoreScatteringOpticalProperties, + lin_computed_layer_properties::linCoreScatteringOpticalProperties, + scattering_interface, + τ_sum, + lin_τ_sum, + m, quad_points, + I_static, architecture, qp_μN, iz) where {FT} + @unpack qp_μ, μ₀ = quad_points + # Just unpack core optical properties from + @unpack τ, ϖ, Z⁺⁺, Z⁻⁺ = computed_layer_properties + @unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺ = lin_computed_layer_properties + + # SUNITI, check? Also, better to write function here + dτ_max = minimum([maximum(τ .* ϖ), FT(0.001) * minimum(qp_μ)]) + _, ndoubl = doubling_number(dτ_max, maximum(τ .* ϖ)) + scatter = true # edit later + arr_type = array_type(architecture) + # Compute dτ vector + dτ = τ ./ 2^ndoubl + expk = arr_type(exp.(-dτ /μ₀)) + + @unpack Z⁺⁺_λ₁λ₀, Z⁻⁺_λ₁λ₀ = RS_type + # If there is scattering, perform the elemental and doubling steps + if scatter + #@show τ, ϖ, RS_type.fscattRayl + + @timeit "elemental_inelastic" elemental_inelastic!(RS_type, + pol_type, SFI, + τ_sum, lin_τ_sum, + dτ, ϖ, + lin_dτ, lin_ϖ, + Z⁺⁺_λ₁λ₀, Z⁻⁺_λ₁λ₀, + m, ndoubl, scatter, + quad_points, + added_layer, lin_added_layer, + I_static, architecture) + #println("Elemental inelastic done...") + @timeit "elemental" elemental!(pol_type, SFI, + τ_sum, lin_τ_sum, + dτ, lin_dτ, + computed_layer_properties, + lin_computed_layer_properties, + m, ndoubl, scatter, + quad_points, + added_layer, lin_added_layer, + architecture) + #println("Elemental done...") + @timeit "doubling_inelastic" doubling_inelastic!(RS_type, pol_type, + SFI, + expk, lin_expk, + ndoubl, + added_layer, lin_added_layer, + I_static, architecture) + #println("Doubling done...") + #@timeit "doubling" doubling!(pol_type, SFI, expk, ndoubl, added_layer, I_static, architecture) + else # This might not work yet on GPU! + # If not, there is no reflectance. Assign r/t appropriately + added_layer.r⁻⁺[:] .= 0; + added_layer.r⁺⁻[:] .= 0; + added_layer.j₀⁻[:] .= 0; + added_layer.ier⁻⁺[:] .= 0; + added_layer.ier⁺⁻[:] .= 0; + added_layer.ieJ₀⁻[:] .= 0; + added_layer.iet⁻⁻[:] .= 0; + added_layer.iet⁺⁺[:] .= 0; + added_layer.ieJ₀⁺[:] .= 0; + temp = Array(exp.(-τ./qp_μN')) + #added_layer.t⁺⁺, added_layer.t⁻⁻ = (Diagonal(exp(-τ_λ / qp_μN)), Diagonal(exp(-τ_λ / qp_μN))) + for iλ = 1:length(τ) + added_layer.t⁺⁺[:,:,iλ] = Diagonal(temp[iλ,:]); + added_layer.t⁻⁻[:,:,iλ] = Diagonal(temp[iλ,:]); + end + end + + # @assert !any(isnan.(added_layer.t⁺⁺)) + + # If this TOA, just copy the added layer into the composite layer + if (iz == 1) + composite_layer.T⁺⁺[:], composite_layer.T⁻⁻[:] = (added_layer.t⁺⁺, added_layer.t⁻⁻) + composite_layer.R⁻⁺[:], composite_layer.R⁺⁻[:] = (added_layer.r⁻⁺, added_layer.r⁺⁻) + composite_layer.J₀⁺[:], composite_layer.J₀⁻[:] = (added_layer.j₀⁺, added_layer.j₀⁻ ) + composite_layer.ieT⁺⁺[:], composite_layer.ieT⁻⁻[:] = (added_layer.iet⁺⁺, added_layer.iet⁻⁻) + composite_layer.ieR⁻⁺[:], composite_layer.ieR⁺⁻[:] = (added_layer.ier⁻⁺, added_layer.ier⁺⁻) + composite_layer.ieJ₀⁺[:], composite_layer.ieJ₀⁻[:] = (added_layer.ieJ₀⁺, added_layer.ieJ₀⁻ ) + + # If this is not the TOA, perform the interaction step + else + @timeit "interaction" interaction!(RS_type, scattering_interface, SFI, + composite_layer, lin_composite_layer, + added_layer, lin_added_layer, I_static) + end +end + + diff --git a/src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl b/src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl index 1c5d633a..6d0d2c28 100644 --- a/src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl +++ b/src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl @@ -1,4 +1,5 @@ -function constructCoreOpticalProperties(RS_type::AbstractRamanType{FT}, iBand, m, model) where FT +function constructCoreOpticalProperties(RS_type::AbstractRamanType{FT}, + iBand, m, model) where FT @unpack τ_rayl, τ_aer, τ_abs, aerosol_optics, greek_rayleigh = model #@show typeof(τ_rayl[1]), typeof(τ_aer[1]), typeof(τ_abs[1]) @assert all(iBand .≤ length(τ_rayl)) "iBand exceeded number of bands" @@ -38,6 +39,16 @@ function constructCoreOpticalProperties(RS_type::AbstractRamanType{FT}, iBand, m pol_type, μ, aerosol_optics[iB][i].greek_coefs, m, arr_type=arr_type) + #dAerZ⁺⁺ = zeros(4, size(AerZ⁺⁺,1), size(AerZ⁺⁺,2)) + #dAerZ⁻⁺ = zeros(4, size(AerZ⁺⁺,1), size(AerZ⁺⁺,2)) + + #for ctr=1:4 + # dAerZ⁺⁺[ctr,:,:], dAerZ⁻⁺[ctr,:,:] = + # Scattering.compute_Z_moments( + # pol_type, μ, + # aerosol_optics[iB][i].greek_coefs, + # m, arr_type=arr_type) + #end #@show typeof(AerZ⁺⁺), typeof(aerosol_optics[iB][i]), typeof(FT.(τ_aer[iB][i,:])) # Generate Core optical properties for Aerosols i aer = createAero.(τ_aer[iB][i,:], @@ -121,7 +132,7 @@ function getG_atSun(lod::CoreDirectionalScatteringOpticalProperties,quad_points: end -function expandOpticalProperties(in::CoreScatteringOpticalProperties, arr_type) +function (in::CoreScatteringOpticalProperties, arr_type) @unpack τ, ϖ, Z⁺⁺, Z⁻⁺ = in @assert length(τ) == length(ϖ) "τ and ϖ sizes need to match" if size(Z⁺⁺,3) == 1 diff --git a/src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl b/src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl new file mode 100644 index 00000000..1796c430 --- /dev/null +++ b/src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl @@ -0,0 +1,213 @@ +function constructCoreOpticalProperties( + RS_type::AbstractRamanType{FT}, + iBand, m, model, lin) where FT + @unpack τ_rayl, τ_aer, τ_abs, aerosol_optics, greek_rayleigh = model + @unpack lin_τ_aer, lin_τ_abs, lin_aerosol_optics = lin + + #@show typeof(τ_rayl[1]), typeof(τ_aer[1]), typeof(τ_abs[1]) + @assert all(iBand .≤ length(τ_rayl)) "iBand exceeded number of bands" + + arr_type = array_type(model.params.architecture) + + pol_type = model.params.polarization_type + + # Quadrature points: + μ = Array(model.quad_points.qp_μ ) + N = length(model.quad_points.qp_μN) + # Number of Aerosols: + nAero = size(τ_aer[iBand[1]],1) + nZ = size(τ_rayl[1],2) + # Rayleigh Z matrix: + Rayl𝐙⁺⁺, Rayl𝐙⁻⁺ = Scattering.compute_Z_moments(pol_type, μ, + greek_rayleigh, m, + arr_type = arr_type); + + band_layer_props = []; + band_fScattRayleigh = []; + # @show arr_type + for iB in iBand + # Define G here as ones: + #G = arr_type(ones(FT,N) + rayl = [CoreScatteringOpticalProperties( + arr_type(τ_rayl[iB][:,i]),RS_type.ϖ_Cabannes[iB], + Rayl𝐙⁺⁺, Rayl𝐙⁻⁺) for i=1:nZ] + # Initiate combined properties with rayleigh + combo = rayl + lin_combo = [] + + # Loop over all aerosol types: + for i=1:nAero + # Precomute Z matrices per type (constant per layer) + #@show typeof(aerosol_optics[iB][i].greek_coefs), typeof(pol_type), typeof(μ) + AerZ⁺⁺, AerZ⁻⁺ = Scattering.compute_Z_moments( + pol_type, μ, + aerosol_optics[iB][i].greek_coefs, + m, arr_type=arr_type) + dAerZ⁺⁺ = zeros(4, size(AerZ⁺⁺,1), size(AerZ⁺⁺,2)) + dAerZ⁻⁺ = zeros(4, size(AerZ⁺⁺,1), size(AerZ⁺⁺,2)) + + for ctr=1:4 + dAerZ⁺⁺[ctr,:,:], dAerZ⁻⁺[ctr,:,:] = + Scattering.compute_Z_moments( + pol_type, μ, + lin_aerosol_optics[iB][i].d_greek_coefs[ctr], + m, arr_type=arr_type) + end + #@show typeof(AerZ⁺⁺), typeof(aerosol_optics[iB][i]), typeof(FT.(τ_aer[iB][i,:])) + # Generate Core optical properties for Aerosols i + aer, lin_aer = createAero.(τ_aer[iB][i,:], + [aerosol_optics[iB][i]], + [AerZ⁺⁺], [AerZ⁻⁺], lin_aerosol_optics[iB][i], dAerZ⁺⁺, dAerZ⁻⁺) + #@show typeof(aer), typeof(combo) + # Mix with previous Core Optical Properties + (combo, lin_combo) = (combo, lin_combo) .+ (aer, lin_aer) + end + #@show typeof(combo) + # TODO Type check τ_abs, τ_aer, rayl[i].τ ./ combo[i].τ + # Somewhere here we can add canopy later as well! + ### + + # fScattRayleigh: + # Assume ϖ of 1 for Rayleight here: + #@show size(combo) + #fScattRayleigh = [FT.(Array(rayl[i].τ ./ combo[i].τ)) for i=1:nZ] + fScattRayleigh = [Array(rayl[i].τ ./ combo[i].τ) for i=1:nZ] + lin_fScattRayleigh = [Array(-rayl[i].τ .* lin_combo[i].lin_τ ./ combo[i].τ^2) for i=1:nZ] + + #@show fScattRayleigh, rayl[1].τ, combo[1].τ + # Create Core Optical Properties merged with trace gas absorptions: + #@show typeof(combo.+ + #[CoreAbsorptionOpticalProperties(arr_type((τ_abs[iB][:,i]))) for i=1:nZ]) + + # Gaseous Absorption + gabs = [CoreAbsorptionOpticalProperties( + arr_type(τ_abs[iB][:,i])) for i=1:nZ] + lin_gabs = [linCoreAbsorptionOpticalProperties( + arr_type(lin,τ_abs[iB][:,:,i])) for i=1:nZ] + (combo, lin_combo) = (combo, lin_combo) .+ (gabs, lin_gabs) + push!(band_layer_props, combo) + push!(lin_band_layer_props, lin_combo) + push!(band_fScattRayleigh,fScattRayleigh) + push!(lin_band_fScattRayleigh,lin_fScattRayleigh) + end + + layer_opt = [] + fscat_opt = [] + lin_layer_opt = [] + lin_fscat_opt = [] + for iz = 1:nZ + push!(layer_opt, prod([band_layer_props[i][iz] for i=1:length(iBand)])); + push!(fscat_opt, [band_fScattRayleigh[i][iz] for i=1:length(iBand)]); + push!(lin_layer_opt, prod([lin_band_layer_props[i][:,iz] for i=1:length(iBand)])); + push!(lin_fscat_opt, [lin_band_fScattRayleigh[i][:,iz] for i=1:length(iBand)]); + end + # For now just one band_fScattRayleigh + return layer_opt, fscat_opt, lin_layer_opt, lin_fscat_opt +end + +function createAero(τAer, aerosol_optics, AerZ⁺⁺, AerZ⁻⁺, lin_aerosol_optics, dAerZ⁺⁺, dAerZ⁻⁺) + @unpack k_ref, k, fᵗ, ω̃ = aerosol_optics + @unpack dk_ref, dk, dfᵗ, dω̃ = lin_aerosol_optics + τ_mod = (1-fᵗ * ω̃ ) * τAer; + ϖ_mod = (1-fᵗ) * ω̃/(1-fᵗ * ω̃) + + lin_τ_mod[1] = (τ_mod/τAer)*(k/k_ref) + lin_ϖ_mod[1] = 0 + for ctr=2:5 #ctr=1 corresponds to the derivative with respect to τ_ref, the rest for the microphysical parameters nᵣ, nᵢ, r₀. and σ₀ + mctr = ctr-1 + lin_τ_mod[ctr] = (τ_mod/k)*(dk[mctr] - dk_ref[mctr]*(k/k_ref)) + - τAer*(fᵗ*dω̃[mctr] + dfᵗ[mctr]*ω̃) + lin_ϖ_mod[ctr] = (dω̃[mctr]*(1-fᵗ) - dfᵗ[mctr]*ω̃*(1-ω̃))/(1-fᵗ * ω̃)^2 + end + CoreScatteringOpticalProperties(τ_mod, ϖ_mod,AerZ⁺⁺, AerZ⁻⁺), + linCoreScatteringOpticalProperties(lin_τ_mod, lin_ϖ_mod, dAerZ⁺⁺, dAerZ⁻⁺) +end + +# Extract scattering definitions and integrated absorptions for the source function! +function extractEffectiveProps( + lods::Array,#{CoreScatteringOpticalProperties{FT},1} + lin_lods::Array, + quad_points::QuadPoints{FT} + ) where FT + + #FT = eltype(lods[1].τ) + nSpec = length(lods[1].τ) + nZ = length(lods) + # First the Scattering Interfaces: + scattering_interface = ScatteringInterface_00() + scattering_interfaces_all = [] + τ_sum_all = similar(lods[1].τ,(nSpec,nZ+1)) #?? + #lin_τ_sum_all = similar(lin_lods[1].τ,(nSpec,nZ+1)) #?? + τ_sum_all[:,1] .= 0 + lin_τ_sum_all[:,1] .= 0 + #@show FT + for iz =1:nZ + # Need to check max entries in Z matrices here as well later! + scatter = maximum(lods[iz].τ .* lods[iz].ϖ) > 2eps(FT) + scattering_interface = get_scattering_interface(scattering_interface, scatter, iz) + push!(scattering_interfaces_all, scattering_interface) + #@show typeof(τ_sum_all[:,iz]), typeof(lods[iz].τ) + @views τ_sum_all[:,iz+1] = τ_sum_all[:,iz] + getG_atSun(lods[iz], quad_points) * lods[iz].τ + for ctr = 1:size(lin_τ_sum_all,1) + @views lin_τ_sum_all[ctr,:,iz+1] = lin_τ_sum_all[ctr,:,iz] + getG_atSun(lods[iz], quad_points) * lods[iz].lin_τ[ctr,:] + end + end + return scattering_interfaces_all, τ_sum_all, lin_τ_sum_all +end +#= +function getG_atSun(lod::CoreScatteringOpticalProperties,quad_points::QuadPoints{FT}) where FT + FT(1) +end + +function getG_atSun(lod::CoreDirectionalScatteringOpticalProperties,quad_points::QuadPoints{FT}) where FT + @unpack iμ₀ = quad_points + gfct = Array(lod.G)[iμ₀] + return gfct +end +=# + +function (in::linCoreScatteringOpticalProperties, arr_type) + @unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺ = in + @assert size(lin_τ) == size(lin_ϖ) "τ and ϖ sizes need to match" + if size(lin_Z⁺⁺,4) == 1 + lin_Z⁺⁺ = _repeat(lin_Z⁺⁺,1,1,1,length(lin_τ[1,:])) + lin_Z⁻⁺ = _repeat(lin_Z⁻⁺,1,1,1,length(lin_τ[1,:])) + return linCoreScatteringOpticalProperties(arr_type(τ), arr_type(ϖ), arr_type(Z⁺⁺), arr_type(Z⁻⁺)) + else + @assert size(lin_Z⁺⁺,4) == length(lin_τ[1,:]) "Z and τ dimensions need to match " + linCoreScatteringOpticalProperties(arr_type(lin_τ), arr_type(lin_ϖ), arr_type(lin_Z⁺⁺), arr_type(lin_Z⁻⁺)) + end +end + +function expandOpticalProperties(in::linCoreDirectionalScatteringOpticalProperties, arr_type) + @unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺, lin_G = in + @assert size(lin_τ) == size(lin_ϖ) "τ and ϖ sizes need to match" + if size(lin_Z⁺⁺,4) == 1 + lin_Z⁺⁺ = _repeat(lin_Z⁺⁺,1,1,1,length(lin_τ[1,:])) + lin_Z⁻⁺ = _repeat(lin_Z⁻⁺,1,1,1,length(lin_τ[1,:])) + return linCoreDirectionalScatteringOpticalProperties( + arr_type(lin_τ), arr_type(lin_ϖ), + arr_type(lin_Z⁺⁺), arr_type(lin_Z⁻⁺), + arr_type(lin_G)) + else + @assert size(lin_Z⁺⁺,4) == length(lin_τ[1,:]) "Z and τ dimensions need to match " + linCoreDirectionalScatteringOpticalProperties( + arr_type(lin_τ), arr_type(lin_ϖ), + arr_type(lin_Z⁺⁺), arr_type(lin_Z⁻⁺), + arr_type(lin_G)) + end +end + +#= +function expandBandScalars(RS_type, x) +#test = [length(RS_type.bandSpecLim[iB]) for iB in RS_type.iBand] +#@show test, sum(test), size(x[1]) +#@show eltype(x[1]),sum([length(RS_type.bandSpecLim[iB]) for iB in RS_type.iBand]) +out = zeros(eltype(x[1]),sum([length(RS_type.bandSpecLim[iB]) for iB in RS_type.iBand])) +for iB in RS_type.iBand +out[RS_type.bandSpecLim[iB]] .= expandScalar(x[iB],length(RS_type.bandSpecLim[iB])) +end +return out +end +=# +#expandScalar(x,n) = x.*ones(n); \ No newline at end of file diff --git a/src/CoreRT/Surfaces/lin_lambertian_surface.jl b/src/CoreRT/Surfaces/lin_lambertian_surface.jl new file mode 100644 index 00000000..c1ff2ae8 --- /dev/null +++ b/src/CoreRT/Surfaces/lin_lambertian_surface.jl @@ -0,0 +1,198 @@ +#= + +This file specifies how to create surface layers, given the surface type, and related info + +=# + +""" + $(FUNCTIONNAME)(lambertian::LambertianSurfaceScalar{FT}) + +Computes (in place) surface optical properties for a (scalar) lambertian albedo as [`AddedLayer`](@ref) + + - `lambertian` a [`LambertianSurfaceScalar`](@ref) struct that defines albedo as scalar + - `SFI` bool if SFI is used + - `m` Fourier moment (starting at 0) + - `pol_type` Polarization type struct + - `quad_points` Quadrature points struct + - `τ_sum` total optical thickness from TOA to the surface + - `architecture` Compute architecture (GPU,CPU) +""" +function create_surface_layer!(lambertian::LambertianSurfaceScalar{FT}, + added_layer::Union{AddedLayer,AddedLayerRS}, + lin_added_layer::Union{linAddedLayer,linAddedLayerRS}, + SFI, + m::Int, + pol_type, + quad_points, + τ_sum, + lin_τ_sum, + Nparams, i_surf, + architecture) where {FT} + + @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart, iμ₀, μ₀ = quad_points + # Get size of added layer + Nquad = size(added_layer.r⁻⁺,1) ÷ pol_type.n + tmp = ones(pol_type.n*Nquad) + arr_type = array_type(architecture) + T_surf = arr_type(Diagonal(tmp)) + if m == 0 + # Albedo normalized by π (and factor 2 for 0th Fourier Moment) + ρ = 2lambertian.albedo#/FT(π) + + + + R_surf = Array(Diagonal(vcat(ρ, zeros(FT,pol_type.n-1)))) + R_surf = repeat(R_surf',Nquad) + R_surf = repeat(R_surf',Nquad) + + # Move to architecture: + R_surf = arr_type(R_surf) + + + # Source function of surface: + if SFI + I₀_NquadN = similar(qp_μN); + I₀_NquadN[:] .=0; + I₀_NquadN[iμ₀Nstart:pol_type.n*iμ₀] = pol_type.I₀; + + added_layer.j₀⁺[:,1,:] .= 0; #I₀_NquadN .* exp.(-τ_sum/μ₀)'; + added_layer.j₀⁻[:,1,:] .= μ₀*(R_surf*I₀_NquadN) .* exp.(-τ_sum/μ₀)'; + for ctr=1:Nparams + if ctr==i_surf + lin_added_layer.dxj₀⁺[ctr,:,1,:] .= 0;#I₀_NquadN .* exp.(-τ_sum/μ₀)'/lambertian.ρ; + lin_added_layer.dxj₀⁻[ctr,:,1,:] .= μ₀*(R_surf*I₀_NquadN) .* exp.(-τ_sum/μ₀)'/lambertian.ρ; + else + lin_added_layer.dxj₀⁺[ctr,:,1,:] .= 0; #I₀_NquadN .* exp.(-τ_sum/μ₀)'*(-1/μ₀).*lin_τ_sum[ctr,:]; + lin_added_layer.dxj₀⁻[ctr,:,1,:] .= μ₀*(R_surf*I₀_NquadN) .* exp.(-τ_sum/μ₀)'*(-1/μ₀).*lin_τ_sum[ctr,:]; + end + end + end + R_surf = R_surf * Diagonal(qp_μN.*wt_μN) + + + #@show size(added_layer.r⁻⁺), size(R_surf), size(added_layer.j₀⁻) + added_layer.r⁻⁺ .= R_surf; + added_layer.r⁺⁻ .= 0; + added_layer.t⁺⁺ .= T_surf; + added_layer.t⁻⁻ .= T_surf; + + lin_added_layer.dxr⁻⁺ .= 0; + lin_added_layer.dxr⁺⁻ .= 0; + lin_added_layer.dxt⁺⁺ .= 0; + lin_added_layer.dxt⁻⁻ .= 0; + lin_added_layer.dxr⁻⁺[i_surf,:,:,:] .= R_surf./lambertian.ρ; + + + else + added_layer.r⁻⁺ .= 0; + added_layer.r⁻⁺ .= 0; + added_layer.t⁺⁺ .= T_surf; + added_layer.t⁻⁻ .= T_surf; + added_layer.j₀⁺ .= 0; + added_layer.j₀⁻ .= 0; + + lin_added_layer.dxr⁻⁺ .= 0; + lin_added_layer.dxr⁺⁻ .= 0; + lin_added_layer.dxt⁺⁺ .= 0; + lin_added_layer.dxt⁻⁻ .= 0; + lin_added_layer.dxj₀⁻ .= 0; + lin_added_layer.dxj₀⁺ .= 0; + end +end + +function create_surface_layer!(lambertian::LambertianSurfaceLegendre{FT}, + added_layer::Union{AddedLayer,AddedLayerRS}, + lin_added_layer::Union{linAddedLayer,linAddedLayerRS}, + SFI, + m::Int, + pol_type, + quad_points, + τ_sum, + lin_τ_sum, + Nparams, i_surf, #here i_surf is an array of indices + architecture) where {FT} + FT2 = Float64 + if m == 0 + @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart, iμ₀, μ₀ = quad_points + legendre_coeff = lambertian.legendre_coeff + arr_type = array_type(architecture) + # Albedo normalized by π (and factor 2 for 0th Fourier Moment) + # a) Define range for legendre polymonial: + x = collect(range(FT2(-1), FT2(1), length=length(τ_sum))); + # Legendre Polynomial basis functions: + P = Scattering.compute_legendre_poly(x,length(legendre_coeff))[1] + # Evaluate Polynomial (as matrix multiplication) + albedo = P * legendre_coeff + ρ = arr_type(2albedo) + # Get size of added layer + dim = size(added_layer.r⁻⁺) + Nquad = dim[1] ÷ pol_type.n + + R_surf = Array(Diagonal(vcat(FT(1), zeros(FT,pol_type.n-1)))) + R_surf = repeat(R_surf',Nquad) + R_surf = repeat(R_surf',Nquad) + + # Move to architecture: + R_surf = arr_type(R_surf) + + # Source function of surface: + if SFI + I₀_NquadN = similar(qp_μN); + I₀_NquadN[:] .=0; + I₀_NquadN[iμ₀Nstart:pol_type.n*iμ₀] = pol_type.I₀; + added_layer.j₀⁺[:] .= 0 + # Suniti double-check + added_layer.j₀⁻[:,1,:] = μ₀*(R_surf*I₀_NquadN) .* (ρ .* exp.(-τ_sum/μ₀))'; + for ctr=1:Nparams + if ctr∈i_surf + lin_added_layer.dxj₀⁺[:,1,:] .= 0;#I₀_NquadN .* exp.(-τ_sum/μ₀)'/lambertian.ρ; + lin_added_layer.dxj₀⁻[:,1,:] .= 2μ₀*(R_surf*I₀_NquadN) .* P[:,ctr-i_surf[1]+1] .* exp.(-τ_sum/μ₀)'; + else + lin_added_layer.dxj₀⁺[ctr,:,1,:] .= 0; #I₀_NquadN .* exp.(-τ_sum/μ₀)'*(-1/μ₀).*lin_τ_sum[ctr,:]; + lin_added_layer.dxj₀⁻[ctr,:,1,:] .= μ₀*(R_surf*I₀_NquadN) .* ρ .* exp.(-τ_sum/μ₀)'*(-1/μ₀).*lin_τ_sum[ctr,:]; + end + end + + end + R_surf = R_surf * Diagonal(qp_μN.*wt_μN) + siz = size(added_layer.r⁻⁺) + R_surf3D = reshape(reduce(hcat,[i*R_surf for i in Array(ρ)]), siz...); + tmp = ones(pol_type.n*Nquad) + T_surf = arr_type(Diagonal(tmp)) + + #@show size(added_layer.r⁻⁺), size(R_surf), size(added_layer.j₀⁻) + added_layer.r⁻⁺ .= R_surf3D; + added_layer.r⁺⁻ .= 0; + added_layer.t⁺⁺ .= T_surf; + added_layer.t⁻⁻ .= T_surf; + + lin_added_layer.dxr⁻⁺ .= 0; + lin_added_layer.dxr⁺⁻ .= 0; + lin_added_layer.dxt⁺⁺ .= 0; + lin_added_layer.dxt⁻⁻ .= 0; + for ctr in i_surf #all elements in i_surf are assumed to be contiguous, e.g. 10, 11, 12 + lin_added_layer.dxr⁻⁺[ctr,:,:,:] .= reshape(reduce(hcat,[i*R_surf for i in Array(P[:,ctr-i_surf[1]+1])]), siz...); + end + else + added_layer.r⁻⁺[:] .= 0; + added_layer.r⁻⁺[:] .= 0; + added_layer.t⁺⁺[:] .= 0; + added_layer.t⁻⁻[:] .= 0; + added_layer.j₀⁺[:] .= 0; + added_layer.j₀⁻[:] .= 0; + + lin_added_layer.dxr⁻⁺ .= 0; + lin_added_layer.dxr⁺⁻ .= 0; + lin_added_layer.dxt⁺⁺ .= 0; + lin_added_layer.dxt⁻⁻ .= 0; + lin_added_layer.dxj₀⁻ .= 0; + lin_added_layer.dxj₀⁺ .= 0; + end +end + +#function reflectance(sur::LambertianSurfaceScalar{FT}, μᵢ::FT, μᵣ::FT, dϕ::FT) where FT +# return sur.albedo +#end + + + diff --git a/src/CoreRT/lin_model_from_parameters.jl b/src/CoreRT/lin_model_from_parameters.jl index 0023966f..716526d6 100644 --- a/src/CoreRT/lin_model_from_parameters.jl +++ b/src/CoreRT/lin_model_from_parameters.jl @@ -149,6 +149,8 @@ function model_from_parameters(params::vSmartMOM_Parameters) @timeit "Mie calc" aerosol_optics_raw, lin_aerosol_optics_raw = compute_aerosol_optical_properties(mie_model, FT); @show aerosol_optics_raw.k + aerosol_optics_raw.k_ref = k_ref + lin_aerosol_optics_raw.dk_ref = dk_ref # Compute truncated aerosol optical properties (phase function and fᵗ), consistent with Ltrunc: #@show i_aer, i_band aerosol_optics[i_band][i_aer] = Scattering.truncate_phase(truncation_type, diff --git a/src/CoreRT/lin_types.jl b/src/CoreRT/lin_types.jl new file mode 100644 index 00000000..9eab8a68 --- /dev/null +++ b/src/CoreRT/lin_types.jl @@ -0,0 +1,786 @@ +#= + +This file contains the linearization of all relevant types that are used in the vSmartMOM module: + +- `AtmosphericProfile` stores all relevant atmospheric profile information +- `AbstractObsGeometry` specifies the RT geometry +- `RT_Aerosol` holds an Aerosol with additional RT parameters +- `AbstractQuadratureType` specifies the quadrature type to use +- `AbstractSourceType` specifies the source type +- `CompositeLayer` and `AddedLayer` specify the layer properties +- `AbstractScatteringInterface` specifies the scattering interface type +- `AbstractSurfaceType` specify the type of surface in the RT simulation +- `AbsorptionParameters`, `ScatteringParameters`, and `vSmartMOM_Model` hold model parameters +- `QuadPoints` holds quadrature points, weights, etc. +- `ComputedAtmosphereProperties` and `ComputedLayerProperties` hold intermediate computed properties + +=# +#= +"Struct for an atmospheric profile" +struct AtmosphericProfile{FT, VMR <: Union{Real, Vector}} + "Temperature Profile" + T::Array{FT,1} + "Pressure Profile (Full)" + p_full::Array{FT,1} + "Specific humidity profile" + q::Array{FT,1} + "Pressure Levels" + p_half::Array{FT,1} + "H2O Volume Mixing Ratio Profile" + vmr_h2o::Array{FT,1} + "Vertical Column Density (Dry)" + vcd_dry::Array{FT,1} + "Vertical Column Density (H2O)" + vcd_h2o::Array{FT,1} + "Volume Mixing Ratio of Constituent Gases" + vmr::Dict{String, VMR} + "Layer height (meters)" + Δz::Array{FT,1} +end + +"Types for describing atmospheric parameters" +abstract type AbstractObsGeometry end + +"Observation Geometry (basics)" +Base.@kwdef struct ObsGeometry{FT} <: AbstractObsGeometry + "Solar Zenith Angle `[Degree]`" + sza::FT + "Viewing Zenith Angle(s) `[Degree]`" + vza::Array{FT,1} + "Viewing Azimuth Angle(s) `[Degree]`" + vaz::Array{FT,1} + "Altitude of Observer `[Pa]`" + obs_alt::FT +end +=# + +mutable struct RT_Aerosol{}#FT<:Union{AbstractFloat, ForwardDiff.Dual}} + "Aerosol" + aerosol::Aerosol#{FT} + "Reference τ" + τ_ref#::FT + "Vertical distribution as function of p (using Distributions.jl)" + profile::Distribution#::FT +end +#= +"Quadrature Types for RT streams" +abstract type AbstractQuadratureType end + +""" + struct RadauQuad + +Use Gauss Radau quadrature scheme, which includes the SZA as point (see Sanghavi vSmartMOM) + +""" +struct RadauQuad <:AbstractQuadratureType end + +""" + struct GaussQuadHemisphere + +Use Gauss quadrature scheme, define interval [-1,1] within an hemisphere (90⁰), repeat for both + +""" +struct GaussQuadHemisphere <:AbstractQuadratureType end + +""" + struct GaussQuadFullSphere + +Use Gauss quadrature scheme, define interval [-1,1] for full sphere (180⁰), take half of it (less points near horizon compared to GaussQuadHemisphere) + +""" +struct GaussQuadFullSphere <:AbstractQuadratureType end + +"Abstract Type for Source Function Integration" +abstract type AbstractSourceType end + +"Use Dummy Node Integration (DNI), SZA has to be a full node, Radau required" +struct DNI <:AbstractSourceType end + +"Use Source Function Integration (SFI), Solar beam embedded in source term, can work with all quadrature schemes (faster)" +struct SFI <:AbstractSourceType end + +"Abstract Type for Layer R,T and J matrices" +abstract type AbstractLayer end +=# +"Composite Layer Matrices (`-/+` defined in τ coordinates, i.e. `-`=outgoing, `+`=incoming" +Base.@kwdef struct linCompositeLayer{FT} <: AbstractLayer + "Composite layer Reflectance matrix R (from + -> -)" + dR⁻⁺::AbstractArray{FT,3} + "Composite layer Reflectance matrix R (from - -> +)" + dR⁺⁻::AbstractArray{FT,3} + "Composite layer transmission matrix T (from + -> +)" + dT⁺⁺::AbstractArray{FT,3} + "Composite layer transmission matrix T (from - -> -)" + dT⁻⁻::AbstractArray{FT,3} + "Composite layer source matrix J (in + direction)" + dJ₀⁺::AbstractArray{FT,3} + "Composite layer source matrix J (in - direction)" + dJ₀⁻::AbstractArray{FT,3} +end + +"Added (Single) Layer Matrices (`-/+` defined in τ coordinates, i.e. `-`=outgoing, `+`=incoming" +Base.@kwdef struct linAddedLayer{FT} <: AbstractLayer + "Added layer Reflectance matrix R (from + -> -)" + dr⁻⁺::AbstractArray{FT,3} + "Added layer transmission matrix T (from + -> +)" + dt⁺⁺::AbstractArray{FT,3} + "Added layer Reflectance matrix R (from - -> +)" + dr⁺⁻::AbstractArray{FT,3} + "Added layer transmission matrix T (from - -> -)" + dt⁻⁻::AbstractArray{FT,3} + "Added layer source matrix J (in + direction)" + dj₀⁺::AbstractArray{FT,3} + "Added layer source matrix J (in - direction)" + dj₀⁻::AbstractArray{FT,3} + "Added layer Reflectance matrix R (from + -> -)" + dxr⁻⁺::AbstractArray{FT,3} + "Added layer transmission matrix T (from + -> +)" + dxt⁺⁺::AbstractArray{FT,3} + "Added layer Reflectance matrix R (from - -> +)" + dxr⁺⁻::AbstractArray{FT,3} + "Added layer transmission matrix T (from - -> -)" + dxt⁻⁻::AbstractArray{FT,3} + "Added layer source matrix J (in + direction)" + dxj₀⁺::AbstractArray{FT,3} + "Added layer source matrix J (in - direction)" + dxj₀⁻::AbstractArray{FT,3} +end + +"Composite Layer Matrices (`-/+` defined in τ coordinates, i.e. `-`=outgoing, `+`=incoming" +struct linCompositeLayerRS{FT} <: AbstractLayer + "Composite layer Reflectance matrix R (from + -> -)" + dR⁻⁺::AbstractArray{FT,3} + "Composite layer Reflectance matrix R (from - -> +)" + dR⁺⁻::AbstractArray{FT,3} + "Composite layer transmission matrix T (from + -> +)" + dT⁺⁺::AbstractArray{FT,3} + "Composite layer transmission matrix T (from - -> -)" + dT⁻⁻::AbstractArray{FT,3} + "Composite layer source matrix J (in + direction)" + dJ₀⁺::AbstractArray{FT,3} + "Composite layer source matrix J (in - direction)" + dJ₀⁻::AbstractArray{FT,3} + + # Additional Arrays for Raman scattering + "Composite layer Reflectance matrix ieR (from + -> -)" + dieR⁻⁺::AbstractArray{FT,4} + "Composite layer Reflectance matrix ieR (from - -> +)" + dieR⁺⁻::AbstractArray{FT,4} + "Composite layer transmission matrix ieT (from + -> +)" + dieT⁺⁺::AbstractArray{FT,4} + "Composite layer transmission matrix ieT (from - -> -)" + dieT⁻⁻::AbstractArray{FT,4} + "Composite layer source matrix ieJ (in + direction)" + dieJ₀⁺::AbstractArray{FT,4} + "Composite layer source matrix ieJ (in - direction)" + dieJ₀⁻::AbstractArray{FT,4} +end + +"Added (Single) Layer Matrices (`-/+` defined in τ coordinates, i.e. `-`=outgoing, `+`=incoming" +struct linAddedLayerRS{FT} <: AbstractLayer + "Added layer Reflectance matrix R (from + -> -)" + dr⁻⁺::AbstractArray{FT,3} + "Added layer transmission matrix T (from + -> +)" + dt⁺⁺::AbstractArray{FT,3} + "Added layer Reflectance matrix R (from - -> +)" + dr⁺⁻::AbstractArray{FT,3} + "Added layer transmission matrix T (from - -> -)" + dt⁻⁻::AbstractArray{FT,3} + "Added layer source matrix J (in + direction)" + dj₀⁺::AbstractArray{FT,3} + "Added layer source matrix J (in - direction)" + dj₀⁻::AbstractArray{FT,3} + + # Additional Arrays for Raman scattering + "Added layer Reflectance matrix ieR (from + -> -)" + dier⁻⁺::AbstractArray{FT,4} + "Added layer transmission matrix ieT (from + -> +)" + diet⁺⁺::AbstractArray{FT,4} + "Added layer Reflectance matrix ieR (from - -> +)" + dier⁺⁻::AbstractArray{FT,4} + "Added layer transmission matrix ieT (from - -> -)" + diet⁻⁻::AbstractArray{FT,4} + "Added layer source matrix ieJ (in + direction)" + diej₀⁺::AbstractArray{FT,4} + "Added layer source matrix ieJ (in - direction)" + diej₀⁻::AbstractArray{FT,4} + + "Added layer Reflectance matrix R (from + -> -)" + dxr⁻⁺::AbstractArray{FT,3} + "Added layer transmission matrix T (from + -> +)" + dxt⁺⁺::AbstractArray{FT,3} + "Added layer Reflectance matrix R (from - -> +)" + dxr⁺⁻::AbstractArray{FT,3} + "Added layer transmission matrix T (from - -> -)" + dxt⁻⁻::AbstractArray{FT,3} + "Added layer source matrix J (in + direction)" + dxj₀⁺::AbstractArray{FT,3} + "Added layer source matrix J (in - direction)" + dxj₀⁻::AbstractArray{FT,3} + + # Additional Arrays for Raman scattering + "Added layer Reflectance matrix ieR (from + -> -)" + diexr⁻⁺::AbstractArray{FT,4} + "Added layer transmission matrix ieT (from + -> +)" + diext⁺⁺::AbstractArray{FT,4} + "Added layer Reflectance matrix ieR (from - -> +)" + diexr⁺⁻::AbstractArray{FT,4} + "Added layer transmission matrix ieT (from - -> -)" + diext⁻⁻::AbstractArray{FT,4} + "Added layer source matrix ieJ (in + direction)" + diexj₀⁺::AbstractArray{FT,4} + "Added layer source matrix ieJ (in - direction)" + diexj₀⁻::AbstractArray{FT,4} +end +# Multisensor Composite layers +# Elastic +"Composite Layer Matrices (`-/+` defined in τ coordinates, i.e. `-`=outgoing, `+`=incoming" +Base.@kwdef struct linCompositeLayerMS{M} <: AbstractLayer + "Composite layer Reflectance matrix R (from + -> -)" + dtopR⁻⁺::M#AbstractArray{FT,4} + "Composite layer Reflectance matrix R (from - -> +)" + dtopR⁺⁻::M + "Composite layer transmission matrix T (from + -> +)" + dtopT⁺⁺::M + "Composite layer transmission matrix T (from - -> -)" + dtopT⁻⁻::M + "Composite layer source matrix J (in + direction)" + dtopJ₀⁺::M + "Composite layer source matrix J (in - direction)" + dtopJ₀⁻::M + "Composite layer Reflectance matrix R (from + -> -)" + dbotR⁻⁺::M + "Composite layer Reflectance matrix R (from - -> +)" + dbotR⁺⁻::M + "Composite layer transmission matrix T (from + -> +)" + dbotT⁺⁺::M + "Composite layer transmission matrix T (from - -> -)" + dbotT⁻⁻::M + "Composite layer source matrix J (in + direction)" + dbotJ₀⁺::M + "Composite layer source matrix J (in - direction)" + dbotJ₀⁻::M +end + +# Inelastic +"Composite Layer Matrices (`-/+` defined in τ coordinates, i.e. `-`=outgoing, `+`=incoming" +Base.@kwdef struct linCompositeLayerMSRS{M1, M2} <: AbstractLayer + "Composite layer Reflectance matrix R (from + -> -)" + dtopR⁻⁺::M1 + "Composite layer Reflectance matrix R (from - -> +)" + dtopR⁺⁻::M1 + "Composite layer transmission matrix T (from + -> +)" + dtopT⁺⁺::M1 + "Composite layer transmission matrix T (from - -> -)" + dtopT⁻⁻::M1 + "Composite layer source matrix J (in + direction)" + dtopJ₀⁺::M1 + "Composite layer source matrix J (in - direction)" + dtopJ₀⁻::M1 + + # Additional Arrays for Raman scattering + "Composite layer Reflectance matrix ieR (from + -> -)" + dtopieR⁻⁺::M2 + "Composite layer Reflectance matrix ieR (from - -> +)" + dtopieR⁺⁻::M2 + "Composite layer transmission matrix ieT (from + -> +)" + dtopieT⁺⁺::M2 + "Composite layer transmission matrix ieT (from - -> -)" + dtopieT⁻⁻::M2 + "Composite layer source matrix ieJ (in + direction)" + dtopieJ₀⁺::M2 + "Composite layer source matrix ieJ (in - direction)" + dtopieJ₀⁻::M2 + + "Composite layer Reflectance matrix R (from + -> -)" + dbotR⁻⁺::M1 + "Composite layer Reflectance matrix R (from - -> +)" + dbotR⁺⁻::M1 + "Composite layer transmission matrix T (from + -> +)" + dbotT⁺⁺::M1 + "Composite layer transmission matrix T (from - -> -)" + dbotT⁻⁻::M1 + "Composite layer source matrix J (in + direction)" + dbotJ₀⁺::M1 + "Composite layer source matrix J (in - direction)" + dbotJ₀⁻::M1 + + # Additional Arrays for Raman scattering + "Composite layer Reflectance matrix ieR (from + -> -)" + dbotieR⁻⁺::M2 + "Composite layer Reflectance matrix ieR (from - -> +)" + dbotieR⁺⁻::M2 + "Composite layer transmission matrix ieT (from + -> +)" + dbotieT⁺⁺::M2 + "Composite layer transmission matrix ieT (from - -> -)" + dbotieT⁻⁻::M2 + "Composite layer source matrix ieJ (in + direction)" + dbotieJ₀⁺::M2 + "Composite layer source matrix ieJ (in - direction)" + dbotieJ₀⁻::M2 +end +#= +"Abstract Type for Scattering Interfaces" +abstract type AbstractScatteringInterface end +"No scattering in either the added layer or the composite layer" +struct ScatteringInterface_00 <: AbstractScatteringInterface end +"No scattering in inhomogeneous composite layer; Scattering in homogeneous layer, added to bottom of the composite layer." +struct ScatteringInterface_01 <: AbstractScatteringInterface end +"Scattering in inhomogeneous composite layer; No scattering in homogeneous layer which is added to bottom of the composite layer." +struct ScatteringInterface_10 <: AbstractScatteringInterface end +"Scattering in inhomogeneous composite layer; Scattering in homogeneous layer, added to bottom of the composite layer." +struct ScatteringInterface_11 <: AbstractScatteringInterface end +"Scattering Interface between Surface and Composite Layer" +struct ScatteringInterface_AtmoSurf <: AbstractScatteringInterface end + +"Abstract Type for Surface Types" +abstract type AbstractSurfaceType end + +"Lambertian Surface (scalar per band)" +struct LambertianSurfaceScalar{FT} <: AbstractSurfaceType + "Albedo (scalar)" + albedo::FT +end + +"Defined as Array (has to have the same length as the band!)" +struct LambertianSurfaceSpectrum{FT} <: AbstractSurfaceType + "Albedo (vector)" + albedo::AbstractArray{FT,1} +end + +"Defined as Array (has to have the same length as the band!)" +struct rpvSurfaceScalar{FT} <: AbstractSurfaceType + "Overall reflectance level parameter (scalar)" + ρ₀::FT + "Hotspot function parameter (1.0 = no hotspot)" + ρ_c::FT + "Anisotropy shape parameter. k < 1.0 (> 1.0) corresponds to a bowl (bell) shape." + k::FT + "Asymmetry parameter, Θ < 0.0 (> 0.0) corresponds to a predominantly backward (forward) scattering." + Θ::FT +end + +struct RossLiSurfaceScalar{FT} <: AbstractSurfaceType + "Volumetric RossThick fraction" + fvol::FT + "Geometric LiSparse fraction" + fgeo::FT + "Isotropic reflectance fraction" + fiso::FT +end + +"Defined by Legendre polynomial terms as function of spectral grid, which is scaled to [-1,1] (degree derived from length of `a_coeff`)" +struct LambertianSurfaceLegendre{FT} <: AbstractSurfaceType + "albedo = legendre_coeff[1] * P₀ + legendre_coeff[2]*P₁ + legendre_coeff[3]*P₂ + ... " + legendre_coeff::AbstractArray{FT,1} +end +=# +#= +""" + struct AbsorptionParameters + +A struct which holds all absorption-related parameters (before any computations) +""" +mutable struct AbsorptionParameters + "Molecules to use for absorption calculations (`nBand, nMolecules`)" + molecules::AbstractArray + "Volume-Mixing Ratios" + vmr::Dict + "Type of broadening function (Doppler/Lorentz/Voigt)" + broadening_function::AbstractBroadeningFunction + "Complex Error Function to use in Voigt calculations" + CEF::AbstractComplexErrorFunction + "Wing cutoff to use in cross-section calculation (cm⁻¹)" + wing_cutoff::Integer + "Lookup table type" + luts::AbstractArray +end +=# +#= +""" + struct ScatteringParameters + +A struct which holds all scattering-related parameters (before any computations) +""" +mutable struct ScatteringParameters{FT<:Union{AbstractFloat, ForwardDiff.Dual}} + "List of scattering aerosols and their properties" + rt_aerosols::Vector{RT_Aerosol} + "Maximum aerosol particle radius for quadrature points/weights (µm)" + r_max::FT + "Number of quadrature points for integration of size distribution" + nquad_radius::Integer + "Reference wavelength (µm)" + λ_ref::FT + "Reference refractive index" + n_ref::Complex{FT} + "Algorithm to use for fourier decomposition (NAI2/PCW)" + decomp_type::AbstractFourierDecompositionType +end +=# +#= +""" + struct vSmartMOM_Parameters + +A struct which holds all initial model parameters (before any computations) + +# Fields +$(DocStringExtensions.FIELDS) +""" +mutable struct vSmartMOM_Parameters{FT<:Union{AbstractFloat, ForwardDiff.Dual}} + + # radiative_transfer group + "Spectral bands (`nBand`)" + spec_bands::AbstractArray + "Surface (Bidirectional Reflectance Distribution Function)" + brdf::AbstractArray + "Quadrature type for RT streams (RadauQuad/GaussQuadHemisphere/GaussQuadFullSphere)" + quadrature_type::AbstractQuadratureType + "Type of polarization (I/IQ/IQU/IQUV)" + polarization_type::AbstractPolarizationType + "Hard cutoff for maximum number of Fourier moments to loop over" + max_m::Integer + "Exclusion angle for forward peak [deg]" + Δ_angle::FT + "Truncation length for legendre terms (scalar for now, can do `nBand` later)" + l_trunc::Integer + "Depolarization factor" + depol::FT + "Float type to use in the RT (Float64/Float32)" + float_type::DataType + "Architecture to use for calculations (CPU/GPU)" + architecture::AbstractArchitecture + + # geometry group + "Solar zenith angle [deg]" + sza::FT + "Viewing zenith angles [deg]" + vza::AbstractArray{FT} + "Viewing azimuthal angles [deg]" + vaz::AbstractArray{FT} + "Altitude of observer [Pa]" + obs_alt::FT + + # atmospheric_profile group + "Temperature Profile [K]" + T::AbstractArray{FT} + "Pressure Profile [hPa]" + p::AbstractArray{FT} + "Specific humidity profile" + q::AbstractArray{FT} + "Length of profile reduction" + profile_reduction_n::Integer + + # absorption group + "Optional struct that holds all absorption-related parameters" + absorption_params::Union{AbsorptionParameters, Nothing} + + # scattering group + "Optional struct that holds all aerosol scattering-related parameters" + scattering_params::Union{ScatteringParameters, Nothing} + +end +=# +#= +""" + struct QuadPoints{FT} + +A struct which holds Quadrature points, weights, etc + +# Fields +$(DocStringExtensions.FIELDS) +""" +struct QuadPoints{FT} + "μ₀, cos(SZA)" + μ₀::FT + "Index in quadrature points with sun" + iμ₀::Int + "Index in quadrature points with sun (in qp_μN)" + iμ₀Nstart::Int + "Quadrature points" + qp_μ::AbstractArray{FT,1} + "Weights of quadrature points" + wt_μ::AbstractArray{FT,1} + "Quadrature points (repeated for polarizations)" + qp_μN::AbstractArray{FT,1} + "Weights of quadrature points (repeated for polarizations)" + wt_μN::AbstractArray{FT,1} + "Number of quadrature points" + Nquad::Int +end +=# +""" + struct vSmartMOM_Model + +A struct which holds all derived model parameters (including any computations) + +# Fields +$(DocStringExtensions.FIELDS) +""" +mutable struct vSmartMOM_lin{AE, TAB, TR, TAE, Ogeom, PRO} + + #"Struct with all individual parameters" + #params::PA # vSmartMOM_Parameters + + "Truncated aerosol optics" + lin_aerosol_optics::AE # AbstractArray{AbstractArray{AerosolOptics}} + #"Greek coefs in Rayleigh calculations" + #greek_rayleigh::GR # GreekCoefs + #"Quadrature points/weights, etc" + #quad_points::QP # QuadPoints + + "Array to hold cross-sections over entire atmospheric profile" + lin_τ_abs::TAB # AbstractArray{AbstractArray} + #"Rayleigh optical thickness" + lin_τ_rayl::TR # AbstractArray{AbstractArray} + "Aerosol optical thickness" + lin_τ_aer::TAE # AbstractArray{AbstractArray} + + #"Observational Geometry (includes sza, vza, vaz)" + #obs_geom::Ogeom # ObsGeometry + "Atmospheric profile to use" + lin_profile::PRO #AtmosphericProfile +end + +""" + struct ComputedAtmosphereProperties + +A struct which holds (for the entire atmosphere) all key layer optical properties required for the RT core solver + +# Fields +$(DocStringExtensions.FIELDS) +""" +Base.@kwdef struct linComputedAtmosphereProperties + + "Absorption optical depth vectors (wavelength dependent)" + lin_τ_λ_all + "Albedo vectors (wavelength dependent)" + lin_ϖ_λ_all + "Absorption optical depth scalars (not wavelength dependent)" + lin_τ_all + "Albedo scalars (not wavelength dependent)" + lin_ϖ_all + "Combined Z moments (forward)" + lin_Z⁺⁺_all + "Combined Z moments (backward)" + lin_Z⁻⁺_all + "Maximum dτs" + lin_dτ_max_all + "dτs" + lin_dτ_all + #"Number of doublings (for all layers)" + #ndoubl_all + "dτs (wavelength dependent)" + lin_dτ_λ_all + "All expk" + lin_expk_all + #"Scattering flags" + #scatter_all + "Sum of optical thicknesses of all layers above the current layer" + lin_τ_sum_all + #"elastic (Cabannes) scattering fraction of Rayleigh (Cabannes+Raman) scattering per layer" + #ϖ_Cabannes_all + "Rayleigh fraction of scattering cross section per layer" + lin_fscattRayl_all + #"Scattering interface type for each layer" + #scattering_interfaces_all +end + + + +""" + struct ComputedLayerProperties + +A struct which holds all key layer optical properties required for the RT core solver + +# Fields +$(DocStringExtensions.FIELDS) +""" +Base.@kwdef struct linComputedLayerProperties + + "Absorption optical depth vector (wavelength dependent)" + lin_τ_λ + "Albedo vector (wavelength dependent)" + lin_ϖ_λ + "Absorption optical depth scalar (not wavelength dependent)" + lin_τ + "Albedo scalar (not wavelength dependent)" + lin_ϖ + "Combined Z moment (forward)" + lin_Z⁺⁺ + "Combined Z moment (backward)" + lin_Z⁻⁺ + "Maximum dτ" + lin_dτ_max + "dτ" + lin_dτ + #"Number of doublings" + #ndoubl + "dτ (wavelength dependent)" + lin_dτ_λ + "expk" + lin_expk + #"Scattering flag" + #scatter + "Sum of optical thicknesses of all layers above the current layer" + lin_τ_sum + "Fraction of scattering caused by Rayleigh" + lin_fscattRayl + #"Elastic fraction (Cabannes) of Rayleigh (Cabannes+Raman) scattering" + #ϖ_Cabannes + #"Scattering interface type for current layer" + #scattering_interface +end + +abstract type AbstractOpticalProperties end + +# Core optical Properties COP +Base.@kwdef struct linCoreScatteringOpticalProperties{FT,FT2,FT3} <: AbstractOpticalProperties + "Absorption optical depth (scalar or wavelength dependent)" + lin_τ::FT + "Single scattering albedo" + lin_ϖ::FT2 + "Z scattering matrix (forward)" + lin_Z⁺⁺::FT3 + "Z scattering matrix (backward)" + lin_Z⁻⁺::FT3 +end + +# Core optical Properties COP with directional cross section +Base.@kwdef struct linCoreDirectionalScatteringOpticalProperties{FT,FT2,FT3,FT4} <: AbstractOpticalProperties + "Absorption optical depth (scalar or wavelength dependent)" + lin_τ::FT + "Single scattering albedo" + lin_ϖ::FT2 + "Z scattering matrix (forward)" + lin_Z⁺⁺::FT3 + "Z scattering matrix (backward)" + lin_Z⁻⁺::FT3 + "Ross kernel; cross section projection factor along µ (G ∈ [0,1], 1 for isotropic σ)" + lin_G::FT4 +end + +Base.@kwdef struct linCoreAbsorptionOpticalProperties{FT} <: AbstractOpticalProperties + "Absorption optical depth (scalar or wavelength dependent)" + lin_τ::FT +end + +# Adding Core Optical Properties, can have mixed dimensions! +function Base.:+( x::CoreScatteringOpticalProperties{xFT, xFT2, xFT3}, + y::CoreScatteringOpticalProperties{yFT, yFT2, yFT3} + dx::linCoreScatteringOpticalProperties{xFT, xFT2, xFT3}, + dy::linCoreScatteringOpticalProperties{yFT, yFT2, yFT3} + ) where {xFT, xFT2, xFT3, yFT, yFT2, yFT3} + # Predefine some arrays: + xZ⁺⁺ = x.Z⁺⁺ + xZ⁻⁺ = x.Z⁻⁺ + yZ⁺⁺ = y.Z⁺⁺ + yZ⁻⁺ = y.Z⁻⁺ + + τ = x.τ .+ y.τ + wx = x.τ .* x.ϖ + wy = y.τ .* y.ϖ + w = wx .+ wy + ϖ = w ./ τ + + #xlin_Z⁺⁺ = dx.lin_Z⁺⁺ + #xlin_Z⁻⁺ = dx.lin_Z⁻⁺ + #ylin_Z⁺⁺ = dy.lin_Z⁺⁺ + #ylin_Z⁻⁺ = dy.lin_Z⁻⁺ + + lin_τ = [dx.lin_τ, dy.lin_τ] + lin_wx = (dx.lin_τ .* (x.ϖ - ϖ) .+ x.τ .* dx.lin_ϖ)./τ + lin_wy = (dy.lin_τ .* (y.ϖ - ϖ) .+ y.τ .* dy.lin_ϖ)./τ + lin_ϖ = [lin_wx, lin_wy] + + #@show xFT, xFT2, xFT3 + all(wx .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, y.Z⁺⁺, y.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dy.lin_Z⁺⁺, dy.lin_Z⁻⁺)) : nothing + all(wy .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dx.lin_Z⁺⁺, dx.lin_Z⁻⁺)) : nothing + + n = length(w); + + wy = wy ./ w + wx = wx ./ w + wx = reshape(wx,1,1,n) + wy = reshape(wy,1,1,n) + + Z⁺⁺ = (wx .* xZ⁺⁺ .+ wy .* yZ⁺⁺) + Z⁻⁺ = (wx .* xZ⁻⁺ .+ wy .* yZ⁻⁺) + + tmpx1 = (dx.lin_τ.*x.ϖ .+ x.τ.*dx.lin_ϖ)./w + tmpy1 = (dy.lin_τ.*y.ϖ .+ y.τ.*dy.lin_ϖ)./w + tmpx2 = x.τ.*x.ϖ./w + tmpy2 = y.τ.*y.ϖ./w + xlin_Z⁺⁺ = tmpx1.*(x.Z⁺⁺.-Z⁺⁺) .+ tmpx2.*dx.dZ⁺⁺ + ylin_Z⁺⁺ = tmpy1.*(y.Z⁺⁺.-Z⁺⁺) .+ tmpy2.*dy.dZ⁺⁺ + xlin_Z⁻⁺ = tmpx1.*(x.Z⁻⁺.-Z⁻⁺) .+ tmpx2.*dx.dZ⁻⁺ + ylin_Z⁻⁺ = tmpy1.*(y.Z⁻⁺.-Z⁻⁺) .+ tmpy2.*dy.dZ⁻⁺ + lin_Z⁺⁺ = [xlin_Z⁺⁺, ylin_Z⁺⁺] + lin_Z⁻⁺ = [xlin_Z⁻⁺, ylin_Z⁻⁺] + + CoreScatteringOpticalProperties(τ, ϖ, Z⁺⁺, Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺)) +end + +# Concatenate Core Optical Properties, can have mixed dimensions! +function Base.:*( x::CoreScatteringOpticalProperties, y::CoreScatteringOpticalProperties, dx::linCoreScatteringOpticalProperties, dy::linCoreScatteringOpticalProperties ) + arr_type = array_type(architecture(x.τ)) + x = expandOpticalProperties(x,arr_type); + y = expandOpticalProperties(y,arr_type); + + dx = expandOpticalProperties(dx,arr_type); + dy = expandOpticalProperties(dy,arr_type); + + CoreScatteringOpticalProperties([x.τ; y.τ],[x.ϖ; y.ϖ],cat(x.Z⁺⁺,y.Z⁺⁺, dims=3), cat(x.Z⁻⁺,y.Z⁻⁺, dims=3) ), + linCoreScatteringOpticalProperties(cat(dx.dτ, dy.dτ, dims=2),cat(dx.dϖ, dy.dϖ, dims=2), cat(dx.dZ⁺⁺,dy.dZ⁺⁺, dims=4), cat(dx.dZ⁻⁺,dy.dZ⁻⁺, dims=4) ) +end + +function Base.:+( x::CoreScatteringOpticalProperties, y::CoreAbsorptionOpticalProperties, dx::linCoreScatteringOpticalProperties, dy::linCoreAbsorptionOpticalProperties ) + τ = x.τ .+ y.τ + wx = x.τ .* x.ϖ + ϖ = (wx) ./ τ + + lin_τ = [dx.lin_τ, dy.lin_τ] + + lin_wx = (dx.lin_τ .* (x.ϖ - ϖ) + x.τ .* dx.lin_ϖ)./τ + lin_wy = -(dy.lin_τ .* ϖ)./τ + lin_ϖ = [lin_wx, lin_wy] + + tmpx1 = (dx.lin_τ.*x.ϖ .+ x.τ.*dx.lin_ϖ)./w + tmpx2 = x.τ.*x.ϖ./w + + xlin_Z⁺⁺ = tmpx1.*(x.Z⁺⁺.-Z⁺⁺) .+ tmpx2.*dx.dZ⁺⁺ + xlin_Z⁻⁺ = tmpx1.*(x.Z⁻⁺.-Z⁻⁺) .+ tmpx2.*dx.dZ⁻⁺ + + lin_Z⁺⁺ = [xlin_Z⁺⁺, zeros(shape(Z⁺⁺))] #Check for size + lin_Z⁻⁺ = [xlin_Z⁻⁺, zeros(shape(Z⁻⁺))] + + CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺) +end + +function Base.:+( y::CoreAbsorptionOpticalProperties, x::CoreScatteringOpticalProperties, dy::linCoreAbsorptionOpticalProperties, dx::linCoreScatteringOpticalProperties ) + τ = x.τ .+ y.τ + wx = x.τ .* x.ϖ + ϖ = (wx) ./ τ + + lin_τ = [dy.lin_τ, dx.lin_τ] + + lin_wx = (dx.lin_τ .* (x.ϖ - ϖ) + x.τ .* dx.lin_ϖ)./τ + lin_wy = -(dy.lin_τ .* ϖ)./τ + lin_ϖ = [lin_wy, lin_wx] + + tmpx1 = (dx.lin_τ.*x.ϖ .+ x.τ.*dx.lin_ϖ)./w + tmpx2 = x.τ.*x.ϖ./w + + xlin_Z⁺⁺ = tmpx1.*(x.Z⁺⁺.-Z⁺⁺) .+ tmpx2.*dx.dZ⁺⁺ + xlin_Z⁻⁺ = tmpx1.*(x.Z⁻⁺.-Z⁻⁺) .+ tmpx2.*dx.dZ⁻⁺ + + lin_Z⁺⁺ = [zeros(shape(Z⁺⁺)), xlin_Z⁺⁺] #Check for size + lin_Z⁻⁺ = [zeros(shape(Z⁺⁺)), xlin_Z⁻⁺] + + CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺) +end + + +function Base.:*( x::FT, y::CoreScatteringOpticalProperties{FT}, dy::linCoreScatteringOpticalProperties{FT} ) where FT + CoreScatteringOpticalProperties(y.τ * x, y.ϖ, y.Z⁺⁺, y.Z⁻⁺, y.G), linCoreScatteringOpticalProperties(dy.lin_τ * x, dy.lin_ϖ, dy.linZ⁺⁺, dy.linZ⁻⁺, dy.linG) +end + diff --git a/src/CoreRT/rt_helper_functions.jl b/src/CoreRT/rt_helper_functions.jl index e0dff6d7..730641b1 100644 --- a/src/CoreRT/rt_helper_functions.jl +++ b/src/CoreRT/rt_helper_functions.jl @@ -209,7 +209,7 @@ make_composite_layer(RS_type::Union{RRS, RRS_plus, VS_0to1_plus, VS_1to0_plus}, default_J_matrix_ie(FT, arr_type, NSens, dims, nSpec, RS_type.n_Raman) ) # With linearization -#= + "Make an added layer, supplying all default matrices" make_added_layer(RS_type::Union{noRS, noRS_plus}, lin::vSmartMOM_lin, FT, arr_type, dims, nSpec) = AddedLayer( default_matrix(FT, arr_type, dims, nSpec), @@ -396,7 +396,7 @@ make_composite_layer(RS_type::Union{RRS, RRS_plus, VS_0to1_plus, VS_1to0_plus}, default_J_matrix_ie(FT, arr_type, lin.nparams, NSens, dims, nSpec, RS_type.n_Raman), default_J_matrix_ie(FT, arr_type, lin.nparams, NSens, dims, nSpec, RS_type.n_Raman) ) -=# + #ending linearization block "Given a ComputedAtmosphereProperties object, extract a ComputedLayerProperties object using data from the iz index of all arrays in the ComputedAtmosphereProperties" diff --git a/src/CoreRT/rt_run.jl b/src/CoreRT/rt_run.jl index 71f59a90..4ff9d84b 100644 --- a/src/CoreRT/rt_run.jl +++ b/src/CoreRT/rt_run.jl @@ -41,8 +41,6 @@ function rt_run(RS_type::AbstractRamanType, #FT = eltype(sza) # Get the float-type to use Nz = length(model.profile.p_full) # Number of vertical slices - - RS_type.bandSpecLim = [] # (1:τ_abs[iB])#zeros(Int64, iBand, 2) #Suniti: how to do this? #Suniti: make bandSpecLim a part of RS_type (including noRS) so that it can be passed into rt_kernel and elemental/doubling/interaction and postprocessing_vza without major syntax changes #put this code in model_from_parameters @@ -209,7 +207,7 @@ function rt_run(RS_type::AbstractRamanType, #end end -#= + function rt_run(RS_type::AbstractRamanType, model::vSmartMOM_Model, lin::vSmartMOM_lin, iBand) @unpack obs_alt, sza, vza, vaz = model.obs_geom # Observational geometry properties @@ -230,8 +228,6 @@ function rt_run(RS_type::AbstractRamanType, #FT = eltype(sza) # Get the float-type to use Nz = length(model.profile.p_full) # Number of vertical slices - - RS_type.bandSpecLim = [] # (1:τ_abs[iB])#zeros(Int64, iBand, 2) #Suniti: how to do this? #Suniti: make bandSpecLim a part of RS_type (including noRS) so that it can be passed into rt_kernel and elemental/doubling/interaction and postprocessing_vza without major syntax changes #put this code in model_from_parameters @@ -313,14 +309,14 @@ function rt_run(RS_type::AbstractRamanType, # Create arrays @timeit "Creating layers" added_layer = - make_added_layer(RS_type, FT_dual, arr_type, lin.nparams, dims, nSpec) + make_added_layer(RS_type, FT_dual, arr_type, lin, dims, nSpec) # Just for now, only use noRS here @timeit "Creating layers" added_layer_surface = - make_added_layer(RS_type, FT_dual, arr_type, lin.nparams, dims, nSpec) + make_added_layer(RS_type, FT_dual, arr_type, lin, dims, nSpec) @timeit "Creating layers" composite_layer = - make_composite_layer(RS_type, FT_dual, arr_type, lin.nparams, dims, nSpec) + make_composite_layer(RS_type, FT_dual, arr_type, lin, dims, nSpec) @timeit "Creating arrays" I_static = - Diagonal(arr_type(Diagonal{FT}(ones(dims[1])))); + Diagonal(arr_type(Diagonal{FT}(ones(dims[1])))); #TODO: if RS_type!=noRS, create ϖ_λ₁λ₀, i_λ₁λ₀, fscattRayl, Z⁺⁺_λ₁λ₀, Z⁻⁺_λ₁λ₀ (for input), and ieJ₀⁺, ieJ₀⁻, ieR⁺⁻, ieR⁻⁺, ieT⁻⁻, ieT⁺⁺, ier⁺⁻, ier⁻⁺, iet⁻⁻, iet⁺⁺ (for output) #getRamanSSProp(RS_type, λ, grid_in) @@ -330,95 +326,108 @@ function rt_run(RS_type::AbstractRamanType, # Loop over fourier moments for m = 0:max_m - 1 - println("Fourier Moment: ", m, "/", max_m-1) - - # Azimuthal weighting - weight = m == 0 ? FT(0.5) : FT(1.0) - # Set the Zλᵢλₒ interaction parameters for Raman (or nothing for noRS) - InelasticScattering.computeRamanZλ!(RS_type, pol_type,Array(qp_μ), m, arr_type) - # Compute the core layer optical properties: - @timeit "OpticalProps" layer_opt_props, fScattRayleigh = - constructCoreOpticalProperties(RS_type,iBand,m,model); - # Determine the scattering interface definitions: - scattering_interfaces_all, τ_sum_all = - extractEffectiveProps(layer_opt_props,quad_points); - #@show typeof(layer_opt_props) - - # Loop over vertical layers: - @showprogress 1 "Looping over layers ..." for iz = 1:Nz # Count from TOA to BOA - - # Construct the atmospheric layer - # From Rayleigh and aerosol τ, ϖ, compute overall layer τ, ϖ - # Suniti: modified to return fscattRayl as the last element of computed_atmosphere_properties - if !(typeof(RS_type) <: noRS) - RS_type.fscattRayl = expandBandScalars(RS_type, fScattRayleigh[iz]) - end + println("Fourier Moment: ", m, "/", max_m-1) + + # Azimuthal weighting + weight = m == 0 ? FT(0.5) : FT(1.0) + # Set the Zλᵢλₒ interaction parameters for Raman (or nothing for noRS) + InelasticScattering.computeRamanZλ!(RS_type, pol_type,Array(qp_μ), m, arr_type) + # Compute the core layer optical properties: + @timeit "OpticalProps" layer_opt_props, fScattRayleigh, + lin_layer_opt_props, lin_fScattRayleigh = + constructCoreOpticalProperties(RS_type, iBand, m, model, lin); + # Determine the scattering interface definitions: + scattering_interfaces_all, τ_sum_all, lin_τ_sum_all = + extractEffectiveProps(layer_opt_props, lin_layer_opt_props, quad_points); + #@show typeof(layer_opt_props) + + # Loop over vertical layers: + @showprogress 1 "Looping over layers ..." for iz = 1:Nz # Count from TOA to BOA + + # Construct the atmospheric layer + # From Rayleigh and aerosol τ, ϖ, compute overall layer τ, ϖ + # Suniti: modified to return fscattRayl as the last element of computed_atmosphere_properties + if !(typeof(RS_type) <: noRS) + RS_type.fscattRayl = expandBandScalars(RS_type, fScattRayleigh[iz]) + + nparams = length(lin_fScattRayleigh[:,1]) + speclen = length(RS_type.fscattRayl[1,:]) + RS_type.lin_fscattRayl = zeros(nparams,speclen) + for ctr = 1:nparams + RS_type.lin_fScattRayleigh[ctr,:] = + expandBandScalars(RS_type, lin_fScattRayleigh[ctr,iz]) + end + end - # Expand all layer optical properties to their full dimension: - @timeit "OpticalProps" layer_opt = - expandOpticalProperties(layer_opt_props[iz], arr_type) - - # Perform Core RT (doubling/elemental/interaction) - rt_kernel!(RS_type, pol_type, SFI, - #bandSpecLim, - added_layer, composite_layer, - layer_opt, - scattering_interfaces_all[iz], - τ_sum_all[:,iz], - m, quad_points, - I_static, - model.params.architecture, - qp_μN, iz) - end - - # Create surface matrices: - create_surface_layer!(brdf, - added_layer_surface, - SFI, m, - pol_type, - quad_points, - arr_type(τ_sum_all[:,end]), - model.params.architecture); - - #@show composite_layer.J₀⁺[iμ₀,1,1:3] - # One last interaction with surface: - @timeit "interaction" interaction!(RS_type, - #bandSpecLim, - scattering_interfaces_all[end], - SFI, - composite_layer, - added_layer_surface, - I_static) - #@show composite_layer.J₀⁺[iμ₀,1,1:3] - hdr_J₀⁻ = similar(composite_layer.J₀⁻) - # One last interaction with surface: - @timeit "interaction_HDRF" interaction_hdrf!(#RS_type, - #bandSpecLim, - #scattering_interfaces_all[end], - SFI, - composite_layer, - added_layer_surface, - m, pol_type, quad_points, - hdr_J₀⁻, bhr_uw, bhr_dw) - - # Postprocess and weight according to vza - @timeit "Postprocessing VZA" postprocessing_vza!(RS_type, - iμ₀, pol_type, - composite_layer, - vza, qp_μ, m, vaz, μ₀, - weight, nSpec, - SFI, - R, R_SFI, - T, T_SFI, - ieR_SFI, ieT_SFI) - - @timeit "Postprocessing HDRF" postprocessing_vza_hdrf!(RS_type, - iμ₀, pol_type, - hdr_J₀⁻, - vza, qp_μ, m, vaz, μ₀, - weight, nSpec, - hdr) + # Expand all layer optical properties to their full dimension: + @timeit "OpticalProps" layer_opt = + expandOpticalProperties(layer_opt_props[iz], arr_type) + @timeit "OpticalProps" lin_layer_opt = + expandOpticalProperties(lin_layer_opt_props[iz], arr_type) + + # Perform Core RT (doubling/elemental/interaction) + rt_kernel!(RS_type, pol_type, SFI, + #bandSpecLim, + added_layer, lin_added_layer, + composite_layer, lin_composite_layer, + layer_opt, + lin_layer_opt, + scattering_interfaces_all[iz], + τ_sum_all[:,iz], + lin_τ_sum_all[:,:,iz], + m, quad_points, + I_static, + model.params.architecture, + qp_μN, iz) + end + + # Create surface matrices: + create_surface_layer!(brdf, + added_layer_surface, + SFI, m, + pol_type, + quad_points, + arr_type(τ_sum_all[:,end]), + model.params.architecture); + + #@show composite_layer.J₀⁺[iμ₀,1,1:3] + # One last interaction with surface: + @timeit "interaction" interaction!(RS_type, + #bandSpecLim, + scattering_interfaces_all[end], + SFI, + composite_layer, + added_layer_surface, + I_static) + #@show composite_layer.J₀⁺[iμ₀,1,1:3] + hdr_J₀⁻ = similar(composite_layer.J₀⁻) + # One last interaction with surface: + @timeit "interaction_HDRF" interaction_hdrf!(#RS_type, + #bandSpecLim, + #scattering_interfaces_all[end], + SFI, + composite_layer, + added_layer_surface, + m, pol_type, quad_points, + hdr_J₀⁻, bhr_uw, bhr_dw) + # Postprocess and weight according to vza + @timeit "Postprocessing VZA" postprocessing_vza!(RS_type, + iμ₀, pol_type, + composite_layer, + vza, qp_μ, m, vaz, μ₀, + weight, nSpec, + SFI, + R, R_SFI, + T, T_SFI, + ieR_SFI, ieT_SFI) + + @timeit "Postprocessing HDRF" postprocessing_vza_hdrf!(RS_type, + iμ₀, pol_type, + hdr_J₀⁻, + vza, qp_μ, m, vaz, μ₀, + weight, nSpec, + hdr) end # Show timing statistics @@ -434,4 +443,3 @@ function rt_run(RS_type::AbstractRamanType, #return SFI ? (R_SFI, T_SFI, ieR_SFI, ieT_SFI) : (R, T) #end end -=# \ No newline at end of file diff --git a/src/Scattering/lin_compute_NAI2.jl b/src/Scattering/lin_compute_NAI2.jl index d57a7d9e..994d30d4 100644 --- a/src/Scattering/lin_compute_NAI2.jl +++ b/src/Scattering/lin_compute_NAI2.jl @@ -318,7 +318,8 @@ function compute_aerosol_optical_properties(model::MieModel{FDT}, FT2::Type=Floa end #@show typeof(convert.(FT2, β)), typeof(greek_coefs) # Return the packaged AerosolOptics object - return AerosolOptics(greek_coefs=greek_coefs, ω̃=FT2(ω̃), k=FT2(bulk_C_ext), fᵗ=FT2(1)), linAerosolOptics(d_greek_coefs=d_greek_coefs, dω̃=FT2(dω̃), dk=FT2(d_bulk_C_ext), dfᵗ=FT2(1)) + return AerosolOptics(greek_coefs=greek_coefs, ω̃=FT2(ω̃), k=FT2(bulk_C_ext), k_ref=FT(0), fᵗ=FT2(1)), + dAerosolOptics(d_greek_coefs=d_greek_coefs, dω̃=FT2(dω̃), dk=FT2(d_bulk_C_ext), dk_ref=zeros(FT, 4), dfᵗ=FT2(1)) else greek_coefs = GreekCoefs(α,β,γ,δ,ϵ,ζ) @@ -332,7 +333,8 @@ function compute_aerosol_optical_properties(model::MieModel{FDT}, FT2::Type=Floa dζ[i,:]) push!(d_greek_coefs, tmp_greek_coefs); end - return AerosolOptics(greek_coefs=greek_coefs, ω̃=(bulk_C_sca / bulk_C_ext), k=(bulk_C_ext), fᵗ=FT(1)), , linAerosolOptics(d_greek_coefs=d_greek_coefs, dω̃=dω̃, dk=d_bulk_C_ext, dfᵗ=FT(1)) + return AerosolOptics(greek_coefs=greek_coefs, ω̃=(bulk_C_sca / bulk_C_ext), k=(bulk_C_ext), k_ref=FT(0), fᵗ=FT(1)), + dAerosolOptics(d_greek_coefs=d_greek_coefs, dω̃=dω̃, dk=d_bulk_C_ext, dk_ref=zeros(FT, 4), dfᵗ=FT(1)) end end diff --git a/src/Scattering/lin_truncate_phase.jl b/src/Scattering/lin_truncate_phase.jl index fef989fa..052125a9 100644 --- a/src/Scattering/lin_truncate_phase.jl +++ b/src/Scattering/lin_truncate_phase.jl @@ -17,7 +17,7 @@ Returns the truncated aerosol optical properties as [`AerosolOptics`](@ref) - `mod` a [`δBGE`](@ref) struct that defines the truncation order (new length of greek parameters) and exclusion angle - `aero` a [`AerosolOptics`](@ref) set of aerosol optical properties that is to be truncated """ -function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}, lin::LinAerosolOptics{FT}; reportFit=false) where {FT} +function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}, lin::dAerosolOptics{FT}; reportFit=false) where {FT} @unpack greek_coefs, ω̃, k = aero @unpack d_greek_coefs, dω̃, dk = lin @unpack α, β, γ, δ, ϵ, ζ = greek_coefs @@ -27,10 +27,10 @@ function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}, lin::LinAerosolOpti l_tr = l_max # Obtain Gauss-Legendre quadrature points and weights for phase function: - μ, w_μ = gausslegendre(length(β)); + μ, w_μ = gausslegendre(Int((length(β)-1)/2)); # Reconstruct phase matrix elements: - scattering_matrix, P, P² = reconstruct_phase(greek_coefs, μ; returnLeg=true) + scattering_matrix, dscattering_matrix, P, P² = reconstruct_phase(greek_coefs, d_greek_coefs, μ; returnLeg=true) @unpack f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄ = scattering_matrix @unpack df₁₁, df₁₂, df₂₂, df₃₃, df₃₄, df₄₄ = dscattering_matrix @@ -71,40 +71,36 @@ function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}, lin::LinAerosolOpti if reportFit println("Errors in δ-BGE fits:") mod_y = convert.(FT, A * cl) - @show StatsBase.rmsd(mod_y, y₁₁; normalize=true) + @show StatsBase.rmsd(mod_y, b; normalize=true) end #= for dβ Ax=b, where - Aᵢⱼ = ∑ₖ w_μₖ Pᵢ(μₖ)Pⱼ(μₖ)/f₁₁²(μₖ), xᵢ=dcᵢ (as in Sanghavi & Stephens 2015), and + Aᵢⱼ = ∑ₖ w_μₖ Pᵢ(μₖ)Pⱼ(μₖ)/f₁₁²(μₖ), xᵢ=c'ᵢ (as in Sanghavi & Stephens 2015), and bᵢ = ∑ₖ w_μₖ ((2/f₁₁(μₖ))∑ₗcₗPₗ(μₖ) - 1) Pᵢ(μₖ)f₁₁'(μₖ)/(f₁₁(μₖ))^2 =# #A = zeros(l_tr, l_tr) x = zeros(l_tr) b = zeros(l_tr) + dcl = zeros(4, l_tr) - for i = 1:l_tr - b[i] = sum(((2./f₁₁).*sum(cl.*P[:,:]) .- 1)w_μ.*P[:,i].*df₁₁[ctr,:]./f₁₁) - #A[i,i] = sum(w_μ.*(P[:,i]./f₁₁).^2) - for j = i+1:l_tr - A[i,j] = sum(w_μ.*P[:,i].*P[:,j]./(f₁₁.^2)) - A[j,i] = A[i,j] + for ctr=1:4 + for i = 1:l_tr + b[i] = sum(((2. / f₁₁).*sum(cl.*(P[:,1:l_tr])', dims=1) .- 1).*w_μ.*P[:,i].*df₁₁[ctr,:]./f₁₁.^2) + end + dcl[ctr,:] = A\b # Julia backslash operator for SVD (just like Matlab); + # B in δ-BGR (β) + if reportFit + println("Errors in δ-BGE fits:") + mod_y = convert.(FT, A * dcl[ctr,:]) + @show StatsBase.rmsd(mod_y, b; normalize=true) end end - cl = A\b # Julia backslash operator for SVD (just like Matlab); - # B in δ-BGR (β) - if reportFit - println("Errors in δ-BGE fits:") - mod_y = convert.(FT, A * cl) - @show StatsBase.rmsd(mod_y, y₁₁; normalize=true) - end - - #= for γ Ax=b, where - Aᵢⱼ = ∑ₖ w_μₖ facᵢP²ᵢ(μₖ)facⱼP²ⱼ(μₖ)/f₁₂²(μₖ), xᵢ=gᵢ (as in Sanghavi & Stephens 2015), and + Aᵢⱼ = ∑ₖ w_μₖ facᵢP²ᵢ(μₖ)facⱼP²ⱼ(μₖ)/f₁₂²(μₖ), xᵢ=γᵢ (as in Sanghavi & Stephens 2015), and bᵢ = ∑ₖ w_μₖ facᵢP²ᵢ(μₖ)/f₁₂(μₖ) =# A = zeros(l_tr, l_tr) @@ -125,7 +121,28 @@ function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}, lin::LinAerosolOpti if reportFit println("Errors in δ-BGE fits:") mod_γ = convert.(FT, B * γᵗ[3:end]) - @show StatsBase.rmsd(mod_γ, y₁₂; normalize=true) + @show StatsBase.rmsd(mod_γ, b; normalize=true) + end + + #= for dγ + Ax=b, where + Aᵢⱼ = ∑ₖ w_μₖ facᵢP²ᵢ(μₖ)facⱼP²ⱼ(μₖ)/f₁₂²(μₖ), xᵢ=γ'ᵢ (as in Sanghavi & Stephens 2015), and + bᵢ = ∑ₖ w_μₖ facᵢP²ᵢ(μₖ)f₁₂'(μₖ)/f₁₂²(μₖ) {2∑ₗP²ₗ(μₖ)γₗ/f₁₂(μₖ) - 1} + =# + x = zeros(l_tr) + b = zeros(l_tr) + dγᵗ = zeros(4, l_tr) + for ctr=1:4 + for i = 3:l_tr + b[i] = fac[i]*sum(((2. / f₁₂).*sum(γᵗ.*(P²[:,1:l_tr])', dims=1) .- 1).*w_μ.*P²[:,i].*df₁₂[ctr,:]./f₁₂.^2) + end + dγᵗ[ctr,1:2] .=0 + dγᵗ[ctr,3:end] = A[3:end,3:end] \ b[3:end] # G in δ-BGE (γ) + if reportFit + println("Errors in δ-BGE fits:") + mod_γ = convert.(FT, A * dγᵗ[ctr,:]) + @show StatsBase.rmsd(mod_γ, b; normalize=true) + end end #= for ϵ @@ -151,8 +168,29 @@ function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}, lin::LinAerosolOpti if reportFit println("Errors in δ-BGE fits:") - mod_ϵ = convert.(FT, B * ϵᵗ[3:end]) - @show StatsBase.rmsd(mod_ϵ, y₃₄; normalize=true) + mod_ϵ = convert.(FT, A * ϵᵗ) + @show StatsBase.rmsd(mod_ϵ, b; normalize=true) + end + + #= for dϵ + Ax=b, where + Aᵢⱼ = ∑ₖ w_μₖ facᵢP²ᵢ(μₖ)facⱼP²ⱼ(μₖ)/f₃₄²(μₖ), xᵢ=ϵ'ᵢ (as in Sanghavi & Stephens 2015), and + bᵢ = ∑ₖ w_μₖ facᵢP²ᵢ(μₖ)f₃₄'(μₖ)/f₃₄²(μₖ) {2∑ₗP²ₗ(μₖ)ϵₗ/f₃₄(μₖ) - 1} + =# + x = zeros(l_tr) + b = zeros(l_tr) + dϵᵗ = zeros(4, l_tr) + for ctr=1:4 + for i = 3:l_tr + b[i] = fac[i]*sum(((2. / f₃₄).*sum(γᵗ.*(P²[:,1:l_tr])', dims=1) .- 1).*w_μ.*P²[:,i].*df₃₄[ctr,:]./f₃₄.^2) + end + dϵᵗ[ctr,1:2] .=0 + dϵᵗ[ctr,3:end] = A[3:end,3:end] \ b[3:end] # G in δ-BGE (γ) + if reportFit + println("Errors in δ-BGE fits:") + mod_ϵ = convert.(FT, A * dϵᵗ[ctr,:]) + @show StatsBase.rmsd(mod_ϵ, b; normalize=true) + end end # Integrate truncated function for later renormalization (here: fraction that IS still scattered): @@ -164,100 +202,25 @@ function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}, lin::LinAerosolOpti αᵗ = (α[1:l_tr] .- (β[1:l_tr] .- cl)) / c₀ # Eq. 38c, derived from β ζᵗ = (ζ[1:l_tr] .- (β[1:l_tr] .- cl)) / c₀ # Eq. 38d, derived from β + for ctr = 1:4 + #dβᵗ[ctr,1] = 0. + for i = 1:l_tr + dβᵗ[ctr,:] = (dcl[ctr,:] - βᵗ*dcl[ctr,1])/c₀ + dδᵗ[ctr,:] = (dδ[ctr,1:l_tr] .- (dβ[ctr,1:l_tr] .- dcl[ctr,:]) - δᵗ*dcl[ctr,1]) / c₀ # Eq. 38b, derived from β + dαᵗ[ctr,:] = (dα[ctr,1:l_tr] .- (dβ[ctr,1:l_tr] .- dcl[ctr,:]) - αᵗ*dcl[ctr,1]) / c₀ # Eq. 38c, derived from β + dζᵗ[ctr,:] = (dζ[ctr,1:l_tr] .- (dβ[ctr,1:l_tr] .- dcl[ctr,:]) - ζᵗ*dcl[ctr,1]) / c₀ # Eq. 38d, derived from β + end + end # Adjust scattering and extinction cross section! - greek_coefs = GreekCoefs(αᵗ, βᵗ, γᵗ, δᵗ, ϵᵗ, ζᵗ ) - - # C_sca = (ω̃ * k); - # C_scaᵗ = C_sca * c₀; - # C_ext = k - (C_sca - C_scaᵗ); - #@show typeof(ω̃), typeof(k),typeof(c₀) - # return AerosolOptics(greek_coefs = greek_coefs, ω̃=C_scaᵗ / C_ext, k=C_ext, fᵗ = 1-c₀) - return AerosolOptics(greek_coefs=greek_coefs, ω̃=ω̃, k=k, fᵗ=(FT(1) - c₀)) -end - -""" - $(FUNCTIONNAME)(mod::δBGE, aero::AerosolOptics)) - -Returns the truncated aerosol optical properties as [`AerosolOptics`](@ref) -- `mod` a [`δBGE`](@ref) struct that defines the truncation order (new length of greek parameters) and exclusion angle -- `aero` a [`AerosolOptics`](@ref) set of aerosol optical properties that is to be truncated -""" -function truncate_phase(mod::δBGE, aero::AerosolOptics{FT}; lin::LinAerosolOptics{FT}, reportFit=false) where {FT} - @unpack greek_coefs, ω̃, k = aero - @unpack d_greek_coefs, dω̃, dk = lin - @unpack α, β, γ, δ, ϵ, ζ = greek_coefs - @unpack dα, dβ, dγ, dδ, dϵ, dζ = d_greek_coefs - @unpack l_max, Δ_angle = mod - - - # Obtain Gauss-Legendre quadrature points and weights for phase function: - μ, w_μ = gausslegendre(length(β)); - μ = convert.(FT,μ) - w_μ = convert.(FT,w_μ) - # Reconstruct phase matrix elements: - scattering_matrix, dscattering_matrix, P, P² = reconstruct_phase(greek_coefs, d_greek_coefs, μ; returnLeg=true) - - @unpack f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄ = scattering_matrix - @unpack df₁₁, df₁₂, df₂₂, df₃₃, df₃₄, df₄₄ = dscattering_matrix - - # Find elements that exclude the peak (if wanted!) - iμ = findall(x -> x < cosd(Δ_angle), μ) - - # Prefactor for P2: - fac = zeros(FT,l_max); - for l = 2:l_max - 1 - fac[l + 1] = sqrt(FT(1) / ( ( l - FT(1)) * l * (l + FT(1)) * (l + FT(2)) )); - end - - # Create subsets (for Ax=y weighted least-squares fits): - y₁₁ = view(f₁₁, iμ) - y₁₂ = view(f₁₂, iμ) - y₃₄ = view(f₃₄, iμ) - A = view(P, iμ, 1:l_max) - B = fac[3:end]' .* view(P², iμ, 3:l_max) - - # Weights (also avoid division by 0) - minY = zeros(length(iμ)) .+ FT(1e-8); - W₁₁ = Diagonal(w_μ[iμ] ./ max(abs.(y₁₁), minY)); - W₁₂ = Diagonal(w_μ[iμ] ./ max(abs.(y₁₂), minY)); - W₃₄ = Diagonal(w_μ[iμ] ./ max(abs.(y₃₄), minY)); - # W₁₂ = Diagonal(w_μ[iμ]); - # W₃₄ = Diagonal(w_μ[iμ]); - # Julia backslash operator for SVD (just like Matlab); - cl = ((W₁₁ * A) \ (W₁₁ * y₁₁)) # B in δ-BGE (β) - γᵗ = similar(cl); γᵗ[1:2] .=0 - ϵᵗ = similar(cl); ϵᵗ[1:2] .=0 - γᵗ[3:end] = ((W₁₂ * B) \ (W₁₂ * y₁₂)) # G in δ-BGE (γ) - ϵᵗ[3:end] = ((W₃₄ * B) \ (W₃₄ * y₃₄)) # E in δ-BGE (ϵ) - - if reportFit - println("Errors in δ-BGE fits:") - mod_y = convert.(FT, A * cl) - mod_γ = convert.(FT, B * γᵗ[3:end]) - mod_ϵ = convert.(FT, B * ϵᵗ[3:end]) - @show StatsBase.rmsd(mod_y, y₁₁; normalize=true) - @show StatsBase.rmsd(mod_γ, y₁₂; normalize=true) - @show StatsBase.rmsd(mod_ϵ, y₃₄; normalize=true) - end - - # Integrate truncated function for later renormalization (here: fraction that IS still scattered): - c₀ = FT(cl[1]) # ( w_μ' * (P[:,1:l_max] * cl) ) / 2 - - # Compute truncated greek coefficients: - βᵗ = cl / c₀ # Eq. 38a, B in δ-BGE (β) - δᵗ = (δ[1:l_max] .- (β[1:l_max] .- cl)) / c₀ # Eq. 38b, derived from β - αᵗ = (α[1:l_max] .- (β[1:l_max] .- cl)) / c₀ # Eq. 38c, derived from β - ζᵗ = (ζ[1:l_max] .- (β[1:l_max] .- cl)) / c₀ # Eq. 38d, derived from β - - # Adjust scattering and extinction cross section! - greek_coefs = GreekCoefs(αᵗ, βᵗ, γᵗ, δᵗ, ϵᵗ, ζᵗ ) - #@show typeof(greek_coefs) - + greek_coefs = GreekCoefs(αᵗ, βᵗ, γᵗ, δᵗ, ϵᵗ, ζᵗ) + d_greek_coefs = GreekCoefs(dαᵗ, dβᵗ, dγᵗ, dδᵗ, dϵᵗ, dζᵗ) + dfᵗ = zeros(4) + dfᵗ = (FT(1) .- dcl[:,1]) # C_sca = (ω̃ * k); # C_scaᵗ = C_sca * c₀; # C_ext = k - (C_sca - C_scaᵗ); #@show typeof(ω̃), typeof(k),typeof(c₀) # return AerosolOptics(greek_coefs = greek_coefs, ω̃=C_scaᵗ / C_ext, k=C_ext, fᵗ = 1-c₀) - return AerosolOptics(greek_coefs=greek_coefs, ω̃=ω̃, k=k, fᵗ=(FT(1) - c₀)) + return AerosolOptics(greek_coefs=greek_coefs, ω̃=ω̃, k=k, k_ref=0, fᵗ=(FT(1) - c₀)), dAerosolOptics(d_greek_coefs, dω̃, dk, dk_ref=zeros(4), dfᵗ) end diff --git a/src/Scattering/lin_types.jl b/src/Scattering/lin_types.jl new file mode 100644 index 00000000..b25ea956 --- /dev/null +++ b/src/Scattering/lin_types.jl @@ -0,0 +1,92 @@ +#= + +This file contains all types that are used in the Scattering module: + +- `GreekCoefs` holds all greek coefficients +- `ScatteringMatrix` holds all computed phase function elements +- `AerosolOptics` holds all computed aerosol optical properties + +=# +#= + +Types that are needed for the output of the Fourier decomposition + +=# + +""" + struct GreekCoefs{FT} + +A struct which holds all Greek coefficient lists (over l) in one object. +See eq 16 in Sanghavi 2014 for details. + +# Fields +$(DocStringExtensions.FIELDS) +""" +struct dGreekCoefs{FT<:Union{AbstractFloat, ForwardDiff.Dual}} + "Greek matrix coefficient α, is in B[2,2]" + dα::Array{FT,1} + "Greek matrix coefficient β, is in B[1,1] (only important one for scalar!)" + dβ::Array{FT,1} + "Greek matrix coefficient γ, is in B[2,1],B[1,2]" + dγ::Array{FT,1} + "Greek matrix coefficient δ, is in B[4,4]" + dδ::Array{FT,1} + "Greek matrix coefficient ϵ, is in B[3,4] and - in B[4,3]" + dϵ::Array{FT,1} + "Greek matrix coefficient ζ, is in B[3,3]" + dζ::Array{FT,1} +end + +""" Extend Base.isapprox (≈) to compare two GreekCoefs """ +function Base.:isapprox(d_greek_coefs_a::dGreekCoefs, d_greek_coefs_b::dGreekCoefs) + field_names = fieldnames(dGreekCoefs) + return all([getproperty(dgreek_coefs_a, field) ≈ getproperty(dgreek_coefs_b, field) for field in field_names]) +end + +""" + struct ScatteringMatrix + +A struct which holds all computed phase function elements. +f₁₁ represents the phase function p for the Intensity (first Stokes Vector element) and is normalized as follows: +1/4π ∫₀²⁽ᵖⁱ⁾ dϕ ∫₋₁¹ p(μ) dμ = 1 + +# Fields +$(DocStringExtensions.FIELDS) +""" +struct dScatteringMatrix{FT} + df₁₁::FT + df₁₂::FT + df₂₂::FT + df₃₃::FT + df₃₄::FT + df₄₄::FT +end + +""" + struct AerosolOptics + +A struct which holds all computed aerosol optics + +# Fields +$(DocStringExtensions.FIELDS) +""" +Base.@kwdef struct dAerosolOptics{FT<:Union{AbstractFloat, ForwardDiff.Dual}} + "derivatives of Greek matrix w.r.t nᵣ, nᵢ, r₀ and σ₀" + d_greek_coefs::GreekCoefs + "derivatives of Single Scattering Albedo w.r.t nᵣ, nᵢ, r₀ and σ₀" + dω̃::FT + "derivatives of Extinction cross-section w.r.t nᵣ, nᵢ, r₀ and σ₀" + dk::FT + "derivatives of Extinction cross-section at reference wavelength w.r.t nᵣ, nᵢ, r₀ and σ₀" + dk_ref::FT + "derivatives of Truncation factor w.r.t nᵣ, nᵢ, r₀ and σ₀" + dfᵗ::FT + "Derivatives" + derivs = zeros(1) +end + +""" Extend Base.isapprox (≈) to compare two AerosolOptics """ +function Base.:isapprox(daerosol_optics_a::dAerosolOptics, daerosol_optics_b::dAerosolOptics) + field_names = fieldnames(dAerosolOptics) + return all([getproperty(daerosol_optics_a, field) ≈ getproperty(daerosol_optics_b, field) for field in field_names]) +end \ No newline at end of file diff --git a/src/Scattering/types.jl b/src/Scattering/types.jl index b428fbba..3f224452 100644 --- a/src/Scattering/types.jl +++ b/src/Scattering/types.jl @@ -250,6 +250,8 @@ Base.@kwdef struct AerosolOptics{FT<:Union{AbstractFloat, ForwardDiff.Dual}} ω̃::FT "Extinction cross-section" k::FT + "Extinction cross-section at reference wavelength" + k_ref::FT "Truncation factor" fᵗ::FT "Derivatives" diff --git a/test/benchmarks/prototype_inelastic_OCO2.jl b/test/benchmarks/prototype_inelastic_OCO2.jl index 3895b903..e584a21e 100644 --- a/test/benchmarks/prototype_inelastic_OCO2.jl +++ b/test/benchmarks/prototype_inelastic_OCO2.jl @@ -1,16 +1,16 @@ ## -using Revise +using InstrumentOperator +using Interpolations using Plots +using Revise using Statistics using vSmartMOM -using vSmartMOM.Architectures using vSmartMOM.Absorption -using vSmartMOM.Scattering +using vSmartMOM.Architectures using vSmartMOM.CoreRT -using vSmartMOM.SolarModel using vSmartMOM.InelasticScattering -using InstrumentOperator -using Interpolations +using vSmartMOM.Scattering +using vSmartMOM.SolarModel # Load OCO Data: # File names: diff --git a/test/benchmarks/testAerosol.jl b/test/benchmarks/testAerosol.jl index c228e5f3..86c6c59b 100644 --- a/test/benchmarks/testAerosol.jl +++ b/test/benchmarks/testAerosol.jl @@ -2,7 +2,7 @@ using CUDA device!(1) using DelimitedFiles using Distributions -using ImageFiltering +#using ImageFiltering using InstrumentOperator using Interpolations using JLD2 @@ -28,16 +28,16 @@ n_bands = length(parameters.spec_bands) n_aer = isnothing(parameters.scattering_params) ? 0 : length(parameters.scattering_params.rt_aerosols) truncation_type = vSmartMOM.Scattering.δBGE{parameters.float_type}(parameters.l_trunc, parameters.Δ_angle) -vmr = isnothing(parameters.absorption_params) ? Dict() : parameters.absorption_params.vmr -p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr = - vSmartMOM.CoreRT.compute_atmos_profile_fields(parameters.T, parameters.p, parameters.q, vmr) +#vmr = isnothing(parameters.absorption_params) ? Dict() : parameters.absorption_params.vmr +#p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr = +# vSmartMOM.CoreRT.compute_atmos_profile_fields(parameters.T, parameters.p, parameters.q, vmr) -profile = vSmartMOM.CoreRT.AtmosphericProfile(parameters.T, p_full, parameters.q, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr) +#profile = vSmartMOM.CoreRT.AtmosphericProfile(parameters.T, p_full, parameters.q, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr) # Reduce the profile to the number of target layers (if specified) -if parameters.profile_reduction_n != -1 - profile = vSmartMOM.CoreRT.reduce_profile(parameters.profile_reduction_n, profile); -end +#if parameters.profile_reduction_n != -1 +# profile = vSmartMOM.CoreRT.reduce_profile(parameters.profile_reduction_n, profile); +#end # aerosol_optics[iBand][iAer] aerosol_optics = [Array{vSmartMOM.Scattering.AerosolOptics}(undef, (n_aer)) for i=1:n_bands]; @@ -45,7 +45,7 @@ FT2 = isnothing(parameters.scattering_params) ? parameters.float_type : typeof(p #FT2 = parameters.float_type # τ_aer[iBand][iAer,iZ] -τ_aer = [zeros(FT2, n_aer, length(profile.p_full)) for i=1:n_bands]; +τ_aer = [zeros(FT2, n_aer, length(model.profile.p_full)) for i=1:n_bands]; # Loop over aerosol type for i_aer=1:n_aer @@ -101,8 +101,8 @@ for i_aer=1:n_aer τ_aer[i_band][i_aer,:] = parameters.scattering_params.rt_aerosols[i_aer].τ_ref * (aerosol_optics[i_band][i_aer].k/k_ref) * - vSmartMOM.CoreRT.getAerosolLayerOptProp(1, c_aero.p₀, c_aero.σp, profile.p_half) - @show vSmartMOM.CoreRT.getAerosolLayerOptProp(1, c_aero.p₀, c_aero.σp, profile.p_half) + vSmartMOM.CoreRT.getAerosolLayerOptProp(1, c_aero.profile.μ, c_aero.profile.σ, model.profile.p_half) + @show vSmartMOM.CoreRT.getAerosolLayerOptProp(1, c_aero.profile.μ, c_aero.profile.σ, model.profile.p_half) @show aerosol_optics[i_band][i_aer].k, k_ref @show τ_aer[i_band][i_aer,:] @show parameters.scattering_params.rt_aerosols[i_aer].τ_ref