In [155]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams.update({'font.size': 14})  # Make the labels larger
In [156]:
# 2(a) A crude model
arr=np.loadtxt("extrasolar_signal.txt")  
# I have downloaded the file from http://sam-dolan.staff.shef.ac.uk/mas212/data/extrasolar_signal.txt
omega = 2*np.pi / 15
tdata = arr[:,0]
vdata = arr[:,1]
ts = tdata; vs = vdata;
N = len(ts)
XX = np.zeros((N,2))
for i in range(N):
    XX[i,:] = [np.sin(omega*ts[i]), np.cos(omega*ts[i])]
XT = XX.transpose()
b = np.dot(XT, vs)
A = np.dot(XT, XX)
beta = np.linalg.solve(A, b)
print("The best-fit parameters are:")
print("\tbeta[0] = %.10f\n\tbeta[1] = %.10f" % (beta[0], beta[1]))

tlin = np.linspace(0, 60, 200)
vfit = beta[0]*np.sin(omega*tlin) + beta[1]*np.cos(omega*tlin)

# Plot the image
plt.plot(ts, vs, ".", label="Data");
plt.plot(tlin, vfit, "-", label="Model");
#plt.legend()
plt.xlabel("Time in days");
plt.ylabel("Relative velocity in m/s"); 
plt.tight_layout(pad=0.5)
plt.savefig("2a_crude_fit.pdf")
The best-fit parameters are:
	beta[0] = -0.7065417249
	beta[1] = -0.1298608881
In [157]:
# 2(b) A good fit?
# First, work out the residuals
vfit = beta[0]*np.sin(omega*ts) + beta[1]*np.cos(omega*ts)
vres = vs - vfit
plt.plot(ts, vres, '.', label="Residuals")
plt.xlabel("Time in days");
plt.ylabel("Residual (m/s)"); 
plt.tight_layout(pad=0.5)
plt.savefig("2b_residuals.pdf")

def RMSD(arr):
    """Calculate the root mean square deviation (RMSD) from the residuals."""
    return np.sqrt(np.sum(arr**2)/len(arr))

print("It does not look 'random': there is evidence of residual structure here (e.g. oscillations).")
print("The residuals are correlated: the n-th residual is not independent of the (n+1)-th and (n-1)-th residual.")
print("This suggests that a better model can be found.")
print("One measure of the 'goodness of fit' is the root mean square deviation (RMSD):")
rmsd0 = RMSD(vres)
print("\tRMSD = %.10f\n" % rmsd0)
It does not look 'random': there is evidence of residual structure here (e.g. oscillations).
The residuals are correlated: the n-th residual is not independent of the (n+1)-th and (n-1)-th residual.
This suggests that a better model can be found.
One measure of the 'goodness of fit' is the root mean square deviation (RMSD):
	RMSD = 0.2341844744

In [158]:
# 2(c) A better linear model
# First try the model suggested in the brief, with four parameters
XX = np.zeros((N,4))
for i in range(N):
    XX[i,:] = [np.sin(omega*ts[i]), np.cos(omega*ts[i]), np.sin(2*omega*ts[i]), np.cos(2*omega*ts[i])]
XT = XX.transpose()
b = np.dot(XT, vs)
A = np.dot(XT, XX)
beta = np.linalg.solve(A, b)
print("The best-fit parameters:")
print("\t".join(["{:.10f}".format(b) for b in beta]))
# Plot the image
vfit = beta[0]*np.sin(omega*tlin) + beta[1]*np.cos(omega*tlin) \
     + beta[2]*np.sin(2*omega*tlin) + beta[3]*np.cos(2*omega*tlin)
plt.plot(ts, vs, ".", label="Data");
plt.plot(tlin, vfit, "-", label="Model");
#plt.legend()
plt.xlabel("Time in days");
plt.ylabel("Relative velocity in m/s"); 
plt.tight_layout(pad=0.5)
plt.savefig("2c_linear_model.pdf")
The best-fit parameters:
-0.7065417249	-0.1288601032	-0.2327004709	-0.1010792683
In [159]:
# 2(c) continued
# Now look at the residuals
vfit = beta[0]*np.sin(omega*ts) + beta[1]*np.cos(omega*ts) \
     + beta[2]*np.sin(2*omega*ts) + beta[3]*np.cos(2*omega*ts)
