Skip to content

Commit

Permalink
Remove unimplemented analytical NSI matrix calculation from nsi_param…
Browse files Browse the repository at this point in the history 
…s.py (#786)
  • Loading branch information
@finnmayhew
finnmayhew authored Jul 30, 2024
1 parent 6ad2b2b commit 081606c
Showing 1 changed file with 1 addition and 154 deletions.
155 changes: 1 addition & 154 deletions pisa/stages/osc/nsi_params.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -384,133 +384,6 @@ def eps_matrix(self):

return nsi_eps

@property
def eps_matrix_analytical(self):
"""Effective NSI coupling matrix calculated analytically."""
# Analytical relations. These are wrong right now! #FIXME
nsi_eps = np.zeros((3, 3, 2), dtype=FTYPE)

sp12 = np.sin(self.phi12)
sp13 = np.sin(self.phi13)
sp23 = np.sin(self.phi23)
cp12 = np.sqrt(1. - sp12**2)
cp13 = np.sqrt(1. - sp13**2)
cp23 = np.sqrt(1. - sp23**2)

sdnsi = np.sin(self.deltansi)
cdnsi = np.cos(self.deltansi)

# eps_ee - eps_mumu (real)
nsi_eps[0, 0, 0] = (
self.eps_scale * cp13**2 * (cp12**2 - sp12**2) +
self.eps_prime * (
(cp12**2 - sp12**2) * (sp13**2 * sp23**2 - cp23**2) -
4 * cp12 * sp12 * sp13 * cp23 * sp23 * cdnsi
)
) - 1
nsi_eps[0, 0, 1] = 0.
# eps_emu (complex)
nsi_eps[0, 1, 0] = (
self.eps_scale * cp12 * sp12 * cp13**2 * np.cos(self.alpha1 - self.alpha2) +
self.eps_prime * (
(
cp12 * sp12 * (sp13**2 * sp23**2 - cp23**2) +
sp13 * cp23 * sp23 * cdnsi * (cp12**2 - sp12**2)
) * np.cos(self.alpha1 - self.alpha2) -
(
sp13 * cp23 * sp23 * sdnsi
) * np.sin(self.alpha1 - self.alpha2)
)
)
nsi_eps[0, 1, 1] = (
self.eps_scale * cp12 * sp12 * cp13**2 * np.sin(self.alpha1 - self.alpha2) +
self.eps_prime * (
(
cp12 * sp12 * (sp13**2 * sp23**2 - cp23**2) +
sp13 * cp23* sp23 * cdnsi * (cp12**2 - sp12**2)
) * np.sin(self.alpha1 - self.alpha2) +
(
sp13 * cp23 * sp23 * sdnsi
) * np.cos(self.alpha1 - self.alpha2)
)
)
# eps_etau (complex)
nsi_eps[0, 2, 0] = (
-self.eps_scale * cp12 * sp13 * cp13 * np.cos(2 * self.alpha1 + self.alpha2) +
self.eps_prime * (
(
cp13 * sp23 * (cp12 * sp13 * sp23 - sp12 * cp23 * cdnsi)
) * np.cos(2 * self.alpha1 + self.alpha2) -
(
cp13 * sp12 * cp23 * sp23 * sdnsi
) * np.sin(2 * self.alpha1 + self.alpha2)
)
)
nsi_eps[0, 2, 1] = (
-self.eps_scale * cp12* sp13 * cp13 * np.sin(2 * self.alpha1 + self.alpha2) +
self.eps_prime * (
(
cp13 * sp23 * (cp12 * sp13 * sp23 - sp12 * cp23 * cdnsi)
) * np.sin(2 * self.alpha1 + self.alpha2) +
(
cp13 * sp23 * sp12 * cp23 * sdnsi
) * np.cos(2 * self.alpha1 + self.alpha2)
)
)
# eps_emu* (complex)
nsi_eps[1, 0, 0] = nsi_eps[0, 1, 0]
nsi_eps[1, 0, 1] = -nsi_eps[0, 1, 1]
# eps_etau* (complex)
nsi_eps[2, 0, 0] = nsi_eps[0, 2, 0]
nsi_eps[2, 0, 1] = -nsi_eps[0, 2, 1]
# eps_mumu - eps_mumu (0 by definition)
nsi_eps[1, 1, 0] = 0.
nsi_eps[1, 1, 1] = 0.
# eps_mutau (complex)
nsi_eps[1, 2, 0] = (
-self.eps_scale * sp12 * cp13 * sp13 * np.cos(self.alpha1 + 2 * self.alpha2) +
self.eps_prime * (
(
cp13 * sp23 * (sp12 * sp13 * sp23 + cp12 * cp23 * cdnsi)
) * np.cos(self.alpha1 + 2 * self.alpha2) +
(
cp12 * cp13 * cp23 * sp23 * sdnsi
) * np.sin(self.alpha1 + 2 * self.alpha2)
)
)
nsi_eps[1, 2, 1] = (
-self.eps_scale * sp12 * cp13 * sp13 * np.sin(self.alpha1 + 2 * self.alpha2) +
self.eps_prime * (
(
-cp12 * cp13 * cp23 * sp23 * sdnsi
) * np.cos(self.alpha1 + 2 * self.alpha2) +
(
cp13 * sp23 * (sp12 * sp13 * sp23 + cp12 * cp23 * cdnsi)
) * np.sin(self.alpha1 + 2 * self.alpha2)
)
)
# eps_mutau* (complex)
nsi_eps[2, 1, 0] = nsi_eps[1, 2, 0]
nsi_eps[2, 1, 1] = -nsi_eps[1, 2, 1]
# eps_tautau - eps_mumu (real)
nsi_eps[2, 2, 0] = (
self.eps_scale * (sp13**2 - cp13**2 * sp12**2) +
self.eps_prime *(
sp23**2 * (cp13**2 - sp12**2 * sp13**2) -
2 * cp12 * sp12 * sp13 * cp23 * sp23 * cdnsi -
cp12**2 * cp23**2
)
)
nsi_eps[2, 2, 1] = 0.

# make this into a complex 2d array
nsi_eps = nsi_eps[:, :, 0] + nsi_eps[:, :, 1] * 1.j
# make sure this is a valid Hermitian potential matrix
# before returning anything
assert np.allclose(nsi_eps, nsi_eps.conj().T, **ALLCLOSE_KW)

return nsi_eps

def test_nsi_params():
"""Unit tests for subclasses of `NSIParams`."""
# TODO: these have to be extended
Expand Down Expand Up @@ -563,35 +436,9 @@ def test_nsi_parameterization():
nsi_params.alpha2 = alpha2
nsi_params.deltansi = deltansi

logging.trace('Checking agreement between numerical & analytical NSI matrix...')

eps_mat_numerical = nsi_params.eps_matrix
eps_mat_analytical = nsi_params.eps_matrix_analytical

try:
close = np.isclose(eps_mat_numerical, eps_mat_analytical, **ALLCLOSE_KW)
if not np.all(close):
logging.debug(
"Numerical NSI matrix:\n%s",
np.array2string(eps_mat_numerical, **ARY2STR_KW)
)
logging.debug(
"Analytical expansion (by hand):\n%s",
np.array2string(eps_mat_analytical, **ARY2STR_KW)
)
raise ValueError(
'Evaluating analytical expressions for NSI matrix elements'
' does not give agreement with numerical calculation!'
' Elementwise agreement:\n%s'
% close
)
except ValueError as err:
logging.warning(
"%s\nThis is expected."
" Going ahead with numerical calculation for now.", err
)

logging.trace('Now checking agreement with sympy calculation...')
logging.trace('Checking agreement with sympy calculation...')

eps_mat_sympy = nsi_sympy_mat_mult(
eps_scale_val=eps_scale,
Expand Down

0 comments on commit 081606c

Please sign in to comment.