Program Listing for File SED.h

↰ Return to documentation for file (src/lib/SED.h)

/*
  17/11/2014
  Class to store one SED
*/

// avoid multiple def of the same class
#ifndef SED_H  // check that this keyword has been set already
#define SED_H  // define the keyword to be checked

#include <string>
#include <vector>

#include "ext.h"
#include "flt.h"
#include "globals.h"
#include "oneElLambda.h"
#include "onesource.h"
#include "opa.h"

enum object_type {
  GAL,
  QSO,
  STAR
};

class SED {
 protected:
  int idAge;

 public:
  vector<oneElLambda> lamb_flux;
  vector<double> kcorr,
      mag_z0;
  vector<double> mag;
  string name;
  bool has_emlines;

  object_type nlib;

  int nummod,
      index,
      index_z0;

  double red,
      chi2 = HIGH_CHI2,
      dm;

  double
      luv,
      lopt,
      lnir,
      ltir;

  double mass,
      age,
      sfr,
      ssfr;

  double ebv,
      mag0,
      distMod;

  int extlawId;

  double qi[4];

  vector<oneElLambda> fac_line;

  SED(const string name, int nummod = 0, string type = "G");
  SED(const string nameC, double tauC, double ageC, int nummodC, string typeC,
      int idAgeC);
  SED(SED const &p) {
    idAge = p.idAge;
    lamb_flux = p.lamb_flux;
    kcorr = p.kcorr;
    mag = p.mag;
    name = p.name;
    has_emlines = p.has_emlines;
    nummod = p.nummod;
    nlib = p.nlib;
    index = p.index;
    index_z0 = p.index_z0;
    red = p.red;
    chi2 = p.chi2;
    dm = p.dm;
    mass = p.mass;
    age = p.age;
    sfr = p.sfr;
    ltir = p.ltir;
    ebv = p.ebv;
    extlawId = p.extlawId;
    fac_line = p.fac_line;
    distMod = p.distMod;
  };


  inline static object_type string_to_object(const string &type) {
    char t = toupper(type[0]);
    if (t == 'S') {
      return STAR;
    } else if (t == 'Q') {
      return QSO;
    } else if (t == 'G') {
      return GAL;
    } else {
      throw invalid_argument("Object type not recognized: " + type);
    }
  }

  object_type get_object_type() const { return nlib; }

  bool is_gal() { return nlib == GAL ? true : false; }

  bool is_qso() { return nlib == QSO ? true : false; }

  bool is_star() { return nlib == STAR ? true : false; }

  virtual ~SED();
  void read(const string &sedFile);
  void warning_integrateSED(const vector<flt> &filters, bool verbose = false);
  vector<double> integrateSED(const flt &filter);
  static void resample(vector<oneElLambda> &lamb_all,
                       vector<oneElLambda> &lamb_new, const int origine,
                       const double lmin, const double lmax);

  void generateCalib(double lmin, double lmax, int Nsteps, int calib);
  int size() { return lamb_flux.size(); }
  void rescale(double scaleFac);
  void compute_magnitudes(const vector<flt> &filters);
  vector<double> compute_fluxes(const vector<flt> &filters);
  double trapzd();
  void sumSpectra(SED addSED, const double rescal);
  void reduce_memory(vector<flt> allFlt);

  /*
   * These functions are different depending on the type of SED
   */
  virtual void writeSED(ofstream &ofs, ofstream &ofsPhys, ofstream &ofsDoc);
  inline void writeSED(const string &binFile, const string &physFile,
                       const string &docFile) {
    ofstream sdocOut, sphysOut, sbinOut;
    sdocOut.open(docFile.c_str());
    if (!sdocOut) {
      throw invalid_argument("Can't open doc file compiling the SED " +
                             docFile);
    }
    sbinOut.open(binFile.c_str(), ios::binary | ios::out);
    if (!sbinOut) {
      throw invalid_argument("Can't open binary file compiling the SED " +
                             binFile);
    }
    if (nlib == GAL) {
      sphysOut.open(physFile.c_str());
      if (!sphysOut) {
        throw invalid_argument(
            "Can't open physical para file associated to SED " + physFile);
      }
    }
    writeSED(sbinOut, sphysOut, sdocOut);
  };

  virtual void writeMag(bool outasc, ofstream &ofsBin, ofstream &ofsDat,
                        vector<flt> allFilters, string magtyp) {};

