import numpy as np
import sys
import os
from scipy.interpolate import interp1d
from astropy.io import ascii
from . import utils
from . import constants as cc
from . import photometry as phot
from . import camera as cam
from .preliminary_computations import A_tel,system_response

#######################################

def set_object(info_dict):
    """ Compute the number of electrons coming from the object per second

    Parameters
    ----------
    info_dict: dictionary

    wavelength : array 
                 wavelengths in angstrom

    Returns
    ---------
    F_e_s : float
            countrate of the object in e-/s

    """
    if info_dict['object_type'] == 'magnitude':
         mag =  info_dict['object_magnitude']
         F_e_s = 10**(0.4*(info_dict['zeropoint'] - mag))        # e/s
         fmag = mag*np.ones(len(info_dict['wavelength_ang']))
         fJy = phot.mag2Jy(info_dict, fmag)   # Jy
         flam = utils.fJy_to_flambda(info_dict['wavelength_ang'], fJy)  # erg/s/cm2/A
         fph = utils.flambda_to_fph(info_dict['wavelength_ang'], flam)        # ph/s/cm2/A   

 
    elif info_dict['object_type'] == 'spectrum':
  
        object_path=info_dict['MainDirectory']+info_dict['object_folder']+info_dict['object_file']

        inFile=open(object_path, "r")
        lines=inFile.readlines()

        x=[] # wavelength
        y=[] # flux in erg/s/cm2/A
        for line in lines:
            if line[0] != "#" and len(line) > 3:
                bits=line.split()
                y.append(float(bits[1]))
                x.append(float(bits[0]))

        x=np.array(x)
        y=np.array(y, dtype=np.float64)

        f = interp1d(x,y,kind='linear')
        if min(info_dict['wavelength_ang']) < min(x) or max(info_dict['wavelength_ang']) > max(x):
            print ('The wavelength coverage must be smaller or equal to the one of the input spectrum. Please adapt it in the configuration file.\nSpectrum wavelength coverage: %.2f-%.2f Angstroms\nCurrent chosen wavelength coverage: %.2f-%.2f Angstroms' % (min(x),max(x),min(info_dict['wavelength_ang']),max(info_dict['wavelength_ang'])))
        flam = f(info_dict['wavelength_ang'])      # erg/s/cm2/A
        fph = utils.flambda_to_fph(info_dict['wavelength_ang'],flam)   # ph/s/cm2/A
        fmag = phot.Jy2Mag(info_dict,utils.flambda_to_fJy(info_dict['wavelength_ang'], flam))              
        F_e_s = np.trapz(fph * info_dict['system_response'] * info_dict['Trans_atmosphere'],info_dict['wavelength_ang']) * info_dict['A_tel'] # e/s 


    elif info_dict['object_type'] == 'grb_sim':
        #-----------------------------
        # Compute the emission spectra
        #-----------------------------
         

        if info_dict['grb_model'] != 'LightCurve':
            if info_dict['grb_model'] == 'gs02':
                try:
                    from pyGRBaglow.synchrotron_model import fireball_afterglow as grb
                except ValueError:
                    print ('Package pyGRBaglow not found. Need to be installed')

                td=info_dict['t_sinceBurst'] # in days
                DIT = info_dict['exptime']   # in second
                #print (info_dict['n0'],info_dict['eps_b'],info_dict['eps_e'],info_dict['E_iso'],info_dict['eta'],info_dict['p'],info_dict['Y'],info_dict['grb_redshift'],info_dict['ism_type'])
                afterglow=grb(n0=info_dict['n0'],eps_b=info_dict['eps_b'],eps_e=info_dict['eps_e'],E_iso=info_dict['E_iso'],eta=info_dict['eta'],p=info_dict['p'],Y=info_dict['Y'],z=info_dict['grb_redshift'],ism_type=info_dict['ism_type'],disp=0)
                time_grb=np.linspace(td,td+DIT/86400,5)     # divides exposure time in 5
                frequencies = cc.c_light_m_s / (info_dict['wavelength_ang']*1e-10)

                afterglow_lc=afterglow.light_curve(time_grb,frequencies)   #in mJy
                #afterglow_lc2=afterglow.light_curve(td+DIT/2/86400,frequencies)
                #factor to convert in Jy
                factor_Jy=1e-3
                factor_time=86400

            elif info_dict['grb_model']== 'SPL': 
                try:
                    from pyGRBaglow.template_models import Templates as grb
                except ValueError:
                    print ('Package grb_afterglow not found. Need to be installed')

                td=info_dict['t_sinceBurst']*86400 # in second
                DIT = info_dict['exptime']   # in second

                time_grb=np.linspace(td,td+DIT,5)
                afterglow=grb(F0=info_dict['F0'],t0=info_dict['t0'],wvl0=info_dict['wvl0'])
                afterglow_lc = afterglow.light_curve(info_dict['wavelength_ang'],time_grb,[info_dict['alpha'],info_dict['beta']],model='SPL')  #in Jy
                #conversion factor set to 1 (already in Jy)
                factor_Jy=1
                factor_time=1

            elif info_dict['grb_model']== 'BPL': 
                try:
                    from pyGRBaglow.template_models import Templates as grb
                except ValueError:
                    print ('Package grb_afterglow not found. Need to be installed')

                td=info_dict['t_sinceBurst']*86400 # in second
                DIT = info_dict['exptime']   # in second

                time_grb=np.linspace(td,td+DIT,5)
                afterglow=grb(F0=info_dict['F0'],t0=info_dict['t0'],wvl0=info_dict['wvl0'])
                afterglow_lc = afterglow.light_curve(info_dict['wavelength_ang'],time_grb,[info_dict['alpha1'],info_dict['alpha2'],info_dict['beta'],info_dict['s']],model='BPL')  #in Jy
                #conversion factor set to 1 (already in Jy)
                factor_Jy=1
                factor_time=1

            #Sum over the time
            #dtime=np.diff(time_series)
            sed_stacked=np.zeros(len(info_dict['wavelength_ang']))
            #sed_stacked2=np.zeros(len(info_dict['wavelength_ang']))
            for i in range(len(info_dict['wavelength_ang'])):
                sed_stacked[i]=np.trapz(afterglow_lc[:,i],time_grb)
                #sed_stacked2[i]=afterglow_lc2[:,i]

            #free memory
            afterglow_lc=None

            #Convert to Jy
            grb_fJy = sed_stacked*factor_Jy /(DIT/factor_time) # GRB SED in Jy in the observer frame
            # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula. It is basically the mean value in Jy in the time interval corresponding to the exposure time. This assumption is valid for rather short exposures and probably not for long exposures, it has to be tested.
            #grb_fJy2 = sed_stacked2*1e-3

            #print (grb_fJy,grb_fJy2)

            #Apply Host galaxy and IGM extinction
            try:
                from pyGRBaglow.igm import meiksin,madau
                from pyGRBaglow.reddening import reddening
            except ValueError:
                print ('Package los_extinction not found. Need to be installed')

            if info_dict['IGM_extinction_model'] =='meiksin':
                grb_fJy *= meiksin(info_dict['wavelength_ang']/10,info_dict['grb_redshift'])
            elif info_dict['IGM_extinction_model'] =='madau':
                grb_fJy *= madau(info_dict['wavelength_ang'],info_dict['grb_redshift'])

            if info_dict['host_extinction_law'] in ['mw','lmc','smc']:
                grb_fJy *= reddening(info_dict['wavelength_ang'],info_dict['grb_redshift'],Av=info_dict['Av_Host']).Pei92(law=info_dict['host_extinction_law'])[1]

            if info_dict['galactic_extinction_law'].lower() != 'none':
                grb_fJy *= reddening(info_dict['wavelength_ang'],0,Av=info_dict['Av_galactic']).Pei92(law='mw')[1]


            """
            #Integration over the exposure time for each wavelength because the GRB spectrum varies with the time of observation
            int_grb_sed = []
            def integrand(x,f):
                return afterglow.light_curve(x,f)
            for index in range(len(info_dict['wavelength_ang'])):
                int_grb_sed.append(quad(integrand,td,td+DIT/86400.,args=(frequencies[index]))[0])     # Flux in mJy.day
         
            # 1e-3 to convert from mJy --> Jy. 
            grb_fJy = np.array(int_grb_sed) *1e-3 /(DIT/86400.) # GRB SED in Jy in the GRB restframe
            # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula
         
            grb_fJy *= meiksin(info_dict['wavelength_ang']/10,info_dict['grb_redshift'])
            grb_fJy *= reddening(info_dict['wavelength_ang'],info_dict['grb_redshift'],Av=info_dict['Av_Host']).Li07(3)
            """


        elif info_dict['grb_model'] == 'LightCurve':
            #-----------------------------
            # Loading light curve
            #-----------------------------
            object_path=info_dict['MainDirectory']+info_dict['object_folder']+info_dict['object_file']

            LC_data = ascii.read(object_path)

            #print (LC_data)
            # Fluxes are assumed to be in microJansky

            # get wavelength
            wvl_list=[]
            for dat in LC_data.group_by(['wvl']).groups.keys:
                wvl_list.append(dat[0])

            td=info_dict['t_sinceBurst'] # in second
            DIT = info_dict['exptime']   # in second

            #print (td)

            mask_time = (LC_data['Time'] >= td) & (LC_data['Time'] <= td+DIT)
            time_list=[]
            for dat in LC_data[mask_time].group_by(['Time']).groups.keys:
                time_list.append(dat[0])
            sed_stacked=np.zeros(len(wvl_list))

            #Sum over the time
            for i,wv in enumerate(wvl_list):
                mask_int = LC_data['wvl'][mask_time] == wv    
                sed_stacked[i]=np.trapz(LC_data['flux'][mask_time][mask_int],time_list)

            # Resample the wavelength
            flux_interp = interp1d(wvl_list,sed_stacked)
            sed_stacked_resampled=flux_interp(info_dict['wavelength_ang'])

            #print (sed_stacked_resampled)
            factor_Jy=1e-6
            factor_time=1

            #Convert to Jy
            grb_fJy = sed_stacked_resampled*factor_Jy /(DIT/factor_time) # GRB SED in Jy in the observer frame
            # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula. It is basically the mean value in Jy in the time interval corresponding to the exposure time. This assumption is valid for rather short exposures and probably not for long exposures, it has to be tested.
            #grb_fJy2 = sed_stacked2*1e-3

            #print (grb_fJy,grb_fJy2)



        flam = utils.fJy_to_flambda(info_dict['wavelength_ang'], grb_fJy)   # erg/s/cm2/A
        fph = utils.flambda_to_fph(info_dict['wavelength_ang'], flam)    # ph/s/cm2/A
        F_e_s = np.trapz(fph * info_dict['system_response'] * info_dict['Trans_atmosphere'],info_dict['wavelength_ang']) * info_dict['A_tel'] # e/s



    info_dict['Object_fph']=fph
    info_dict['Object_fes']=F_e_s
    info_dict['Object_mag']=phot.Jy2Mag(info_dict,utils.flambda_to_fJy(info_dict['wavelength_ang'], flam)) 
    return info_dict

