text:
module Plots

using CSV
using DataFrames
using DelimitedFiles
using Distributed
using Interpolations
using LaTeXStrings
using LinearAlgebra
using FFTW
using FITSIO
using JLD2
using LsqFit
using NPZ
using OffsetArrays
using Polynomials
using Printf
using ProgressMeter
using PyCall
using PyPlot
using SphericalHarmonicModes
using SphericalHarmonicArrays

@everywhere begin
	using ParallelUtilities
	using DistributedArrays
end

pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)

const ticker = PyNULL()
const colors = PyNULL()
const cm = PyNULL()
const axes_grid1 = PyNULL()

function __init__()
	copy!(ticker,  pyimport("matplotlib.ticker"))
	copy!(colors, pyimport("matplotlib.colors"))
	copy!(cm, pyimport("matplotlib.cm"))
	copy!(axes_grid1, pyimport("mpl_toolkits.axes_grid1"))
end
# ho_endpt_so = pyimport("ho_endpt_so")

using Main.HelioseismicKernels: SCRATCH_kerneldir, SCRATCH
using Main.HelioseismicKernels
using Main.HelioseismicKernels: Rsun, r, c, ρ, nr, r_src_default, r_obs_default
using Main.HelioseismicKernels: SCRATCH, SCRATCH_kerneldir
using Main.HelioseismicKernels: Point2D, Point3D
const r_frac = r./Rsun;

import Main.HelioseismicKernels: los_radial, los_earth

los_tag_string(::los_radial) = ""
los_tag_string(::los_earth) = "_los"

"""
	scipercentage(x)

Return the number as a percentage in scientific notation
"""
function scipercentage(x, xmax=100)
	# x = a * 10^b
	# log10(x) = log10(a) + b
	if iszero(x)
		return "0%"
	end
	b = floor(Int, log10(abs(x)))
	a = x / 10.0^b
	a *= 100/xmax
	string(a)*L"\times10^{%$b}"*"%"
end

function text_from_offsettext(ax)
	t = ax.get_offset_text().get_text()
	if isempty(t)
		return ""
	end
	s = split(t, 'e')
	m = parse(Int, s[1])
	ex = parse(Int, replace(s[2], '−'=>'-'))
	txt = L"×10^{%$ex}"
end

function replace_yaxis_offsettext(ax, y=1.02)
	t = text_from_offsettext(ax.yaxis)
	ax.yaxis.get_offset_text().set_visible(false)
	ax.text(0, y, t, transform=ax.transAxes, fontsize=10)
end

# traveltimes

function plot_δτ_validation_v(;kwargs...)

	r_src = get(kwargs,:r_src,r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") Nt dt dν Nν dω

	filename = joinpath(SCRATCH_kerneldir[],"δτ_v_CmC0.fits")
	if isfile(filename)
		δτ_v_CmC0 = FITS(filename,"r") do f
							read(f[1])
						end
	else
		throw(DimensionMismatch("Run δτ_v_compare before plotting"))
	end

	header = FITS(filename,"r") do f
				read_header(f[1])
			end

	ϕ_low = header["PHI_LOW"]
	ϕ_high = header["PHI_HIGH"]
	nϕ = header["NPHI"]

	filename = joinpath(SCRATCH_kerneldir[],"δτ_v_FB.fits")
	if isfile(filename)
		δτ_v_FB = FITS(filename,"r") do f
					read(f[1])
				end
	else
		throw(DimensionMismatch("Run δτ_v_compare before plotting"))
	end

	ϕ2_deg = collect(LinRange(ϕ_low,ϕ_high,nϕ))
	ϕ2_arr = deg2rad.(ϕ2_deg)
	n1 = Point2D(π/2,0)
	n2_arr = [Point2D(π/2,ϕ2) for ϕ2 in ϕ2_arr]

	percentage_diff = (δτ_v_CmC0-δτ_v_FB)./δτ_v_CmC0.*100

	ϕ2_arr = rad2deg.(ϕ2_arr)

	ax_δτ = get(kwargs,:ax_δτ) do
		subplot2grid((3,1),(0,0),rowspan=2)
	end
	ax_δτ.plot(ϕ2_arr,δτ_v_CmC0,label="C - C₀",color="black",ls="dashed")
	ax_δτ.plot(ϕ2_arr,δτ_v_FB,"^",label="First Born",ms=6,ls="None",color="black")
	ax_δτ.set_ylabel("δτ [sec]",fontsize=15)
	ax_δτ.legend(loc="best", fontsize=14)
	ax_δτ.xaxis.set_major_formatter(ticker.NullFormatter())
	ax_δτ.xaxis.set_major_locator(ticker.MaxNLocator(5))
	ax_δτ.yaxis.set_major_locator(ticker.MaxNLocator(4))

	ax_δτ_diff = get(kwargs,:ax_δτ_diff) do
		subplot2grid((3,1),(2,0))
	end
	ax_δτ_diff.plot(ϕ2_arr,percentage_diff,"o-",ms=4,zorder=2,color="black")
	ax_δτ_diff.set_ylabel("% Diff.",fontsize=15)
	ax_δτ_diff.set_xlabel("Angular separation [degrees]",fontsize=15)
	ax_δτ_diff.axhline(0,ls="dotted",color="black",zorder=0)
	ax_δτ_diff.xaxis.set_major_locator(ticker.MaxNLocator(5))
	ax_δτ_diff.yaxis.set_major_locator(ticker.MaxNLocator(3))
	ax_δτ_diff.margins(y=0.2)
	ax_δτ_diff.set_yticklabels([(@sprintf "%.2f" x)*"%"
				for x in ax_δτ_diff.get_yticks()])

	gcf().set_size_inches(6,4)
	tight_layout()

	dτ_arr = DataFrame(dist=ϕ2_deg,dt_FB=δτ_v_FB,dt_C=δτ_v_CmC0,
		percentage_diff=percentage_diff)

	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"dt_v.eps"))
	save && CSV.write(joinpath(SCRATCH[],"dt_v"),dτ_arr,delim=' ')

	dτ_arr
end

function plot_δτ_validation_δc(los = los_radial();kwargs...)

	r_src = get(kwargs,:r_src,r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") Nt dt dν Nν dω

	los_tag = los_tag_string(los)

	filename = joinpath(SCRATCH_kerneldir[],"δτ_δc_CmC0$los_tag.fits")
	if isfile(filename)
		δτ_δc_CmC0 = FITS(filename,"r") do f
							read(f[1])
						end
	else
		throw(DimensionMismatch("Run δτ_δc_compare before plotting"))
	end

	header = FITS(filename,"r") do f
				read_header(f[1])
			end

	ϕ_low = header["PHI_LOW"]
	ϕ_high = header["PHI_HIGH"]

	filename = joinpath(SCRATCH_kerneldir[],"δτ_δc_FB$los_tag.fits")
	if isfile(filename)
		δτ_δc_FB = FITS(filename,"r") do f
					read(f[1])
				end
	else
		throw(DimensionMismatch("Run δτ_δc_compare before plotting"))
	end

	nϕ=get(kwargs,:nϕ,length(δτ_δc_CmC0))

	ϕ2_deg = collect(LinRange(ϕ_low,ϕ_high,nϕ))
	ϕ2_arr = deg2rad.(ϕ2_deg)
	n1 = Point2D(π/2,0)
	n2_arr = [Point2D(π/2,ϕ2) for ϕ2 in ϕ2_arr]

	percentage_diff = (δτ_δc_CmC0-δτ_δc_FB)./δτ_δc_CmC0.*100

	ϕ2_arr = rad2deg.(ϕ2_arr)

	ax_δτ = get(kwargs,:ax_δτ) do
		subplot2grid((3,1),(0,0),rowspan=2)
	end
	ax_δτ.plot(ϕ2_arr,δτ_δc_CmC0,label="C - C₀",color="black",ls="dashed")
	ax_δτ.plot(ϕ2_arr,δτ_δc_FB,"^",label="First Born",ms=6,ls="None",color="black")
	ax_δτ.set_ylabel("Travel time shift [sec]",fontsize=12)
	ax_δτ.legend(loc="best")
	ax_δτ.xaxis.set_major_formatter(ticker.NullFormatter())
	ax_δτ.xaxis.set_major_locator(ticker.MaxNLocator(5))

	ax_δτ_diff = get(kwargs,:ax_δτ_diff) do
		subplot2grid((3,1),(2,0))
	end
	ax_δτ_diff.plot(ϕ2_arr,percentage_diff,"o-",ms=4,zorder=2,color="black")
	ax_δτ_diff.set_ylabel("Difference",fontsize=12)
	ax_δτ_diff.set_xlabel("Angular separation [degrees]",fontsize=12)
	ax_δτ_diff.axhline(0,ls="dotted",color="black",zorder=0)
	ax_δτ_diff.xaxis.set_major_locator(ticker.MaxNLocator(5))
	ax_δτ_diff.yaxis.set_major_locator(ticker.MaxNLocator(3))
	ax_δτ_diff.margins(y=0.2)
	ax_δτ_diff.set_yticklabels([(@sprintf "%.2f" x)*"%"
				for x in ax_δτ_diff.get_yticks()])

	tight_layout()

	dτ_arr = DataFrame(dist=ϕ2_deg,dt_FB=δτ_δc_FB,dt_C=δτ_δc_CmC0,
		percentage_diff=percentage_diff)

	save = get(kwargs,:save,false)
	save && savefig(joinpath(SCRATCH[],"dt_dc$los_tag.eps"))
	save && CSV.write(joinpath(SCRATCH[],"dt_dc$los_tag"),dτ_arr,delim=' ')

	dτ_arr
end

function plot_δτ_validation_δc_with_and_without_los(;kwargs...)
	ax_δτ_radial = subplot2grid((3,2),(0,0),rowspan=2)
	ax_δτ_radial_diff = subplot2grid((3,2),(2,0))
	ax_δτ_los = subplot2grid((3,2),(0,1),rowspan=2)
	ax_δτ_los_diff = subplot2grid((3,2),(2,1))

	plot_δτ_validation_δc(los_radial();kwargs...,
		ax_δτ=ax_δτ_radial,
		ax_δτ_diff=ax_δτ_radial_diff,
		save=false)

	plot_δτ_validation_δc(los_earth();kwargs...,
		ax_δτ=ax_δτ_los,
		ax_δτ_diff=ax_δτ_los_diff,
		save=false)

	ax_δτ_radial.set_title("Radial")
	ax_δτ_los.set_title("Line of sight")

	gcf().set_size_inches(8,4)
	tight_layout()
	savefig(joinpath(SCRATCH[],"dt_dc_radial_and_los.eps"))
end

# HelioseismicKernels

function plot_time_distance_with_and_without_los(;kwargs...) # time in hours

	r_src=get(kwargs,:r_src,r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") dt Nt

	t_max = get(kwargs,:t_max,5)
	t_min = get(kwargs,:t_min,1)
	ℓ_range = get(kwargs,:ℓ_range,20:100)
	nϕ = get(kwargs,:nϕ,maximum(ℓ_range))

	t_max_ind = floor(Int,t_max*60^2/dt) + 1
	t_min_ind = floor(Int,t_min*60^2/dt) + 1
	t_inds = t_min_ind:t_max_ind

	t = t_inds.*dt./60^2
	ϕ = LinRange(5,90,nϕ).*π/180

	nobs1 = Point2D(π/2,0)

	filename = joinpath(SCRATCH_kerneldir[],"C_t_phi.fits")
	if isfile(filename)
		C =	FITS(filename,"r") do f
				read(f[1])
			end
		if size(C,1) == Nt
			C = C[t_inds,:]
		end
	else
		error("file not found")
	end

	filename = joinpath(SCRATCH_kerneldir[],"Ct.fits")
	if isfile(filename)
		C60 =	FITS(filename,"r") do f
					read(f[1])
				end
		if size(C60,1) == Nt
			C60 = C60[t_inds]
		end
	else
		error("file not found")
	end

	filename = joinpath(SCRATCH_kerneldir[],"Ct_los.fits")
	if isfile(filename)
		C60_los =	FITS(filename,"r") do f
						read(f[1])
					end
		if size(C60_los,1) == Nt
			C60_los = C60_los[t_inds]
		end
	else
		error("file not found")
	end

	ax1 = get(kwargs,:ax1) do
		subplot2grid((1,3),(0,0),colspan=2)
	end
	ax2 = get(kwargs,:ax2) do
		subplot2grid((1,3),(0,2),colspan=1)
	end

	C_norm = maximum(abs,C)
	ax1.pcolormesh(ϕ.*180/π,t,C./C_norm,
		cmap="Greys",vmax=0.3,vmin=-0.3,rasterized=true)
	ax1.set_xlabel("Angular separation [degrees]",fontsize=12);
	ax1.set_ylabel("Time [hours]",fontsize=12);
	ax1.axvline(60,color="black",ls="dashed",lw=0.6)
	ax1.set_title("Time-distance diagram",fontsize=12);

	ax2.plot(C60,t,lw=0.7,color="grey",label="Radial")
	ax2.plot(C60_los,t,lw=0.7,color="black",label="LoS")
	ax2.yaxis.set_major_formatter(ticker.NullFormatter())
	ax2.set_ylim(ax1.get_ylim())
	ax2.set_ylabel("Time [hours]",fontsize=12);
	ax2.legend(loc="lower left")

	tight_layout()

	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"C_timedistance.eps"))
end

function plot_time_distance_bounce_filtered(;kwargs...)
	r_src = get(kwargs,:r_src,HelioseismicKernels.r_src_default)
	t_max = get(kwargs,:t_max,6)
	t_min = get(kwargs,:t_min,0.3)
	ℓ_range = get(kwargs,:ℓ_range,20:100)
	nϕ = get(kwargs,:nϕ,maximum(ℓ_range))
	bounce_no = get(kwargs,:bounce_no,1)

	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") dt Nt

	t_max_ind = floor(Int,t_max*60^2/dt)+1
	t_min_ind = floor(Int,t_min*60^2/dt)+1
	t_inds = t_min_ind:t_max_ind

	t = t_inds.*dt./60^2

	ϕ = LinRange(5,90,nϕ).*π/180

	filename = joinpath(SCRATCH_kerneldir[], "C_t_phi.fits")
	if isfile(filename)
		C = FITS(filename,"r") do f
				read(f[1])
			end
		C = C[t_inds, :]
	else
		error("File not found")
	end

	C_filt = copy(C)
	f_t = zeros(size(C_filt))
	for (ind,Δϕ_i) in enumerate(ϕ)
		τ_low_ind,τ_high_ind = HelioseismicKernels.time_window_bounce_filter(Δϕ_i,dt,bounce_no)
		f_t[τ_low_ind:τ_high_ind,ind] .= 1
	end

	C_filt .*= f_t

	ax1 = get(kwargs,:ax1) do
		subplot(121)
	end
	ax2 = get(kwargs,:ax2) do
		subplot(122)
	end

	ax1.pcolormesh(rad2deg.(ϕ), t, C./maximum(abs,C),
		cmap="Greys", vmax=1, vmin=-1, rasterized=true)
	ax1.set_xlabel("Angular separation [degrees]", fontsize=12);
	ax1.set_ylabel("Time [hours]", fontsize=12);
	ax1.set_title("Time-distance diagram", fontsize=12);

	ax2.pcolormesh(rad2deg.(ϕ), t, C_filt./maximum(abs,C_filt),
		cmap="Greys", vmax=1, vmin=-1, rasterized=true)
	ax2.set_xlabel("Angular separation [degrees]", fontsize=12);
	ax2.set_ylabel("Time [hours]", fontsize=12);
	ax2.set_title("Filtered bounce", fontsize=12);
end

