Bugzilla for GALPROP – Bug 39
injection spectrum parameter beta_rig
Last modified: 2010-11-05 07:50:16 PDT
From forum user tomassetti, reposted here. the Galdef option "inj_spectrum_type = beta_rig" is considered in line #64 of "create_transport_array.cc". Here's the code: //if(strcmp(galdef.inj_spectrum_type,"beta_rig")==0) spec_shape*=particle.beta[ip]; // IMOS20000615 if(strcmp(galdef.inj_spectrum_type,"beta_rig")==0) spec_shape/=sqrt(1.+pow(2.e3/particle.rigidity[ip],2)); // IMOS20011210 If I don't make a misunderstanding, the first line (commented) corresponds to a real rigidity power law x beta (i.e. what the user expects from the choice "beta_rig"). Second line (adopted) gives the same rigidity shape for all charged CRs (at injection). But it is "beta_rig" only assuming M/Z = 2GeV for all particles (included protons, on which M/Z~1 GeV). It seems to me a crude approximation. So I wonder what was wrong with the first line... ?
thanks for finding this. I don't know what is going on here offhand, unfortunately the change was not documented in the code. It also needs to be put into the Explanatory Supplement Section 5.5.1. Normally we use only rigidity without beta, so this was probably not tested much.
original forum message: http://galprop.stanford.edu/forum/viewtopic.php?f=19&t=112
(In reply to comment #2) > original forum message: > http://galprop.stanford.edu/forum/viewtopic.php?f=19&t=112 I think this code is most likely an AWS or IMOS `special', since only you guys persist in using cstdio functions. What was the idea behind the code in the first place?
it bears the stamp of IMOS as you see (at least he leaves his mark when making a change). cstdio was all that existed when galprop c++ was invented. .... but now it is anyway a string (undocumented who did that, can guess) anyway it is important to fix this since the beta form is of interest from theory.
NB the current version r585 has if(galdef.inj_spectrum_type == "beta_rig") spec_shape/=sqrt(1.+pow(2.e3/particle.rigidity[ip],2)); // IMOS20011210 which is what I was referring to above.
svn blame create_transport_arrays.cc 548 gudlaugu if(galdef.inj_spectrum_type == "beta_rig") spec_shape/=sqrt(1.+pow(2.e3/particle.rigidity[ip],2)); // IMOS20011210 (which is OK except the formula is from IMOS, the latest change from Gulli).
The idea behind the code: some theories for shock acceleration predict a beta-dependence in addition to the rigidity powerlaw. It is quite commonly used (not by us). Anyway we should get it right.
see e.g. sec 2.2 of http://arxiv.org/abs/1011.0989
(copied from forum, from tomassetti) I did not think to it as a bug. Such spectral shape, only rigidity-dependent, is physically plausible as well as the broken one, often used. However the user should be aware of this, and it can be interested in a real "beta x rig" injection. I also suggest to consider powers of beta (beta^W x rig), with W free parameter to be used instead of nuc_g_1 (that can be fixed =nuc_g_2).
yes we need the more general form, as well as documenting where the current form comes from.
