import numpy as np
import math

def imf(m, mL=0.15, mU=100.):
    return m**(-2.35) * 0.35 / (mL**(-0.35) - mU**(-0.35))

def sspmag(logtsel, logt, m_ini, magiso, magsun=4.80):

    wt   = np.where(abs(logt-logtsel)<0.01) # Select rows with logt == logtsel
    mag_t = magiso[wt]                      # New array with magnitudes for the selected age
    m_ini_t = m_ini[wt]                     # New array with masses for the selected age

    dm   = m_ini_t[1:] - m_ini_t[:-1]       # Calculate mass steps (Delta M)
    lum_t = 10**(-0.4*(mag_t[:-1]-magsun))  # Convert from Magnitudes to luminosities
    nm   = imf(m_ini_t[:-1])                # Calculate IMF weight for each mass 

    lumssp = sum(lum_t*nm*dm)               # Calculate integrated luminosity
    magssp = -2.5*math.log10(lumssp) +magsun # Convert back to magnitude

    return magssp

# Name of file with isochrone data
fn = "isoc_z019.dat"

# Read in data from the file
logt, m_ini, mv = np.loadtxt(fn, usecols=(0,1,9), unpack=True)

# Calculate the integrated V-band magnitude for log(age) = 8.0
mvssp = sspmag(8.0, logt, m_ini, mv)
print(mvssp)