function plot_C_spectrum(;kwargs...)

	r_src=get(kwargs,:r_src,r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") ν_arr

	Nν_Gfn = length(ν_arr)

	filename = joinpath(SCRATCH_kerneldir[],"Cωℓ_in_range.fits")

	Cωℓ,header = FITS(filename,"r") do f
				read(f[1]),read_header(f[1])
			end

	ν_ind_range = get(kwargs,:ν_ind_range,axes(Cωℓ,1))
	ℓ_range = get(kwargs,:ℓ_range,axes(Cωℓ,2))

	ℓ_edges = collect(first(ℓ_range):last(ℓ_range)+1)
	ν_edges = collect(ν_arr[first(ν_ind_range):last(ν_ind_range)]).*1e3
	spec = parent(Cωℓ)./maximum(abs,Cωℓ)

	observed_freq = readdlm("Schou_freq181qr.1288")
	l = observed_freq[:,1]; ν=observed_freq[:,3]

	ax = get(kwargs,:ax) do
		subplot(111)
	end # default is a single subplot

	ax.pcolormesh(ℓ_edges,ν_edges,spec,
		cmap="Greys",vmax=0.1,vmin=0,rasterized=true);
	ax.set_xlabel(L"$\ell$",fontsize=12);
	ax.set_ylabel("Frequency [mHz]",fontsize=12);
	modes_filter = @. (2e3 < ν < 4.5e3) & (1 < l < 100) & ( mod(l,5)==0)
	ax.plot(l[modes_filter],ν[modes_filter]./1e3,
		label="MDI frequencies",marker="o",ls="None",mfc="0.8",mec="0.3",ms=4)
	ax.legend(loc="lower left")
	ax.set_title("Spectrum of C",fontsize=12)
	tight_layout()

	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"C_spectrum.eps"))
end

function plot_C_spectrum_and_time_distance(;kwargs...)

	ax1 = subplot2grid((1,5),(0,0),colspan=2)
	ax2 = subplot2grid((1,5),(0,2),colspan=2)
	ax3 = subplot2grid((1,5),(0,4),colspan=1)

	plot_C_spectrum(;kwargs...,ax=ax1,save=false)
	plot_time_distance_with_and_without_los(;kwargs...,ax1=ax2,ax2=ax3,save=false)

	plt.gcf().set_size_inches(12,4)
	tight_layout()

	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"C_spectrum_timedistance.eps"))
end

function plot_Ct_groups_of_modes(;t_max=6,t_min=0.3)
	r_src=HelioseismicKernels.Rsun-75e5
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") dt Nt

	t_max_ind = floor(Int,t_max*60^2/dt)+1
	t_min_ind = floor(Int,t_min*60^2/dt)+1

	t_inds = t_min_ind:t_max_ind
	t = t_inds.*dt/60^2
	nt = length(t)

	n1 = Point2D(π/2,0)
	n2 = Point2D(π/2,π/3)

	ℓ_ranges = Iterators.partition(20:99,20)

	C = zeros(nt,length(ℓ_ranges)+1)

	for (ind,ℓ_range) in enumerate(ℓ_ranges)
		C[:,ind] = HelioseismicKernels.Ct(n1,n2,ℓ_range=ℓ_range)[t_inds]
	end

	C[:,end] = HelioseismicKernels.Ct(n1,n2,ℓ_range=20:99)[t_inds]

	fig,ax = subplots(nrows=size(C,2),ncols=1)
	for axis in ax
		axis.set_xlim(first(t),last(t))
	end

	ax[1].set_title("Cross covariance (Δϕ=60 degrees)")

	for (ind,ℓ_range) in enumerate(ℓ_ranges)
		ax[ind].plot(t,C[:,ind],label="ℓ=$(ℓ_range[1]:ℓ_range[end])")
		ax[ind].xaxis.set_major_formatter(ticker.NullFormatter())
		ax[ind].yaxis.set_major_locator(ticker.MaxNLocator(3))
		ax[ind].legend(loc="upper right",bbox_to_anchor=(1.01,1.5))
	end

	ax[end].plot(t,C[:,end])
	ax[end].set_title("Sum over ℓ")
	ax[end].yaxis.set_major_locator(ticker.MaxNLocator(3))

	xlabel("Time [hours]",fontsize=12)
	tight_layout()
	savefig("Ct_lranges.eps")
end

function plot_Ct_stacked_meridionalflow(;t_max=6, t_min=0.3, startind=1, nrows=4, kwargs...)
	filename = joinpath(SCRATCH_kerneldir[], "Ct_n2arr.fits")
	Ct, header = FITS(filename, "r") do f
		read(f[1]), read_header(f[1])
	end

	bounce_no = get(kwargs, :bounce_no, 1)

	r_src = r_src_default
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") dt Nt

	t_max_ind = floor(Int,t_max*60^2/dt)+1
	t_min_ind = floor(Int,t_min*60^2/dt)+1

	t_inds = t_min_ind:t_max_ind
	t = t_inds*dt/60^2
	nt = length(t)

	endind = startind + nrows - 1

	# Points
	n1 = Point2D(header["X1THETA"], header["X1PHI"])
	θmin = header["THMIN"]
	θmax = header["THMAX"]
	Nθ = header["NTHPHI"]
	n2_arr = [Point2D(θ,header["X1PHI"]) for θ in LinRange(θmin, θmax, Nθ)]

	fig, axlist = subplots(nrows=nrows, ncols=1)
	for (ind, n2ind) in enumerate(startind:endind)

		axlist[ind].plot(t, Ct[t_inds,n2ind])

		t_inds_low, t_inds_high = HelioseismicKernels.time_window_bounce_filter(n1, n2_arr[n2ind], dt, bounce_no)

		axlist[ind].axvline(t_inds_low * dt/60^2, color="firebrick", ls="dotted")
		axlist[ind].axvline(t_inds_high * dt/60^2, color="firebrick", ls="dotted")

		t_inds_range = HelioseismicKernels.time_window_indices_by_fitting_bounce_peak(Ct[:,n2ind], n1, n2_arr[n2ind],
			dt=dt, bounce_no=bounce_no)

		axlist[ind].axvline(minimum(t_inds_range) * dt/60^2, color="black", ls="dotted")
		axlist[ind].axvline(maximum(t_inds_range) * dt/60^2, color="black", ls="dotted")
	end

	fig.set_size_inches(6,4)
	fig.tight_layout()
end

function plot_Ct_bounce_fit(;t_max=6, t_min=0.3, kwargs...)
	r_src = get(kwargs, :r_src, r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") ν_start_zeros ν_arr Nt dt dν Nν_Gfn

	t_max_ind = floor(Int,t_max*60^2/dt)+1
	t_min_ind = floor(Int,t_min*60^2/dt)+1

	t_inds = t_min_ind:t_max_ind
	t = t_inds*dt/60^2
	nt = length(t)

	filename = joinpath(SCRATCH_kerneldir[], "Ct.fits")
	C_t, header = FITS(filename, "r") do f
		read(f[1]), read_header(f[1])
	end
	Nt = size(C_t,1)
	C_env = HelioseismicKernels.hilbertenvelope(C_t)

	n1 = Point2D(header["X1THETA"], header["X1PHI"])
	n2 = Point2D(header["X2THETA"], header["X2PHI"])

	t_ind_min, t_ind_max = HelioseismicKernels.time_window_bounce_filter(n1, n2, dt)
	peak_center = argmax(C_env[t_ind_min:t_ind_max]) + t_ind_min - 1

	points_around_max = C_env[peak_center-2:peak_center+2]
	amp = C_env[peak_center]
	A,t0,σt = HelioseismicKernels.gaussian_fit(peak_center-2:peak_center+2, points_around_max)

	t_inds_range = max(1,floor(Int64,t0 - 2σt)):min(ceil(Int64,t0 + 2σt),Nt)

	plot(t, C_t[t_inds])
	plot(t, C_env[t_inds])
	axvline(t_ind_min*dt/60^2, ls="dotted", color="black")
	axvline(t_ind_max*dt/60^2, ls="dotted", color="black")
	axvline(peak_center*dt/60^2, ls="dotted", color="orange")

	axvline(minimum(t_inds_range)*dt/60^2, ls="dotted", color="firebrick")
	axvline(maximum(t_inds_range)*dt/60^2, ls="dotted", color="firebrick")

	# τ_ind_range = HelioseismicKernels.time_window_indices_by_fitting_bounce_peak(
	# 					C_t, x1, x2; dt=dt, kwargs...)
end

function plot_h(x1, x2; kwargs...)

	r_src = get(kwargs, :r_src, r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") ν_start_zeros ν_arr Nt dt dν Nν_Gfn

	filename = joinpath(SCRATCH_kerneldir[], "Ct.fits")
	C_t = FITS(filename, "r") do f
		read(f[1])
	end

	# C_t = HelioseismicKernels.fft_ω_to_t(Cω_x1x2,dν)
	τ_ind_range = HelioseismicKernels.time_window_indices_by_fitting_bounce_peak(
						C_t, x1, x2; dt=dt, kwargs...)
	ht = HelioseismicKernels.ht(Cω_x1x2, x1, x2; τ_ind_range=τ_ind_range, kwargs...)
	hω = HelioseismicKernels.fft_t_to_ω(ht, dt)

	# subplot(411)
	# plot(ν_arr,real(Cω_x1x2[ν_start_zeros .+ (1:Nν_Gfn)]))
	# title("C(x₁,x₂,ω)")

	ax2=subplot(311)
	plot((1:Nt).*dt,C_t,color="black")
	axvline(τ_ind_range[1]*dt,ls="solid")
	axvline(τ_ind_range[end]*dt,ls="solid")
	title("C(x₁,x₂,t)")

	subplot(312,sharex=ax2)
	plot((1:Nt).*dt,ht,color="black")
	xlim(0,60^2*6)

	title("h(x₁,x₂,t)")

	subplot(313)
	plot(ν_arr,imag(hω[ν_start_zeros .+ (1:Nν_Gfn)]),label="imag")
	plot(ν_arr,real(hω[ν_start_zeros .+ (1:Nν_Gfn)]),label="real")
	legend(loc="best")
	title("h(x₁,x₂,ω)")

	tight_layout()

	# gcf().savefig(joinpath(SCRATCH[],"h$los_tag.eps"))
end

function plot_δCω_isotropicδc_with_and_without_los(args...;kwargs...)

	r_src=get(kwargs,:r_src,r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") ν_arr ν_start_zeros Nν_Gfn

	ν_zoom_min = get(kwargs,:ν_zoom_min,2.985e-3)
	ν_zoom_max = get(kwargs,:ν_zoom_max,3.05e-3)
	ν_zoom_inds = @. ν_zoom_min < ν_arr < ν_zoom_max
	ν_arr_zoom = ν_arr[ν_zoom_inds]

	ax_radial = subplot(121)
	ax_los = subplot(122)
	inset_coords = [0.7, 0.1, 0.25, 0.25]
	ax_radial_zoom = ax_radial.inset_axes(inset_coords)
	ax_los_zoom = ax_los.inset_axes(inset_coords)

	filename = joinpath(SCRATCH_kerneldir[],"δCω_isotropicδc_FB.fits")
	if isfile(filename)
		δCω_FB = FITS(filename,"r") do f
					read(f[1])
				end
		δCω_FB = reinterpret(ComplexF64,δCω_FB)
	else
		δCω_FB = HelioseismicKernels.δCω_isotropicδc_firstborn_integrated_over_angle(
			args...;kwargs...)
	end
	δCω_FB = δCω_FB[ν_start_zeros .+ (1:Nν_Gfn)]

	filename = joinpath(SCRATCH_kerneldir[],"δCω_isotropicδc_CmC0.fits")
	if isfile(filename)
		δCω_CmC₀ = FITS(filename,"r") do f
					read(f[1])
				end
		δCω_CmC₀ = reinterpret(ComplexF64,δCω_CmC₀)
	else
		δCω_CmC₀ = HelioseismicKernels.δCω_isotropicδc_C_minus_C0(
			args...;kwargs...)
	end
	δCω_CmC₀ = δCω_CmC₀[ν_start_zeros .+ (1:Nν_Gfn)]

	ax_radial.plot(ν_arr.*1e3,real(δCω_FB),label="First Born",color="grey")
	ax_radial.plot(ν_arr.*1e3,real(δCω_CmC₀),label="C - C₀",
			ls="None",marker=".",color="black")
	ax_radial.set_xlabel("Frequency [mHz]",fontsize=12)
	ax_radial.set_ylabel("ℜ[C]",fontsize=12)
	ax_radial.legend(loc="best")
	ax_radial.set_title("Radial",fontsize=12)

	ax_radial_zoom.plot(ν_arr_zoom.*1e3,real(δCω_FB[ν_zoom_inds]),
		label="First Born",color="grey")
	ax_radial_zoom.plot(ν_arr_zoom.*1e3,real(δCω_CmC₀[ν_zoom_inds]),
		label="C - C₀",ls="None",marker=".",color="black")
	ax_radial_zoom.xaxis.set_major_locator(ticker.MaxNLocator(2))
	ax_radial.indicate_inset_zoom(ax_radial_zoom)

	filename = joinpath(SCRATCH_kerneldir[],"δCω_isotropicδc_FB_los.fits")
	if isfile(filename)
		δCω_FB_los = FITS(filename,"r") do f
						read(f[1])
					end
		δCω_FB_los = reinterpret(ComplexF64,δCω_FB_los)
	else
		δCω_FB_los = HelioseismicKernels.δCω_isotropicδc_firstborn_integrated_over_angle_los(
			args...;kwargs...)
	end
	δCω_FB_los = δCω_FB_los[ν_start_zeros .+ (1:Nν_Gfn)]

	filename = joinpath(SCRATCH_kerneldir[],"δCω_isotropicδc_CmC0_los.fits")
	if isfile(filename)
		δCω_CmC₀_los = FITS(filename,"r") do f
						read(f[1])
					end
		δCω_CmC₀_los = reinterpret(ComplexF64,δCω_CmC₀_los)
	else
		δCω_CmC₀_los = HelioseismicKernels.δCω_isotropicδc_C_minus_C0(
			args...;kwargs...)
	end
	δCω_CmC₀_los = δCω_CmC₀_los[ν_start_zeros .+ (1:Nν_Gfn)]

	ax_los.plot(ν_arr.*1e3,real(δCω_FB_los),label="First Born",color="grey")
	ax_los.plot(ν_arr.*1e3,real(δCω_CmC₀_los),label="C - C₀",
				ls="None",marker=".",color="black")
	ax_los.set_xlabel("Frequency [mHz]",fontsize=12)
	ax_los.set_ylabel("ℜ[C]",fontsize=12)
	ax_los.legend(loc="best")
	ax_los.set_title("Line-of-sight",fontsize=12)

	ax_los_zoom.plot(ν_arr_zoom.*1e3,real(δCω_FB_los[ν_zoom_inds]),
		label="First Born",color="grey")
	ax_los_zoom.plot(ν_arr_zoom.*1e3,real(δCω_CmC₀_los[ν_zoom_inds]),
		label="C - C₀",ls="None",marker=".",color="black")
	ax_los_zoom.xaxis.set_major_locator(ticker.MaxNLocator(2))
	ax_los.indicate_inset_zoom(ax_los_zoom)

	gcf().set_size_inches(8,4)
	tight_layout()
	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"dCw_dc_with_and_without_los.eps"))
end

function plot_Ct_with_and_without_los(;t_min=1.5, t_max=2, kwargs...)
	r_src = r_src_default
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load(joinpath(Gfn_path_src,"parameters.jld2"),Nt,dt)

	t_max_ind = floor(Int,t_max*60^2/dt)+1
	t_min_ind = floor(Int,t_min*60^2/dt)+1

	t_inds = t_min_ind:t_max_ind
	t = t_inds*dt/60^2
	nt = length(t)

	C_rad_200_200 = FITS(joinpath(SCRATCH_kerneldir[], "Ct_r1_200_r2_200.fits"), "r") do f
		read(f[1], t_inds)
	end
	C_rad_200_100 = FITS(joinpath(SCRATCH_kerneldir[], "Ct_r1_200_r2_100.fits"), "r") do f
		read(f[1], t_inds)
	end
	C_los_200_200 = FITS(joinpath(SCRATCH_kerneldir[], "Ct_r1_200_r2_200_los.fits"), "r") do f
		read(f[1], t_inds)
	end
	C_los_200_100 = FITS(joinpath(SCRATCH_kerneldir[], "Ct_r1_200_r2_100_los.fits"), "r") do f
		read(f[1], t_inds)
	end

	fig, ax = subplots(nrows=1, ncols=1)

	scale = max(maximum(abs, C_rad_200_200), maximum(abs, C_rad_200_100))

	plot(t, C_rad_200_200/scale, label=L"r_2=200\,"*"km,\nradial", color="black")
	plot(t, C_rad_200_100/scale, label=L"r_2=100\,"*"km,\nradial", color="0.5", marker=".", markevery=3)
	plot(t, C_los_200_200/scale, label=L"r_2=200\,"*"km,\nprojected", color="0.5", marker="^", markevery=3, ls="dashed")
	plot(t, C_los_200_100/scale, label=L"r_2=100\,"*"km,\nprojected", color="0.2", ls="dashed", marker=".", markevery=3)
	mf = ticker.ScalarFormatter(useMathText=true)
	mf.set_powerlimits((-2,2))
	ax.yaxis.set_major_formatter(mf)

	ax.xaxis.set_major_locator(ticker.MaxNLocator(4))
	ax.yaxis.set_major_locator(ticker.MaxNLocator(3))
	xlim(t[1], t[end])

	xlabel("Time [hours]",fontsize=12)
	ylabel("Cross covariance\n[Arbitrary units]",fontsize=12)
	legend(loc="lower left",fontsize=10, ncol=2)

	fig.set_size_inches(5,3)
	fig.tight_layout()
	fig.savefig(joinpath(SCRATCH[], "C_heights_los.eps"))
end

function plot_δCω_uniform_rotation_FB_rotatedframe(;kwargs...)

	r_src=get(kwargs,:r_src,r_src_default)
	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src)
	@load joinpath(Gfn_path_src,"parameters.jld2") ν_arr ν_start_zeros Nν_Gfn

	ν_zoom_min = 2.985e-3
	ν_zoom_max = 3.05e-3
	ν_zoom_inds = @. ν_zoom_min < ν_arr < ν_zoom_max
	ν_arr_zoom = ν_arr[ν_zoom_inds]

	filename = joinpath(SCRATCH_kerneldir[],"δCω_flows_FB.fits")
	δCω_FB = FITS(filename,"r") do f
		copy(reinterpret(ComplexF64,read(f[1])))
	end

	δCω_FB = δCω_FB[ν_start_zeros .+ (1:Nν_Gfn)]

	filename = joinpath(SCRATCH_kerneldir[],"δCω_flows_rotated.fits")
	δCω_rot = FITS(filename,"r") do f
		copy(reinterpret(ComplexF64,read(f[1])))
	end

	δCω_rot = δCω_rot[ν_start_zeros .+ (1:Nν_Gfn)]

	ax_radial = subplot(111)

	ax_radial.plot(ν_arr.*1e3,imag(δCω_FB),
		color="gray",label="First Born")
	ax_radial.plot(ν_arr.*1e3,imag(δCω_rot),
		color="black",ls="None",marker=".",label="Solid-body rotation")
	ax_radial.set_xlabel("Frequency [mHz]",fontsize=12)
	ax_radial.set_ylabel("ℑ[C]",fontsize=12)
	ax_radial.legend(loc="best")
	ax_radial.set_title("Radial",fontsize=12)

	inset_coords = [0.7, 0.1, 0.25, 0.25]
	ax_radial_zoom = ax_radial.inset_axes(inset_coords)

	ax_radial_zoom.plot(ν_arr_zoom.*1e3,imag(δCω_FB[ν_zoom_inds]),
		color="grey")
	ax_radial_zoom.plot(ν_arr_zoom.*1e3,imag(δCω_rot[ν_zoom_inds]),
		ls="None",marker=".",color="black")
	ax_radial_zoom.xaxis.set_major_locator(ticker.MaxNLocator(2))
	ax_radial.indicate_inset_zoom(ax_radial_zoom)

	tight_layout()
	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"dCw_v_FB_rotated.eps"))
