In [103]:

```
%matplotlib inline
```

Now we have a grasp of the basics of python, but the whole reason for downloading python in the first place was to analyze solar data. Let's take a closer look at examples of solar data analysis.

We will be using SunPy to access solar data. SunPy is a python package designed to interface between the powerful tools that exist in other Python Libraries with current repositories of solar data. With SunPy we will show how to: download solar data sets from the VSO, calibrate to industry standards, plot and overlay a time series.

One of the most common data types in solar data processing is a time series. A time series is a measurement of how one physical parameter changes as a function of time. This example shows how to fit a gaussian to a spectral line. In this example, it will be as "real world" as possible.

First, let's import some useful libraries.

In [110]:

```
from datetime import datetime, timedelta #we saw these in the last tutorial
import numpy as np
from astropy.io import fits #how to read .fits files
from astropy.modeling import models, fitting #some handy fitting tools from astropy
import matplotlib.pyplot as plt
from scipy.integrate import trapz #numerical itegration tool
import astropy.units as u #units!!
```

Next we need to load in the data set we want to work with:

In [111]:

```
filename = '/Users/mskirk/Documents/Conferences/SDO 2016/EVS_L2_2011045_01_005_01.fits'
```

In [112]:

```
hdulist = fits.open(filename)
```

So what did we get when we opened the file? Let's take a look:

In [113]:

```
len(hdulist)
```

Out[113]:

We got 4 items in the list. Lets take a look at the first one:

In [114]:

```
hdulist[0].header
```

Out[114]:

This doesn't contain any useful information. Let's take a look at the second item:

In [115]:

```
hdulist[1].header
```

Out[115]:

In [116]:

```
hdulist[2].header
```

Out[116]:

In [117]:

```
hdulist[3].header
```

Out[117]:

Let's take a look at some of the data we've got.

In [118]:

```
len(hdulist[1].data)
```

Out[118]:

In [119]:

```
hdulist[1].data["WAVELENGTH"]
```

Out[119]:

In [120]:

```
len(hdulist[3].data)
```

Out[120]:

In [121]:

```
hdulist[3].data[0]["IRRADIANCE"]
```

Out[121]:

In [122]:

```
plt.plot(hdulist[1].data["WAVELENGTH"], hdulist[3].data[0]["IRRADIANCE"])
plt.ylim(0, 0.007)
plt.xlim(5, 40)
```

Out[122]:

In [123]:

```
w = np.logical_and(hdulist[1].data["WAVELENGTH"] > 30.2, hdulist[1].data["WAVELENGTH"] < 30.6)
```

In [124]:

```
w
```

Out[124]:

In [125]:

```
wavelength = hdulist[1].data["WAVELENGTH"][w]
```

In [126]:

```
wavelength
```

Out[126]:

Now, we need to add some units to this data. The header of the file tells us what the units are:

In [127]:

```
wavelength = u.Quantity(wavelength, unit="nm")
```

Let's do the same thing to the irradiance data.

In [128]:

```
irr = hdulist[3].data[0]["IRRADIANCE"][w]
```

In [129]:

```
irr = u.Quantity(irr, unit="W/m**2")
```

What have we got? Let's plot it and take a look:

In [130]:

```
plt.plot(wavelength, irr, 'o')
```

Out[130]:

Now that we have extracted the He II line from our total spectrum, we want to fit it with a gaussian. Do do this we will make use of a couple of packages in in astropy. We will initialize the gaussian fit with some approximations (max, center, FWHM):

In [131]:

```
g_init = models.Gaussian1D(0.007, 30.4, .2)
```

Now let's define a fitting method and produce a fit:

In [132]:

```
fit_g = fitting.LevMarLSQFitter()
```

In [133]:

```
g = fit_g(g_init, wavelength, irr)
```

Let's take a look at some of the qualities of our fitted gaussian:

In [134]:

```
g.mean
```

Out[134]:

In [135]:

```
g.stddev
```

Out[135]:

In [136]:

```
g.amplitude
```

Out[136]:

In [137]:

```
g
```

Out[137]:

In [138]:

```
x = np.linspace(30.2, 30.6, 1000)
```

To find the values of our fit at each location, it is easy:

In [139]:

```
y = g(x)
```

Now we can plot it:

In [140]:

```
plt.plot(wavelength, irr, 'o')
plt.plot(x, y)
```

Out[140]:

Let's find the area under the curve we just created. We can numerically integrate it easily:

In [141]:

```
intensity = trapz(y,x)
```

In [142]:

```
intensity
```

Out[142]:

The file we downloaded is a hyper-spectrum. This means that the spectrum changes over time. Can we find how the intensity of the line changes over time using the same fitting tools that we just showcased? Sure, we just need to put everything into a loop.

In [143]:

```
x = np.linspace(30.2, 30.6, 1000) #define our wavelength array
intensity = [] # initialize the intensity array
w = np.logical_and(hdulist[1].data["WAVELENGTH"] > 30.2, hdulist[1].data["WAVELENGTH"] < 30.6) # isolate the 304 line
wavelength = hdulist[1].data["WAVELENGTH"][w] # filter all of the spectra
fit_g = fitting.LevMarLSQFitter() #fitting algorithm
for spectrum in hdulist[3].data: #looping over all the spectra (360)
irr = spectrum["IRRADIANCE"][w] #pulling out the irradiance
g_init = models.Gaussian1D(0.007, 30.4, .2) #initialize the gaussian fit
g = fit_g(g_init, wavelength, irr) #do the fitting
y = g(x) # create the curve
intensity.append(trapz(y,x)) # aggregate the results
intensity = u.Quantity(intensity, unit="nm") # adding units.
```

In [144]:

```
plt.plot(intensity)
```

Out[144]:

In [ ]:

```
```