Per-orbit irrep projectors.

Former-commit-id: 7d66eb721b7f75a99c29a054f0fd60520eb28722

qpms/scatsystem.c

@@ -562,3 +562,36 @@ complex double *qpms_orbit_action_matrix(complex double *target,
   return target;
 }
 
+complex double *qpms_orbit_irrep_projector_matrix(complex double *target,
+    const qpms_ss_orbit_type_t *ot, const qpms_vswf_set_spec_t *bspec,
+    const qpms_finite_group_t *sym, const qpms_iri_t iri) {
+  assert(sym);
+  assert(sym->rep3d);
+  assert(ot); assert(ot->size > 0);
+  assert(iri < sym->nirreps); assert(sym->irreps);
+  check_norm_compat(bspec);
+  const size_t n = bspec->n;
+  const qpms_gmi_t N = ot->size;
+  if (target == NULL)
+    target = malloc(n*n*N*N*sizeof(complex double));
+  if (target == NULL) abort();
+  memset(target, 0, n*n*N*N*sizeof(complex double));
+  // Workspace for the orbit group action matrices
+  complex double *tmp = malloc(n*n*N*N*sizeof(complex double));
+  const int d = sym->irreps[iri].dim;
+  double prefac = d / (double) sym->order;
+  for(int partner = 0; partner < d; ++partner) {
+    for(qpms_gmi_t g = 0; g < sym->order; ++g) {
+      // We use the diagonal elements of D_g
+      complex double D_g_conj = sym->irreps[iri].m[g*d*d + partner*d + partner];
+      qpms_orbit_action_matrix(tmp, ot, bspec, sym, g);
+      // TODO kahan sum?
+      for(size_t i = 0; i < n*n*N*N; ++i)
+        target[i] += prefac * D_g_conj * tmp[i];
+    }
+  }
+  free(tmp);
+  return target;
+}
+
+

qpms/scatsystem.h

@@ -327,7 +327,23 @@ complex double *qpms_orbit_action_matrix(
 		/// The index of the operation in sym to represent.
 		const qpms_gmi_t g);
 
-
+/// Constructs a dense matrix representation of a irrep projector for an orbit type.
+/** TODO detailed doc */
+complex double *qpms_orbit_irrep_projector_matrix(
+		/// Target array. If NULL, a new one is allocated.
+		/** The size of the array is (orbit->size * bspec->n)**2
+		 * (it makes sense to assume all the T-matrices share their spec).
+		 */
+		complex double *target,
+		/// The orbit (type).
+		const qpms_ss_orbit_type_t *orbit,
+		/// Base spec of the t-matrices (we don't know it from orbit, as it has 
+		/// only T-matrix indices.
+		const qpms_vswf_set_spec_t *bspec,
+		/// The symmetry group used to generate the orbit (must have rep3d filled).
+		const struct qpms_finite_group_t *sym,
+		/// The index of the irreducible representation of sym.
+		const qpms_iri_t iri);
 
 #if 0
 // Abstract types that describe T-matrix/particle/scatsystem symmetries