This program calculates the pressure P and total mass M as a function of radius r in a spherically symmetric, static mass, in the context of general relativity. The ODE for the pressure is known as the Tolman-Oppenheimer-Volkoff (TOV) equation, and the ODE for the mass-energy comes from performing the Leibniz integral rule on the definition of m(r) during the TOV derivation. (This is shown in pretty much every GR book.) The TOV equation + mass-energy equation lead to 2 ODEs --- one for m(r) and the other for P(r) --- but in these equations we have a third variable, the mass-energy density rho(r). Thus we close the system by providing rho(r), or more commonly, rho(P, T), which is known in thermodynamics as an equation of state (EOS). Currently this program can use polytropes or any tabulated EOS data. The required form of such data is: First line: number of data entries Remaining lines: 4 columns of data 1.) (energy density)/c^2 in units of g/cm^3 2.) pressure in units of erg/cm^3 3.) (currently unused) 4.) (currently unused) If you use a polytrope, remember you have to set the linear parameter K as well as the exponent of rho, called Gamma. For neutron-degenerate matter, Gamma ~ 2.5 is a good choice. As for K, I prefer to use a "real" (tabulated) EOS, find the value of the pressure at the center of the star, and set K such that the polytrope also matches that pressure. After all, the central pressure sets the scale for most of the thermodynamic/mechanical variables in the TOV system. Setting K this way is an arbitrary choice, however, and you can do with it whatever you like. Since the TOV equation contains a rho(P, T) term along with a P(rho, T) term, I use a root-finder to invert the EOS. This method is predicated on the assumption that T(r) is known everywhere in the star. Since the TOV system does not provide any constaints on T(r), we have some freedom here: lots of people set T=0 in neutron stars, others make it isothermal, and if you want to be really fancy you could calculate T using the diffusion approximation (see any stellar evolution textbook for the derivation). The last option requires in addition a total luminosity parameter for the star. The way the program exists currently, it assumes T is wrapped up in the tabulated EOS data, but this restriction would not be difficult to lift. Among other things, this program can tell you how massive is a relativistic star in hydrostatic equilibrium. To find this, just pick some final radius for the integration and watch P(r) drop. The surface of the star R is at the radius where P(r = R) = 0. Of course the integrator has no idea that P(r) can't be negative, but the ODEs certainly do, and will result in lots of NaNs once the integrator reaches the surface. So pick a large enough 'final' radius such that P(r) at least reaches zero, and ignore everything past that point where it goes negative/NaN. Remember also that you need to specify rho(r = 0) as an initial condition. By construction M(r = 0) = 0 so that part's done. If you vary rho(r = 0) you should find that TOV predicts a single maximum mass for a given EOS. That's called the TOV limit (just like the Chandrasekhar limit for white dwarfs). REFERENCES ---------- http://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation