#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;
}

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

double FluxHI(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)
//      J.Peterson & K.Bandura, astroph-0606104  (eq 7)
//      F.Abdalla & Rawlings, astroph-0411342 (eq 7 may pb de d'unites?)
{
  return  m / (2.35e+5 * d*d )
          * Jansky2Watt_cst
          * (Fr_HyperFin_Par*1e+9)/SpeedOfLight_Cst;
}