Source code for common.sky.grb

'''
Created on 12 févr. 2018
@author: Colley Jean-Marc, APC/IN2P3/CNRS
'''

import numpy as np
import scipy.interpolate as spi
import matplotlib.pyplot as plt



# TODO à supprimer ???


class SimuGRB(object):
    def __init__(self):
        self.lc = np.empty((1,2), dtype=np.float32)       
        
        
    def _update_interpol(self):
        self.lc_itp = spi.interp1d(self.lc[:,0], self.lc[:,1])
        
        
    def read_grb_lc_from_swift(self, p_file):
        """
        """
        a = np.loadtxt(p_file, comments="!")
        self.lc = np.hstack((a[:,0], a[:,3]))
        print(self.lc)
        self._update_interpol()

                
    def get_nb_total_photon(self):
        """
        """        
        nb_total = np.trapz(self.lc[:,1], self.lc[:,0])
        return nb_total
    
    
    def simu_time_lc(self, p_ntime):
        # 1-tirage couple (t_i, v_i) uniforme
        # 2- interpol en t_i de la courbe de lumière
        # 3- test de rejet avec where
        # ai-je assez de couple sinon retour à 1
        again = True
        n_target = p_ntime
        max_v = np.max(self.lc[:,1])
        t_lc = None
        while(again):
            print(n_target)            
            ns = n_target*3
            t_i = np.random.uniform(0, self.lc[-1,0], ns)
            v_i = np.random.uniform(0, max_v, ns)
            lc_i = self.lc_itp(t_i)
            idx_ok = np.where(v_i < lc_i)[0]
            if t_lc is None:                
                t_lc = t_i[idx_ok]
            else:
                t_lc = np.append(t_lc, t_i[idx_ok])                
            n_target -= idx_ok.size
            if (n_target <= 0) : 
                again = False            
        return t_lc[0:p_ntime+1]

    
    def plot_lc(self):
        plt.figure()
        plt.plot(self.lc[:,0], self.lc[:,1])
        plt.title("Courbe de lumière")
        
        
    def plot_lci(self):
        plt.figure()
        plt.plot(self.lc[:,0], self.lc[:,1],"*")
        x = np.linspace(self.lc[0,0], self.lc[-1,0])       
        plt.plot(x, self.lc_itp(x))
        plt.title("Courbe de lumière interpol")