From c4de6fb06eb90514e967df14c754b64e59018082 Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 16 Jun 2021 09:55:57 +0100
Subject: [PATCH] LISA angles tunable, block-diagonal correlation matrices

---
 examples/Nell_example.ipynb |  2 +-
 schnell/__init__.py         |  5 +++--
 schnell/correlation.py      | 35 ++++++++++++++++++++++++++++++++++-
 schnell/detector.py         |  8 +++++---
 4 files changed, 43 insertions(+), 7 deletions(-)

diff --git a/examples/Nell_example.ipynb b/examples/Nell_example.ipynb
index dac59be..f507ac1 100644
--- a/examples/Nell_example.ipynb
+++ b/examples/Nell_example.ipynb
@@ -537,7 +537,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.8.5"
   }
  },
  "nbformat": 4,
diff --git a/schnell/__init__.py b/schnell/__init__.py
index 4d450ee..8928c9f 100644
--- a/schnell/__init__.py
+++ b/schnell/__init__.py
@@ -5,8 +5,9 @@
     NoiseCorrelationConstantIdentity,
     NoiseCorrelationConstantR,
     NoiseCorrelationFromFunctions,
-    NoiseCorrelationLISA)
+    NoiseCorrelationLISA,
+    NoiseCorrelationBoxDiagonal)
 from .detector import (  # noqa
     Detector, GroundDetectorTriangle,
     GroundDetector, LISADetector)
-__version__ = '0.2.0'
+__version__ = '0.3.0'
diff --git a/schnell/correlation.py b/schnell/correlation.py
index c872125..8288576 100644
--- a/schnell/correlation.py
+++ b/schnell/correlation.py
@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.linalg import block_diag
 from .detector import LISADetector
 
 
@@ -9,7 +10,7 @@ class NoiseCorrelationBase(object):
     Do not use the bare class.
     """
     def __init__(self, ndet):
-        self.ndet
+        self.ndet = ndet
 
     def _get_corrmat(self, f):
         raise NotImplementedError("Don't use the NoiseCorrelationBase class")
@@ -81,6 +82,38 @@ def __init__(self, ndet, r):
                     np.full([self.ndet, self.ndet], r))
 
 
+class NoiseCorrelationBoxDiagonal(NoiseCorrelationBase):
+    """ This implements a correlation matrix from a set of
+    noise correlation objects. The resulting matrix will be
+    the combination of all correlation matrices in a
+    block-diagonal way.
+
+    Args:
+        corrs: list of `NoiseCorrelationBase` objects.
+    """
+    def __init__(self, corrs):
+        if not isinstance(corrs, list):
+            raise TypeError("`corr` must be a list of "
+                            "`NoiseCorrelation` objects.")
+        self.ncorrs = len(corrs)
+        for c in corrs:
+            if not isinstance(c, NoiseCorrelationBase):
+                raise TypeError("`corr` must be a list of "
+                                "`NoiseCorrelation` objects.")
+        self.corrs = corrs
+        self.ndet = np.sum([c.ndet for c in self.corrs])
+
+    def _get_corrmat(self, f):
+        f_use = np.atleast_1d(f)
+        mats = np.array([c.get_corrmat(f_use)
+                         for c in self.corrs])
+        mat = np.zeros([len(f_use), self.ndet, self.ndet])
+        for i, freq in enumerate(f_use):
+            ms = mats[:, i, :, :]
+            mat[i, :, :] = block_diag(*ms)
+        return mat
+
+
 class NoiseCorrelationFromFunctions(NoiseCorrelationBase):
     """ This implements a correlation matrix that has
     the same auto-correlation PSD for all detectors and
diff --git a/schnell/detector.py b/schnell/detector.py
index 5b90a6b..edbc040 100644
--- a/schnell/detector.py
+++ b/schnell/detector.py
@@ -365,15 +365,15 @@ class LISADetector(Detector):
         static (bool): if `True`, a static configuration corresponding
             to a perfect equilateral triangle in the x-y plane will
             be assumed.
+        kappa (float): initial longitude (degrees).
+        lambd (float): initial orientation angle (degrees).
     """
     trans_freq_earth = 2 * np.pi / (365 * 24 * 3600)
     R_AU = 1.496E11
-    kap = 0  # initial longitude
-    lam = 0  # initial orientation
     clight = 299792458.
 
     def __init__(self, detector_id, is_L5Gm=False,
-                 static=False):
+                 static=False, kappa=0, lambd=0):
         self.i_d = detector_id % 3
         self.name = 'LISA_%d' % self.i_d
         self.get_transfer = self._get_transfer_LISA
@@ -384,6 +384,8 @@ def __init__(self, detector_id, is_L5Gm=False,
             self.L = 2.5E9
             self.e = 0.00482419
         self.static = static
+        self.kap = np.radians(kappa)  # initial longitude
+        self.lam = np.radians(lambd)  # Orientation angle
 
     def _get_transfer_LISA(self, u, f, nv):
         # Eq. 48 in astro-ph/0105374