From beef0ea9b8c7a4e277cf88ab6ab099c3c1d2a1cc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Thu, 8 Aug 2019 22:18:42 +0300
Subject: [PATCH] WIP T-matrix of axially symmetric T-matrix.

Former-commit-id: b699e21f6693607b1a790cbdf7a0d7986c6ddff9
---
 qpms/tmatrices.c | 36 +++++++++++++++++++++++++++
 qpms/tmatrices.h | 64 ++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 100 insertions(+)

diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c
index 48fdb4d..e99ed98 100644
--- a/qpms/tmatrices.c
+++ b/qpms/tmatrices.c
@@ -582,4 +582,40 @@ complex double *qpms_apply_tmatrix(
   return f;
 }
 
+qpms_arc_function_retval_t qpms_arc_sphere(double theta, const void *R) {
+  qpms_arc_function_retval_t retval = {*(const double*)R, 0};
+  return retval;
+}
 
+qpms_arc_function_retval_t qpms_arc_cylinder(double theta, const void *param) {
+  const qpms_arc_cylinder_params_t *p = param;
+  double thresh = atan(2 * p->R / p->h);
+  QPMS_ENSURE(theta >= 0 && theta <= M_PI, 
+      "theta = %g, but it must lie in interval [0, M_PI]", theta);
+  qpms_arc_function_retval_t res;
+  if (theta < thresh) {
+    // Upper base
+    res.r = 0.5 * p->h / cos(theta);
+    res.beta = -theta;
+  } else if (theta <= M_PI - thresh) {
+    // Side
+    res.r = p->R / cos(theta - M_PI_2);
+    res.beta = -theta + M_PI_2;
+  } else {
+    // Lower base
+    res.r = 0.5 * p->h / cos(theta - M_PI);
+    res.beta = -theta + M_PI;
+  }
+  return res;
+}
+
+#if 0
+qpms_errno_t qpms_tmatrix_axialsym_fill(
+		qpms_tmatrix_t *t, complex double omega, qpms_epsmu_generator_t outside,
+		qpms_epsmu_generator_t inside,qpms_arc_function_t shape)
+{
+
+
+}
+
+#endif
diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h
index e0a9bf9..3bba2e0 100644
--- a/qpms/tmatrices.h
+++ b/qpms/tmatrices.h
@@ -345,8 +345,72 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(
 		) {
 	qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
 	qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg);
+	return t;
 };
 
+/// Return value type for qpms_arc_function_t.
+typedef struct qpms_arc_function_retval_t {
+	double r; ///< Distance from the origin.
+	double beta; ///< Angle between surface normal and radial direction.
+} qpms_arc_function_retval_t;
+
+/// Prototype for general parametrisation of \f$ C_\infty \f$-symmetric particle's surface.
+typedef struct qpms_arc_function_t {
+	/// Arc parametrisation function.
+	/** TODO link to notes.
+	 *
+	 * Implemented by:
+	 * qpms_arc_cylinder(),
+	 * qpms_arc_sphere().
+	 */
+	qpms_arc_function_retval_t (*function) (
+			double theta, ///< Polar angle from interval \f$ [0, \pi] \f$. 
+			const void *params ///< Pointer to implementation specific parameters.
+			);
+	const void *params;
+} qpms_arc_function_t;
+
+/// Parameter structure for qpms_arc_cylinder().
+typedef struct qpms_arc_cylinder_params_t {
+	double R; ///< Cylinder radius.
+	double h; ///< Cylinder height.
+} qpms_arc_cylinder_params_t;
+
+/// Arc parametrisation of cylindrical particle; for qpms_arc_function_t.
+qpms_arc_function_retval_t qpms_arc_cylinder(double theta,
+		const void *params; ///< Points to qpms_arc_cylinder_params_t
+		);
+
+/// Arc parametrisation of spherical particle; for qpms_arc_function_t.
+/** Useful mostly only for benchmarks, as one can use the Mie-Lorentz solution. */
+qpms_arc_function_retval_t qpms_arc_sphere(double theta,
+		const void *R; ///< Points to double containing particle's radius.
+		);
+
+
+/// Replaces T-matrix contents with those of a particle with \f$ C_\infty \f$ symmetry.
+qpms_errno_t qpms_tmatrix_axialsym_fill(
+		qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
+		complex double omega, ///< Angular frequency.
+		qpms_epsmu_generator_t outside, ///< Optical properties of the outside medium.
+		qpms_epsmu_generator_t inside, ///< Optical properties of the particle's material.
+		qpms_arc_function_t shape ///< Particle surface parametrisation.
+		);
+
+/// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry.
+static inline qpms_tmatrix_t *qpms_tmatrix_axialsym(
+		const qpms_vswf_set_spec_t *bspec,
+		complex double omega, ///< Angular frequency.
+		qpms_epsmu_generator_t outside, ///< Optical properties of the outside medium.
+		qpms_epsmu_generator_t inside, ///< Optical properties of the particle's material.
+		qpms_arc_function_t shape ///< Particle surface parametrisation.
+		) {
+	qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
+	qpms_tmatrix_axialsym_fill(t, omega, outside, inside, shape);
+	return t;
+}
+
+
 
 #if 0
 // Abstract types that describe T-matrix/particle/scatsystem symmetries