0

Essentially, what i want to do is set up the initial conditions for an n-body simulation of a galaxy. The paper i tried to go by is found here http://arxiv.org/abs/1204.0513 .

The paper describes two density functions:

  • a surface density function of form rho(r) for density at distance r from the center

  • and a density function of form rho(r, z) with z being the z component in (x,y,z) coordinates.

What i did so far for each particle (star) was:

  • use a monte carlo method with the first function to get a distance from the galaxy center
  • use a monte carlo method with the second function to get a z component
  • generate x,y coords randomly on a circle of radius r, at height z

But perhaps the code for this would be more helpful:

/* generate 3 random numbers x,y,z having a domain from 0 to 1 (not including 0.0) */

// The distance from the galaxy center for the current star.
// Rs is a constant parameter in the surface density function.
float r = -Rs * log(1.0f - x);

// The z coordinate for the current star
float Sz = -0.5f * (0.1f * Rs) * log(-((z-1)/z));

// The x,y coordinates for the current star
float Sx = sqrt(r*r - Sz*Sz) * cos(2.0f * PI * y);
float Sy = sqrt(r*r - Sz*Sz) * sin(2.0f * PI * y);

Visually, this produces the expected shape, until the obvious situation arises, namely if the generated Sz value is higher than the generated r value.

So my question is firstly if i'm even on the right track to begin with, and if so if anyone can either suggest a correction mechanism for the case i described above or an alternate method of generating these coordinates.

andrei
  • 339
  • 3
  • 12

1 Answers1

0

I think you have two basic options here. Option 1 would be acceptance/rejection - throw away the current attempt and try again if the generated Sz value pushes you over the r threshold. Option 2 would be to scale the result - normalize the three dimensions to r, then scale by a random number so the vector's length is some proportion <= r. The distribution of the scaling value would determine how dense outcomes are at different radii, a modeling decision you'd have to make.

pjs
  • 18,696
  • 4
  • 27
  • 56
  • Sscaling the result is something i'm currently looking into. Rejecting a given Sz or generating a new random z until one fits are both undesirable due to the logic behind this code. – andrei Jan 15 '14 at 08:07
  • Generating a new z until one fits would skew the distribution, the results would yield the conditional distribution of z given x and y rather than the joint distribution of x, y, and z. – pjs Jan 15 '14 at 14:58