import numpy as np
import sys
import os
from scipy.interpolate import interp1d
from astropy.io import ascii
from . import utils
from . import constants as cc
from . import photometry as phot
from . import camera as cam
from .preliminary_computations import A_tel,system_response

#######################################

def set_object(info_dict):
    """ Compute the number of electrons coming from the object per second

    Parameters
    ----------
    info_dict: dictionary

    wavelength : array 
                 wavelengths in angstrom

    Returns
    ---------
    F_e_s : float
            countrate of the object in e-/s

    """
    if info_dict['object_type'] == 'magnitude':
         mag =  info_dict['object_magnitude']
         F_e_s = 10**(0.4*(info_dict['zeropoint'] - mag))        # e/s
         fmag = mag*np.ones(len(info_dict['wavelength_ang']))
         fJy = phot.mag2Jy(info_dict, fmag)   # Jy
         flam = utils.fJy_to_flambda(info_dict['wavelength_ang'], fJy)  # erg/s/cm2/A
         fph = utils.flambda_to_fph(info_dict['wavelength_ang'], flam)        # ph/s/cm2/A   

 
    elif info_dict['object_type'] == 'spectrum':
  
        object_path=info_dict['MainDirectory']+info_dict['object_folder']+info_dict['object_file']

        inFile=open(object_path, "r")
        lines=inFile.readlines()

        x=[] # wavelength
        y=[] # flux in erg/s/cm2/A
        for line in lines:
            if line[0] != "#" and len(line) > 3:
                bits=line.split()
                y.append(float(bits[1]))
                x.append(float(bits[0]))

        x=np.array(x)
        y=np.array(y, dtype=np.float64)

        f = interp1d(x,y,kind='linear')
        if min(info_dict['wavelength_ang']) < min(x) or max(info_dict['wavelength_ang']) > max(x):
            print ('The wavelength coverage must be smaller or equal to the one of the input spectrum. Please adapt it in the configuration file.\nSpectrum wavelength coverage: %.2f-%.2f Angstroms\nCurrent chosen wavelength coverage: %.2f-%.2f Angstroms' % (min(x),max(x),min(info_dict['wavelength_ang']),max(info_dict['wavelength_ang'])))
        flam = f(info_dict['wavelength_ang'])      # erg/s/cm2/A
        fph = utils.flambda_to_fph(info_dict['wavelength_ang'],flam)   # ph/s/cm2/A
        fmag = phot.Jy2Mag(info_dict,utils.flambda_to_fJy(info_dict['wavelength_ang'], flam))              
        F_e_s = np.trapz(fph * info_dict['system_response'] * info_dict['Trans_atmosphere'],info_dict['wavelength_ang']) * info_dict['A_tel'] # e/s 


    elif info_dict['object_type'] == 'grb_sim':
        #-----------------------------
        # Compute the emission spectra
        #-----------------------------
         

        if info_dict['grb_model'] != 'LightCurve':
            if info_dict['grb_model'] == 'gs02':
                try:
                    from pyGRBaglow.synchrotron_model import fireball_afterglow as grb
                except ValueError:
                    print ('Package pyGRBaglow not found. Need to be installed')

                td=info_dict['t_sinceBurst'] # in days
                DIT = info_dict['exptime']   # in second
                #print (info_dict['n0'],info_dict['eps_b'],info_dict['eps_e'],info_dict['E_iso'],info_dict['eta'],info_dict['p'],info_dict['Y'],info_dict['grb_redshift'],info_dict['ism_type'])
                afterglow=grb(n0=info_dict['n0'],eps_b=info_dict['eps_b'],eps_e=info_dict['eps_e'],E_iso=info_dict['E_iso'],eta=info_dict['eta'],p=info_dict['p'],Y=info_dict['Y'],z=info_dict['grb_redshift'],ism_type=info_dict['ism_type'],disp=0)
                time_grb=np.linspace(td,td+DIT/86400,5)     # divides exposure time in 5
                frequencies = cc.c_light_m_s / (info_dict['wavelength_ang']*1e-10)

                afterglow_lc=afterglow.light_curve(time_grb,frequencies)   #in mJy
                #afterglow_lc2=afterglow.light_curve(td+DIT/2/86400,frequencies)
                #factor to convert in Jy
                factor_Jy=1e-3
                factor_time=86400

            elif info_dict['grb_model']== 'SPL': 
                try:
                    from pyGRBaglow.template_models import Templates as grb
                except ValueError:
                    print ('Package grb_afterglow not found. Need to be installed')

                td=info_dict['t_sinceBurst']*86400 # in second
                DIT = info_dict['exptime']   # in second

                time_grb=np.linspace(td,td+DIT,5)
                afterglow=grb(F0=info_dict['F0'],t0=info_dict['t0'],wvl0=info_dict['wvl0'])
                afterglow_lc = afterglow.light_curve(info_dict['wavelength_ang'],time_grb,[info_dict['alpha'],info_dict['beta']],model='SPL')  #in Jy
                #conversion factor set to 1 (already in Jy)
                factor_Jy=1
                factor_time=1

            elif info_dict['grb_model']== 'BPL': 
                try:
                    from pyGRBaglow.template_models import Templates as grb
                except ValueError:
                    print ('Package grb_afterglow not found. Need to be installed')

                td=info_dict['t_sinceBurst']*86400 # in second
                DIT = info_dict['exptime']   # in second

                time_grb=np.linspace(td,td+DIT,5)
                afterglow=grb(F0=info_dict['F0'],t0=info_dict['t0'],wvl0=info_dict['wvl0'])
                afterglow_lc = afterglow.light_curve(info_dict['wavelength_ang'],time_grb,[info_dict['alpha1'],info_dict['alpha2'],info_dict['beta'],info_dict['s']],model='BPL')  #in Jy
                #conversion factor set to 1 (already in Jy)
                factor_Jy=1
                factor_time=1

            #Sum over the time
            #dtime=np.diff(time_series)
            sed_stacked=np.zeros(len(info_dict['wavelength_ang']))
            #sed_stacked2=np.zeros(len(info_dict['wavelength_ang']))
            for i in range(len(info_dict['wavelength_ang'])):
                sed_stacked[i]=np.trapz(afterglow_lc[:,i],time_grb)
                #sed_stacked2[i]=afterglow_lc2[:,i]

            #free memory
            afterglow_lc=None

            #Convert to Jy
            grb_fJy = sed_stacked*factor_Jy /(DIT/factor_time) # GRB SED in Jy in the observer frame
            # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula. It is basically the mean value in Jy in the time interval corresponding to the exposure time. This assumption is valid for rather short exposures and probably not for long exposures, it has to be tested.
            #grb_fJy2 = sed_stacked2*1e-3

            #print (grb_fJy,grb_fJy2)

            #Apply Host galaxy and IGM extinction
            try:
                from pyGRBaglow.igm import meiksin,madau
                from pyGRBaglow.reddening import reddening
            except ValueError:
                print ('Package los_extinction not found. Need to be installed')

            if info_dict['IGM_extinction_model'] =='meiksin':
                grb_fJy *= meiksin(info_dict['wavelength_ang']/10,info_dict['grb_redshift'])
            elif info_dict['IGM_extinction_model'] =='madau':
                grb_fJy *= madau(info_dict['wavelength_ang'],info_dict['grb_redshift'])

            if info_dict['host_extinction_law'] in ['mw','lmc','smc']:
                grb_fJy *= reddening(info_dict['wavelength_ang'],info_dict['grb_redshift'],Av=info_dict['Av_Host']).Pei92(law=info_dict['host_extinction_law'])[1]

            if info_dict['galactic_extinction_law'].lower() != 'none':
                grb_fJy *= reddening(info_dict['wavelength_ang'],0,Av=info_dict['Av_galactic']).Pei92(law='mw')[1]


            """
            #Integration over the exposure time for each wavelength because the GRB spectrum varies with the time of observation
            int_grb_sed = []
            def integrand(x,f):
                return afterglow.light_curve(x,f)
            for index in range(len(info_dict['wavelength_ang'])):
                int_grb_sed.append(quad(integrand,td,td+DIT/86400.,args=(frequencies[index]))[0])     # Flux in mJy.day
         
            # 1e-3 to convert from mJy --> Jy. 
            grb_fJy = np.array(int_grb_sed) *1e-3 /(DIT/86400.) # GRB SED in Jy in the GRB restframe
            # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula
         
            grb_fJy *= meiksin(info_dict['wavelength_ang']/10,info_dict['grb_redshift'])
            grb_fJy *= reddening(info_dict['wavelength_ang'],info_dict['grb_redshift'],Av=info_dict['Av_Host']).Li07(3)
            """


        elif info_dict['grb_model'] == 'LightCurve':
            #-----------------------------
            # Loading light curve
            #-----------------------------
            object_path=info_dict['MainDirectory']+info_dict['object_folder']+info_dict['object_file']

            LC_data = ascii.read(object_path)

            #print (LC_data)
            # Fluxes are assumed to be in microJansky

            # get wavelength
            wvl_list=[]
            for dat in LC_data.group_by(['wvl']).groups.keys:
                wvl_list.append(dat[0])

            td=info_dict['t_sinceBurst'] # in second
            DIT = info_dict['exptime']   # in second

            #print (td)

            mask_time = (LC_data['Time'] >= td) & (LC_data['Time'] <= td+DIT)
            time_list=[]
            for dat in LC_data[mask_time].group_by(['Time']).groups.keys:
                time_list.append(dat[0])
            sed_stacked=np.zeros(len(wvl_list))

            #Sum over the time
            for i,wv in enumerate(wvl_list):
                mask_int = LC_data['wvl'][mask_time] == wv    
                sed_stacked[i]=np.trapz(LC_data['flux'][mask_time][mask_int],time_list)

            # Resample the wavelength
            flux_interp = interp1d(wvl_list,sed_stacked)
            sed_stacked_resampled=flux_interp(info_dict['wavelength_ang'])

            #print (sed_stacked_resampled)
            factor_Jy=1e-6
            factor_time=1

            #Convert to Jy
            grb_fJy = sed_stacked_resampled*factor_Jy /(DIT/factor_time) # GRB SED in Jy in the observer frame
            # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula. It is basically the mean value in Jy in the time interval corresponding to the exposure time. This assumption is valid for rather short exposures and probably not for long exposures, it has to be tested.
            #grb_fJy2 = sed_stacked2*1e-3

            #print (grb_fJy,grb_fJy2)



        flam = utils.fJy_to_flambda(info_dict['wavelength_ang'], grb_fJy)   # erg/s/cm2/A
        fph = utils.flambda_to_fph(info_dict['wavelength_ang'], flam)    # ph/s/cm2/A
        F_e_s = np.trapz(fph * info_dict['system_response'] * info_dict['Trans_atmosphere'],info_dict['wavelength_ang']) * info_dict['A_tel'] # e/s



    info_dict['Object_fph']=fph
    info_dict['Object_fes']=F_e_s
    info_dict['Object_mag']=phot.Jy2Mag(info_dict,utils.flambda_to_fJy(info_dict['wavelength_ang'], flam)) 
    return info_dict

