codeeee

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import ascii
import glob

c = 2.99792458e8 #m/sec

star1 = np.loadtxt('Ab15f.dat.fix', delimiter=',',\
               unpack=True, skiprows=1)
star2 = np.loadtxt('Bb14.dat.fix', delimiter=',',\
               unpack=True, skiprows=1)
star3 = np.loadtxt('Bb42.dat.fix', delimiter=',',\
               unpack=True, skiprows=1)

plt.plot(star1[0],star1[1])
plt.xlim([3000,9000])
plt.xlabel('Wavelength (Ångström)')
plt.ylabel('Flux density (erg/cm²/s/Ångström)')
plt.title('Spectrum of star Ab15f, type A2V')
plt.show()

plt.plot(star2[0],star2[1])
plt.xlim([3000,9000])
plt.xlabel('Wavelength (Ångström)')
plt.ylabel('Flux density (erg/cm²/s/Ångström)')
plt.title('Spectrum of star Bb14, type A0V')
plt.show()

plt.plot(star3[0],star3[1])
plt.xlim([3000,9000])
plt.xlabel('Wavelength (Ångström)')
plt.ylabel('Flux density (erg/cm²/s/Ångström)')
plt.title('Spectrum of star Bb42, type F8V')
plt.show()

#3A

vega = ascii.read('vega.fnu.csv')
vega_wavelength = 10*vega[0][:]
vega_flux = vega[1][:]*(c/((1e-9*vega[0][:])**2))*1e-10

plt.plot(vega_wavelength,vega_flux)
plt.xlim([3000,9000])
plt.xlabel('Wavelength (Ångström)')
plt.ylabel('Flux density (erg/cm²/s/Ångström)')
plt.title('Vega')
plt.show()

lambda0_U = 3659
dlambda_U = 660

lowerboundU = lambda0_U-(dlambda_U/2)
upperboundU = lambda0_U+(dlambda_U/2)

wavelengthsU = np.linspace(lowerboundU,upperboundU,100)
interpFluxU = np.interp(wavelengthsU,vega_wavelength,vega_flux)

C_U = 2.5 * np.log10(np.trapz(interpFluxU*wavelengthsU)/np.trapz(wavelengthsU))
print(C_U)

#plt.plot(wavelengthsU,interpFluxU,'-x',label='interpolant')
#plt.plot(vega_wavelength,vega_flux, 'o',label='data')
#plt.xlim([lowerboundU,upperboundU])
#plt.xlabel('Wavelength (Ångström)')
#plt.ylabel('Flux density (erg/cm²/s/Ångström)')
#plt.title('Vega in the U-band range ($\lambda$ = 3659)')
#plt.legend(loc='lower right')
#plt.show()

lambda0_B = 4582
dlambda_B = 940

lowerboundB = lambda0_B-(dlambda_B/2)
upperboundB = lambda0_B+(dlambda_B/2)

wavelengthsB = np.linspace(lowerboundB,upperboundB,100)
interpFluxB = np.interp(wavelengthsB,vega_wavelength,vega_flux)

C_B = 2.5 * np.log10(np.trapz(interpFluxB*wavelengthsB)/np.trapz(wavelengthsB))
print(C_B)

lambda0_V = 5448
dlambda_V = 880

lowerboundV = lambda0_V-(dlambda_V/2)
upperboundV = lambda0_V+(dlambda_V/2)

wavelengthsV = np.linspace(lowerboundV,upperboundV,100)
interpFluxV = np.interp(wavelengthsV,vega_wavelength,vega_flux)

C_V = 2.5 * np.log10(np.trapz(interpFluxV*wavelengthsV)/np.trapz(wavelengthsV))
print(C_V)

#Next, we read in all the files:

filenames = sorted(glob.glob('*.dat.fix'))
magnitude_U = np.zeros(92) #Empty arrays to fill with magnitudes of stars
magnitude_B = np.zeros(92)
magnitude_V = np.zeros(92)

for i in range(len(filenames)):
    star = np.loadtxt(filenames[i], delimiter=',',\
                            unpack=True, skiprows=1)

    interpFluxU = np.interp(wavelengthsU,star[0],star[1])
    m_U = -2.5 * \
    np.log10(np.trapz(interpFluxU*wavelengthsU)/np.trapz(wavelengthsU)) + C_U
    magnitude_U[i] = m_U

    interpFluxB = np.interp(wavelengthsB,star[0],star[1])
    m_B = -2.5 * \
    np.log10(np.trapz(interpFluxB*wavelengthsB)/np.trapz(wavelengthsB)) + C_B
    magnitude_B[i] = m_B

    interpFluxV = np.interp(wavelengthsV,star[0],star[1])
    m_V = -2.5 * \
    np.log10(np.trapz(interpFluxV*wavelengthsV)/np.trapz(wavelengthsV)) + C_V
    magnitude_V[i] = m_V