text:
l_min(SHModes::SphericalHarmonicModes.SHModeRange) = minimum(l_range(SHModes))
l_max(SHModes::SphericalHarmonicModes.SHModeRange) = maximum(l_range(SHModes))

m_min(SHModes::SphericalHarmonicModes.SHModeRange) = minimum(m_range(SHModes))
m_max(SHModes::SphericalHarmonicModes.SHModeRange) = maximum(m_range(SHModes))

########################################################################################
# Utility function
########################################################################################

function getkernelmodes(; kwargs...)
	s_min = get(kwargs, :s_min, 0) :: Int
	s_max = get(kwargs, :s_max, s_min) :: Int
	t_min = get(kwargs, :t_min, 0) :: Int
	t_max = get(kwargs, :t_max, t_min) :: Int

	m = get(kwargs, :SHModes, LM(s_min:s_max, t_min:t_max))
	LM{UnitRange{Int}, UnitRange{Int}}(m)
end

function unpackGfnparams(p_Gsrc, r_src, p_Gobs1, r_obs1, p_Gobs2, r_obs2)
	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src)
	p_Gobs1 = read_all_parameters(p_Gobs1, r_src = r_obs1)
	p_Gobs2 = read_all_parameters(p_Gobs2, r_src = r_obs2)

	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	Gfn_path_obs1, NGfn_files_obs1 =  p_Gobs1.path, p_Gobs1.num_procs
	Gfn_path_obs2, NGfn_files_obs2 = p_Gobs2.path, p_Gobs2.num_procs

	(; Gfn_path_src, NGfn_files_src,
		Gfn_path_obs1, NGfn_files_obs1,
		Gfn_path_obs2, NGfn_files_obs2,
		p_Gsrc, p_Gobs1, p_Gobs2)
end

VSHType(::los_radial) = SH()
VSHType(::los_earth) = GSH()

function los_projected_spheroidal_biposh_flippoints(xobs1::SphericalPoint, xobs2::SphericalPoint, los, SHModes, jₒjₛ_allmodes)
	l1, l2 = line_of_sight_covariant.((xobs1, xobs2), los)
	_Y12, _Y21 = biposh_flippoints(VSHType(los), Point2D(xobs1)..., Point2D(xobs2)..., SHModes, jₒjₛ_allmodes)
	Y12 = biposh_spheroidal(los_projected_biposh!(los, _Y12, l1, l2))
	Y21 = biposh_spheroidal(los_projected_biposh!(los, _Y21, l2, l1))
	return Y12, Y21
end

function trim_m(Y::SHArray{<:Any,1, <:Vector, <:Tuple{LM}}, SHModes)
	Yp = parent(Y)
	modes = only(SphericalHarmonicArrays.modes(Y))
	modes_new = LM(l_range(modes), intersect(m_range(SHModes), m_range(modes)))
	inds = [ind for (ind,(l,m)) in enumerate(modes) if m in m_range(SHModes, l)]
	Yp = Yp[inds]
	SHArray(Yp, modes_new)
end
function trim_m(Y, SHModes)
	SHArray([trim_m(Y[ind], SHModes) for ind in eachindex(Y)],
		SphericalHarmonicArrays.modes(Y))
end

function rotate_biposh!(Y1′2′_j1j2, Y12_j1j2, M, D)
	modes = only(SphericalHarmonicArrays.modes(Y1′2′_j1j2))
	modes_allm = only(SphericalHarmonicArrays.modes(Y12_j1j2))
	for l in l_range(modes)
		Dl = D[l]
		for m in m_range(modes, l)
			Y1′2′_j1j2[(l,m)] = M * sum(Dl[m′,m] * Y12_j1j2[(l,m′)] for m′ in m_range(modes_allm, l))
		end
	end
	return Y1′2′_j1j2
end
function rotate_biposh(Y12, M, D, SHModes)
	Y1′2′ = trim_m(Y12, SHModes)
	for ind in eachindex(Y12)
		Y12_j1j2 = Y12[ind]
		Y1′2′_j1j2 = Y1′2′[ind]
		rotate_biposh!(Y1′2′_j1j2, Y12_j1j2, M, D)
	end
	return Y1′2′
end

_basisrotationmatrix(::los_earth, M1, M2) = Diagonal(LinearAlgebra.kron(M1, M2))
_basisrotationmatrix(::los_radial, M1, M2) = 1

function los_projected_spheroidal_biposh_flippoints((xobs1, xobs1′)::NTuple{2,SphericalPoint},
	(xobs2, xobs2′)::NTuple{2,SphericalPoint}, los, SHModes, jₒjₛ_allmodes)

	SHModes_allm = LM(l_range(SHModes))

	_Y12_allm, _Y21_allm = biposh_flippoints(VSHType(los), Point2D(xobs1)..., Point2D(xobs2)..., SHModes_allm, jₒjₛ_allmodes)
	_Y12 = trim_m(_Y12_allm, SHModes)
	_Y21 = trim_m(_Y21_allm, SHModes)
	Y12 = los_projected_biposh_spheroidal(_Y12, xobs1, xobs2, los)
	Y21 = los_projected_biposh_spheroidal(_Y21, xobs2, xobs1, los)
	R, M11′, M22′ = rotation_points_12(HelicityCovariant(), xobs1, xobs2, xobs1′, xobs2′)
	M11′22′ = _basisrotationmatrix(los, M11′, M22′)
	M22′11′ = _basisrotationmatrix(los, M22′, M11′)
	Rinv = inv(R)
	α, β, γ = Rinv.theta1, Rinv.theta2, Rinv.theta3
	lmax = maximum(l_range(SHModes))
	D = OffsetArray([OffsetArray(wignerD(j, α, β, γ), -j:j, -j:j) for j in 0:lmax], 0:lmax)
	_Y1′2′ = rotate_biposh(_Y12_allm, M11′22′, D, SHModes)
	_Y2′1′ = rotate_biposh(_Y21_allm, M22′11′, D, SHModes)
	Y1′2′ = los_projected_biposh_spheroidal(_Y1′2′, xobs1′, xobs2′, los)
	Y2′1′ = los_projected_biposh_spheroidal(_Y2′1′, xobs2′, xobs1′, los)
	return Y12, Y21, Y1′2′, Y2′1′
end

########################################################################################
# Validation for uniform rotation
########################################################################################

function populatekernelvalidation!(::Flow, ::los_radial, K::AbstractVector{<:Real},
	(j, ω), ∂ϕ₂Pⱼ_cosχ, conjhω, H¹₁jj_r₁r₂, H¹₁jj_r₂r₁, dω)

	pre = 2*√(3/π) * dω/2π * ω^3 * Powspec(ω) * ∂ϕ₂Pⱼ_cosχ
	for r_ind in UnitRange(axes(K, 1))
		K[r_ind] += pre * imag( conjhω * ( H¹₁jj_r₁r₂[r_ind] + conj(H¹₁jj_r₂r₁[r_ind]) ) )
	end
end

function populatekernelvalidation!(::Flow, ::los_earth, K::AbstractVector{<:Real},
	(j, ω), Y12j, conjhω, H¹₁jj_r₁r₂, H¹₁jj_r₂r₁, dω)

	pre = 2 * dω/2π * ω^3 * Powspec(ω) * √((j*(j+1))/π)

	for r_ind in UnitRange(axes(K, 1))
		for α₂ in 0:1, α₁ in 0:1
			Pʲʲ₁₀_α₁α₂ = pre * imag(Y12j[α₁, α₂])
			iszero(Pʲʲ₁₀_α₁α₂) && continue

			K[r_ind] += Pʲʲ₁₀_α₁α₂ * imag(conjhω *
				( H¹₁jj_r₁r₂[α₁, α₂, r_ind] + conj(H¹₁jj_r₂r₁[α₂, α₁, r_ind] ) ) )
		end
	end
end

function populatekernelvalidation!(::Flow, ::los_radial, K::AbstractMatrix{<:Real},
	(j, ω), ∂ϕ₂Pⱼ_cosχ::AbstractVector{<:Real},
	conjhω::AbstractVector{ComplexF64}, H¹₁jj_r₁r₂, H¹₁jj_r₂r₁, dω)

	pre = 2*√(3/4π) * dω/2π * ω^3 * Powspec(ω)
	for n2ind in axes(K, 2)
		pre2 = pre * ∂ϕ₂Pⱼ_cosχ[n2ind]
		conjhω_n2 = conjhω[n2ind]
		for r_ind in axes(K, 1)
			K[r_ind, n2ind] +=  pre2 * 2imag( conjhω_n2 * ( H¹₁jj_r₁r₂[r_ind] + conj(H¹₁jj_r₂r₁[r_ind]) ) )
		end
	end
end

function populatekernelvalidation!(::Flow, ::los_earth, K::AbstractMatrix{<:Real},
	(j, ω), Y12j::AbstractVector, conjhω, H¹₁jj_r₁r₂, H¹₁jj_r₂r₁, dω)

	pre = 2 * dω/2π * ω^3 * Powspec(ω) * √((j*(j+1))/π)

	for n2ind in axes(K, 2)
		Y12j_n₂ = Y12j[n2ind]
		conjhω_n₂ = conjhω[n2ind]

		for r_ind in UnitRange(axes(K, 1))
			for α₂ in 0:1, α₁ in 0:1
				Pʲʲ₁₀_α₁α₂ = pre * imag(Y12j_n₂[α₁, α₂])
				iszero(Pʲʲ₁₀_α₁α₂) && continue

				K[r_ind, n2ind] += Pʲʲ₁₀_α₁α₂ * imag(conjhω_n₂ *
					(H¹₁jj_r₁r₂[α₁, α₂, r_ind] + conj(H¹₁jj_r₂r₁[α₂, α₁, r_ind] ) ) )
			end
		end
	end
end

function kernel_ℑu⁺₁₀_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction, hω,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	@unpack p_Gsrc, p_Gobs1, p_Gobs2, Gfn_path_src, NGfn_files_src,
	Gfn_path_obs1, NGfn_files_obs1, Gfn_path_obs2, NGfn_files_obs2 =
		unpackGfnparams(p_Gsrc, r_src, p_Gobs1, radius(xobs1), p_Gobs2, radius(xobs2))

	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2 =
			Gfn_fits_files.((Gfn_path_src, Gfn_path_obs1, Gfn_path_obs2),
				((ℓ_arr, 1:Nν_Gfn),), (ℓ_ωind_iter_on_proc,),
				(NGfn_files_src, NGfn_files_obs1, NGfn_files_obs2))

	K = zeros(nr)

	arrs = allocatearrays(Flow(), los, r₁_ind == r₂_ind)
	@unpack Gsrc, Gobs1, Gobs2, G¹₁jj_r₁, G¹₁jj_r₂, H¹₁jj_r₁r₂, H¹₁jj_r₂r₁ = arrs

	Y12 = los_projected_biposh_spheroidal(computeY₁₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)

	for (j, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		conjhω = conj(hω[ω_ind])

		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (j, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		Gγr_r₁_rsrc = αrcomp(Gsrc, r₁_ind, los)
		Gγr_r₂_rsrc = αrcomp(Gsrc, r₂_ind, los)

		# Green function about receiver location
		read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
			(ℓ_arr, 1:Nν_Gfn), (j, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)

		radial_fn_uniform_rotation_firstborn!(G¹₁jj_r₁, Gsrc, Gobs1, j, los)
		Hjₒjₛω!(H¹₁jj_r₁r₂, G¹₁jj_r₁, Gγr_r₂_rsrc)

		if r₁_ind != r₂_ind
			# Green function about receiver location
			read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
			(ℓ_arr, 1:Nν_Gfn), (j, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)

			radial_fn_uniform_rotation_firstborn!(G¹₁jj_r₂, Gsrc, Gobs2, j, los)
			Hjₒjₛω!(H¹₁jj_r₂r₁, G¹₁jj_r₂, Gγr_r₁_rsrc)
		end

		populatekernelvalidation!(Flow(), los, K, (j, ω),
			Y12[j], conjhω, H¹₁jj_r₁r₂, H¹₁jj_r₂r₁, dω)
	end

	foreach(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))
	return @. (ρ/r)*K
end

function kernel_ℑu⁺₁₀_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	nobs1::Point2D, nobs2_arr::Vector{<:Point2D}, los::los_direction, hω,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = p_Gobs1,
	r_src = r_src_default)

	hω = permutedims(hω) #  convert to (n2, ω)

	r_obs_ind = radial_grid_index(nobs1)

	@unpack p_Gsrc, p_Gobs1, Gfn_path_src, NGfn_files_src =
		unpackGfnparams(p_Gsrc, r_src, p_Gobs1, r_obs_default, p_Gobs2, r_obs_default)

	Gfn_path_obs, NGfn_files_obs = p_Gobs1.path, p_Gobs1.num_procs

	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	Gfn_fits_files_src, Gfn_fits_files_obs =
			Gfn_fits_files.((Gfn_path_src, Gfn_path_obs),
				((ℓ_arr, 1:Nν_Gfn),), (ℓ_ωind_iter_on_proc,),
				(NGfn_files_src, NGfn_files_obs))

	K = zeros(nr, length(nobs2_arr))

	arrs = allocatearrays(Flow(), los, true)
	@unpack Gsrc, H¹₁jj_r₁r₂, G¹₁jj_r₁, Gobs1 = arrs

	Y12 = los_projected_biposh_spheroidal(computeY₁₀, nobs1, nobs2_arr, los, ℓ_ωind_iter_on_proc)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		conjhω = conj(@view hω[:, ω_ind])

		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		G_r₂_rsrc = αrcomp(Gsrc, r_obs_ind, los)

		# Green function about receiver location
		read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs, 1:2, srcindFITS(los), :, 1)

		radial_fn_uniform_rotation_firstborn!(G¹₁jj_r₁, Gsrc, Gobs1, ℓ, los)
		Hjₒjₛω!(H¹₁jj_r₁r₂, G¹₁jj_r₁, G_r₂_rsrc)

		populatekernelvalidation!(Flow(), los, K, (ℓ, ω),
			view(Y12, :, ℓ), conjhω, H¹₁jj_r₁r₂, H¹₁jj_r₁r₂, dω)
	end

	map(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs))
	return @. (ρ/r)*K
end

# Traveltimes

kernelfilenameℑu⁺₁₀(::TravelTimes, ::los_radial) = "K_δτ_ℑu⁺₁₀.fits"
kernelfilenameℑu⁺₁₀(::TravelTimes, ::los_earth) = "K_δτ_ℑu⁺₁₀_los.fits"
kernelfilenameℑu⁺₁₀(::Amplitudes, ::los_radial) = "K_A_ℑu⁺₁₀.fits"
kernelfilenameℑu⁺₁₀(::Amplitudes, ::los_earth) = "K_A_ℑu⁺₁₀_los.fits"

function kernel_ℑu⁺₁₀(comm, m::SeismicMeasurement, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	r_src, = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs...)
	@unpack ν_arr, ℓ_arr = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)

	# Compute the appropriate window function depending on the parameter of interest
	hω_arr = get(kwargs, :hω) do
		h = hω(comm, m, xobs1, xobs2, los; kwargs..., print_timings = false)
		_broadcast(h, 0, comm)
	end

	# Use the window function to compute the kernel
	K = pmapsum(comm, kernel_ℑu⁺₁₀_partial, iters,
		xobs1, xobs2, los, hω_arr, p_Gsrc, p_Gobs1, p_Gobs2, r_src)

	K === nothing && return nothing

	if get(kwargs, :save, true)
		filename = joinpath(SCRATCH_kerneldir[], kernelfilenameℑu⁺₁₀(m, los))
		FITS(filename,"w") do f
			write(f, reinterpret_as_float(K))
		end
	end

	return K
