begin
    min_period = 0.75
    max_period = 1*time_span
    f_nyquist = 0.5*nobs/time_span
    P_grid = reverse(1 ./ range(1/max_period, 1/min_period, step = f_nyquist/(2π*oversamp)))
end;

t_daily = collect(1:nobs);

t_regular = collect(range(0,time_span; step=floor(time_span/nobs))[1:nobs]);

t_perturbed = collect(t_regular .+ (2.0 .*rand(length(t_regular)).-1.0).*max_Δtobs);

t_uniform = collect(sort(time_span * rand(nobs)));

t_plt = 1:0.5:time_span;

begin
    regenerate_residuals
    rv_residuals = σ_rv .* randn(nobs);
end;

rv_obs_daily = rv_true.(period,K,phase,t_daily) .+ rv_residuals;

rv_obs_regular = rv_true.(period,K,phase,t_regular) .+ rv_residuals;

rv_obs_perturbed = rv_true.(period,K,phase,t_perturbed) .+ rv_residuals;

rv_obs_uniform = rv_true.(period,K,phase,t_uniform) .+ rv_residuals;

rv_plt = rv_true.(period,K,phase,t_plt);

output_daily = periodogram(rv_obs_daily, t_daily, P_grid);

output_regular = periodogram(rv_obs_regular, t_regular, P_grid);

output_perturbed = periodogram(rv_obs_perturbed, t_perturbed, P_grid);

output_uniform = periodogram(rv_obs_uniform, t_uniform, P_grid);

rv_true(period::Real,K::Real,ϕ::Real, t::Real) = K * cos(2π*t/period + ϕ)

rv_true (generic function with 1 method)

function build_design_matrix_sinusoidal!(A::AbstractMatrix, P::Real, t::AbstractVector{<:Real})
    @assert size(A,1) == length(t)
    @assert size(A,2) == 3
    A[:,1] .= cos.(2π/P .* t)
    A[:,2] .= sin.(2π/P .* t)
    A[:,3] .= 1
    return A
end

build_design_matrix_sinusoidal! (generic function with 1 method)

function periodogram(y::AbstractVector, t::AbstractVector,
                     P_grid::AbstractVector)
    @assert length(t) == length(y)
    n = length(y)
    A = Matrix{Float64}(undef,n,3)
    predict = Vector{Float64}(undef,n)
    periodogram_rss = Vector{Float64}(undef,length(P_grid))
    periodogram_K = Vector{Float64}(undef,length(P_grid))
    
    for (i,P) in enumerate(P_grid)
        build_design_matrix_sinusoidal!(A,P,t)
        sol = solve(LinearProblem(A, y))
        predict .= A*sol
        periodogram_K[i] = sqrt(sol[1]^2+sol[2]^2)
        periodogram_rss[i] = sum((predict.-y).^2)
    end
    rss0 = sum(y.^2)
    periodogram_rss_norm = (rss0.-periodogram_rss)./rss0
    idx_bf = argmax(periodogram_rss)
    P_bf = P_grid[idx_bf]
    build_design_matrix_sinusoidal!(A,P_bf,t)
    sol_bf = solve(LinearProblem(A, y))
    K_bf = sqrt(sol_bf[1]^2+sol_bf[2]^2)
    pred_bf = A*sol_bf
    return (;periodogram_rss, periodogram_rss_norm, periodogram_K, P_bf, K_bf, pred_bf, idx_bf, sol_bf)
end

periodogram (generic function with 1 method)

begin
    sched_list = [:daily => "Daily", :regular => "Regular", :perturbed=> "Perturbed", :uniform => "Uniform", :none=>"None"]
    sched_dict = Dict(first.(sched_list) .=> last.(sched_list))
end;

begin
    t_sched = (;daily=t_daily, regular=t_regular, perturbed=t_perturbed, uniform=t_uniform)
    rv_obs_sched = (;daily=rv_obs_daily, regular=rv_obs_regular, perturbed=rv_obs_perturbed, uniform=rv_obs_uniform)
    output_sched = (;daily=output_daily, regular=output_regular, perturbed=	output_perturbed, uniform=output_uniform)
end;

begin
xmin = exp(log_xmin)
xmax = exp(log_xmax)
end;

phase = deg2rad(phase_deg);

begin
    #σ_tobs = σ_tobs_hr/24
    max_Δtobs = max_Δtobs_hr/24
end;

time_span = time_span_yr * 365.2425;

begin
    using LinearAlgebra, LinearSolve
    using Random, Distributions
    using Plots, LaTeXStrings
    using PlutoUI, PlutoTeachingTools
end

TableOfContents()

To run this notebook locally, download this file and open it with Pluto.jl.