From f01e9ba34b37dde057f57c532b10e5629b5f6c75 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Thu, 21 Feb 2019 09:40:30 +0000
Subject: [PATCH] scatsystem: Comparing T-matrices

Former-commit-id: 3f1f012a8c518926ac8515884de898764807d516
---
 qpms/scatsystem.c | 18 ++++++++++++++++++
 qpms/scatsystem.h |  8 +++++++-
 qpms/vswf.c       | 10 ++++++++++
 qpms/vswf.h       |  7 +++++++
 4 files changed, 42 insertions(+), 1 deletion(-)

diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index 4021eed..5001b05 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -128,4 +128,22 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) {
   return T;
 }
 
+bool qpms_tmatrix_isnear(const qpms_tmatrix *A, const qpms_tmatrix *B,
+    const double rtol, const double atol)
+{
+  if (!qpms_vswf_set_specisidentical(A->spec, B->spec))
+    return false;
+  if (A->m == B->m)
+    return true;
+  const size_t n = A->spec->n;
+  for (size_t i = 0; i < n*n; ++i) {
+    const double tol = atol + rtol * (cabs(B->m[i]));
+    if ( cabs(B->m[i] - A->m[i]) > tol )
+      return false;
+  }
+  return true;
+}
+    
+
+
 
diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h
index 1c6ed17..4ff4588 100644
--- a/qpms/scatsystem.h
+++ b/qpms/scatsystem.h
@@ -46,7 +46,13 @@ typedef struct qpms_particle_t {
 	const qpms_tmatrix_t *tmatrix; ///< T-matrix; not owned by qpms_particle_t.
 } qpms_particle_t;
 
-
+/// Check T-matrix equality/similarity.
+/** 
+ *  This function actually checks for identical vswf specs.
+ * TODO define constants with "default" atol, rtol for this function.
+ */
+bool qpms_tmatrix_isnear(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
+		const double rtol, const double atol);
 
 /// Creates a T-matrix from another matrix and a symmetry operation.
 /** The symmetry operation is expected to be a unitary (square) 
diff --git a/qpms/vswf.c b/qpms/vswf.c
index bdd8589..da237b0 100644
--- a/qpms/vswf.c
+++ b/qpms/vswf.c
@@ -49,6 +49,16 @@ qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *s, const qpms_uvswf
   return QPMS_SUCCESS;
 }
 
+bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
+                    const qpms_vswf_set_spec_t *b) {
+  if (a == b) return true;
+  if (a->n != b->n) return false;
+  for (size_t i = 0; i < a->n; ++i)
+    if (a->ilist[i] != b->ilist[i])
+      return false;
+  return true;
+}
+
 void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *s) {
   if(s) free(s->ilist);
   free(s);
diff --git a/qpms/vswf.h b/qpms/vswf.h
index 1389aa8..3bf4d4c 100644
--- a/qpms/vswf.h
+++ b/qpms/vswf.h
@@ -34,6 +34,13 @@ qpms_vswf_set_spec_t *qpms_vswf_set_spec_init();
 qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *self, qpms_uvswfi_t u);
 /// Destroys a \ref qpms_vswf_set_spec_t.
 void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *);
+/// Compares two vswf basis specs.
+/**
+ * Checks whether ilist is the same and of the same length.
+ * If yes, returns true, else returns false.
+ */
+bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
+		const qpms_vswf_set_spec_t *b);
 
 /// NOT IMPLEMENTED Evaluates a set of VSWF basis functions at a given point.
 /** The list of basis wave indices is specified in \a setspec;