25 #include <gsl/gsl_complex.h> 26 #include <gsl/gsl_complex_math.h> 27 #include <gsl/gsl_matrix.h> 28 #include <gsl/gsl_eigen.h> 29 #include <gsl/gsl_math.h> 30 #include <gsl/gsl_const_num.h> 31 #include <gsl/gsl_const_mksa.h> 32 #include <gsl/gsl_const_cgsm.h> 41 fEval = gsl_vector_alloc(numNu);
42 fEvec = gsl_matrix_complex_alloc(numNu, numNu);
43 H_GSL = gsl_matrix_complex_alloc(numNu, numNu);
44 W_GSL = gsl_eigen_hermv_alloc(numNu);
49 const int numNu = p.
fEval->size;
50 fEval = gsl_vector_alloc(numNu);
51 fEvec = gsl_matrix_complex_alloc(numNu, numNu);
52 H_GSL = gsl_matrix_complex_alloc(numNu, numNu);
53 W_GSL = gsl_eigen_hermv_alloc(numNu);
64 gsl_matrix_complex_free(
H_GSL);
65 gsl_eigen_hermv_free(
W_GSL);
66 gsl_matrix_complex_free(
fEvec);
67 gsl_vector_free(
fEval);
73 gsl_eigen_hermv_workspace*
W_GSL;
79 static std::complex<double>
zero(0.0,0.0);
80 static std::complex<double>
one (1.0,0.0);
83 static const double kKm2eV = 5.06773103202e+09;
84 static const double kK2 = 4.62711492217e-09;
88 static const double kGf = 1.166371E-5;
152 cout <<
"First argument should be smaller than second argument" <<
endl;
153 cout <<
"Setting reverse order (Theta" << j << i <<
"). " <<
endl;
159 cout <<
"Theta" << i << j <<
" not valid for " <<
fNumNus;
160 cout <<
" neutrinos. Doing nothing." <<
endl;
178 cout <<
"First argument should be smaller than second argument" <<
endl;
179 cout <<
"Setting reverse order (Delta" << j << i <<
"). " <<
endl;
185 cout <<
"Delta" << i << j <<
" not valid for " <<
fNumNus;
186 cout <<
" neutrinos. Doing nothing." <<
endl;
190 cout <<
"Rotation " << i << j <<
" is real. Doing nothing." <<
endl;
208 cout <<
"Dm" << j <<
"1 not valid for " <<
fNumNus;
209 cout <<
" neutrinos. Doing nothing." <<
endl;
229 if(
fTheta[i][j]==0)
return;
244 for(
int k=0; k<
i; k++){
245 fHmsBufferC =
fHms[k][
i];
247 fHms[k][
i] *= fCosBuffer;
248 fHms[k][
i] +=
fHms[k][
j] * fSinBuffer * conj(fExpBuffer);
250 fHms[k][
j] *= fCosBuffer;
251 fHms[k][
j] -= fHmsBufferC * fSinBuffer * fExpBuffer;
255 for(
int k=i+1; k<
j; k++){
256 fHmsBufferC =
fHms[k][
j];
258 fHms[k][
j] *= fCosBuffer;
259 fHms[k][
j] -= conj(
fHms[i][k]) * fSinBuffer * fExpBuffer;
261 fHms[
i][k] *= fCosBuffer;
262 fHms[
i][k] += fSinBuffer * fExpBuffer * conj(fHmsBufferC);
267 fHmsBufferD = real(
fHms[j][j]);
269 fHms[
i][
i] *= fCosBuffer * fCosBuffer;
270 fHms[
i][
i] += 2 * fSinBuffer * fCosBuffer * real(
fHms[i][j] * conj(fExpBuffer));
273 fHms[
j][
j] *= fCosBuffer * fCosBuffer;
274 fHms[
j][
j] += fSinBuffer * fHmsBufferC * fSinBuffer;
275 fHms[
j][
j] -= 2 * fSinBuffer * fCosBuffer * real(
fHms[i][j] * conj(fExpBuffer));
277 fHms[
i][
j] -= 2 * fSinBuffer * real(
fHms[i][j] * conj(fExpBuffer)) * fSinBuffer * fExpBuffer;
278 fHms[
i][
j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC) * fExpBuffer;
284 for(
int k=i+1; k<
j; k++){
285 fHms[k][
j] = -conj(
fHms[i][k]) * fSinBuffer * fExpBuffer;
287 fHms[
i][k] *= fCosBuffer;
291 fHmsBufferD = real(
fHms[i][i]);
293 fHms[
i][
j] = fSinBuffer * fCosBuffer * (
fHms[
j][
j] - fHmsBufferD) * fExpBuffer;
295 fHms[
i][
i] *= fCosBuffer * fCosBuffer;
296 fHms[
i][
i] += fSinBuffer * fHms[
j][
j] * fSinBuffer;
298 fHms[
j][
j] *= fCosBuffer * fCosBuffer;
299 fHms[
j][
j] += fSinBuffer * fHmsBufferD * fSinBuffer;
308 for(
int k=0; k<
i; k++){
309 fHmsBufferC =
fHms[k][
i];
311 fHms[k][
i] *= fCosBuffer;
314 fHms[k][
j] *= fCosBuffer;
315 fHms[k][
j] -= fHmsBufferC * fSinBuffer;
320 fHmsBufferD = real(
fHms[j][j]);
322 fHms[
i][
i] *= fCosBuffer * fCosBuffer;
323 fHms[
i][
i] += 2 * fSinBuffer * fCosBuffer * real(
fHms[i][j]);
326 fHms[
j][
j] *= fCosBuffer * fCosBuffer;
327 fHms[
j][
j] += fSinBuffer * fHmsBufferC * fSinBuffer;
328 fHms[
j][
j] -= 2 * fSinBuffer * fCosBuffer * real(
fHms[i][j]);
330 fHms[
i][
j] -= 2 * fSinBuffer * real(
fHms[i][j]) * fSinBuffer;
331 fHms[
i][
j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC);
339 fHms[
i][
i] = fSinBuffer * fHms[
j][
j] * fSinBuffer;
341 fHms[
j][
j] *= fCosBuffer * fCosBuffer;
371 for(
int i=0;
i<
j;
i++){
375 for(
int i=0;
i<
j;
i++){
405 double kr2GNe =
kK2*M_SQRT2*
kGf*Ne;
411 *gsl_matrix_complex_ptr(
d->
H_GSL,
i,
i) = gsl_complex_rect(real(buf),0);
414 if(anti>0) *gsl_matrix_complex_ptr(
d->
H_GSL,
j,
i) = gsl_complex_rect(real(buf),-imag(buf));
415 else *gsl_matrix_complex_ptr(
d->
H_GSL,
j,
i) = gsl_complex_rect(real(buf), imag(buf));
419 if(anti>0) *gsl_matrix_complex_ptr(
d->
H_GSL,
i,
i) = gsl_complex_add_real(gsl_matrix_complex_get(
d->
H_GSL,
i,
i) , kr2GNe/2);
420 else *gsl_matrix_complex_ptr(
d->
H_GSL,
i,
i) = gsl_complex_add_real(gsl_matrix_complex_get(
d->
H_GSL,
i,
i) , -kr2GNe/2);;
424 if(anti>0) *gsl_matrix_complex_ptr(
d->
H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(
d->
H_GSL,0,0) , kr2GNe);
425 else *gsl_matrix_complex_ptr(
d->
H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(
d->
H_GSL,0,0) , -kr2GNe);
450 gsl_complex buf = gsl_matrix_complex_get(
d->
fEvec,
j,
i);
460 gsl_complex buf = gsl_matrix_complex_get(
d->
fEvec,
i,
j);
475 const std::list<double>& Ne,
478 if (L.size()!=Ne.size()) abort();
479 std::list<double>::const_iterator Li (L.begin());
480 std::list<double>::const_iterator Lend(L.end());
481 std::list<double>::const_iterator Ni (Ne.begin());
482 for (; Li!=Lend; ++Li, ++Ni) {
std::vector< double > fDm
m^2_i - m^2_1 in vacuum
virtual void SetDelta(int i, int j, double delta)
virtual void SetAngle(int i, int j, double th)
static const double kK2
mole/GeV^2/cm^3 to eV
virtual double P(int flv) const
virtual void PropMatter(double L, double E, double Ne, int anti=1)
virtual void SetDm(int i, double dm)
gsl_matrix_complex * fEvec
double fCachedNe
Cached electron density.
gsl_matrix_complex * H_GSL
static const double kGeV2eV
GeV to eV.
virtual void RotateH(int i, int j)
Rotate the Hamiltonian by theta_ij and delta_ij.
static constexpr double L
gsl_eigen_hermv_workspace * W_GSL
std::vector< complex > fNuState
The neutrino current state.
fvar< T > exp(const fvar< T > &x)
std::complex< double > complex
virtual void SetStdPars()
Set standard 3-flavor parameters.
virtual void ResetToFlavour(int flv=1)
Oscillation probability calculators.
int fCachedAnti
Cached nu/nubar selector.
bool fBuiltHms
Tag to avoid rebuilding Hms.
std::vector< std::vector< double > > fTheta
theta[i][j] mixing angle
assert(nhit_max >=nhit_nbins)
virtual void SolveHam(double E, double Ne, int anti)
double fCachedE
Cached neutrino energy.
std::vector< std::vector< complex > > fHms
matrix H*2E in eV^2
static std::complex< double > one(1.0, 0.0)
std::vector< std::vector< double > > fDelta
delta[i][j] CP violating phase
virtual void InitializeVectors()
static const double kKm2eV
km to eV^-1
static std::complex< double > zero(0.0, 0.0)