####################################### Compute the number of electrons coming from the object per second      Parameters     ----------     info_dict: dictionary      wavelength : array                   wavelengths in angstrom      Returns     ---------     F_e_s : float             countrate of the object in e-/s # e/s # Jy # erg/s/cm2/A # ph/s/cm2/A # wavelength # flux in erg/s/cm2/A # erg/s/cm2/A # ph/s/cm2/A # e/s #----------------------------- # Compute the emission spectra #----------------------------- # in days # in second #print (info_dict['n0'],info_dict['eps_b'],info_dict['eps_e'],info_dict['E_iso'],info_dict['eta'],info_dict['p'],info_dict['Y'],info_dict['grb_redshift'],info_dict['ism_type']) # divides exposure time in 5 #in mJy #afterglow_lc2=afterglow.light_curve(td+DIT/2/86400,frequencies) #factor to convert in Jy # in second # in second #in Jy #conversion factor set to 1 (already in Jy) # in second # in second #in Jy #conversion factor set to 1 (already in Jy) #Sum over the time #dtime=np.diff(time_series) #sed_stacked2=np.zeros(len(info_dict['wavelength_ang'])) #sed_stacked2[i]=afterglow_lc2[:,i] #free memory #Convert to Jy # GRB SED in Jy in the observer frame # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula. It is basically the mean value in Jy in the time interval corresponding to the exposure time. This assumption is valid for rather short exposures and probably not for long exposures, it has to be tested. #grb_fJy2 = sed_stacked2*1e-3 #print (grb_fJy,grb_fJy2) #Apply Host galaxy and IGM extinction #Integration over the exposure time for each wavelength because the GRB spectrum varies with the time of observation             int_grb_sed = []             def integrand(x,f):                 return afterglow.light_curve(x,f)             for index in range(len(info_dict['wavelength_ang'])):                 int_grb_sed.append(quad(integrand,td,td+DIT/86400.,args=(frequencies[index]))[0])     # Flux in mJy.day                       # 1e-3 to convert from mJy --> Jy.              grb_fJy = np.array(int_grb_sed) *1e-3 /(DIT/86400.) # GRB SED in Jy in the GRB restframe             # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula                       grb_fJy *= meiksin(info_dict['wavelength_ang']/10,info_dict['grb_redshift'])             grb_fJy *= reddening(info_dict['wavelength_ang'],info_dict['grb_redshift'],Av=info_dict['Av_Host']).Li07(3) #----------------------------- # Loading light curve #----------------------------- #print (LC_data) # Fluxes are assumed to be in microJansky # get wavelength # in second # in second #print (td) #Sum over the time # Resample the wavelength #print (sed_stacked_resampled) #Convert to Jy # GRB SED in Jy in the observer frame # /(DIT/86400.) to go from Jy.day --> Jy because of the integration in time due to the time dependancy. it has to be divided by the exposure time to recover an e-/s unit for the SNR formula. It is basically the mean value in Jy in the time interval corresponding to the exposure time. This assumption is valid for rather short exposures and probably not for long exposures, it has to be tested. #grb_fJy2 = sed_stacked2*1e-3 #print (grb_fJy,grb_fJy2) # erg/s/cm2/A # ph/s/cm2/A # e/s