.. _custom_emitters_guide: Custom emitters distribution ============================ .. code:: ipython3 import numpy as np .. code:: ipython3 import jetset print('tested on jetset',jetset.__version__) .. parsed-literal:: tested on jetset 1.2.2 The user can bulid custom emitters distributions using the :class:`.EmittersDistribution` class. The following examples show how to implement it .. code:: ipython3 from jetset.jet_emitters import EmittersDistribution you need to define a function that describes your functional form (use numpy functions to make the code more performant) .. code:: ipython3 def distr_func_super_exp(gamma,gamma_cut,s,a): return np.power(gamma,-s)*np.exp(-(1/a)*(gamma/gamma_cut)**a) then you have to link the parmeters in your funtcion to a paramters of the ``EmittersDistribution`` class. .. note:: It is important to note that each parameter has the proper `par_type` string, and it is importat to properly set this `par_type` if you want to have good results using the phenomenological constraining. Have a look at the predefined model to understand how to set the `par_type` string. As a general rule: 1) the slope of the power law that starts from 'gmin' should be defined 'LE_spectral_slope' 2) any slope of the power law above the break should be defined as 'HE_spectral_slope' 3) values of gamma defining the transition from a power law trend toward a curved trend (including cut-off) should be defined as 'turn-over-energy' 4) curvature in log-parabolic trends and/or super-exp exponent, should be defined as 'spectral_curvature' 4) the spectral_type='user_defined' *will not work* for phenomenological constraining: .. code:: ipython3 print('allowed_sepctral types for phenomenological constraining',EmittersDistribution.spectral_types_obs_constrain()) .. parsed-literal:: allowed_sepctral types for phenomenological constraining ['bkn', 'plc', 'lp', 'lppl', 'pl', 'lpep', 'array'] .. code:: ipython3 n_e_super_exp=EmittersDistribution('super_exp',spectral_type='user_defined') n_e_super_exp.add_par('gamma_cut',par_type='turn-over-energy',val=50000.,vmin=1., vmax=None, unit='lorentz-factor') n_e_super_exp.add_par('s',par_type='LE_spectral_slope',val=2.3,vmin=-10., vmax=10, unit='') n_e_super_exp.add_par('a',par_type='spectral_curvature',val=1.8,vmin=0., vmax=100., unit='') now you have to link your defined functional form to the ``EmittersDistribution`` class. .. code:: ipython3 n_e_super_exp.set_distr_func(distr_func_super_exp) parameters can be easily set .. code:: ipython3 n_e_super_exp.parameters.s.val=.4 n_e_super_exp.parameters.s.val=2.0 n_e_super_exp.parameters.gamma_cut.val=1E5 n_e_super_exp.normalize=True n_e_super_exp.parameters.gmax.val=1E6 .. code:: ipython3 n_e_super_exp.parameters.show_pars() .. raw:: html Table length=6
namepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
gminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
gmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
Nemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
gamma_cutturn-over-energylorentz-factor*1.000000e+051.000000e+00--FalseFalse
sLE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
aspectral_curvature1.800000e+000.000000e+001.000000e+02FalseFalse
.. code:: ipython3 p=n_e_super_exp.plot() .. image:: custom_emitters_files/custom_emitters_17_0.png .. code:: ipython3 p=n_e_super_exp.plot(energy_unit='eV') .. image:: custom_emitters_files/custom_emitters_18_0.png here we define a bkn power-law .. code:: ipython3 def distr_func_bkn(gamma_break,gamma,s1,s2): return np.power(gamma,-s1)*(1.+(gamma/gamma_break))**(-(s2-s1)) n_e_bkn=EmittersDistribution('bkn',spectral_type='bkn') n_e_bkn.add_par('gamma_break',par_type='turn-over-energy',val=1E3,vmin=1., vmax=None, unit='lorentz-factor') n_e_bkn.add_par('s1',par_type='LE_spectral_slope',val=2.5,vmin=-10., vmax=10, unit='') n_e_bkn.add_par('s2',par_type='HE_spectral_slope',val=3.2,vmin=-10., vmax=10, unit='') n_e_bkn.set_distr_func(distr_func_bkn) n_e_bkn.parameters.show_pars() n_e_bkn.parameters.s1.val=2.0 n_e_bkn.parameters.s2.val=3.5 p=n_e_bkn.plot() .. raw:: html Table length=6
namepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
gminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
gmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
Nemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
gamma_breakturn-over-energylorentz-factor*1.000000e+031.000000e+00--FalseFalse
s1LE_spectral_slope2.500000e+00-1.000000e+011.000000e+01FalseFalse
s2HE_spectral_slope3.200000e+00-1.000000e+011.000000e+01FalseFalse
.. image:: custom_emitters_files/custom_emitters_20_1.png Passing the custom distribution to the Jet class ------------------------------------------------ .. code:: ipython3 from jetset.jet_model import Jet my_jet=Jet(electron_distribution=n_e_bkn) **now the ``n_e_bkn`` will be deep copyed, so changes applied to the one passed to the model will not affect the original one** .. code:: ipython3 my_jet.parameters.N.val=5E4 my_jet.show_model() my_jet.IC_nu_size=100 my_jet.eval() .. parsed-literal:: -------------------------------------------------------------------------------- model description: -------------------------------------------------------------------------------- type: Jet name: jet_leptonic electrons distribution: type: bkn gamma energy grid size: 201 gmin grid : 2.000000e+00 gmax grid : 1.000000e+06 normalization: False log-values: False ratio of cold protons to relativistic electrons: 1.000000e-01 radiative fields: seed photons grid size: 100 IC emission grid size: 100 source emissivity lower bound : 1.000000e-120 spectral components: name:Sum, state: on name:Sync, state: self-abs name:SSC, state: on external fields transformation method: blob SED info: nu grid size jetkernel: 1000 nu size: 500 nu mix (Hz): 1.000000e+06 nu max (Hz): 1.000000e+30 flux plot lower bound : 1.000000e-30 -------------------------------------------------------------------------------- .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm35.000000e+040.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*1.000000e+031.000000e+00--FalseFalse
jet_leptonics1LE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonics2HE_spectral_slope3.500000e+00-1.000000e+011.000000e+01FalseFalse
.. parsed-literal:: -------------------------------------------------------------------------------- Since as default, the ``Nomralization`` is false, let’s check the actual number density of particles and conpare it to the parameter ``N`` .. code:: ipython3 print('N_particle=',my_jet.emitters_distribution.eval_N(),'N parameter=',my_jet.parameters.N.val) .. parsed-literal:: N_particle= 24608.46344775512 N parameter= 50000.0 .. note:: N_particle is different from N, because the distribution is not normalized .. code:: ipython3 my_jet.eval() .. code:: ipython3 p=my_jet.plot_model() p.setlim(y_min=1E-16,y_max=1E-13) .. image:: custom_emitters_files/custom_emitters_29_0.png Now we switch on the normalization for the emetters distribtuion, and we keep all the parameters unchanged, including N .. code:: ipython3 my_jet.Norm_distr = True my_jet.parameters.N.val=5E4 my_jet.show_model() my_jet.IC_nu_size=100 my_jet.eval() .. parsed-literal:: -------------------------------------------------------------------------------- model description: -------------------------------------------------------------------------------- type: Jet name: jet_leptonic electrons distribution: type: bkn gamma energy grid size: 201 gmin grid : 2.000000e+00 gmax grid : 1.000000e+06 normalization: True log-values: False ratio of cold protons to relativistic electrons: 1.000000e-01 radiative fields: seed photons grid size: 100 IC emission grid size: 100 source emissivity lower bound : 1.000000e-120 spectral components: name:Sum, state: on name:Sync, state: self-abs name:SSC, state: on external fields transformation method: blob SED info: nu grid size jetkernel: 1000 nu size: 500 nu mix (Hz): 1.000000e+06 nu max (Hz): 1.000000e+30 flux plot lower bound : 1.000000e-30 -------------------------------------------------------------------------------- .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm35.000000e+040.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*1.000000e+031.000000e+00--FalseFalse
jet_leptonics1LE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonics2HE_spectral_slope3.500000e+00-1.000000e+011.000000e+01FalseFalse
.. parsed-literal:: -------------------------------------------------------------------------------- and we check again the actual number density of particles and conpare it to the parameter N .. code:: ipython3 print('N_particle=',my_jet.emitters_distribution.eval_N(),'N parameter=',my_jet.parameters.N.val) .. parsed-literal:: N_particle= 50000.0 N parameter= 50000.0 .. note:: N_particle and N now are the same, because the distribution is normalized .. code:: ipython3 p=my_jet.plot_model() p.setlim(y_min=1E-16,y_max=1E-13) .. image:: custom_emitters_files/custom_emitters_35_0.png Building a distribution from an external array ---------------------------------------------- Here we just build two arrays, but you can pass any ``n_gamma`` and ``gamma`` array wit the same size, and with ``gamma>1`` and ``n_gamma>0`` .. code:: ipython3 from jetset.jet_emitters import EmittersArrayDistribution import numpy as np # gamma array gamma = np.logspace(1, 8, 500) # gamma array this is n(\gamma) in 1/cm^3/gamma n_gamma = gamma ** -2 * 1E-5 * np.exp(-gamma / 1E5) N1 = np.trapz(n_gamma, gamma) n_distr = EmittersArrayDistribution(name='array_distr', emitters_type='electrons', gamma_array=gamma, n_gamma_array=n_gamma,normalize=False) N2 = np.trapz(n_distr._array_n_gamma, n_distr._array_gamma) ``N1`` and ``N2`` are used only for the purpose of checking, you can skip them .. code:: ipython3 p=n_distr.plot() .. image:: custom_emitters_files/custom_emitters_40_0.png .. code:: ipython3 my_jet = Jet(emitters_distribution=n_distr, verbose=False) my_jet.show_model() .. parsed-literal:: -------------------------------------------------------------------------------- model description: -------------------------------------------------------------------------------- type: Jet name: jet_leptonic electrons distribution: type: array_distr gamma energy grid size: 501 gmin grid : 1.000000e+01 gmax grid : 1.000000e+08 normalization: False log-values: False ratio of cold protons to relativistic electrons: 1.000000e-01 radiative fields: seed photons grid size: 100 IC emission grid size: 100 source emissivity lower bound : 1.000000e-120 spectral components: name:Sum, state: on name:Sync, state: self-abs name:SSC, state: on external fields transformation method: blob SED info: nu grid size jetkernel: 1000 nu size: 500 nu mix (Hz): 1.000000e+06 nu max (Hz): 1.000000e+30 flux plot lower bound : 1.000000e-30 -------------------------------------------------------------------------------- .. raw:: html Table length=9
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*1.000000e+011.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+081.000000e+001.000000e+15FalseFalse
jet_leptonicNscaling_factor1.000000e+000.000000e+00--FalseFalse
.. parsed-literal:: -------------------------------------------------------------------------------- you can also skip the next cell, it is just to check .. code:: ipython3 N3 = np.trapz(my_jet.emitters_distribution.n_gamma_e, my_jet.emitters_distribution.gamma_e) np.testing.assert_allclose(N1, N2, rtol=1E-5) np.testing.assert_allclose(N1, N3, rtol=1E-2) np.testing.assert_allclose(N1, my_jet.emitters_distribution.eval_N(), rtol=1E-2) ``N`` will act as a scaling factor for the array when normalization is set to ``False`` .. code:: ipython3 my_jet.parameters.N.val=1E9 print('this is the actual number of emitters dendisty %2.2f'%my_jet.emitters_distribution.eval_N(),'this the scaling factor',my_jet.parameters.N.val) .. parsed-literal:: this is the actual number of emitters dendisty 999.56 this the scaling factor 1000000000.0 .. code:: ipython3 my_jet.eval() p=my_jet.plot_model() .. image:: custom_emitters_files/custom_emitters_46_0.png you can still normalize the distribution .. code:: ipython3 my_jet.Norm_distr = True my_jet.parameters.N.val=2000 print('this is the actaul number of emitters dendisty %2.2f'%my_jet.emitters_distribution.eval_N(),'this the scaling factor',my_jet.parameters.N.val) .. parsed-literal:: this is the actaul number of emitters dendisty 2000.00 this the scaling factor 2000 .. code:: ipython3 my_jet.eval() p=my_jet.plot_model() .. image:: custom_emitters_files/custom_emitters_49_0.png .. code:: ipython3 my_jet.save_model('test_jet_custom_emitters_array.pkl') new_jet = Jet.load_model('test_jet_custom_emitters_array.pkl') .. raw:: html Table length=9
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*1.000000e+011.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+081.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm32.000000e+030.000000e+00--FalseFalse
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
.. code:: ipython3 new_jet.eval() p=new_jet.plot_model() .. image:: custom_emitters_files/custom_emitters_51_0.png