In [457]:

    

import numpy as np

import scipy.integrate as integ
import scipy.interpolate as interp

import matplotlib.pyplot as plt
import matplotlib.patches as patch


    

In [458]:

    

# Loading the wavelength and flux of the spectrum
mat = np.loadtxt('NEWNGC1850SOAR.txt')
wavelength, flux = mat.transpose()


    

In [459]:

    

# Spectral plot
plt.plot(wavelength, flux)


    

    
Out[459]:





[<matplotlib.lines.Line2D at 0x7f1a5082c390>]






    

In [460]:

    

# G Band
wl_count = 4000
wl_start, wl_end = 4318, 4364
wl_chosen_start, wl_chosen_end = 4324, 4364
wl_center = wl_chosen_start + (wl_chosen_end - wl_chosen_start) / 2


    

In [461]:

    

# Isolate G Band
gstart_index, gend_index = np.searchsorted(wavelength, [wl_start-2, wl_end+2])
wl_gband = wavelength[gstart_index:gend_index]
flux_gband = flux[gstart_index:gend_index]


    

In [462]:

    

# Plot G Band
plt.plot(wl_gband, flux_gband)


    

    
Out[462]:





[<matplotlib.lines.Line2D at 0x7f1a50759d10>]






    

In [463]:

    

# Cubic spline interpolation of the spectrum within the chosen band
spline = interp.interp1d(wl_gband, flux_gband, kind='cubic')
wl_gband_fine = np.linspace(wl_start, wl_end, wl_count)
flux_gband_fine = spline(wl_gband_fine)
plt.plot(wl_gband_fine, flux_gband_fine)


    

    
Out[463]:





[<matplotlib.lines.Line2D at 0x7f1a506a3350>]






    

In [464]:

    

# Fit a line through the chosen points to determine the continuum level
wl_chosen = [wl_chosen_start, wl_chosen_end]
flux_chosen = spline(wl_chosen)
plt.plot(wl_gband_fine, flux_gband_fine, '-',
         wl_chosen, flux_chosen, 'o-')


    

    
Out[464]:





[<matplotlib.lines.Line2D at 0x7f1a505de750>,
 <matplotlib.lines.Line2D at 0x7f1a505de9d0>]






    

In [465]:

    

# Compute the area between the continuum level and the curve
lower_limit, upper_limit = wl_chosen
inverse_area, err = integ.quad(spline, lower_limit, upper_limit)
total_area = height * (upper_limit - lower_limit)
area = total_area - inverse_area


    

In [466]:

    

# The equivalent width is the width of a rectangle with an area equal to the calculated area
# and a height equal to the continuum level
height = np.average(flux_chosen)
width = area / height
wl_width_start, wl_width_end = wl_center-width, wl_center+width


    

In [467]:

    

# Plot the chosen band and the equivalent rectangle (shaded)
axis = plt.gca()
axis.add_patch(
    patch.Rectangle(
        (wl_width_start, 0),
        width, height,
        facecolor='grey',
        alpha=0.5
    )
)

# Plot the spectrum restricted to the points chosen
chosen_index_start, chosen_index_end = np.searchsorted(wl_gband_fine, wl_chosen)
wl_chosen_res = wl_gband_fine[chosen_index_start:chosen_index_end]
flux_chosen_res = spline(wl_chosen_res)

axis.set_title("EQW = " + str(width), fontsize=20)
axis.plot(
    wl_chosen_res, flux_chosen_res, '-',
    wl_chosen, flux_chosen, 'o-'
)


    

    
Out[467]:





[<matplotlib.lines.Line2D at 0x7f1a5051d750>,
 <matplotlib.lines.Line2D at 0x7f1a5052a150>]






    

In [467]: