#!/usr/bin/env python

# Authors of this code:

# Daniel A. Brandt, Ph.D., Michigan Tech Research Institute, daabrand@mtu.edu

# Aaron L. Bukowski, Ph.D., University of Michigan, abukowski@umich.edu

# Aaron J. Ridley, Ph.D., University of Michigan, ridley@umich.edu

# This file contains a suite of tools that do the following:

# 1. Download FISM2 data for a time period the user desires.

# 2. Rebin that data into the binning scheme the user desires (i.e. EUVAC-37, NEUVAC-59, or SOLOMON).

# 3. Outputs a FISM2 file with the rebinnined irradiances in the desired bins (for use by euv.cpp)

# Top-level imports:

import argparse

import numpy as np

from datetime import datetime, timedelta

import pathlib

import os, sys

import pooch

from netCDF4 import Dataset

import scipy.integrate as integ

# Directory management:

here = pathlib.Path(__file__).parent.resolve()

euvDir = here.parent.joinpath('share/run/UA/inputs')

# Physical constants:

h = 6.62607015e-34 # Planck's constant in SI units of J s

c = 299792458 # Speed of light in m s^-1

# Helper Functions:

def getFism2(dateStart, dateEnd, source, downloadDir=None):

"""

Given a starting date and an ending date, automatically download irradiance data from LISIRD for a specific source,

including FISM2 daily or FISM2 in the Standard Bands.

:param dateStart: str

The starting date for the data in YYYY-MM-DD format.

:param dateEnd: str

The ending date for the data in YYYY-MM-DD format.

:param source: str

The type of data to be obtained. Valid inputs are:

- FISM2 (for daily averages of FISM2 data)

- FISM2S (for daily averages of FISM2 standard bands, according to Solomon and Qian 2005)

:return times: ndarray

Datetime values for each spectrum.

:return wavelengths: ndarray

Wavelength bins (bin boundaries) for the spectral data.

:return irradiance: ndarray

A 2D array where each row is a spectrum at a particular time, and the columns are wavelength bands.

"""

# Converting the input time strings to datetimes:

try:

dateStartDatetime = datetime.strptime(dateStart, "%Y-%m-%d")

dateEndDatetime = datetime.strptime(dateEnd, "%Y-%m-%d")

except:

dateStartDatetime = datetime.strptime(dateStart, "%Y%m%d")

dateEndDatetime = datetime.strptime(dateEnd, "%Y%m%d")

# Check if the user has asked for a source that can be obtained:

validSources = ['FISM2', 'FISM2S']

if source not in validSources:

raise ValueError("Variable 'source' must be either 'FISM2' or 'FISM2S.")

# If the download directory is not specified, set it to the top directory that the package is in:

if downloadDir is None:

downloadDir = os.getcwd()

# Download the most recent file for the corresponding source and read it in:

if source == 'FISM2':

url = 'https://lasp.colorado.edu/eve/data_access/eve_data/fism/daily_hr_data/daily_data.nc'

fname = 'FISM2_daily_data.nc'

urlObtain(url, loc=downloadDir, fname=fname) # hash='dbee404e1c75689b47691b8a4a733236bb66abbdc0f01b8cbd8236f69fe9d469'

datetimes, wavelengths, irradiance, uncertainties = obtainFism2(os.path.join(downloadDir, fname))

else:

url = 'https://lasp.colorado.edu/eve/data_access/eve_data/fism/daily_bands/daily_bands.nc'

fname = 'FISM2_daily_bands.nc'

urlObtain(url, loc=downloadDir, fname=fname) # hash='27e3183f8ad6b289de191a63d3feada64c9d3f6b2973315ceda4a42c41638465'

datetimes, wavelengths, irradiance, uncertainties = obtainFism2(os.path.join(downloadDir, fname), bands=True)

# Subset the data according to user demands:

validInds = np.where((datetimes >= dateStartDatetime) & (datetimes <= dateEndDatetime))[0]

times = datetimes[validInds]

if source == 'FISM2S':

irradiance = irradiance[-1, validInds, :]

else:

irradiance = irradiance[validInds, :]

# Return the resulting data:

return times, wavelengths, irradiance

def obtainFism2(myFism2File, bands=False):

"""

Load in spectrum data from a FISM2 file.

:param myFism2File: str

The location of the NETCDF4 file.

:param bands: bool

If True, loads in the data segmented into the Solomon and Qian 2005 standard bands.

:return datetimes: ndarray

An array of datetimes for each TIMED/SEE spectra.

:return wavelengths: ndarray

A one-dimensional array of wavelengths at which there are irradiance values.

:return irradiances: ndarray

A two-dimensional array of irradiance values at each time.

:return uncertainties: ndarray

A two-dimensional array of irradiance uncertainty values at each time.

"""

fism2Data = Dataset(myFism2File)