end

########################################################################################
# Kernels for flow velocity
########################################################################################

@inline function mulprefactor!(K::AbstractArray{<:Real}, vindskip = 0)
	@turbo @. K *= ρ
	@assert axes(K,1) == eachindex(r) "first axis of K does not match that of r"
	for trail_ind in CartesianIndices(UnitRange.(axes(K)[3:end])), vind in UnitRange(axes(K, 2))
		vind == vindskip && continue # radial component
		for r_ind in eachindex(r)
			K[r_ind, vind, trail_ind] /= r[r_ind]
		end
	end
end
@inline function mulprefactor!(K::StructArray{<:Complex}, vindskip = 0)
	mulprefactor!(K.re, vindskip)
	mulprefactor!(K.im, vindskip)
end

function populatekernel!(::Flow, ::los_radial, K::AbstractArray{<:Complex,3},
	SHModes, Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
	(l, m), (jₛ, jₒ, ω),
	twoimagconjhωHγℓjₒjₛ_r₁r₂,
	twoimagconjhωconjHγℓjₒjₛ_r₂r₁,
	pre⁰, pre⁺, phase)

	mode_ind = modeindex(firstshmodes(Yjₒjₛ_lm_n₁n₂), (l, m))

	Pʲₒʲₛₗₘ_₀₀_n₁n₂ = Yjₒjₛ_lm_n₁n₂[mode_ind]
	Pʲₒʲₛₗₘ_₀₀_n₂n₁ = Yjₒjₛ_lm_n₂n₁[mode_ind]

	iszero(Pʲₒʲₛₗₘ_₀₀_n₁n₂) && iszero(Pʲₒʲₛₗₘ_₀₀_n₂n₁) && return

	anyzeromom = l == 0 || jₛ == 0
	evenmom = iseven(jₒ + jₛ + l)

	zeroindT, plusindT = axes(twoimagconjhωconjHγℓjₒjₛ_r₂r₁, 2)
	minusindK, zeroindK, plusindK = axes(K, 1)
	lm_ind_K = modeindex(SHModes, (l, m))

	if evenmom
		@turbo for r_ind in eachindex(r)
			TP₁₂⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₁n₂ * twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, zeroindT]
			TP₂₁⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₂n₁ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, zeroindT]

			K[zeroindK, r_ind, lm_ind_K] += TP₂₁⁰ - TP₁₂⁰
		end
	end
	anyzeromom && return K
	@turbo for r_ind in eachindex(r)
		TP₁₂¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₁n₂ * twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, plusindT]
		TP₂₁¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₂n₁ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, plusindT]

		K[plusindK, r_ind, lm_ind_K] += TP₂₁¹ - TP₁₂¹
		K[minusindK, r_ind, lm_ind_K] += phase * (TP₂₁¹ - TP₁₂¹)
	end

	return K
end
function populatekernel!(::Flow, ::los_radial, K::AbstractArray{<:Complex,4},
	SHModes, (Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₁′n₂′)::Tuple{Any, Any},
	(Yjₒjₛ_lm_n₂n₁, Yjₒjₛ_lm_n₂′n₁′)::Tuple{Any, Any},
	(l, m), (jₛ, jₒ, ω),
	twoimagconjhωHγℓjₒjₛ_r₁r₂,
	twoimagconjhωconjHγℓjₒjₛ_r₂r₁,
	pre⁰, pre⁺, phase)

	mode_ind = modeindex(firstshmodes(Yjₒjₛ_lm_n₁n₂), (l, m))

	Pʲₒʲₛₗₘ_₀₀_n₁n₂ = Yjₒjₛ_lm_n₁n₂[mode_ind]
	Pʲₒʲₛₗₘ_₀₀_n₂n₁ = Yjₒjₛ_lm_n₂n₁[mode_ind]
	Pʲₒʲₛₗₘ_₀₀_n₁′n₂′ = Yjₒjₛ_lm_n₁′n₂′[mode_ind]
	Pʲₒʲₛₗₘ_₀₀_n₂′n₁′ = Yjₒjₛ_lm_n₂′n₁′[mode_ind]

	iszero(Pʲₒʲₛₗₘ_₀₀_n₁n₂) && iszero(Pʲₒʲₛₗₘ_₀₀_n₂n₁) &&
	iszero(Pʲₒʲₛₗₘ_₀₀_n₁′n₂′) && iszero(Pʲₒʲₛₗₘ_₀₀_n₂′n₁′) && return

	anyzeromom = l == 0 || jₛ == 0
	evenmom = iseven(jₒ + jₛ + l)

	zeroindT, plusindT = axes(twoimagconjhωconjHγℓjₒjₛ_r₂r₁, 3)
	minusindK, zeroindK, plusindK = axes(K, 2)
	lm_ind_K = modeindex(SHModes, (l, m))

	if evenmom
		@turbo for r_ind in eachindex(r)
			TP₁₂⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₁n₂ * twoimagconjhωHγℓjₒjₛ_r₁r₂[1, r_ind, zeroindT]
			TP₂₁⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₂n₁ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[1, r_ind, zeroindT]
			TP₁′₂′⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₁′n₂′ * twoimagconjhωHγℓjₒjₛ_r₁r₂[2, r_ind, zeroindT]
			TP₂′₁′⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₂′n₁′ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[2, r_ind, zeroindT]

			K[1, zeroindK, r_ind, lm_ind_K] += TP₂₁⁰ - TP₁₂⁰
			K[2, zeroindK, r_ind, lm_ind_K] += TP₂′₁′⁰ - TP₁′₂′⁰
		end
	end
	anyzeromom && return K
	@turbo for r_ind in eachindex(r)
		TP₁₂¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₁n₂ * twoimagconjhωHγℓjₒjₛ_r₁r₂[1, r_ind, plusindT]
		TP₂₁¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₂n₁ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[1, r_ind, plusindT]
		TP₁′₂′¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₁′n₂′ * twoimagconjhωHγℓjₒjₛ_r₁r₂[2, r_ind, plusindT]
		TP₂′₁′¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₂′n₁′ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[2, r_ind, plusindT]

		K[1, plusindK, r_ind, lm_ind_K] += TP₂₁¹ - TP₁₂¹
		K[2, plusindK, r_ind, lm_ind_K] += TP₂′₁′¹ - TP₁′₂′¹
		K[1, minusindK, r_ind, lm_ind_K] += phase * (TP₂₁¹ - TP₁₂¹)
		K[2, minusindK, r_ind, lm_ind_K] += phase * (TP₂′₁′¹ - TP₁′₂′¹)
	end

	return K
end
function populatekernel!(::Flow, ::los_earth, K::AbstractArray{<:Complex,3},
	SHModes, Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
	(l, m), (jₛ, jₒ, ω),
	twoimagconjhωHγℓjₒjₛ_r₁r₂,
	twoimagconjhωconjHγℓjₒjₛ_r₂r₁,
	pre⁰, pre⁺, phase)

	mode_ind = modeindex(firstshmodes(Yjₒjₛ_lm_n₁n₂), (l, m))
	lm_ind_K = modeindex(SHModes, (l, m))

	_Yjₒjₛ_lm_n₁n₂ = parent(Yjₒjₛ_lm_n₁n₂[mode_ind])
	_Yjₒjₛ_lm_n₂n₁ = parent(Yjₒjₛ_lm_n₂n₁[mode_ind])

	minusindK, zeroindK, plusindK = axes(K, 1)
	zeroindT, plusindT = axes(twoimagconjhωconjHγℓjₒjₛ_r₂r₁, 2)

	anyzeromom = l == 0 || jₛ == 0
	evenmom = iseven(jₒ + jₛ + l)
	if evenmom
		@turbo for r_ind in eachindex(r)
			TP₁₂⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₁n₂ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, zeroindT])
			TP₂₁⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₂n₁ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, zeroindT])

			K[zeroindK, r_ind, lm_ind_K] += TP₂₁⁰ - TP₁₂⁰
		end
	end
	anyzeromom && return K
	@turbo for r_ind in eachindex(r)
		TP₁₂¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₁n₂ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, plusindT])
		TP₂₁¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₂n₁ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, plusindT])

		K[plusindK, r_ind, lm_ind_K] += TP₂₁¹ - TP₁₂¹
		K[minusindK, r_ind, lm_ind_K] += phase * (TP₂₁¹ - TP₁₂¹)
	end

	return K
end
function populatekernel!(::Flow, ::los_earth, K::AbstractArray{<:Complex,4},
	SHModes,
	(Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₁′n₂′)::NTuple{2,Any},
	(Yjₒjₛ_lm_n₂n₁, Yjₒjₛ_lm_n₂′n₁′)::NTuple{2,Any},
	(l, m), (jₛ, jₒ, ω),
	twoimagconjhωHγℓjₒjₛ_r₁r₂,
	twoimagconjhωconjHγℓjₒjₛ_r₂r₁,
	pre⁰, pre⁺, phase)

	mode_ind = modeindex(firstshmodes(Yjₒjₛ_lm_n₁n₂), (l, m))
	lm_ind_K = modeindex(SHModes, (l, m))

	_Yjₒjₛ_lm_n₁n₂ = parent(Yjₒjₛ_lm_n₁n₂[mode_ind])
	_Yjₒjₛ_lm_n₂n₁ = parent(Yjₒjₛ_lm_n₂n₁[mode_ind])
	_Yjₒjₛ_lm_n₁′n₂′ = parent(Yjₒjₛ_lm_n₁′n₂′[mode_ind])
	_Yjₒjₛ_lm_n₂′n₁′ = parent(Yjₒjₛ_lm_n₂′n₁′[mode_ind])

	minusindK, zeroindK, plusindK = axes(K, 2)
	zeroindT, plusindT = axes(twoimagconjhωconjHγℓjₒjₛ_r₂r₁, 3)

	anyzeromom = l == 0 || jₛ == 0
	evenmom = iseven(jₒ + jₛ + l)
	if evenmom
		@turbo for r_ind in eachindex(r)
			TP₁₂⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₁n₂ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[1, r_ind, zeroindT])
			TP₂₁⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₂n₁ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[1, r_ind, zeroindT])
			TP₁′₂′⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₁′n₂′ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[2, r_ind, zeroindT])
			TP₂′₁′⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₂′n₁′ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[2, r_ind, zeroindT])

			K[1, zeroindK, r_ind, lm_ind_K] += TP₂₁⁰ - TP₁₂⁰
			K[2, zeroindK, r_ind, lm_ind_K] += TP₂′₁′⁰ - TP₁′₂′⁰
		end
	end
	anyzeromom && return K
	@turbo for r_ind in eachindex(r)
		TP₁₂¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₁n₂ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[1, r_ind, plusindT])
		TP₂₁¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₂n₁ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[1, r_ind, plusindT])
		TP₁′₂′¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₁′n₂′ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[2, r_ind, plusindT])
		TP₂′₁′¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₂′n₁′ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[2, r_ind, plusindT])

		K[1, plusindK, r_ind, lm_ind_K] += TP₂₁¹ - TP₁₂¹
		K[2, plusindK, r_ind, lm_ind_K] += TP₂′₁′¹ - TP₁′₂′¹
		K[1, minusindK, r_ind, lm_ind_K] += phase * (TP₂₁¹ - TP₁₂¹)
		K[2, minusindK, r_ind, lm_ind_K] += phase * (TP₂′₁′¹ - TP₁′₂′¹)
	end

	return K
end

function populatekernelrθϕl0!(::Flow, ::los_radial, K::AbstractArray{<:Real, 3},
	SHModes, Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
	(l,_), (jₛ, jₒ, ω), twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂,
	pre⁰, pre⁺, phase)

	lm_ind_K = modeindex(SHModes, (l, 0))
	lm_ind = modeindex(firstshmodes(Yjₒjₛ_lm_n₁n₂), (l, 0))
	Pʲₒʲₛₗₘ_₀₀_n₁n₂ = Yjₒjₛ_lm_n₁n₂[lm_ind]
	Pʲₒʲₛₗₘ_₀₀_n₂n₁ = Yjₒjₛ_lm_n₂n₁[lm_ind]

	rindK, θindK, ϕindK = axes(K, 1)
	zeroindT, plusindT = axes(twoimagconjhωconjHγℓjₒjₛ_r₂r₁, 2)

	anyzeromom = l == 0 || jₛ == 0
	evenmom = iseven(jₒ + jₛ + l)

	if evenmom
		@turbo for r_ind in eachindex(r)
			TP₁₂⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₁n₂ * twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, zeroindT]
			TP₂₁⁰ = pre⁰ * Pʲₒʲₛₗₘ_₀₀_n₂n₁ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, zeroindT]

			K[rindK, r_ind, lm_ind_K] += real(TP₂₁⁰) - real(TP₁₂⁰)
		end
	end

	# l = 0 only has an r component
	anyzeromom && return K
	@turbo for r_ind in eachindex(r)
		TP₁₂¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₁n₂ * twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, plusindT]
		TP₂₁¹ = pre⁺ * Pʲₒʲₛₗₘ_₀₀_n₂n₁ * twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, plusindT]

		tempre = real(TP₂₁¹) - real(TP₁₂¹)
		tempim = imag(TP₂₁¹) - imag(TP₁₂¹)

		K[θindK, r_ind, lm_ind_K] += evenmom * 2tempre
		K[ϕindK, r_ind, lm_ind_K] += !evenmom * (-2tempim)
	end

	return K
end
function populatekernelrθϕl0!(::Flow, ::los_earth, K::AbstractArray{<:Real, 3},
	SHModes, Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
	(l,_), (jₛ, jₒ, ω), twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂,
	pre⁰, pre⁺, phase)

	lm_ind_K = modeindex(SHModes, (l, 0))
	lm_ind = modeindex(firstshmodes(Yjₒjₛ_lm_n₁n₂), (l, 0))

	_Yjₒjₛ_lm_n₁n₂ = parent(Yjₒjₛ_lm_n₁n₂[lm_ind])
	_Yjₒjₛ_lm_n₂n₁ = parent(Yjₒjₛ_lm_n₂n₁[lm_ind])

	rindK, θindK, ϕindK = axes(K, 1)
	zeroindT, plusindT = axes(twoimagconjhωconjHγℓjₒjₛ_r₂r₁, 2)

	anyzeromom = l == 0 || jₛ == 0
	evenmom = iseven(jₒ + jₛ + l)

	if evenmom
		@turbo for r_ind in eachindex(r)
			TP₁₂⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₁n₂ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, zeroindT])
			TP₂₁⁰ = pre⁰ * sum(_Yjₒjₛ_lm_n₂n₁ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, zeroindT])

			K[rindK, r_ind, lm_ind_K] += real(TP₂₁⁰) - real(TP₁₂⁰)
		end
	end
	anyzeromom && return K
	@turbo for r_ind in eachindex(r)
		TP₁₂¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₁n₂ .* twoimagconjhωHγℓjₒjₛ_r₁r₂[r_ind, plusindT])
		TP₂₁¹ = pre⁺ * sum(_Yjₒjₛ_lm_n₂n₁ .* twoimagconjhωconjHγℓjₒjₛ_r₂r₁[r_ind, plusindT])

		tempre = real(TP₂₁¹) - real(TP₁₂¹)
		tempim = imag(TP₂₁¹) - imag(TP₁₂¹)

		K[θindK, r_ind, lm_ind_K] += evenmom * 2tempre
		K[ϕindK, r_ind, lm_ind_K] += !evenmom * (-2tempim)
	end

	return K
end

