notebook.community

Edit and run 





In [1]:



    


import scipy as sp
from matplotlib import pyplot as plt
%pylab inline
import yaml



    


















    






Populating the interactive namespace from numpy and matplotlib

















In [2]:



    


document='''
- crab_pulsar:
    spectrum : 
        powerlaw:
             K : 1
             gamma : 2
    spatial  : 
        point :
            ra : 83.63
            dec : 22.01
    polar :
        angle : 190
        degree : 17.8
- crab_nebula:
'''
print yaml.load(document)



    


















    






[{'crab_pulsar': {'polar': {'angle': 190, 'degree': 17.8}, 'spectrum': {'powerlaw': {'K': 1, 'gamma': 2}}, 'spatial': {'point': {'dec': 22.01, 'ra': 83.63}}}}, {'crab_nebula': None}]

















In [3]:



    


documet='''
- a
- b
'''
print yaml.load(documet)



    


















    






['a', 'b']

















In [78]:



    


from ximpol.srcmodel import xGenerator,xSpectralComponent,xSource



    











In [5]:



    


class xLightCurve(object):
    def __init__(self):
        self.f = None
        pass
    def Steady(self,rate):  
        self.f = rate
    def FunctionalLightCurve(self,function):
        self.f = function
        pass
    def TemplateLightCurve(self,times,rates):
        self.f = interpolate.UnivariateSpline(times,rates,k=1,s=0)
        pass
    def Phase(self,rate,period):
        pass
    def __call__(self,x):
        if isinstance(self.f,float): return sp.ones_like(x)*self.f
        return self.f(x)



    











In [12]:



    






    


















    






Component (pl1+pl2)*(co1) added










    
























In [75]:



    


tmin=0
tmax=100
emin=1
emax=5

aeff=xAeff()
aeff.load_arf('../irf/fits/xipe_baseline.arf')


C=lambda t: 10.0*(1.0+sp.cos(t))
gamma=lambda t: -2+0.01*t

mySource=xSource('test',resolve_name=False)
mySource.setRADec(83.63,22.01)
print mySource.getLB()
spectrum=xSpectralComponent('spectrum')

times=sp.linspace(tmin,tmax,100)
flux=[]
events=[]
for t in times:
    spectrum.powerlaw(C(t),gamma(t))
    x,y = aeff.convolve(spectrum)
    f   = interpolate.UnivariateSpline(x,y,k=1,s=0)
    flux.append(f.integral(emin,emax))
    pass

lc   = interpolate.UnivariateSpline(times,flux,k=1,s=0)

S=xSimulator(lc,lc.integral)
S.setMinMax(tmin,tmax)
events_times=S.generate()
events_energies=[]
for evt in events_times:
    spectrum.powerlaw(C(evt),gamma(evt))
    x,y = aeff.convolve(spectrum)
    f   = interpolate.UnivariateSpline(x,y,k=1,s=0)
    S=xSimulator(f,f.integral)
    S.setMinMax(emin,emax)
    events_energies.append(S.generate(1)[0])
    pass



plt.plot(times,flux)


#spectrum.plot(sp.linspace(1,10,100))
#plt.xscale('log')
#plt.yscale('log')
#spectrum.polarization(0.1,89)
#mySource.addComponent(spectrum)


#temporalProfile=xLightCurve()
#temporalProfile.Steady(100.0)




#plt.plot(x,Crab.components[0](x))

#plt.plot(x,y)
#plt.plot(x,f(x))
#S=xSimulator(f,f.integral)
#S.setMinMax(1,10)

#S.generate()



    


















    






(184.55973191463613, -5.789183329711647)










    
Out[75]:








[<matplotlib.lines.Line2D at 0x10cae01d0>]










    
























In [76]:



    


plt.plot(events_times,events_energies,'o')
print len(events_energies)



    


















    






171901

1 case of a simple constant function:

