def sfhlum(logt, lumt, t0=1e10, taustar=1e15)
# logt = array of log(t) values
# lumt = array of SSP luminosities at log(t). Delta log(t) assumed constant!

    ltot = 0.
    for tt, ll in zip(logt, lumt):    # Loop over logt, lumt
        ti = 10**tt                   # Age of SSP
        tau = t0 - ti                 # Time from t0
        psi = math.exp(-tau/taustar)  # Star formation rate
        ltot = ltot + ll*ti*psi       # Add contribution to total luminosity
    return ltot