Astrophysics

% A planet around a star at the origin, with GM = 1
f = @(tt, u) [u(3); u(4); -u(1)/norm(u(1:2))^3; -u(2)/norm(u(1:2))^3];
[tt, u] = ode45(f, [0 12], [1; 0; 0; 1.25]);
plot(u(:, 1), u(:, 2), 'LineWidth', 1.8), hold on
scatter(0, 0, 160, [0.95 0.65 0.1], 'filled')
scatter(u(1, 1), u(1, 2), 40, 'k', 'filled')
hold off
axis equal
title('Kepler orbit from ode45')
subtitle(sprintf('energy drift %.1e', abs(0.5*sum(u(end, 3:4).^2) - 1/norm(u(end, 1:2)) - (0.5*1.25^2 - 1))))

% Planck's law for four stellar temperatures, with Wien's peaks
h = 6.626e-34; c = 2.998e8; k = 1.381e-23;
lam = linspace(100e-9, 2500e-9, 600);           % wavelength (m)
for T = [3000 4500 5800 8000]
  B = 2*h*c^2 ./ lam.^5 ./ (exp(h*c ./ (lam*k*T)) - 1);
  plot(lam*1e9, B/1e13, 'LineWidth', 2, 'DisplayName', sprintf('%d K', T))
  hold on
  lp = 2.898e-3 / T;                              % Wien's displacement law
  scatter(lp*1e9, 2*h*c^2/lp^5/(exp(h*c/(lp*k*T)) - 1)/1e13, 30, 'k', 'filled')
end
xline(380, 'k:', 'visible'), xline(750, 'k:')
hold off
legend show
xlabel('wavelength (nm)'), ylabel('spectral radiance (10¹³ W sr⁻¹ m⁻³)')
title('Hotter stars are brighter and bluer')

% A planet with 10% of the star's radius crosses its disc
rng(6)
tt = linspace(-4, 4, 600);                       % hours from mid-transit
depth = 0.1^2; dur = 3; ing = 0.35;              % (Rp/R*)², duration and ingress (h)
flux = ones(size(tt));
inT = abs(tt) < dur/2;
limb = 0.6 + 0.5*sqrt(max(0, 1 - (2*tt(inT)/dur).^2));   % limb darkening
flux(inT) = 1 - depth * min(1, (dur/2 - abs(tt(inT))) / ing) .* limb;
obs = flux + 0.0015*randn(size(tt));

scatter(tt, obs, 6, [0.65 0.65 0.65], 'filled'), hold on
plot(tt, flux, 'LineWidth', 2)
edges = -4:0.25:4; bc = []; bm = [];              % 15-minute bins
for k = 1:numel(edges) - 1
  sel = tt >= edges(k) & tt < edges(k+1);
  bc = [bc mean(tt(sel))]; bm = [bm mean(obs(sel))];
end
errorbar(bc, bm, 0.0015/sqrt(18) * ones(size(bc)), 'o')
hold off
xlabel('hours from mid-transit'), ylabel('relative flux')
title(sprintf('Transit depth %.1f%%  →  planet radius ≈ %.2f × star', 100*depth, sqrt(depth)))

rng(12)
d = 5 + 395*rand(1, 40);                        % distance (Mpc)
v = 70*d + 600*randn(size(d));                  % velocity (km/s), with peculiar motions
p = polyfit(d, v, 1);
scatter(d, v, 40, 'filled'), hold on
x = linspace(0, 420, 10);
plot(x, polyval(p, x), 'LineWidth', 2), hold off
xlabel('distance (Mpc)'), ylabel('recession velocity (km/s)')
title(sprintf('H₀ ≈ %.1f km/s/Mpc  →  Hubble time ≈ %.1f billion years', p(1), 977.8/p(1)))

r = linspace(0.5, 30, 300);                      % radius (kpc)
G = 4.30e-6; M = 5e10;                            % G in kpc (km/s)² / M☉, disc mass in M☉
vStars = sqrt(G*M ./ r .* (1 - exp(-r/3)));      % visible matter only
vHalo = 160 * sqrt(1 - atan(r/4) ./ (r/4));      % dark-matter halo
vTotal = sqrt(vStars.^2 + vHalo.^2);
rng(1)
ro = 2:2:28;
vo = interp1(r, vTotal, ro) + 6*randn(size(ro));
plot(r, vStars, '--', 'LineWidth', 2), hold on
plot(r, vHalo, ':', 'LineWidth', 2)
plot(r, vTotal, 'LineWidth', 2.5)
errorbar(ro, vo, 8*ones(size(ro)), 'ko')
hold off
legend('stars and gas', 'dark-matter halo', 'total', 'observed')
xlabel('radius (kpc)'), ylabel('rotation speed (km/s)')
title('Flat rotation curves point to dark matter')