Increased clarity in input/output

I have been playing around with Synchrad for a little while now trying to apply it to a non-standard use-case, where there is no obvious "fundamental frequency" like in an undulator or betatron radiation. The tutorials are based around these schemes, and while it is possible for me to extract some methodology from these tutorials, it would be of much help to be able to see what outputs are available from the different functions.

Maybe it is because I am new to Python, but it is not at all obvious to me how to see what is output by the different functions in the code. In Spyder, I can check the main Synchrad module, but it doesn't tell me anything about what information is available after calc = SynchRad( calc_input ) has been run. Moreover, if I then try to look into calc in the variable browser, it says that "the variable is not picklable". I would therefore very much appreciate that something is added at least to the documentation stating what output is available, as well as what units it is given in.

Something which I am currently trying to do is to simply calculate the total photon number / emitted energy, but in for example the undulator tutorial there are some unclear inputs that are (at least to me) not obvious, such as the lambda0_um and phot_num variables. phot_num evidently changes the unit and/or scaling of the output to something I don't understand. Or rather: I guess that calc.get_energy( phot_num=True, normalize_to_weights=True ) should give me the emitted average number of photons per electron, but if I set phot_num=False, I get a very large value which seems to be in neither J nor eV.

I am also interested in getting the radiated power as a function of propagation distance/time (which I guess I can set with the variable nSnaps in the calculate_spectrum?) as well as the physical near-field source size (to be able to calculate the peak brilliance), but again it is not obvious to me how to handle the outputs.

Any clarification here would be much appreciated! Best regards, Jonas

Hi @jonasbjorklundsvensson

Thank you for your interest in Synchrad. The code's documentation is indeed being a WIP, and in particular the utility part, but i'll try to add it asap. Meanwhile, I can answer some of your questions here.

As you may see in tutorials/Introduction, the actual expressions which we calculate are dimensionless, so all it defined by the units in which you measure space/time in your tracks. For a non-normalized case, e.g. space units being meters (as well as timeStep=c*dt), the frequency axis in the grid is normalized to a frequency of a photon with lambda_0 = 1 m, so if you wish to define a grid with photon energies from 1 eV to 1 keV, you can do it, for example, as:

calc_input = {
    "grid": [
        (1/energy_1m_eV , 1e3/energy_1m_eV),  
        (0, theta_max),
        (0.0, 2 * np.pi),
        (1024, 32, 32),
    ],  
}

where the energy of a photon with 1m wavelength energy_1m_eV = 2*pi*c*hbar/e can be imported from synchrad.utils.

In this non-normalized case, in the utilities methods you may pass the lambda0_um=1e6 which corresponds to the 1m unit length being measured in microns. In that case get_energy should return the total energy in J and if phot_num=True is selected the returned value should be a number of photons. Note, that if you wish to obtain the total photon number / emitted energy you have to ensure that all the emission is contained in the 3D region (omega, theta, phi) defined by your grid -- so it is recommended to plot the 2D (omega, theta) maps for different phi and be sure noting is truncated.

Concerning structure of calc and the calculation result. The calc is an object with many different data, and after calc.calculate_spectrum has been executed, the radiation is contained in the dictionary calc.Data['radiation'], however, I do not recommend working with it directly. Instead the utils methods should be used, and I would also recommend to save the simulation to the file, by giving file_spectrum='./spectrum.h5' argument to calculate_spectrum -- then you'll always re-initilze the executed simulation with calc = SynchRad(file_spectrum='./spectrum.h5') and work with the data.

The main utilities:

get_full_spectrum: returns a total calculated spectrum as an array defined on the grid
get_energy_spectrum: returns a 1D spectrum overs frequencies integrated over all angles
get_energy: further integrates the spectrum over the frequencies to provide total emission energy
get_spot_cartesian gets an angular spot (total spectrum integrated over frequencies)

You are correct that using nSnaps gets the propagation divided into a number of intervals. In order to work with such data you will only need to use the iteration= argument for the utility methods, and the corresponding track iterations are stored as calc.snap_iterations.

Let me know if this helps, and do not hesitate to ask further -- it will help me to decide how to put these explanations into the docs/tutorials.

If you have difficulties with python, i can post some examples from the analysis I do myself. For an example, you can use the following script to plot the 2D (omega, theta) maps with radiation truncated for emission energies <50 eV:

calc = SynchRad(file_spectrum='./spectrum.h5')
spect_ax_keV = calc.get_spectral_axis()  * energy_1m_eV * 1e-3

