profiles.cc

Go to the documentation of this file.

#include "../include/clumps.h"
#include "../include/params.h"
#include "../include/profiles.h"

using namespace std;
#include <iostream>
#include <math.h>
#include <stdlib.h>

//______________________________________________________________________________
double factor_r_2_to_rscale(double par_prof[4])
{
   //--- Returns F = r_2 / rscale
   //  par_prof[0]    Shape parameter #1
   //  par_prof[1]    Shape parameter #2
   //  par_prof[2]    Shape parameter #3
   //  par_prof[3]    card_profile [gENUM_PROFILE]
   //
   // N.B.: r_2 corresponds to d/dr(r^2 rho(r)) |(r=r_2) = 0
   // For a generic profile (alpha,beta,gamma), we obtain
   //    EINASTO  : rscale=r_2 => F = 1
   //    EINASTO_N: rscale=re  =>  F = pow(2n/dn,n) [Eq. 8 Graham et al., AJ 132, 2701 (2006)]
   //    ZHAO     : rscale=rs  =>  F = pow((beta-2)/(2-gamma),-1/alpha)
   // N.B.: r_2 is not defined for kZHAO if  beta<=2, i.e. outer slope (par_prof[1]) => abort

   int card_profile = (int)par_prof[3];
   switch (card_profile) {
      case kEINASTO:
         return 1.;
      case kEINASTO_N:
         return pow(2. * par_prof[0] / (3.*par_prof[0] - 0.3333333333 + 0.0079 / par_prof[0]), par_prof[0]);
      case kZHAO:
         if (par_prof[1] > 2.)
            return pow((par_prof[1] - 2.) / (2. - par_prof[2]), -1. / par_prof[0]);
         else {
            cout << ">>> r_{-2} is not defined for kZHAO if  beta<=2 => abort() " << endl;
            abort();
         }
      default :
         cout << ">>> card_profile chosen does not correspond to any profile => abort() " << endl;
         abort();
   }
   return 0.;
}

//______________________________________________________________________________
string get_name_prof(double par_prof[4])
{
   //--- Returns a name as "profile name (#1,#2,#3)"
   //  par_prof[0]    Shape parameter #1
   //  par_prof[1]    Shape parameter #2
   //  par_prof[2]    Shape parameter #3
   //  par_prof[3]    card_profile [gENUM_PROFILE]
   //
   int card_profile = (int)par_prof[3];
   char tmp[100];
   string name = gNAMES_PROFILE[card_profile];

   switch (card_profile) {
      case kEINASTO:
         sprintf(tmp, "%s (#alpha=%.2lf)", name.c_str(), par_prof[0]);
         return tmp;
      case kEINASTO_N:
         sprintf(tmp, "%s (n=%d)", name.c_str(), (int)par_prof[0]);
         return tmp;
      case kZHAO:
         sprintf(tmp, "%s (%.1lf,%.1lf,%.1lf)",
                 name.c_str(), par_prof[0], par_prof[1], par_prof[2]);
         return tmp;
      default :
         cout << ">>> card_profile chosen does not correspond to any profile => abort() " << endl;
         abort();
   }
   return "";
}

//______________________________________________________________________________
double get_rsat(double par[6])
{
   //--- Returns r_sat [kpc] for which rho(r<r_sat)=rho_sat
   //  par[0]    Dark matter density normalisation [Msol/kpc^3]
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]

   int profile = (int)par[5];
   switch (profile) {
      case kEINASTO:
         return 0.;
      case kEINASTO_N:
         return 0.;
      case kZHAO:
         if (par[4] > 0.)
            return par[1] * pow(par[0] / gDM_RHOSAT, 1. / par[4]);
         else
            return 0.;
      default :
         cout << ">>> get_rsat: card_profile " << par[5]
              << " does not correspond to any profile => abort() " << endl;
         abort();
   }
   return 0.;
}

//______________________________________________________________________________
double get_rvir(double par[6], double const& rho_vir)
{
   //--- Returns Rvir, such that rho(Rvir)=rho_vir
   //  par[0]    Dark matter density normalisation [Msol/kpc^3]
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]
   //  rho_vir   Alleged DM density at the virial radius [Msol/kpc^3]

   double rmin = par[1] / 10.;
   double rmax = 1.e3;
   int n = 1000;
   double r = rmin;
   double dr = pow(rmax / rmin, 1. / (double)n);
   for (int i = 0; i < n; ++i) {
      r *= dr;
      double res = 0.;
      rho(r, par, res);
      if (res < rho_vir)
         return r;
   }
   return 0.;
}
//______________________________________________________________________________
double rho_EINASTO(double &r, double par[3])
{
   //--- Returns rho(r) for Einasto profile (Navarro et al. 2004 or Springel et al. 2008) [Msol/kpc^3]
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    rho_{-2}: density normalisation [Msol/kpc^3]
   //  par[1]    r_{-2}: radius at which the slope in -2 (replaces rs for this profile) [kpc]
   //  par[2]    alpha: shape parameter of the profile
   //
   //    => rho(r)= rho_{-2} * exp{-2/alpha [(r/r_{-2})^{alpha} - 1] }

   double res = par[0] * exp(-2. / par[2] * (pow(r / par[1], par[2]) - 1.));
   if (res < gDM_RHOSAT)
      return res;
   else
      return gDM_RHOSAT;
}

//______________________________________________________________________________
double rho_EINASTO_N(double &r, double par[3])
{
   //--- Returns rho(r) for Einasto r^{1/n} profile (Merritt et al. 2006) [Msol/kpc^3]
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    rho_e: density normalisation [Msol/kpc^3]
   //  par[1]    r_e: radius containing half the mass of the halo (replaces rs for this profile) [kpc]
   //  par[2]    n: integer index of the profile (shape parameter)
   //
   //    => rho(r)= rho_e * exp{-dn [(r/r_e)^{1/n} - 1] }



   if (par[2] < 0.5) {
      cout << "WARNING in rho_EINASTO_N: n < 0.5 not authorised in CLUMPY" << endl;
      abort();
   }
   double dn = 3.*par[2] - 0.3333333333 + 0.0079 / par[2];

   double res = par[0] * exp(-dn * (pow((r / par[1]), (1. / par[2])) - 1.));
   if (res < gDM_RHOSAT)
      return res;
   else
      return gDM_RHOSAT;
}

//______________________________________________________________________________
double rho_ZHAO(double &r, double par[5])
{
   //--- Returns rho(r) for a ZHAO (alpha,beta,gamma) [Msol/kpc^3]
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    rho_s: density normalisation [Msol/kpc^3]
   //  par[1]    r_s: scale radius [kpc]
   //  par[2]    alpha: transition slope
   //  par[3]    beta: outer slope
   //  par[4]    gamma: inner slope
   //
   //    => rho(r)= rho_s/[(r/r_s)^{gamma} [1+(r/r_s)^{alpha}]^{(beta-gamma)/alpha}]


   // Required to avoid divergences at r=0
   if (par[4] > 1.e-5) {
      if (r > par[1] * pow(par[0] / gDM_RHOSAT, 1. / par[4])) {
         return par[0] * pow(par[1] / r, par[4])
                / pow(1. + pow(r / par[1], par[2]), (par[3] - par[4]) / par[2]);
      } else
         return gDM_RHOSAT;
   } else
      return par[0] / pow(1. + pow(r / par[1], par[2]), par[3] / par[2]);
}

//______________________________________________________________________________
void rho(double &r, double par[6], double& res)
{
   //--- Sets res = Density profile [Msol/kpc^3]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    Normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]
   // Output:
   //  res       Value of rho(r)

   int profile1 = (int)par[5];

   // Calculate rho(r)
   switch (profile1) {
      case kEINASTO:
         res = rho_EINASTO(r, par);
         break;
      case kEINASTO_N:
         res = rho_EINASTO_N(r, par);
         break;
      case kZHAO:
         res = rho_ZHAO(r, par);
         break;
      default :
         cout << ">>> rho1: card_profile " << par[5]
              << " does not correspond to any profile => abort() " << endl;
         abort();
   }
   return;
}