lc=xLightCurve() lc.Steady(100.0) S=xSimulator(lc) S.setMinMax(0,100) S.generate() S.plot(100)





In [90]:



    


#1 case of a template (points in x and y):
tmin=0
tmax=10
x=sp.linspace(tmin,tmax,100)
y=1000.*sp.ones_like(x)
lc=xLightCurve()
lc.TemplateLightCurve(x,y)
S=xSimulator(lc)
S.SetMinMax(tmin,tmax)
S.generate()
S.plot(100)



    


















    






Done










    
























In [16]:



    


aeff=xAeff()
aeff.load_arf('../irf/fits/xipe_baseline.arf')


x,y=aeff.convolve(Crab.components[0])
f = interpolate.UnivariateSpline(x,y,k=1,s=0)
#plt.plot(x,Crab.components[0](x))

#plt.plot(x,y)
#plt.plot(x,f(x))
S=xSimulator(f,f.integral)
S.setMinMax(1,10)

S.generate()
S.plot()
plt.xscale('log')
plt.yscale('log')
S.nEvents



    


















    






Done










    
Out[16]:








342










    
























In [210]:



    


S.plot()
S.events



    


















    
Out[210]:








array([ 1.00548119,  1.00897353,  1.01402756, ...,  9.84914435,
        9.86744118,  9.88262894])










    
























In [133]:



    


def dNdE(energy): 
    return 1.0*sp.power(energy,-4.8)

x,y=aeff.convolve(dNdE)
f = interpolate.UnivariateSpline(x,y,k=1,s=0)


#plt.plot(x,y)
#plt.plot(x,f(x))
S=xSimulator(f,f.integral)
S.SetMinMax(0.1,100)
S.generate()
S.plot()
plt.xscale('log')
plt.yscale('log')



    


















    






Done










    
























In [20]:



    


import argparse
parser = argparse.ArgumentParser(description='ximpol simulator for X-ray polarimeter')
parser.add_argument('-c',, nargs='+', type=int, help='an integer to be summed')
    parser.add_argument(
        '--log', default=sys.stdout, type=argparse.FileType('w'),
        help='the file where the sum should be written')
    args = parser.parse_args()
    args.log.write('%s' % sum(args.integers))
    args.log.close()



    











In [83]:



    


sp.random.normal??



    











In [101]:



    


class xPsf():
    def __init__(self):
        self.psf=sp.random.normal
        pass
    def __call__(self,energy=0,theta=0):        
        return self.psf(0.0,15.0)



    











In [102]:



    


psf=xPsf()



    











In [105]:



    


psf()



    


















    
Out[105]:








-7.44764831069462

















In [106]:



    


plt.hist2d??



    











In [40]:



    


import scipy as sp
from astropy.io import fits
import numpy as np
class event():
    def __init__(self):
        self.energy=0
        self.time=0
        self.ra=0
        self.dec=0
        self.x=0
        self.y=0
        self.l=0
        self.b=0
        self.angle=0
    def __str__(self):
        return 'Time=%10.3f Energy = %10.1f, Ra=%10.3f Dec=%10.3f Angle=%10.3f' %(self.time,self.energy,self.ra,self.dec,self.angle)
    