k = calc.Args['omega'][:,None,None] * energy_1m_eV
spect_filter = (k>50)

extent = [spect_ax_keV.min(), 
          spect_ax_keV.max(), 
          -calc.Args['theta'].max()* 1e3,
          calc.Args['theta'].max()* 1e3 ]

spect_val = calc.get_full_spectrum(spect_filter=spect_filter)
i_th = spect_val.shape[2] // 4                                   ## looking at phi=pi/2 slice
i_th_diag = i_th + spect_val.shape[2] // 2
spect_val = np.c_[spect_val[:, :,  i_th_diag][:, ::-1], 
                  spect_val[:, 1:, i_th     ]]
plt.figure()
plt.imshow( spect_val.T,
            extent = extent,
            cmap=plt.cm.nipy_spectral,
            origin='lower', 
            aspect='auto' )

plt.colorbar()

Hi @hightower8083 !

Sorry for the slow reply... have been swamped with stuff lately.

Thanks for the input, I'm going to try your suggestions out, but as this project is a bit on the backburner it might be a while before I make more progress.

Some additional questions: the spectrum is given in photons / 0.1%BW; I want to calculate the critical energy from usual definitions (i.e., half the energy contained in the spectrum up until the critical energy) and I'm not entirely sure how I should approach this. Getting the energy spectrum should simply be multiplying the photon number in each energy bin by the photon energy in the corresponding bin, i.e. dW/dE = E*dN/dE, but I'm not sure how to integrate the spectrum when the energy is not absolute but given in relative units (0.1%BW).

Related to this; if I get the angle-integrated spectrum by for example spectrum = calc.get_energy_spectrum( normalize_to_weights=True ), should I scale this by the same factor as I scale my energy axis, for example by a factor of 1e3 if I rescale from eV to keV? It seems to me that since the BW-bit in the spectrum is relative, I should not do that, right?

Best regards, Jonas

Hi @jonasbjorklundsvensson

Yes, the units may seem a bit confusing. If you take Jackson's formula for the far-field, and divide it by \hbar you get it dimensionless:

$\cfrac{\mathrm{d}\mathcal{W}}{\mathrm{d}\Omega\, \mathrm{d}\varepsilon_{ph}} = \frac{\alpha}{4\pi^2}\left| \int_{-\infty}^{\infty} \mathrm{d}t\,\cfrac{\bf n\!\times\! (\bf n\!-\!\beta)\!\times\!\dot{\beta}}{(1-{\bf n\,\beta})^2}\;\mathrm{e}^{i \omega (t-\mathbf{n}\mathbf{r}_e/c )}\right|^2\,,$

and the code actually calculates the last term |..|^2. The factor \alpha / 4 \pi^2 is added by the utilities methods, and when you call get_.. methods without lambda0_um (i.e. lambda0_um=None), you get just this number.

When you use some value for lambda0_um it adds another factor that fixes Joule for the units of the integral over the photon frequencies, which are normalized in the code to the \omega_0=2\pi c/\lambda_0. This is only meaningful for get_energy and get_spot_cartesian, called with phot_num=False.

When you consider "0.1%b.w." binning, it means that you take d\epsilon = 1e-3 * \epsilon and dW = \epsilon * dN, and this dN will be your number of photons per 0.1%b.w. So you can easily get it, by call ingget_.. methods without lambda0_um, and multiplyinh the result by 1e-3 (it comes from 0.1%). You can see how it is done in this tutorial, where the line spect_val *= 1e-3 appears (cell In[8]). Note that the data is defined on the regular grid and the binning notion only concerns the units.

About the critical energy. I'm not sure what's its usual definition -- more often we fit it with classical synchrotron radiation spectra. If you want to get the median energy which divides full emitted energy equally you may try looking for it like this

spectrum = calc.get_energy_spectrum( normalize_to_weights=True )
spect_ax_keV = calc.get_spectral_axis()  * energy_1m_eV * 1e-3

spect_val_int = spect_val.cumsum()
energy_crit = spect_ax_keV[ ( spect_val_int <= 0.5 * spect_val_int[-1] ).sum() ]

# estimate relative error
Energy_above = calc.get_energy(normalize_to_weights=True, spect_filter=(spect_ax_keV>energy_crit)[:,None, None])
Energy_below = calc.get_energy(normalize_to_weights=True, spect_filter=(spect_ax_keV<=energy_crit)[:,None, None])
print( 2* np.abs(Energy_above-Energy_below)/ (Energy_above+Energy_below))

Let me know if it helps

hightower8083 / synchrad

Increased clarity in input/output #20