//______________________________________________________________________________
void rho_mix(double &r, double par[21], double& res)
{
   //--- Sets res = Density profile [Msol/kpc^3]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    rho_1(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    rho_1(r) scale radius [kpc]
   //  par[2]    rho_1(r) shape parameter #1
   //  par[3]    rho_1(r) shape parameter #2
   //  par[4]    rho_1(r) shape parameter #3
   //  par[5]    rho_1(r) card_profile [gENUM_PROFILE]
   //  par[6]    UNUSED
   //  par[7]    UNUSED
   //  par[8]    UNUSED
   //  par[9]    UNUSED
   //  par[10]   switch_rho - selects which combination of rho_1 and rho_2 to use
   //                        0 -> rho(r) = rho1(r)
   //                        1 -> rho(r) = rho1(r)-rho2(r)
   //                        2 -> rho(r) = rho1(r)*rho2(r)
   //  par[11]   rho_2(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[12]   rho_2(r) scale radius [kpc]
   //  par[13]   rho_2(r) shape parameter #1
   //  par[14]   rho_2(r) shape parameter #2
   //  par[15]   rho_2(r) shape parameter #3
   //  par[16]   rho_2(r) card_profile [gENUM_PROFILE]
   //  par[17]   UNUSED
   //  par[18]   UNUSED
   //  par[19]   UNUSED
   //  par[20]   UNUSED
   // Output:
   //  res       Value of rho(r)

   int switch_rho = (int)par[10];
   int profile1 = (int)par[5];

   // Calculate rho1(r)
   switch (profile1) {
      case kEINASTO:
         res = rho_EINASTO(r, par);
         break;
      case kEINASTO_N:
         res = rho_EINASTO_N(r, par);
         break;
      case kZHAO:
         res = rho_ZHAO(r, par);
         break;
      default :
         cout << ">>> rho_1: card_profile " << par[5]
              << " does not correspond to any profile => abort() " << endl;
         abort();
   }


   // Check combi and calculate rho2 if necessary
   if (switch_rho == 0) {
      if (res > gDM_RHOSAT)
         res = gDM_RHOSAT;
      return;
   } else {
      int profile2 = (int)par[16];
      double rho_2 = 0.;
      switch (profile2) {
         case kEINASTO:
            rho_2 = rho_EINASTO(r, &par[11]);
            break;
         case kEINASTO_N:
            rho_2 = rho_EINASTO_N(r, &par[11]);
            break;
         case kZHAO:
            rho_2 = rho_ZHAO(r, &par[11]);
            break;
         default :
            cout << ">>> rho_2: card_profile " << par[16]
                 << " does not correspond to any profile => abort() " << endl;
            abort();
      }

      if (switch_rho == 1) {
         if (rho_2 < gDM_RHOSAT)
            res  -= rho_2;
         return;
      } else if (switch_rho == 2) {
         res  *= rho_2;
         return;
      } else {
         cout << ">>> switch_rho:  " << par[1]
              << " does not correspond to any combi for rho_1 and rho_2 => abort() " << endl;
         abort();
      }
   }
   return;
}

//______________________________________________________________________________
void rho2(double &r, double par[6], double& res)
{
   //--- Sets res = Density square of the profile [Msol^2/kpc^6]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    Normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]
   // Output:
   //  res       Value of rho^2(r)

   rho(r, par, res);
   res *= res;
   return;
}