end

# In this function we compute the cross covariance in two ways
# firstly we compute C(x1, x2) and rotate the coordinate frame to obtain C′(x1′, x2′)
# secondly we compute C(x1′, x2′) directly
# these two should be equal
function plot_δCt_rotatedpoints_WignerD(;t_max=2.2, t_min=1, kwargs...)

	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src_default)
	@load joinpath(Gfn_path_src,"parameters.jld2") dt Nt

	t_max_ind = floor(Int,t_max*60^2/dt)+1
	t_min_ind = floor(Int,t_min*60^2/dt)+1

	t_inds = t_min_ind:t_max_ind
	t = t_inds*dt/60^2
	nt = length(t)

	C_direct, hdr = FITS(joinpath(SCRATCH_kerneldir[], "Ct_radial.fits"), "r") do f
		read(f[1]), read_header(f[1])
	end
	x1′ = Point3D(hdr["X1R"], hdr["X1THETA"], hdr["X1PHI"])
	x2′ = Point3D(hdr["X2R"], hdr["X2THETA"], hdr["X2PHI"])

	C_rot, hdr = FITS(joinpath(SCRATCH_kerneldir[], "Ct_rot_radial.fits"), "r") do f
		read(f[1]), read_header(f[1])
	end
	x1 = Point3D(hdr["X1R"], hdr["X1THETA"], hdr["X1PHI"])
	x2 = Point3D(hdr["X2R"], hdr["X2THETA"], hdr["X2PHI"])

	scale = maximum(abs, C_direct[t_inds])

	fig = figure()
	axcov = subplot2grid((3,1), (0,0), rowspan=2)
	axcov.plot(t, C_direct[t_inds] / scale, label="direct", color="grey")
	axcov.plot(t, C_rot[t_inds] / scale, label="rotated", ls="None", marker=".", color="black")
	axcov.legend(loc="best", fontsize=12)
	axcov.set_ylabel("Cross covariance\n[Arbitrary units]", fontsize=12)
	axcov.xaxis.set_major_formatter(ticker.NullFormatter())
	axcov.xaxis.set_major_locator(ticker.MaxNLocator(4))
	axcov.yaxis.set_major_locator(ticker.MaxNLocator(3))
	axcov.yaxis.get_offset_text().set_visible(false)

	axdiff = subplot2grid((3,1), (2,0))
	axdiff.plot(t, (C_direct - C_rot)[t_inds] / scale, color="black")
	axdiff.set_xlabel("Time [hours]", fontsize=12)
	axdiff.set_ylabel("Difference", fontsize=12)
	axdiff.yaxis.set_major_locator(ticker.MaxNLocator(3))
	axdiff.xaxis.set_major_locator(ticker.MaxNLocator(4))
	mf = ticker.ScalarFormatter(useMathText=true)
	mf.set_powerlimits((-2,2))
	axdiff.yaxis.set_major_formatter(mf)

	fig.set_size_inches(6,4)
	fig.tight_layout()

	fig.savefig(joinpath(SCRATCH[], "Ct_flows_rot.eps"))
end

# Kernels

function compute_Kv_10_different_arrival_bounces()
	n1 = Point2D(π/2,0); n2 = Point2D(π/2,π/3)
	kernels = zeros(nr,3)

	for (ind,bounce_no) in enumerate((1,2,4))
		kernels[:,ind] .= imag.(HelioseismicKernels.kernel_uniform_rotation_uplus(n1,n2,ℓ_range=20:100,bounce_no=bounce_no))
	end

	mkpath("$SCRATCH[]/kernels")
	f = FITS("$SCRATCH[]/kernels/K10_different_bounces.fits","w")
	write(f,kernels)
	close(f)

	return kernels
end

function plot_Kv_10_different_arrival_bounces()

	kernel_file = "$SCRATCH[]/kernels/K10_different_bounces.fits"
	if isfile(kernel_file)
		f = FITS(kernel_file)
		kernels = read(f[1])
		close(f)
	else
		kernels = compute_K10_different_arrival_bounces()
	end

	subplot(211)
	plot(r_frac,normalize(kernels[:,1],Inf),label="first bounce")
	xlim(0.8,r[end])
	legend(loc="best")
	title(L"Normalized $K_{10}(r)$",fontsize=12)

	subplot(212)
	plot(r_frac,kernels[:,3],label="fourth bounce")
	xlim(0.8,r[end])
	legend(loc="best")

	xlabel(L"$r/R_\odot$",fontsize=12)

	tight_layout()
end

function plot_Kvrₗ₀(;kwargs...)
	kernel_file = joinpath(SCRATCH_kerneldir[],"Kₗ₀_δτ_u_rθϕ.fits")

	kernel,header = FITS(kernel_file,"r") do f
				read(f[1],:,1,:),read_header(f[1])
	end
	kernel_modes = LM(header["L_MIN"]:header["L_MAX"],0:0)
	s_min = kernel_modes.l_min
	s_max = kernel_modes.l_max

	ls_arr = ["solid","dashed"]
	markers = ["None","."]
	lw = [2,1]
	colors = ["red","blue","green","orange"]

	ax1 = get(kwargs,:ax) do
		subplot(111)
	end

	inset_coords = [0.1, 0.1, 0.5, 0.3]
	ax1_zoom = ax1.inset_axes(inset_coords)

	for s=s_min:s_max
		s_ind = modeindex(kernel_modes,s,0)
		for ax in [ax1,ax1_zoom]
			ax.plot(r_frac,kernel[:,s_ind],label="ℓ=$s",
					ls=ls_arr[div(s-1,4)+1],
					color=colors[mod(div(s-1,2),4)+1],
					marker=markers[div(s-1,4)+1],
					markevery=40,linewidth=lw[div(s-1,4)+1])
		end
	end

	ax1.set_ylabel(L"\mathrm{K}_{v_r,ℓ0}",fontsize=14)
	ax1.legend(loc="best")
	ax1.set_xlim(0.2,r_frac[end])
	ax1.yaxis.set_major_locator(ticker.MaxNLocator(4))

	ax1_zoom.set_xlim(0.97,r_frac[end])
	ax1_zoom.xaxis.set_major_locator(ticker.MaxNLocator(4))

	ax1.indicate_inset_zoom(ax1_zoom)
	gcf().set_size_inches(8,4)
	tight_layout()
	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"Kvr.eps"),dpi=320)
end

function plot_Kvθₗ₀(;kwargs...)

	kernel_file = joinpath(SCRATCH_kerneldir[],"Kₗ₀_δτ_u_rθϕ.fits")

	kernel,header = FITS(kernel_file,"r") do f
				read(f[1],:,2,:),read_header(f[1])
	end
	kernel_modes = LM(header["L_MIN"]:header["L_MAX"],0:0)
	s_min = kernel_modes.l_min
	s_max = kernel_modes.l_max

	ls_arr = ["solid","dashed"]
	markers = ["None","."]
	lw = [2,1]
	colors = ["red","blue","green","orange"]

	ax2 = get(kwargs,:ax) do
		subplot(111)
	end

	inset_coords = [0.1, 0.1, 0.5, 0.3]
	ax2_zoom = ax2.inset_axes(inset_coords)

	for s=s_min:s_max
		s_ind = modeindex(kernel_modes,s,0)
		for ax in [ax2,ax2_zoom]
			ax.plot(r_frac,kernel[:,s_ind],label="ℓ=$s")
		end
	end
	ax2.legend(loc="upper left")
	ax2.set_xlim(0.2,r_frac[end])
	ax2.set_ylabel(L"$\mathrm{K}_{v_θ,ℓ0}$",fontsize=14)
	ax2.set_xlabel(L"$r/R_\odot$",fontsize=12)
	ax2.yaxis.set_major_locator(ticker.MaxNLocator(4))

	ax2_zoom.set_xlim(0.97,r_frac[end])
	ax2_zoom.xaxis.set_major_locator(ticker.MaxNLocator(4))

	ax2.indicate_inset_zoom(ax2_zoom)

	gcf().set_size_inches(8,4)
	tight_layout()
	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"Kvtheta.eps"),dpi=320)
end

function plot_Kvϕₗ₀(;kwargs...)
	kernel_file = joinpath(SCRATCH_kerneldir[],"Kₗ₀_δτ_u_rθϕ.fits")

	kernel,header = FITS(kernel_file,"r") do f
				read(f[1],:,3,:),read_header(f[1])
			end
	kernel_modes = LM(header["L_MIN"]:header["L_MAX"],0:0)
	s_min = kernel_modes.l_min
	s_max = kernel_modes.l_max

	ax1 = get(kwargs,:ax) do
		subplot(111)
	end

	inset_coords = [0.1, 0.1, 0.5, 0.3]
	ax1_zoom = ax1.inset_axes(inset_coords)

	ls_arr = ["solid","dashed"]
	markers = ["None","."]
	lw = [2,1]
	colors = ["red","blue","green","orange"]

	for s=s_min:s_max
		s_ind = modeindex(kernel_modes,s,0)
		for ax in [ax1,ax1_zoom]
			ax.plot(r_frac,kernel[:,s_ind],label="ℓ=$s")
		end
	end
	ax1.set_ylabel(L"$\mathrm{K}_{vϕ,ℓ0}$",fontsize=14)
	ax1.legend(loc="upper left")
	ax1.set_xlim(0.2,r_frac[end])
	ax1.set_xlabel(L"$r/R_\odot$",fontsize=12)
	ax1.axhline(0,color="black",ls="dotted",zorder=1)
	ax1.yaxis.set_major_locator(ticker.MaxNLocator(4))

	ax1_zoom.set_xlim(0.95,r_frac[end])

	ax1.indicate_inset_zoom(ax1_zoom)

	gcf().set_size_inches(8,4)
	tight_layout()
	save = get(kwargs,:save,true)
	save && savefig(joinpath(SCRATCH[],"Kvphi.eps"),dpi=320)
end

function plot_Kvrθₗ₀(;kwargs...)
	ax1 = subplot(211)
	ax2 = subplot(212)

	plot_Kvrₗ₀(ax=ax1,save=false)
	plot_Kvθₗ₀(ax=ax2,save=false)

	gcf().set_size_inches(8,8)
	tight_layout()
	savefig(joinpath(SCRATCH[],"Kvrtheta.eps"),dpi=320)
end

function plot_Kψϕ_s0(;kwargs...)
	kernel_file = joinpath(SCRATCH_kerneldir[],"Kψ_imag.fits")

	kernel = FITS(kernel_file,"r") do f
				read(f[1])
			end

	s_max = size(kernel,2)

	ax = get(kwargs,:ax) do
		subplot(111)
	end

	inset_coords = [0.3, 0.65, 0.5, 0.3]
	ax_zoom = ax.inset_axes(inset_coords)

	ls_arr = ["solid","dashed"]
	markers = ["None","."]
	lw = [2,1]
	colors = ["red","blue","green","orange"]

	for s=1:s_max
		for axis in [ax,ax_zoom]
			axis.plot(r_frac,ρ.*kernel[:,s],label="ℓ=$s",
				ls=ls_arr[div(s-1,4)+1],
				color=colors[mod(div(s-1,2),4)+1],
				marker=markers[div(s-1,4)+1],
				markevery=40,linewidth=lw[div(s-1,4)+1])
		end
	end
	ax.set_ylabel(L"\rho\,K_{{ψ_ϕ},ℓ0}",fontsize=14)
	ax.legend(loc="best")
	ax.set_xlim(0.2,r_frac[end])
	ax.set_xlabel(L"$r/R_\odot$",fontsize=12)
	ax.xaxis.set_major_locator(ticker.MaxNLocator(4))

	ax_zoom.set_xlim(0.97,r_frac[end])
	ax_zoom.xaxis.set_major_locator(ticker.MaxNLocator(4))

	ax.indicate_inset_zoom(ax_zoom)
	tight_layout()

	gcf().set_size_inches(8,4)
	save = get(kwargs,:save,false)
	save && savefig(joinpath(SCRATCH[],"Kpsiphi.eps"),dpi=320)
end