class eventList():
    def __init__(self):
        self.energy_array = sp.array([],dtype='double')
        self.time_array   = sp.array([],dtype='double')
        self.ra_array     = sp.array([],dtype='double')
        self.dec_array    = sp.array([],dtype='double')
        self.x_array      = sp.array([],dtype='int')
        self.y_array      = sp.array([],dtype='int')
        self.l_array      = sp.array([],dtype='double')
        self.b_array      = sp.array([],dtype='double')
        self.angle_array  = sp.array([],dtype='double')
        pass
    def fill(self,event):
        self.energy_array=sp.append(self.energy_array,event.energy)
        self.time_array=sp.append(self.time_array,event.time)
        self.ra_array=sp.append(self.ra_array,event.ra)
        self.dec_array=sp.append(self.dec_array,event.dec)
        self.x_array=sp.append(self.x_array,event.x)
        self.y_array=sp.append(self.y_array,event.y)
        self.l_array=sp.append(self.l_array,event.l)
        self.b_array=sp.append(self.b_array,event.b)
        self.angle_array=sp.append(self.angle_array,event.angle)
        pass
    def write_fits(self,file_path):       
        col_time = fits.Column(name='time', format='E', array=self.time_array)
        col_energy = fits.Column(name='energy', format='E', array=self.energy_array)
        col_ra = fits.Column(name='ra', format='E', array=self.ra_array)
        col_dec = fits.Column(name='dec', format='E', array=self.dec_array)
        col_l = fits.Column(name='l', format='E', array=self.l_array)
        col_b = fits.Column(name='b', format='E', array=self.b_array)
        col_x = fits.Column(name='x', format='D', array=self.x_array)
        col_y = fits.Column(name='y', format='D', array=self.y_array)
        col_angle = fits.Column(name='angle', format='E', array=self.angle_array)
        cols = fits.ColDefs([col_time,col_energy,col_ra,col_dec,col_l,col_b,col_x,col_y,col_angle])
        tbhdu = fits.BinTableHDU.from_columns(cols)
        tbhdu.writeto(file_path,clobber=True)



    











In [41]:



    


#import sp.random
import scipy as sp
nevents=100
el=eventList()
for i in range(nevents):
    e=event()
    e.time=1.0*i
    e.energy=sp.random.uniform()
    el.fill(e)
    print e
el.write_fits('fava.fits')



    


















    






Time=     0.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     1.000 Energy =        1.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     2.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     3.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     4.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     5.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     6.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     7.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     8.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=     9.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    10.000 Energy =        1.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    11.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    12.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    13.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    14.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    15.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    16.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    17.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    18.000 Energy =        0.6, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    19.000 Energy =        0.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    20.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    21.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    22.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    23.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    24.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    25.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    26.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    27.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    28.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    29.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    30.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    31.000 Energy =        0.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    32.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    33.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    34.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    35.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    36.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    37.000 Energy =        1.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    38.000 Energy =        0.6, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    39.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    40.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    41.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    42.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    43.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    44.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    45.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    46.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    47.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    48.000 Energy =        0.6, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    49.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    50.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    51.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    52.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    53.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    54.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    55.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    56.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    57.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    58.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    59.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    60.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    61.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    62.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    63.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    64.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    65.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    66.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    67.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    68.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    69.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    70.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    71.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    72.000 Energy =        0.6, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    73.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    74.000 Energy =        1.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    75.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    76.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    77.000 Energy =        0.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    78.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    79.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    80.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    81.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    82.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    83.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    84.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    85.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    86.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    87.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    88.000 Energy =        0.5, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    89.000 Energy =        0.6, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    90.000 Energy =        0.3, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    91.000 Energy =        0.8, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    92.000 Energy =        0.7, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    93.000 Energy =        0.1, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    94.000 Energy =        0.0, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    95.000 Energy =        0.6, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    96.000 Energy =        0.6, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    97.000 Energy =        0.2, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    98.000 Energy =        0.4, Ra=     0.000 Dec=     0.000 Angle=     0.000
Time=    99.000 Energy =        0.9, Ra=     0.000 Dec=     0.000 Angle=     0.000

















In [24]:



    


print b[0]



    


















    






Time=    10.000 Energy =        0.0, Ra=     0.000 Dec=     0.000 Angle=     0.000

















In [35]:



    


sp.array([],dtype='double')



    


















    
Out[35]:








array([], dtype=float64)

















In [27]:



    


a



    


















    
Out[27]:








array([], dtype=float64)

















In [ ]:

Content source: lucabaldini/ximpol

Similar notebooks:

notebook.community | gallery | about