!
! Survival probability calculation for varying electron number density along neutrino path
!
! Based 'flavour basis' equations, section II, in
!
! Three flavour neutrino oscillations in matter: Flavour diagonal potentials, the adiabatic basis, and the CP phase
! James P. Kneller, Gail C. McLaughlin
! PHYSICAL REVIEW D 80, 053002 (2009)
!
! conventions
!
! mass eigenstates (vacuum) index imass = 1, 2 , 3
!
! flavour eigenstates (e,mu,tau) index ialpha = 1, 2, 3
!
! pseudo code outline
!
! acquire mixing parameters : mixing angles (3), phase angle (1) and known differences of squared mass (2)
!
! set flag to indicate normal / inverted mass hierachy and/or set angles omega(alpha) (to be used in equation 12)
!
! set values for arbritrary phases alpha, beta and epsilon (CP)
!
! acquire path information: distance travelled (n segments) and electron number density (n segments)
!
! calculate all variables which depend on mixing parameters: sin12, sin23, sin13, cos12, cos13, cos23, msq12, msq13, msq23
!
! for all path length segments
!
! calculate elements of Hamiltonian matrix HF in flavour base (equations 6a to 6 f)
!
! calculate values of q(alpha) and r(alpha) (equations 10a to 10f)
!
! calculate values of TF, QF, RF (equations 8a to 8c or alternatively equations 9a to 10f)
!
! calculate eigenvalues eigen_value(alpha) (equation 12)
!
! find three eigenvectors eigen_vector(alpha) of Hamiltonian HF for given eigenvalues (use LINPACK, MKL, GSL, other ? )
!
! Note that by definition : HF * eigen_vector(alpha) = eigen_value(alpha) * eigen_vector(alpha)
! Solution of Schroedinger equation : sum_over_alpha: ( a(alpha) * exp( -j * eigen_value(alpha) * time / hbar ) * eigen_vector(alpha) ) = sum_over_alpha: state(alpha,time)
! Boundary condition at start time: sum_over_alpha: ( a(alpha) * eigen_vector(alpha) ) = sum_over_alpha: previous_state(alpha) ( in practice a 3 x 3 system of equations )
!
! find coefficients a(alpha) by numerically solving ( 3 x 3 ) linear system (question: are the coefficients real numbers?)
!
! calculate state at end of step: sum_over_alpha: state(alpha) = sum_over_alpha: a(alpha) * exp( -j * eigen_value(alpha) * time / hbar ) * eigen_vector(alpha)
!
! calculate survival probabilities P(beta) = ||**2 where alpha is state at end of last step
!
! enter next points P(beta,t) or P(beta,L) with L=c*t into a plot (histogram? points?)
!
! check: 0 < P(beta,L) < 1 and sum_over_beta: P(beta,L) = 1
!
! if this is the end of the last step then exit loop
!
! generate plot, display, save plot
!
! Note(s):
!
! alternative algorythm: start with matrix HF and calculate both eigenvalues and eigenvectors numerically? Would this be any slower?
!
! TF, QF, RF or else eigen_value(alpha) plus eigen_vector(alpha) (total of 12 values) could be tabulated, and then recovered from a spline function
!
! test case: results should reproduce the vacuum case (zero matter density)
!
! test case: results should reproduce the constant matter density case (see journal paper)
!
! test case: results should reproduce the exponential matter density case (see journal paper)
!
!
! Questions:
!
! Implement as mixed language code? Fortran for complex numbers, vectors, matrices and if needed 128 bit precision? C++ to the use CERN root libraries for plots?
!
! Implement twice as two independent pieces of code as cross check? Is this needed as there are several test cases?
!
! Can precise results be approximated by average matter density over all path length? Or a low order polynominal mass density?
!
! Under which conditions can several neutrino path segments be combined into a single average density segment?
!
! What is the 'oscillation length scale' over which variations in abundance become significant?
program neutrino_flux_calculation
print*," Neutrino flux calculation"
end program