Getting Started

Fetching an image with Fido

For this demo, we’ll be using sunpy and matplotlib. Let’s do some imports:

import matplotlib.pyplot as plt
from sunpy.net import Fido, attrs as a
from sunpy.map import Map

Great! Now we’ll use Fido to search between two times for an image taken by SDO/HMI in continuum intensity mode:

# Use fido to search for an SDO HMI continuum intensity image
# between the two times below
result = Fido.search(a.Time('2013/5/13 15:55', '2013/5/13 15:55:30'),
                     a.Instrument('HMI'), a.vso.Physobs('intensity'))

The variable result stores the information about the images that match the search criteria. if we feed it to Fido.fetch, sunpy will download the image:

# Download the file(s) found to the `data` directory
downloaded_files = Fido.fetch(result, path='data/.')

We can open the downloaded image using Map, like so:

# Use Sunpy to open up the image
m = Map(downloaded_files[0])

and finally we can see what the map looks like with:

# and use sunpy's built-in features to plot it
m.peek()
plt.show()

which gives us the following awesome plot!

(Source code)

Simulating a transit

Now that we have an image of the Sun to work with, let’s simulate a transit of an Earth-like planet transiting the Sun.

As with before, we’ll have a bunch of packages to import things from:

import astropy.units as u
from astropy.constants import R_earth, R_sun
from sunpy.map import Map
from stash import simulate_lightcurve
import matplotlib.pyplot as plt

The first thing we’ll do is open up the image that we downloaded in the previous example:

# Load image from the `fetch_and_plot_example.py` script
image = Map('data/hmi_ic_45s_2013_05_13_15_56_15_tai_continuum.fits').data

That’s the image of the Sun that we’ll simulate a transit on top of. Now let’s define some characteristics of the planet that will be doing the transiting:

# Set up planet parameters
orbital_period = 365 * u.day
semimajor_axis = 1 * u.AU
impact_parameter = 0.2
R_planet = R_earth
R_star = R_sun

The impact parameter b is defined on [-1, 1], and describes how far from the center of the solar disk, in units of solar radii, the planet appears to cross over. A planet that transits the solar equator has b=0. A planet that just barely grazes the tip of the Sun has b=1. (For well-aligned planets, you can convert between the solar latitude being occulted and the impact parameter by noting that \(b = \sin \theta\), where \(\theta\) is the latitude.)

Now we call the simulate_lightcurve function to simulate a light curve of a transit of this planet on our image, image:

# Simulate a light curve for that system, return a `LightCurve` object
lc = simulate_lightcurve(image, orbital_period, semimajor_axis,
                         impact_parameter, R_planet, R_star)

We can plot the LightCurve object using the built-in plot function:

# Plot the resulting light curve
lc.plot()

# Show me the plot!
plt.show()

and we’ll see something like this:

(Source code)

Would you look at that –– the planet occulted a starspot, causing the apparent brightness to temporarily increase during the transit, because less flux was missing when the planet was over the spot, compared to when the planet is over the typically bright photosphere. Now we’re cooking!

Simulating a bunch of transits

Now this time, let’s iterate over impact parameter and see all of the different transit light curves we could get as we vary \(b \in [-1, 0]\):

# Iterate over a range of impact parameters:
for impact_parameter in np.arange(-0.8, 0, 0.05):
    # Simulate a light curve for that system, return a `LightCurve` object
    lc = simulate_lightcurve(image, orbital_period, semimajor_axis,
                             impact_parameter, R_planet, R_star)

    # Plot the resulting light curve
    lc.plot()

(Source code)

You can see that the planet occulted the big starspot at one of the impact parameters that we swept through in the for loop.