function plot_Kcₗₘ_reim(;kwargs...)
	lmax = get(kwargs,:lmax,50)
	smax = get(kwargs,:smax,2lmax)
	tmax = get(kwargs,:tmax,smax)

	Kst,header = FITS(joinpath(SCRATCH[],"kernels",
		"Kst_δτ_c_lmax$(lmax)_smax$(smax)_tmax$(tmax).fits"),"r") do f
		reinterpret(ComplexF64,read(f[1])),read_header(f[1])
	end

	SHModes = st(0:smax,0:tmax)

	r_frac_cutoff = 0.8
	r_frac_cutoff_zoom = 0.99
	r_frac_trimmed = r_frac[r_frac .> r_frac_cutoff]
	r_frac_zoomed = r_frac[r_frac .> r_frac_cutoff_zoom]
	K_plot = (c.*Kst)[r_frac .> r_frac_cutoff,:]
	K_plot_zoom = (c.*Kst)[r_frac .> r_frac_cutoff_zoom,:]

	modes = [(4,2),(5,3),(7,5)]
	linecolors = ["0.2","0.4","0.6"]
	linestyles = ["solid","dashed","solid"]

	ax_real = get(kwargs,:ax1) do
		subplot(121)
	end
	ax_imag = get(kwargs,:ax2) do
		subplot(122)
	end
	inset_coords = [0.15, 0.1, 0.35, 0.35]
	ax_real_zoom = ax_real.inset_axes(inset_coords)
	ax_imag_zoom = ax_imag.inset_axes(inset_coords)

	ymax = 0
	for ((s,t),lc,ls) in zip(modes,linecolors,linestyles)
		SHind = modeindex(SHModes,(s,t))
		if SHind > size(Kst,2)
			continue
		end
		ax_real.plot(r_frac_trimmed,real(K_plot[:,SHind]),
			label="(l=$s,m=$t)",color=lc,ls=ls)
		ax_real_zoom.plot(r_frac_zoomed,real(K_plot_zoom[:,SHind]),color=lc,ls=ls)
		ax_real_zoom.axhline(0,ls="dotted",color="black",zorder=0)
		ymax = max(ymax,maximum(abs,real(K_plot[:,SHind])))
	end
	ax_real.set_xlabel(L"r/R_\odot",fontsize=12)
	ax_real.set_ylabel(L"\Re [cK]",fontsize=12)
	ax_real.xaxis.set_major_locator(ticker.MaxNLocator(4))
	ax_real.legend(loc="best",fontsize=12)

	ax_real.set_ylim(-ymax*1.1,ymax*1.1)
	ax_real_zoom.set_ylim(ax_real.get_ylim())
	ax_real_zoom.yaxis.set_major_locator(ticker.NullLocator())
	ax_real.indicate_inset_zoom(ax_real_zoom)
	ax_real.set_title("Real part",fontsize=12)

	ymax = 0
	for ((s,t),lc,ls) in zip(modes,linecolors,linestyles)
		SHind = modeindex(SHModes,(s,t))
		if SHind > size(Kst,2)
			continue
		end
		ax_imag.plot(r_frac_trimmed,imag(K_plot[:,SHind]),
			label="(l=$s,m=$t)",color=lc,ls=ls)
		ax_imag_zoom.plot(r_frac_zoomed,imag(K_plot_zoom[:,SHind]),color=lc,ls=ls)
		ax_imag_zoom.axhline(0,ls="dotted",color="black",zorder=0)
		ymax = max(ymax,maximum(abs,imag(K_plot[:,SHind])))
	end
	ax_imag.set_xlabel(L"r/R_\odot",fontsize=12)
	ax_imag.set_ylabel(L"\Im [cK]",fontsize=12)
	ax_imag.xaxis.set_major_locator(ticker.MaxNLocator(4))
	ax_imag.legend(loc="best",fontsize=12)

	ax_imag.set_ylim(-ymax*1.1,ymax*1.1)
	ax_imag_zoom.set_ylim(ax_imag.get_ylim())
	ax_imag_zoom.yaxis.set_major_locator(ticker.NullLocator())
	ax_imag.indicate_inset_zoom(ax_imag_zoom)
	ax_imag.set_title("Imaginary part",fontsize=12)

	tight_layout()
	save = get(kwargs,:save,false)
	save && savefig(joinpath(SCRATCH[],"Kst_somemodes.eps"))
end

function plot_Kcₗₘ_spectrum(;kwargs...)
	lmax = get(kwargs,:lmax,50)
	smax = get(kwargs,:smax,2lmax)
	tmax = get(kwargs,:tmax,smax)

	Kst,header = FITS(joinpath(SCRATCH[],"kernels",
		"Kst_δτ_c_lmax$(lmax)_smax$(smax)_tmax$(tmax).fits"),"r") do f
		reinterpret(ComplexF64,read(f[1])),read_header(f[1])
	end

	SHModes = st(0:smax,0:tmax)
	ax = get(kwargs,:ax) do
		subplot(111)
	end

	r₁ = header["R1"]

	K_spectrum = zeros(0:smax,0:tmax)

	r₁_ind = HelioseismicKernels.radial_grid_index(r₁)

	for (SHind,(s,t)) in enumerate(SHModes)
		K_spectrum[s,t] = abs2(Kst[r₁_ind,SHind])
	end

	p = ax.imshow(permutedims(K_spectrum),cmap="Greys",origin="lower",
		vmax=maximum(K_spectrum)/10,rasterized=true)
	ax.set_xlabel(L"ℓ",fontsize=12)
	ax.set_ylabel(L"m",fontsize=12)
	ax.set_title("Spectrum",fontsize=12)

	ax.yaxis.set_major_locator(ticker.MaxNLocator(5,integer=true))
	ax.xaxis.set_major_locator(ticker.MaxNLocator(5,integer=true))
end

function plot_Kcₗₘ_reim_spectrum(;kwargs...)

	ax1 = subplot(131)
	ax2 = subplot(132)
	ax3 = subplot(133)

	plot_Kcₗₘ_reim(;kwargs...,ax1=ax1,ax2=ax2)
	plot_Kcₗₘ_spectrum(;kwargs...,ax=ax3)

	gcf().set_size_inches(10,4)
	tight_layout()

	savefig(joinpath(SCRATCH[],"Kst_reim_spectrum.eps"),dpi=160)
end

function plot_kernel_timing_scaling_benchmark(n1=Point2D(π/2,0),n2=Point2D(π/2,π/3);
	s_max=10,ℓ_range=20:30:100)

	ns = 3; chunksize = max(1,div(s_max-1,ns))
	s_range = 1:chunksize:s_max
	evaltime = zeros(length(ℓ_range),length(s_range))
	p = Progress(length(evaltime), 1)
	for (s_ind,s_max) in enumerate(s_range), (ℓind,ℓ) in enumerate(ℓ_range)
		evaltime[ℓind,s_ind] = @elapsed HelioseismicKernels.flow_axisymmetric_without_los(
									n1,n2,s_max,ℓ_range=20:ℓ);
		next!(p)
	end

	subplot(211)
	for (s_ind,s_max) in enumerate(s_range)
		plot(ℓ_range,evaltime[:,s_ind],label="s_max = $s_max",marker="o",ms=4,ls="solid")
	end
	xlabel("Maximum ℓ",fontsize=12)
	ylabel("Runtime [sec]",fontsize=12)
	legend(loc="best",fontsize=12)

	subplot(212)
	plot(s_range,evaltime[end,:],marker="o",ms=4,ls="solid",label="ℓ = $(last(ℓ_range))")
	xlabel("Maximum s",fontsize=12)
	ylabel("Runtime [sec]",fontsize=12)

	tight_layout()
end

function plot_2points_same_longitude(n1,n2,ax;kwargs...)
	θ = kwargs[:θ]
	ϕ = get(kwargs,:ϕ,(n1.ϕ + n2.ϕ)/2)
	r_obs_1 = isa(n1,Point2D) ? r_obs_default/Rsun : n1.r/Rsun
	r_obs_2 = isa(n2,Point2D) ? r_obs_default/Rsun : n2.r/Rsun

	x1 = r_obs_1*sin(n1.ϕ)
	z1 = r_obs_1*cos(n1.ϕ)

	x2 = r_obs_2*sin(n2.ϕ)
	z2 = r_obs_2*cos(n2.ϕ)

	ax.plot([x1],[z1],marker="o",ls="None",ms=6,color="black")
	ax.plot([x2],[z2],marker="o",ls="None",ms=6,color="black")

	xmin = all(x->x>0,cos.(θ)) ? -0.1 : -1.1
	xmax = all(x->x<0,cos.(θ)) ?  0.1 :  1.1
	ymin = all(x->x>0,sin.(θ)) ? -0.1 : -1.1
	ymax = all(x->x<0,sin.(θ)) ?  0.1 :  1.1
	ax.set_xlim(xmin,xmax)
	ax.set_ylim(ymin,ymax)

	ax.xaxis.set_major_locator(ticker.MaxNLocator(4))
	ax.yaxis.set_major_locator(ticker.MaxNLocator(4))
end

function plot_K_longitudinal_slice(n1,n2,K_plot::Matrix{Float64};kwargs...)

	θ = get(kwargs,:θ) do
		LinRange(0,π,size(K_plot,2))
	end
	ϕ = get(kwargs,:ϕ) do
		(n1.ϕ + n2.ϕ)/2
	end
	dθ = θ[2] - θ[1]

	vmax = maximum(abs.(K_plot))
	vmax_scale = get(kwargs,:vmax_scale,1)
	vmax /= vmax_scale

	fig = get(kwargs,:fig) do
		figure()
	end

	ax = get(kwargs,:ax) do
		subplot(111,aspect="equal")
	end

	z = r_frac*cos.(θ)'
	x = r_frac*sin.(θ)'

	p = ax.pcolormesh(x,z,K_plot,cmap="RdBu_r",rasterized=true,vmax=vmax,vmin=-vmax)

	if n1.ϕ == ϕ && n2.ϕ == ϕ
		plot_2points_same_longitude(n1,n2,ax;kwargs...,θ=θ)
	end

	ax.set_xlim(-0.1,1.1)
	ax.set_ylim(-1.1,1.1)

	ax.set_xlabel(L"r/R_\odot",fontsize=12)
	ax.set_ylabel(L"r/R_\odot",fontsize=12)

	ax.yaxis.set_major_locator(ticker.MaxNLocator(4))
	ax.xaxis.set_major_locator(ticker.MaxNLocator(4))

	divider = axes_grid1.make_axes_locatable(ax)
	cax = divider.new_horizontal(size="5%",pad=0.05)
	fig.add_axes(cax)
	fig.colorbar(p, cax=cax, extend= (vmax_scale==1) ? "neither" : "both")

	mf = ticker.ScalarFormatter(useMathText=true)
	mf.set_powerlimits((-2,2))
	cax.yaxis.set_major_formatter(mf)
	cax.yaxis.set_offset_position("left")
	cax.yaxis.set_major_locator(ticker.MaxNLocator(3))
end

function plot_Kc_longitudinal_slice(K_fits_filename;kwargs...)
	K,header = FITS(K_fits_filename,"r") do f
		c.*read(f[1]),read_header(f[1])
	end

	x1 = Point3D(header["R1"],header["TH1"],header["PHI1"])
	x2 = Point3D(header["R2"],header["TH2"],header["PHI2"])

	θ_slice = LinRange(header["THSLMN"],header["THSLMX"],header["NTH"])
	ϕ_slice = header["PHSL"]

	plot_K_longitudinal_slice(x1,x2,K;θ=θ_slice,ϕ=ϕ_slice,kwargs...)
end

function plot_Kc_3D_longitudinal_slice_twopanels(;kwargs...)
	fig = get(kwargs,:fig) do
		figure()
	end
	ax1 = get(kwargs,:ax_long_1) do
		subplot(121,aspect="equal")
	end
	ax2 = get(kwargs,:ax_long_2) do
		subplot(122,aspect="equal")
	end

	filename = joinpath(SCRATCH_kerneldir[],"Krθ_from_Kst.fits")
	plot_Kc_longitudinal_slice(filename;kwargs...,ax=ax1,fig = fig)

	header = FITS(filename,"r") do f
		read_header(f[1])
	end

	lmax = header["LMAX"]

	filename = joinpath(SCRATCH_kerneldir[],
			"Kc_3D_long_slice_jcutoff$(lmax).fits")

	plot_Kc_longitudinal_slice(filename;kwargs...,ax=ax2,fig = fig)

	ax1.set_title("This work",fontsize=12)
	ax2.set_title("3D kernel",fontsize=12)

	save = get(kwargs,:save,true)
	if save
		fig.set_size_inches(8,4)
		tight_layout()
		fig.savefig(joinpath(SCRATCH[],"compare_3D_kernel_longslice.eps"),
		dpi=get(kwargs,:dpi,160))
	end
end

function plot_Ku_longitudinal_slice(;kwargs...)

	K_fits_filename = joinpath(SCRATCH_kerneldir[],"Ku_rθslice_from_Kₗₘ.fits")

	K,header = FITS(K_fits_filename,"r") do f
		c.*read(f[1]),read_header(f[1])
	end

	x1 = Point3D(header["R1"],header["TH1"],header["PHI1"])
	x2 = Point3D(header["R2"],header["TH2"],header["PHI2"])

	θ_slice = LinRange(header["THSLMN"],header["THSLMX"],header["NTH"])
	ϕ = header["PHSL"]

	fig = figure()

	ax1 = subplot(131,aspect="equal")
	vr_scale = get(kwargs,:vr_scale,1)
	plot_K_longitudinal_slice(x1,x2,K[:,:,1];θ=θ_slice,ϕ=ϕ,kwargs...,
		fig=fig,ax=ax1,vmax_scale=vr_scale)

	ax2 = subplot(132,aspect="equal")
	vθ_scale = get(kwargs,:vθ_scale,1)
	plot_K_longitudinal_slice(x1,x2,K[:,:,2];θ=θ_slice,ϕ=ϕ,kwargs...,
		fig=fig,ax=ax2,vmax_scale=vθ_scale)

	ax3 = subplot(133,aspect="equal")
	vϕ_scale = get(kwargs,:vϕ_scale,1)
	plot_K_longitudinal_slice(x1,x2,K[:,:,3];θ=θ_slice,ϕ=ϕ,kwargs...,
		fig=fig,ax=ax3,vmax_scale=vϕ_scale)

	fig.set_size_inches(12,4)
	fig.tight_layout()
end

function plot_2points_same_latitude(n1,n2,ax;kwargs...)
	θ = get(kwargs,:θ,(n1.θ + n2.θ)/2)
	ϕ = kwargs[:ϕ]
	r_obs_1 = isa(n1,Point2D) ? r_obs_default/Rsun : n1.r/Rsun
	r_obs_2 = isa(n2,Point2D) ? r_obs_default/Rsun : n2.r/Rsun

	x1 = r_obs_1*cos(n1.ϕ)
	y1 = r_obs_1*sin(n1.ϕ)

	x2 = r_obs_2*cos(n2.ϕ)
	y2 = r_obs_2*sin(n2.ϕ)

	ax.plot([x1],[y1],marker="o",ls="None",ms=6,mfc="yellow",mec="black")
	ax.plot([x2],[y2],marker="o",ls="None",ms=6,mfc="yellow",mec="black")

	xmin = all(x->x>0,cos.(ϕ)) ? -0.1 : minimum(cos.(ϕ))*1.1
	xmax = all(x->x<0,cos.(ϕ)) ?  0.1 : maximum(cos.(ϕ))*1.1
	ymin = all(x->x>0,sin.(ϕ)) ? -0.1 : minimum(sin.(ϕ))*1.1
	ymax = all(x->x<0,sin.(ϕ)) ?  0.1 : maximum(sin.(ϕ))*1.1

	ax.xaxis.set_major_locator(ticker.MaxNLocator(4))
	ax.yaxis.set_major_locator(ticker.MaxNLocator(4))
end

