Justine Devin / gammalib

Commit 0be65377cdbbd02b1ab9400ce7c9436dccbd4a81

Authored by knodlseder
1 parent 7b177e08

implement node loading

Inline Side-by-side
Showing 1 changed file with 39 additions and 8 deletions
src/model/GModelSpectralFunc.cpp
  View file @0be6537
1 1 /***************************************************************************
2 2 * GModelSpectralFunc.cpp - Spectral function model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2010 by Jurgen Knodlseder *
  4 + * copyright (C) 2010 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software; you can redistribute it and/or modify *
... ... @@ -12,7 +12,7 @@
12 12 ***************************************************************************/
13 13 /**
14 14 * @file GModelSpectralFunc.cpp
15   - * @brief GModelSpectralFunc class implementation.
  15 + * @brief Spectral function model class implementation
16 16 * @author J. Knodlseder
17 17 */
18 18  
... ... @@ -23,12 +23,14 @@
23 23 #include <math.h>
24 24 #include "GException.hpp"
25 25 #include "GTools.hpp"
  26 +#include "GCsv.hpp"
26 27 #include "GModelSpectralFunc.hpp"
27 28  
28 29 /* __ Method name definitions ____________________________________________ */
29 30 #define G_PAR "GModelSpectralFunc::par(int)"
30 31 #define G_READ "GModelSpectralFunc::read(GXmlElement&)"
31 32 #define G_WRITE "GModelSpectralFunc::write(GXmlElement&)"
  33 +#define G_LOAD_NODES "GModelSpectralFunc::load_nodes(std::string&)"
32 34  
33 35 /* __ Macros _____________________________________________________________ */
34 36  
... ... @@ -208,11 +210,13 @@ GModelSpectralFunc* GModelSpectralFunc::clone(void) const
208 210 * \f[I(E)=norm f(E)\f]
209 211 * where
210 212 * \f$norm\f$ is the normalization of the function.
  213 + * Note that the node energies are stored as log10 of energy in units of
  214 + * MeV.
211 215 ***************************************************************************/
212 216 double GModelSpectralFunc::eval(const GEnergy& srcEng)
213 217 {
214 218 // Interpolate function
215   - double func = m_nodes.interpolate(srcEng.MeV(), m_values);
  219 + double func = m_nodes.interpolate(srcEng.log10MeV(), m_values);
216 220  
217 221 // Compute function value
218 222 double value = norm() * func;
... ... @@ -240,8 +244,8 @@ double GModelSpectralFunc::eval(const GEnergy&amp; srcEng)
240 244 ***************************************************************************/
241 245 double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng)
242 246 {
243   - // Interpolate function (dummy so far)
244   - double func = m_nodes.interpolate(srcEng.MeV(), m_values);
  247 + // Interpolate function
  248 + double func = m_nodes.interpolate(srcEng.log10MeV(), m_values);
245 249  
246 250 // Compute function value
247 251 double value = norm() * func;
... ... @@ -363,6 +367,7 @@ std::string GModelSpectralFunc::print(void) const
363 367 // Append header
364 368 result.append("=== GModelSpectralFunc ===\n");
365 369 result.append(parformat("Function file")+m_filename);
  370 + result.append(parformat("Number of nodes")+str(m_nodes.size()));
366 371 result.append(parformat("Number of parameters")+str(size()));
367 372 for (int i = 0; i < size(); ++i) {
368 373 result.append("\n"+parformat("Parameter "+str(i+1)));
... ... @@ -442,10 +447,17 @@ void GModelSpectralFunc::free_members(void)
442 447 *
443 448 * @param[in] filename File name.
444 449 *
445   - * Node energies are stored in units of MeV.
446   - * At least 2 nodes are required.
  450 + * @exception GException::not_enough_data
  451 + * File contains less than 2 nodes
  452 + * @exception GException::not_enough_columns
  453 + * File contains less than 2 columns
447 454 *
448   - * @todo Not yet implemented.
  455 + * The file function is stored as a column separated value table (CSV) in an
  456 + * ASCII file with (at least) 2 columns. The first column specifies the
  457 + * energy in MeV while the second column specifies the intensity at this
  458 + * energy in units of ph/cm2/s/MeV.
  459 + * The node energies will be stored as log10 of the energy given in units
  460 + * of MeV. At least 2 nodes and 2 columns are required.
449 461 ***************************************************************************/
450 462 void GModelSpectralFunc::load_nodes(const std::string& filename)
451 463 {
... ... @@ -456,6 +468,25 @@ void GModelSpectralFunc::load_nodes(const std::string&amp; filename)
456 468 // Set filename
457 469 m_filename = filename;
458 470  
  471 + // Load file
  472 + GCsv csv = GCsv(filename);
  473 +
  474 + // Check if there are at least 2 nodes
  475 + if (csv.nrows() < 2)
  476 + throw GException::not_enough_data(G_LOAD_NODES, filename,
  477 + csv.nrows());
  478 +
  479 + // Check if there are at least 2 columns
  480 + if (csv.ncols() < 2)
  481 + throw GException::not_enough_columns(G_LOAD_NODES, filename,
  482 + csv.ncols());
  483 +
  484 + // Setup nodes
  485 + for (int i = 0; i < csv.nrows(); ++i) {
  486 + m_nodes.append(log10(csv.real(i,0)));
  487 + m_values.push_back(csv.real(i,1));
  488 + }
  489 +
459 490 // Return
460 491 return;
461 492 }
... ...