source: Sophya/trunk/Cosmo/SimLSS/schechter.cc@ 3267

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "pexceptions.h"

#include "constcosmo.h"
#include "schechter.h"

///////////////////////////////////////////////////////////
//***************** Schechter Functions *****************//
///////////////////////////////////////////////////////////

// HI mass function:
//    see M.Zwaan astroph-0502257 (MNRAS Letters, Volume 359, Issue 1, pp. L30-L34.)
//    see M.Zwaan astroph-9707109 (ApJ, Volume 509, Issue 2, pp. 687-702.)

Schechter::Schechter(double nstar,double mstar,double alpha)
  : nstar_(nstar) , mstar_(mstar) , alpha_(alpha) , outvalue_(0)
{
}

Schechter::Schechter(Schechter& f)
  : nstar_(f.nstar_) , mstar_(f.mstar_) , alpha_(f.alpha_) , outvalue_(f.outvalue_)
{
}

Schechter::~Schechter(void)
{
}

void Schechter::SetOutValue(unsigned short outvalue)
// outvalue = 0 : give dn/dm
//          = 1 : give m*dn/dm
{
  if(outvalue>1) {
    cout<<"Schechter::SetOutValue: Error bad outvalue: "<<outvalue<<endl;
    throw ParmError("Schechter::SetOutValue: Error bad outvalue");
  }
  outvalue_ = outvalue;
}

double Schechter::operator() (double m)
// Return : "dn/dm = f(m)" or "m*dn/dm = f(m)"
{
  double x = m/mstar_;
  double dndm = nstar_/mstar_ * pow(x,alpha_) * exp(-x);
  if(outvalue_) dndm *= m;  // on veut m*dn/dm
  return dndm;
}

void Schechter::Print(void)
{
  cout<<"Schechter::Print: nstar="<<nstar_<<" Mpc^-3"
      <<"  mstar="<<mstar_<<" MSol"
      <<"  alpha="<<alpha_
      <<"  (outvalue="<<outvalue_<<" -> return ";
  if(outvalue_) cout<<"m*dn/dm)"; else cout<<"dn/dm)";
  cout<<endl;
}

///////////////////////////////////////////////////////////
//******************* Le Flux a 21cm ********************//
///////////////////////////////////////////////////////////

double Msol2FluxHI(double m,double d)
// Input:
//    m : masse de HI en "Msol"
//    d : distance en "Mpc"  (si cosmo d=DLum(z))
// Return:
//    le flux total emis a 21 cm en W/m^2
// Ref:
// -- Binney & Merrifield, Galactic Astronomy p474 (ed:1998)
// S(W/m^2) = 1e-26 * Nu_21cm(Hz) * m(en masse solaire) /(2.35e5 * C(km/s) * Dlum^2)
//          = 2.0e-28 * m / Dlum^2
// -- J.Peterson & K.Bandura, astroph-0606104  (eq 7)
//    F.Abdalla & Rawlings, astroph-0411342 (eq 7 mais pb de d'unites?)
//          (A_21cm = 2.86888e-15 s^-1)
// S(W/m^2) = 3/4 * h(J.s) * Nu_21cm(Hz) * A_21cm(s^-1) * Msol(kg)/m_proton(kg)
//                         / (4 Pi) / (Dlum(Mpc) * 3.0857e22(m/Mpc))^2
//          = 2.0e-28 * m / Dlum^2
//-------------
{
  return  2.0e-28 * m / (d*d);
}

double FluxHI2Msol(double f,double d)
// Input:
//    f : flux total emis a 21 cm en W/m^2
//    d : distance en "Mpc"  (si cosmo d=DLum(z))
// Return:
//    m : masse de HI en "Msol"
{
  return f *d*d / 2.0e-28;
}