vres = vs - vfit
plt.plot(ts, vres, '.', label="Residuals")
plt.xlabel("Time in days");
plt.ylabel("Residual (m/s)"); 
plt.tight_layout(pad=0.5)
plt.savefig("2c_residuals.pdf")

print("The residuals are still somewhat correlated.")
print("The root mean square deviation (RMSD) is:")
rmsd1 = RMSD(vres)
print("\t%.10f\n" % rmsd1)
print("The new model has a lower RMSD (%.3f vs %.3f)." % (rmsd0, rmsd1))
The residuals are still somewhat correlated.
The root mean square deviation (RMSD) is:
	0.1508936456

The new model has a lower RMSD (0.234 vs 0.151).
In [160]:
# 2(d) A non-linear model
from scipy.optimize import curve_fit
def f2(t, b0, b1, b2):
    return b0*np.sin(b1*t + b2)

# Choosing a good starting estimate for the parameters is essential
params0 = (0.6, omega, np.pi)  
lowerbounds = (0,0,0)
upperbounds = (2, 2*omega, 2*np.pi)
popt, pcov = curve_fit(f2, ts, vs, params0, bounds=(lowerbounds, upperbounds))
print("The best-fit parameters:")
print("\t".join(["{:.10f}".format(p) for p in popt]))
vres = vs - f2(ts, popt[0], popt[1], popt[2])
rmsd2 = RMSD(vres)
print("The RMSD is %.10f" % rmsd2)
print("This suggests that this non-linear model is not significantly better than the three-parameter linear model (%.3f)" % rmsd0)
print("and worse than the six-parameter linear model (%.3f)." % rmsd1)

# An alternative model:
def f3(t, b0, b1, b2, b3, b4, b5):
    return b0*np.sin(b1*t + b3*np.sin(b4*t + b5) + b2)

params0 = (0.6, omega, np.pi, 0.1, omega, 0)
lowerbounds = (0,0,0,0,0,0)
upperbounds = (2, 2*omega, 2*np.pi, omega, 2*omega, 2*np.pi)
popt3, pcov3 = curve_fit(f3, ts, vs, params0, bounds=(lowerbounds, upperbounds))
print("\n -- New model --")
print("The best-fit parameters:")
print("\t".join(["{:.10f}".format(p) for p in popt3]))
p = popt3
vres = vs - f3(ts, p[0], p[1], p[2], p[3], p[4], p[5])
rmsd3 = RMSD(vres)
print("The RMSD is %.10f" % rmsd3)

p = popt
vlin2 = f2(tlin, p[0], p[1], p[2])
p = popt3
vlin3 = f3(tlin, p[0], p[1], p[2], p[3], p[4], p[5])
plt.plot(ts, vs, ".", label="Data");
plt.plot(tlin, vlin2, "b-", label="Non-linear Model #1");
plt.plot(tlin, vlin3, "g-", label="Non-linear Model #2");
#plt.legend()
plt.xlabel("Time in days");
plt.ylabel("Relative velocity in m/s"); 
plt.tight_layout(pad=0.5)
plt.savefig("2d_nonlinear_model.pdf")
The best-fit parameters:
0.7203971778	0.4208663592	3.2630519183
The RMSD is 0.2334931187
This suggests that this non-linear model is not significantly better than the three-parameter linear model (0.234)
and worse than the six-parameter linear model (0.151).

 -- New model --
The best-fit parameters:
0.7977009754	0.4230453734	3.1962860431	0.4188790205	0.4225317719	0.1002164046
The RMSD is 0.1519714419
In [161]:
# 3(a) Numerical solution of ODEs with midpoint method
def dXdt(X, t, p, e):
    r = X[0]; rdot = X[1]; phi = X[2];
    return np.array([rdot, -1/r**2 + p/r**3, p**0.5/r**2])

