Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Springel Model - Initial conditions of galaxy for simulation

  1. May 1, 2013 #1

    I'am studying a code (called "starscream") which allows to make initial conditions (positions and velocities) for NBody simulation. They are based on Springel and White (1999) model. I have several problems about the underlying equations used in this code.

    Firstly, See attachment "fig1.png", for computing the disk scale length, into the "disk_scale_length_func" routine, the author seems to have forgotten to multiply the line 128 by the line 129, i.e the
    factor "sqrt(1.0/f_c_func(gal->halo_concentration))".

    Indeed, the definition of the disk scale length is :


    where [tex] f_{R}(\lambda,c,m_{d},j_{d})=2\bigg[\int_{0}^{\infty}\,e^{-u}\,u^{2}\,\dfrac{v_{c}(R_{d}\,u)}{v_{200}}\,du\bigg]^{-1}[/tex]

    and [tex] f_{c}=\dfrac{c}{2}\,\dfrac{1-1/(+c)^{2}-2\,(ln(1+c))/(1+c)}{[c/(1+c)-ln(1+c)]^{2}} [/tex]

    Here is the original part of this code :

    Code (Text):

    // This function determines the scale length of the disk to within 15% or so using
    // the fitting formula provided in MMW97.
    double disk_scale_length_func(galaxy *gal) {

      int i;
      double base, power, c, f_conc, f, scale, old_scale, v_c;

      c = gal->halo_concentration;
      f_conc = 2.0/3.0+ pow((c/21.5),0.7);
      base = (gal->j_d*gal->lambda)/(0.1*gal->m_d);
      power = (-0.06+2.71*gal->m_d+0.0047*gal->m_d/(gal->j_d*gal->lambda));
      f = pow(base,power)*(1.0-3.0*gal->m_d+5.2*gal->m_d*gal->m_d) *
      gal->disk_scale_length = (1.0/sqrt(2.0))*(gal->j_d/gal->m_d)*gal->lambda*

      // Now try to iterate for the correct value.
      for (i = 0; i < 50; ++i) {
        scale = sqrt(1.0/2.0)*(gal->j_d/gal->m_d)*gal->lambda*gal->r200*
        [B]2.0*(gal->v200/j_d_func(gal))             // product missing "*" with the next line
        gal->disk_scale_length = scale;
        printf("scale = %f\n",scale);

      return scale;

    Is anyone can confirm this error ?


    Secondly, See attachment "fig2.png", into the "dj_d_func" routine, I don't understand the lines 199 to 201. Actually, from this equation in the doc :

    [/tex] (eq 1)

    this should be :

    v_{c}^{2}=\dfrac{GM_{d}(r)}{r} \rightarrow M_{d}(r)=\dfrac{v_{c}^{2}\,r}{G}

    with [tex]v_{c}^{2}[/tex] computed from the (eq 1) equation.

    The mass density of the disk is :

    [/tex] (eq 2)

    With the expression of [tex]\Sigma_{0}[/tex] on (eq 2), we can have :


    So into the code, on line 199, with [tex]y=R/(2\,R_{d})[/tex] and we should have :


    and not only one product "y*y".

    Here is the original part of this code (with in bold my correction) :

    Code (Text):

    double dj_d_func(double radius, void *params) {

      int status;
      galaxy *gal = (galaxy *) params;
      double I_0, K_0, I_1, K_1, halo_mass, disk_mass, y, a;
      double total_disk_mass, v_c, scale;
      gsl_sf_result result;

      scale = gal->disk_scale_length;
      total_disk_mass = gal->disk_mass;
      a = gal->halo_a_value;
      halo_mass = gal->halo_mass;

      y = radius/(2.0*scale);
      status = gsl_sf_bessel_I0_e(y,&result);
      I_0 = result.val;
      status = gsl_sf_bessel_K0_e(y,&result);
      K_0 = result.val;
      status = gsl_sf_bessel_I1_e(y,&result);
      I_1 = result.val;
      status = gsl_sf_bessel_K1_e(y,&result);
      K_1 = result.val;

      disk_mass = 2.0*total_disk_mass*y*y*(I_0*K_0 - I_1*K_1);

      //disk_mass = 4.0*total_disk_mass*y*y*y*(I_0*K_0 + I_1*K_1);

      halo_mass = halo_mass*radius*radius/((radius+a)*(radius+a));
      v_c = sqrt(G*(halo_mass/(radius*kpc) + disk_mass/(scale*kpc)));

      return v_c*(radius/scale)*(radius/scale)*exp(-radius/scale);
    Is my modification correct ?

    I make you notice too that there is a minus sign on the original code at line 199 : (I_{0}*K_{0} - I_{1}*K_{1}) : is this good ?

    Thanks for your help

    Attached Files:

  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted