diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 2a63a3f..d90239c 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -3,6 +3,9 @@ find_package(GSL 2.0 REQUIRED) find_package(BLAS REQUIRED) find_package(LAPACK REQUIRED) +# disable an annoying warning that gives false positives probably due to a bug in gcc +# and other not very relevant warnings +set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wno-int-in-bool-context -Wno-comment") #includes set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS}) diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 602ac4a..5c1f23f 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -122,9 +122,10 @@ typedef struct qpms_ss_particle_orbitinfo { } qpms_ss_particle_orbitinfo_t; struct qpms_trans_calculator; +struct qpms_epsmu_generator_t; typedef struct qpms_scatsys_t { - qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. + struct qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. qpms_abstract_tmatrix_t **tm; ///< T-matrices in the system qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[]. diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 1e2c529..eed05e9 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1007,7 +1007,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { if(f->op.irot3arr.owns_ops) free(f->op.irot3arr.ops); break; - case QPMS_TMATRIX_OPERATION_FINITE_GROUP: // Group never owned + case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: // Group never owned break; case QPMS_TMATRIX_OPERATION_COMPOSE_SUM: { @@ -1015,7 +1015,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { &(f->op.compose_sum); if(o->opmem) { for(size_t i = 0; i < o->n; ++i) - qpms_tmatrix_operation_clear(o->ops[i]); + qpms_tmatrix_operation_clear(&(o->opmem[i])); free(o->opmem); } free(o->ops); @@ -1027,7 +1027,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { &(f->op.compose_chain); if(o->opmem) { for(size_t i = 0; i < o->n; ++i) - qpms_tmatrix_operation_clear(o->ops[i]); + qpms_tmatrix_operation_clear(&(o->opmem[i])); free(o->opmem); } free(o->ops); @@ -1037,3 +1037,62 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { QPMS_WTF; } } + +qpms_tmatrix_t *qpms_tmatrix_apply_operation( + const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) { + // Certain operations could be optimized, but the effect would be marginal. + qpms_tmatrix_t *res = qpms_tmatrix_copy(orig); + return qpms_tmatrix_apply_operation_inplace(f, res); +} + +static qpms_tmatrix_t *qtao_compose_sum_inplace(qpms_tmatrix_t *T, + const struct qpms_tmatrix_operation_compose_sum *cs) { + qpms_tmatrix_t *tmp_target = qpms_tmatrix_init(T->spec); + qpms_tmatrix_t *sum = qpms_tmatrix_init(T->spec); + for (size_t i = 0; i < cs->n; ++i) { + memcpy(tmp_target->m, T->m, SQ(T->spec->n) * sizeof(complex double)); + QPMS_ENSURE(qpms_tmatrix_apply_operation_inplace(cs->ops[i] , tmp_target), + "Got NULL pointer from qpms_tmatrix_apply_operation_inplace, hupsis!"); + for (size_t j = 0; j < SQ(T->spec->n); ++j) + sum->m[j] += tmp_target->m[j]; + } + for(size_t j = 0; j < SQ(T->spec->n); ++j) + T->m[j] = sum->m[j] * cs->factor; + qpms_tmatrix_free(sum); + qpms_tmatrix_free(tmp_target); + return T; +} + +static qpms_tmatrix_t *qtao_compose_chain_inplace(qpms_tmatrix_t *T, + const struct qpms_tmatrix_operation_compose_chain *cc) { + for(size_t i = 0; i < cc->n; ++i) + qpms_tmatrix_apply_operation_inplace(cc->ops[i], T); + return T; +} + +static qpms_tmatrix_t *qtao_scmulz_inplace(qpms_tmatrix_t *T, + const struct qpms_tmatrix_operation_scmulz *s) { + for(size_t i = 0; i < SQ(T->spec->n); ++i) + T->m[i] *= s->m[i]; + return T; +} + +qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace( + const qpms_tmatrix_operation_t *f, qpms_tmatrix_t *T) { + switch(f->typ) { + case QPMS_TMATRIX_OPERATION_LRMATRIX: + return qpms_tmatrix_apply_symop_inplace(T, f->op.lrmatrix.m); + case QPMS_TMATRIX_OPERATION_IROT3: + return qpms_tmatrix_symmetrise_irot3arr_inplace(T, 1, &(f->op.irot3)); + case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: + return qpms_tmatrix_symmetrise_finite_group_inplace(T, f->op.finite_group); + case QPMS_TMATRIX_OPERATION_COMPOSE_SUM: + return qtao_compose_sum_inplace(T, &(f->op.compose_sum)); + case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN: + return qtao_compose_chain_inplace(T, &(f->op.compose_chain)); + case QPMS_TMATRIX_OPERATION_SCMULZ: + return qtao_scmulz_inplace(T, &(f->op.scmulz)); + default: + QPMS_WTF; + } +} diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index 83aa71b..163c6d3 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -235,7 +235,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( complex double *qpms_apply_tmatrix( complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated. const complex double *a, ///< Incident field coefficient array of size T->spec->n. - const qpms_tmatrix_t *T + const qpms_tmatrix_t *T ///< T-matrix \a T to apply. ); /// Generic T-matrix generator function that fills a pre-initialised qpms_tmatrix_t given a frequency. @@ -429,7 +429,7 @@ qpms_arc_function_retval_t qpms_arc_cylinder(double theta, ); /// Arc parametrisation of spherical particle; for qpms_arc_function_t. -/** Useful mostly only for benchmarks, as one can use the Mie-Lorentz solution. */ +/** Useful mostly only for benchmarks or debugging, 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. ); @@ -520,7 +520,7 @@ typedef struct qpms_tmatrix_function_t { /// Specifies different kinds of operations done on T-matrices typedef enum { - QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f[ T' = MTM^\dagger \f]; see @ref qpms_tmatrix_operation_lrmatrix. + QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f$ T' = MTM^\dagger \f$; see @ref qpms_tmatrix_operation_lrmatrix. QPMS_TMATRIX_OPERATION_IROT3, ///< Single rotoreflection specified by a qpms_irot3_t. QPMS_TMATRIX_OPERATION_IROT3ARR, ///< Symmetrise using an array of rotoreflection; see @ref qpms_tmatrix_operation_irot3arr. QPMS_TMATRIX_OPERATION_COMPOSE_SUM, ///< Apply several transformations and sum the results, see @ref qpms_tmatrix_operation_compose_sum. @@ -534,13 +534,13 @@ typedef enum { struct qpms_tmatrix_operation_lrmatrix { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ - const complex double *m; + complex double *m; bool owns_m; ///< Whether \a m is owned by this; }; struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations. -/// Specifies a composed operation of type \f$ T' = c\Sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t. +/// Specifies a composed operation of type \f$ T' = c\sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t. struct qpms_tmatrix_operation_compose_sum { size_t n; ///< Number of operations in ops; const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers array \a ops[] always owned by this.) @@ -550,21 +550,21 @@ struct qpms_tmatrix_operation_compose_sum { * the memory held by \a opmem and to be properly initialised. * Each \a ops member has to point to _different_ elements of \a opmem. */ - qpms_tmatrix_operation_t *opmem; + struct qpms_tmatrix_operation_t *opmem; }; /// Specifies a composed operation of type \f$ T' = f_{n-1}(f_{n-2}(\dots f_0(T)\dots))) \f$ for qpms_tmatrix_operation_t. struct qpms_tmatrix_operation_compose_chain { size_t n; ///< Number of operations in ops; const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers owned by this.) - qpms_tmatrix_operation_t *opmem; ///< (Optional) operations buffer into which elements of \a ops point. (Owned by this or NULL.) + struct qpms_tmatrix_operation_t *opmem; ///< (Optional) operations buffer into which elements of \a ops point. (Owned by this or NULL.) }; /// Specifies an elementwise complex multiplication of type \f$ T'_{ij} = M_{ij}T_{ij} \f$ for qpms_tmatrix_operation_t. struct qpms_tmatrix_operation_scmulz { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ - const complex double *m; + complex double *m; bool owns_m; ///< Whether \a m is owned by this. }; @@ -573,7 +573,7 @@ struct qpms_tmatrix_operation_scmulz { * or qpms_symmetrise_tmdata_irot3arr_inplace(). */ struct qpms_tmatrix_operation_irot3arr { size_t n; ///< Number of rotoreflections; - const qpms_irot3_t *ops; ///< Rotoreflection array of size \a n. + qpms_irot3_t *ops; ///< Rotoreflection array of size \a n. bool owns_ops; ///< Whether \a ops array is owned by this. }; @@ -584,14 +584,13 @@ typedef struct qpms_tmatrix_operation_t { union { struct qpms_tmatrix_operation_lrmatrix lrmatrix; struct qpms_tmatrix_operation_scmulz scmulz; - qpms_irot3 irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3 - struct qpms_tmatrix_operation_finite_group finite_group; + qpms_irot3_t irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3 struct qpms_tmatrix_operation_irot3arr irot3arr; struct qpms_tmatrix_operation_compose_sum compose_sum; struct qpms_tmatrix_operation_compose_chain compose_chain; /// Finite group for QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE. /** Not owned by this; \a rep3d must be filled. */ - const qpms_finite_group_t *finitegroup; + const qpms_finite_group_t *finite_group; } op; ///< Operation data; actual type is determined by \a typ. } qpms_tmatrix_operation_t;