GammaLib

#! /usr/bin/env python

# ==========================================================================

# Create and test table model

# -------------------------------------------------------------------------

#

# Copyright (C) 2019 Juergen Knoedlseder

#

# This program is free software: you can redistribute it and/or modify

# it under the terms of the GNU General Public License as published by

# the Free Software Foundation, either version 3 of the License, or

# (at your option) any later version.

#

# This program is distributed in the hope that it will be useful,

# but WITHOUT ANY WARRANTY; without even the implied warranty of

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

# GNU General Public License for more details.

#

# You should have received a copy of the GNU General Public License

# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#

# ==========================================================================

import math

import gammalib

# ================== #

# Create table model #

# ================== #

def create_table_model(nindex=100, ncutoff=50, nebins=50):

    """

    Create an exponentially cutoff power law table model

    """

    # Set energy boundaries

    emin  = gammalib.GEnergy(10.0, 'GeV')

    emax  = gammalib.GEnergy(300.0, 'TeV')

    ebins = gammalib.GEbounds(nebins, emin, emax)

    # Set table model parameters

    val1 = [-3.0+i*0.02 for i in range(nindex)]

    val2 = [math.pow(10.0,-1.0+i*0.05+6.0) for i in range(ncutoff)]

    par1 = gammalib.GModelSpectralTablePar(gammalib.GModelPar('Index',-2.0), val1)

    par2 = gammalib.GModelSpectralTablePar(gammalib.GModelPar('Cutoff',1.0e6), val2)

    pars = gammalib.GModelSpectralTablePars()

    pars.append(par1)

    pars.append(par2)

    # Set table model spectra

    spectra = gammalib.GNdarray(nindex, ncutoff, nebins)

    pivot   = gammalib.GEnergy(0.3, 'TeV')

    for ii, index in enumerate(val1):

        for ic, cutoff in enumerate(val2):

            #print(index,cutoff)

            ecut = gammalib.GEnergy(cutoff, 'MeV')

            spec = gammalib.GModelSpectralExpPlaw(1.0, index, pivot, ecut)

            for ie in range(nebins):

                spectra[ii, ic, ie] = spec.eval(ebins.elogmean(ie))

    # Create table model

    model = gammalib.GModelSpectralTable(ebins, pars, spectra)

    # Set interpolation method

    model.table_par('Index').method(0)   # Linear for index

    #model.table_par('Cutoff').method(1)  # Logarithmic for cut-off

    model.table_par('Cutoff').method(0)  # Linear for cut-off

    # Return model

    return model

# ========== #

# Check flux #

# ========== #

def check_flux(model):

    """

    Create an exponentially cutoff power law table model

    """

    # Exponential cutoff

    pivot = gammalib.GEnergy(0.3, 'TeV')

    ecut  = gammalib.GEnergy(1.0, 'TeV')

    spec  = gammalib.GModelSpectralExpPlaw(1.0, -2.0, pivot, ecut)

    emin  = gammalib.GEnergy(100, 'GeV')

    emax  = gammalib.GEnergy(10, 'TeV')

    # Print flux

    print('GModelSpectralExpPlaw.flux():  %f' % spec.flux(emin, emax))

    print('GModelSpectralTable.flux():    %f' % model.flux(emin, emax))

    print('GModelSpectralExpPlaw.eflux(): %f' % spec.eflux(emin, emax))

    print('GModelSpectralTable.eflux():   %f' % model.eflux(emin, emax))

    # Return

    return

# ======================== #

# Main routine entry point #

# ======================== #

if __name__ == '__main__':

    # Write header

    print('*******************************')

    print('* Create and test table model *')

    print('*******************************')

    # Create table model

    model = create_table_model()

    # Save table model

    model.save('model_table.fits', True)

    # Print table model

    print(model)

    # Check flux

    check_flux(model)

    # Load table model

    model.load('model_table.fits')

    # Re-print table model

    print(model)

    # Re-save table model

    model.save('model_table2.fits', True)

    # Update

    #model.update()

    #model.update()