function plot_K_latitudinal_slice(n1,n2,K_plot::Matrix{Float64};kwargs...)

	ϕ = get(kwargs,:ϕ) do
		LinRange(0,2π,4size(K_plot,2))
	end
	θ = get(kwargs,:θ) do
		(n1.θ + n2.θ)/2
	end

	vmax = maximum(abs.(K_plot))
	vmax_scale = get(kwargs,:vmax_scale,1)
	vmax /= vmax_scale

	fig = get(kwargs,:fig) do
		figure()
	end

	ax = get(kwargs,:ax) do
		subplot(111,aspect="equal")
	end

	x = r_frac*cos.(ϕ)';
	y = r_frac*sin.(ϕ)';

	p = ax.pcolormesh(x,y,K_plot,
		cmap="RdBu_r",rasterized=true,vmax=vmax,vmin=-vmax)

	if n1.θ == θ && n2.θ == θ
		plot_2points_same_latitude(n1,n2,ax;kwargs...,ϕ=ϕ)
	end

	ax.set_xlabel(L"r/R_\odot",fontsize=12)
	ax.set_ylabel(L"r/R_\odot",fontsize=12)

	ax.set_xlim(-1.1,1.1)
	ax.set_ylim(-1.1,1.1)

	ax.yaxis.set_major_locator(ticker.MaxNLocator(4))
	ax.xaxis.set_major_locator(ticker.MaxNLocator(4))

	divider = axes_grid1.make_axes_locatable(ax)
	cax = divider.new_horizontal(size="5%",pad=0.05)
	fig.add_axes(cax)
	fig.colorbar(p, cax=cax, extend= (vmax_scale==1) ? "neither" : "both")

	mf = ticker.ScalarFormatter(useMathText=true)
	mf.set_powerlimits((-2,2))
	cax.yaxis.set_major_formatter(mf)
	cax.yaxis.set_offset_position("left")
	cax.yaxis.set_major_locator(ticker.MaxNLocator(3))
end

function plot_Kc_latitudinal_slice(K_fits_filename;kwargs...)
	K,header = FITS(K_fits_filename,"r") do f
		c.*read(f[1]),read_header(f[1])
	end

	x1 = Point3D(header["R1"],header["TH1"],header["PHI1"])
	x2 = Point3D(header["R2"],header["TH2"],header["PHI2"])

	ϕ = LinRange(header["PHSLMN"],header["PHSLMX"],header["NPHI"])

	plot_K_latitudinal_slice(x1,x2,K;ϕ=ϕ,kwargs...)
end

function plot_Kc_3D_latitudinal_slice_twopanels(;kwargs...)
	fig = get(kwargs,:fig) do
		figure()
	end
	ax1 = get(kwargs,:ax_lat_1) do
		subplot(121,aspect="equal")
	end
	ax2 = get(kwargs,:ax_lat_2) do
		subplot(122,aspect="equal")
	end

	filename = joinpath(SCRATCH_kerneldir[],"Krϕ_from_Kst.fits")
	plot_Kc_latitudinal_slice(filename;kwargs...,ax=ax1,fig=fig)

	header = FITS(filename,"r") do f
		read_header(f[1])
	end

	lmax = header["LMAX"]

	filename = joinpath(SCRATCH_kerneldir[],
		"Kc_3D_lat_slice_jcutoff$(lmax).fits")

	plot_Kc_latitudinal_slice(filename;kwargs...,ax=ax2,fig=fig)

	ax1.set_title("This work",fontsize=12)
	ax2.set_title("3D kernel",fontsize=12)

	save = get(kwargs,:save,true)
	if save
		fig.set_size_inches(8,4)
		tight_layout()
		fig.savefig(joinpath(SCRATCH[],"compare_3D_kernel_latslice.eps"),
		dpi=get(kwargs,:dpi,160))
	end
end

function plot_Ku_latitudinal_slice(;kwargs...)

	K_fits_filename = joinpath(SCRATCH_kerneldir[],"Ku_rϕslice_from_Kₗₘ.fits")

	K,header = FITS(K_fits_filename,"r") do f
		c.*read(f[1]),read_header(f[1])
	end

	x1 = Point3D(header["R1"],header["TH1"],header["PHI1"])
	x2 = Point3D(header["R2"],header["TH2"],header["PHI2"])

	θ = header["THSL"]
	ϕ_slice = LinRange(header["PHSLMN"],header["PHSLMX"],header["NPHI"])

	fig = figure()

	ax1 = subplot(131,aspect="equal")
	vr_scale = get(kwargs,:vr_scale,1)
	plot_K_latitudinal_slice(x1,x2,K[:,:,1];θ=θ,ϕ=ϕ_slice,fig=fig,ax=ax1,
		kwargs...,vmax_scale=vr_scale)

	ax2 = subplot(132,aspect="equal")
	vθ_scale = get(kwargs,:vθ_scale,1)
	plot_K_latitudinal_slice(x1,x2,K[:,:,2];θ=θ,ϕ=ϕ_slice,fig=fig,ax=ax2,
		kwargs...,vmax_scale=vθ_scale)

	ax3 = subplot(133,aspect="equal")
	vϕ_scale = get(kwargs,:vϕ_scale,1)
	plot_K_latitudinal_slice(x1,x2,K[:,:,3];θ=θ,ϕ=ϕ_slice,fig=fig,ax=ax3,
		kwargs...,vmax_scale=vϕ_scale)

	fig.set_size_inches(12,4)
	fig.tight_layout()
end

function plot_Kc_3D_lat_long_slice_twopanels(;kwargs...)
	fig = figure()
	ax_lat_1 = subplot(221,aspect="equal")
	ax_lat_2 = subplot(222,aspect="equal")
	ax_long_1 = subplot(223,aspect="equal")
	ax_long_2 = subplot(224,aspect="equal")

	vmax_scale_lat = get(kwargs,:vmax_scale_lat) do
		200
	end
	vmax_scale_long = get(kwargs,:vmax_scale_long) do
		20
	end

	plot_Kc_3D_latitudinal_slice_twopanels(fig=fig,ax_lat_1=ax_lat_1,
		ax_lat_2=ax_lat_2,vmax_scale=vmax_scale_lat,save=false)
	plot_Kc_3D_longitudinal_slice_twopanels(fig=fig,ax_long_1=ax_long_1,
		ax_long_2=ax_long_2,vmax_scale=vmax_scale_long,save=false)

	fig.set_size_inches(8,8)
	tight_layout()

	fig.savefig(joinpath(SCRATCH[],"compare_3D_kernel_latlongslice.eps"),
		dpi=get(kwargs,:dpi,160))
end

function plot_Kv_10_different_n1n2()
	filename = joinpath(SCRATCH_kerneldir[],"K_δτ_v10_equator.fits")
	K1 = FITS(filename,"r") do f
		read(f[1])
	end

	filename = joinpath(SCRATCH_kerneldir[],"K_δτ_v10_diagonal.fits")
	K2 = FITS(filename,"r") do f
		read(f[1])
	end

	filename = joinpath(SCRATCH_kerneldir[],"K_δτ_v10_meridian.fits")
	K3 = FITS(filename,"r") do f
		read(f[1])
	end

	scale = max(maximum(abs,K1),maximum(abs,K2),maximum(abs,K3))

	ax1 = subplot(131)
	ax2 = subplot(132)
	ax3 = subplot(133)

	ax1.plot(r_frac,K1./scale,color="black",lw=2)
	ax2.plot(r_frac,K2./scale,color="black",lw=2)
	ax3.plot(r_frac,K3./scale,color="black",lw=2)

	ax2.set_ylim(ax1.get_ylim())
	ax3.set_ylim(ax1.get_ylim())

	for ax in [ax1,ax2,ax3]
		ax.yaxis.set_major_locator(ticker.MaxNLocator(3))
		ax.xaxis.set_major_locator(ticker.MaxNLocator(3))
		ax.set_xlabel(L"$R/R_\odot$",fontsize=12)
	end

	ax1.set_ylabel("Kvϕ_10",fontsize=12)

	gcf().set_size_inches(8,2)
	tight_layout()
	savefig(joinpath(SCRATCH[],"kernel_pairs_of_points.eps"),dpi=320)
end

function plot_kernel_with_and_without_los_sameheight(args...;kwargs...)

	K2D = HelioseismicKernels.flow_axisymmetric_without_los(args...;kwargs...)
	Base.GC.gc()
	K2D_los = HelioseismicKernels.flow_axisymmetric_with_los(args...;kwargs...)

	plot(r_frac,K2D[:,-1,1],label="without LoS")
	plot(r_frac,K2D_los[:,-1,1],label="with LoS")
	xlim(0.7,r_frac[end])
	legend(loc="best",fontsize=12)
	xlabel(L"$r/R_\odot$",fontsize=12)
	ylabel("Kernel [s/(cm/s Mm)]",fontsize=12)
	title("ϕ component of kernel")
	tight_layout()
end

function plot_Ku_rotatedpoints(;kwargs...)
	f1 = joinpath(SCRATCH_kerneldir[], "Kst_δτ_u_jmax100_lmax1_mmax0_los_twopoints.fits")
	f2 = joinpath(SCRATCH_kerneldir[], "Kst_δτ_u_jmax100_lmax1_mmax0_los.fits")
	header = FITSIO.read_header(f1)
	lmin = header["L_MIN"]::Int
	lmax = header["L_MAX"]::Int
	mmin = header["M_MIN"]::Int
	mmax = header["M_MAX"]::Int
	SHModes = LM(lmin:lmax, mmin:mmax)
	K1 = reinterpret(ComplexF64, FITSIO.fitsread(f1, 1, :, :, modeindex(SHModes, (1,0)), 2))
	K2 = reinterpret(ComplexF64, FITSIO.fitsread(f2, 1, :, :, modeindex(SHModes, (1,0))))
	f, ax = subplots(nrows = 1, ncols = 1)
	ax.plot(r_frac, -2imag(K1[:,3]), color = "grey")
	ax.plot(r_frac, -2imag(K2[:,3]), marker = ".", color = "black", ls = "None", markevery = 20)
	xlim(0.5, 1.01)
	xlabel(L"r/R_\odot", fontsize = 11)
	ylabel(L"\Im\left[K^{1}_{10}\left(r\right)\right] " * L"\left[s^2/cm^4\right]", fontsize = 11)
	mf = ticker.ScalarFormatter(useMathText=true)
	mf.set_powerlimits((-2,2))
	ax.xaxis.set_major_locator(ticker.MaxNLocator(4))
	ax.yaxis.set_major_locator(ticker.MaxNLocator(3))
	ax.yaxis.set_major_formatter(mf)
	ax2 = ax.inset_axes([0.15, 0.18, 0.6, 0.5])
	ax2.plot(r_frac, -2imag(K1[:,3]), color = "grey")
	ax2.plot(r_frac, -2imag(K2[:,3]), marker = ".", color = "black", ls = "None", markevery = 20)
	ax2.set_xlim(0.96, 1.003)
	ax2.xaxis.set_major_locator(ticker.MaxNLocator(3))
	ax2.yaxis.set_major_locator(ticker.MaxNLocator(3))
	ax2.yaxis.set_major_formatter(mf)
	ax.indicate_inset_zoom(ax2)
	f.set_size_inches(4.5,3.2)
	f.tight_layout()
	f.savefig("plots/K_n1n2_los_direct_rotated.eps")
end

function plot_Kl0_rθ()
	filename = "Kl0_δτ_u_rθϕ.fits"
	filepath = joinpath(SCRATCH_kerneldir[],filename)

	header = FITS(filepath,"r") do f
		read_header(f[1])
	end;
	SHModes = LM(header["L_MIN"]:header["L_MAX"],
		header["M_MIN"]:header["M_MAX"]);

	r_frac_top_cutoff = 1+100e5/Rsun
	r_top_cutoffind = searchsortedfirst(r_frac,r_frac_top_cutoff)
	r_frac_bot_cutoff = 0.97
	r_bot_cutoffind = searchsortedfirst(r_frac,r_frac_bot_cutoff)
	r_inds_range = r_bot_cutoffind:r_top_cutoffind

	K = FITS(filepath,"r") do f
		arr = read(f[1],:,1:2,:)
		SHArray(arr,(axes(arr)[1:2]...,SHModes))
	end

	# Meridional flows
	# We choose even ℓ
	modes = [(2,0), (4,0), (6,0)]
	fig, axrθ = subplots(nrows=length(modes), ncols=2)
	for (ind,mode) in enumerate(modes)
		l, = mode
		axrθ[ind,1].plot(r_frac[r_inds_range], K[r_inds_range,1,mode],
			label="ℓ=$l",color="black")

		axrθ[ind,2].plot(r_frac[r_inds_range], K[r_inds_range,2,mode],
			label="ℓ=$l",color="black")
	end

	for ax in axrθ
		ax.margins(0.05,0.1)
		ax.legend(loc="best",fontsize=12)

		ax.xaxis.set_major_locator(ticker.MaxNLocator(3))
	end
	for ax in (@view axrθ[end,:])
		ax.set_xlabel(L"r/R_\odot",fontsize=12)
	end
	for ax in (@view axrθ[:,1])
		ax.set_ylabel(L"\Re[K]",fontsize=12)
	end
	axrθ[1,1].set_title(L"K_r",fontsize=12)
	axrθ[1,2].set_title(L"K_\theta",fontsize=12)

	fig.set_size_inches(6,4)
	tight_layout()

	for ax in axrθ
		t = text_from_offsettext(ax.yaxis)
		ax.yaxis.get_offset_text().set_visible(false)
		ax.text(0, 1.02, t, transform=ax.transAxes, fontsize=10)
	end

	savefig(joinpath(SCRATCH[],"Kl0rtheta.eps"))
end

function plot_Kl0_ϕ()
	filename = "Kl0_δτ_u_rθϕ.fits"
	filepath = joinpath(SCRATCH_kerneldir[],filename)

	header = FITS(filepath,"r") do f
		read_header(f[1])
	end
	SHModes = LM(header["L_MIN"]:header["L_MAX"],
		header["M_MIN"]:header["M_MAX"])

	K = FITS(filepath,"r") do f
		arr = read(f[1],:,3,:)
		SHArray(arr,(axes(arr,1),SHModes))
	end

	r_frac_top_cutoff = 1+100e5/Rsun
	r_top_cutoffind = searchsortedfirst(r_frac,r_frac_top_cutoff)
	r_frac_bot_cutoff = 0.97
	r_bot_cutoffind = searchsortedfirst(r_frac,r_frac_bot_cutoff)
	r_inds_range = r_bot_cutoffind:r_top_cutoffind

	# Rotation
	# We choose odd ℓ
	fig,axϕ = subplots(nrows=1,ncols=3)
	modes = [(1,0), (3,0), (5,0)]
	for (ind,mode) in enumerate(modes)
		l, = mode
		axϕ[ind].plot(r_frac[r_inds_range], K[r_inds_range,mode],
			label="ℓ = $l",color="black")
	end
	# mf = ticker.ScalarFormatter(useMathText=true)
	# mf.set_powerlimits((-2,2))
	for ax in axϕ
		ax.margins(0.05,0.1)
		ax.legend(loc="best",fontsize=12)
		ax.set_xlabel(L"r/R_\odot",fontsize=12)
		ax.yaxis.set_major_locator(ticker.MaxNLocator(3))
		# ax.yaxis.set_major_formatter(mf)
	end
	axϕ[1].set_ylabel(L"\Im[K_\phi]",fontsize=12)

	fig.set_size_inches(7,2)
	tight_layout()

	for ax in axϕ
		t = text_from_offsettext(ax.yaxis)
		ax.yaxis.get_offset_text().set_visible(false)
		ax.text(0, 1.02, t, transform=ax.transAxes, fontsize=10)
	end

	savefig(joinpath(SCRATCH[],"Kl0phi.eps"))
end

# Runtimes

