tmatrices.h

Go to the documentation of this file.

 
#ifndef TMATRICES_H
#define TMATRICES_H
// #include "qpms_types.h" // included via materials.h
// #include <gsl/gsl_spline.h> // included via materials.h
#include "materials.h"
#include <stdio.h>
 
 
 
struct qpms_finite_group_t;
typedef struct qpms_finite_group_t qpms_finite_group_t;
 
static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
    return t->m + (t->spec->n * rowno);
}
 
 
qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec);
 
 
qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T);
 
 
qpms_tmatrix_t *qpms_tmatrix_mv(qpms_tmatrix_t *dest, const qpms_tmatrix_t *orig);
 
void qpms_tmatrix_free(qpms_tmatrix_t *t);
 
 
bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
        const double rtol, const double atol);
 
 
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
        const qpms_tmatrix_t *T, 
        const complex double *M 
        );
 
 
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
        qpms_tmatrix_t *T, 
        const complex double *M 
        );
 
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
        const qpms_tmatrix_t *T, 
        const complex double *M 
        );
 
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(
        const qpms_tmatrix_t *T 
        );
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
        const qpms_tmatrix_t *T, 
        int N 
        );
 
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
        qpms_tmatrix_t *T, 
        const complex double *M 
        );
 
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(
        qpms_tmatrix_t *T 
        );
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(
        qpms_tmatrix_t *T, 
        int N 
        );
 
 
qpms_errno_t qpms_load_scuff_tmatrix(
        const char *path, 
        const qpms_vswf_set_spec_t *bspec, 
        size_t *n, 
        double **freqs, 
        double **freqs_su, 
        qpms_tmatrix_t **tmatrices_array, 
        complex double **tmdata 
        );
 
 
extern bool qpms_load_scuff_tmatrix_crash_on_failure;
 
 
qpms_errno_t qpms_read_scuff_tmatrix(
        FILE *f, 
        const qpms_vswf_set_spec_t *bspec, 
        size_t *n, 
        double **freqs, 
        double **freqs_su, 
        qpms_tmatrix_t **tmatrices_array, 
 
        complex double ** tmdata 
        );
 
 
qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
        complex double *tmdata, const size_t tmcount,
        const qpms_vswf_set_spec_t *bspec,
        size_t n_symops,
        const qpms_irot3_t *symops
        );
 
 
qpms_errno_t qpms_symmetrise_tmdata_finite_group(
        complex double *tmdata, const size_t tmcount,
        const qpms_vswf_set_spec_t *bspec,
        const qpms_finite_group_t *pointgroup
        );
 
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
        qpms_tmatrix_t *T,
        size_t n_symops,
        const qpms_irot3_t *symops
        );
 
 
qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
        qpms_tmatrix_t *T,
        const qpms_finite_group_t *pointgroup
        );
 
complex double *qpms_apply_tmatrix(
        complex double *f_target, 
        const complex double *a, 
        const qpms_tmatrix_t *T 
        );
 
 
typedef struct qpms_tmatrix_generator_t {
    qpms_errno_t (*function) (qpms_tmatrix_t *t, 
            complex double omega, 
            const void *params 
            );
    const void *params; 
} qpms_tmatrix_generator_t;
 
qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
        const qpms_vswf_set_spec_t *bspec,
        qpms_tmatrix_generator_t gen,
        complex double omega);
 
 
 
qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
        complex double omega,
        const void *tmatrix_orig
        );
 
/* Fuck this, include the whole <gsl/gsl_spline.h>
typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct
typedef struct gsl_interp_type gsl_interp_type;
extern const gsl_interp_type * gsl_interp_linear;
extern const gsl_interp_type * gsl_interp_polynomial;
extern const gsl_interp_type * gsl_interp_cspline;
extern const gsl_interp_type * gsl_interp_cspline_periodic;
extern const gsl_interp_type * gsl_interp_akima;
extern const gsl_interp_type * gsl_interp_akima_periodic;
extern const gsl_interp_type * gsl_interp_steffen;
*/
 
// struct gsl_interp_accel; // use if lookup proves to be too slow
typedef struct qpms_tmatrix_interpolator_t {
    const qpms_vswf_set_spec_t *bspec;
    //bool owns_bspec;
    gsl_spline **splines_real; 
    gsl_spline **splines_imag; 
    // gsl_interp_accel **accel_real;
    // gsl_interp_accel **accel_imag;
} qpms_tmatrix_interpolator_t;
 
