test.py

Knödlseder Jürgen, 06/23/2020 04:18 PM

Download (3.93 KB)

 
1
#! /usr/bin/env python
2
# ==========================================================================
3
# Create and test table model
4
# -------------------------------------------------------------------------
5
#
6
# Copyright (C) 2019 Juergen Knoedlseder
7
#
8
# This program is free software: you can redistribute it and/or modify
9
# it under the terms of the GNU General Public License as published by
10
# the Free Software Foundation, either version 3 of the License, or
11
# (at your option) any later version.
12
#
13
# This program is distributed in the hope that it will be useful,
14
# but WITHOUT ANY WARRANTY; without even the implied warranty of
15
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16
# GNU General Public License for more details.
17
#
18
# You should have received a copy of the GNU General Public License
19
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
20
#
21
# ==========================================================================
22
import math
23
import gammalib
24

    
25

    
26
# ================== #
27
# Create table model #
28
# ================== #
29
def create_table_model(nindex=100, ncutoff=50, nebins=50):
30
    """
31
    Create an exponentially cutoff power law table model
32
    """
33
    # Set energy boundaries
34
    emin  = gammalib.GEnergy(10.0, 'GeV')
35
    emax  = gammalib.GEnergy(300.0, 'TeV')
36
    ebins = gammalib.GEbounds(nebins, emin, emax)
37

    
38
    # Set table model parameters
39
    val1 = [-3.0+i*0.02 for i in range(nindex)]
40
    val2 = [math.pow(10.0,-1.0+i*0.05+6.0) for i in range(ncutoff)]
41
    par1 = gammalib.GModelSpectralTablePar(gammalib.GModelPar('Index',-2.0), val1)
42
    par2 = gammalib.GModelSpectralTablePar(gammalib.GModelPar('Cutoff',1.0e6), val2)
43
    pars = gammalib.GModelSpectralTablePars()
44
    pars.append(par1)
45
    pars.append(par2)
46

    
47
    # Set table model spectra
48
    spectra = gammalib.GNdarray(nindex, ncutoff, nebins)
49
    pivot   = gammalib.GEnergy(0.3, 'TeV')
50
    for ii, index in enumerate(val1):
51
        for ic, cutoff in enumerate(val2):
52
            #print(index,cutoff)
53
            ecut = gammalib.GEnergy(cutoff, 'MeV')
54
            spec = gammalib.GModelSpectralExpPlaw(1.0, index, pivot, ecut)
55
            for ie in range(nebins):
56
                spectra[ii, ic, ie] = spec.eval(ebins.elogmean(ie))
57

    
58
    # Create table model
59
    model = gammalib.GModelSpectralTable(ebins, pars, spectra)
60

    
61
    # Set interpolation method
62
    model.table_par('Index').method(0)   # Linear for index
63
    #model.table_par('Cutoff').method(1)  # Logarithmic for cut-off
64
    model.table_par('Cutoff').method(0)  # Linear for cut-off
65

    
66
    # Return model
67
    return model
68

    
69

    
70
# ========== #
71
# Check flux #
72
# ========== #
73
def check_flux(model):
74
    """
75
    Create an exponentially cutoff power law table model
76
    """
77
    # Exponential cutoff
78
    pivot = gammalib.GEnergy(0.3, 'TeV')
79
    ecut  = gammalib.GEnergy(1.0, 'TeV')
80
    spec  = gammalib.GModelSpectralExpPlaw(1.0, -2.0, pivot, ecut)
81
    emin  = gammalib.GEnergy(100, 'GeV')
82
    emax  = gammalib.GEnergy(10, 'TeV')
83

    
84
    # Print flux
85
    print('GModelSpectralExpPlaw.flux():  %f' % spec.flux(emin, emax))
86
    print('GModelSpectralTable.flux():    %f' % model.flux(emin, emax))
87
    print('GModelSpectralExpPlaw.eflux(): %f' % spec.eflux(emin, emax))
88
    print('GModelSpectralTable.eflux():   %f' % model.eflux(emin, emax))
89

    
90
    # Return
91
    return
92

    
93

    
94
# ======================== #
95
# Main routine entry point #
96
# ======================== #
97
if __name__ == '__main__':
98

    
99
    # Write header
100
    print('*******************************')
101
    print('* Create and test table model *')
102
    print('*******************************')
103

    
104
    # Create table model
105
    model = create_table_model()
106

    
107
    # Save table model
108
    model.save('model_table.fits', True)
109

    
110
    # Print table model
111
    print(model)
112

    
113
    # Check flux
114
    check_flux(model)
115

    
116
    # Load table model
117
    model.load('model_table.fits')
118

    
119
    # Re-print table model
120
    print(model)
121

    
122
    # Re-save table model
123
    model.save('model_table2.fits', True)
124

    
125
    # Update
126
    #model.update()
127
    #model.update()
128