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:

**Code:**

// 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

As you can see

**particle.secondary_source_function** is the same for both 2D and 3D. Function

**pbarqqFF** 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.

The only difference between 2D and 3D is the following:

**Code:**

n_spatial_dimensions = 3

z_min = -10.0 min z

z_max = +10.0 max z

dz = 1.0 delta z

x_min = -10.0 min x

x_max = +10.0 max x

dx = 1.0 delta x

y_min = -10.0 min y

y_max = +10.0 max y

dy = 1.0 delta y

prop_x = 1 1=propagate in x (3D)

prop_y = 1 1=propagate in y (3D)

prop_z = 1 1=propagate in z (2D, 3D)

prop_p = 1 1=propagate in p (2D, 3D)

p_Ekin_grid = Ekin p||Ekin alignment

p_min = 100 min momentum (MV)

p_max = 40000 max momentum (MV)

p_factor = 1.5 momentum factor

Ekin_min = 1.0e1 min kinetic energy per nucleon (MeV)

Ekin_max = 1.0e6 max kinetic energy per nucleon (MeV)

Ekin_factor = 1.5 kinetic energy per nucleon factor

Am I missing something?