function plotruntimes_soundspeed()
	t = readdlm("runtimes",comments=true)
	sizes = t[:,1]
	runtimes_k3D = t[:,2]/60^2
	runtimes_k1pt5d_rad = t[:,3]/60^2
	runtimes_k1pt5d_los = t[:,4]/60^2
	ax = subplot(121)
	ax.plot(sizes,runtimes_k3D,color="grey",marker=".",
		label=L"\mathrm{K}\left(r,\theta,\phi \right)")
	ax.plot(sizes,runtimes_k1pt5d_rad,color="black",
		label=L"\mathrm{K}_{\ell m}\left(r\right)"*", radial",marker="s",ms=4)
	ax.plot(sizes,runtimes_k1pt5d_los,color="black",
		ls="dashed",label=L"\mathrm{K}_{\ell m}\left(r\right)"*
		", line-of-sight",marker="^",ms=4)
	legend(loc="best",fontsize=12)
	xlabel("Maximum ℓ",fontsize=12)
	ylabel("Net computing time [hours]",fontsize=12)
	ax.yaxis.set_major_locator(ticker.MaxNLocator(5))

	axtwin = twinx()
	scale = 60/28 # minutes per core
	axtwin.plot(sizes,runtimes_k3D .* scale,color="grey",marker=".")
	axtwin.plot(sizes,runtimes_k1pt5d_rad .* scale,color="black",
		marker="s",ms=4)
	axtwin.plot(sizes,runtimes_k1pt5d_los .* scale,color="black",
		ls="dashed",marker="^",ms=4)
	axtwin.yaxis.set_major_locator(ticker.MaxNLocator(5))
	ylabel("Run-time on 28 cores [minutes]",fontsize=12,labelpad=10)
	axtwin.yaxis.set_major_locator(ticker.MaxNLocator(5))

	ax = subplot(122)
	np = [28,56,78,112]
	runtime_rad_ℓ20 = t[end,3:2:end]
	runtime_los_ℓ20 = t[end,4:2:end]

	plot(np,runtime_rad_ℓ20./np./60,color="black",ls="solid",
		marker="o",label="radial",ms=4)
	plot(np,runtime_los_ℓ20./np./60,color="black",ls="dashed",
		marker="o",label="line-of-sight",ms=4)
	ylabel("Run-time [minutes]",fontsize=12)
	xlabel("Number of cores",fontsize=12)
	legend(loc="best",fontsize=12)
	ylim(0,runtime_los_ℓ20[1]/28/60*1.05)
	xticks(np)

	gcf().set_size_inches(9,4)
	tight_layout()
	savefig(joinpath(SCRATCH[],"runtimes.eps"))
end

function plotruntimes_flows()
	t = readdlm("runtimesflows56p",comments=true)
	lmax = t[:,1]
	trad_hr = t[:,2]/60^2
	trad56p_hr = t[:,3]/60^2
	tlos_hr = t[:,4]/60^2
	tlos56p_hr = t[:,5]/60^2

	fig, axlist = subplots(nrows=2, ncols=1)
	ax = axlist[1]
	ax.plot(lmax, trad_hr, color="0.5",marker=".",
		label="radial")
	ax.plot(lmax, tlos_hr, color="black",
		label="line of sight",marker="s",ms=4)

	ax.set_yscale("log", basey=10)
	ax.legend(loc="best",fontsize=12)
	# ax.yaxis.set_major_locator(ticker.MaxNLocator(5))
	ax.xaxis.set_major_locator(ticker.MaxNLocator(4,integer=true))
	# ax.yaxis.set_major_formatter(ticker.StrMethodFormatter("{x:.0f}"))
	# ax.yaxis.set_minor_formatter(ticker.NullFormatter())
	tk = [1,10,100]
	ax.set_yticks(tk)
	ax.set_yticklabels(string.(tk))

	# Linear fit
	p = polyfit(lmax[end-2:end], log.(trad_hr[end-2:end]), 1)
	lmax_dense = LinRange(lmax[end-2], lmax[end], 10)
	trend = (lmax_dense/lmax_dense[1]).^2 * trad_hr[end-2] * 3
	ax.plot(lmax_dense, trend, color="black", ls="dashed")
	ax.text(lmax_dense[1], trend[1]*2, L"\ell_{max}^2", fontsize=12)

	ax2 = axlist[2]
	t = readdlm("runtimes",comments=true)
	lmax_ss = t[:,1]
	runtimes_ss_rad = t[:,3]/60^2
	plot(lmax_ss, runtimes_ss_rad, color="black", ls="dotted", label="sound speed")
	plot(lmax, trad_hr, color="black", label="flows")
	ax2.legend(loc="best", fontsize=11)
	ax2.set_xlabel("Maximum angular degree "*L"\ell_{max}",fontsize=12)
	ax2.xaxis.set_major_locator(ticker.MaxNLocator(4,integer=true))

	# Trick to add a common y label to both subplots
	fig.add_subplot(111, frameon=false)
	plt.tick_params(labelcolor="none", top=false, bottom=false, left=false, right=false)
	ylabel("Net computing time [hours]", fontsize=12)

	gcf().set_size_inches(5,4)
	fig.tight_layout()
	savefig(joinpath(SCRATCH[],"runtimesflows.eps"))
end

function plotruntime_flows_2()
	radout = readdlm("kernel_radial224.out")
	Tradlmax = Int.(radout[2:end,2])
	Tradtime = Int.(radout[2:end,3]) * 224 / 60^2

	losout = readdlm("kernel_los224.out")
	Tlosmax = Int.(losout[2:end,2])
	Tlostime = Int.(losout[2:end,3]) * 224 / 60^2

	Δt_los_rad = Tlostime - Tradtime
	perc_Δt_los_rad = Δt_los_rad ./ Tradtime
	percentage_to_abs = LinearInterpolation(reverse(perc_Δt_los_rad), reverse(Δt_los_rad), extrapolation_bc = Line())

	f, axlist = subplots(nrows = 2, ncols = 1)
	axlist[1].semilogy(Tradlmax, Tradtime, color= "black", label = "radial")
	axlist[1].semilogy(Tlosmax, Tlostime, color = "grey", marker = ".", ls = "dashed", label = "line-of-sight")
	axlist[1].legend(loc="best", fontsize = 11)
	axlist[1].xaxis.set_major_formatter(ticker.NullFormatter())
	axlist[1].yaxis.set_major_locator(ticker.LogLocator())
	axlist[1].yaxis.set_major_formatter(ticker.ScalarFormatter())
	axlist[1].set_ylabel("Computation\n time [hr]", fontsize = 11)

	axlist[2].plot(Tradlmax, perc_Δt_los_rad, color = "black", marker = ".")
	axlist[2].set_ylabel("Difference", fontsize = 11)
	axlist[2].set_xlabel("Maximum angular degree", fontsize = 11)
	axlist[2].yaxis.set_major_formatter(ticker.PercentFormatter(xmax = 1))
	ax2_2 = axlist[2].twinx()
	ax2_2.plot(Tradlmax, perc_Δt_los_rad, color = "black", marker = ".")
	ax2_2.yaxis.set_major_locator(ticker.LinearLocator(4))
	axlist[2].yaxis.set_major_locator(ticker.LinearLocator(4))
	yt = axlist[2].get_yticks()
	ax2_2.set_yticklabels(round.(Int, percentage_to_abs.(yt)))
	ax2_2.set_ylabel("Hours", fontsize = 11)

	f.subplots_adjust(hspace = 0)
	f.set_size_inches(4,4)
	f.tight_layout()
	f.savefig("plots/runtimes.eps")
end

function compare_radial_los_kernel(; s_max = 15)
	Kradfile = joinpath(SCRATCH_kerneldir[], "Kst_δτ_u_jmax100_lmax$(s_max)_mmax$(s_max).fits")
	_Krad = reinterpret(ComplexF64, FITSIO.fitsread(Kradfile, 1, :, 2:2, :)::Array{Float64, 3})
	Krad = SHArray(_Krad, (:, :, LM(0:s_max, 0:s_max)))
	Klosfile = joinpath(SCRATCH_kerneldir[], "Kst_δτ_u_jmax100_lmax$(s_max)_mmax$(s_max)_los.fits")
	_Klos = reinterpret(ComplexF64, FITSIO.fitsread(Klosfile, 1, :, 2:2, :)::Array{Float64, 3})
	Klos = SHArray(_Klos, (:, :, LM(0:s_max, 0:s_max)))

	nrows = 4; ncols = 4
	f, axlist = subplots(nrows = nrows, ncols = ncols)
	r_in_cutoff = 0.99
	r_out_cutoff = 1
	r_inds = searchsortedfirst(r_frac, r_in_cutoff):searchsortedlast(r_frac, r_out_cutoff)

	for ax in axlist
		ax.set_visible(false)
	end

	for (row,l) in enumerate(0:8:s_max), (col,m) in enumerate(0:8:l)
		axlist[row,col].plot(r_frac[r_inds], real(Krad[r_inds, 1, (l,m)]),
			color = "grey", label = "radial")
		axlist[row,col].plot(r_frac[r_inds], real(Klos[r_inds, 1, (l,m)]),
			ls = "dashed", marker = ".", markevery = 20,
			ms = 6, color = "black", label = "line-of-sight")
		axlist[row,col].xaxis.set_major_locator(ticker.MaxNLocator(2))
		axlist[row,col].set_visible(true)
		axlist[row,col].yaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=true))
		axlist[row,col].set_title("($l, $m)", loc = "right", fontsize = 11)
		if row != nrows
			axlist[row,col].xaxis.set_major_formatter(ticker.NullFormatter())
		else
			axlist[row,col].set_xlabel(L"r/R_\odot", fontsize = 11)
		end
	end

	lg = axlist[1,1].get_legend_handles_labels()
	axlist[1,ncols].legend(lg..., loc="best", fontsize = 11)
	axlist[1,ncols].set_visible(true)
	axlist[1,ncols].axis("off")

	axl = f.add_subplot(111, frame_on=false)
	tick_params(labelcolor="none", bottom=false, left = false)
	axl.set_ylabel("Sensitivity kernel "* L"\left[s^2/cm^4\right]", fontsize = 11)

	f.set_size_inches(8,5)
	f.tight_layout()
	f.savefig("plots/kernel_radial_los_limb.eps")
end

# Green function

# function plot_Gfn_with_krishnendu(ℓ=40,ν=3e-3;δc=0,kwargs...)
# 	# Gfn_krishnendu_path = "/scratch/jb6888/Greenfn_krishnendu"

# 	r_src = get(kwargs,:r_src,r_src_default)

# 	c_scale=1+δc

# 	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(r_src,c_scale=c_scale)
# 	println("Loading Greenfn from $Gfn_path_src")
# 	@load joinpath(Gfn_path_src,"parameters.jld2") ν_arr

# 	ν_ind = argmin(abs.(ν_arr .- ν))

# 	ν_ongrid = ν_arr[ν_ind]
# 	ω = 2π*ν_ongrid
# 	γ = HelioseismicKernels.γ_damping
# 	ω_γ = ω-im*γ

# 	ρ = HelioseismicKernels.ρ
# 	c = HelioseismicKernels.c
# 	g = HelioseismicKernels.g
# 	N² = HelioseismicKernels.N2

# 	@. N²+= g^2/c^2*(1-1/c_scale^2)
# 	@. c*=c_scale

# 	a = @. g/c^2
# 	Sₗ² = @. ℓ*(ℓ+1)*c^2/r^2
#     b = @. (Sₗ² /ω_γ^2 -1)/ (ρ*c^2)
#     c = @. ρ*(ω_γ^2 - N²)

#     matrix_kri = ho_endpt_so.genmatrix(a,b,c,r)
#     src = HelioseismicKernels.source(ω,ℓ,r_src=r_src)[1]
#     src_kri = vcat(src[nr:end],[0,0],src[1:nr-1])
#     sol_kri = matrix_kri\src_kri
#     ξ_kri = sol_kri[1:nr]
# 	p_kri = sol_kri[nr+1:2nr]

# 	G1_kri = @. HelioseismicKernels.Ω(ℓ,0)/(ρ*r*ω^2) * p_kri

# 	# Gfn_krishnendu = npzread(joinpath(Gfn_krishnendu_path,
# 	# 	"Gr_nuind_$(ν_ind-1)_ell$(ℓ)_dc$(@sprintf "%g" δc).npy"))[1,:]
# 	# r_krishnendu = readdlm("krishnendu_solar_model/model0.2to1")[:,1][end:-1:1]

# 	@load joinpath(Gfn_path_src,"parameters.jld2") Nν_Gfn ℓ_arr num_procs

# 	proc_id = get_processor_id_from_split_array(ℓ_arr,1:Nν_Gfn,(ℓ,ν_ind),num_procs)
# 	Gfn_fits = HelioseismicKernels.Gfn_fits_files(Gfn_path_src,(proc_id,))
# 	G0_jishnu = HelioseismicKernels.read_Gfn_file_at_index(Gfn_fits,
# 	   ℓ_arr,1:Nν_Gfn,(ℓ,ν_ind),num_procs,:,1,1,1)
# 	G1_jishnu = HelioseismicKernels.read_Gfn_file_at_index(Gfn_fits,
# 	   ℓ_arr,1:Nν_Gfn,(ℓ,ν_ind),num_procs,:,2,1,1)

# 	# display(hcat(ξ_kri,G0_jishnu))

# 	ax1=subplot(221)
# 	plot(r_frac,real(ξ_kri),label="Mandal et al. [2017]",color="grey")
# 	plot(r_frac,real(G0_jishnu),markevery=20,marker=".",
# 		ls="None",label="This work",color="black")
# 	xlim(0.995,r_frac[end]);
# 	ylabel(L"\Re[G^{(-1)}_{(-1),\ell\omega}]",fontsize=12)
# 	legend(loc="best")
# 	xlabel(L"r/R_\odot",fontsize=12)
# 	ax1.xaxis.set_major_locator(ticker.MaxNLocator(3))

# 	ax2=subplot(222)
# 	plot(r_frac,imag(ξ_kri),label="Mandal et al. [2017]",color="grey")
# 	plot(r_frac,imag(G0_jishnu),markevery=20,marker=".",
# 		ls="None",label="This work",color="black")
# 	xlim(0.995,r_frac[end]);
# 	legend(loc="best")
# 	ax2.xaxis.set_major_locator(ticker.MaxNLocator(3))
# 	xlabel(L"r/R_\odot",fontsize=12)
# 	ylabel(L"\Im[G^{(-1)}_{(-1),\ell\omega}]",fontsize=12)

# 	ax3=subplot(223)
# 	plot(r_frac,real(G1_kri),label="Mandal et al. [2017]",color="grey")
# 	plot(r_frac,real(G1_jishnu),markevery=20,marker=".",
# 		ls="None",label="This work",color="black")
# 	xlim(0.995,r_frac[end]);
# 	ylabel(L"\Re[G^{(+1)}_{(-1),\ell\omega}]",fontsize=12)
# 	legend(loc="best")
# 	xlabel(L"r/R_\odot",fontsize=12)
# 	ax3.xaxis.set_major_locator(ticker.MaxNLocator(3))

# 	ax4=subplot(224)
# 	plot(r_frac,imag(G1_kri),label="Mandal et al. [2017]",color="grey")
# 	plot(r_frac,imag(G1_jishnu),markevery=20,marker=".",
# 		ls="None",label="This work",color="black")
# 	xlim(0.995,r_frac[end]);
# 	legend(loc="best")
# 	xlabel(L"r/R_\odot",fontsize=12)
# 	ylabel(L"\Im[G^{(+1)}_{(-1),\ell\omega}]",fontsize=12)
# 	ax4.xaxis.set_major_locator(ticker.MaxNLocator(3))

# 	tight_layout()
# end