reinterpretSMatrix(A::AbstractArray, ::los_earth) = dropdims(reinterpret(SMatrix{2,2,eltype(A),4}, reshape(A, 4, size(A)[3:end]...)), dims = 1)
reinterpretSMatrix(A::AbstractArray, ::los_radial) = A

function kernel_uₛₜ_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction, SHModes, hω_arr,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	@unpack p_Gsrc, p_Gobs1, p_Gobs2, Gfn_path_src, NGfn_files_src, Gfn_path_obs1, NGfn_files_obs1,
	Gfn_path_obs2, NGfn_files_obs2 = unpackGfnparams(p_Gsrc, r_src, p_Gobs1, radius(xobs1), p_Gobs2, radius(xobs2));

	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	s_range = l_range(SHModes)
	s_max = maximum(s_range)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))
	obs_at_same_height = r₁_ind == r₂_ind

	ℓ_range_proc = _extremaelementrange(ℓ_ωind_iter_on_proc, dims = 1)

	# Get a list of all modes that will be accessed.
	# This can be used to open the fits files before the loops begin.
	# This will cut down on FITS IO costs
	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs1 = Gfn_fits_files(Gfn_path_obs1,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs1)
	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs2 = Gfn_fits_files(Gfn_path_obs2,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs2)

	Kre = zeros(3, nr, length(SHModes)); # first axis is for vector indices (-1,0,1 in the helicity basis)
	K = StructArray{ComplexF64}((Kre, zero(Kre))); # Kγₗₘ(r, x₁, x₂)

	ind⁰ = 2 # index of the vector component 0

	arrs = allocatearrays(Flow(), los, obs_at_same_height);
	@unpack Gsrc, drGsrc, Gparts_r₁, Gparts_r₂,
	Gγℓjₒjₛ_r₁, Gγℓjₒjₛ_r₂, Hγℓjₒjₛ_r₁r₂, Hγℓjₒjₛ_r₂r₁ = arrs;
	@unpack Gparts_r₁_2, Gparts_r₂_2 = arrs;
	@unpack twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂ = arrs;
	twoimagconjhωconjHγℓjₒjₛ_r₂r₁P, twoimagconjhωHγℓjₒjₛ_r₁r₂P = map(parent, (twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂));

	HT = HybridArray{Tuple{sizeindG(los)...,sizeindG(los)...,StaticArrays.Dynamic(),2}}
	H₂₁ = _structarrayparent(Hγℓjₒjₛ_r₂r₁, HT)
	H₁₂ = _structarrayparent(Hγℓjₒjₛ_r₁r₂, HT)

	twoimagconjhωconjHγℓjₒjₛ_r₂r₁S = reinterpretSMatrix(twoimagconjhωconjHγℓjₒjₛ_r₂r₁P, los);
	twoimagconjhωHγℓjₒjₛ_r₁r₂S = reinterpretSMatrix(twoimagconjhωHγℓjₒjₛ_r₁r₂P, los);

	jₒjₛ_allmodes = L2L1Triangle(ℓ_range_proc, s_max, ℓ_arr)
	jₒrange = l2_range(jₒjₛ_allmodes)
	@unpack Gobs1_cache, Gobs2_cache = arrs
	for jₒ in jₒrange
		Gobs1_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
		if !obs_at_same_height
			Gobs2_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
		end
	end

	@timeit localtimer "biposh" begin
		Y12, Y21 = los_projected_spheroidal_biposh_flippoints(xobs1, xobs2, los, SHModes, jₒjₛ_allmodes)
	end

	C = zeros(0:1, 0:s_max, l2_range(jₒjₛ_allmodes), ℓ_range_proc)

	for jₛ in axes(C, 4), jₒ in axes(C, 3), ℓ in axes(C, 2)

		C[0, ℓ, jₒ, jₛ] = if iseven(jₒ+jₛ+ℓ)
							clebschgordan(Float64, jₒ, 0, jₛ, 0, ℓ, 0)
						else
							zero(Float64)
						end

		C[1, ℓ, jₒ, jₛ] = if ℓ > 0 && jₛ > 0
							clebschgordan(Float64, jₒ, 0, jₛ, 1, ℓ, 1)
						else
							zero(Float64)
						end
	end

	# for some reason the types of these variables were not being inferred
	C⁰::Float64, C⁺::Float64 = 0, 0
	phase::Int = 1
	pre⁰::Float64, pre⁺::Float64 = 0, 0

	ω_ind_prev = -1 # something unrealistic to start off with

	# Loop over the Greenfn files
	for (jₛ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		h_ω = hω_arr[ω_ind]
		dωω³Pω = dω/2π * ω^3 * Powspec(ω)
		Ωjₛ0 = Ω(jₛ, 0)

		@timeit localtimer "FITS" begin
		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about the source radius
		read_Gfn_file_at_index!(drGsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 2)
		end # timer

		Gα₂r_r₁_rₛ = αrcomp(Gsrc, r₁_ind, los)
		Gα₂r_r₂_rₛ = αrcomp(Gsrc, r₂_ind, los)

		for jₒ in l2_range(jₒjₛ_allmodes, jₛ)

			jₒjₛ_ind = modeindex(jₒjₛ_allmodes, (jₒ, jₛ))

			Yjₒjₛ_lm_n₁n₂ = Y12[jₒjₛ_ind]
			Yjₒjₛ_lm_n₂n₁ = Y21[jₒjₛ_ind]

			(Yjₒjₛ_lm_n₁n₂ === nothing || Yjₒjₛ_lm_n₂n₁ === nothing) && continue

			# Change shmodes to ML to avoid recomputing the radial term for same l
			modesML = intersect(ML(firstshmodes(Yjₒjₛ_lm_n₁n₂)), ML(SHModes))
			modesML === nothing && continue

			Gobs1 = Gobs1_cache[jₒ]
			@timeit localtimer "FITS" begin
			# Green functions based at the observation point
			if ω_ind != ω_ind_prev || jₒ ∉ l2_range(jₒjₛ_allmodes, jₛ-1)
				read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
				(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)
			end
			end # timer

			@timeit localtimer "radial term 1" begin
			Gⱼₒⱼₛω_u⁰⁺_firstborn!(Gparts_r₁, Gparts_r₁_2, Gsrc, drGsrc, jₛ, Gobs1, jₒ, los)
			end # timer

			if !obs_at_same_height
				Gobs2 = Gobs2_cache[jₒ]
				@timeit localtimer "FITS" begin
				if ω_ind != ω_ind_prev || jₒ ∉ l2_range(jₒjₛ_allmodes, jₛ-1)
					read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
					(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)
				end
				end # timer

				@timeit localtimer "radial term 1" begin
				Gⱼₒⱼₛω_u⁰⁺_firstborn!(Gparts_r₂, Gparts_r₂_2, Gsrc, drGsrc, jₛ, Gobs2, jₒ, los)
				end #timer
			end

			l_prev = -1 # something unrealistic to start off

			for (l, m) in modesML
				# Check if triangle condition is satisfied
				# The outer loop is over all possible jₛ and jₒ for a given l_max
				# Not all combinations would contribute towards a specific l
				δ(jₛ, jₒ, l) || continue

				if l != l_prev
					@timeit localtimer "radial term 2" begin
						@timeit localtimer "H" begin
							Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(Hγℓjₒjₛ_r₁r₂, Gparts_r₁_2, Gγℓjₒjₛ_r₁, Gα₂r_r₂_rₛ, jₒ, jₛ, l)

							if !obs_at_same_height
								Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(Hγℓjₒjₛ_r₂r₁, Gparts_r₂_2, Gγℓjₒjₛ_r₂, Gα₂r_r₁_rₛ, jₒ, jₛ, l)
							end
						end

						@timeit localtimer "T" begin
							@turbo for I in CartesianIndices(H₂₁)
								twoimagconjhωconjHγℓjₒjₛ_r₂r₁P[I] = -2imag(h_ω * H₂₁[I])
								twoimagconjhωHγℓjₒjₛ_r₁r₂P[I] = 2imag(conj(h_ω) * H₁₂[I])
							end
						end

						C⁰ = C[0, l, jₒ, jₛ]
						C⁺ = C[1, l, jₒ, jₛ]
						phase = (-1)^(jₒ + jₛ + l)
						coeffj = √((2jₒ+1)*(2jₛ+1)/π/(2l+1))

						pre⁰ = dωω³Pω * coeffj * C⁰
						pre⁺ = dωω³Pω * coeffj * Ωjₛ0 * C⁺
						l_prev = l
					end
				end

				@timeit localtimer "kernel" begin
				populatekernel!(Flow(), los, K,
					SHModes,
					Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
					(l, m), (jₛ, jₒ, ω),
					twoimagconjhωHγℓjₒjₛ_r₁r₂S,
					twoimagconjhωconjHγℓjₒjₛ_r₂r₁S,
					pre⁰, pre⁺, phase)
				end #timer
			end
		end

		ω_ind_prev = ω_ind
	end

	foreach(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))

	_K = permutedims(K, [2,1,3])
	mulprefactor!(_K, ind⁰)
	AK = Array(no_offset_view(_K))

	return AK
end

function kernel_uₛₜ_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	(xobs1, xobs1′)::NTuple{2,SphericalPoint},
	(xobs2, xobs2′)::NTuple{2,SphericalPoint}, los::los_direction, SHModes, hω_arr,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	@unpack p_Gsrc, p_Gobs1, p_Gobs2, Gfn_path_src, NGfn_files_src, Gfn_path_obs1, NGfn_files_obs1,
	Gfn_path_obs2, NGfn_files_obs2 = unpackGfnparams(p_Gsrc, r_src, p_Gobs1, radius(xobs1), p_Gobs2, radius(xobs2));

	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	s_range = l_range(SHModes)
	s_max = maximum(s_range)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))
	obs_at_same_height = r₁_ind == r₂_ind

	ℓ_range_proc = _extremaelementrange(ℓ_ωind_iter_on_proc, dims = 1)

	# Get a list of all modes that will be accessed.
	# This can be used to open the fits files before the loops begin.
	# This will cut down on FITS IO costs
	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs1 = Gfn_fits_files(Gfn_path_obs1,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs1)
	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs2 = Gfn_fits_files(Gfn_path_obs2,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs2)

	Kre = zeros(2, 3, nr, length(SHModes)); # first axis is for vector indices (-1,0,1 in the helicity basis)
	K = StructArray{ComplexF64}((Kre, zero(Kre))); # Kγₗₘ(r, x₁, x₂)

	ind⁰ = 2 # index of the vector component 0

	arrs = allocatearrays(Flow(), los, obs_at_same_height, Val(2));
	@unpack Gsrc, drGsrc, Gparts_r₁, Gparts_r₂,
	Gγℓjₒjₛ_r₁, Gγℓjₒjₛ_r₂, Hγℓjₒjₛ_r₁r₂, Hγℓjₒjₛ_r₂r₁ = arrs;
	@unpack Gparts_r₁_2, Gparts_r₂_2 = arrs;
	@unpack twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂ = arrs;
	twoimagconjhωconjHγℓjₒjₛ_r₂r₁P, twoimagconjhωHγℓjₒjₛ_r₁r₂P = map(parent, (twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂));

	HT = HybridArray{Tuple{sizeindG(los)...,sizeindG(los)...,StaticArrays.Dynamic(),2}}
	H₂₁ = _structarrayparent(Hγℓjₒjₛ_r₂r₁, HT)
	H₁₂ = _structarrayparent(Hγℓjₒjₛ_r₁r₂, HT)

	twoimagconjhωconjHγℓjₒjₛ_r₂r₁S = reinterpretSMatrix(twoimagconjhωconjHγℓjₒjₛ_r₂r₁P, los);
	twoimagconjhωHγℓjₒjₛ_r₁r₂S = reinterpretSMatrix(twoimagconjhωHγℓjₒjₛ_r₁r₂P, los);

	inds_lead_T = CartesianIndices((sizeindG(los)..., sizeindG(los)...))
	inds_trail_T = CartesianIndices(axes(H₂₁)[ndims(inds_lead_T) .+ (1:2)])

	jₒjₛ_allmodes = L2L1Triangle(ℓ_range_proc, s_max, ℓ_arr)
	jₒrange = l2_range(jₒjₛ_allmodes)
	@unpack Gobs1_cache, Gobs2_cache = arrs
	for jₒ in jₒrange
		Gobs1_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
		if !obs_at_same_height
			Gobs2_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
		end
	end

	@timeit localtimer "biposh" begin
		Y12, Y21, Y1′2′, Y2′1′ =
			los_projected_spheroidal_biposh_flippoints((xobs1, xobs1′),
			(xobs2, xobs2′), los, SHModes, jₒjₛ_allmodes)
	end

	C = zeros(0:1, 0:s_max, l2_range(jₒjₛ_allmodes), ℓ_range_proc)

	for jₛ in axes(C, 4), jₒ in axes(C, 3), ℓ in axes(C, 2)

		C[0, ℓ, jₒ, jₛ] = if iseven(jₒ+jₛ+ℓ)
							clebschgordan(Float64, jₒ, 0, jₛ, 0, ℓ, 0)
						else
							zero(Float64)
						end

		C[1, ℓ, jₒ, jₛ] = if ℓ > 0 && jₛ > 0
							clebschgordan(Float64, jₒ, 0, jₛ, 1, ℓ, 1)
						else
							zero(Float64)
						end
	end

	# for some reason the types of these variables were not being inferred
	C⁰::Float64, C⁺::Float64 = 0, 0
	phase::Int = 1
	pre⁰::Float64, pre⁺::Float64 = 0, 0

	ω_ind_prev = -1 # something unrealistic to start off with

	# Loop over the Greenfn files
	for (jₛ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		h_ω_1 = hω_arr[ω_ind, 1]
		h_ω_2 = hω_arr[ω_ind, 2]
		h_ω = SVector{2}(h_ω_1, h_ω_2)
		dωω³Pω = dω/2π * ω^3 * Powspec(ω)
		Ωjₛ0 = Ω(jₛ, 0)

		@timeit localtimer "FITS" begin
		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about the source radius
		read_Gfn_file_at_index!(drGsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 2)
		end # timer

		Gα₂r_r₁_rₛ = αrcomp(Gsrc, r₁_ind, los)
		Gα₂r_r₂_rₛ = αrcomp(Gsrc, r₂_ind, los)

		for jₒ in l2_range(jₒjₛ_allmodes, jₛ)

			jₒjₛ_ind = modeindex(jₒjₛ_allmodes, (jₒ, jₛ))

			Yjₒjₛ_lm_n₁n₂ = Y12[jₒjₛ_ind]
			Yjₒjₛ_lm_n₂n₁ = Y21[jₒjₛ_ind]
			Yjₒjₛ_lm_n₁′n₂′ = Y1′2′[jₒjₛ_ind]
			Yjₒjₛ_lm_n₂′n₁′ = Y2′1′[jₒjₛ_ind]

			(Yjₒjₛ_lm_n₁n₂ === nothing || Yjₒjₛ_lm_n₂n₁ === nothing) && continue
			(Yjₒjₛ_lm_n₁′n₂′ === nothing || Yjₒjₛ_lm_n₂′n₁′ === nothing) && continue

			# Change shmodes to ML to avoid recomputing the radial term for same l
			modesML = intersect(ML(firstshmodes(Yjₒjₛ_lm_n₁n₂)), ML(SHModes))
			modesML === nothing && continue

			Gobs1 = Gobs1_cache[jₒ]
			@timeit localtimer "FITS" begin
			# Green functions based at the observation point
			if ω_ind != ω_ind_prev || jₒ ∉ l2_range(jₒjₛ_allmodes, jₛ-1)
				read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
				(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)
			end
			end # timer

			@timeit localtimer "radial term 1" begin
			Gⱼₒⱼₛω_u⁰⁺_firstborn!(Gparts_r₁, Gparts_r₁_2, Gsrc, drGsrc, jₛ, Gobs1, jₒ, los)
			end # timer

			if !obs_at_same_height
				Gobs2 = Gobs2_cache[jₒ]
				@timeit localtimer "FITS" begin
				if ω_ind != ω_ind_prev || jₒ ∉ l2_range(jₒjₛ_allmodes, jₛ-1)
					read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
					(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)
				end
				end # timer

				@timeit localtimer "radial term 1" begin
				Gⱼₒⱼₛω_u⁰⁺_firstborn!(Gparts_r₂, Gparts_r₂_2, Gsrc, drGsrc, jₛ, Gobs2, jₒ, los)
				end #timer
			end

			l_prev = -1 # something unrealistic to start off

			for (l, m) in modesML
				# Check if triangle condition is satisfied
				# The outer loop is over all possible jₛ and jₒ for a given l_max
				# Not all combinations would contribute towards a specific l
				δ(jₛ, jₒ, l) || continue

				if l != l_prev
					@timeit localtimer "radial term 2" begin
						@timeit localtimer "H" begin
							Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(Hγℓjₒjₛ_r₁r₂, Gparts_r₁_2, Gγℓjₒjₛ_r₁, Gα₂r_r₂_rₛ, jₒ, jₛ, l)

							if !obs_at_same_height
								Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(Hγℓjₒjₛ_r₂r₁, Gparts_r₂_2, Gγℓjₒjₛ_r₂, Gα₂r_r₁_rₛ, jₒ, jₛ, l)
							end
						end

						@timeit localtimer "T" begin
							for Itrail in inds_trail_T, nind in eachindex(h_ω), Ilead in inds_lead_T
								twoimagconjhωconjHγℓjₒjₛ_r₂r₁P[Ilead, nind, Itrail] = -2imag(h_ω[nind] * H₂₁[Ilead, Itrail])
								twoimagconjhωHγℓjₒjₛ_r₁r₂P[Ilead, nind, Itrail] = 2imag(conj(h_ω[nind]) * H₁₂[Ilead, Itrail])
							end
						end

						C⁰ = C[0, l, jₒ, jₛ]
						C⁺ = C[1, l, jₒ, jₛ]
						phase = (-1)^(jₒ + jₛ + l)
						coeffj = √((2jₒ+1)*(2jₛ+1)/π/(2l+1))

						pre⁰ = dωω³Pω * coeffj * C⁰
						pre⁺ = dωω³Pω * coeffj * Ωjₛ0 * C⁺
						l_prev = l
					end
				end

				@timeit localtimer "kernel" begin
				populatekernel!(Flow(), los, K,
					SHModes,
					(Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₁′n₂′),
					(Yjₒjₛ_lm_n₂n₁, Yjₒjₛ_lm_n₂′n₁′),
					(l, m), (jₛ, jₒ, ω),
					twoimagconjhωHγℓjₒjₛ_r₁r₂S,
					twoimagconjhωconjHγℓjₒjₛ_r₂r₁S,
					pre⁰, pre⁺, phase)
				end #timer
			end
		end

		ω_ind_prev = ω_ind
	end

	foreach(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))

	_K = permutedims(K, [3,2,4,1])
	mulprefactor!(_K, ind⁰)
	AK = Array(no_offset_view(_K))

	return AK
end

# Compute Kₛₜ first and then compute Kₛ₀_rθϕ from that
function kernel_uₛ₀_rθϕ_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction, SHModes, hω_arr,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	K = kernel_uₛₜ_partial(localtimer, ℓ_ωind_iter_on_proc,
		xobs1, xobs2, los, SHModes, hω_arr,
		p_Gsrc, p_Gobs1, p_Gobs2, r_src)

	K_rθϕ = zeros(size(K, 1), 3, size(K, 3))

	ind⁰, ind⁺ = UnitRange(axes(K,2))[2:3]

	@turbo for st_ind in UnitRange(axes(K, 3))
		for r_ind in UnitRange(axes(K, 1))
			# r-component
			K_rθϕ[r_ind, 1, st_ind] = real(K[r_ind, ind⁰, st_ind])
		end
		for r_ind in UnitRange(axes(K, 1))
			# θ-component
			K_rθϕ[r_ind, 2, st_ind] = 2real(K[r_ind, ind⁺, st_ind])
		end
		for r_ind in UnitRange(axes(K, 1))
			# ϕ-component
			K_rθϕ[r_ind, 3, st_ind] = -2imag(K[r_ind, ind⁺, st_ind])
		end
	end

	K_rθϕ
end

# Compute Kₛ₀_rθϕ directly
function kernel_uₛ₀_rθϕ_partial_2(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction, SHModes, hω_arr,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	@unpack p_Gsrc, p_Gobs1, p_Gobs2, Gfn_path_src, NGfn_files_src, Gfn_path_obs1, NGfn_files_obs1,
	Gfn_path_obs2, NGfn_files_obs2 = unpackGfnparams(p_Gsrc, r_src, p_Gobs1, radius(xobs1), p_Gobs2, radius(xobs2))

	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	s_max = l_max(SHModes)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))
	obs_at_same_height = r₁_ind == r₂_ind

	ℓ_range_proc = _extremaelementrange(ℓ_ωind_iter_on_proc, dims = 1)

	# Get a list of all modes that will be accessed.
	# This can be used to open the fits files before the loops begin.
	# This will cut down on FITS IO costs
	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs1 = Gfn_fits_files(Gfn_path_obs1,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs1)
	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs2 = Gfn_fits_files(Gfn_path_obs2,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs2)

	K = zeros(3, nr, length(SHModes)) # Kγₗₘ(r, x₁, x₂)
	indr = 1

	arrs = allocatearrays(Flow(), los, obs_at_same_height)
	@unpack Gsrc, drGsrc, Gobs1, Gobs2, Gparts_r₁, Gparts_r₂, Gparts_r₁_2, Gparts_r₂_2,
	Gγℓjₒjₛ_r₁, Gγℓjₒjₛ_r₂, Hγℓjₒjₛ_r₁r₂, Hγℓjₒjₛ_r₂r₁ = arrs
	@unpack twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂ = arrs
	twoimagconjhωconjHγℓjₒjₛ_r₂r₁P, twoimagconjhωHγℓjₒjₛ_r₁r₂P = map(parent, (twoimagconjhωconjHγℓjₒjₛ_r₂r₁, twoimagconjhωHγℓjₒjₛ_r₁r₂))

	HT = HybridArray{Tuple{sizeindG(los)...,sizeindG(los)...,StaticArrays.Dynamic(),2}}
	H₂₁R = HT(parent(Hγℓjₒjₛ_r₂r₁.re))
	H₂₁I = HT(parent(Hγℓjₒjₛ_r₂r₁.im))
	H₁₂R = HT(parent(Hγℓjₒjₛ_r₁r₂.re))
	H₁₂I = HT(parent(Hγℓjₒjₛ_r₁r₂.im))

	twoimagconjhωconjHγℓjₒjₛ_r₂r₁S = reinterpretSMatrix(twoimagconjhωconjHγℓjₒjₛ_r₂r₁P, los)
	twoimagconjhωHγℓjₒjₛ_r₁r₂S = reinterpretSMatrix(twoimagconjhωHγℓjₒjₛ_r₁r₂P, los)

	jₒjₛ_allmodes = L2L1Triangle(ℓ_range_proc, s_max, ℓ_arr)

	@timeit localtimer "biposh" begin
		Y12, Y21 = los_projected_spheroidal_biposh_flippoints(xobs1, xobs2, los, SHModes, jₒjₛ_allmodes)
	end

	C = zeros(0:1, 0:s_max, l2_range(jₒjₛ_allmodes), ℓ_range_proc)

	@timeit localtimer "CG" begin
	for jₛ in axes(C, 4), jₒ in axes(C, 3), ℓ in axes(C, 2)

		C[0, ℓ, jₒ, jₛ] = if iseven(jₒ+jₛ+ℓ)
							clebschgordan(Float64, jₒ, 0, jₛ, 0, ℓ, 0)
						else
							zero(Float64)
						end

		C[1, ℓ, jₒ, jₛ] = if ℓ > 0 && jₛ > 0
							clebschgordan(Float64, jₒ, 0, jₛ, 1, ℓ, 1)
						else
							zero(Float64)
						end
	end
	end # timer

	C⁰, C⁺ = zero(Float64), zero(Float64)
	phase = 1
	pre⁰, pre⁺ = zero(Float64), zero(Float64)

	# Loop over the Greenfn files
	for (jₛ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		h_ω = hω_arr[ω_ind]
		dωω³Pω = dω/2π * ω^3 * Powspec(ω)
		Ωjₛ0 = Ω(jₛ, 0)

		@timeit localtimer "FITS" begin
		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about the source radius
		read_Gfn_file_at_index!(drGsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 2)
		end # timer

		Gα₂r_r₁_rₛ = αrcomp(Gsrc, r₁_ind, los)
		Gα₂r_r₂_rₛ = αrcomp(Gsrc, r₂_ind, los)

		for jₒ in l2_range(jₒjₛ_allmodes, jₛ)

			Yjₒjₛ_lm_n₁n₂ = Y12[(jₒ, jₛ)]
			Yjₒjₛ_lm_n₂n₁ = Y21[(jₒ, jₛ)]

			(Yjₒjₛ_lm_n₁n₂ === nothing || Yjₒjₛ_lm_n₂n₁ === nothing) && continue

			@timeit localtimer "FITS" begin
			# Green functions based at the observation point
			read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
			(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)
			end # timer

			@timeit localtimer "radial term 1" begin
			Gⱼₒⱼₛω_u⁰⁺_firstborn!(Gparts_r₁, Gparts_r₁_2, Gsrc, drGsrc, jₛ, Gobs1, jₒ, los)
			end # timer

			if !obs_at_same_height
				@timeit localtimer "FITS" begin
					read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
					(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)
				end # timer

				@timeit localtimer "radial term 1" begin
				Gⱼₒⱼₛω_u⁰⁺_firstborn!(Gparts_r₂, Gparts_r₂_2, Gsrc, drGsrc, jₛ, Gobs2, jₒ, los)
				end #timer
			end

			for (l, m) in firstshmodes(Yjₒjₛ_lm_n₁n₂)
				# Check if triangle condition is satisfied
				# The outer loop is over all possible jₛ and jₒ for a given s_max
				# Not all combinations would contribute towards a specific l
				δ(jₛ, jₒ, l) || continue

				@timeit localtimer "radial term 2" begin
					@timeit localtimer "H" begin
						Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(Hγℓjₒjₛ_r₁r₂, Gparts_r₁_2, Gγℓjₒjₛ_r₁, Gα₂r_r₂_rₛ, jₒ, jₛ, l)

						if !obs_at_same_height
							Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(Hγℓjₒjₛ_r₂r₁, Gparts_r₂_2, Gγℓjₒjₛ_r₂, Gα₂r_r₁_rₛ, jₒ, jₛ, l)
						end
					end
					# @. twoimagconjhωconjHγℓjₒjₛ_r₂r₁ = -2imag(hω * Hγℓjₒjₛ_r₂r₁)
					# @. twoimagconjhωHγℓjₒjₛ_r₁r₂ = 2imag(conjhω * Hγℓjₒjₛ_r₁r₂)
					@timeit localtimer "T" begin
						@turbo for I in CartesianIndices(H₂₁R)
							twoimagconjhωconjHγℓjₒjₛ_r₂r₁P[I] = -2(real(h_ω) * H₂₁I[I] + imag(h_ω) * H₂₁R[I])
							twoimagconjhωHγℓjₒjₛ_r₁r₂P[I] = 2(real(h_ω) * H₁₂I[I] - imag(h_ω) * H₁₂R[I])
						end
					end
				end #timer

				C⁰ = C[0, l, jₒ, jₛ]
				C⁺ = C[1, l, jₒ, jₛ]
				phase = (-1)^(jₒ + jₛ + l)
				coeffj = √((2jₒ+1)*(2jₛ+1)/π/(2l+1))

				pre⁰ = dωω³Pω * coeffj * C⁰
				pre⁺ = dωω³Pω * coeffj * Ωjₛ0 * C⁺

				@timeit localtimer "kernel" begin
				populatekernelrθϕl0!(Flow(), los, K, SHModes, Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
				(l, m), (jₛ, jₒ, ω), twoimagconjhωconjHγℓjₒjₛ_r₂r₁S, twoimagconjhωHγℓjₒjₛ_r₁r₂S,
				pre⁰, pre⁺, phase)
				end #timer
			end
		end
	end

	foreach(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))
	_K = Base.PermutedDimsArray(K, (2,1,3))
	mulprefactor!(_K, indr)

	return Array(_K)
end

function generatefitsheader(xobs1::SphericalPoint, xobs2::SphericalPoint, SHModes, j_range, ν_ind_range)
	FITSHeader(["r1","th1","phi1","r2","th2","phi2",
		"l_min","m_min","l_max","m_max",
		"j_min","j_max","nui_min","nui_max"],
		Any[float(radius(xobs1)), float(xobs1.θ), float(xobs1.ϕ),
		float(radius(xobs2)), float(xobs2.θ), float(xobs2.ϕ),
		l_min(SHModes), m_min(SHModes), l_max(SHModes), m_max(SHModes),
		minimum(j_range), maximum(j_range),
		minimum(ν_ind_range), maximum(ν_ind_range)],
		["Radius of the first observation point",
		"Colatitude of the first observation point",
		"Azimuth of the first observation point",
		"Radius of the second observation point",
		"Colatitude of the second observation point",
		"Azimuth of the second observation point",
		"Minimum angular degree of the flow",
		"Minimum azimuthal order of the flow",
		"Maximum angular degree of the flow",
		"Maximum azimuthal order of the flow",
		"Minimum wave mode","Maximum wave mode",
		"Minimum wave frequency index","Maximum wave frequency index"])
end
function generatefitsheader((xobs1, xobs1′)::NTuple{2,SphericalPoint},
		(xobs2, xobs2′)::NTuple{2,SphericalPoint}, SHModes, j_range, ν_ind_range)

	f = generatefitsheader(xobs1, xobs2, SHModes, j_range, ν_ind_range)

	for (k, v, c) in zip(["r1p","th1p","phi1p","r2p","th2p","phi2p"],
		Any[radius(xobs1′), Point2D(xobs1′)..., radius(xobs2′), Point2D(xobs2′)...],
		["Radius of the first observation point rotated",
		"Colatitude of the first observation point rotated",
		"Azimuth of the first observation point rotated",
		"Radius of the second observation point rotated",
		"Colatitude of the second observation point rotated",
		"Azimuth of the second observation point rotated"])

		f[k] = v
		set_comment!(f, k, c)
	end

	return f
end

function modetag(j_range, SHModes)
   "jmax$(maximum(j_range))_lmax$(l_max(SHModes))_mmax$(m_max(SHModes))"
end
function kernelfilenameuₛₜ(m::SeismicMeasurement, los::los_direction, j_range, SHModes, tag="")
	mtag = m isa TravelTimes ? "δτ" : "A"
	lostag = los isa los_radial ? "" : "_los"
	"Kst_$(mtag)_u_$(modetag(j_range, SHModes))$(lostag)"*string(tag)*".fits"
end

function kernel_uₛₜ(comm, m::SeismicMeasurement, xobs1, xobs2,
	los::los_direction = los_radial(); kwargs...)

	SHModes = getkernelmodes(; kwargs...)

	r_src, = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs...)
	@unpack ν_arr, ℓ_arr = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)
	j_range, ν_ind_range = iters

	hω_arr = get(kwargs, :hω) do
		h = hω(comm, m, xobs1, xobs2, los; kwargs..., print_timings = false)
		_broadcast(h, 0, comm)
	end

	K_δτ_uₛₜ = pmapsum(comm, kernel_uₛₜ_partial, iters,
		xobs1, xobs2, los, SHModes, hω_arr, p_Gsrc, p_Gobs1, p_Gobs2, r_src)

	K_δτ_uₛₜ === nothing  && return nothing
	K = OffsetArray(K_δτ_uₛₜ, :, -1:1, ..)

	if get(kwargs, :save, true)

		filepath = joinpath(SCRATCH_kerneldir[], kernelfilenameuₛₜ(m, los, j_range, SHModes))

		header = generatefitsheader(xobs1, xobs2, SHModes, j_range, ν_ind_range)

		FITS(filepath,"w") do f
			FITSIO.write(f, reinterpret_as_float(K_δτ_uₛₜ), header = header)
		end
	end

	SHArray(K, (axes(K)[1:2]..., SHModes, axes(K)[4:end]...))
end

kernelfilenameuₛ₀rθϕ(::TravelTimes, ::los_radial) = "Kl0_δτ_u_rθϕ.fits"
kernelfilenameuₛ₀rθϕ(::TravelTimes, ::los_earth) = "Kl0_δτ_u_rθϕ_los.fits"
kernelfilenameuₛ₀rθϕ(::Amplitudes, ::los_radial) = "Kl0_A_u_rθϕ.fits"
kernelfilenameuₛ₀rθϕ(::Amplitudes, ::los_earth) = "Kl0_A_u_rθϕ_los.fits"

# Compute Kₛₜ first and then compute Kₛ₀_rθϕ from that
function kernel_uₛ₀_rθϕ(comm, m::SeismicMeasurement, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	_SHModes = getkernelmodes(; kwargs...)
	SHModes = LM(l_range(_SHModes), 0)

	r_src, = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs...)
	@unpack ν_arr, ℓ_arr = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)
	j_range, ν_ind_range = iters

	hω_arr = get(kwargs, :hω) do
		h = hω(comm, m, xobs1, xobs2, los; kwargs..., print_timings = false)
		_broadcast(h, 0, comm)
	end

	K_δτ_uₛ₀ = pmapsum(comm, kernel_uₛ₀_rθϕ_partial, iters,
		xobs1, xobs2, los, SHModes, hω_arr, p_Gsrc, p_Gobs1, p_Gobs2, r_src)

	K_δτ_uₛ₀ === nothing  && return nothing

	if get(kwargs, :save, true)
		filename = joinpath(SCRATCH_kerneldir[], kernelfilenameuₛ₀rθϕ(m, los))

		header = generatefitsheader(xobs1, xobs2, SHModes, j_range, ν_ind_range)

		FITS(filename,"w") do f
			FITSIO.write(f, reinterpret_as_float(K_δτ_uₛ₀), header = header)
		end
	end

	SHArray(K_δτ_uₛ₀, (axes(K_δτ_uₛ₀)[1:2]..., SHModes))
end

# Compute Kₛ₀_rθϕ directly
function kernel_uₛ₀_rθϕ_2(comm, m::SeismicMeasurement, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	_SHModes = getkernelmodes(; kwargs...)
	SHModes = LM(l_range(_SHModes), 0)

	r_src, = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs...)
	@unpack ν_arr, ℓ_arr = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)
	j_range, ν_ind_range = iters

	hω_arr = get(kwargs, :hω) do
		h = hω(comm, m, xobs1, xobs2, los; kwargs..., print_timings = false)
		_broadcast(h, 0, comm)
	end

	K_δτ_uₛ₀ = pmapsum(comm, kernel_uₛ₀_rθϕ_partial_2, iters,
		xobs1, xobs2, los, SHModes, hω_arr, p_Gsrc, p_Gobs1, p_Gobs2, r_src)

	K_δτ_uₛ₀ === nothing  && return nothing
	if get(kwargs, :save, true)
		filename = joinpath(SCRATCH_kerneldir[], kernelfilenameuₛ₀rθϕ(m, los))

		header = generatefitsheader(xobs1, xobs2, SHModes, j_range, ν_ind_range)

		FITS(filename,"w") do f
			FITSIO.write(f, reinterpret_as_float(K_δτ_uₛ₀), header = header)
		end
	end

	SHArray(K_δτ_uₛ₀, (axes(K_δτ_uₛ₀)[1:2]..., SHModes))
end

function kernel_ψϕₛ₀(comm, m::SeismicMeasurement, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	Kv = kernel_uₛ₀_rθϕ(comm, m, xobs1, xobs2, los; kwargs..., save = false)
	Kv === nothing && return nothing
	kernel_ψϕₛ₀(Kv; kwargs...)
end

function kernel_ψϕₛ₀(Kv::AbstractArray; kwargs...)
	mode_range = firstshmodes(Kv);
	Kψϕ = zeros(axes(Kv, 1), mode_range);
	kernel_ψϕₛ₀!(Kψϕ, Kv; save = get(kwargs, :save, true))

	if get(kwargs, :save, true)
		filename = joinpath(SCRATCH_kerneldir[], "Kpsi_imag.fits")
		FITS(filename,"w") do f
			write(f, reinterpret_as_float(Kψϕ))
		end
	end

	return Kψϕ
end

function kernel_ψϕₛ₀!(Kψϕ::SHArray, Kv::SHArray; kwargs...)
	mode_range = firstshmodes(Kv)

	for l in l_range(mode_range)
		@views Kψϕ[:, (l, 0)] .= -(ddr*(Kv[:, 2, (l, 0)]./ρ) .+
			@. (Kv[:, 2, (l, 0)] - 2Ω(l, 0)*Kv[:, 1, (l, 0)])/(ρ*r) )
	end

	return Kψϕ
end

#################################################################################################################
# Validation for isotropic sound speed perturbation
#################################################################################################################

function populatekernelvalidation!(::SoundSpeed, ::los_radial, K::AbstractVector,
	(ℓ, ω), Y12ℓ, conjhω, H¹₁jjω_r₁r₂, H¹₁jj_r₂r₁, dω)

	pre = 1/√(4π) * dω/2π * ω^2 * Powspec(ω) * Y12ℓ
	@. K += pre * 2real(conjhω * ( H¹₁jjω_r₁r₂ + conj(H¹₁jj_r₂r₁) ) )
end

function populatekernelvalidation!(::SoundSpeed, ::los_radial, K::AbstractMatrix,
	(ℓ, ω), Y12ℓ, conjhω, H¹₁jjω_r₁r₂::AbstractVector{<:Complex},
	H¹₁jj_r₂r₁::AbstractVector{<:Complex}, dω)

	pre = 1/√(4π) * dω/2π * ω^2 * Powspec(ω)

	for n2ind in UnitRange(axes(K, 2))
		temp_ωℓ_n2 = pre * Y12ℓ[n2ind]
		conjhω_n2 = conjhω[n2ind]
		for r_ind in UnitRange(axes(K, 1))
			K[r_ind, n2ind] += temp_ωℓ_n2 *
				2real(conjhω_n2 * ( H¹₁jjω_r₁r₂[r_ind] + conj(H¹₁jj_r₂r₁[r_ind]) ) )
		end
	end
end

function populatekernelvalidation!(::SoundSpeed, ::los_earth, K::AbstractVector,
	(j, ω), Y12j, conjhω, Hjj_r₁r₂, Hjj_r₂r₁, dω)

	pre = dω/2π * ω^2 * Powspec(ω) * √(1/π)

	for r_ind in UnitRange(axes(K, 1))
		for α₂ in 0:1, α₁ in 0:1

			llY = pre * real(Y12j[α₁, α₂])
			iszero(llY) && continue

			K[r_ind] += real(conjhω *
				( Hjj_r₁r₂[α₁, α₂, r_ind] + conj(Hjj_r₂r₁[α₂, α₁, r_ind]) ) ) * llY
		end
	end
end

function populatekernelvalidation!(::SoundSpeed, ::los_earth, K::AbstractMatrix,
	(j, ω), Y12j, conjhω, Hjjω_r₁r₂, Hjj_r₂r₁, dω)

	pre = dω/2π * ω^2 * Powspec(ω) * √(1/π)

	for n2ind in axes(K, 2)
		Yjn₂ = Y12j[n2ind]
		conjhω_n2 = conjhω[n2ind]

		for r_ind in UnitRange(axes(K, 1))
			for α₂ in 0:1, α₁ in 0:1

				llY = pre * real(Yjn₂[α₁, α₂])
				iszero(llY) && continue

				K[r_ind, n2ind] += real(conjhω_n2 *
					( Hjjω_r₁r₂[α₁, α₂, r_ind] + conj(Hjj_r₂r₁[α₂, α₁, r_ind]) ) ) * llY
			end
		end
	end
end

function kernel_δc₀₀_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction, hω,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	@unpack p_Gsrc, p_Gobs1, p_Gobs2, Gfn_path_src, NGfn_files_src,
	Gfn_path_obs1, NGfn_files_obs1, Gfn_path_obs2, NGfn_files_obs2 =
		unpackGfnparams(p_Gsrc, r_src, p_Gobs1, radius(xobs1), p_Gobs2, radius(xobs2))

	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2 =
		Gfn_fits_files.((Gfn_path_src, Gfn_path_obs1, Gfn_path_obs2),
			((ℓ_arr, 1:Nν_Gfn),), (ℓ_ωind_iter_on_proc,),
			(NGfn_files_src, NGfn_files_obs1, NGfn_files_obs2))

	K = zeros(nr)

	arrs = allocatearrays(SoundSpeed(), los, r₁_ind == r₂_ind);
	@unpack Gsrc, drGsrc, divGsrc, Gobs1, drGobs1, Gobs2, divGobs = arrs;
	fjj_r₁_rsrc, fjj_r₂_rsrc = arrs.fjₒjₛ_r₁_rsrc, arrs.fjₒjₛ_r₂_rsrc;
	H¹₁jj_r₁r₂, H¹₁jj_r₂r₁ = arrs.Hjₒjₛω_r₁r₂, arrs.Hjₒjₛω_r₂r₁;

	Y12 = los_projected_biposh_spheroidal(computeY₀₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		conjhω = conj(hω[ω_ind])

		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about the source radius
		read_Gfn_file_at_index!(drGsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:1, srcindFITS(los), :, 2)

		Gγr_r₁_rsrc = αrcomp(Gsrc, r₁_ind, los)
		Gγr_r₂_rsrc = αrcomp(Gsrc, r₂_ind, los)

		# Green function about receiver location
		read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about receiver location
		read_Gfn_file_at_index!(drGobs1, Gfn_fits_files_obs1,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs1, 1:1, srcindFITS(los), :, 2)

		radial_fn_isotropic_δc_firstborn!(fjj_r₁_rsrc,
			Gsrc, drGsrc, divGsrc, Gobs1, drGobs1, divGobs, ℓ)

		Hjₒjₛω!(H¹₁jj_r₁r₂, fjj_r₁_rsrc, Gγr_r₂_rsrc)

		if r₁_ind != r₂_ind
			# Green function about receiver location
			read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
				(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)

			# Derivative of Green function about receiver location
			read_Gfn_file_at_index!(drGobs2, Gfn_fits_files_obs2,
				(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs2, 1:1, srcindFITS(los), :, 2)

			radial_fn_isotropic_δc_firstborn!(fjj_r₂_rsrc,
				Gsrc, drGsrc, divGsrc, Gobs2, drGobs2, divGobs, ℓ)

			Hjₒjₛω!(H¹₁jj_r₂r₁, fjj_r₂_rsrc, Gγr_r₁_rsrc)
		end

		populatekernelvalidation!(SoundSpeed(), los, K,
		(ℓ, ω), Y12[ℓ], conjhω, H¹₁jj_r₁r₂, H¹₁jj_r₂r₁, dω)
	end

	map(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))
	return K
end

function kernel_δc₀₀_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	nobs1::Point2D, nobs2_arr::Vector{<:Point2D}, los::los_direction, hω_arr,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = p_Gobs1,
	r_src = r_src_default)

	hω_arr = permutedims(hω_arr) #  convert to (n2, ω)

	r_obs_ind = radial_grid_index(nobs1)

	@unpack p_Gsrc, p_Gobs1, Gfn_path_src, NGfn_files_src =
		unpackGfnparams(p_Gsrc, r_src, p_Gobs1, r_obs_default, p_Gobs2, r_obs_default)

	Gfn_path_obs, NGfn_files_obs = p_Gobs1.path, p_Gobs1.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	Gfn_fits_files_src, Gfn_fits_files_obs =
		Gfn_fits_files.((Gfn_path_src, Gfn_path_obs),
			((ℓ_arr, 1:Nν_Gfn),), (ℓ_ωind_iter_on_proc,),
			(NGfn_files_src, NGfn_files_obs))

	K = zeros(nr, length(nobs2_arr))

	arrs = allocatearrays(SoundSpeed(), los, true)
	@unpack Gsrc, drGsrc, divGsrc, divGobs = arrs
	Gobs, drGobs = arrs.Gobs1, arrs.drGobs1
	f_robs_rsrc, Hjj_robs_rsrc = arrs.fjₒjₛ_r₁_rsrc, arrs.Hjₒjₛω_r₁r₂

	Y12 = los_projected_biposh_spheroidal(computeY₀₀, nobs1, nobs2_arr, los, ℓ_ωind_iter_on_proc)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		conjhω = conj(@view hω_arr[:, ω_ind])

		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about the source radius
		read_Gfn_file_at_index!(drGsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:1, srcindFITS(los), :, 2)

		Gγr_robs_rsrc = αrcomp(Gsrc, r_obs_ind, los)

		# Green function about receiver location
		read_Gfn_file_at_index!(Gobs, Gfn_fits_files_obs,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about receiver location
		read_Gfn_file_at_index!(drGobs, Gfn_fits_files_obs,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs, 1:1, srcindFITS(los), :, 2)

		radial_fn_isotropic_δc_firstborn!(f_robs_rsrc,
			Gsrc, drGsrc, divGsrc, Gobs, drGobs, divGobs, ℓ)

		Hjₒjₛω!(Hjj_robs_rsrc, f_robs_rsrc, Gγr_robs_rsrc)

		populatekernelvalidation!(SoundSpeed(), los, K,
			(ℓ, ω), (@view Y12[:, ℓ]), conjhω, Hjj_robs_rsrc, Hjj_robs_rsrc, dω)
	end

	map(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs))
	return K
end

kernelfilenameδc₀₀(::TravelTimes, ::los_radial) = "K_δτ_δc₀₀.fits"
kernelfilenameδc₀₀(::TravelTimes, ::los_earth) = "K_δτ_δc₀₀_los.fits"
kernelfilenameδc₀₀(::Amplitudes, ::los_radial) = "K_A_δc₀₀.fits"
kernelfilenameδc₀₀(::Amplitudes, ::los_earth) = "K_A_δc₀₀_los.fits"

function kernel_δc₀₀(comm, m::SeismicMeasurement, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	r_src, = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs..., c_scale = 1)
	@unpack ν_arr, ℓ_arr = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)

	hω_arr = get(kwargs, :hω) do
		h = hω(comm, m, xobs1, xobs2, los; kwargs..., print_timings = false)
		_broadcast(h, 0, comm)
	end

	K = pmapsum(comm, kernel_δc₀₀_partial, iters,
		xobs1, xobs2, los, hω_arr, p_Gsrc, p_Gobs1, p_Gobs2, r_src)
	K === nothing && return nothing

	if get(kwargs, :save, true)
		filename = joinpath(SCRATCH_kerneldir[], kernelfilenameδc₀₀(m, los))
		FITS(filename,"w") do f
			write(f, reinterpret_as_float(K))
		end
	end

	return K
