import astropy.constants as constants
import astropy.units as u
import numpy  as np
from scipy.integrate import simpson

def synPhot(wave, flux, flux_err, Tw, Tf):
    '''
    Usage:   
                        - Synthetic photometry
    Input:           
        wave            - Wavelength pixels of input spectrum in unit of angstrom.
        flux            - Spectral flux in unit of erg/s/cm^2/angstrom.
        flux_err        - Flux uncertainty.
        Tw              - filter response function: wavelength in angstrom.
        Tf              - filter response function: transmission at each wavelength.
            
    Output:
        band_mag        - synthetic magnitude
        band_mag_err    - magnitude uncertainty

    Author:             He Zhao
    Email:              he.zhao@oca.eu, hzhao@pmo.ac.cn
    Finished:           2024.10.23
    Update Log:
    
    '''

    
    light_speed = constants.c.to_value(u.angstrom/u.s)
    
    Fnu = flux * (wave**2/light_speed)
    Fnu_err = flux_err * (wave**2/light_speed)
    # interpolate response function to have the same samplings as spectrum
    Tff = np.interp(wave, Tw, Tf)
    up  = simpson(Fnu*Tff/wave, wave)
    low = simpson(Tff/wave, wave)
    band_flux = up/low * 1e23
    band_flux_err = np.sqrt(np.sum(((Fnu_err*Tff/wave)/low)**2*400.)) * 1e23
    if np.isnan(band_flux) or band_flux < 0:
        band_mag = np.nan
        band_mag_err = np.nan
    else:
        band_mag = -2.5*np.log10(band_flux) + 8.90007
        band_mag_err = band_flux_err / band_flux * 2.5 / np.log(10.)
    return band_mag, band_mag_err, band_flux, band_flux_err