def midpoint(dXfn, tmax, nsteps, ic, par):
    dt = tmax / nsteps;
    Xs = np.zeros((3, nsteps+1))
    ts = dt*np.arange(nsteps+1)
    Xs[:,0] = ic
    for k in range(nsteps):
        X = Xs[:,k]
        Xhalf = X + 0.5*dt*dXfn(X, ts[k], par[0], par[1])
        Xs[:,k+1] = X + dt*dXfn(Xhalf, ts[k]+dt/2.0, par[0], par[1])
    return (ts,Xs)
In [137]:
# 3(b) Orbits and signals  
tmax = 60;
nsteps = 5000;
ics = [(1.0, 0.0, 0.0), (1.0, 0.4, np.pi/6), (1.0,0.8,-np.pi/4)]

plt.figure(figsize=(5,5))
for ic in ics:
    p0, e0, phi0 = ic;
    ts, Xs = midpoint(dXdt, tmax, nsteps, [p0/(1+e0), 0.0, phi0], (p0,e0))
    rs = Xs[0,::4]; phis = Xs[2,::4]
    plt.plot(rs*np.cos(phis), rs*np.sin(phis), '-', label="e = %.1f" % e0);
plt.axis('equal');
plt.ylim(-1.5,4.5); plt.xlim(-4.5,1.5);
plt.xlabel('$r \, \cos \phi$', fontsize=16); plt.ylabel('$r \, \sin \phi$', fontsize=16);
plt.tight_layout(pad=0.5)
plt.savefig("3b_orbits.pdf")
In [162]:
# 3(b) [continued]
plt.figure(figsize=(6,5))
for ic in ics:
    p0, e0, phi0 = ic;
    ts, Xs = midpoint(dXdt, tmax, nsteps, [p0/(1+e0), 0.0, phi0], (p0,e0))
    rs = Xs[0,:]; rdots = Xs[1,:]; phis = Xs[2,:]
    vs= rdots*np.cos(phis) - p0**0.5/rs*np.sin(phis)
    plt.plot(ts, vs, '-', label="e = %.1f" % e0);

plt.xlim(0,60)
plt.xlabel('Time in days', fontsize=16); 
plt.ylabel('Relative velocity in m/s', fontsize=16);
plt.tight_layout(pad=0.5)
plt.savefig("3b_example_signals.pdf")
In [163]:
# 3(c) Fitting the physical model
from scipy.optimize import leastsq

print("Time spacing in data = %.4f" % (ts[1] - ts[0])) # The data is linearly-spaced in t.

dT = tdata[1] - tdata[0]  # The data is linearly-spaced
# This step size is too large for 'midpoint', so we will use intermediate steps

def calcV(ic, dt, tmax=60.0):
    p0, e0, phi0 = ic
    nsteps = int(tmax / dt)
    ts, Xs = midpoint(dXdt, tmax, nsteps, [p0/(1+e0), 0.0, phi0], (p0,e0))
    rs = Xs[0,:]; rdots = Xs[1,:]; phis = Xs[2,:]
    vs = rdots*np.cos(phis) - p0**0.5/rs*np.sin(phis)
    return vs

def calcVres(ic, dT = 0.3, tmax=60.0):
    n2 = 6
    dt = dT / n2
    vs = calcV(ic, dt)
    return vs[::n2] - vdata

# Before attempting to minimize, it is important to find good initial estimates.
# After some experimentation, I came to the following:
ic0 = np.array((1.5, 0.4, 0))
tfit = np.linspace(0, 60, 2400)
vfit = calcV(ic0, tfit[1]-tfit[0])
plt.plot(tfit, vfit , 'r-')
plt.plot(tdata, vdata, 'b.')
Time spacing in data = 0.0120
Out[163]:
[<matplotlib.lines.Line2D at 0x115594cf8>]
In [164]:
# 3(c). Finding the best-fit parameters.
# This may take some seconds to run. 
bestfit = leastsq(calcVres, ic0)
ic1 = bestfit[0]
print("Best fit parameters:\n\tp = %.6f,\te = %.6f,\tphi = %.6f" % tuple(ic1))
tfit = np.linspace(0, 60, 2400)
vfit = calcV(bestfit[0], tfit[1]-tfit[0])
plt.plot(tfit, vfit , 'r-')
plt.plot(tdata, vdata, 'b.')
plt.xlabel('Time in days', fontsize=16); 
plt.ylabel('Relative velocity in m/s', fontsize=16);
plt.tight_layout(pad=0.5)
plt.title("Exoplanet model: best fit")
plt.savefig("3b_example_signals.pdf")
Best fit parameters:
	p = 1.490787,	e = 0.396944,	phi = -0.002900
