diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000..29f3d36
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,8 @@
+language: python
+
+python:
+ - "2.6.0"
+ - "2.7.0"
+ - "3.6.0"
+
+script: pytest
\ No newline at end of file
diff --git a/README.md b/README.md
deleted file mode 100644
index efd1af1..0000000
--- a/README.md
+++ /dev/null
@@ -1 +0,0 @@
-# Exoplanet_MCMC
diff --git a/corner_plot.png b/corner_plot.png
new file mode 100644
index 0000000..434b3c9
Binary files /dev/null and b/corner_plot.png differ
diff --git a/data/__init__.py b/data/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/data/__pycache__/__init__.cpython-35.pyc b/data/__pycache__/__init__.cpython-35.pyc
new file mode 100644
index 0000000..f6e7459
Binary files /dev/null and b/data/__pycache__/__init__.cpython-35.pyc differ
diff --git a/data/__pycache__/beta_pic_b.cpython-35.pyc b/data/__pycache__/beta_pic_b.cpython-35.pyc
new file mode 100644
index 0000000..4d91ec2
Binary files /dev/null and b/data/__pycache__/beta_pic_b.cpython-35.pyc differ
diff --git a/exoplanet_mcmc/.cache/v/cache/lastfailed b/exoplanet_mcmc/.cache/v/cache/lastfailed
new file mode 100644
index 0000000..faa52c9
--- /dev/null
+++ b/exoplanet_mcmc/.cache/v/cache/lastfailed
@@ -0,0 +1,3 @@
+{
+ "plot_test.py": true
+}
\ No newline at end of file
diff --git a/exoplanet_mcmc/__init__.py b/exoplanet_mcmc/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/exoplanet_mcmc/__pycache__/__init__.cpython-35.pyc b/exoplanet_mcmc/__pycache__/__init__.cpython-35.pyc
new file mode 100644
index 0000000..a1e5b27
Binary files /dev/null and b/exoplanet_mcmc/__pycache__/__init__.cpython-35.pyc differ
diff --git a/exoplanet_mcmc/__pycache__/exoplanet_mcmc.cpython-35.pyc b/exoplanet_mcmc/__pycache__/exoplanet_mcmc.cpython-35.pyc
new file mode 100644
index 0000000..f7321c8
Binary files /dev/null and b/exoplanet_mcmc/__pycache__/exoplanet_mcmc.cpython-35.pyc differ
diff --git a/exoplanet_mcmc/__pycache__/plot_test.cpython-35-PYTEST.pyc b/exoplanet_mcmc/__pycache__/plot_test.cpython-35-PYTEST.pyc
new file mode 100644
index 0000000..115c5c2
Binary files /dev/null and b/exoplanet_mcmc/__pycache__/plot_test.cpython-35-PYTEST.pyc differ
diff --git a/exoplanet_mcmc/__pycache__/test_case.cpython-35-PYTEST.pyc b/exoplanet_mcmc/__pycache__/test_case.cpython-35-PYTEST.pyc
new file mode 100644
index 0000000..3c5b745
Binary files /dev/null and b/exoplanet_mcmc/__pycache__/test_case.cpython-35-PYTEST.pyc differ
diff --git a/exoplanet_mcmc/exoplanet_mcmc.py b/exoplanet_mcmc/exoplanet_mcmc.py
new file mode 100644
index 0000000..54dcbd7
--- /dev/null
+++ b/exoplanet_mcmc/exoplanet_mcmc.py
@@ -0,0 +1,306 @@
+import numpy as np
+import astropy.constants as cons
+import astropy.units as u
+import emcee
+import corner
+import matplotlib.pyplot as plt
+import scipy.optimize as op
+from matplotlib.ticker import FormatStrFormatter
+
+
+def mas_unc_to_m_unc(distance, mas_unc):
+ """ Converts mas uncertainties to meters at the plane of the system.
+ Arguments:
+ distance - distance to the exoplanet in meters
+ mas_unc - the uncertainty on the observations in mas"""
+ radians_unc = mas_unc.to('rad') #convert milli-arcseconds to radians
+ m_unc = np.tan(radians_unc) * distance #convert radians to parsec in the plane of the system
+ return m_unc.to('m') #convert parsecs to meters
+
+
+def lnlike(pars, data, info):
+ """ Likelihood function that returns the chi-squared distance between the model, using a set of orbital parameters
+ passed by MCMC, and the observed values.
+ Arguments:
+ pars - passed during MCMC; tuple of unitless parameters in order:
+ total_mass [in kg]
+ eccentricity
+ semimajor_axis [in m]
+ time_of_periastron [between -1 and 1, units of period]
+ argument_of_periastron [-2pi to 2pi]
+ data - data to run mcmc on; tuple of data with units in order:
+ true_anomalies [rad]
+ radii [m]
+ info - tuple of info in order of:
+ time corresponding to data [sec]
+ sigma_theta [rad]
+ sigma_r [m]"""
+ mpm, e, a, time_periastron, arg_periastron = pars # Unpack parameters
+ times, sigma_theta, sigma_r = info # Unpacking info
+ thetas, rs = data # Unpacking data
+
+ mpm = mpm * u.kg # Total mass
+ a = a * u.m # Semimajor axis
+ T = np.sqrt(4. * np.pi ** 2. * a ** 3. / (cons.G * mpm)) # Period
+ mean_anomalies = (2. * np.pi) * (times - time_periastron * T) / T # Calculating mean anomalies
+ chis = np.zeros(len(times) * 2) # For each time, two distances to calculate
+
+ #Iterative scheme to solve Kepler's Eq. as a function of time
+ #described on Solar System Dynamics pg 36(Murray, Dermott)
+ for i in range(len(times)):
+ k = 0.85 # Suggested by Murray&Dermott
+ mean_anomaly = mean_anomalies[i]
+ eccentric_anomaly_0 = mean_anomaly + np.sign(np.sin(mean_anomaly * u.rad)) * k * e #First guess
+ eccentric_anomaly_i = eccentric_anomaly_0
+
+ for j in range(3): # For low eccentricities, takes very few (~1) iterations to converge
+ fppp = e * np.cos(eccentric_anomaly_i * u.rad)
+ fpp = e * np.sin(eccentric_anomaly_i * u.rad)
+ fp = 1. - e * np.cos(eccentric_anomaly_i * u.rad)
+ f = eccentric_anomaly_i - e * np.sin(eccentric_anomaly_i * u.rad) - mean_anomaly
+ d1 = -f / fp
+ d2 = -f / (fp + .5 * d1 * fpp)
+ d3 = -f / (fp + .5 * d2 * fpp + (1 / 6.) * d2 ** 2. * fppp)
+ eccentric_anomaly_next = eccentric_anomaly_i + d3
+ eccentric_anomaly_i = eccentric_anomaly_next
+
+ true_anomaly = 2. * np.arctan(np.sqrt((1 + e) / (1 - e)) * np.tan(eccentric_anomaly_next * u.rad / 2.)) #Angle around orbit
+ true_r = a * (1. - e ** 2.) / (1. + e * np.cos(true_anomaly + arg_periastron * u.rad)) #Radial distance
+ chis[2 * i] = (thetas[i] - true_anomaly) / sigma_theta #distance between model and observations
+ chis[2 * i + 1] = (rs[i] - true_r) / sigma_r #distance between model and observations
+ chis = np.array(chis)
+ return -0.5 * (np.sum(chis * chis)) #chi-sruared distance
+
+def random_orbit():
+ """ Instantiates a system with random argument values."""
+ sys=System(mstar=np.random.randint(1,5), mplanet=np.random.randint(1, 10), distance=5., semimajor_axis=np.random.randint(1,25), eccentricity=np.random.random(), argument_of_periastron=2.*np.pi*np.random.random(), time_periastron=0, timesteps=25)
+ return sys
+
+class System(object):
+ def __init__(self, mstar, mplanet, distance, semimajor_axis, eccentricity, time_periastron, argument_of_periastron, timesteps):
+ """ Instantiates a System
+ Arguments [all without units attached]:
+ mstar - stellar mass in units of Msun
+ mplanet - planetary mass in units of Mjup
+ distance - distance to system in parsecs
+ semimajor_axis - in au
+ eccentricity
+ time_of_periastron [between -1 and 1, units of period]
+ argument_of_periastron [-2pi to 2pi]
+ timesteps - number of timesteps in an orbit"""
+ self.mstar = mstar * cons.M_sun #Attach units
+ self.mplanet = mplanet * cons.M_jup #Attach units
+ self.distance = distance * cons.pc #Attach units
+ self.a = semimajor_axis * cons.au #Attach units
+ self.e = eccentricity
+ self.tau = time_periastron
+ self.arg = argument_of_periastron * u.rad #Attach units
+
+ mu = cons.G * (self.mstar + self.mplanet)
+ self.T = np.sqrt(4. * np.pi ** 2. * self.a ** 3. / mu)
+ n = 2. * np.pi / self.T #Mean motion
+ self.times = np.linspace(0, self.T, timesteps)
+ mean_anomalies = n * (self.times - self.tau)
+ rs = []
+ true_anomalies = []
+
+ #Iterative scheme to solve Kepler's Eq. as a function of time
+ #described on Solar System Dynamics pg 36(Murray, Dermott)
+ for i in range(timesteps):
+ k = 0.85 # Suggested by Murray&Dermott
+ mean_anomaly = mean_anomalies[i]
+ eccentric_anomaly_0 = mean_anomaly + np.sign(np.sin(mean_anomaly * u.rad)) * k * self.e # First guess
+ eccentric_anomaly_i = eccentric_anomaly_0
+
+ for i in range(3): # For low eccentricities, takes very few (~1) iterations to converge
+ fppp = self.e * np.cos(eccentric_anomaly_i * u.rad)
+ fpp = self.e * np.sin(eccentric_anomaly_i * u.rad)
+ fp = 1. - self.e * np.cos(eccentric_anomaly_i * u.rad)
+ f = eccentric_anomaly_i - self.e * np.sin(eccentric_anomaly_i * u.rad) - mean_anomaly
+ d1 = -f / fp
+ d2 = -f / (fp + .5 * d1 * fpp)
+ d3 = -f / (fp + .5 * d2 * fpp + (1 / 6.) * d2 ** 2. * fppp)
+
+ eccentric_anomaly_next = eccentric_anomaly_i + d3
+ eccentric_anomaly_i = eccentric_anomaly_next
+
+ cos_true_anomaly = (np.cos(eccentric_anomaly_next * u.rad) - self.e) / (
+ 1. - self.e * np.cos(eccentric_anomaly_next * u.rad))
+ r = self.a * (1. - self.e ** 2.) / (1. + self.e * cos_true_anomaly)
+ true_anomaly = 2. * np.arctan(
+ np.sqrt((1 + self.e) / (1 - self.e)) * np.tan(eccentric_anomaly_next * u.rad / 2.))
+
+ rs.append(r.value)
+ true_anomalies.append(true_anomaly.value)
+
+ self.rs = np.array(rs) * u.m
+ self.true_anomalies = np.array(true_anomalies) * u.rad
+
+ def plot_orbit(self):
+ """ Method to plot the orbit, with no noise added. Called on an instance of System."""
+ fig, ax = plt.subplots(1, 1)
+ ax.set_aspect('equal')
+ ax.scatter(0, 0, marker='*', color='k', s=100, label='Host Star') #Plot host star at 0,0
+ ax.scatter(self.rs.to('au') * np.cos(self.true_anomalies + self.arg),
+ self.rs.to('au') * np.sin(self.true_anomalies + self.arg), color='#ff7119', s=25, label='Exoplanet') #Plot the exoplanet position at each timestep
+ lims=ax.dataLim.get_points() #Get the extrema of the data
+ ymin=lims[0][1]
+ ymax=lims[1][1]
+ xmin=lims[0][0]
+ xmax=lims[1][0]
+ ax.set_xticks([xmin, 0, xmax]) #Only plot the extrema of the orbit and 0
+ ax.set_yticks([ymin, 0, ymax]) #Only plot the extrema of the orbit and 0
+ ax.xaxis.set_major_formatter(FormatStrFormatter('%.1f'+' AU')) #Set 1 decimal precison and add AU to labels
+ ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'+' AU')) #Set 1 decimal precison and add AU to labels
+ handles, labels=ax.get_legend_handles_labels()
+ for key in ax.spines:
+ ax.spines[key].set_visible(False) #Remove borders
+ circle=plt.Circle((0,0), radius=1., color='g', fill=False, linewidth=2) #Plot circle of 1AU for reference
+ ax.add_artist(circle)
+ ax.grid(True, alpha=0.4) #Show grid
+ ax.legend(handles=handles+[circle], labels=['Host Star', 'Exoplanet', "Earth's Orbit"], loc=1)
+ return fig
+
+ def generate_mcmc_sample(self, mas_unc, sigma_true_anomaly, sigma_mass, indices):
+ """ For an instance of system, generates a sample to run MCMC on.
+ Arguments:
+ mas_unc - uncertainty in observations in mas, no units attached
+ sigma_true_anomaly - uncertainty in observations of position angle in degrees
+ sigma_mass - uncertainty in the stellar mass as a fraction of true stellar mass
+ indices - slice of the time, true_anomaly, and radius arrays to use
+ """
+ self.sigma_r = mas_unc_to_m_unc(self.distance, mas_unc * u.mas)
+ self.sigma_true_anom = (sigma_true_anomaly * u.deg).to('rad')
+ sigma_mass = sigma_mass * self.mstar
+ self.noisy_rs = self.rs[indices[0]: indices[1]] + self.sigma_r * np.random.normal(loc=0,
+ size=(indices[1] - indices[0])) #Add gaussian noise to radial distance
+ self.noisy_true_anomalies = self.true_anomalies[indices[0]: indices[1]] + self.sigma_true_anom * np.random.normal(
+ loc=0, size=(indices[1] - indices[0])) #Add gaussian noise to angles
+ self.noisy_mass = self.mstar + sigma_mass * np.random.randn() #Add gaussian noise to stellar mass
+ self.mcmc_times = self.times[indices[0]: indices[1]] #Select times of mcmc data
+
+ def plot_mcmc_sample(self):
+ """ Method to plot the data which the MCMC will run on. Called on an instance of System for which
+ generate_mcmc_sample has been run. """
+ fig, ax = plt.subplots(1, 1)
+ ax.set_aspect('equal')
+ ax.scatter(0, 0, marker='*', color='k', s=100, label='Host Star')
+ keep=np.logical_and(self.times>=self.mcmc_times.min(),self.times<=self.mcmc_times.max())
+ ax.scatter(self.noisy_rs.to('au') * np.cos(self.noisy_true_anomalies + self.arg),
+ self.noisy_rs.to('au') * np.sin(self.noisy_true_anomalies + self.arg), color='#ff7119', s=25, label='Noise Added (Observations)')
+ ax.scatter(self.rs[keep].to('au') * np.cos(self.true_anomalies[keep] + self.arg),
+ self.rs[keep].to('au') * np.sin(self.true_anomalies[keep] + self.arg), color='k', s=25, alpha=.3, label='True Values')
+ handles, labels=ax.get_legend_handles_labels()
+ lims=ax.dataLim.get_points() #Get the extrema of the data
+ ymin=lims[0][1]
+ ymax=lims[1][1]
+ xmin=lims[0][0]
+ xmax=lims[1][0]
+ ax.set_xticks([xmin, 0, xmax]) #Only plot the extrema of the orbit and 0
+ ax.set_yticks([ymin, 0, ymax])
+ ax.xaxis.set_major_formatter(FormatStrFormatter('%.1f'+' AU')) #Set 1 decimal precison and add AU to labels
+ ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'+' AU'))
+ for key in ax.spines:
+ ax.spines[key].set_visible(False) #Remove borders
+ circle=plt.Circle((0,0), radius=1., color='g', fill=False, linewidth=2) #Plot circle of 1AU for reference
+ ax.add_artist(circle)
+ ax.grid(True, alpha=0.4) #Show grid
+ ax.legend(handles=handles+[circle], labels=['Host Star', 'Noise Added (Observations)', 'True Values', "Earth's Orbit"], loc=1)
+ return fig
+
+
+ def generate_first_guess(self):
+ """ Method to generate the first guess [initial position in parameter space] for the MCMC, by minimizing the
+ likelihood function. """
+ self.data = [self.noisy_true_anomalies, self.noisy_rs] #To be passed to probability functions
+ self.info = [self.mcmc_times, self.sigma_true_anom, self.sigma_r] #To be passed to probability functions
+ nll = lambda *args: -lnlike(*args) #Dummy function to pass to minimize
+ result = op.minimize(fun=nll, x0=[(self.noisy_mass + self.mplanet).value, self.e, self.a.value, self.tau, self.arg.value],
+ args=(self.data, self.info), bounds=[(0, 100.*self.noisy_mass.value), (.001, .999), (0, 100*self.a.value), (-1,1),(-2.*np.pi, 2*np.pi)])
+ #Minimize within reasonable bounds
+ minv = result["x"] #Call the resulting parameters
+ minv = np.array(minv, dtype=float)
+ self.pos0 = minv
+
+ def lnprior(self, pars):
+ """ Flat priors in each parameter, inspired by those used for Beta Pic b in (Wang et al. 2016: 1607.05272)
+ Arguments:
+ pars - passed during MCMC; tuple of UNITLESS parameters in order:
+ total_mass [in kg]
+ eccentricity
+ semimajor_axis [in m]
+ time_of_periastron [between -1 and 1]
+ argument_of_periastron [-2pi to 2pi]"""
+ mpm, e, a, time_periastron, arg_periastron = pars #Unpack parameters
+ #Flat priors in all 5 parameters
+ if ((self.noisy_mass * .5).value < mpm < (self.noisy_mass * 2.).value) and (0.0001 < e < 0.99) and (
+ (self.a / 2.).value < a < (self.a * 2.).value) and (-1 < time_periastron < 1) and (
+ -2 * np.pi < arg_periastron < 2 * np.pi):
+ return 0.0
+ return -np.inf
+
+ def lnprob(self, pars, data, info):
+ """ Likelihood function that returns the chi-squared distance between the model, using a set of orbital parameters passed by MCMC, and the observed values.
+ Arguments:
+ pars - passed during MCMC; tuple of unitless parameters in order:
+ total_mass [in kg]
+ eccentricity
+ semimajor_axis [in m]
+ time_of_periastron [between -1 and 1, units of period]
+ argument_of_periastron [-2pi to 2pi]
+ data - data to run mcmc on; tuple of data with units in order:
+ true_anomalies [rad]
+ radii [m]
+ info - tuple of info in order of:
+ time corresponding to data [sec]
+ sigma_theta [rad]
+ sigma_r [m]"""
+ lp = self.lnprior(pars)
+ if ~(np.isfinite(lp)):
+ return -np.inf
+ return lp + lnlike(pars, data, info)
+
+ def runmcmc(self, p0spread, nwalker, nburn, nsteps):
+ """ Method to run emcee on an instance of System for which generate_sample_data has been run.
+ Arguments:
+ p0spread - fraction of range of parameter space which is used as variance of normal distribution to sample from
+ nwalker - number of MCMC walkers
+ nburn - number of burn-in steps
+ nsteps - number of MCMC steps
+ Note: For the EnsembleSampler default affine invariant algorithm, it is good to use many walkers and fewer steps."""
+ ndim = len(self.pos0) #Dimensions are length of initial parameter space position vector
+ pos = np.array([self.pos0 + p0spread * np.array(
+ [1.5 * self.noisy_mass.value, 1., 1.5 * self.a.value, 2., 4. * np.pi]) * np.random.randn(ndim) for i in
+ range(nwalker)]) #Define initial position for each walker using gaussian noise with variance of p0spread * (range of prior in each dimension)
+ for i in range(ndim):
+ print(pos[:,i].min(), pos[:,i].max()) #Print the min and max of each parameter
+ sampler = emcee.EnsembleSampler(nwalkers=nwalker, dim=ndim, lnpostfn=self.lnprob, args=(self.data, self.info))
+ for i, result in enumerate(sampler.sample(pos, iterations=nburn)):
+ if i==nburn-1: #Save end result of burn-in steps
+ newpos, prob, state = result
+ if (i+1) % (nburn/10.) == 0: #Print progress every 10 percent
+ print("Burn-in progress: {0:5.1%}".format(float(i) / nburn)) #With a total of 5 characters and 1 decimal precision, print percent progress
+ sampler.reset() #Dump burn-in information
+ for i, result in enumerate(sampler.sample(newpos, iterations=nsteps)): #Start from where the walkers ended up after burn-in
+ if (i+1) % (nsteps/10.) == 0: #Print progress every 10 percent
+ print("MCMC progress: {0:5.1%}".format(float(i) / nsteps)) #With a total of 5 characters and 1 decimal precision, print percent progress
+ self.mcmc = sampler
+ if (np.mean(sampler.acceptance_fraction) < .25) or (np.mean(sampler.acceptance_fraction) > .5):
+ print("The mean acceptance fraction was {0:.3f}. This value suggests the MCMC was not tuned well.".format(np.mean(sampler.acceptance_fraction)))
+
+
+ def walker_plot(self):
+ """ Method to plot the walkers in each parameter space. Called on an instance of System for which runmcmc has been run."""
+ f, axs = plt.subplots(self.pos0.shape[0], 1, figsize=(4, 12)) #a subplot for each dimension
+ for i, ax in enumerate(f.axes):
+ for j in range(self.mcmc.chain.shape[0]): #Plot each walker
+ ax.plot(self.mcmc.chain[j, :, i], color='k', alpha=0.2) #As indices are chain[dimension, steps, walker]
+ ax.set_ylabel(['M', 'e', 'a', 'tau', 'arg'][i])
+ return f
+
+ def corner_plot(self):
+ """ Method to plot the corner plot for the MCMC. Called on an instance of System for which runmcmc has been run."""
+ samples = self.mcmc.chain[:, :, :].reshape((-1, self.pos0.shape[0])) #Add all of the walkers together, so you have achain for each dimension
+ fig = corner.corner(samples[:, :], labels=["M", "e", "a", "tau", "arg"],
+ truths=[(self.mstar + self.mplanet).value, self.e, self.a.value, self.tau, self.arg.value])
+ return fig
diff --git a/mcmc_sample_plot.png b/mcmc_sample_plot.png
new file mode 100644
index 0000000..084d5e7
Binary files /dev/null and b/mcmc_sample_plot.png differ
diff --git a/orbit_plot.png b/orbit_plot.png
new file mode 100644
index 0000000..58fde12
Binary files /dev/null and b/orbit_plot.png differ
diff --git a/readme.md b/readme.md
new file mode 100644
index 0000000..d890eef
--- /dev/null
+++ b/readme.md
@@ -0,0 +1,114 @@
+# Exoplanet_MCMC Demonstration
+
+The goal of this notebook is to demonstrate the effect of different properties of astrometric data on the posterior distributions resulting from the MCMC package I have written. I will show the effect of having more or less densely sampled data, greater or fewer points included, and smaller or larger uncertainties on the observations.
+
+## Using the module exoplanet_mcmc
+
+To get to a corner plot of the orbital parameters of a system, only four functions must be used:
+
+### Create an instance of System
+
+```python
+print(exmc.System.__init__.__doc__)
+# Instantiates a System
+# Arguments [all without units attached]:
+# mstar - stellar mass in units of Msun
+# mplanet - planetary mass in units of Mjup
+# distance - distance to system in parsecs
+# semimajor_axis - in au
+# eccentricity
+# time_of_periastron [between -1 and 1, units of period]
+# argument_of_periastron [-2pi to 2pi]
+# timesteps - number of timesteps in an orbit
+
+sys=exmc.System(mstar=2., mplanet=2., semimajor_axis=2., eccentricity=.2, distance=5., argument_of_periastron=np.pi/4., time_periastron=0, timesteps=40)
+
+sys.plot_orbit()
+```
+
+
+
+### Pick a sample of the orbit to add noise to and use as observations for the MCMC
+
+```python
+print(exmc.System.generate_mcmc_sample.__doc__)
+# For an instance of system, generates a sample to run MCMC on.
+# Arguments:
+# mas_unc - uncertainty in observations in mas, no units attached
+# sigma_true_anomaly - uncertainty in observations of position angle in degrees
+# sigma_mass - uncertainty in the stellar mass as a fraction of true stellar mass
+# indices - slice of the time, true_anomaly, and radius arrays to use
+
+sys.generate_mcmc_sample(indices=[10,20], mas_unc=5., sigma_true_anomaly=3., sigma_mass=.03)
+
+sys.plot_mcmc_sample()
+```
+
+
+
+### Generate the first guess of the parameters for the sample observations
+
+```python
+print(exmc.System.generate_first_guess.__doc__)
+# Method to generate the first guess [initial position in parameter space] for the MCMC, by minimizing the
+# likelihood function.
+
+sys.generate_first_guess()
+```
+
+### Run the MCMC on the sample of the data
+
+```python
+print(exmc.System.runmcmc.__doc__)
+# Method to run emcee on an instance of System for which generate_sample_data has been run.
+# Arguments:
+# p0spread - fraction of range of parameter space which is used as variance of normal distribution to sample from
+# nwalker - number of MCMC walkers
+# nburn - number of burn-in steps
+# nsteps - number of MCMC steps
+# Note: For the EnsembleSampler default affine invariant algorithm, it is good to use many walkers and fewer steps.
+
+sys.runmcmc(p0spread=.01, nwalker=500, nburn=200, nsteps=400)
+
+#3.8032890891e+30 4.17956644614e+30
+#0.172149097695 0.230489784341
+#287480747958.0 313270933742.0
+#-0.0728754135833 0.0707064628404
+#-0.344134986222 0.42200270733
+#Burn in progress: 9.5%
+#Burn in progress: 19.5%
+#Burn in progress: 29.5%
+#Burn in progress: 39.5%
+#Burn in progress: 49.5%
+#Burn in progress: 59.5%
+#Burn in progress: 69.5%
+#Burn in progress: 79.5%
+#Burn in progress: 89.5%
+#Burn in progress: 99.5%
+#MCMC progress: 9.8%
+#MCMC progress: 19.8%
+#MCMC progress: 29.8%
+#MCMC progress: 39.8%
+#MCMC progress: 49.8%
+#MCMC progress: 59.8%
+#MCMC progress: 69.8%
+#MCMC progress: 79.8%
+#MCMC progress: 89.8%
+#MCMC progress: 99.8%
+```
+
+## Plot the results
+
+```sys.walker_plot()```
+
+
+
+If the MCMC was run with enough steps, the plot of the walkers in each parameter space should look like "white noise", as above.
+
+```sys.corner_plot()```
+
+
+
+The true values of the parameters are overplotted in blue. We can see that they fall near the regions of maximum probability.
+
+
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..c5a0262
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,6 @@
+matplotlib
+astropy
+scipy
+emcee
+corner
+numpy
\ No newline at end of file
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/__pycache__/__init__.cpython-35.pyc b/tests/__pycache__/__init__.cpython-35.pyc
new file mode 100644
index 0000000..2540769
Binary files /dev/null and b/tests/__pycache__/__init__.cpython-35.pyc differ
diff --git a/tests/__pycache__/test_exoplanet_mcmc.cpython-35-PYTEST.pyc b/tests/__pycache__/test_exoplanet_mcmc.cpython-35-PYTEST.pyc
new file mode 100644
index 0000000..cabda75
Binary files /dev/null and b/tests/__pycache__/test_exoplanet_mcmc.cpython-35-PYTEST.pyc differ
diff --git a/tests/test_exoplanet_mcmc.py b/tests/test_exoplanet_mcmc.py
new file mode 100644
index 0000000..b67f2f1
--- /dev/null
+++ b/tests/test_exoplanet_mcmc.py
@@ -0,0 +1,26 @@
+from exoplanet_mcmc.exoplanet_mcmc import *
+import numpy as np
+
+def test_system1():
+ sys=random_orbit()
+ assert sys.mstar.value>0
+
+def test_system2():
+ sys=random_orbit()
+ assert sys.mplanet.value>0
+
+def test_system3():
+ sys=random_orbit()
+ assert (sys.e>0)&(sys.e<1)
+
+def test_system4():
+ sys=random_orbit()
+ assert sys.a.value>0
+
+def test_system5():
+ sys=random_orbit()
+ assert (sys.tau>-1)&(sys.tau<1)
+
+def test_system6():
+ sys=random_orbit()
+ assert (sys.arg.value>-2*np.pi)&(sys.arg.value<2*np.pi)
\ No newline at end of file
diff --git a/walker_plot.png b/walker_plot.png
new file mode 100644
index 0000000..22b7ec0
Binary files /dev/null and b/walker_plot.png differ