Bug 39 - injection spectrum parameter beta_rig
injection spectrum parameter beta_rig
Status: NEW
Product: GALPROP Code
Classification: Unclassified
Component: Code
54 (in WebRun)
PC Linux
: --- enhancement
Assigned To: Troy Porter
Depends on:
Blocks:
  Show dependency treegraph
 
Reported: 2010-11-05 06:58 PDT by Andy Strong
Modified: 2010-11-05 07:50 PDT (History)
3 users (show)

See Also:
Web Browser: ---


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Andy Strong 2010-11-05 06:58:47 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... ?
Comment 1 Andy Strong 2010-11-05 07:00:02 PDT
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.
Comment 2 Andy Strong 2010-11-05 07:01:31 PDT
original forum message:
http://galprop.stanford.edu/forum/viewtopic.php?f=19&t=112
Comment 3 Troy Porter 2010-11-05 07:05:42 PDT
(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?
Comment 4 Andy Strong 2010-11-05 07:11:18 PDT
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.
Comment 5 Andy Strong 2010-11-05 07:14:14 PDT
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.
Comment 6 Andy Strong 2010-11-05 07:19:48 PDT
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).
Comment 7 Andy Strong 2010-11-05 07:22:23 PDT
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.
Comment 8 Andy Strong 2010-11-05 07:31:58 PDT
see e.g. sec 2.2 of 
http://arxiv.org/abs/1011.0989
Comment 9 Andy Strong 2010-11-05 07:49:22 PDT
(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).
Comment 10 Andy Strong 2010-11-05 07:50:16 PDT
yes we need the more general form, as well as documenting where the current form comes from.