In [141]:
# 3(b) continued
# Now look at the residuals
vres = calcVres(ic1)
plt.plot(tdata, vres, '.', label="Residuals")
plt.xlabel("Time in days");
plt.ylabel("Residual (m/s)"); 
plt.tight_layout(pad=0.5)
plt.savefig("2b_residuals.pdf")

print("The residuals do not look correlated any more.")
print("The root mean square deviation (RMSD) is now:")
rmsdbest = RMSD(vres)
print("\t%.10f" % rmsd1)
print("The new model has a lower RMSD than the crude model (%.3f vs %.3f)." % (rmsdbest, rmsd0))
The residuals do not look correlated any more.
The root mean square deviation (RMSD) is now:
	0.1508936456
The new model has a lower RMSD than the crude model (0.102 vs 0.234).

3(c) Can we infer the planet's mass and semi-major axis? This is a little tricky, because the equations I gave you have had dimensionful constants removed. Let's start with Kepler's third law, in the form $$ a^3 = \frac{G (M+m)}{\Omega^2} $$ where $a$ is the semi-major axis, $\Omega$ is the orbital angular frequency ($\Omega = 2\pi / T$, where $T$ is the period), and $M$ and $m$ are the masses of the star and planet, respectively; $G$ is Newton's gravitational constant. We may assume that $m \ll M$.

In [165]:
G = 6.674e-11  # in SI units
M = 2e30 # The mass of the Sun in kg (we assume the star has this mass)
Omega = omega / (24*3600)  # the orbital frequency is encoded in the frequency of the signal
     # Here I have converted from 'per day' to Hz = s^{-1}. 
a3 = G*M / Omega**2
a = a3**(1/3.0)
astunit = 1.496e11 # This is the Earth-Sun distance in metres. 
print("The semi-major axis 'a' is approximately:\n\t%.3e metres, or\n\t%.3f astronomical units." % (a, a/astunit))
print("So the planet is %.1f times closer to its star than Earth is." % (astunit/a))
The semi-major axis 'a' is approximately:
	1.784e+10 metres, or
	0.119 astronomical units.
So the planet is 8.4 times closer to its star than Earth is.

To estimate the planet's mass, we can use the amplitude of the relative velocity of the star. The star is pulled by the gravitational attraction of the planet. The gravitational acceleration is $\approx G m / a^2$. The star is orbiting the centre of mass of the system, at a radius $r \approx \frac{m}{M} \times a$. Using the circular-motion formula, acceleration = $v^2 / r$. Equating these two results for acceleration, and rearranging leads to $$ m \approx v \sqrt{\frac{a M}{G}} . $$

In [166]:
v = 0.8  # approximately, in metres per second
mEarth = 5.972e24  # mass of the Earth in kg.
mass = v*np.sqrt(a*M / G)
print("The mass of the planet is approximately:")
print("\t%.3e kg, or\n\t%.3f times the mass of the Earth." % (mass, mass/mEarth))
print("\nThis is a crude estimate, based on assuming a circular orbit.")
print("With some extra 'physics', it is possible to get more accurate estimates from the semi-latus rectum 'p' and eccentricity 'e' parameters")
The mass of the planet is approximately:
	1.850e+25 kg, or
	3.097 times the mass of the Earth.

This is a crude estimate, based on assuming a circular orbit.
With some extra 'physics', it is possible to get more accurate estimates from the semi-latus rectum 'p' and eccentricity 'e' parameters
In [ ]: