Add Bortfeld single BP model

[grzanka]

See libamtrack/library#22

[grzanka]

See https://github.com/khpeek/spread-out-Bragg-peak

[grzanka]

https://github.com/afg1/SOBPtools/blob/566835a0790166c045ff26555cc5c713055d8dd3/README.md

[grzanka]

See https://github.com/jasqs/PTClib/blob/master/lib/analyseDose.py

[grzanka]

Also https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/63405/versions/3/previews/BraggPeakAnalysis/analyseBP.m/index.html

http://iopscience.iop.org/article/10.1088/0031-9155/56/24/003/meta

https://aapm.onlinelibrary.wiley.com/doi/abs/10.1118/1.598116

[grzanka]

[grzanka]

Some prototype code to start with:

import mpmath as mp  # its part of sympy, offers high-precision floating-point arithmetic

# lets use again models.py file and display the implementation of Wilkens model
from models import WilkensLET

def fun(energy_MeV, Phi_0_cm2, sigma_energy_MeV, z_cm, eps = 0.1):
  
  rho_g_cm3 = 1.0
  
  beta_cm = 0.012  # 1/cm
  
  gamma = 0.6

  # range
  range_cm = ERSCalc.range_cm(energy_MeV)

  # range straggling of monoenergetical protons, see Appendix [1]
  sigma_mono_cm = 0.012 * range_cm ** 0.935

  # range equivalent of energy straggline, equation (A2) in [1]
  sigma_r_cm = sigma_energy_MeV
  sigma_r_cm *= ERSCalc.alpha_cm_MeV
  sigma_r_cm *= ERSCalc.p
  sigma_r_cm *= (energy_MeV ** (ERSCalc.p - 1.0))

  # total sigma, equation (A3) in [1]
  sigma_cm = (sigma_mono_cm ** 2 + sigma_r_cm ** 2) ** 0.5
  
  A = sigma_cm ** (1.0 / ERSCalc.p)
  A *= mp.gamma(1.0 / ERSCalc.p)
  A /= mp.sqrt( 2.0 * np.pi)
  A /= rho_g_cm3
  A /= ERSCalc.p
  A /= (ERSCalc.alpha_cm_MeV ** (1.0 / ERSCalc.p))
  A /= (1.0 + beta_cm * range_cm)
  
  # zeta variable introduced for equation (A8) in [1]
  zeta = (range_cm - z_cm) / sigma_cm
  
  B = beta_cm / ERSCalc.p + gamma * beta_cm + eps / range_cm
  
  result = Phi_0_cm2 * A
  bracket_part1 = WilkensLET.parabolic_integral(1.0 / ERSCalc.p, -zeta) / sigma_cm
  bracket_part2 = B * WilkensLET.parabolic_integral(1.0 + (1.0 / ERSCalc.p), -zeta)
  result *= ( bracket_part1 + bracket_part2 )
  
  result *= 1.6021766e-10
  
  return result

[grzanka]

Basic implementation is in 9469f2f

TODO:

• better description
• test directory with validation notebooks