wavelengths = np.asarray(fism2Data.variables["wavelength"])

if bands == True: # STANDARD BANDS

flux = np.asarray(fism2Data.variables["ssi"]) # photons/cm2/second

# bandwidths = np.asarray(fism2Data.variables['band_width'])

pFlux = flux * 1.0e4 # photons/m2/second

# Convert fluxes to irradiances:

irr = np.zeros_like(flux)

for i in range(flux.shape[1]):

irr[:, i] = spectralIrradiance(pFlux[:, i], wavelengths[i] * 10.0) # W/m^2

irradiance = np.array([flux, irr])

uncertainties = np.full_like(irradiance, fill_value=np.nan) # TODO: Replace with an estimation of uncertainty

else: # NATIVE DATA

irradiance = np.asarray(fism2Data.variables["irradiance"]) # W/m^2/nm

uncertainties = np.asarray(fism2Data.variables["uncertainty"])

dates = fism2Data.variables["date"]

datetimes = []

for i in range(len(dates)):

year = dates[i][:4]

day = dates[i][4:]

currentDatetime = (

datetime(int(year), 1, 1)

+ timedelta(int(day) - 1)

+ timedelta(hours=12)

)

datetimes.append(currentDatetime)

datetimes = np.asarray(datetimes)

return datetimes, wavelengths, irradiance, uncertainties

def rebin(fism_out, saveLoc=os.getcwd(), binning_scheme='EUVAC', zero=True):

"""

Takes the output of getFism and rebins the data into whatever format the user desires.

Args:

fism_out: arraylike

The output of getFism2. Contains 4 elements: (1) datetime values for the FISM2 spectra, (2) the wavelengths

of the spectrum, (3) the actual FISM2 irradiance spectra.

saveLoc: path

Path to save data files to. Defaults to the current working directory.

binning_scheme: str

Determines the binning scheme to be used. Valid arguments include the following:

'EUVAC' or 'Euvac' or 'euvac': Uses the 37 wavelength band scheme described in Richards, et al. 1994; doi.org/10.1029/94JA00518

'NEUVAC' or 'Neuvac' or 'neuvac': Uses the 59 wavelength band scheme described in Brandt and Ridley, 2024; doi.org/10.1029/2024SW004043

'HFG': Uses the 23 wavelength band scheme described in Solomon and Qian, 2005; https://doi.org/10.1029/2005JA011160

'SOLOMON' or 'Solomon' or 'solomon': Same situation as for argument 'HFG'.

NOTE: If 'HFG' or 'SOLOMON' is chosen, the values of fism_out must correspond to getFism2 being run with the

argument source='FISM2'. If this IS NOT the case, an error will be thrown.

zero: bool

Controls whether singular (bright) wavelength lines are set to a value of zero after they are extracted.

Default is True.

Returns:

fism2_file: str

The location of the rebinned FISM data.

fism2_data: arraylike

Contains 3 elements: (a) a list of datetimes for the data and (b) the rebinned FISM2 data.

"""

# Unpack the contents of fism_out:

datetimes, wavelengths, irradiance = fism_out

# Get the native wavelength resolution of the input data:

# nativeResolution = np.concatenate((np.diff(wavelengths), np.array([np.diff(wavelengths)[-1]])), axis=0)

nativeWavelengths = wavelengths.copy()

if binning_scheme != 'HFG' and binning_scheme != 'SOLOMON' and binning_scheme != 'Solomon' and binning_scheme != 'solomon':

if binning_scheme == 'EUVAC' or binning_scheme == 'Euvac' or binning_scheme == 'euvac':

# Grab the euv_37.csv file:

fileStr = str(euvDir.joinpath('euv.csv'))

bin_bounds = read_euv_csv_file(fileStr)

tag = '_37'

elif binning_scheme == 'NEUVAC' or binning_scheme == 'Neuvac' or binning_scheme == 'neuvac':

# Grab the euv_59.csv file:

fileStr = str(euvDir.joinpath('euv_59.csv'))

bin_bounds = read_euv_csv_file(fileStr)

tag = '_59'

else:

raise FileNotFoundError('The .csv files for specifying bin boundaries cannot be found!!')

# Perform the rebinning!

shorts = bin_bounds['short'] / 10.

longs = bin_bounds['long'] / 10.

newWaves = 0.5 * (shorts + longs)

# Instantiate the new data array:

if len(irradiance.shape) < 2:

fism2_data = np.zeros((1, newWaves.shape[0]))

else:

fism2_data = np.zeros((irradiance.shape[0], newWaves.shape[0]))

# First go through all the wavelengths that are singular

myData = irradiance

for iWave, short in enumerate(shorts):

long = longs[iWave]

if (long == short):

i = np.argmin(np.abs(wavelengths - short))

i2 = np.argmin(np.abs(nativeWavelengths - short))

