linearize for oco - 2

src/CoreRT/Surfaces/lin_lambertian_surface.jl

@@ -63,7 +63,7 @@ function create_surface_layer!(lambertian::LambertianSurfaceScalar{FT},
                     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,:];
+                    lin_added_layer.dxj₀⁻[ctr,:,1,:] .= -(R_surf*I₀_NquadN) .* exp.(-τ_sum/μ₀)'.*lin_τ_sum[ctr,:];
                 end
             end
         end
@@ -149,7 +149,7 @@ function create_surface_layer!(lambertian::LambertianSurfaceLegendre{FT},
                     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,:];
+                    lin_added_layer.dxj₀⁻[ctr,:,1,:] .= -(R_surf*I₀_NquadN) .* ρ .* exp.(-τ_sum/μ₀)'.*lin_τ_sum[ctr,:];
                 end
             end
 

src/CoreRT/lin_postprocessing_vza.jl

@@ -0,0 +1,168 @@
+#=
+
+This file contains the function to perform azimuthal-weighting to the RT matrices after all 
+kernel calculations. 
+
+=#
+
+"Perform post-processing to azimuthally-weight RT matrices"
+function lin_postprocessing_vza!(RS_type::noRS, iμ₀, pol_type, 
+        composite_layer, 
+        lin_composite_layer,
+        vza, qp_μ, m, vaz, μ₀, weight, 
+        nSpec, nparams, SFI, R_SFI, dR_SFI, T_SFI, dT_SFI) #,ieR_SFI, ieT_SFI)
+
+    # idx of μ0 = cos(sza)
+    st_iμ0, istart0, iend0 = get_indices(iμ₀, pol_type);
+
+    # Convert these to Arrays (if CuArrays), so they can be accessed by index
+    R⁻⁺ = Array(composite_layer.R⁻⁺);
+    T⁺⁺ = Array(composite_layer.T⁺⁺);
+    J₀⁺ = Array(composite_layer.J₀⁺);
+    J₀⁻ = Array(composite_layer.J₀⁻);
+
+    dR⁻⁺ = Array(lin_composite_layer.dR⁻⁺);
+    dT⁺⁺ = Array(lin_composite_layer.dT⁺⁺);
+    dJ₀⁺ = Array(lin_composite_layer.dJ₀⁺);
+    dJ₀⁻ = Array(lin_composite_layer.dJ₀⁻);
+
+    minUp = 0.0;
+    maxUp = 0.0;
+    # Loop over all viewing zenith angles
+    for i = 1:length(vza)
+
+        # Find the nearest quadrature point idx
+        iμ = nearest_point(qp_μ, cosd(vza[i]));
+        st_iμ, istart, iend = get_indices(iμ, pol_type);
+
+        # Compute bigCS
+        cos_m_phi, sin_m_phi = (cosd(m * vaz[i]), sind(m * vaz[i]));
+        bigCS = weight * Diagonal([cos_m_phi, cos_m_phi, sin_m_phi, sin_m_phi][1:pol_type.n]);
+
+        #@show J₀⁻[istart:iend,1, 1], J₀⁺[istart:iend,1, 1]
+        # Accumulate Fourier moments after azimuthal weighting
+        for s = 1:nSpec
+
+            #if SFI
+                dR  = bigCS * J₀⁻[istart:iend,1, s]
+                dRR = dR ./ R_SFI[i,:,s]
+                if dRR[1] < minUp 
+                    minUp = dRR[1]
+                elseif dRR[1] > maxUp 
+                    maxUp = dRR[1]
+                end
+
+                R_SFI[i,:,s] .+= dR;
+                T_SFI[i,:,s] .+= bigCS * J₀⁺[istart:iend,1, s];
+
+                for ctr=1:nparams
+                    dR = bigCS * dJ₀⁻[ctr,istart:iend,1, s]
+
+                    dR_SFI[ctr, i,:,s] .+= dR;
+                    dT_SFI[ctr, i,:,s] .+= bigCS * dJ₀⁺[ctr,istart:iend,1, s];
+                end
+            #else
+            #    R[i,:,s] .+= bigCS * (R⁻⁺[istart:iend, istart0:iend0, s] / μ₀) * pol_type.I₀;
+            #    T[i,:,s] .+= bigCS * (T⁺⁺[istart:iend, istart0:iend0, s] / μ₀) * pol_type.I₀;
+            #end
+        end
+    end
+
+    if m>0
+        @info "Radiance Δ range in %: from $(minUp*100) to $(maxUp*100) for m=$(m)" 
+    end
+end
+#=
+"RAMI: Perform post-processing to azimuthally-weight hdr matrices"
+function postprocessing_vza_hdrf!(RS_type::noRS, iμ₀, pol_type, 
+        hdr_J₀⁻,  vza, qp_μ, m, vaz, μ₀, weight, 
+        nSpec, hdr)
+    
+    # idx of μ0 = cos(sza)
+    st_iμ0, istart0, iend0 = get_indices(iμ₀, pol_type);
+
+    # Convert these to Arrays (if CuArrays), so they can be accessed by index
+    hdr_J₀⁻ = Array(hdr_J₀⁻);
+    
+    minUp = 0.0;
+    maxUp = 0.0;
+    # Loop over all viewing zenith angles
+    for i = 1:length(vza)
+
+        # Find the nearest quadrature point idx
+        iμ = nearest_point(qp_μ, cosd(vza[i]));
+        st_iμ, istart, iend = get_indices(iμ, pol_type);
+        
+        # Compute bigCS
+        cos_m_phi, sin_m_phi = (cosd(m * vaz[i]), sind(m * vaz[i]));
+        bigCS = weight * Diagonal([cos_m_phi, cos_m_phi, sin_m_phi, sin_m_phi][1:pol_type.n]);
+
+        #@show J₀⁻[istart:iend,1, 1], J₀⁺[istart:iend,1, 1]
+        # Accumulate Fourier moments after azimuthal weighting
+        for s = 1:nSpec
+            dR  = bigCS * hdr_J₀⁻[istart:iend,1, s]
+            hdr[i,:,s] .+= dR;                
+        end
+    end
+end
+
+
+function postprocessing_vza!(RS_type::Union{RRS, VS_0to1_plus, VS_1to0_plus}, 
+        iμ₀, pol_type, 
+        composite_layer, 
+        lin_composite_layer,
+        vza, qp_μ, m, vaz, μ₀, weight, 
+        nSpec, nparams, 
+        SFI, 
+        R_SFI, dR_SFI, T_SFI, dT_SFI, ieR_SFI, ieT_SFI)
+    
+    # idx of μ0 = cos(sza)
+    st_iμ0, istart0, iend0 = get_indices(iμ₀, pol_type);
+
+    # Convert these to Arrays (if CuArrays), so they can be accessed by index
+    R⁻⁺ = Array(composite_layer.R⁻⁺);
+    T⁺⁺ = Array(composite_layer.T⁺⁺);
+    J₀⁺ = Array(composite_layer.J₀⁺);
+    J₀⁻ = Array(composite_layer.J₀⁻);
+
+    #ieR⁻⁺ = Array(composite_layer.ieR⁻⁺);
+    #ieT⁺⁺ = Array(composite_layer.ieT⁺⁺);
+    ieJ₀⁺ = Array(composite_layer.ieJ₀⁺);
+    ieJ₀⁻ = Array(composite_layer.ieJ₀⁻);
+    # Loop over all viewing zenith angles
+    
+
+    for i = 1:length(vza)
+
+        # Find the nearest quadrature point idx
+        iμ = nearest_point(qp_μ, cosd(vza[i]));
+        st_iμ, istart, iend = get_indices(iμ, pol_type);
+        
+        # Compute bigCS
+        cos_m_phi, sin_m_phi = (cosd(m * vaz[i]), sind(m * vaz[i]));
+        bigCS = weight * Diagonal([cos_m_phi, cos_m_phi, sin_m_phi, sin_m_phi][1:pol_type.n]);
+
+        # Accumulate Fourier moments after azimuthal weighting
+        for s = 1:nSpec
+            
+            if SFI
+                R_SFI[i,:,s] += bigCS * J₀⁻[istart:iend,1, s];
+                T_SFI[i,:,s] += bigCS * J₀⁺[istart:iend,1, s];
+                #@show i, s, R_SFI[i,:,s]
+                #@show i, s, ieR_SFI[i,:,s]
+                for t =1:size(ieJ₀⁺,4)# in eachindex ieJ₀⁺[1,1,1,:]
+                    ieR_SFI[i,:,s] += bigCS * ieJ₀⁻[istart:iend,1, s, t];
+                    ieT_SFI[i,:,s] += bigCS * ieJ₀⁺[istart:iend,1, s, t];
+                    #@show i, s, t, ieR_SFI[i,:,s]
+                end
+                #@show i, s, ieR_SFI[i,:,s]
+            else
+                R[i,:,s] += bigCS * (R⁻⁺[istart:iend, istart0:iend0, s] / μ₀) * pol_type.I₀;
+                T[i,:,s] += bigCS * (T⁺⁺[istart:iend, istart0:iend0, s] / μ₀) * pol_type.I₀;
+                #no modification for Raman because SFI will be the only mode used for inelastic computations
+            end
+            
+        end
+    end
+end
+=#

src/CoreRT/rt_run.jl

@@ -340,7 +340,8 @@ function rt_run(RS_type::AbstractRamanType,
         scattering_interfaces_all, τ_sum_all, lin_τ_sum_all = 
             extractEffectiveProps(layer_opt_props, lin_layer_opt_props, quad_points);
         #@show typeof(layer_opt_props)
-
+        nparams = length(lin_fScattRayleigh[:,1])
+        speclen = length(RS_type.fscattRayl[1,:])
         # Loop over vertical layers: 
         @showprogress 1 "Looping over layers ..." for iz = 1:Nz  # Count from TOA to BOA
 
@@ -350,8 +351,7 @@ function rt_run(RS_type::AbstractRamanType,
             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,:] = 
@@ -384,10 +384,13 @@ function rt_run(RS_type::AbstractRamanType,
         # Create surface matrices:
         create_surface_layer!(brdf, 
                     added_layer_surface, 
+                    lin_added_layer_surface,
                     SFI, m, 
                     pol_type, 
                     quad_points, 
-                    arr_type(τ_sum_all[:,end]), 
+                    arr_type(τ_sum_all[:,end]),
+                    arr_type(lin_τ_sum_all[:,end]), 
+                    Nparams, i_surf,
                     model.params.architecture);
 
         #@show composite_layer.J₀⁺[iμ₀,1,1:3]
@@ -397,37 +400,42 @@ function rt_run(RS_type::AbstractRamanType,
                             scattering_interfaces_all[end], 
                             SFI, 
                             composite_layer, 
-                            added_layer_surface, 
+                            lin_composite_layer,
+                            added_layer_surface,
+                            lin_added_layer_surface,
+                            nparams, 
                             I_static)
         #@show composite_layer.J₀⁺[iμ₀,1,1:3]
-        hdr_J₀⁻ = similar(composite_layer.J₀⁻)
+        #hdr_J₀⁻ = similar(composite_layer.J₀⁻)
         # One last interaction with surface:
-        @timeit "interaction_HDRF" interaction_hdrf!(#RS_type,
+        #=@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, 
+        @timeit "Postprocessing VZA" postprocessing_vza!(#RS_type, 
                     iμ₀, pol_type, 
                     composite_layer, 
+                    lin_composite_layer,
                     vza, qp_μ, m, vaz, μ₀, 
-                    weight, nSpec, 
+                    weight, nSpec, nparams,
                     SFI, 
-                    R, R_SFI, 
-                    T, T_SFI,
-                    ieR_SFI, ieT_SFI)
-
+                    R_SFI, dR_SFI, 
+                    T_SFI, dT_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
@@ -438,7 +446,7 @@ function rt_run(RS_type::AbstractRamanType,
     #if RAMI
     #@show size(hdr), size(bhr_dw)
     #hdr = hdr[:,1,:] ./ bhr_dw[1,:]
-    return SFI ? (R_SFI, T_SFI, ieR_SFI, ieT_SFI, hdr, bhr_uw[1,:], bhr_dw[1,:]) : (R, T)
+    return R_SFI, T_SFI, dR_SFI, dT_SFI
     #else
     #return SFI ? (R_SFI, T_SFI, ieR_SFI, ieT_SFI) : (R, T)
     #end