TableOfContents()

question(md"""
What models other than Gaussian distributions are used for noise?
""")

let 
    plt = plot(xlabel="x", ylabel="pdf(x)", xlims=(-6,6))
    plot!(plt, Normal(0,1), label="Normal")
    plot!(plt, TDist(dof_t), label="Student-t dof=" * string(dof_t))
end

let
    plt = plot(xlabel="x", ylabel="pdf(x)", xlims=(-5,5))
    plot!(plt, Normal(0,1), label="Normal")
        
    mix = MixtureModel(Normal[Normal(0.0, 1), Normal(0.0, σ_outliers)],
                        [frac_in_std_normal, 1-frac_in_std_normal])
    x = -5:0.05:5
    plot!(plt, x, pdf.(mix,x), label="α N(0,1) + (1-α) N(0,σ)")
end

question(md"""
Do traditional methods and Bayesian methods yield the same results? Or, is one preferred over the other in terms of reliability and accuracy?
""")

let
    max_xplt = 10
    min_x_plt = μ_true-5*σ
    max_x_plt = μ_true+5*σ
    plt = plot(xlims=(min_x_plt,max_x_plt))
    xplt = -max_xplt:0.05:max_xplt
    plot!(plt, xplt, pdf.(prior,xplt), label="Prior")
    plot!(plt, xplt, pdf.(posterior,xplt), label="Posterior")
    plot!(plt, xplt, pdf.(true_dist,xplt), label="True distribution")
    plot!(plt,fill(sum_xs/n_bayes_plt,2), [0,ylims(plt)[2]], label="Mean xᵢ's")
end

@bind n_bayes_plt NumberField(1:ndraw,default=1)

begin
    ndraw = 100
    x_bayes_sample = μ_true .+ σ .* randn(ndraw)
end	;

begin
    sum_xs = sum(view(x_bayes_sample,1:n_bayes_plt))
    prefac = 1/(1/σ₀^2+n_bayes_plt/σ^2) 
    μ_post =  prefac*(μ₀/σ₀^2+sum_xs/σ^2)
    σ_post =  sqrt(prefac)
    prior = Normal(μ₀,σ₀)
    true_dist = Normal(μ_true,σ)
    posterior = Normal(μ_post,σ_post)
end;

question(md"""From my understanding of periodograms, you are able to detect half periods, double periods, and even triple or quad periods.  Explain.""")

let
    x = x_alias_plt
    plt = plot(x, y1, lweight=2, label= "P")
    if show_2P
        y2 = cos.(4π.*x)
        plot!(plt, x, y2, linestyle=:dot, linecolor=2, label= "2P")
        ys = cos.(4π.*xs)
        scatter!(plt, xs[1:num_obs], ys[1:num_obs],  markercolor=2,  label= "2P")
    end
    if show_3P
        y3 = cos.(6π.*x)
        plot!(plt, x, y3, linestyle=:dot, linecolor=3, label= "3P")
        ys = cos.(6π.*xs)
        scatter!(plt, xs[1:num_obs], ys[1:num_obs],  markercolor=3,  label= "3P")
    end
    if show_Po2
        y1o2 = cos.(π.*x)
        plot!(plt, x, y1o2, linestyle=:dot, linecolor=4, label= "P/2")
        ys = cos.(π.*xs)
        scatter!(plt, xs[1:num_obs], ys[1:num_obs],  markercolor=4,  label= "P/2")
    end
    plt
end

begin
    xmax = 10
    max_num_points = 100
end;

begin
    x_alias_plt = range(0,stop=xmax, step=0.05)
    y1 = cos.(2π.*x_alias_plt)
    xs = xmax.*rand(max_num_points)
    redraw_points 
end;

question(md"""Is there a difference between a periodogram and a fourier transform, and if there is, what is the advantage of finding planets using a periodogram as opposed to a fourier transform?""")

question(md"""
- How can you remove the effects of stellar processes effects on RV on the data?
- What are ancillary indicators? 
""")

TwoColumn(md"""
A true Doppler shift due to reflex motion from a single planet:
- Results in a simple shift of the spectrum
- Causes all lines to shift by an equal $\Delta\lambda/\lambda$
- Is strictly periodic with constant amplitude.
- Is strictly periodic with constant period & phase.
""",
md"""
Spurious RV measurements due to stellar variability:
- Affects line depths & shapes
- Affects some lines more than others
- Changes in magnitude with time
- May be quaesiperiodic (e.g., rotation), but not strictly periodic.
""")

question(md"""Would using space telescopes significantly improve data analysis?""")

heatmap(x_plt,y_plt,grid_I',clims=(0,1), size=(bigfigsize,bigfigsize), title="Intensity") 

heatmap(x_plt,y_plt,grid_rv',seriescolor=cgrad(ColorSchemes.balance), size=(bigfigsize,bigfigsize), title="RV") 

heatmap(x_plt,y_plt,grid_Irv'.*gridsize^2,seriescolor=cgrad(ColorSchemes.balance), size=(bigfigsize,bigfigsize), title="Contribution to RV") 

begin
    period_in_days = 3
    orbit = KeplerianOrbit(period=period_in_days, t0=0, b=b_plt, Omega=Ω_plt )
    rv_plt = calc_all(t_plt, orbit)
    fluxes = @. ld(orbit, t_anim, r_ratio_plt)

    plt_star = plot_star_weighted_rv(t_plt)
    plot!(plt_star,map(tt->Orbits.relative_position(orbit,tt)[1],t_anim),map(tt->Orbits.relative_position(orbit,tt)[2],t_anim), size=(figsize,figsize), color=:black)
    title!(plt_star,"Weighted RV w/ Planet")
    
    plt_spectrum = plot_spectrum_frame(t_plt)
    
    path_rv = map(tt->calc_rv_inplace(tt,grid_Irv,orbit,r_ratio_plt),t_anim);
    plt_rm = plot_rm_frame(t_plt,rv_plt,path_rv)
end;

plt_star

title!(plt_spectrum,"Absorption Line Profile")

title!(plt_rm,"Rossiter-McLaughlin Effect\nOn Measured RV")

begin
    plot(t_anim.*24, fluxes, label=:none, size=(figsize,figsize//2))
    scatter!([t_plt.*24],[ld(orbit, t_plt, r_ratio_plt)], color=:blue, label=:none)
    xlabel!("Time (hr)")
    ylabel!("Flux")
    title!("Transit light curve")
end

begin
    radius_in_rearth = (r_ratio_plt/0.00916794)
    mass_in_mearth = min(radius_in_rearth^3,317.9)
    rv_amp = 0.09*mass_in_mearth * (365.25/period_in_days)^(1/3)
end;

path_rv_movie = map(tt->calc_rv_inplace(tt,grid_Irv,orbit,r_ratio_plt),t_anim);

begin
    t_frame = t_anim[mod(i_frame,length(t_anim))+1]
    rv_frame = calc_all(t_frame, orbit)

    plt_star_frame = plot_star_weighted_rv(t_plt)
    plot!(plt_star_frame,map(tt->Orbits.relative_position(orbit,tt)[1],t_anim),map(tt->Orbits.relative_position(orbit,tt)[2],t_anim), size=(figsize,figsize), color=:black)
    title!(plt_star_frame,"Weighted RV w/ Planet")
    
    plt_spectrum_frame = plot_spectrum_frame(t_plt)

    plt_rm_frame = plot_rm_frame(t_frame,rv_frame,path_rv_movie)
end;

@bind i_frame Clock(0.25, false)

speed:

plt_star_frame

plt_spectrum_frame

plt_rm_frame

LocalResource("../_assets/week5/figures/geometry.png")

LocalResource("../_assets/week5/figures/Albrecht+2021_fig4.png")

LocalResource("../_assets/week5/figures/misalign_cartoon_darker.png")

ChooseDisplayMode()

Full Width Mode       Present Mode

begin
    using Transits
    using Distributions, StatsPlots # For questions
    using Plots, ColorSchemes
    using PlutoUI, PlutoTeachingTools
end

#LocalResource("../_assets/week5/WinnFig8a.png")

#LocalResource("../_assets/week5/WinnFig8b.png")

question(str; invite="Question") = Markdown.MD(Markdown.Admonition("tip", invite, [str]))

question (generic function with 1 method)

begin
    gridsize = min(100,round(Int64,8*(1/r_ratio_plt)))
    figsize = 400
    bigfigsize = 800
end

800

begin
    num_times_anim = 100
    num_times_movie = 200
    t_anim = range(-0.06, 0.06, length=num_times_anim) # days from t0
    t_movie = range(-0.06, 0.06, length=num_times_movie) # days from t0
end

-0.06:0.0006030150753768845:0.06

begin
    u = [0.4, 0.26] # quad limb dark
    ld = PolynomialLimbDark(u)
end

PolynomialLimbDark u_n: [-1.0, 0.4, 0.26]

begin
    x_plt = -1:(1//gridsize):1
    y_plt = -1:(1//gridsize):1
    grid_I = [ I(x,y) for x in x_plt, y in y_plt ]
    grid_rv = [rv_patch(x,y) for x in x_plt, y in y_plt ]
    sum_grid_I = sum(grid_I)
    grid_Irv = grid_I .* grid_rv ./ sum_grid_I
end;

v_grid = rv_patch.(x_plt,0.);

begin
    workspace_map = zeros(size(grid_Irv))
    workspace_map_spectrum = zeros(size(grid_Irv))
    workspace_spectra = zeros(length(x_plt))
end;

function I(x::Real, y::Real; radius_ratio=1)
    r = sqrt(x^2+y^2)
    if r>1 return 0 end #missing end
    μ = cos(asin(r))
    ld(μ, radius_ratio)	
end

I (generic function with 1 method)

function rv_patch(x::Real, y::Real)
    r = sqrt(x^2+y^2)
    if r>1 return 0 end # missing end
    Rsol = 696_340_000
    sec_per_day = 24*60^2
    A = 14.713 
    B = -2.396
    C = -1.787
    vsini = (π/180) / sec_per_day * Rsol
    #sinφ = y
    #vsini *= (A + B * sinφ^2 + C * sinφ^4)
    vsini *= A # + B * sinφ^2 + C * sinφ^4)
    x*vsini
end

rv_patch (generic function with 1 method)

function get_iter_to_change(x,y,r, grid)
    x_lo = searchsortedfirst(grid,x-r-0.5*grid.step)
    x_hi = searchsortedlast(grid,x+r+0.5*grid.step)
    y_lo = searchsortedfirst(grid,y-r-0.5*grid.step)
    y_hi = searchsortedlast(grid,y+r+0.5*grid.step)
    Iterators.product(x_lo:x_hi, y_lo:y_hi)
end

get_iter_to_change (generic function with 1 method)

begin
function planet_mask(x::Real, y::Real, xx::Real,yy::Real,r::Real)
    (xx-x)^2+(yy-y)^2<r^2 ? 0 : 1
end
function planet_mask(t, orbit::Transits.Orbits.AbstractOrbit, r; grid::AbstractVector = x_plt)
    x,y,z = Orbits.relative_position(orbit,t)
    [ planet_mask(x,y, xx,yy,r) for xx in grid, yy in grid ]
end
end

planet_mask (generic function with 2 methods)

function calc_map_inplace!(map, t, grid, orb, r_ratio)
    #map .= grid.*planet_mask(t,orb,r_ratio)
    x,y,z = Orbits.relative_position(orb,t)
    iter = get_iter_to_change(x,y,r_ratio, x_plt)
    map .= grid
    for (i,j) in iter
        map[i,j] = grid[i,j]*planet_mask(x,y,x_plt[i],y_plt[j], r_ratio)
    end
    map
end

calc_map_inplace! (generic function with 1 method)

function calc_rv_inplace(t, grid_Irv, orb, r_ratio, map=workspace_map)
    calc_map_inplace!(map, t, grid_Irv, orb, r_ratio)
    sum(map)
end

calc_rv_inplace (generic function with 2 methods)

function calc_spectra_inplace!(spectrum, t, grid_I, orb, r_ratio, map=workspace_map_spectrum) 
    calc_map_inplace!(map, t, grid_I, orb, r_ratio)
    spectrum .= sum(map,dims=2)
end

calc_spectra_inplace! (generic function with 2 methods)

function calc_all(t, orbit)
    rv = calc_rv_inplace(t, grid_Irv, orbit, r_ratio_plt)
    calc_spectra_inplace!(workspace_spectra, t, grid_I, orbit, r_ratio_plt, workspace_map_spectrum) 
    rv
end

calc_all (generic function with 1 method)

function plot_star_weighted_rv(t)
    plt = plot(xlims=(-1,1), ylims=(-1,1), size=(figsize,figsize), legend=:none);
    plt = heatmap!(plt,x_plt,y_plt,workspace_map',seriescolor=cgrad(ColorSchemes.balance))
end

plot_star_weighted_rv (generic function with 1 method)

function plot_spectrum_frame(t)
    plot(v_grid, 1 .-workspace_spectra./sum(workspace_spectra), xlabel="v (m/s)", ylabel="Flux", legend=:none, size=(figsize,figsize//2))
end

plot_spectrum_frame (generic function with 1 method)

function plot_rm_frame(t, rv, path_rv)
    plt = plot(t_anim.*24,path_rv, size=(figsize,figsize//2), label=:none)
    xlabel!(plt,"Time (hr)")
    ylabel!(plt,"ΔRV (m/s)")
    scatter!(plt,[t*24],[rv], mc=:blue, label=:none)
end

plot_rm_frame (generic function with 1 method)

if make_anim
    anim_rm_rv = @animate for tt ∈ t_movie
        #plot_rm_frame(tt)
        rv_frame = calc_all(tt, orbit)
        plt_rm_frame = plot_rm_frame(tt,rv_frame,path_rv_movie)
    end
    gif(anim_rm_rv, "anim_rm_rv_curve.gif", fps = 15)
    LocalResource("anim_rm_rv_curve.gif")
end

if make_anim
    anim_rm_star = @animate for tt ∈ t_movie
        rv_frame = calc_all(tt, orbit)
        plot_star_weighted_rv(tt)
    end
    gif(anim_rm_star, "anim_rm_disk.gif", fps = 15)
    LocalResource("anim_rm_disk.gif")
end

if make_anim
    anim_rm_spectra = @animate for tt ∈ t_movie
        rv_frame = calc_all(tt, orbit)
        plot_spectrum_frame(tt)
    end
    gif(anim_rm_spectra, "anim_rm_spectra.gif", fps = 15)
    LocalResource("anim_rm_spectra.gif")
end

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

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

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

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