try:

fism2_data[:, iWave] = myData[:, i] * (nativeWavelengths[i2 + 1] - nativeWavelengths[i2])

except:

fism2_data[:, iWave] = myData[i] * (nativeWavelengths[i2 + 1] - nativeWavelengths[i2])

if zero == True:

# Zero out bin so we don't double count it.

try:

myData[:, i] = np.zeros_like(myData[:, i])

except:

myData[i] = 0.0

# Then go through the ranges

for iWave, short in enumerate(shorts):

long = longs[iWave]

if (long != short):

d1 = np.abs(wavelengths - short)

iStart = np.argmin(d1)

d2 = np.abs(wavelengths - long)

iEnd = np.argmin(d2)

wave_int = 0.0

# For wavelengths at or below 0.2 nm, just compute the sum:

if long <= 0.2:

for i in range(iStart + 1, iEnd + 1):

fism2_data[:, iWave] += myData[:, i] * \

(wavelengths[i + 1] - wavelengths[i])

wave_int += (wavelengths[i + 1] - wavelengths[i])

else:

# For issues computing the sum, integrate instead:

try:

fism2_data[:, iWave] = integ.trapezoid(myData[:, iStart:iEnd], wavelengths[iStart:iEnd], axis=1)

except:

fism2_data[:, iWave] = integ.trapezoid(myData[iStart:iEnd], wavelengths[iStart:iEnd])

elif binning_scheme == 'HFG' or binning_scheme == 'SOLOMON' or binning_scheme == 'Solomon' or binning_scheme == 'solomon':

# Determine whether the supplied data already conforms to the Solomon and Qian binning scheme.

tag = '_solomon'

if fism_out[2].shape[1] != 23:

raise ValueError("Incorrect dimensions for element 3 of argument 'fism_out'. Dimensions must be (n,23), "

"resulting from running function 'getFism' with argument 'stanBands'=True. ")

# Should the data confirm to the proper dimensions, there is no rebinning step that needs to be done. Simply

# continue.

fism2_data = fism_out[2]

else:

# If the input irradiance data DOES NOT conform to the Solomon and Qian binning scheme, throw an error.

raise ValueError("Invalid value for argument 'binning_scheme'. Must be 'EUVAC', 'NEUVAC', 'HFG', or 'SOLOMON'.")

# Save the rebinned data to a relative path (outside the package directory) in the form of a .txt file:

fism2_file = pathlib.Path(os.getcwd()).joinpath('fism2_file'+tag+'.txt')

saveFism(fism2_data, datetimes, fism2_file)

return fism2_file, fism2_data

def saveFism(data, times, filename):

"""

Takes (rebinned) FISM2 data and saves it a .txt file at a user-defined location.

Args:

data: numpy.ndarray

Irradiance data in a nxm array where the first dimension corresponds to observations (the spectrum number)

and the second dimension corresponds to wavelengths.

times: numpy.ndarray

The time values at which each spectrum is recorded.

filename: str

The desired location where the data will be saved.

Returns:

Nothing. Simply saves a file.

"""

# A helper function for working with integers:

def numStr(num):

return ',' + str(int(num))

# Define a helper function for opening a file to write the data, in such a way as to include parent directories if

# needed:

def safe_open_w(path):

''' Open "path" for writing, creating any parent directories as needed.

(https://stackoverflow.com/questions/23793987/write-a-file-to-a-directory-that-doesnt-exist)

'''

os.makedirs(os.path.dirname(path), exist_ok=True)

return open(path, 'w')

# Open the new file and begin writing, line by line:

with safe_open_w(str(filename)) as output:

# Write the header information:

# output.write("#START\n")

# Write the irradiances themselves:

firstLine = ['%.6g' % (element) for element in data[0, :]]

firstLine_joined = ','.join(firstLine)

# The first line should always be a duplicate of the first line of data, but starting at UTC=00:00 of the first date:

output.write(str(times[0].year) + numStr(times[0].month) + numStr(

times[0].day) + ',0,0,0,' + firstLine_joined + '\n')

# The rest of the lines can be straight from the data:

for i in range(data.shape[0]):

currentLine_joined = ','.join(['%.6g' % (element) for element in data[i, :]])

output.writelines(str(times[i].year) + numStr(times[i].month) + numStr(

times[i].day) + numStr(times[i].hour) + ',0,0,' + currentLine_joined + '\n')

# The last line should occur 12 hours from the last datapoint, but have duplicate values there:

lastLine_joined = ','.join(['%.6g' % (element) for element in data[-1, :]])

lastTime = times[-1] + timedelta(hours=12)

output.write(str(lastTime.year) + numStr(lastTime.month) + numStr(

lastTime.day) + ',0,0,0,' + lastLine_joined + '\n')

print('Irradiance data saved to: ')

os.system('readlink -f '+str(filename))

return

def spectralIrradiance(photonFlux, wavelength):

"""

Convert the photon flux to the corresponding spectral irradiance, given a specific wavelength.

Args:

photonFlux: numpy.ndarray, float, or int

Photon flux in units of photons s^-1 m^-2. For a singular wavelength, units are in photons m^-2.

wavelength: wavelength: float

A specific wavelength in Angstroms.

Returns:

irradiance: numpy.ndarray or float

The corresponding spectral irradiance in units of W/m^2/nm.

"""

photonEnergy = (h*c) / (wavelength*1e-10) # Convert the wavelength in the denominator to meters.

irradiance= photonFlux * photonEnergy

return irradiance

def read_euv_csv_file(file):

"""

Originally written by Aaron J. Ridley, within the file 'fism2_process.py':

https://github.com/aaronjridley/EUV/blob/main/fism2_process.py

This file reads in binning data from a CSV file that specifies bin boundaries and cross sections for either the

EUVAC model or the NEUVAC model.

Args:

file: str

The location of the .csv file to be read.

Returns:

wavelengths: numpy.ndarray

The wavelength bin boundaries for either the EUVAC model or the NEUVAC model.

"""

fpin = open(file, 'r')

iFound = 0

afac = []

f74113 = []

for line in fpin:

aline = line.split(',')

s = aline[-1].strip().split('.')[0]

if (aline[0].strip() == "Short"):

if (s.isnumeric()):

short = np.asarray(aline[5:], dtype=float)

else:

short = np.asarray(aline[5:-1], dtype=float)

iFound += 1

if (aline[0].strip() == "Long"):

if (s.isnumeric()):

long = np.asarray(aline[5:], dtype=float)

else:

long = np.asarray(aline[5:-1], dtype=float)

if (aline[0].strip() == "F74113"):

if (s.isnumeric()):

f74113 = np.asarray(aline[5:], dtype=float)

else:

f74113 = np.asarray(aline[5:-1], dtype=float)

iFound += 1

if (aline[0].strip() == "AFAC"):

if (s.isnumeric()):

afac = np.asarray(aline[5:], dtype=float)

else:

afac = np.asarray(aline[5:-1], dtype=float)

iFound += 1

# Save and convert from Angstroms to nm (FISM is in nm)

wavelengths = {'short': short / 10.0,

'long': long / 10.0,

'afac': afac,

'f74113': f74113}

return wavelengths

def urlObtain(URL, loc=None, fname=None, hash=None):

"""

Helper function that uses Pooch to download files to a location specified by the user.

:param URL: str

The location of a file to be downloaded.

:param loc: str

The place the file will be downloaded.

:param fname: str

The name the file will have once it is downloaded.

:param hash: str

A known hash (checksum) of the file. Will be used to verify the download or check if an existing file needs to

be updated.

:return:

"""

if loc is None:

loc = os.getcwd()

if os.path.isfile(str(loc) + '/' + fname) is False:

fname_loc = pooch.retrieve(url=URL, known_hash=hash, fname=fname, path=loc)

else:

fname_loc = str(loc) + '/' + fname

return fname_loc

def get_args():

parser = argparse.ArgumentParser(description = 'Create FISM input data')

parser.add_argument('start',

help='Start date (format YYYYMMDD)',

type=str)

parser.add_argument('end',

help='End date (format YYYYMMDD)',

type=str)

parser.add_argument('-b', '--binning',

help="Binning scheme to use. Can be [solomon,neuvac,euvac] "

"(case insensitive)",

type=str, default="neuvac")

args = parser.parse_args()

return args

# Execution (testing):

if __name__ == '__main__':

# Download some FISM2 data for the time period stated by the user.

args = get_args()

dateStart = args.start

dateEnd = args.end

binning_scheme = args.binning

if binning_scheme == 'HFG' or binning_scheme == 'SOLOMON' or binning_scheme == 'Solomon' or binning_scheme == 'solomon':

# SOLOMON (STAN BANDS; b23)

fism2_out_23 = getFism2(dateStart, dateEnd, 'FISM2S', downloadDir=here)

fism2_file_23, fism2_data_23 = rebin(fism2_out_23, binning_scheme=binning_scheme)

else:

fism2_out_raw = getFism2(dateStart, dateEnd, 'FISM2', downloadDir=here)

if binning_scheme == 'NEUVAC' or binning_scheme == 'Neuvac' or binning_scheme == 'neuvac':

# NEUVAC BINS (b59)

fism2_file_59, fism2_data_59 = rebin(fism2_out_raw, binning_scheme=binning_scheme, zero=True)

else:

# EUVAC BINS (b37)

fism2_file_37, fism2_data_37 = rebin(fism2_out_raw, binning_scheme=binning_scheme, zero=True)

# Exit with a zero error code:

sys.exit(0)