After experiencing the same problem I traced the error back to line 44 in D_xx.cc. The problem is that regardless of the value of n_spatial_dimensions the array Dxx.d2[ir][iz].s[ip] is filled.
If the user specified n_spatial_dimensions = 3, Dxx.d3[ix][iy][iz] will be created rather than Dxx.d2[ir][iz]. But D_xx.cc tries to fill d2[...] which obviously can not work. What I did was insert if (n_spatial_dimensions == 3) ... and if (n_spatial_dimensions == 2) statements into the code.
Here's my updated version of the file:
Code:
//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****|
// * D_xx.cc * galprop package * 02/13/2003
//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****|
//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****|
// Wave damping formalism is described in:
//
// Ptuskin, V.S., et al. 2006, ApJ 642, 902
// Ptuskin, V.S., et al. 2005, Adv. Space Res. 35, 162
//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****|
using namespace std;//AWS20050624
#include<cstdio>
#include<cstdlib>
#include"galprop_classes.h"
#include"galprop.h"
int iprotons,ir,ix,iy,iz,ip;
int damping_min_ip; // IMOS20060330
//this is to avoid problems of using Galprop class members in static function "fu" IMOS20060322
Particle *protons;
double damping_p0;
int n_spatial_dimensions, diff_reacc;
int Galprop::D_xx(Particle &particle,int iprotons_,int ir_,int ix_,int iy_,int iz_,int ip_)
{
iprotons=iprotons_; ir=ir_; ix=ix_; iy=iy_; iz=iz_; ip=ip_;
double L_cm, Lp_cm, tmp;
// integration parameters
double a=particle.rigidity[ip],ai;
//this is to avoid problems of using Galprop class members in static function "fu" IMOS20060322
protons=&gcr[iprotons];
damping_p0=galdef.damping_p0;
n_spatial_dimensions=galdef.n_spatial_dimensions;
diff_reacc=galdef.diff_reacc;
// STANDARD DIFFUSION COEFFICIENT (galdef.diff_reacc =0, 1, 2, -1==beta^3 Dxx)
if(galdef.diff_reacc < 3)
{
if (n_spatial_dimensions == 2) {
particle.Dxx.d2[ir][iz].s[ip] = particle.beta[ip] *galdef.D0_xx;
if(galdef.diff_reacc<0) particle.Dxx.d2[ir][iz].s[ip] = pow(particle.beta[ip],galdef.diff_reacc) *galdef.D0_xx;
if(particle.rigidity[ip]<galdef>=galdef.D_rigid_br)
particle.Dxx.d2[ir][iz].s[ip]*= pow(particle.rigidity[ip]/galdef.D_rigid_br, galdef.D_g_2);
return 0;
}
if (n_spatial_dimensions == 3) {
particle.Dxx.d3[ix][iy][iz].s[ip] = particle.beta[ip] *galdef.D0_xx;
if(galdef.diff_reacc<0) particle.Dxx.d3[ix][iy][iz].s[ip] = pow(particle.beta[ip],galdef.diff_reacc) *galdef.D0_xx;
if(particle.rigidity[ip]<galdef>=galdef.D_rigid_br)
particle.Dxx.d3[ix][iy][iz].s[ip]*= pow(particle.rigidity[ip]/galdef.D_rigid_br, galdef.D_g_2);
return 0;
}
}
// WAVE DAMPING (see Ptuskin et al. astro-ph/0301420)
if(ip==particle.n_pgrid-1) damping_min_ip=0; // IMOS20060330
L_cm = galdef.damping_max_path_L; // max free path
if(ip<damping_min_ip) // IMOS20060330
{
if (n_spatial_dimensions == 2) {
particle.Dxx.d2[ir][iz].s[ip] = particle.beta[ip] *C*L_cm/3.;
// printout at the solar system position
// if(iz==particle.n_zgrid/2+1 && ir==9) cout<<D_xx>>>> "<<particle.rigidity[ip]<<" "<<particle.Ekin[ip]<<" "<<particle.beta[ip] *C*L_cm/3.<<" "<<-ai<<endl;
return 0;
}
if (n_spatial_dimensions == 3) {
particle.Dxx.d3[ix][iy][iz].s[ip] = particle.beta[ip] *C*L_cm/3.;
// printout at the solar system position
// if(iz==particle.n_zgrid/2+1 && ir==9) cout<<D_xx>>>> "<<particle.rigidity[ip]<<" "<<particle.Ekin[ip]<<" "<<particle.beta[ip] *C*L_cm/3.<<" "<<-ai<<endl;
return 0;
}
}
Lp_cm = 3./C *galdef.D0_xx*pow(particle.rigidity[ip]/galdef.D_rigid_br,galdef.D_g_1);
for(int i=1;i<gcr> galdef.damping_p0)
{
damping_p0 = gcr[iprotons].rigidity[i]; // re-definition of galdef.damping_p0
break;
}
/*******************TEST
ir=0; iz=0;
for(int i=0; i<particle.n_pgrid; i++)
{
gcr[iprotons].cr_density.d2[0][0].s[i] = pow(gcr[iprotons].p[i],-2.);
cout<<p>"<<damping_p0<<" "<<a<<" "<<ai<<endl>0.) L_cm = Lp_cm/pow(tmp, 2);
}
// Kraichnan diffusion with wave damping ##
if(galdef.diff_reacc==12)
{
tmp = 1. -galdef.damping_const_G*(-ai); // /pow(particle.Z, 3./2.)
if(tmp>0.) L_cm = Lp_cm/tmp;
}
if (L_cm>galdef.damping_max_path_L && gcr[iprotons].rigidity[ip]<1.e4) // IMOS20050907
{
L_cm = galdef.damping_max_path_L;
damping_min_ip=ip;
}
if (n_spatial_dimensions == 2 ) particle.Dxx.d2[ir][iz].s[ip] = particle.beta[ip] *C*L_cm/3.;
if (n_spatial_dimensions == 3 ) particle.Dxx.d3[ix][iy][iz].s[ip] = particle.beta[ip] *C*L_cm/3.;
if(iz==particle.n_zgrid/2+1 && ir==9)
cout<<D_xx>>>> "<<particle.rigidity[ip]<<" "<<particle.Ekin[ip]
<<" "<<particle.beta[ip] *C*L_cm/3.<<" "<<-ai<<endl;
return 0;
}
//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****|
double Galprop::fu(double x)
{
int i,n,m;
double int_psi=0., y;
// cout<<ir>>>> x,n= "<<x<<" "<<n<<" gcr[iprotons].p[n]= "<<gcr[iprotons].p[n]<<endl;
// integration over rigidity (= momentum for protons)
if(n_spatial_dimensions==2)
{
// fit each interval with power-law and integrate analytically
for(int_psi=0., i=n-1; i<m>cr_density.d2[ir] [iz].s[i]/protons->cr_density.d2[ir] [iz].s[i+1])
/log(protons-> p[i]/protons-> p[i+1]);
// integrate (psi/p) dp
if(i>n-1) int_psi +=protons->cr_density.d2[ir][iz].s[i ] // fit norm.
*pow(protons-> p[i ],-y)/y
*(pow(protons-> p[i+1], y)
-pow(protons-> p[i ], y));
else int_psi +=protons->cr_density.d2[ir][iz].s[i ] // fit norm.
*pow(protons-> p[i ],-y)/y
*(pow(protons-> p[i+1], y)
-pow(x, y));
}
}
if(n_spatial_dimensions==3)
{
// fit each interval with power-law and integrate analytically
for(int_psi=0., i=n-1; i<m>cr_density.d3[ix][iy][iz].s[i]/protons->cr_density.d3[ix][iy][iz].s[i+1])
/log(protons-> p[i]/protons-> p[i+1]);
// integrate (psi/p) dp
if(i>n-1) int_psi +=protons->cr_density.d3[ix][iy][iz].s[i ] // fit norm.
*pow(protons-> p[i ],-y)/y
*(pow(protons-> p[i+1], y)
-pow(protons-> p[i ], y));
else int_psi +=protons->cr_density.d3[ix][iy][iz].s[i ] // fit norm.
*pow(protons-> p[i ],-y)/y
*(pow(protons-> p[i+1], y)
-pow(x, y));
}
}
//if(ir==0 && iz==80 && ip==40) cout<<" integral = "<<n<<" "<<y<<" "<<int_psi<<" "<<protons>p[n]<<" "<<x<<" "<<damping_p0 <<endl; //exit(1);
if(diff_reacc==11) return ( pow(x,2./3.)*int_psi ); // Kolmogorov
if(diff_reacc==12) return (sqrt(x) *int_psi ); // Kraichnan
return 0; // in case of an error return 0.
}
Regards
Bastian