//______________________________________________________________________________
void rho2_mix(double &r, double par[21], double& res)
{
   //--- Sets res = Density square of the profile [Msol^2/kpc^6]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    rho_1(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    rho_1(r) scale radius [kpc]
   //  par[2]    rho_1(r) shape parameter #1
   //  par[3]    rho_1(r) shape parameter #2
   //  par[4]    rho_1(r) shape parameter #3
   //  par[5]    rho_1(r) card_profile [gENUM_PROFILE]
   //  par[6]    UNUSED
   //  par[7]    UNUSED
   //  par[8]    UNUSED
   //  par[9]    UNUSED
   //  par[10]   switch_rho - selects which combination of rho_1 and rho_2 to use
   //                        0 -> rho(r) = rho1(r)
   //                        1 -> rho(r) = rho1(r)-rho2(r)
   //                        2 -> rho(r) = rho1(r)*rho2(r)
   //  par[11]   rho_2(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[12]   rho_2(r) scale radius [kpc]
   //  par[13]   rho_2(r) shape parameter #1
   //  par[14]   rho_2(r) shape parameter #2
   //  par[15]   rho_2(r) shape parameter #3
   //  par[16]   rho_2(r) card_profile [gENUM_PROFILE]
   //  par[17]   UNUSED
   //  par[18]   UNUSED
   //  par[19]   UNUSED
   //  par[20]   UNUSED
   // Output:
   //  res       Value of rho^2(r)

   rho_mix(r, par, res);
   res *= res;
   return;
}

//______________________________________________________________________________
void r2rho(double &r, double par[6], double& res)
{
   //--- Sets res = r^2 * Density of the profile [Msol/kpc]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    Normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]
   // Output:
   //  res       Value of r^2*rho(r)

   rho(r, par, res);
   res *= r * r;
   return;
}

//______________________________________________________________________________
void r2rho_mix(double &r, double par[21], double& res)
{
   //--- Sets res = r^2 * Density of the profile [Msol/kpc]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    rho_1(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    rho_1(r) scale radius [kpc]
   //  par[2]    rho_1(r) shape parameter #1
   //  par[3]    rho_1(r) shape parameter #2
   //  par[4]    rho_1(r) shape parameter #3
   //  par[5]    rho_1(r) card_profile [gENUM_PROFILE]
   //  par[6]    UNUSED
   //  par[7]    UNUSED
   //  par[8]    UNUSED
   //  par[9]    UNUSED
   //  par[10]   switch_rho - selects which combination of rho_1 and rho_2 to use
   //                        0 -> rho(r) = rho1(r)
   //                        1 -> rho(r) = rho1(r)-rho2(r)
   //                        2 -> rho(r) = rho1(r)*rho2(r)
   //  par[11]   rho_2(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[12]   rho_2(r) scale radius [kpc]
   //  par[13]   rho_2(r) shape parameter #1
   //  par[14]   rho_2(r) shape parameter #2
   //  par[15]   rho_2(r) shape parameter #3
   //  par[16]   rho_2(r) card_profile [gENUM_PROFILE]
   //  par[17]   UNUSED
   //  par[18]   UNUSED
   //  par[19]   UNUSED
   //  par[20]   UNUSED
   // Output:
   //  res       Value of r^2*rho(r)

   rho_mix(r, par, res);
   res *= r * r;
   return;
}

//______________________________________________________________________________
void r2rho2(double &r, double par[6], double& res)
{
   //--- Sets res = r^2 * Density square of the profile [Msol^2/kpc^4]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    Normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]
   // Output:
   //  res       Value of r^2*rho^2(r)


   rho(r, par, res);
   res *= res * r * r;
   return;
}

//______________________________________________________________________________
void r2rho2_mix(double &r, double par[21], double& res)
{
   //--- Sets res = r^2 * Density square of the profile [Msol^2/kpc^4]
   // Input:
   //  r         Distance from the centre of the profile [kpc]
   //  par[0]    rho_1(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[1]    rho_1(r) scale radius [kpc]
   //  par[2]    rho_1(r) shape parameter #1
   //  par[3]    rho_1(r) shape parameter #2
   //  par[4]    rho_1(r) shape parameter #3
   //  par[5]    rho_1(r) card_profile [gENUM_PROFILE]
   //  par[6]    UNUSED
   //  par[7]    UNUSED
   //  par[8]    UNUSED
   //  par[9]    UNUSED
   //  par[10]   switch_rho - selects which combination of rho_1 and rho_2 to use
   //                        0 -> rho(r) = rho1(r)
   //                        1 -> rho(r) = rho1(r)-rho2(r)
   //                        2 -> rho(r) = rho1(r)*rho2(r)
   //  par[11]   rho_2(r) normalisation: density [Msol/kpc^3] or proxy for dPdV [kpc^{-3}]
   //  par[12]   rho_2(r) scale radius [kpc]
   //  par[13]   rho_2(r) shape parameter #1
   //  par[14]   rho_2(r) shape parameter #2
   //  par[15]   rho_2(r) shape parameter #3
   //  par[16]   rho_2(r) card_profile [gENUM_PROFILE]
   //  par[17]   UNUSED
   //  par[18]   UNUSED
   //  par[19]   UNUSED
   //  par[20]   UNUSED
   // Output:
   //  res       Value of r^2*rho^2(r)

   rho_mix(r, par, res);
   res *= res * r * r;
   return;
}

//______________________________________________________________________________
double r_2_to_rscale(double const& r_2, double par_prof[4])
{
   //--- Returns scale radius rs [kpc] from r_2 [kpc] for a given profile.
   //  r_2            Radius [kpc] for which the DM profile  slope = -2
   //  par_prof[0]    Shape parameter #1
   //  par_prof[1]    Shape parameter #2
   //  par_prof[2]    Shape parameter #3
   //  par_prof[3]    card_profile [gENUM_PROFILE]

   return r_2 / factor_r_2_to_rscale(par_prof);
}

//______________________________________________________________________________
double rscale_to_r_2(double const& rs, double par_prof[4])
{
   //--- Returns r_2 [kpc] from the scale radius rs [kpc] for a given profile
   //  r_s            Scale radius [kpc] (rs, re or r{-2}, depending on card_profile)
   //  par_prof[0]    Shape parameter #1
   //  par_prof[1]    Shape parameter #2
   //  par_prof[2]    Shape parameter #3
   //  par_prof[3]    card_profile [gENUM_PROFILE]

   return rs * factor_r_2_to_rscale(par_prof);
}

//______________________________________________________________________________
void set_par0_given_mref(double par[6], double const& r_ref, double const& m_ref, double const& eps)
{
   //--- Returns DM density normalisation [Msol/kpc^3] matching the constraint M(r<r_ref)=m_ref
   // Input:
   //  par[0]    UNUSED
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]
   //  r_ref     Reference radius [kpc]
   //  m_ref     Mass [Msol] within the radius r_ref
   //  eps       Relative precision sought for integration
   // Output:
   //  par[0]    Dark matter density normalisation [Msol/kpc^3]

   const int npar = 7;
   double par_tmp[npar] = {1., par[1], par[2], par[3], par[4], par[5], r_ref};
   par[0] = m_ref / mass_singlehalo(par_tmp, eps);
}

//______________________________________________________________________________
void set_par0_given_rhoref(double par[6], double& r_ref, double const& rho_ref)
{
   //--- Sets normalisation par[0] such that rho(r_ref)=rho_ref
   // Input:
   //  par[0]    UNUSED
   //  par[1]    Scale radius [kpc]
   //  par[2]    Shape parameter #1
   //  par[3]    Shape parameter #2
   //  par[4]    Shape parameter #3
   //  par[5]    card_profile [gENUM_PROFILE]
   //  r_ref     Reference radius [kpc]
   //  rho_ref   Density at r=r_ref [Msol/kpc^3]
   // Output:
   //  par[0]    Dark matter density normalisation [Msol/kpc^3]

   // Calculate normalisation
   double res = 0.;
   par[0] = 1.;
   rho(r_ref, par, res);
   par[0] = rho_ref / res;
   return;
}