scatsystem: Comparing T-matrices

Former-commit-id: 3f1f012a8c518926ac8515884de898764807d516

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;
+}
+    
+
+
 

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) 

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);

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;