function plot_reciprocity(G_reciprocity;kwargs...)

	Gfn_path_src = HelioseismicKernels.Gfn_path_from_source_radius(
					get(kwargs,:r_src,r_src_default))

	@load(joinpath(Gfn_path_src,"parameters.jld2"),
		ν_arr,Nν_Gfn,ℓ_arr,num_procs)

	ℓ_range = get(kwargs,:ℓ_range,ℓ_arr)
	ν_ind_range = get(kwargs,:ν_ind_range,1:Nν_Gfn)

	vmax = maximum(G_reciprocity)
	norm= colors.LogNorm(vmin=vmax/1e4, vmax=vmax)

	subplot(131)
	pcolormesh(ℓ_range,ν_arr[ν_ind_range].*1e3,
		G_reciprocity[:,:,1],cmap="Greys",rasterized=true,norm=norm)
	title(L"|G^1_0\,(r_2,r_1)|",fontsize=12)
	xlabel("ℓ",fontsize=12)
	ylabel("Frequency [mHz]",fontsize=12)
	gca().xaxis.set_major_locator(ticker.MaxNLocator(4,integer=true))
	gca().yaxis.set_major_locator(ticker.MaxNLocator(5))
	colorbar(extend="min")

	subplot(132)
	pcolormesh(ℓ_range,ν_arr[ν_ind_range].*1e3,
		G_reciprocity[:,:,2],cmap="Greys",rasterized=true,norm=norm)
	title(L"|G^0_1\,(r_1,r_2)|",fontsize=12)
	xlabel("ℓ",fontsize=12)
	gca().xaxis.set_major_locator(ticker.MaxNLocator(4,integer=true))
	gca().yaxis.set_major_locator(ticker.MaxNLocator(5))
	colorbar(extend="min")

	subplot(133)
	Gdiff = G_reciprocity[:,:,1] .- G_reciprocity[:,:,2]
	vmax = maximum(abs.(Gdiff))
	norm= colors.LogNorm(vmin=vmax/1e4, vmax=vmax)
	pcolormesh(ℓ_range,ν_arr[ν_ind_range].*1e3,
		Gdiff,cmap="Greys",rasterized=true,norm=norm)
	title("Difference",fontsize=12)
	xlabel("ℓ",fontsize=12)
	gca().xaxis.set_major_locator(ticker.MaxNLocator(4,integer=true))
	gca().yaxis.set_major_locator(ticker.MaxNLocator(5))
	colorbar(extend="min")

	tight_layout()
end

function plot_reciprocity(;kwargs...)
	G_reciprocity = HelioseismicKernels.Gfn_reciprocity(;kwargs...)
	plot_reciprocity(G_reciprocity;kwargs...)
end

function plot_solar_damping(;fit_order=5)
	HMIfreq=readdlm("m181q.1216");
	ℓ_HMI,ν_HMI,γ_HMI,Δγ_HMI = HMIfreq[:,1],HMIfreq[:,3],HMIfreq[:,5],HMIfreq[:,10];

	# Fit modes above ℓ=11 and 2mHz<ν<4mHz
	mode_filter = (ℓ_HMI.>11) .& (ν_HMI .> 2e3) .& (ν_HMI .< 4.5e3);

	ν_HMI = ν_HMI[mode_filter];
	ℓ_HMI = ℓ_HMI[mode_filter];
	γ_HMI = γ_HMI[mode_filter];
	Δγ_HMI = Δγ_HMI[mode_filter];

	# Fit γ(ω) in Hz, the HMI frequencies are in μHz
	γ_damping = polyfit(ν_HMI.*(2π*1e-6),γ_HMI.*(2π*1e-6),fit_order);

	γ_fit = @. γ_damping(ν_HMI.*(2π*1e-6))*1e6/2π;

	ν_HMI_mHz = ν_HMI./1e3;

	ax_ℓν = subplot2grid((3,3),(1,0),colspan=2,rowspan=2)
	norm = colors.LogNorm(vmin=minimum(γ_HMI),vmax=maximum(γ_HMI))
	scatter(ℓ_HMI,ν_HMI_mHz,c=γ_HMI,cmap="Greys",s=10,norm=norm)
	sm = cm.ScalarMappable(norm=norm, cmap="Greys")
	sm.set_array(γ_HMI)
	colorbar(sm, ax=ax_ℓν,format=ticker.LogFormatter(base=10,
						labelOnlyBase=false, minor_thresholds=(10,2)),
						orientation="horizontal")
	xlabel("ℓ",fontsize=14)
	ylabel("Frequency [mHz]",fontsize=14)

	ax_ℓ = subplot2grid((3,3),(0,0),colspan=2,rowspan=1)
	plot(ℓ_HMI,γ_HMI,"o",ls="None",ms=4,mfc="0.7",mec="0.3")
	ylabel(L"$\gamma$ [$\mu$Hz]",fontsize=14);

	ax_ν = subplot2grid((3,3),(1,2),rowspan=2,colspan=1)
	plot(γ_HMI,ν_HMI_mHz,"o",ls="None",ms=4,mfc="0.3",mec="0.6")

	p = sortperm(ν_HMI_mHz)
	ν_HMI_mHz = ν_HMI_mHz[p]
	γ_fit = γ_fit[p]

	suffix = (fit_order==1) ? "st" : (fit_order==2) ? "nd" :
			(fit_order==3) ? "rd" : "th"

	plot(γ_fit,ν_HMI_mHz,color="black",label="fit")
	legend(loc="lower right",fontsize=16)

	xlabel(L"$\gamma$ [$\mu$Hz]",fontsize=16);
	cb_dummy=colorbar(sm, ax=ax_ν,orientation="horizontal")
	cb_dummy.ax.set_visible(false)

	ax_ℓν.yaxis.set_major_locator(ticker.MaxNLocator(4))
	ax_ℓν.xaxis.set_major_locator(ticker.MaxNLocator(4))

	ax_ν.yaxis.set_major_locator(ticker.NullLocator())
	ax_ℓ.xaxis.set_major_locator(ticker.NullLocator())

	gcf().set_size_inches(8,6)
	tight_layout()
	gcf().subplots_adjust(hspace=0,wspace=0)

	savefig(joinpath(SCRATCH[],"solar_damping.eps"),dpi=160)

	γ_damping
end

function plot_δG_isotropicδc(xobs,xsrc;kwargs...)
	Gfn_path_src = Greenfn_3D.Gfn_path_from_source_radius(xsrc)
	@load(joinpath(Gfn_path_src,"parameters.jld2"),ν_start_zeros,Nν_Gfn,ν_arr,ℓ_arr)

	δGFB = Greenfn_3D.δGrr_isotropicδc_firstborn_integrated_over_angle(xobs,xsrc;kwargs...)
	δGfull = Greenfn_3D.δGrr_isotropicδc_GminusG0(xobs,xsrc;kwargs...)

	ℓ_range = get(kwargs,:ℓ_range,ℓ_arr)

	subplot(211)
	plot(ν_arr.*1e3,real(δGFB[ν_start_zeros .+ (1:Nν_Gfn)]),
		label="First Born",color="grey")
	plot(ν_arr.*1e3,real(δGfull[ν_start_zeros .+ (1:Nν_Gfn)]),
		marker=".",ls="None",label="G-G₀",color="black")
	title("Change in Gᵣᵣ for ℓ=$(length(ℓ_range)==1 ? first(ℓ_range) : ℓ_range)",
		fontsize=15)
	legend(loc="best",fontsize=12)
	ylabel(L"\Re \left[\delta G_{rr} \right]",fontsize=15)
	gca().xaxis.set_major_locator(ticker.MaxNLocator(5))
	gca().yaxis.set_major_locator(ticker.MaxNLocator(4))

	subplot(212)
	plot(ν_arr.*1e3,imag(δGFB[ν_start_zeros .+ (1:Nν_Gfn)]),
		label="First Born",color="grey")
	plot(ν_arr.*1e3,imag(δGfull[ν_start_zeros .+ (1:Nν_Gfn)]),
		marker=".",ls="None",label="G-G₀",color="black")
	xlabel("Frequency [mHz]",fontsize=15)
	legend(loc="best",fontsize=12)
	ylabel(L"\Im \left[\delta G_{rr} \right]",fontsize=15)
	gca().xaxis.set_major_locator(ticker.MaxNLocator(5))
	gca().yaxis.set_major_locator(ticker.MaxNLocator(4))

	tight_layout()

	savefig(joinpath(SCRATCH[],"dGrr.eps"))
end

function plot_streamfunction_meridionalflow(;kwargs...)
	fig = figure()
	axflow = subplot2grid((3,2), (0,0), rowspan=3)
	θ = LinRange(0, π, 100)
	rmin = 0.6Rsun
	rmax = 0.9Rsun
	rinds = rmin .< r .< rmax
	r_cutoff = r[rinds]./Rsun
	ψ = TravelTimeShifts.meridional_flow_ψ₂₀rθ(θ'; kwargs...)
	ψ = ψ[rinds, :]
	x = [ri*sin(θi) for ri in r_cutoff, θi in θ]
	z = [ri*cos(θi) for ri in r_cutoff, θi in θ]
	p = axflow.pcolormesh(x, z, ψ, cmap="RdBu_r", rasterized=true)
	cb = colorbar(mappable = p, ax = axflow)
	cb.ax.yaxis.set_major_locator(ticker.MaxNLocator(3))
	cb.ax.yaxis.offsetText.set_visible(false)

	axflow.xaxis.set_major_locator(ticker.NullLocator())
	axflow.yaxis.set_major_locator(ticker.NullLocator())

	θ = LinRange(0, π, 20)
	rskipinds = 40
	x = [ri*sin(θi) for ri in r_cutoff, θi in θ][1:rskipinds:end, :]
	z = [ri*cos(θi) for ri in r_cutoff, θi in θ][1:rskipinds:end, :]

	ur, uθ = TravelTimeShifts.meridional_flow_u₂₀rθ(θ'; kwargs...)
	ur = ur[rinds, :][1:rskipinds:end, :]
	uθ = uθ[rinds, :][1:rskipinds:end, :]
	ux = ur .* sin.(θ') .+ uθ .* cos.(θ')
	uz = ur .* cos.(θ') .- uθ .* sin.(θ')

	ukey = maximum(sqrt, ux.^2 .+ ux.^2)

	q = axflow.quiver(x, z, ux, uz)
	axflow.quiverkey(q, 0.8, 0.92, ukey, "$(round(ukey/100, sigdigits=2)) m/s", labelsep=0.05)
	axflow.set_ylabel("The model of meridional flow", fontsize=12)

	axtime = subplot2grid((3,2), (0,1), rowspan=2)

	axtime.xaxis.set_major_locator(ticker.MaxNLocator(4))
	axtime.yaxis.set_major_locator(ticker.MaxNLocator(3))

	arr, header = FITS(joinpath(SCRATCH_kerneldir[], "streamfunction_traveltimes.fits"), "r") do f
		read(f[1]), read_header(f[1])
	end
	θ_arr = rad2deg.(LinRange(header["THMIN"], header["THMAX"], Int(header["NTH"])))
	δτ_u = arr[:,2]
	δτ_ψ = arr[:,4]
	axtime.plot(θ_arr, δτ_u, color="black", label="Flow")
	axtime.plot(θ_arr, δτ_ψ, marker="s", ms=4, color="black", ls="None", label="Stream\nfunction")
	axtime.legend(loc="upper right", fontsize=11)

	axtime.set_ylabel("Travel-time shift [sec]", fontsize=11)
	# axtime.set_xlabel("Latitude [degrees]", fontsize=11)

	axdiff = subplot2grid((3,2), (2,1), sharex=axtime)
	δτ_reldiff = arr[:,5]
	axdiff.plot(θ_arr, δτ_reldiff, color="black")
	axdiff.set_ylabel("% Diff.", fontsize=11)
	axdiff.set_xlabel("Latitude [degrees]", fontsize=11)
	# axdiff.yaxis.set_major_formatter(ticker.PercentFormatter())
	# axdiff.yaxis.set_major_formatter(ticker.FuncFormatter((x,pos)->scipercentage(x)))
	axdiff.ticklabel_format(axis="y", style="sci", scilimits=(-2,2))
	axdiff.axhline(0, ls="dotted", color="black", zorder=0)

	fig.set_size_inches(6, 4)

	fig.tight_layout()

	replace_yaxis_offsettext(axdiff, 1.03)

	fig.tight_layout()
	fig.savefig(joinpath(SCRATCH[], "meridionalflow_traveltimes.eps"))
end

end


text:
print("Testing acqr_simula: ")
nt = 1000
seed = 99
ssm = acqr(0.8,1.5,2,3,0.5,2)
m = acqr_simula(ssm,nt,seed)
y = m["y"]

if sum(y.^2) > 1e-6
	println("OK")
end

ext: jl


text:
using WriteVTK

#=
This is the 1D WriteMesh routine
=#
function writemesh(base_name, x1; fields = (), realelems = 1:size(x1)[end])
    (Nqr, _) = size(x1)
    Nsubcells = (Nqr - 1)

    cells =
        Array{MeshCell{Array{Int, 1}}, 1}(undef, Nsubcells * length(realelems))
    for e in realelems
        offset = (e - 1) * Nqr
        for i in 1:(Nqr - 1)
            cells[i + (e - 1) * Nsubcells] =
                MeshCell(VTKCellTypes.VTK_LINE, offset .+ [i, i + 1])
        end
    end

    vtkfile = vtk_grid("$(base_name)", @view(x1[:]), cells; compress = false)
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 2D WriteMesh routine
=#
function writemesh(
    base_name,
    x1,
    x2;
    x3 = nothing,
    fields = (),
    realelems = 1:size(x1)[end],
)
    @assert size(x1) == size(x2)
    (Nqr, Nqs, _) = size(x1)
    Nsubcells = (Nqr - 1) * (Nqs - 1)

    cells =
        Array{MeshCell{Array{Int, 1}}, 1}(undef, Nsubcells * length(realelems))
    ind = LinearIndices((1:Nqr, 1:Nqs))
    for e in realelems
        offset = (e - 1) * Nqr * Nqs
        for j in 1:(Nqs - 1)
            for i in 1:(Nqr - 1)
                cells[i + (j - 1) * (Nqr - 1) + (e - 1) * Nsubcells] = MeshCell(
                    VTKCellTypes.VTK_PIXEL,
                    offset .+ ind[i:(i + 1), j:(j + 1)][:],
                )
            end
        end
    end

    if x3 == nothing
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            cells;
            compress = false,
        )
    else
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            @view(x3[:]),
            cells;
            compress = false,
        )
    end
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 3D WriteMesh routine
=#
function writemesh(
    base_name,
    x1,
    x2,
    x3;
    fields = (),
    realelems = 1:size(x1)[end],
)
    (Nqr, Nqs, Nqt, _) = size(x1)
    (Nr, Ns, Nt) = (Nqr - 1, Nqs - 1, Nqt - 1)
    Nsubcells = Nr * Ns * Nt
    cells =
        Array{MeshCell{Array{Int, 1}}, 1}(undef, Nsubcells * length(realelems))
    ind = LinearIndices((1:Nqr, 1:Nqs, 1:Nqt))
    for e in realelems
        offset = (e - 1) * Nqr * Nqs * Nqt
        for k in 1:Nt
            for j in 1:Ns
                for i in 1:Nr
                    cells[i + (j - 1) * Nr + (k - 1) * Nr * Ns + (e - 1) *
                                                                 Nsubcells] =
                        MeshCell(
                            VTKCellTypes.VTK_VOXEL,
                            offset .+ ind[i:(i + 1), j:(j + 1), k:(k + 1)][:],
                        )
                end
            end
        end
    end

    vtkfile = vtk_grid(
        "$(base_name)",
        @view(x1[:]),
        @view(x2[:]),
        @view(x3[:]),
        cells;
        compress = false,
    )
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end


text:
# TODO place this doc
# """
# Type assertions are supported in the left hand of rules
# to match and access literal values both when using classic
# rewriting and EGraph based rewriting.
# To use a type assertion pattern, add `::T` after
# a pattern variable in the `left_hand` of a rule.
# """
using Parameters

import Base.==

abstract type Rule end
# Must override
==(a::Rule, b::Rule) = false
==(a::Rule, b::Rule) = false
canprune(r::Type{<:Rule}) = false
canprune(r::T) where {T<:Rule}= canprune(T)

abstract type SymbolicRule <: Rule end

