00001 /* This file is part of Cloudy and is copyright (C)1978-2010 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "physconst.h" 00005 #include "dense.h" 00006 #include "pressure.h" 00007 #include "radius.h" 00008 #include "colden.h" 00009 #include "dark_matter.h" 00010 #include "cosmology.h" 00011 00012 void GravitationalPressure( void ) 00013 { 00014 DEBUG_ENTRY( "GravitationalPressure()" ); 00015 00016 double M_dark, M_self, M_external, g_dark, g_self, g_external; 00017 double R = radius.Radius - radius.dRadSign * radius.drad / 2; 00018 00019 M_dark = M_self = M_external = g_dark = g_self = g_external = 0.; 00020 00021 if( dark.lgNFW_Set ) 00022 { 00023 double rho_crit = 3.*POW2(cosmology.H_0*1e5/MEGAPARSEC)/(8.*PI*GRAV_CONST); 00024 double c_200 = dark.r_200/dark.r_s; 00025 ASSERT( c_200 > 0. ); 00026 double delta_c = (200./3.) * POW3( c_200 ) / ( log(1.+c_200) - c_200/(1.+c_200) ); 00027 M_dark = 4*PI*rho_crit*delta_c*POW3(dark.r_s); 00028 M_dark *= 1./(1.+R/dark.r_s) + log(1.+R/dark.r_s) - 1.; 00029 00030 g_dark = -GRAV_CONST * M_dark / (R*R); 00031 } 00032 00033 /* (self-)Gravity forces: Yago Ascasibar (UAM, Spring 2009) */ 00034 00035 // add all external mass components 00036 for( unsigned int i=0; i < pressure.external_mass[0].size(); i++ ) 00037 { 00038 double M_i = pressure.external_mass[0][i]; 00039 if( R < pressure.external_mass[1][i] ) 00040 M_i *= pow( R / pressure.external_mass[1][i], pressure.external_mass[2][i] ); 00041 M_external += M_i; 00042 } 00043 00044 // evaluate gravitational acceleration 00045 switch( pressure.gravity_symmetry ) 00046 { 00047 case -1: // no self-gravity 00048 break; 00049 00050 case 0: // spherical symmetry 00051 M_self = dense.xMassTotal - dense.xMassDensity * radius.dVeffVol; 00052 // dense.xMassTotal has already been updated in radius_increment.cpp 00053 M_self *= 4*PI* radius.rinner*radius.rinner; 00054 M_self *= pressure.self_mass_factor; 00055 00056 g_self = - GRAV_CONST * M_self / (R*R); 00057 g_external = - GRAV_CONST * M_external * SOLAR_MASS / (R*R); 00058 break; 00059 00060 case 1: // mid-plane symmetry 00061 M_self = colden.TotMassColl + dense.xMassDensity * radius.drad_x_fillfac / 2.; 00062 // colden.TotMassColl will be updated later on in radius_increment.cpp 00063 M_self *= pressure.self_mass_factor; 00064 M_self *= 2.; // mid-plane symmetry 00065 00066 // Gravitational acceleration due to an infinite slab is independent of distance from 00067 // the slab and is given by g = - 2 PI G * sigma. 00068 g_self = - 2*PI* GRAV_CONST * M_self; 00069 g_external = - 2*PI* GRAV_CONST * M_external * SOLAR_MASS/PARSEC/PARSEC; 00070 00071 if( dark.lgNFW_Set ) 00072 fprintf( ioQQQ, " WARNING: Setting both mid-plane baryonic gravity symmetry and an NFW dark matter halo is almost certainly unphysical!\n" ); 00073 break; 00074 00075 default: 00076 fprintf( ioQQQ, " Unknown gravitational symmetry = %d !!!\n", pressure.gravity_symmetry ); 00077 TotalInsanity(); 00078 }; 00079 00080 pressure.RhoGravity_dark = dense.xMassDensity * g_dark * radius.drad_x_fillfac; 00081 pressure.RhoGravity_self = dense.xMassDensity * g_self * radius.drad_x_fillfac; 00082 pressure.RhoGravity_external = dense.xMassDensity * g_external * radius.drad_x_fillfac; 00083 pressure.RhoGravity = pressure.RhoGravity_dark + pressure.RhoGravity_self + pressure.RhoGravity_external; 00084 00085 return; 00086 }