end

########################################################################################
# Sound-speed kernels
########################################################################################

function populatekernel!(::SoundSpeed, ::los_radial, K::StructArray{<:Complex}, SHModes, Yjₒjₛ_n1n2, Yjₒjₛ_n2n1,
	(l, m), (jₛ, jₒ, ω), pre, tworealconjhωHjₒjₛω_r₁r₂, tworealconjhωconjHjₒjₛω_r₂r₁)

	lm_ind = modeindex(firstshmodes(Yjₒjₛ_n1n2), (l, m))
	lm_ind_K = modeindex(SHModes, (l, m))

	_Yjₒjₛ_lm_n₁n₂ = Yjₒjₛ_n1n2[lm_ind]
	_Yjₒjₛ_lm_n₂n₁ = Yjₒjₛ_n2n1[lm_ind]

	@turbo for r_ind in eachindex(r)
		T₁₂ = pre * conj(_Yjₒjₛ_lm_n₁n₂) * tworealconjhωHjₒjₛω_r₁r₂[r_ind]
		T₂₁ = pre * conj(_Yjₒjₛ_lm_n₂n₁) * tworealconjhωconjHjₒjₛω_r₂r₁[r_ind]

		K[r_ind, lm_ind_K] += T₂₁ + T₁₂
	end

	return K
end

function populatekernel!(::SoundSpeed, ::los_earth, K::StructArray{<:Complex}, SHModes, Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
	(l, m), (jₛ, jₒ, ω), pre, tworealconjhωHjₒjₛω_r₁r₂, tworealconjhωconjHjₒjₛω_r₂r₁)

	lm_ind = modeindex(firstshmodes(Yjₒjₛ_lm_n₁n₂), (l, m))
	lm_ind_K = modeindex(SHModes, (l, m))

	_Yjₒjₛ_lm_n₁n₂ = parent(Yjₒjₛ_lm_n₁n₂[lm_ind])
	_Yjₒjₛ_lm_n₂n₁ = parent(Yjₒjₛ_lm_n₂n₁[lm_ind])

	@turbo for r_ind in eachindex(r)
		T₁₂ = pre * sum(conj(_Yjₒjₛ_lm_n₁n₂) .* tworealconjhωHjₒjₛω_r₁r₂[r_ind])
		T₂₁ = pre * sum(conj(_Yjₒjₛ_lm_n₂n₁) .* tworealconjhωconjHjₒjₛω_r₂r₁[r_ind])

		K[r_ind, lm_ind_K] += T₂₁ + T₁₂
	end
	return K
end

function kernel_δcₛₜ_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction, SHModes, hω_arr,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	@unpack p_Gsrc, p_Gobs1, p_Gobs2, Gfn_path_src, NGfn_files_src,
	Gfn_path_obs1, NGfn_files_obs1, Gfn_path_obs2, NGfn_files_obs2 =
		unpackGfnparams(p_Gsrc, r_src, p_Gobs1, radius(xobs1), p_Gobs2, radius(xobs2))

	@unpack ℓ_arr, ω_arr, Nν_Gfn, dω = p_Gsrc

	s_max = l_max(SHModes)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))
	obs_at_same_height = r₁_ind == r₂_ind

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs1 = Gfn_fits_files(Gfn_path_obs1,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs1)

	Gfn_fits_files_obs2 = Gfn_fits_files(Gfn_path_obs2,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, s_max, NGfn_files_obs2)

	arrs = allocatearrays(SoundSpeed(), los, obs_at_same_height)
	@unpack Gsrc, drGsrc, Gobs1, drGobs1, Gobs2, drGobs2, divGsrc, divGobs = arrs;
	@unpack fjₒjₛ_r₁_rsrc, fjₒjₛ_r₂_rsrc, Hjₒjₛω_r₁r₂, Hjₒjₛω_r₂r₁ = arrs;
	@unpack Gobs1_cache, drGobs1_cache, Gobs2_cache, drGobs2_cache = arrs;
	@unpack tworealconjhωHjₒjₛω_r₁r₂, tworealconjhωconjHjₒjₛω_r₂r₁ = arrs;
	tworealconjhωconjHjₒjₛω_r₂r₁P, tworealconjhωHjₒjₛω_r₁r₂P = map(parent, (tworealconjhωconjHjₒjₛω_r₂r₁, tworealconjhωHjₒjₛω_r₁r₂));

	HT = HybridArray{Tuple{sizeindG(los)...,sizeindG(los)...,StaticArrays.Dynamic()}}
	H₂₁ = _structarrayparent(Hjₒjₛω_r₂r₁, HT)
	H₁₂ = _structarrayparent(Hjₒjₛω_r₁r₂, HT)

	tworealconjhωconjHjₒjₛω_r₂r₁S = reinterpretSMatrix(tworealconjhωconjHjₒjₛω_r₂r₁P, los)
	tworealconjhωHjₒjₛω_r₁r₂S = reinterpretSMatrix(tworealconjhωHjₒjₛω_r₁r₂P, los)

	ℓ_range_proc = UnitRange(extremaelement(ℓ_ωind_iter_on_proc, dims = 1)...)

	jₒjₛ_allmodes = L2L1Triangle(ℓ_range_proc, s_max, ℓ_arr)
	jₒrange = l2_range(jₒjₛ_allmodes)
	@unpack Gobs1_cache, Gobs2_cache = arrs
	for jₒ in jₒrange
		Gobs1_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
		drGobs1_cache[jₒ] = zeros(ComplexF64, axes(drGsrc))
		if !obs_at_same_height
			Gobs2_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
			drGobs2_cache[jₒ] = zeros(ComplexF64, axes(drGsrc))
		end
	end

	@timeit localtimer "CG" begin
		C = SHVector([zeros(abs(jₒ-jₛ):jₒ+jₛ) for (jₒ, jₛ) in jₒjₛ_allmodes], jₒjₛ_allmodes)
		for (jₒ, jₛ) in jₒjₛ_allmodes
			C_jₒjₛ = C[(jₒ, jₛ)]
			for l in abs(jₒ-jₛ):jₒ+jₛ
				if isodd(jₒ+jₛ+l)
					continue
				end
				C_jₒjₛ[l] = clebschgordan(Float64, jₒ, 0, jₛ, 0, l, 0)
			end
		end
	end

	_K = zeros(nr, length(SHModes))
	K = StructArray{ComplexF64}((_K, zero(_K)))

	@timeit localtimer "biposh" begin
		Y12, Y21 = los_projected_spheroidal_biposh_flippoints(xobs1, xobs2, los, SHModes, jₒjₛ_allmodes)
	end

	ω_ind_prev = -1
	# Loop over the Greenfn files
	for (jₛ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]
		conjhω = conj(hω_arr[ω_ind])
		dωω²Pω = dω/2π * ω^2 * Powspec(ω)

		@timeit localtimer "FITS" begin

		# Green function about the source radius
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		Gαr_r₁_rsrc = αrcomp(Gsrc, r₁_ind, los)
		Gαr_r₂_rsrc = αrcomp(Gsrc, r₂_ind, los)

		# Derivative of Green function about the source radius
		read_Gfn_file_at_index!(drGsrc, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (jₛ, ω_ind), NGfn_files_src, 1:1, srcindFITS(los), :, 2)
		end #localtimer

		for jₒ in l2_range(jₒjₛ_allmodes, jₛ)

			jₒjₛ_ind = modeindex(jₒjₛ_allmodes, (jₒ, jₛ))

			Yjₒjₛ_lm_n₁n₂ = Y12[jₒjₛ_ind]
			Yjₒjₛ_lm_n₂n₁ = Y21[jₒjₛ_ind]

			(Yjₒjₛ_lm_n₁n₂ === nothing || Yjₒjₛ_lm_n₂n₁ === nothing) && continue

			Cljₒjₛ = C[(jₒ, jₛ)]

			Gobs1 = Gobs1_cache[jₒ]
			drGobs1 = drGobs1_cache[jₒ]
			Gobs2 = Gobs2_cache[jₒ]
			drGobs2 = drGobs2_cache[jₒ]

			@timeit localtimer "FITS" begin

			if ω_ind != ω_ind_prev || jₒ ∉ l2_range(jₒjₛ_allmodes, jₛ-1)
				# Green functions based at the observation point for ℓ′
				read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
				(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)

				read_Gfn_file_at_index!(drGobs1, Gfn_fits_files_obs1,
				(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs1, 1:1, srcindFITS(los), :, 2)
			end

			end # timer

			@timeit localtimer "radial term" begin

			# precompute the radial term in f
			radial_fn_δc_firstborn!(fjₒjₛ_r₁_rsrc, Gsrc, drGsrc, jₛ, divGsrc,
				Gobs1, drGobs1, jₒ, divGobs)

			end # timer

			@timeit localtimer "radial term 2" begin
				Hjₒjₛω!(Hjₒjₛω_r₁r₂, fjₒjₛ_r₁_rsrc, Gαr_r₂_rsrc)
			end
			@timeit localtimer "radial term 3" begin
				@turbo for I in eachindex(H₁₂)
					tworealconjhωHjₒjₛω_r₁r₂P[I] = 2real(conjhω * H₁₂[I])
				end
			end

			if !obs_at_same_height

				@timeit localtimer "FITS" begin
				if ω_ind != ω_ind_prev || jₒ ∉ l2_range(jₒjₛ_allmodes, jₛ-1)
					read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
					(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)

					read_Gfn_file_at_index!(drGobs2, Gfn_fits_files_obs2,
					(ℓ_arr, 1:Nν_Gfn), (jₒ, ω_ind), NGfn_files_obs2, 1:1, srcindFITS(los), :, 2)
				end
				end # timer

				@timeit localtimer "radial term" begin
				radial_fn_δc_firstborn!(fjₒjₛ_r₂_rsrc, Gsrc, drGsrc, jₛ, divGsrc,
					Gobs2, drGobs2, jₒ, divGobs)
				end # timer

				@timeit localtimer "radial term 2" begin
				Hjₒjₛω!(Hjₒjₛω_r₂r₁, fjₒjₛ_r₂_rsrc, Gαr_r₁_rsrc)
				end
			end
			@timeit localtimer "radial term 3" begin
				@turbo for I in eachindex(H₂₁)
					tworealconjhωconjHjₒjₛω_r₂r₁P[I] = 2real(conjhω * conj(H₂₁[I]))
				end
			end

			for (l, m) in firstshmodes(Yjₒjₛ_lm_n₁n₂)

				# The Clebsch-Gordan coefficients imply the selection
				# rule that only even ℓ+ℓ′+s modes contribute
				isodd(jₒ + jₛ + l) && continue

				# Check if triangle condition is satisfied
				# The loop is over all possible ℓ and ℓ′ for a given s_max
				# Not all combinations would contribute towards a specific s
				δ(jₛ, jₒ, l) || continue

				pre = dωω²Pω * Njₒjₛs(jₒ, jₛ, l) * Cljₒjₛ[l]

				@timeit localtimer "kernel" begin

				populatekernel!(SoundSpeed(), los, K, SHModes,
					Yjₒjₛ_lm_n₁n₂, Yjₒjₛ_lm_n₂n₁,
					(l, m), (jₛ, jₒ, ω), pre,
					tworealconjhωHjₒjₛω_r₁r₂S,
					tworealconjhωconjHjₒjₛω_r₂r₁S)
				end
			end # lm
		end # jₒ

		ω_ind_prev = ω_ind
	end # (jₛ, ω)

	map(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))
	Array(K)
end

function kernelfilenameδcₛₜ(m::SeismicMeasurement, los::los_direction, j_range, SHModes)
	mtag = m isa TravelTimes ? "δτ" : "A"
	lostag = los isa los_radial ? "" : "_los"
	"Kst_$(mtag)_c_$(modetag(j_range, SHModes))$(lostag).fits"
end

function kernel_δcₛₜ(comm, m::SeismicMeasurement, xobs1, xobs2,
	los::los_direction = los_radial(); kwargs...)

	SHModes = getkernelmodes(; kwargs...)

	r_src, = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs..., c_scale = 1)
	@unpack ν_arr, ℓ_arr = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)
	j_range, ν_ind_range = iters

	hω_arr = get(kwargs, :hω) do
		h = hω(comm, m, xobs1, xobs2, los; kwargs..., print_timings = false)
		_broadcast(h, 0, comm)
	end

	K_δcₛₜ = pmapsum(comm, kernel_δcₛₜ_partial, iters,
		xobs1, xobs2, los, SHModes, hω_arr, p_Gsrc, p_Gobs1, p_Gobs2, r_src)

	K_δcₛₜ === nothing && return nothing

	if get(kwargs, :save, true)
		filename = joinpath(SCRATCH_kerneldir[], kernelfilenameδcₛₜ(m, los, j_range, SHModes))
		header = generatefitsheader(xobs1, xobs2, SHModes, j_range, ν_ind_range)
		FITS(filename,"w") do f
			FITSIO.write(f, reinterpret_as_float(K_δcₛₜ), header = header)
		end
	end

	SHArray(K_δcₛₜ, (axes(K_δcₛₜ, 1), SHModes))
end


text:
module DataTransform

# For loading pipelines
using JSON # TODO: vs JSON3?
using YAML

# For manipulating data
using Query
using DataFrames
using FeatureEng

greet() = print("Hello World!")

end # module

ext: jl


text:
using Main.MyModule

f()

max_stars_repo_name: oist/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
module Charges

# Flat brokerage in Rs. per executed order. So a buy and sell will become Rs. 40
FLAT_BROKERAGE = 20.00
# in percent
BROKERAGE = 0.03 * 0.01

# STT (Securities Transaction Tax)
# in percent
STT = 0.025 * 0.01
STT_COMMODITY = 0.01 * 0.01

# Exchange Transaction Charges
# percent
EXCHANGE_TRANSACTION_CHARGES_EQUITY_INTRADAY = 0.00325 * 0.01
EXCHANGE_TRANSACTION_CHARGES_CURRENCY_FUTURES = 0.0009 * 0.01
EXCHANGE_TRANSACTION_CHARGES_COMMODITY_FUTURES = 0.0026 * 0.01

# percent
GST = 18 * 0.01

# SEBI charges Rs.5 / crore
SEBI_CHARGES = 1_00_00_000 / 10

# Stamp duty
# percent
STAMP_DUTY = 0.002 * 0.01
STAMP_DUTY_COMMODITY_FUTURES = 0.002 * 0.01

function commodity_futures(commodity, buy_price, sell_price, quantity)
    unit_change =
        if commodity == "crudeoil"
            100
        elseif commodity == "gold"
            100
        elseif commodity == "goldm"
            10
        elseif commodity == "naturalgas"
            1250
        elseif commodity == "silver"
            30
        elseif commodity == "silverm"
            5
        elseif commodity == "silvermic"
            1
        end

    quantity = quantity * unit_change
    turnover = (buy_price + sell_price) * quantity

    brokerage = turnover * BROKERAGE
    @debug brokerage

    brokerage =
      if brokerage > FLAT_BROKERAGE * 2
        FLAT_BROKERAGE * 2
      else
        brokerage
      end
    @debug brokerage

    stt = sell_price * quantity * STT_COMMODITY

    exchange_transaction_charges = turnover * EXCHANGE_TRANSACTION_CHARGES_COMMODITY_FUTURES

    gst = (brokerage + exchange_transaction_charges) * GST
    sebi_charges = turnover / SEBI_CHARGES
    stamp_duty = (buy_price * quantity) * STAMP_DUTY_COMMODITY_FUTURES

    total_charges = brokerage + stt + exchange_transaction_charges + gst + sebi_charges + stamp_duty

    gross_profit = (sell_price - buy_price) * quantity
    net_profit = gross_profit - total_charges

    points_to_breakeven = 0
    pips_to_breakeven = 0

    return Dict(
      "turnover" => turnover,
      "gross_profit" => gross_profit,
      "brokerage" => brokerage,
      "stt" => stt,
      "exchange_transaction_charges" => exchange_transaction_charges,
      "gst" => gst,
      "sebi_charges" => sebi_charges,
      "stamp_duty" => stamp_duty,
      "total_charges" => total_charges,
      "net_profit" => net_profit,
      "points_to_breakeven" => points_to_breakeven,
      "pips_to_breakeven" => pips_to_breakeven
    )
