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!
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:
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()
You can see that the planet occulted the big starspot at one of the impact
parameters that we swept through in the for
loop.