  inline void readSEDBin(const string &fname) {
    ifstream sbinIn;
    sbinIn.open(fname.c_str(), ios::binary);
    if (!sbinIn) {
      throw invalid_argument("Can't open the binary file compiling the SED " +
                             fname);
    }
    readSEDBin(sbinIn);
  }
  virtual void readSEDBin(ifstream &ins);
  virtual void readMagBin(ifstream &ins) {};
  virtual void sumEmLines() {};
  virtual void kcorrec(const vector<double> &magz0) {};
  virtual void add_neb_cont() {};  // Add continuum
  virtual void calc_ph() {};

  virtual void SEDproperties() {};

  void generate_spectra(double zin = 0.0, double dmin = 1.0,
                        vector<opa> opaAll = {});

  virtual void clean() {
    lamb_flux.clear();
    mag.clear();
    kcorr.clear();
    fac_line.clear();
  };

  void fit_normalization(const onesource &source, const int imagm);

  inline bool is_same_model(const SED &other) {
    return ((*this).nummod == other.nummod && (*this).ebv == other.ebv &&
            (*this).age == other.age);
  }

  pair<vector<double>, vector<double>> get_data_vector(double minl, double maxl,
                                                       bool mag,
                                                       double offset = 0.0);

  void redshift();

  void applyExt(const double ebv, const ext &obj);

  void applyExtLines(const ext &obj);

  void applyOpa(const vector<opa> &opaAll);

  inline void emplace_back(const double lambda, const double value) {
    lamb_flux.emplace_back(lambda, value, 1);
  }

  inline void set_vector(const vector<double> &x, const vector<double> &y) {
    if (x.size() != y.size()) throw runtime_error("vector sizes are different");
    lamb_flux.clear();
    for (size_t k = 0; k < x.size(); k++) {
      emplace_back(x[k], y[k]);
    }
  }
};

class GalSED : public SED {
 public:
  vector<double> flEm;
  string format;
  double tau, zmet, d4000,
      fracEm;  //< fraction of the emmission line considered

  GalSED(SED const &p) : SED(p) { nlib = GAL; };
  GalSED(GalSED const &p) : SED(p) {
    flEm = p.flEm;
    format = p.format;
    tau = p.tau;
    zmet = p.zmet;
    lnir = p.lnir;
    luv = p.luv;
    lopt = p.lopt;
    d4000 = p.d4000;
    fracEm = p.fracEm;
  };

  GalSED(const string nameC, int nummodC = 0);
  GalSED(const string nameC, double tauC, double ageC, string formatC,
         int nummodC, string typeC, int idAgeC);
  ~GalSED() { flEm.clear(); }

  void SEDproperties();
  void add_neb_cont();
  GalSED generateEmSED(const string &emtype);
  void generateEmEmpUV(double MNUV_int, double NUVR);
  void generateEmEmpSFR(double MNUV_int, double NUVR);
  void generateEmPhys(double zmet, double qi);
  void generateEmSpectra(int nstep);
  void sumEmLines();

  void kcorrec(const vector<double> &magz0);
  void rescaleEmLines();
  void zdepEmLines(int flag);
  void calc_ph();

  void writeSED(ofstream &ofs, ofstream &ofsPhys, ofstream &ofsDoc);
  void readSEDBin(ifstream &ins);

  void writeMag(bool outasc, ofstream &ofsBin, ofstream &ofsDat,
                vector<flt> allFilters, string magtyp) const;
  void readMagBin(ifstream &ins);

  void clean() {
    SED::clean();
    flEm.clear();
  }
};

class QSOSED : public SED {
 public:
  QSOSED(SED const &p) : SED(p) { nlib = QSO; };
  QSOSED(QSOSED const &p) : SED(p){};
  QSOSED(const string nameC, int nummodC = 0) : SED(nameC, nummodC, "QSO"){};
  ~QSOSED(){};

  void writeMag(bool outasc, ofstream &ofsBin, ofstream &ofsDat,
                vector<flt> allFilters, string magtyp) const;

  void readMagBin(ifstream &ins);
};

class StarSED : public SED {
 public:
  StarSED(SED const &p) : SED(p) { nlib = STAR; };
  StarSED(StarSED const &p) : SED(p){};
  StarSED(const string name, int nummod = 0) : SED(name, nummod, "STAR") { ; }
  ~StarSED() { ; }

  void writeMag(bool outasc, ofstream &ofsBin, ofstream &ofsDat,
                vector<flt> allFilters, string magtyp) const;
  void readMagBin(ifstream &ins);
};

#endif