"""
Rules defined as `left_hand => right_hand` are
called *symbolic rewrite* rules. Application of a *rewrite* Rule
is a replacement of the `left_hand` pattern with
the `right_hand` substitution, with the correct instantiation
of pattern variables. Function call symbols are not treated as pattern
variables, all other identifiers are treated as pattern variables.
Literals such as `5, :e, "hello"` are not treated as pattern
variables.

```julia
Rule(:(a * b => b * a))
```
"""
struct RewriteRule <: SymbolicRule
    left::Pattern
    right::Pattern
    prune::Bool
    RewriteRule(l, r) = new(l,r,false)
    RewriteRule(l,r,p) = new(l,r,p)
end
canprune(t::Type{RewriteRule}) = true

# =============================================================================

# Only the last LHS is rewritten
struct MultiPatRewriteRule <: SymbolicRule
    left::Pattern
    right::Pattern
    # additional lhs patterns
    pats::Vector{Pattern}
end
==(a::MultiPatRewriteRule, b::MultiPatRewriteRule) = a.left == b.left &&
    all(a.pats .== b.pats) && (a.right == b.right)

abstract type BidirRule <: SymbolicRule end
==(a::BidirRule, b::BidirRule) = (a.left == b.left) && (a.right == b.right)

"""
This type of *anti*-rules is used for checking contradictions in the EGraph
backend. If two terms, corresponding to the left and right hand side of an
*anti-rule* are found in an [`EGraph`], saturation is halted immediately.
"""
struct UnequalRule <: BidirRule
    left::Pattern
    right::Pattern
end

"""
```julia
Rule(:(a * b == b * a))
```
"""
struct EqualityRule <: BidirRule
    left::Pattern
    right::Pattern
end

"""
Rules defined as `left_hand |> right_hand` are
called `dynamic` rules. Dynamic rules behave like anonymous functions.
Instead of a symbolic substitution, the right hand of
a dynamic `|>` rule is evaluated during rewriting:
matched values are bound to pattern variables as in a
regular function call. This allows for dynamic computation
of right hand sides.

Dynamic rule
```julia
Rule(:(a::Number * b::Number |> a*b))
```
"""
struct DynamicRule <: Rule
    left::Pattern
    right::Any
    patvars::Vector{PatVar} # useful set of pattern variables
    prune::Bool
    DynamicRule(l::Pattern, r, prune) = new(l, r, unique(patvars(l)), prune)
    DynamicRule(l, r) = new(l,r,false)
end
canprune(t::Type{DynamicRule}) = true

==(a::DynamicRule, b::DynamicRule) = (a.left == b.left) && (a.right == b.right)


text:
module MyPkg1

using ForwardDiff

# greet() = print("Hello World!")
include("Foo.jl")

export f1

end # module


text:
geosum(r::Real, n::Integer) = (1 - r^(n + 1)) / (1 - r)
geosum(a::Real, r::Real, n::Integer) = a * geosum(r, n)

monthly_rate(annual_rate::Real) = (1 + annual_rate) ^ (1/12) - 1

"""
Future value. Assumes payments at the beginning of the month.
"""
function fv(monthly_payment::Real, annual_rate::Real; months::Integer=0,
            init::Real=0, years::Integer=0)
  months == 0 && (months = years * 12)

  r = monthly_rate(annual_rate)
  init_future = init * (1 + r) ^ months
  return geosum(monthly_payment, 1 + r, months) - monthly_payment + init_future
end

"""
Computes monthly payments. Assumes payments at the beginning of the month.
"""
function monthly_payment(goal::Real, annual_rate::Real; months::Integer=0,
                         years::Integer=0, init::Real=0)
  months == 0 && (months = years * 12)
  r = monthly_rate(annual_rate)
  init_future = init * (1 + r) ^ months
  return (goal - init_future) / (geosum(1 + r, months) - 1)
end


text:
# Thank you Distributions.jl!
using Distributions
# normal distribution methods
function distribution_derivative(w::Normal, x::Float64)
  _sigma = std(w)
  xn = (x - mean(w)) / _sigma
  return pdf(w, xn) / _sigma
end

function peizer_pratt_method_2_inversion(z::Float64, n::Int)
  isodd(n) || error("n must be an odd number, $n not allowed")

  res = (z / (n + 1.0 / 3.0 + 0.1 / (n + 1.0)))
  res *= res
  res = exp(-res * (n + 1.0 / 6.0))
  res = 0.5 + (z > 0.0 ? 1.0 : -1.0) * sqrt(0.25 * (1.0 - res))
  return res
end


text:
using StaticArrays

struct PBIntegratedDiffusion <: ContinuousTimeProcess{Float64}
    γ::Float64
end
struct PBIntegratedDiffusionAux <: ContinuousTimeProcess{Float64}
    γ::Float64
end

PorPtilde = Union{PBIntegratedDiffusion, PBIntegratedDiffusionAux}

βu(t, x::Float64, P::PBIntegratedDiffusion) = - (x+sin(x)) + 1/2
βu(t, x::Float64, P::PBIntegratedDiffusionAux) = -x + 1/2
# not really a 'beta'

Bridge.b(t::Float64, x, P::PorPtilde) = Bridge.b!(t, x, copy(x), P)
function Bridge.b!(t, x, out, P::PorPtilde)
    out[1] = x[2]
    out[2] = βu(t, x[2], P)
    out
end

function Bridge.σ!(t, x, dm, out, P::PorPtilde)
    out[1] = 0.0
    out[2] = dm*P.γ
    out
end
function Bridge.a(t, P::PorPtilde)
    [0.0 0.0; 0.0 P.γ^2]
end
Bridge.a(t, x, P::PorPtilde) = Bridge.a(t, P::PorPtilde)

Bridge.constdiff(::PorPtilde) = true

function Bridge.B!(t, arg, out, P::PBIntegratedDiffusionAux)
    B = [0.0 1.0; 0.0 -1.0]
    out .= (B*arg)
    out
end
function BBt!(t, arg, out, P::PBIntegratedDiffusionAux)
    B = [0.0 1.0; 0.0 -1.0]
    out .= (B*arg + arg*B')
    out
end

function Bridge.dP!(t, p, out, P)
    BBt!(t, p, out, P)
    out[2,2] -= P.γ^2
    out
end

Bridge.constdiff(::PBIntegratedDiffusion) = true
Bridge.b(t::Float64, x, P::PBIntegratedDiffusion) = ℝ{2}(x[2], βu(t, x[2], P))
Bridge.σ(t, x, P::PBIntegratedDiffusion) = ℝ{2}(0.0, P.γ)
Bridge.b(t::Float64, x, P::PBIntegratedDiffusionAux) = ℝ{2}(x[2], βu(t, x[2], P))

Bridge.σ(t, P::PBIntegratedDiffusionAux) =  ℝ{2}(0.0, P.γ)
Bridge.σ(t, x, P::PBIntegratedDiffusionAux) = Bridge.σ(t, P)

Bridge.B(t, P::PBIntegratedDiffusionAux) = @SMatrix [0.0 1.0; 0.0 -1.0]
Bridge.β(t, P::PBIntegratedDiffusionAux) = ℝ{2}(0, 1/2)
Bridge.a(t, P::PBIntegratedDiffusionAux) = @SMatrix [0.0 0.0; 0.0 P.γ^2]

SKIP! = true
if !SKIP!

include("partialbridge!.jl")

@test Bridge.b!(0.0, X1.yy[1], copy( X1.yy[1]),  P) == Bridge.b(0.0, X1.yy[1], P)

@test Bridge.σ!(0.0, X1.yy[1], 1.0, copy( X1.yy[1]),  P) == Bridge.σ(0.0, X1.yy[1], P)

@test Bridge.b!(0.0, X1.yy[1], copy( X1.yy[1]),  Pt) == Bridge.b(0.0, X1.yy[1], Pt)

end

include("partialparam.jl")

# Generate Data
Random.seed!(1)

W2 = sample(tt, Wiener())

X2 = solve(Euler(), x0, W2, P)

if !SKIP!
    @test W2.yy == W1.yy
    @test X2.yy == X1.yy

end

Po2 = Bridge.PartialBridgeνH(tt, P, Pt, L, v, ϵ, Σ)

Ft = copy(Po2.ν)
Ht = copy(Po2.H)

ν, H⁺, C_ = Bridge.updateνH⁺C(L, Σ, v, ϵ)

F, H, C = Bridge.updateFHC(L, Σ, v, zero(Ft[end]), zero(Ht[end]), ϵ)
@test C ≈ C_
@test F ≈ H*ν
@test H⁺ ≈ inv(H)

Ft, Ht, C = Bridge.partialbridgeodeHνH!(Bridge.R3(), tt, Ft, Ht, Pt, (F, H, C))

@test abs(C - Po2.C) < 0.03
@test maximum(norm.(Ht .- Po2.H)./norm(Ht)) < 1e-5
@test maximum(norm.(Po2.H .* Po2.ν .- Ft)) < 0.015 #not very precise, update if cond number becomes better
@test_broken cond(Ht[1]) < 1.e7 #

#LP2 = -0.5*(x0'*Po2.H[1]*x0 - 2*x0'*Po2.F)[] - Po2.C
LP2 = -0.5*(x0'*Po2.H[1]*x0 - 2*x0'*Po2.H[1]*Po2.ν[1])[] - Po2.C

@show LP, LP2
@test abs(LP - LP2) < 0.01

# no epsilon
F, H, C = Bridge.updateFHC(L, Σ, v, zero(Ft[end]), zero(Ht[end]), 0.0)
Ft, Ht, C = Bridge.partialbridgeodeHνH!(Bridge.R3(), tt, Ft, Ht, Pt, (F, H, C))
LP3 = -0.5*(x0'*Ht[1]*x0  - 2*x0'*Ft[1])[] - C
@test abs(LP - LP2) < 0.01

Xo2 = copy(X2)
solve!(Euler(), Xo2, x0, W2, Po2)

if !SKIP!
    @test norm(Xo1.yy - Xo2.yy) <sqrt(eps())
end
# Likelihood

ll2 = llikelihood(Bridge.LeftRule(), Xo2, Po2)

#@test ll1 ≈ ll2
if !SKIP!
    @test abs(ll1 - ll2) < 0.0002
end

lls2 = Float64[]

    # MCMC parameter
function mcmc2(x0, tt, Po2)
    iterations = 1000
    subsamples = 0:100:iterations
    ρ = 0.9

    # initalization
    W = sample(tt, Wiener())
    X = solve(Euler(), x0, W, Po2)
    ll = llikelihood(Bridge.LeftRule(), X, Po2)

    acc = 0

    Wo = copy(W)
    Xo = copy(X)
    Wrho = copy(W)

    XX = Any[]
    if 0 in subsamples
        push!(XX, copy(X))
    end
    for iter in 1:iterations
        # Proposal
        sample!(Wrho, Wiener())
        Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*Wrho.yy

        solve!(Euler(), Xo, x0, Wo, Po2)
        llo = llikelihood(Bridge.LeftRule(), Xo, Po2)
        push!(lls2, llo)

        if log(rand()) <= llo - ll
            X, Xo = Xo, X
            W, Wo = Wo, W

            ll = llo
            acc += 1
        end
        if iter in subsamples
            push!(XX, copy(X))
        end
    end

end

@time @testset "MCMC" begin
    mcmc2(ℝ{2}(x0), tt, Po2)
    if !SKIP!
        @test norm(lls1 - lls2) < sqrt(eps())
    end
    #@test  norm(lls1 - lls2) < 0.005

end


text:
export prove_implication

function equivalence_rule(rule::Pair{U, V}) where {U <: AbstractExpression, V <: AbstractExpression}
    [
        rule,
        rule.second => rule.first
    ]
end

p, q = LogicalSymbol.([:p, :q])
default_replacement_rules = collect(Iterators.flatten([
    # material implication
    equivalence_rule(p → q => ¬p ∨ q),
    # material equivalence
    equivalence_rule(p ⟷ q => (p ∧ q) ∨ (¬p ∧ ¬q)),
    # demorgan's laws
    equivalence_rule(¬(p ∧ q) => ¬p ∨ ¬q),
    equivalence_rule(¬(p ∨ q) => ¬p ∧ ¬q),
    # double negation
    equivalence_rule(¬¬p => p),
]))

struct RuleApplication
    rule::Pair{AbstractExpression, AbstractExpression}
    expression::AbstractExpression
end
expression(ra::RuleApplication) = ra.expression

function prove_implication(imp_expr::LogicalExpression)
    if operation(imp_expr) != →
        throw(ErrorException("Single argument `prove_implication` requires that the provided expression is an implication (→)"))
    end

    prove_implication(left(imp_expr), right(imp_expr))
end
function prove_implication(left::AbstractExpression, right::AbstractExpression; replacement_rules=default_replacement_rules)
    visited_expressions = Set{AbstractExpression}([left])
    expression_queue = AbstractExpression[left]
    while length(expression_queue) > 0
        current_expression = first(expression_queue)
        if current_expression == right
            return true
        end
        deleteat!(expression_queue, 1)

        for rr ∈ replacement_rules
            # @info replacement_metarules rr.first hash.(keys(replacement_metarules)) |> first hash(rr.first)
            replaced_expression = recursivereplace(current_expression, rr)
            has_redundant_negations = recursivematches(replaced_expression, ¬¬¬p)
            if replaced_expression != current_expression && !has_redundant_negations && replaced_expression ∉ visited_expressions
                push!(expression_queue, replaced_expression)
                push!(visited_expressions, replaced_expression)
            end
        end
    end

    false
end


text:
module DenseTests

using Photon, Test
using Photon.Layers
import Knet
import Photon.Layers: Dense, Sequential
import Photon: hasgpu, KorA

function dense_model_cpu()

    function run_model(data)
        model = Sequential(Dense(10, :relu), Dense(1))
        model(data)
    end

    data = randn(10, 16)
    pred = run_model(data)
    @test size(pred) == (1, 16)
    @test typeof(pred) == Array{Float64,2}

    data = randn(Float32, 10, 16)
    pred = run_model(data)
    @test size(pred) == (1, 16)
    @test typeof(pred) == Array{Float32,2}

end

function dense_model_gpu()
    function run_model(data)
        model = Sequential(Dense(10), Dense(1))
        model(data)
    end

    data = KorA(randn(10, 16))
    pred = run_model(data)
    @test size(pred) == (1, 16)
    @test typeof(pred) == Knet.KnetArray{Float32,2}
end

@testset "Dense" begin
    dense_model_cpu()
    hasgpu() && dense_model_gpu()
end

end


text:
using SparseArrays, Test

#-------------------------------------------------------------------------------
#   MODEL STRUCT DEFINITION
#-------------------------------------------------------------------------------

@testset "Model Generation" begin
    G = [0.0 0.0 1.0 1.0;  0.0 1.0 0.0 1.0]
    T1 = SparseArrays.sparse(Array{Int8,2}([1 1 0 0;  0 1 1 0;  0 0 1 1]))
    T2 = SparseArrays.sparse(Array{Int8,2}([1 1 1]))

    @test isnothing(Lar.Model())

    m = Lar.Model(G)
    @test size(m, 0, 1) == 2
    @test size(m, 0, 2) == 4
    @test size(m, 1, 1) == 0
    @test size(m, 1, 2) == 4
    @test size(m, 2, 1) == 0
    @test size(m, 2, 2) == 0
    @test_throws BoundsError size(m, 3, 1)

    m = Lar.Model(G, [T1, T2])
    @test size(m, 0, 1) == 2
    @test size(m, 0, 2) == 4
    @test size(m, 1, 1) == 3
    @test size(m, 1, 2) == 4
    @test size(m, 2, 1) == 1
    @test size(m, 2, 2) == 3
    @test_throws BoundsError size(m, 3, 1)

    @test_throws ArgumentError Lar.Model(G, [T1])
    @test_throws ArgumentError Lar.Model(G, [T1, T1])
    @test_throws ArgumentError Lar.Model(G, [T1, T1, T1])

end

#-------------------------------------------------------------------------------
#   OTHER BASE REIMPLEMENTATIONS
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#   GEOMETRY MANIPULATION
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#   TOPOLOGY MANIPULATION
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#   MODEL MANIPULATION
#-------------------------------------------------------------------------------