Source Code
Louis Marmet
2006 November 1st
You can copy, run and distribute this code as long as the name
of the original author is always included.
/* Static galactic model II
Copyright © Louis Marmet,
2005 December 28.
This code assumes the following:
- Matter distribution in a galaxy is given by the sum of two functions:
- rho_2d[ri]: mass per unit area in the disk at
radius ri,
- rho_3d[ri]: mass per unit volume in the bulge at
radius ri,
- The galaxy is in a steady state of motion,
- Gravity interaction is instantaneous
************************************* EXAMPLES
********************************
When trying these examples, watch the calculated values of chi.
// Example 1: DISK ONLY
Take file r_density_raw.txt
USE:
#define
SMOOTH_STR
1000 // smoothing strength (smaller number means stronger smoothing)
#define
LOOPS
(300) // the calculation will be
repeated this number of times
#define
BULGE
(1==0) // true is there is a Bulge
component
#define DISK
(1==1) // true
if there is a Disk component
The resulting file r_density.txt will contain the density distribution.
*************************************
// Example 2: BULGE ONLY
Take file s_density_raw.txt
USE:
#define
SMOOTH_STR
1000 // smoothing strength (smaller number means stronger smoothing)
#define
LOOPS
(300) // the calculation will be
repeated this number of times
#define
BULGE
(1==1) // true is there is a Bulge
component
#define DISK
(1==0) // true
if there is a Disk component
The resulting file s_density.txt
will contain the density distribution.
*************************************
When trying this example,
watch the calculated values of chi and the Angular momentum.
I have to stress that to
obtain a solution, there is no given parameter
that will work from the beginning to the end of the calculation.
You
MUST change these parameters manually many times, while running the
program for a few thousand iterations between changes.
Look at the program to figure that the
fluctuations happen because the program is randomly moving small
amounts of mass up and down the radius of the galaxy. If an
improvement is seen (smaller chi or total angular momentum closer to
target), the mass is left at its new location, otherwise it is returned
to the old location. The amount of mass that is displaced is
proportional to RANDOM_GAIN. If RANDOM_GAIN is too large, a lot
of
mass will be displaced, and this might lead to divergence. If
RANDOM_GAIN is too small, very little mass is moved and convergence is
slow. You have to play with the value of RANDOM_GAIN to get
reasonable
variations in chi but avoid divergence.
Also, there is some mass
transfered between the bulge and the disk. If the new
distribution is
better (smaller chi or total angular momentum closer to target), the
mass is left at its new location, otherwise it is returned to the
old location. The amount of mass that is displaced is
proportional to MASS_TRANSFER . If MASS_TRANSFER is too large, a lot of mass will be
displaced, and this might lead to
divergence. If MASS_TRANSFER is too
small, very little mass is moved and
convergence is slow. You have to play with the value of MASS_TRANSFER to
get minimal variations in the angular momentum but avoid divergence.
Finally, once in a while, the distribution functions are
smoothed
(otherwise, they become very noisy). STRONG smoothing is obtained
with
a SMALL value of SMOOTH_STR. WEAK smoothing is obtained with a
LARGE
value of SMOOTH_STR. Look at the program after
// ************** Smoothing
of density functions
to see why this is so. Strong smoothing is good except that it
prevents the mass distribution function to have sharp features.
These
sharp features are necessary to get a good solution, so eventually,
when the program converges towards a solution, you must increase the
value of SMOOTH_STR otherwise chi won't decrease. So follow the
procedure given in example 3 to get the distribution function for an
angular momentum of 85. It takes some time if the computer is not
fast.
// Example
3: DISK AND BULGE
take file sr_density_raw.txt
which is the solution for 78.0 angular momentum
USE:
#define ANGULARM_TARGET
78.0 // [M_galaxy r_galaxy km/s]
#define SMOOTH_STR
1000 // smoothing
strength (smaller number means stronger smoothing)
#define LOOPS
(300) // the
calculation will be repeated this number of times
#define BULGE
(1==1) // true is
there is a Bulge component
#define DISK
(1==1) // true
if there is a Disk component
#define RANDOM_GAIN
1
// gain on the amplitude of the test variations
#define MASS_TRANSFER
(1.0/2) // gain on the mass transfer to balance disk and
bulge angular momentum
and run the program.
One gets a value of chi about equal to 3.
Now, we want to see what the galaxy is like with an angular momentum of
85. Keep same file sr_density_raw.txt, and
USE:
#define ANGULARM_TARGET
85.0 // [M_galaxy r_galaxy km/s]
#define SMOOTH_STR
10 //
smoothing strength : MORE smoothing
#define LOOPS
(600) // the
calculation will be repeated this number of times
and run the program.
If chi oscillates a lot, use a smaller MASS_TRANSFER
#define MASS_TRANSFER
(1.0/10) // gain on the mass transfer to balance disk and bulge
angular momentum
and run the program.
If chi doesn't get smaller, it may be because the smoothing is too
strong. The smoother curve is actually not a good
solution, so decrease smoothing strength by increasing SMOOTH_STR
#define SMOOTH_STR
50 //
smoothing strength : LESS smoothing
and run the program.
It may take a lot of loops...
#define LOOPS
(3000) // the
calculation will be repeated this number of times
and run the program.
The angular momentum varies a lot, try not to transfer too much:
#define SMOOTH_STR
100
#define LOOPS
(1500) // the
calculation will be repeated this number of times
#define MASS_TRANSFER
(1.0/10) // gain on the mass transfer to balance disk and bulge
angular momentum
and run the program.
This is better, now the angular momentum is 84.956 (close to target)
and chi=2.16 [km/s].
Reduce smoothing:
#define SMOOTH_STR
300
and run the program.
Now chi=1.73 [km/s]. sr_density.txt has the
density distributions.
To reach the best solution, continue to increase SMOOTH_STR,
decrease RANDOM_GAIN and decrease MASS_TRANSFER. The proper
choice of which one has to be changed depends on what is happening with
the quantities chi and angular momentum.
The value of SMOOTH_STR should be small (3 to 10) when
attempting to find a new solution. This strong smoothing allows
the mass to distribute radially. Monitor the value of chi and if
it doesn't decrease, there is
probably too much smoothing, so increase the value of SMOOTH_STR to
reduce the smoothing strength.
If the angular momentum oscillates each time it is reported,
there is too much mass transferred between the disk and bulge. Reduce
MASS_TRANSFER to correct that.
If the value of chi varies a lot, possibly RANDOM_GAIN is too
large.
*/
#include <ansi_c.h>
// Physical constants
#define
C_PI
3.141592653589793
#define
C_G
6.673e-11 // [m^3/kg/s^2]
#define
C_C
299792458 // [m/s]
#define
C_LIGHT_YEAR
9.4605284e15 // [m]
#define
C_KPC
3.08e19
//
[m]
#define C_MSOL
1.99e30
//
[kg]
// Physical parameters
#define
GALAXY_RADIUS
7.5e20 // [m]
#define
DISK_THICKNESS
4.5e19 // [m]
#define
VELOCITY_TARGET
175000 // [m/s]
#define
ANGULARM_TARGET
92.69 // [M_galaxy
r_galaxy km/s]
// Program parameters
#define
N_RADIAL_SAMPLES
100 //
Number of sampling points, from 0 to GALAXY_RADIUS
#define
RADIUS_STEP
(GALAXY_RADIUS/N_RADIAL_SAMPLES) // step size
#define
PHI_STEP
(C_PI/19)
// smoothing strength; (smaller number means stronger smoothing)
// 3
is
good value for chi<1% of VELOCITY_TARGET, but is too strong for disk
functions
#define
SMOOTH_STR
1000 // Start with 3 and then increase to allow faster
varying functions
#define
LOOPS
(335*2*1) // the
calculation will be repeated this number of times
#define
BULGE
(1==0) // true is there is a Bulge component
#define
DISK
(1==1) // true if there is a Disk component
// gain on the amplitude of
the test variations. A too large value of the gain may slow down
convergence
#define
RANDOM_GAIN
1
#define
MASS_TRANSFER
(1.0/2) // gain on the mass transfer
to balance disk and bulge angular momentum
double
*rho_2d; //
mass per unit area (per unit radius per unit phi/r)
double
*rho_3d; //
mass per unit volume
double
*v_target; // target
velocity
// Calculates the mass of the ring ri
double mass_ring(int ri)
{
if (ri==0)
return(rho_2d[0]*C_PI*DISK_THICKNESS*RADIUS_STEP*RADIUS_STEP/4);
return(rho_2d[ri]*ri*2*C_PI*RADIUS_STEP*DISK_THICKNESS*RADIUS_STEP);
}
// Calculates the mass of the shell ri
double mass_shell(int ri)
{
if (ri==0)
return(rho_3d[0]*C_PI*RADIUS_STEP*RADIUS_STEP*RADIUS_STEP/6);
return(rho_3d[ri]*(12*ri*ri+1)*C_PI*RADIUS_STEP*RADIUS_STEP*RADIUS_STEP/3);
}
// X component of gravitational acceleration on a
point located at distance ri_m
// This will integrate the field from a matter distribution in a disk
and in the bulge
double accG_x(int
ri_m) //
[m/s^2]
{
int
ri;
double
ax_r, ax_s, phi, s, c;
double
r, r_m, dx, dy, distance2,
distance;
if (ri_m==0)
return(0);
r_m=ri_m*RADIUS_STEP;
ax_r=0;
if (DISK)
//
********** 2-D integral over phi (0 to PI, then multiply by 2)
{
for
(phi=PHI_STEP/2; phi<C_PI;
phi+=PHI_STEP)
{
s=sin(phi);
c=cos(phi);
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++) //
2-D integral over all r
{
r=ri*RADIUS_STEP;
dx=r*c-r_m;
dy=r*s;
distance2=dx*dx+dy*dy;
distance=sqrt(distance2);
ax_r+=dx/(distance*distance2)*mass_ring(ri);
}
}
ax_r*=PHI_STEP/C_PI;
}
ax_s=0;
if (BULGE)
//
********** 3-D integral over r<r_m
{
for
(ri=0; ri<ri_m; ri++)
ax_s-=mass_shell(ri);
ax_s-=mass_shell(ri)/2; // add
only half of last shell
}
ax_s/=r_m*r_m;
return(C_G*(ax_r+ax_s));
}
// Find the deviation from the target velocity
double chi(void)
{
int
ri;
double
r, a, dv, chi;
chi=0;
for (ri=1;
ri<N_RADIAL_SAMPLES; ri++)
{
r=ri*RADIUS_STEP;
a=accG_x(ri);
if (a<0)
dv=sqrt(-a*r)-v_target[ri]; //
a=v^2/r
else
dv=-sqrt(a*r)-v_target[ri];
chi+=dv*dv;
}
return(sqrt(chi/N_RADIAL_SAMPLES));
}
//
**************************************************************
main()
{
int
ri, loop,
rv, zi;
double
r, acc, chi_old, chi_new,
size, datum, intensity;
double
mass, mass_2D, mass_3D;
double
ang_mom2D, ang_mom3D, c1,
c2;
time_t t;
FILE *fh;
rho_2d=(double
*)malloc(N_RADIAL_SAMPLES*sizeof(double));
rho_3d=(double
*)malloc(N_RADIAL_SAMPLES*sizeof(double));
v_target=(double
*)malloc(N_RADIAL_SAMPLES*sizeof(double));
// ********** When program is run for the first
time, this will generate an approximate density profile.
// ********** Do not read data from file (see conditional compile
below: set to 1==0)
// Evaluate density in disk
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
{
r=ri*RADIUS_STEP;
if (r<0.1)
rho_2d[ri]=1e-26+2e-24*r/GALAXY_RADIUS; // rho_2d[] mass per unit volume
else
rho_2d[ri]=2e-21*GALAXY_RADIUS/r;
}
// Evaluate density in bulge
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
{
r=ri*RADIUS_STEP;
rho_3d[ri]=1e-22*GALAXY_RADIUS/(RADIUS_STEP+r); // rho_3d[] mass per unit volume
}
// Evaluate target velocity
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
{
r=ri/(double)(N_RADIAL_SAMPLES-1);
v_target[ri]=VELOCITY_TARGET*(1-(1/0.15/0.15)*(r-0.15)*(r-0.15));
// sharp slope
if (r>0.15)
v_target[ri]=VELOCITY_TARGET; //
constant
}
#if 1==1
// ************ read
data from file
if (DISK &&
BULGE)
fh=fopen("sr_density_raw.txt", "r");
if (!DISK &&
BULGE)
fh=fopen("s_density_raw.txt", "r");
if (DISK &&
!BULGE)
fh=fopen("r_density_raw.txt", "r");
setvbuf(fh, NULL, _IOFBF, 5000);
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
{
fscanf(fh, "%i
%lf", &rv,
&datum); if
(ri<N_RADIAL_SAMPLES) rho_2d[ri]=datum;
fscanf(fh, "%lf",
&datum); if
(ri<N_RADIAL_SAMPLES) rho_3d[ri]=datum;
fscanf(fh, "%lf\n",
&datum); if
(ri<N_RADIAL_SAMPLES) v_target[ri]=datum;
}
fclose(fh);
#endif
rho_2d[N_RADIAL_SAMPLES-1]=0;
rho_3d[N_RADIAL_SAMPLES-1]=0;
v_target[0]=0;
if (!DISK) for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
rho_2d[ri]=0;
if (!BULGE) for
(ri=0; ri<N_RADIAL_SAMPLES; ri++)
rho_3d[ri]=0;
// Randomize
time(&t);
srand(t);
chi_new=chi();
// loop
for (loop=0;
loop<LOOPS; loop++)
{
if (DISK)
// Modify disk density at random
{
rv=(N_RADIAL_SAMPLES-1)*(rand()/((double)RAND_MAX+1));
chi_old=chi_new;
size=RANDOM_GAIN*(2*rand()/((double)RAND_MAX)-1.0)*chi_old/VELOCITY_TARGET*rho_2d[rv];
rho_2d[rv]+=size;
if
(rho_2d[rv]<0) rho_2d[rv]=size/10;
chi_new=chi();
if
(chi_new>chi_old)
{
rho_2d[rv]-=size;
chi_new=chi_old;
}
}
if (BULGE)
// Modify bulge density at random
{
rv=(N_RADIAL_SAMPLES-1)*(rand()/((double)RAND_MAX+1));
chi_old=chi_new;
size=RANDOM_GAIN*(2*rand()/((double)RAND_MAX)-1.0)*chi_old/VELOCITY_TARGET*rho_3d[rv];
rho_3d[rv]+=size;
if
(rho_3d[rv]<0) rho_3d[rv]=size/10;
chi_new=chi();
if
(chi_new>chi_old)
{
rho_3d[rv]-=size;
chi_new=chi_old;
}
}
// Calculate mass and angular momentum of disk
and bulge
mass=mass_ring(0)+mass_shell(0);
ang_mom2D=ang_mom3D=0;
for
(ri=1;
ri<N_RADIAL_SAMPLES; ri++)
{
mass_2D=mass_ring(ri);
mass_3D=mass_shell(ri);
mass+=mass_2D+mass_3D;
acc=accG_x(ri);
if (acc<0)
{
r=ri*RADIUS_STEP;
ang_mom2D+=mass_2D*sqrt(-acc*r)*r;
// mass * v * r
ang_mom3D+=2.0/3*mass_3D*sqrt(-acc*r)*r; // 2/3 mass * v * r
}
}
// Once in a while, transfer mass between bulge
and disk, smooth density functions
if ((loop%(N_RADIAL_SAMPLES))==1)
{
// ************* Compensate
for incorrect angular momentum
if (DISK
&& BULGE)
{
for (rv=0;
rv<N_RADIAL_SAMPLES-1; rv++)
{
r=rv*RADIUS_STEP;
mass_2D=mass_ring(rv);
mass_3D=mass_shell(rv);
c2=(ang_mom2D+ang_mom3D)/mass/GALAXY_RADIUS/1000-ANGULARM_TARGET;
if (c2>0)
c1=-mass_2D; // move mass to bulge
else
c1=mass_3D; // move mass to
disk
c2=fabs(N_RADIAL_SAMPLES*c2/VELOCITY_TARGET*1000)*MASS_TRANSFER;
if (c2>0.5)
c2=0.5;
c1*=c2;
rho_2d[rv]+=c1*rho_2d[rv]/mass_ring(rv);
rho_3d[rv]-=c1*rho_3d[rv]/mass_shell(rv);
}
}
// ************** Smoothing
of density functions
if
(DISK) // Smooth 2-D density,
keeping both ends fixed
for (ri=1;
ri<N_RADIAL_SAMPLES-2; ri++)
rho_2d[ri]=(SMOOTH_STR*rho_2d[ri]+rho_2d[ri-1]+rho_2d[ri+1])/(SMOOTH_STR+2);
if
(BULGE) // Smooth 3-D density;
keeping both ends fixed
for (ri=1;
ri<N_RADIAL_SAMPLES-2; ri++)
rho_3d[ri]=(SMOOTH_STR*rho_3d[ri]+rho_3d[ri-1]+rho_3d[ri+1])/(SMOOTH_STR+2);
chi_new=chi();
}
// *** print data
if
((loop%(N_RADIAL_SAMPLES/2))==0) printf("loop=%i:
chi=%lf [km/s]\n",
LOOPS-loop, chi_new/1000);
if
((loop%(N_RADIAL_SAMPLES*3))==0)
{
printf("\n
Radius Density Acceleration
Velocity\n");
printf("
[kpc] [ag/m^3]
[nm/s^2] [km/s]\n");
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
{
r=ri*RADIUS_STEP;
acc=accG_x(ri);
printf("%8.2lf
%8.3lf %9.4lf",
r/C_KPC, (rho_2d[ri]+rho_3d[ri])*1e21, acc*1e9);
if (acc>0) acc=0;
printf("
%8.3lf\n",
sqrt(-acc*r)/1000); // acc=m
v^2/r
ri+=N_RADIAL_SAMPLES/100;
if (r>0.1*GALAXY_RADIUS)
ri+=N_RADIAL_SAMPLES/10;
}
printf("Total
mass: %8.3lf [GM_sol] Angular momentum target: %8.3lf [M_galaxy
r_galaxy km/s]\n",
mass/C_MSOL/1e9, ANGULARM_TARGET);
printf("Angular momentum: %8.3lf (disk) + %8.3lf
(bulge) = %8.3lf
[M_galaxy r_galaxy km/s]\n",
ang_mom2D/mass/GALAXY_RADIUS/1000,
ang_mom3D/mass/GALAXY_RADIUS/1000,
(ang_mom2D+ang_mom3D)/mass/GALAXY_RADIUS/1000);
}
}// END: loop
// Save raw data
if (DISK &&
BULGE)
fh=fopen("sr_density_raw.txt", "w");
if (!DISK &&
BULGE)
fh=fopen("s_density_raw.txt", "w");
if (DISK &&
!BULGE)
fh=fopen("r_density_raw.txt", "w");
setvbuf(fh, NULL, _IOFBF, 5000);
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
fprintf(fh, "%6i
%10.4lg %10.4lg
%10.4lg\n", ri, rho_2d[ri], rho_3d[ri], v_target[ri]);
fclose(fh);
// Save data
if (DISK &&
BULGE)
fh=fopen("sr_density.txt", "w");
if (!DISK &&
BULGE)
fh=fopen("s_density.txt", "w");
if (DISK &&
!BULGE)
fh=fopen("r_density.txt", "w");
setvbuf(fh, NULL, _IOFBF, 5000);
fprintf(fh, "
radius
rho 2-D rho 3-D acceleration
velocity target
Intensity\n");
fprintf(fh, "
[kpc] [ag/m^3]
[ag/m^3] [nm/s^2]
[km/s] [km/s]
[M_sol/pc^2]\n");
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
{
r=ri*RADIUS_STEP;
acc=accG_x(ri);
if (acc>0)
acc=0;
fprintf(fh, "%8.2lf
%9.4lf %9.4lf %9.4lf
%8.3lf %8.3lf",
r/C_KPC,
rho_2d[ri]*1e21, rho_3d[ri]*1e21, acc*1e9, sqrt(-acc*r)/1000,
v_target[ri]/1000);
// find intensity seen from
above
intensity=rho_2d[ri]*DISK_THICKNESS+rho_3d[ri]*RADIUS_STEP;
for
(zi=1;
zi<N_RADIAL_SAMPLES; zi++)
{
rv=sqrt(r*r+zi*RADIUS_STEP*zi*RADIUS_STEP)/RADIUS_STEP;
// total radius
if
(rv<N_RADIAL_SAMPLES) intensity+=2*rho_3d[rv]*RADIUS_STEP;
}
fprintf(fh, "
%8.3lf\n", intensity/C_MSOL*(C_KPC/1000)*(C_KPC/1000));
}
mass_2D=mass_3D=0;
for (ri=0;
ri<N_RADIAL_SAMPLES; ri++)
{
mass_2D+=mass_ring(ri);
mass_3D+=mass_shell(ri);
}
fprintf(fh, "Mass disk %8.3lf
+ mass bulge %8.3lf =
Total mass %8.3lf [GM_sol]\n",
mass_2D/C_MSOL/1e9, mass_3D/C_MSOL/1e9, (mass_2D+mass_3D)/C_MSOL/1e9);
fprintf(fh, "Angular momentum
target: %8.3lf
[M_galaxy r_galaxy km/s]\n", ANGULARM_TARGET);
fprintf(fh, "Angular
momentum: %8.3lf (disk) +
%8.3lf (bulge) = %8.3lf [M_galaxy r_galaxy km/s]\n",
ang_mom2D/(mass_2D+mass_3D)/GALAXY_RADIUS/1000,
ang_mom3D/(mass_2D+mass_3D)/GALAXY_RADIUS/1000,
(ang_mom2D+ang_mom3D)/(mass_2D+mass_3D)/GALAXY_RADIUS/1000);
fclose(fh);
free(v_target);
free(rho_3d);
free(rho_2d);
return;
}