CRPropa / CRPropa3

CRPropa is a public astrophysical simulation framework for propagating extraterrestrial ultra-high energy particles. https://crpropa.github.io/CRPropa3/
https://crpropa.desy.de
GNU General Public License v3.0
65 stars 66 forks source link

Regarding use of the TabularPhotonField class #418

Closed PrantikS closed 1 year ago

PrantikS commented 1 year ago

Hi, I am trying to create a custom photon field using the class TabularPhotonField(). I have created the required files (My_Field_photonEnergy.txt, My_Field_photonDensity.txt) as per these instructions and placed them in the directory share/crpropa/Scaling.

When I tried to use this custom field in the following way: TabularPhotonField(My_Field, False) It gives this error "NameError: name 'My_Field' is not defined." Do I need to define this somewhere?

I also tried calling the inbuilt field IRB_Stecker05, the same way: TabularPhotonField(IRB_Stecker05, True). However, it also shows this error:

TypeError: Wrong number or type of arguments for overloaded function 'new_TabularPhotonField'. Possible C/C++ prototypes are: crpropa::TabularPhotonField::TabularPhotonField(std::string const,bool const) crpropa::TabularPhotonField::TabularPhotonField(std::string const)

Kindly, let me know how to fix this issue!

lukasmerten commented 1 year ago

We are currently working on this. Examples for the custom photon fields will be updated shortly.

To make custom photon fields fully available in all interaction in CRPropa, you will also need to make use of the CRPropa-data repository. Since tabulated data as interaction rates are calculated with function from there. To streamline the whole process also the CRPropa-data repository will be updated soon (see the already opended PR.

This my answer of how to implement custom photon fields will be shortly outdated, I have to ask to wait until the already developed changes are merged. In case this does not solve your problem, you can reopen this issue.

PrantikS commented 1 year ago

@lukasmerten Thanks a lot!

I think, to add a new photon field, I need to first add it to the photonField.py file. Is it correct? Also, how do I modify an existing photon field? For example, if I want to modify the following photon field (URB_Nitu21) defined in the file photonField.py, what is the correct way? Do I have to change the parameters: p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10 ? If so, could you help me to understand what are these parameters and how to obtain them?

`class URB_Nitu21: """ Universal Radio Background from Nitu et al. 2021. Reference: I. C. Nitu, H. T. J. Bevings, J. D. Bray, A. M. M. Scaife Astroparticle Physics 126 (2021) 102532. https://arxiv.org/abs/2004.13596 """ name = 'URB_Nitu21' info = 'URB_Nitu21' redshift = None

def getDensity(self, eps, z=0):
    """
    Comoving spectral number density dn/deps [1/m^3/J] at given photon energy eps [J]
    """
    p0 = -1.9847e1
    p1 = -2.9857e-1
    p2 = -2.6984e-1
    p3 = 9.5394e-2
    p4 = -4.9059e-2
    p5 = 4.4297e-3
    p6 = 7.6038e-3
    p7 = -1.9690e-3
    p8 = -2.2573e-4
    p9 = 1.1762e-4
    p10 = -9.9443e-6
    p = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10]

    eps = np.r_[eps]
    nu = eps / h
    I = 0.
    for k in range(len(p)):
        I += (p[k] * np.power(np.log10(nu / 1e6), k))
    I = 10. ** I
    I = 4 * np.pi / (h * c0) * (I / eps)

    I[eps < self.getEmin()] = 0.
    I[eps > self.getEmax()] = 0.

    return I

def getEmin(self, z=0):
    """Minimum effective photon energy in [J]"""
    return 1e3 * Hz * h

def getEmax(self, z=0):
    """Maximum effective photon energy in [J]"""
    return 1e12 * Hz * h`
JulienDoerner commented 1 year ago

hey @PrantikS, The parameters of the Nitu field are defined in the mentioned paper (see table 2).

To create your own photon field, you can have a look at the file calc_all.py in the mentioned PR. There is a function createPhotonTargetInteractions(fields: list) which creates the data products needed for CRPropa. For this you need a class where the photon density and some other information of your field is given. The easiest way of creating this class, is using inheritance from the bass class in photonField.py.

PrantikS commented 1 year ago

Hi @JulienDoerner, Thanks a lot for the help. I have created the fields that I needed.