end

end


text:
import Distributions.MvNormal

"""
Forecast the state of the process at the next time step.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- x : AbstractMvNormal.  Current state estimate.
- u : Vector, optional.  External input.
- t : Real, optional. Time supplied to the time-varying model matrices (if the
    model is not time-varying, it will be ignored.)
#### Returns
- MvNormal distribution representing the forecast and its associated uncertainty.
"""
function predict( m :: AbstractGaussianSSM,
                  x :: AbstractMvNormal;
                  u :: Vector = zeros(m.nu),
                  t :: Real = 0.0)

    F  = process_matrix(m, x, t)
    CI = control_input(m, u, t)
    return MvNormal(F * mean(x) .+ CI, F * cov(x) * F' + m.V(t))
end

"""
Observe the state process with uncertainty.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- x : AbstractMvNormal.  Current state estimate.
- u : Vector, optional.  External input.
- t : Real, optional. Time supplied to the time-varying model matrices (if the
model is not time-varying, it will be ignored.)
#### Returns
- MvNormal distribution, representing the probability of recording any
particular observation data.
"""
function observe(m::AbstractGaussianSSM, x::AbstractMvNormal, t::Real=0.0)
    G = observation_matrix(m, x, t)
    return MvNormal(G * mean(x), G * cov(x) * G' + m.W(t))
end

function _filter(m      :: AbstractGaussianSSM,
                 y      :: Array{T},
                 x0     :: AbstractMvNormal,
                 u      :: Array{T},
                 times  :: Vector{T},
                 filter :: AbstractKalmanFilter) where T

    x_filtered = Array{AbstractMvNormal}(size(y, 2))
    loglik = 0.0
    x_pred = predict(m, x0, u=u[:, 1], t=times[1])
    x_filtered[1] = update(m, x_pred, y[:, 1], filter, 1)
    loglik = logpdf(x_filtered[1], mean(x_pred))

    if ! any(isnan(y[:, 1]))
        loglik += logpdf(observe(m, x_filtered[1]), y[:,1])
    end

    for i in 2:size(y, 2)
        x_pred = predict(m, x_filtered[i-1], u=u[:, i], t=times[i])
        # Check for missing values in observation
        if any(isnan(y[:, i]))
        x_filtered[i] = x_pred
        else
            x_filtered[i] = update(m, x_pred, y[:, i], filter, i)
            loglik += logpdf(observe(m, x_filtered[i]), y[:, i])
        end
        loglik += logpdf(x_pred, mean(x_filtered[i]))
    end

    return FilteredState(y, x_filtered, u, times, loglik, false)

end

"""
Given a process/observation model and a set of data, estimate the states of the
hidden process at each time step.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- y : Array of data.  Each column is a set of observations at one time, while
each row is one observed variable.  Can contain NaNs to represent missing data.
- x0 : AbstractMvNorm. Multivariate normal distribution representing our estimate
of the hidden state and our uncertainty about that estimate before our first observation.
(i.e., if the first data arrives at t=1, this is where we think the process is at t=0.)
- u : Matrix, optional.  Array of external inputs.
- t : Real, optional. Time supplied to the time-varying model matrices (if the
model is not time-varying, it will be ignored.)
#### Returns
- A FilteredState object holding the estimates of the hidden state.

#### Notes
This function only does forward-pass filtering--that is, the state estimate at time
t incorporates data from 1:t, but not from t+1:T.  For full forward-and-backward
filtering, run `smooth` on the FilteredState produced by this function.
"""
function filter( m      :: LinearGaussianSSM,
                 y      :: Array{T},
                 x0     :: AbstractMvNormal;
                 filter :: LinearKalmanFilter=KF(),
                 u      :: Matrix{T}=zeros(m.nu, size(y, 2)),
                 times  :: Vector{T}=zeros(size(y, 2))) where T

    return _filter(m, y, x0, u, times, filter)

end

function filter( m      :: NonlinearGaussianSSM,
                 y      :: Array{T},
                 x0     :: AbstractMvNormal;
                 filter :: NonlinearKalmanFilter=EKF(),
                 u      :: Matrix{T}=zeros(m.nu, size(y, 2)),
                 times  :: Vector{T}=zeros(size(y, 2))) where T

    return _filter(m, y, x0, u, times, filter)

end

"""
Given a process/observation model and a set of data, estimate the states of the
hidden process at each time step using forward and backward passes.

This function has two methods.  It can take either an AbstractGaussianSSM and a the
FilteredState produced by a previous call to `filter`, or it can take a model, data,
set, and initial state estimate, in which case it does the forward and backward
passes at once.

"""
function _smooth( m  :: AbstractGaussianSSM,
                  fs :: FilteredState{T}) where T

    n = size(fs.observations, 2)
    smooth_dist = Array(AbstractMvNormal, n)
    smooth_dist[end] = fs.state[end]
    if ! any(isnan(fs.observations[:, 1]))
        loglik = logpdf(observe(m, smooth_dist[end], fs.times[end]),
            fs.observations[:, end])
    end

    for i in (n - 1):-1:1
        state_pred = predict(m, fs.state[i], u=fs.input[:, i], t=fs.times[i])
        P = cov(fs.state[i])
        F = process_matrix(m, fs.state[i], fs.times[i])
        J = P * F' * inv(cov(state_pred))
        x_smooth = mean(fs.state[i]) + J *
            (mean(smooth_dist[i+1]) - mean(state_pred))
        P_smooth = P + J * (cov(smooth_dist[i+1]) - cov(state_pred)) * J'
        smooth_dist[i] = MvNormal(x_smooth, P_smooth)
        loglik += logpdf(predict(m, smooth_dist[i], u=fs.input[:, i], t=fs.times[i]),
            mean(smooth_dist[i+1]))
        if ! any(isnan(fs.observations[:, i]))
            loglik += logpdf(observe(m, smooth_dist[i], fs.times[i]), fs.observations[:, i])
        end
    end
    return FilteredState(fs.observations, smooth_dist, fs.input, fs.times, loglik, true)
end

function smooth(m::AbstractGaussianSSM, fs::FilteredState{T}) where T
    return _smooth(m, fs)
end

## Linear Gaussian method
function smooth( m      :: LinearGaussianSSM,
                 y      :: Array{T},
                 x0     :: AbstractMvNormal;
                 u      :: Matrix{T}=zeros(m.nu, size(y, 2)),
                 filter :: LinearKalmanFilter=KF()) where T

    fs = filter(m, y, x0, u=u, filter=filter)

    return _smooth(m, fs)

end

## Nonlinear Gaussian method
function smooth( m  :: NonlinearGaussianSSM,
                 y  :: Array{T},
                 x0 :: AbstractMvNormal;
                 u  :: Matrix{T} = zeros(m.nu, size(y, 2)),
                 filter::NonlinearKalmanFilter=EKF())  where T

    fs = filter(m, y, x0, u=u, filter=filter)

    return _smooth(m, fs)

end

"""
Simulate a realization of a state-space process and observations of it.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- x0 : AbstractMvNorm. Distribution of the state process's starting value.

#### Returns
- (x, y) : Tuple of the process and observation arrays.
"""
function simulate( m  :: AbstractGaussianSSM,
                   n  :: Int64,
                   x0 :: AbstractMvNormal)

    F = process_matrix(m, x0, 1)
    G = observation_matrix(m, x0, 1)
    x = zeros(size(F, 1), n)
    y = zeros(size(G, 1), n)
    x[:, 1] = rand(MvNormal(F * mean(x0), m.V(1)))
    y[:, 1] = rand(MvNormal(G * x[:, 1], m.W(1)))
    for i in 2:n
        F = process_matrix(m, x[:, i-1], i)
        G = observation_matrix(m, x[:, i-1], i)
        x[:, i] = F * x[:, i-1] + rand(MvNormal(m.V(i)))
        y[:, i] = G * x[:, i] + rand(MvNormal(m.W(i)))
    end
    return (x, y)

end

"""
Fit a state-space model to data using maximum likelihood.
"""
function fit(build_func, y, x0, params0, args...)

    function objective(params)
        ssm = build_func(params)
        return -loglikelihood(smooth(m, y, x0))
    end
    fit = optimize(objective, params0, args...)

end


text:
#TODO:
#handle batches
#gradient propagation line

softmax(x) = exp.(x) ./ sum(exp.(x), ndims(x))

"""
Sample from Gumbel(0,1)
"""
function sample_gumbel(shape, eps::Float64=1e-20)
    u = rand(shape)
    return -log.(-log.(u + eps) + eps)
end

"""
Draw a sample from the Gumbel-softmax distribution
"""
function gumbel_softmax_sample(logits, temperature::Float64)
    @assert temperature > 0.0
    y = logits + sample_gumbel(size(logits))
    return softmax(y ./ temperature)
end
using ReverseDiff
"""
Sample from the Gumbel-Softmax distribution and optionally discretize.
Args:
    logits: [batch_size, n_class] unnormalized log-probs
    temperature: non-negative scalar
    hard: if true, take argmax, but differential w.r.t. soft sample y
Returns:
    [batch_size, n_class] sample from the Gumbel-softmax distribution.
    If hard=true, then the returned sample will be one-hot, otherwise it
    will be a probability distribution that sums to 1 across classes.
"""
function gumbel_softmax(logits, temperature::Float64, hard::Bool=false)
    y = gumbel_softmax_sample(logits, temperature)
    if hard
        imax = indmax(y)
        y_hard = [i == imax ? 1.0 : 0.0 for i=1:length(y)]
        #= m = maximum(y) =#
        #= y_hard = [yy == m ? 1.0 : 0.0 for yy in y] =#
        #y_hard = zeros(Float64, length(logits))
        #y_hard[imax] += 1.0
        #= y_hard = zeros(Int, size(y)) =#
        #= for i = 1:length(imax) =#
        #=     y_hard[i, imax[i]] = 1 =#
        #= end =#
        #y_hard = [ j==imax[i] ? 1 : 0 for i=1:size(y,1), j=1:size(y,2) ]
        #return stop_gradient(y - y) + y
        return stop_gradient(y_hard - y) + y
    end
    return y
end

function indmaxcols(X)
    return [indmax(X[i,:]) for i=1:size(X,1)]
end

using ReverseDiff: TrackedReal, TrackedArray, ForwardOptimize
using ForwardDiff: Dual, value
stop_gradient(x) = x
@inline function stop_gradient{N}(n::Dual{N})
    Dual(value(n), 0.0)
end
stop_gradient(t::TrackedArray) = map(ForwardOptimize(stop_gradient), t)
stop_gradient{T<:TrackedReal}(t::Array{T}) = map(ForwardOptimize(stop_gradient), t)

#test
#ReverseDiff.gradient(x->sum(gumbel_softmax(x,0.3,true)),(logits,))
#ReverseDiff.gradient(x->sum(stop_gradient(ones(length(logits))-x)+x), (logits,))


text:
using DFTranspose
using Test

using Test, DataFrames, Random, CategoricalArrays, Dates, PooledArrays
const ≅ = isequal
@testset "general usage" begin

    df = DataFrame(x1 = [1,2,3,4], x2 = [1,4,9,16])
    dft =DataFrame(_variables_ =["x1","x2"], _c1 = [1,1], _c2 = [2,4],
                    _c3 = [3,9], _c4 = [4,16])
    @test df_transpose(df, [:x1, :x2]) == dft

    df = DataFrame(rand(10,5), :auto)
    df2 = df_transpose(df, r"x")

    @test Matrix(df) == permutedims(Matrix(df2[!,r"_c"]))

    df = DataFrame(rand(10,5), :auto)
    allowmissing!(df)
    df2 = df_transpose(df, r"x")

    @test df_transpose(df, r"x") == df_transpose(df, :)
    @test df_transpose(df, r"x") == df_transpose(df, [1,2,3,4,5])
    @test Matrix(df) == permutedims(Matrix(df2[!,2:end]))

    df = DataFrame(rand(10,5), :auto)
    allowmissing!(df)
    df[1, 4] = missing
    df[4, 5] = missing
    df2 = df_transpose(df, r"x")
    @test Matrix(df) ≅ permutedims(Matrix(df2[!,Not(:_variables_)]))

    df = DataFrames.hcat!(DataFrame(rand(Bool, 15,2),[:b1,:b2]), DataFrame(rand(15,5),:auto))
    df2 = df_transpose(df, :)
    @test Matrix(df) == permutedims(Matrix(df2[!,Not(:_variables_)]))

    df = DataFrame(foo = ["one", "one", "one", "two", "two","two"],
                    bar = ['A', 'B', 'C', 'A', 'B', 'C'],
                    baz = [1, 2, 3, 4, 5, 6],
                    zoo = ['x', 'y', 'z', 'q', 'w', 't'])
    df2 = df_transpose(df, :baz, :foo, id = :bar)
    df3 = DataFrame(foo = ["one", "two"], _variables_ = ["baz", "baz"], A = [1, 4], B = [2, 5], C = [3, 6])

    @test df_transpose(df, :baz, :foo, id = :bar) == df_transpose(df, [:baz], :foo, id = :bar)
    @test df_transpose(df, :baz, :foo, id = :bar) == df_transpose(df, [:baz], [:foo], id = :bar)
    @test df_transpose(df, :baz, :foo, id = :bar) == df_transpose(df, :baz, [:foo], id = :bar)
    @test df2 == df3

    df = DataFrame(id = [1, 2, 3, 1], x1 = rand(4), x2 = rand(4))
    @test_throws AssertionError df_transpose(df, r"x", id = :id)

    df = DataFrame(rand(1:100, 1000, 10), :auto)
    insertcols!(df, 1, :g => repeat(1:100, inner = 10))
    insertcols!(df, 2, :id => repeat(1:10, 100))
    # duplicate and id within the last group
    df[1000, :id] = 1
    @test_throws AssertionError df_transpose(df, r"x", [:g], id = :id )

    df = DataFrame(rand(1000, 100), :auto)
    mdf = Matrix(df)

    @test describe(df_transpose(df, r"x"), sum => :sum, cols = r"_c").sum ≈ sum(mdf, dims = 2)

    df = DataFrame([[1, 2], [1.1, 2.0],[1.1, 2.1],[1.1, 2.0]]
                    ,[:person, Symbol("11/2020"), Symbol("12/2020"), Symbol("1/2021")])
    dft = df_transpose(df, Not(:person), :person,
                        renamerowid = x -> Date(x, dateformat"m/y"),
                        variable_name = "Date",
                         renamecolid = x -> "measurement")
    dftm = DataFrame(person = [1,1,1,2,2,2],
                Date = Date.(repeat(["2020-11-01","2020-12-01","2021-01-01"], 2)),
                measurement = [1.1, 1.1, 1.1, 2.0, 2.1, 2.0])
    @test dft == dftm

    df = DataFrame(rand(100,5),:auto)
    insertcols!(df, 1, :id => repeat(1:20, 5))
    insertcols!(df, 1, :g => repeat(1:5, inner = 20))

    dft = df_transpose(df, r"x", [:g, :id], variable_name = "variable", renamecolid = x -> "value")
    dfs = sort(stack(df, r"x"), [:g, :id])

    @test dft == dfs

    df = DataFrame(group = repeat(1:3, inner = 2),
                                 b = repeat(1:2, inner = 3),
                                 c = repeat(1:1, inner = 6),
                                 d = repeat(1:6, inner = 1),
                                 e = string.('a':'f'))
    allowmissing!(df)
    df[2, :b] = missing
    df[4, :d] = missing
    df2 = df_transpose(df, 2:4, [:group], id = :e, default_fill = 0)

    df3 = DataFrame(group = repeat(1:3, inner = 3),
                    _variables_ = string.(repeat('b':'d', 3)),
                    a = [1,1,1,0,0,0,0,0,0],
                    b = [missing,1,2,0,0,0,0,0,0],
                    c = [0,0,0, 1,1,3, 0,0,0],
                    d = [0,0,0, 2,1,missing, 0,0,0],
                    e = [0,0,0,0,0,0,2,1,5],
                    f = [0,0,0,0,0,0,2,1,6] )
    @test df2 ≅ df3

    df = DataFrame(id = ["r3", "r1", "r2" , "r4"], x1 = [1,2,3,4], x2 = [1,4,9,16])
    df2 = df_transpose(df, [:x1,:x2], id = :id)
    df3 = DataFrame([["x1", "x2"],
                        [1, 1],
                        [2, 4],
                        [3, 9],
                        [4, 16]],
                        [:_variables_, :r3, :r1, :r2, :r4])
    @test df3 == df2

    pop = DataFrame(country = ["c1","c1","c2","c2","c3","c3"],
                            sex = ["male", "female", "male", "female", "male", "female"],
                            pop_2000 = [100, 120, 150, 155, 170, 190],
                            pop_2010 = [110, 120, 155, 160, 178, 200],
                            pop_2020 = [115, 130, 161, 165, 180, 203])

    popt = df_transpose(pop, r"pop_", :country,
                            id = :sex, variable_name = "year",
                            renamerowid = x -> match(r"[0-9]+",x).match, renamecolid = x -> x * "_pop")
    poptm = DataFrame([ ["c1", "c1", "c1", "c2", "c2", "c2", "c3", "c3", "c3"],
            SubString{String}["2000", "2010", "2020", "2000", "2010", "2020", "2000", "2010", "2020"],
            Union{Missing, Int64}[100, 110, 115, 150, 155, 161, 170, 178, 180],
            Union{Missing, Int64}[120, 120, 130, 155, 160, 165, 190, 200, 203]],
            [:country,:year,:male_pop,:female_pop])

    @test popt == poptm
    popt = df_transpose(pop, r"pop", r"cou", id = r"sex",  variable_name = "year",
                            renamerowid = x -> match(r"[0-9]+",x).match, renamecolid = x -> x * "_pop")
    @test popt == poptm
    pop.country = PooledArray(pop.country)
    popt = df_transpose(pop, r"pop", r"cou", id = r"sex",  variable_name = "year",
                            renamerowid = x -> match(r"[0-9]+",x).match, renamecolid = x -> x * "_pop")
    @test popt.country == PooledArray(["c1", "c1", "c1", "c2", "c2", "c2", "c3", "c3", "c3"])
    df =  DataFrame(region = repeat(["North","North","South","South"],2),
                 fuel_type = repeat(["gas","coal"],4),
                 load = [.1,.2,.5,.1,6.,4.3,.1,6.],
                 time = [1,1,1,1,2,2,2,2],
                 )

    df2 = df_transpose(df, :load, :time, id = 1:2)
    df3 = DataFrame([[1, 2],
         ["load", "load"],
         Union{Missing, Float64}[0.1, 6.0],
         Union{Missing, Float64}[0.2, 4.3],
         Union{Missing, Float64}[0.5, 0.1],
         Union{Missing, Float64}[0.1, 6.0]], [:time,:_variables_,Symbol("(\"North\", \"gas\")"),Symbol("(\"North\", \"coal\")"),Symbol("(\"South\", \"gas\")"),Symbol("(\"South\", \"coal\")")])

     @test df2 == df3
     df = DataFrame(A_2018=1:4, A_2019=5:8, B_2017=9:12,
                             B_2018=9:12, B_2019 = [missing,13,14,15],
                              ID = [1,2,3,4])
      f(x) =  match(r"[0-9]+",x).match
      dfA = df_transpose(df, r"A", :ID, renamerowid = f, variable_name = "Year", renamecolid = x->"A");
      dfB = df_transpose(df, r"B", :ID, renamerowid = f, variable_name = "Year", renamecolid = x->"B");
      df2 = outerjoin(dfA, dfB, on = [:ID, :Year])
      df3 = DataFrame([[1, 1, 2, 2, 3, 3, 4, 4, 1, 2, 3, 4],
                 SubString{String}["2018", "2019", "2018", "2019", "2018", "2019", "2018", "2019", "2017", "2017", "2017", "2017"],
                 Union{Missing, Int64}[1, 5, 2, 6, 3, 7, 4, 8, missing, missing, missing, missing],
                 Union{Missing, Int64}[9, missing, 10, 13, 11, 14, 12, 15, 9, 10, 11, 12]], [:ID,:Year,:A,:B])
        @test df2 ≅ df3
        df = DataFrame(paddockId= [0, 0, 1, 1, 2, 2],
                                color= ["red", "blue", "red", "blue", "red", "blue"],
                                count= [3, 4, 3, 4, 3, 4],
                                weight= [0.2, 0.3, 0.2, 0.3, 0.2, 0.2])
        df2 = df_transpose( df_transpose(df, [:count,:weight],[:paddockId,:color]),
                             :_c1, :paddockId, id = 2:3)
         df3 = DataFrame([[0, 1, 2],
                         ["_c1", "_c1", "_c1"],
                         Union{Missing, Float64}[3.0, 3.0, 3.0],
                         Union{Missing, Float64}[0.2, 0.2, 0.2],
                         Union{Missing, Float64}[4.0, 4.0, 4.0],
                         Union{Missing, Float64}[0.3, 0.3, 0.2]],[ :paddockId,:_variables_,Symbol("(\"red\", \"count\")"),Symbol("(\"red\", \"weight\")"),Symbol("(\"blue\", \"count\")"),Symbol("(\"blue\", \"weight\")")])

         @test df2 == df3

        df = DataFrame(x1 = [9,2,8,6,8], x2 = [8,1,6,2,3], x3 = [6,5,3,10,8])
        df2 = df_transpose(df, r"x", renamerowid = x -> match(r"[0-9]+",x).match,renamecolid = x -> "_column_" * string(x))
        df3 = DataFrame([SubString{String}["1", "2", "3"],
                         [9, 8, 6],
                         [2, 1, 5],
                         [8, 6, 3],
                         [6, 2, 10],
                         [8, 3, 8]],[:_variables_,:_column_1,:_column_2,:_column_3,:_column_4,:_column_5])
         @test df2 == df3
         df = DataFrame(a=["x", "y"], b=[1, "two"], c=[3, 4], d=[true, false])
         df2 = df_transpose(df, [:b, :c, :d], id = :a, variable_name = "new_col")
         df3 = DataFrame([["b", "c", "d"],
                         Any[1, 3, true],
                         Any["two", 4, false]],[:new_col,:x,:y])
         @test df2 == df3
         df = DataFrame(rand(1:100,100,4),:auto)
         DataFrames.hcat!(df, DataFrame(rand(100,3),:auto), makeunique = true)
         insertcols!(df,1, :id=>Symbol.(100:-1:1))

         df2 = df_transpose(df, r"x", id = :id, variable_name = "id")
         df3 = permutedims(df, :id)

         @test df2 == df3

end

@testset "Outputs - from DataFrames test set" begin
    df = DataFrame(Fish = CategoricalArray{Union{String, Missing}}(["Bob", "Bob", "Batman", "Batman"]),
                   Key = CategoricalArray{Union{String, Missing}}(["Mass", "Color", "Mass", "Color"]),
                   Value = Union{String, Missing}["12 g", "Red", "18 g", "Grey"])
    levels!(df[!, 1], ["XXX", "Bob", "Batman"])
    levels!(df[!, 2], ["YYY", "Color", "Mass"])
#     Not sure if it is relevant, however, we are doing it here
    df2 = df_transpose(df, [:Value], [:Fish], id = :Key)
    @test levels(df[!, 1]) == ["XXX", "Bob", "Batman"] # make sure we did not mess df[!, 1] levels
    @test levels(df[!, 2]) == ["YYY", "Color", "Mass"] # make sure we did not mess df[!, 2] levels
#     duplicates
    @test_throws AssertionError df_transpose(df, [:Value], id = :Key)

    df = DataFrame(Fish = CategoricalArray{Union{String, Missing}}(["Bob", "Bob", "Batman", "Batman"]),
                   Key = CategoricalArray{Union{String, Missing}}(["Mass", "Color", "Mass", "Color"]),
                   Value = Union{String, Missing}["12 g", "Red", "18 g", "Grey"])
    levels!(df[!, 1], ["XXX", "Bob", "Batman"])
    levels!(df[!, 2], ["YYY", "Color", "Mass"])
    df2 = df_transpose(df, [:Value], [:Fish], id = :Key, renamecolid=x->string("_", uppercase(string(x)), "_"))
    df4 = DataFrame(Fish = Union{String, Missing}["Bob", "Batman"],
                _variables_ = String["Value", "Value"],
                _COLOR_ = Union{String, Missing}["Red", "Grey"],
                _MASS_ = Union{String, Missing}["12 g", "18 g"])
    @test df2 == df4
    # without categorical array
    df = DataFrame(Fish = ["Bob", "Bob", "Batman", "Batman"],
                   Key = ["Mass", "Color", "Mass", "Color"],
                   Value = ["12 g", "Red", "18 g", "Grey"])
    df2 = df_transpose(df, [:Value], [:Fish], id = :Key)
    df4 = DataFrame(Fish = ["Bob", "Batman"],
                    _variables_ = ["Value", "Value"],
                    Mass = ["12 g", "18 g"],
                    Color = ["Red", "Grey"])
    @test df2 ≅ df4
    @test typeof(df2[!, :Fish]) <: Vector{String}
    #Make sure df_transpose works with missing values at the start of the value column
    allowmissing!(df, :Value)
    df[1, :Value] = missing
    df2 = df_transpose(df, [:Value], [:Fish], id =  :Key)
    #This changes the expected result
    allowmissing!(df4, :Mass)
    df4[1, :Mass] = missing
    @test df2 ≅ df4

    # test missing value in grouping variable
    mdf = DataFrame(RowID = 1:4, id=[missing, 1, 2, 3], a=1:4, b=1:4)
    @test select(df_transpose(df_transpose(mdf, [:a,:b], [:RowID, :id]), [:_c1], [:RowID, :id], id = :_variables_), :RowID, :id, :a, :b) ≅ mdf
    @test select(df_transpose(df_transpose(mdf, Not(1,2), [:RowID, :id]), [:_c1], [:RowID, :id], id = :_variables_), :RowID, :id, :a, :b) ≅ mdf

    # test more than one grouping column
    wide = DataFrame(id = 1:12,
                     a  = repeat([1:3;], inner = [4]),
                     b  = repeat([1:4;], inner = [3]),
                     c  = randn(12),
                     d  = randn(12))
    w2 = wide[:, [1, 2, 4, 5]]
    rename!(w2, [:id, :a, :_C_, :_D_])
    long = df_transpose(wide, [:c, :d], [:id, :a, :b])
    wide3 = df_transpose(long, [:_c1], [:id, :a, :b], id = :_variables_)
    @test select(wide3, Not(:_variables_)) == wide
    df = DataFrame([repeat(1:2, inner=4), repeat('a':'d', outer=2), collect(1:8)],
                       [:id, :variable, :value])

    udf = df_transpose(df, [:value], [:id], id = :variable)
    @test select(udf, Not(:_variables_)) == DataFrame([Union{Int, Missing}[1, 2], Union{Int, Missing}[1, 5],
                                Union{Int, Missing}[2, 6], Union{Int, Missing}[3, 7],
                                Union{Int, Missing}[4, 8]], [:id, :a, :b, :c, :d])

    @test isa(udf[!, 1], Vector{Int})
    @test isa(udf[!,:a], Vector{Union{Int,Missing}})
    @test isa(udf[!,:b], Vector{Union{Int,Missing}})
    @test isa(udf[!,:c], Vector{Union{Int,Missing}})
    @test isa(udf[!,:d], Vector{Union{Int,Missing}})
    df = DataFrame([categorical(repeat(1:2, inner=4)),
                           categorical(repeat('a':'d', outer=2)), categorical(1:8)],
                       [:id, :variable, :value])

    udf = df_transpose(df, [:value], [:id], id = :variable)
    @test isa(udf[!, 1], CategoricalArray{Int64,1,UInt32})
    @test isa(udf[!, :a], Vector{Union{Missing, CategoricalValue{Int64, UInt32}}})
    @test isa(udf[!, :b], Vector{Union{Missing, CategoricalValue{Int64, UInt32}}})
    @test isa(udf[!, :c], Vector{Union{Missing, CategoricalValue{Int64, UInt32}}})
    @test isa(udf[!, :d], Vector{Union{Missing, CategoricalValue{Int64, UInt32}}})

    df1 = DataFrame(a=["x", "y"], b=rand(2), c=[1, 2], d=rand(Bool, 2))

    @test_throws MethodError df_transpose(df1)
    @test_throws ArgumentError df_transpose(df1, :bar)

    df1_pd = df_transpose(df1, 2:4, id = 1)
    @test size(df1_pd, 1) == ncol(df1) - 1
    @test size(df1_pd, 2) == nrow(df1) + 1
    @test names(df1_pd) == ["_variables_", "x", "y"]
    @test names(df_transpose(df1, 2:4, id = 1, variable_name = "foo")) == ["foo", "x", "y"]

end


text:
using Pkg; Pkg.activate(@__DIR__)

using Tar, Downloads, Dates, ArtifactUtils

url = "https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.0/cytoscape.min.js"

dir = mkpath(joinpath(@__DIR__, "artifacts"))

#-----------------------------------------------------------------------------# start from scratch
rm(dir, recursive=true, force=true)
rm(joinpath(@__DIR__, "cytoscape_jl.tar"), force=true)
rm(joinpath(@__DIR__, "cytoscape_jl.tar.gz"), force=true)

#-----------------------------------------------------------------------------# download
Downloads.download(url, joinpath(mkdir(dir), "cytoscape.min.js"))

tarfile = joinpath(@__DIR__, "cytoscape_jl.tar")

Tar.create(dir, tarfile)

run(`gzip $tarfile`)

#-----------------------------------------------------------------------------# upload
try
    artifacts_today = "artifacts_$(today())"

    run(`gh release create $artifacts_today $(tarfile * ".gz") --title $artifacts_today --notes ""`)

    @info "Sleeping so artifact is ready on GitHub..."
    sleep(10)
    add_artifact!(
        "Artifacts.toml",
        "cytoscape.min.js",
        "https://github.com/joshday/CytoscapeJS.jl/releases/download/$artifacts_today/cytoscape_jl.tar.gz",
        force=true,
    )
catch ex
    @error "Error (probably the release already exists): $ex"
end