Verify the results of stopping power calculations

[grzanka]

Check from where the differences comes from: core libamtrack library or from some GUI features.

[grzanka]

Data inspection

PSTAR from https://physics.nist.gov/cgi-bin/Star/ap_table.pl

[grzanka]

Looking at https://github.com/libamtrack/library/blob/master/src/AT_StoppingPowerDataPSTAR.c
we have:

		  { 3.000E+01 , 1.876E+01 },
		  { 3.500E+01 , 1.656E+01 },
		  { 4.000E+01 , 1.488E+01 },
		  { 4.500E+01 , 1.354E+01 },
		  { 5.000E+01 , 1.245E+01 },

which is consistent with web GUI plot data and PSTAR webpage database.

The only difference may be attributed to interpolation method used to calculate stopping power at 32 MeV and other intermediate points.

[grzanka]

ICRU data, extracted from https://github.com/libamtrack/library/blob/master/src/AT_StoppingPowerDataICRU.c:

30 | 40 | 50
-- | -- | --
0.01876 | 0.01488 | 0.01245

Seems to be consistent with PSTAR data and GUI plot data.

[grzanka]

Interpolation and data extraction:

ICRU:

mass_stopping_power_MeV_cm2_g[i] = AT_get_interpolated_y_from_input_table(
		AT_stopping_power_ICRU_table[material_no-1].energy_and_stopping_power[0],
		AT_stopping_power_ICRU_table[material_no-1].energy_and_stopping_power[Z[i]],
		AT_stopping_power_ICRU_table[material_no-1].number_of_data_points, E_MeV_u[i]);

// TODO: remove after proper scaling of data
if(material_no == Water_Liquid){
    mass_stopping_power_MeV_cm2_g[i] *= 1000;
}

PSTAR:

mass_stopping_power_MeV_cm2_g[i] = AT_get_interpolated_y_from_input_2d_table(
					source_for_given_material->energy_and_stopping_power,
					n_data,
					E_MeV_u[i]);
if( particle_no[i] != PARTICLE_PROTON_NUMBER){
   double Zeff_ion    =  AT_effective_charge_from_E_MeV_u_single(E_MeV_u[i], particle_no[i]);
   double Zeff_proton =  AT_effective_charge_from_E_MeV_u_single(E_MeV_u[i], PARTICLE_PROTON_NUMBER);
   mass_stopping_power_MeV_cm2_g[i] *= gsl_pow_2(Zeff_ion / Zeff_proton);
}

[grzanka]

For energies between 30 and 50 MeV, libamtrack has PSTAR data for every 5 MeV (5 data points) while ICRU data for every 10 MeV (3 data points).

The stopping power for PSTAR for 35 MeV is taken from data points for PSTAR (1.656E+01 MeV/cm2/g) while for ICRU it is being interpolated from 30 and 40 MeV data points.

Note that PSTAR value for 35 MeV (1.656E+01 ) is different from average value of 30 and 40 MeV ( 1.876E+01 + 1.488E+01 / 2 = 1.682E+01). The difference is ~1.57%.

[grzanka]

PSTAR webpage uses cubic spline interpolation, while libamtrack uses linear one:

Python 3.7.0 (default, Oct  9 2018, 10:31:47) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from scipy import interpolate
>>> x = [30.,35.,40.,45.]
>>> y = [1.876E+01,  1.656E+01 , 1.488E+01,  1.354E+01]
>>> fs = interpolate.interp1d(x,y,kind='cubic')
>>> fl = interpolate.interp1d(x,y,kind='linear')
>>> fs(36.)
array(16.18816)
>>> fl(36.)
array(16.224)

Libamtrack stopping power calculated for 36 MeV protons is 1.6224 which is the close to the value obtained with linear interpolation, while PSTAR value (as calculated by https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html) is 16.18816 which is closer to cubic spline interpolation. The difference is ~0.22%.

[grzanka]

As described in 3.5. Methods Used To Evaluate Stopping Powers on https://physics.nist.gov/PhysRefData/Star/Text/programs.html a cubic spline is used.

[grzanka]

The difference in stopping power prediction is not related to any GUI issue, hence closing in favor of libamtrack/library#25 and libamtrack/library#24