Now I'm pretty sure that there is something going on with GALPROP that I'm not aware about.
I'm running two simulations with the same source function: 2D and 3D. Here is the code of gen_DM_source.cc:
// assign the source function (2D)
if(galaxy.n_spatial_dimensions==2)
{
for(int ir=0; ir<gcr[0].n_rgrid; ir++)
{
for(int iz=0; iz<gcr[0].n_zgrid; iz++)
{
for(int ip=0; ip<particle.n_pgrid; ip++)
{
/////////////////////////////////////////////////////////////////////
if (strcmp(particle.name,"DM_antiprotons") == 0 && galdef.DM_int1 == 777)
{
if (particle.Etot[ip] * 1.e-3 <= DMmass)
{
double source = 0.5 * 3 * pow(10, -26) * pow(800 / DMmass, 2) * pow(DM_profile_av(galaxy.r[ir], galaxy.z[iz], galaxy.dr, galaxy.dz, dzz) / DMmass, 2) *
pbarqqFF(particle.Etot[ip] * 0.001 - 0.938, t_a, aSize, t_b, t_c, t_d) / 1000;
particle.secondary_source_function.d2[ir][iz].s[ip] += (SoL/(4 * Pi)) * source;
}
continue;
}
/////////////////////////////////////////////////////////////////////
} // ip
} // iz
} // ir
} // particle.n_spatial_dimensions==2
// assign the source function (3D)
if(galaxy.n_spatial_dimensions==3)
{
for(int ix=0; ix<gcr[0].n_xgrid; ix++)
{
for(int iy=0; iy<gcr[0].n_ygrid; iy++)
{
for(int iz=0; iz<gcr[0].n_zgrid; iz++)
{
for(int ip=0; ip<particle.n_pgrid; ip++)
{
/////////////////////////////////////////////////////////////////////
if (strcmp(particle.name,"DM_antiprotons") == 0 && galdef.DM_int1 == 777)
{
if (particle.Etot[ip] * 1.e-3 <= DMmass)
{
double source = 0.5 * 3 * pow(10, -26) * pow(800 / DMmass, 2) * pow(DM_profile_av(galaxy.x[ix], galaxy.y[ix], galaxy.z[iz], galaxy.dx, galaxy.dy, galaxy.dz, dzz) / DMmass, 2) *
pbarqqFF(particle.Etot[ip] * 0.001 - 0.938, t_a, aSize, t_b, t_c, t_d) / 1000;
particle.secondary_source_function.d3[ix][iy][iz].s[ip] = (SoL/(4 * Pi)) * source;
}
continue;
}
/////////////////////////////////////////////////////////////////////
} //ip
} // iz
} // iy
} // ix
} // particle.n_spatial_dimensions==3
is the same for both 2D and 3D. Function
is my fragmentaion function of quarks into antiprotons. It is also the same in both cases.
But while the 2D results are acceptable, the 3D ones look crazy. I've deliberately drawn the proton's flux which is the same in both cases (as expected). I think something is going on with 3D source functions, something I'm not accounting for.
I've also attached a galdef file for 2D. The spatial grids for 3D are less detailed because it takes more time to compute in 3 dimensions than in 2, but that wouldn't change the results by 20 orders of magnitude.
FAQ
Search
Register
Login
