cloudy trunk
Loading...
Searching...
No Matches
gravity.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3#include "cddefines.h"
4#include "physconst.h"
5#include "dense.h"
6#include "pressure.h"
7#include "radius.h"
8#include "colden.h"
9#include "dark_matter.h"
10#include "cosmology.h"
11
13{
14 DEBUG_ENTRY( "GravitationalPressure()" );
15
16 double M_dark, M_self, M_external, g_dark, g_self, g_external;
17 double R = radius.Radius - radius.dRadSign * radius.drad / 2;
18
19 M_dark = M_self = M_external = g_dark = g_self = g_external = 0.;
20
21 if( dark.lgNFW_Set )
22 {
23 double rho_crit = 3.*POW2(cosmology.H_0*1e5/MEGAPARSEC)/(8.*PI*GRAV_CONST);
24 double c_200 = dark.r_200/dark.r_s;
25 ASSERT( c_200 > 0. );
26 double delta_c = (200./3.) * POW3( c_200 ) / ( log(1.+c_200) - c_200/(1.+c_200) );
27 M_dark = 4*PI*rho_crit*delta_c*POW3(dark.r_s);
28 M_dark *= 1./(1.+R/dark.r_s) + log(1.+R/dark.r_s) - 1.;
29
30 g_dark = -GRAV_CONST * M_dark / (R*R);
31 }
32
33 /* (self-)Gravity forces: Yago Ascasibar (UAM, Spring 2009) */
34
35 // add all external mass components
36 for( unsigned int i=0; i < pressure.external_mass[0].size(); i++ )
37 {
38 double M_i = pressure.external_mass[0][i];
39 if( R < pressure.external_mass[1][i] )
40 M_i *= pow( R / pressure.external_mass[1][i], pressure.external_mass[2][i] );
41 M_external += M_i;
42 }
43
44 // evaluate gravitational acceleration
45 switch( pressure.gravity_symmetry )
46 {
47 case -1: // no self-gravity
48 break;
49
50 case 0: // spherical symmetry
51 M_self = dense.xMassTotal - dense.xMassDensity * radius.dVeffVol;
52 // dense.xMassTotal has already been updated in radius_increment.cpp
53 M_self *= 4*PI* radius.rinner*radius.rinner;
54 M_self *= pressure.self_mass_factor;
55
56 g_self = - GRAV_CONST * M_self / (R*R);
57 g_external = - GRAV_CONST * M_external * SOLAR_MASS / (R*R);
58 break;
59
60 case 1: // mid-plane symmetry
61 M_self = colden.TotMassColl + dense.xMassDensity * radius.drad_x_fillfac / 2.;
62 // colden.TotMassColl will be updated later on in radius_increment.cpp
63 M_self *= pressure.self_mass_factor;
64 M_self *= 2.; // mid-plane symmetry
65
66 // Gravitational acceleration due to an infinite slab is independent of distance from
67 // the slab and is given by g = - 2 PI G * sigma.
68 g_self = - 2*PI* GRAV_CONST * M_self;
69 g_external = - 2*PI* GRAV_CONST * M_external * SOLAR_MASS/PARSEC/PARSEC;
70
71 if( dark.lgNFW_Set )
72 fprintf( ioQQQ, " WARNING: Setting both mid-plane baryonic gravity symmetry and an NFW dark matter halo is almost certainly unphysical!\n" );
73 break;
74
75 default:
76 fprintf( ioQQQ, " Unknown gravitational symmetry = %d !!!\n", pressure.gravity_symmetry );
78 };
79
80 pressure.RhoGravity_dark = dense.xMassDensity * g_dark * radius.drad_x_fillfac;
81 pressure.RhoGravity_self = dense.xMassDensity * g_self * radius.drad_x_fillfac;
82 pressure.RhoGravity_external = dense.xMassDensity * g_external * radius.drad_x_fillfac;
83 pressure.RhoGravity = pressure.RhoGravity_dark + pressure.RhoGravity_self + pressure.RhoGravity_external;
84
85 return;
86}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define POW3
Definition cddefines.h:936
#define POW2
Definition cddefines.h:929
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_colden colden
Definition colden.cpp:5
t_cosmology cosmology
Definition cosmology.cpp:11
t_dark_matter dark
t_dense dense
Definition dense.cpp:24
void GravitationalPressure(void)
Definition gravity.cpp:12
UNUSED const double SOLAR_MASS
Definition physconst.h:71
UNUSED const double PI
Definition physconst.h:29
UNUSED const double MEGAPARSEC
Definition physconst.h:141
UNUSED const double GRAV_CONST
Definition physconst.h:109
UNUSED const double PARSEC
Definition physconst.h:138
t_pressure pressure
Definition pressure.cpp:5
t_radius radius
Definition radius.cpp:5