ehsteve / roentgen

A Python package for the quantitative analysis of the interaction of x-rays with matter.
12 stars 5 forks source link

Providing continuum x-ray emission #20

Closed ehsteve closed 1 year ago

DanRyanIrish commented 4 years ago

Here is a plot comparing the output of thermal_bremsstrahlung() and the f_vth and brem_49 SSW routines. The roentgen curve has been scaled for easier comparison and the T, EM, n_e, and z input parameters for each routine are listed in the plot title and legend. image

hayesla commented 4 years ago

it looks like as @ehsteve mentioned that its difference in gaunt factor calculations - see figure below for comparison of gaunt factor calculation for roentgen mewe_gaunt and IDL acgaunt. The one here mewe_gaunt is the approximation where the IDL acgaunt uses the more full calculation -

mewe_comparison

ehsteve commented 4 years ago

@hayesla i'm confused because these two gaunt factors seem fairly close and so cannot explain the factor of ten difference I see in @DanRyanIrish at low energies but maybe this is a normalization issue? Dan seems to have normalized to match at higher energies so perhaps they deviate at high energies?

ehsteve commented 4 years ago

Could we get the plots on the same x range? say 1 to 200 keV?

hayesla commented 4 years ago

Could we get the plots on the same x range? say 1 to 200 keV?

oh actually I was playing around with it, and if you scale the output by energy (i.e. keV) and set the density to 0.01 it matches pretty well with the offsets more in line with the gaunt factors. I'm a bit confused @ehsteve - what is the density units supposed to be in thermal_bremsstrahlung?

both_comparsion_fixed

ehsteve commented 4 years ago

That's great! Thanks @hayesla. The density value is n_i * n_e. For solar we usually integrate that out into the emission measure. I just wanted to capture the dependencies but did not have the time to calculate the correct constant (which is why all the units are messed up as well). Alright, to me this is good enough. Certainly for the attenuator paper where we only need the shape to be correct. What do you think @DanRyanIrish?

ehsteve commented 4 years ago

@DanRyanIrish could you compare the XrayTubeSpectrum provided here with the spectra provided by the website you found http://solutioinsilico.com/medical-physics/applications/tasmip-app.php?ans=1&n_sp=1&k1=30&m1=1&t1=100&k2=&m2=1&t2=&k3=&m3=1&t3=

DanRyanIrish commented 4 years ago

XraytubeSpectrum fails when energy (keV) given rather than wavelength (nm). May be an issue with CutoffWavelength

In [18]: roentgen_data = XraytubeEmission(web_data["energy"], 100*u.kV)                                                                                          
---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
<ipython-input-18-f72e14b5789e> in <module>
----> 1 roentgen_data = XraytubeEmission(web_data["energy"], 100*u.kV)

~/miniconda3/envs/roentgen-dev/lib/python3.7/site-packages/astropy/units/decorators.py in wrapper(*func_args, **func_kwargs)
    232             # Call the original function with any equivalencies in force.
    233             with add_enabled_equivalencies(self.equivalencies):
--> 234                 return_ = wrapped_function(*func_args, **func_kwargs)
    235             if wrapped_signature.return_annotation not in (inspect.Signature.empty, None):
    236                 return return_.to(wrapped_signature.return_annotation)

~/sunpy_dev/roentgen/roentgen/continuum/emission.py in XraytubeEmission(wavelength, voltage)
     85     emission : The unscaled emission, provided without units
     86     """
---> 87     f = (wavelength / CutoffWavelength(voltage) - 1) * 1 / wavelength ** 2
     88     # remove unphysical negative values
     89     f[f < 0.0] = 0.0

~/miniconda3/envs/roentgen-dev/lib/python3.7/site-packages/astropy/units/quantity.py in __array_ufunc__(self, function, method, *inputs, **kwargs)
    455         # consistent units between two inputs (e.g., in np.add) --
    456         # and the unit of the result (or tuple of units for nout > 1).
--> 457         converters, unit = converters_and_unit(function, method, *inputs)
    458 
    459         out = kwargs.get('out', None)

~/miniconda3/envs/roentgen-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/converters.py in converters_and_unit(function, method, *args)
    187                             "argument is not a quantity (unless the "
    188                             "latter is all zero/infinity/nan)"
--> 189                             .format(function.__name__))
    190             except TypeError:
    191                 # _can_have_arbitrary_unit failed: arg could not be compared

UnitConversionError: Can only apply 'subtract' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)
DanRyanIrish commented 4 years ago

Here is the comparison of the roentgen and web-tool x-ray tube spectra for a tube voltage of 100 kV image

ehsteve commented 1 year ago

Other packages are to provide this. Out of scope.