void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp);
 
qpms_errno_t qpms_tmatrix_interpolator_eval_fill(qpms_tmatrix_t *target, 
        const qpms_tmatrix_interpolator_t *interp, double freq);
 
 
qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq);
 
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, 
        const double *freqs, const qpms_tmatrix_t *tmatrices_array, 
        const gsl_interp_type *iptype
           //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
               /*, ...? */);
 
 
 
qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, 
            complex double omega, 
            const void *interpolator 
);
 
 
complex double *qpms_mie_coefficients_reflection(
                complex double *target, 
                const qpms_vswf_set_spec_t *bspec, 
                double a, 
                complex double k_i, 
                complex double k_e, 
                complex double mu_i, 
                complex double mu_e, 
                qpms_bessel_t J_ext, 
                qpms_bessel_t J_scat 
                );
 
qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, 
        double a, 
        complex double k_i, 
        complex double k_e, 
        complex double mu_i, 
        complex double mu_e 
        );
 
typedef struct qpms_tmatrix_generator_sphere_param_t {
    qpms_epsmu_generator_t outside;
    qpms_epsmu_generator_t inside;
    double radius;
} qpms_tmatrix_generator_sphere_param_t;
 
qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
        complex double omega,
        const void *params 
        );
 
static inline qpms_tmatrix_t *qpms_tmatrix_spherical(
        const qpms_vswf_set_spec_t *bspec,
        double a, 
        complex double k_i, 
        complex double k_e, 
        complex double mu_i, 
        complex double mu_e 
        ) {
    qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
    qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e);
    return t;
}
 
qpms_errno_t qpms_tmatrix_spherical_mu0_fill(
        qpms_tmatrix_t *t, 
        double a, 
        double omega, 
        complex double epsilon_fg, 
        complex double epsilon_bg 
        ); 
 
static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(
        const qpms_vswf_set_spec_t *bspec,
        double a, 
        double omega, 
        complex double epsilon_fg, 
        complex double epsilon_bg 
        ) {
    qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
    qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg);
    return t;
};
 
typedef struct qpms_arc_function_retval_t {
    double r; 
    double beta; 
} qpms_arc_function_retval_t;
 
typedef struct qpms_arc_function_t {
 
    qpms_arc_function_retval_t (*function) (
            double theta, 
            const void *params 
            );
    const void *params;
} qpms_arc_function_t;
 
typedef struct qpms_arc_cylinder_params_t {
    double R; 
    double h; 
} qpms_arc_cylinder_params_t;
 
qpms_arc_function_retval_t qpms_arc_cylinder(double theta,
        const void *params 
        );
 
 
qpms_arc_function_retval_t qpms_arc_sphere(double theta,
        const void *R 
        );
 
 
 
qpms_errno_t qpms_tmatrix_axialsym_fill(
        qpms_tmatrix_t *t, 
        complex double omega, 
        qpms_epsmu_t outside, 
        qpms_epsmu_t inside, 
        qpms_arc_function_t shape, 
        qpms_l_t lMax_extend 
        );
 
static inline qpms_tmatrix_t *qpms_tmatrix_axialsym(
        const qpms_vswf_set_spec_t *bspec,
        complex double omega, 
        qpms_epsmu_t outside, 
        qpms_epsmu_t inside, 
        qpms_arc_function_t shape, 
 
        qpms_l_t lMax_extend 
        ) {
    qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
    qpms_tmatrix_axialsym_fill(t, omega, outside, inside, shape, lMax_extend);
    return t;
}
 
typedef struct qpms_tmatrix_generator_axialsym_param_t {
    qpms_epsmu_generator_t outside;
    qpms_epsmu_generator_t inside;
    qpms_arc_function_t shape;
    qpms_l_t lMax_extend;
} qpms_tmatrix_generator_axialsym_param_t;
 
 
qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, 
            complex double omega, 
            const void *params 
);
 
 
qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(complex double *target,
        complex double omega,
        const qpms_tmatrix_generator_axialsym_param_t *param,
        qpms_normalisation_t norm,
        qpms_bessel_t J
        );
 
qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target,
        complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,
        qpms_arc_function_t shape, qpms_l_t lMaxQR, qpms_normalisation_t norm,
        qpms_bessel_t J 
        );
 
 
typedef struct qpms_tmatrix_function_t {
        const qpms_vswf_set_spec_t *spec;
        const qpms_tmatrix_generator_t *gen; 
} qpms_tmatrix_function_t;
 
// FIXME the name is not very intuitive.
static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, complex double omega) {
    return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega);
}
 
typedef enum {
    QPMS_TMATRIX_OPERATION_NOOP, 
    QPMS_TMATRIX_OPERATION_LRMATRIX, 
    QPMS_TMATRIX_OPERATION_IROT3, 
    QPMS_TMATRIX_OPERATION_IROT3ARR, 
    QPMS_TMATRIX_OPERATION_COMPOSE_SUM, 
    QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN, 
    QPMS_TMATRIX_OPERATION_SCMULZ, 
    //QPMS_TMATRIX_OPERATION_POINTGROUP, ///< TODO the new point group in pointgroup.h
    QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE 
} qpms_tmatrix_operation_kind_t;
 
struct qpms_tmatrix_operation_lrmatrix {
 
    complex double *m; 
        size_t m_size; 
    bool owns_m; 
};
 
struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations.
 
struct qpms_tmatrix_operation_compose_sum {
    size_t n; 
    const struct qpms_tmatrix_operation_t **ops; 
    double factor; 
 
    struct qpms_tmatrix_operation_t *opmem; 
};
 
struct qpms_tmatrix_operation_compose_chain {
    size_t n; 
    const struct qpms_tmatrix_operation_t **ops; 
    struct qpms_tmatrix_operation_t *opmem; 
    size_t opmem_size; 
    _Bool *ops_owned; 
};
 
struct qpms_tmatrix_operation_scmulz {
 
    complex double *m;
        size_t m_size; 
    bool owns_m; 
};
 
 
struct qpms_tmatrix_operation_irot3arr {
    size_t n; 
    qpms_irot3_t *ops; 
    bool owns_ops; 
};
 
typedef struct qpms_tmatrix_operation_t {
    qpms_tmatrix_operation_kind_t typ;
    union {
        struct qpms_tmatrix_operation_lrmatrix lrmatrix;
        struct qpms_tmatrix_operation_scmulz scmulz;
        qpms_irot3_t irot3; 
        struct qpms_tmatrix_operation_irot3arr irot3arr;
        struct qpms_tmatrix_operation_compose_sum compose_sum;
        struct qpms_tmatrix_operation_compose_chain compose_chain;
 
        const qpms_finite_group_t *finite_group; 
    } op; 
} qpms_tmatrix_operation_t;
 
static const qpms_tmatrix_operation_t qpms_tmatrix_operation_noop = {.typ = QPMS_TMATRIX_OPERATION_NOOP};
 
qpms_tmatrix_t *qpms_tmatrix_apply_operation(const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig);
 
qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(const qpms_tmatrix_operation_t *op, qpms_tmatrix_t *orig);
 
qpms_tmatrix_t *qpms_tmatrix_apply_operation_replace(qpms_tmatrix_t *dest, 
        const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig);
 
 
void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f);
 
 
void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *target, const qpms_tmatrix_operation_t *src);
 
void qpms_tmatrix_operation_compose_chain_init(
        qpms_tmatrix_operation_t *target, 
        size_t nops, 
        size_t opmem_size 
        );
 
#if 0
// Abstract types that describe T-matrix/particle/scatsystem symmetries
// To be implemented later. See also the thoughts in the beginning of groups.h.
 
typedef *char qpms_tmatrix_id_t; 
 
 
typedef struct qpms_abstract_tmatrix_t{
    qpms_tmatrix_id_t id;
    qpms_irot3_t *invar_gens;
    qpms_gmi_t invar_gens_size;
 
} qpms_abstract_tmatrix_t;
 
 
typedef struct qpms_abstract_particle_t{
} qpms_abstract_particle_t;
 
typedef struct qpms_abstract_particle_t {
    cart3_t pos; 
    const qpms_abstract_tmatrix_t *tmatrix; 
} qpms_abstract_particle_t;
 
 
typedef qpms_particle_tid_t qpms_abstract_particle_tid_t;
#endif // 0
 
#endif //TMATRICES_H