# This program will read in one of my spectrum files obtained
# with the McMath-Pierce Solar Telescope on Kitt Peak, AZ and
# make a plot of the spectrum.  In the future, I will add a
# line ID function and an equivalent width of a spectral line
# function.  This program will also investigate the use of
# the Matplotlib Pyplot savefig function to make hardcopy
# files of the plots without going through the interactive
# GUI.

# When running this program, one needs to pass the name of
# the spectrum data file (in this example `sun981.wvc' to be
# read and plotted:
#       python3 specplot sun981.wvc

# Always include these next import commands.

import sys
import numpy as np
import matplotlib.pyplot as plt

# Retrieve the data filename and make a root-name for the
# plot output files.  Then make an array of plot output
# filenames.  Assume a maximum of 20 plot files.

nfmax = 20  # Maximum number of encapsulated postscript files.

if (len(sys.argv) == 2):
        spcname = sys.argv[1]

iperiod = spcname.find('.')
rootname = spcname[0:iperiod]
fsuffix = '.eps'    # Assume encapsulated postscript files.

# Make filenames for the postscript plot files.
plname = []

for i in range(0, nfmax):
    plname.append(rootname+str(i).zfill(2)+fsuffix)

# Read in the data from the spectrum file.

rdline = ''
starname = ''

fspc = open(spcname, 'r')

# Look for the name of the object that was observed.

while rdline != '  \n':
    rdline = fspc.readline()
    qstar = rdline.find('Star:')
    if qstar > -1:
        ieql = rdline.find(' = ')
        if ieql > 0:
            starname = rdline[7:ieql]

rdline = ''
sodate = ''
snccd = ''

# Look for the obs date, CCD #, and number of data points in the spectrum.

while rdline != '  \n':
    rdline = fspc.readline()
    qnpix = rdline.find('Number of pixels:')
    qdate = rdline.find('Obs-Date:')
    qccdn = rdline.find('CCD Picture Number:')
    slen = len(rdline)
    if qnpix > -1:
        snp = rdline[19:slen]
        np = int(snp)
    if qccdn > -1:
        snccd = rdline[qccdn+20:slen]
    if qdate > -1:
        sodate = rdline[11:22]

# Look for the number of data points in the spectrum.

rdline = ''
sfocus = ''

while rdline != '  \n':
    rdline = fspc.readline()
    qfocus = rdline.find('Telescope focus:')
    if qfocus > -1:
        sfocus = rdline[24:31]
        fwhm = float(sfocus)

rdline = ''

while rdline != '> Wavelength (Angstroms)\n':
    rdline = fspc.readline()

# Read in the wavelength data.
        
qaduflux = -1
swave = ''

while qaduflux == -1:
    rdline = fspc.readline()
    qaduflux = rdline.find('ADU-Flux:')
    if qaduflux == -1:
        swave = swave + rdline

# Read in the flux data.
        
qston = -1
sflux = ''

while qston == -1:
    rdline = fspc.readline()
    qston = rdline.find('Signal-to-Noise:')
    if qston == -1:
        sflux = sflux + rdline

fspc.close()

print('Object observed: "'+starname+'"')
print('Date of Observation: '+sodate)
print('CCD Picture #: '+snccd)
print('Number of pixels in spectrum: ', np)
print('Telescope FWHM: '+sfocus+' Angstroms')

# Convert wavelengths to floats.

sswave = swave.split()
wvlen = len(sswave)
wave = []

for i in range(0, wvlen):
    wave.append(float(sswave[i]))

# Convert fluxes to floats.

ssflux = sflux.split()
fllen = len(ssflux)
flux = []

for i in range(0, fllen):
    flux.append(float(ssflux[i]))
    
# Show the spectrum.

plt.plot(wave, flux, 'k-')
plt.xlabel('Wavelength (A)', labelpad=10)
plt.ylabel('Flux (ADU)', labelpad=10)
plt.title(starname+' Observation')
plt.figtext(0.18, 0.8, 'CCD Picture #'+snccd)
plt.figtext(0.7, 0.83, sodate)

plt.show()