scatsystem.h

Go to the documentation of this file.

 
#ifndef QPMS_SCATSYSTEM_H
#define QPMS_SCATSYSTEM_H
#include "qpms_types.h"
#include "vswf.h"
#include "tmatrices.h"
#include <stdbool.h>
 
 
void qpms_scatsystem_set_nthreads(long n);
 
 
typedef struct qpms_particle_t {
    // Does it make sense to ever use other than cartesian coords for this?
    cart3_t pos; 
    const qpms_tmatrix_function_t *tmg; 
    qpms_tmatrix_operation_t op; 
} qpms_particle_t;
 
struct qpms_finite_group_t;
typedef struct qpms_finite_group_t qpms_finite_group_t;
 
typedef struct qpms_particle_tid_t {
    // Does it make sense to ever use other than cartesian coords for this?
    cart3_t pos; 
    qpms_ss_tmi_t tmatrix_id; 
} qpms_particle_tid_t;
 
 
typedef qpms_gmi_t qpms_ss_orbit_pi_t; 
typedef qpms_ss_tmi_t qpms_ss_oti_t; 
 
 
typedef struct qpms_ss_orbit_type_t {
    qpms_ss_orbit_pi_t size; 
    size_t bspecn; 
 
    qpms_ss_orbit_pi_t *action;
 
    qpms_ss_tmi_t *tmatrices;
 
    size_t *irbase_sizes;
    //The following are pretty redundant, TODO reduce them at some point.
    size_t *irbase_cumsizes;
    size_t *irbase_offsets;
 
    complex double *irbases;
    size_t instance_count;
    qpms_ss_pi_t p_offset;
} qpms_ss_orbit_type_t;
 
 
typedef ptrdiff_t qpms_ss_osn_t; 
 
typedef struct qpms_ss_particle_orbitinfo {
    qpms_ss_oti_t t; 
#define QPMS_SS_P_ORBITINFO_UNDEF (-1) 
    qpms_ss_osn_t osn; 
    qpms_ss_orbit_pi_t p; 
} qpms_ss_particle_orbitinfo_t;
 
typedef struct qpms_ss_derived_tmatrix_t {
    qpms_ss_tmgi_t tmgi; 
    struct qpms_tmatrix_operation_t op; 
} qpms_ss_derived_tmatrix_t;
 
typedef struct qpms_scatsys_periodic_info_t {   
 
    cart3_t lattice_basis[3];
 
 
    cart3_t reciprocal_basis[3];
 
 
    double unitcell_volume;
 
 
    double eta;
} qpms_scatsys_periodic_info_t;
 
 
struct qpms_trans_calculator;
struct qpms_epsmu_generator_t;
 
 
typedef struct qpms_scatsys_t {
    int lattice_dimension;
    struct qpms_epsmu_generator_t medium; 
    
 
    qpms_tmatrix_function_t *tmg; 
    qpms_ss_tmgi_t tmg_count; 
 
    qpms_ss_derived_tmatrix_t *tm; 
    qpms_ss_tmi_t tm_count; 
    qpms_ss_tmi_t tm_capacity; 
    
    qpms_particle_tid_t *p; 
    qpms_ss_pi_t p_count; 
    qpms_ss_pi_t p_capacity; 
 
    //TODO the index types do not need to be so big.
    const struct qpms_finite_group_t *sym; 
    qpms_ss_pi_t *p_sym_map; 
    qpms_ss_tmi_t *tm_sym_map; 
 
    qpms_ss_oti_t orbit_type_count;
    qpms_ss_orbit_type_t *orbit_types; 
 
    qpms_ss_particle_orbitinfo_t *p_orbitinfo; 
 
    size_t fecv_size; 
    size_t *saecv_sizes; 
 
    size_t *fecv_pstarts; 
    size_t *saecv_ot_offsets; 
    //size_t **saecv_pstarts; ///< NI. Indices of where pi'th particle's excitation coeff start in a symmetry-adjusted e.c.v.
 
    // TODO shifted origin of the symmetry group etc.
    // TODO some indices for fast operations here.
    // private
    size_t max_bspecn; 
    
    qpms_ss_pi_t *p_by_orbit;
 
    // We keep the p_orbitinfo arrays in this chunk in order to avoid memory fragmentation
    char *otspace;
    char *otspace_end;
    double lenscale; // radius of the array, used as a relative tolerance measure
    struct qpms_trans_calculator *c;
 
    qpms_scatsys_periodic_info_t per;
} qpms_scatsys_t;
 
static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) {
    return ss->tmg[ss->tm[tmi].tmgi].spec;
}
 
static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) {
    return ss->tmg[ss->tm[ss->p[pi].tmatrix_id].tmgi].spec;
}
 
typedef struct qpms_scatsys_at_omega_t {
    const qpms_scatsys_t *ss; 
 
    qpms_tmatrix_t **tm;
    complex double omega; 
    qpms_epsmu_t medium; 
    complex double wavenumber; 
} qpms_scatsys_at_omega_t;
 
 
 
qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym,
        complex double omega, const struct qpms_tolerance_spec_t *tol);
 
void qpms_scatsys_free(qpms_scatsys_t *s);
 
 
void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw);
 
 
qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss,
        complex double omega);
 
 
complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U,
        const qpms_scatsys_t *ss, qpms_iri_t iri);
 
 
complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed,
        const complex double *orig_full, const qpms_scatsys_t *ss,
        qpms_iri_t iri);
 
 
complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full,
        const complex double *orig_packed, const qpms_scatsys_t *ss,
        qpms_iri_t iri, bool add);
 
 
complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
        const complex double *orig_full, const qpms_scatsys_t *ss,
        qpms_iri_t iri);
 
 
complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
        const complex double *orig_packed, const qpms_scatsys_t *ss,
        qpms_iri_t iri, bool add);
 
 
complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
        const complex double *orig_full, const qpms_scatsys_t *ss,
        qpms_iri_t iri);
 
 
complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
        const complex double *orig_packed, const qpms_scatsys_t *ss,
        qpms_iri_t iri, bool add);
 
 
complex double *qpms_scatsys_build_translation_matrix_full(
        complex double *target,
        const qpms_scatsys_t *ss,
        complex double k 
        );
 
 
complex double *qpms_scatsysw_build_modeproblem_matrix_full(
        complex double *target,
        const qpms_scatsys_at_omega_t *ssw
        );
 
 
complex double *qpms_scatsys_build_translation_matrix_e_full(
        complex double *target,
        const qpms_scatsys_t *ss,
        complex double k, 
        qpms_bessel_t J
        );
 
 
complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
        complex double *target,
        const qpms_scatsys_t *ss,
        qpms_iri_t iri,
        complex double k, 
        qpms_bessel_t J
        );
 
 
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
        complex double *target,
        const qpms_scatsys_at_omega_t *ssw,
        qpms_iri_t iri 
        );
 
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
        complex double *target,
        const qpms_scatsys_at_omega_t *ssw,
        qpms_iri_t iri 
        );
 
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
        complex double *target,
        const qpms_scatsys_at_omega_t *ssw,
        qpms_iri_t iri 
        );
 
struct qpms_scatsys_at_omega_k_t; // Defined below.
typedef struct qpms_ss_LU {
    const qpms_scatsys_at_omega_t *ssw;
    const struct qpms_scatsys_at_omega_k_t *sswk; 
    bool full; 
    qpms_iri_t iri; 
    complex double *a;
    int *ipiv; 
} qpms_ss_LU;
void qpms_ss_LU_free(qpms_ss_LU);
 
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
        complex double *target, 
        int *target_piv, 
        const qpms_scatsys_at_omega_t *ssw
        );
 
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
        complex double *target, 
        int *target_piv, 
        const qpms_scatsys_at_omega_t *ssw,
        qpms_iri_t iri
        );
 
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
        complex double *modeproblem_matrix_full, 
        int *target_piv, 
        const qpms_scatsys_at_omega_t *ssw, 
        const struct qpms_scatsys_at_omega_k_t *sswk 
        );
 
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(
        complex double *modeproblem_matrix_irrep_packed, 
        int *target_piv, 
        const qpms_scatsys_at_omega_t *ssw,
        qpms_iri_t iri
        );
 
complex double *qpms_scatsys_scatter_solve(
        complex double *target_f, 
        const complex double *a_inc, 
        qpms_ss_LU ludata 
        );
 
// ======================= Periodic system -only related stuff =============================
 
 
typedef struct qpms_scatsys_at_omega_k_t {
    const qpms_scatsys_at_omega_t *ssw;
    double k[3]; 
    double eta; 
} qpms_scatsys_at_omega_k_t;
 
 
complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
        complex double *target,
        const qpms_scatsys_at_omega_k_t *sswk
        );
 
complex double *qpms_scatsys_periodic_build_translation_matrix_full(
        complex double *target,
        const qpms_scatsys_t *ss,
        complex double wavenumber, 
        const cart3_t *wavevector, 
        double eta 
        );
 
complex double *qpms_scatsyswk_build_translation_matrix_full(
        complex double *target,
        const qpms_scatsys_at_omega_k_t *sswk
        );
 
 
qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
        complex double *target, 
        int *target_piv, 
        const qpms_scatsys_at_omega_k_t *sswk
        );
 
 
struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
        const qpms_scatsys_t *ss, 
        const double k[3],
        complex double omega_centre, 
        double omega_rr, 
        double omega_ri, 
        size_t contour_npoints, 
        double rank_tol, 
        size_t rank_sel_min, 
        double res_tol 
        );
 
// ======================= Periodic system -only related stuff end =========================
 
 
qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path);
 
qpms_scatsys_t *qpms_scatsys_load(char *path);
 
struct qpms_finite_group_t;
 
 
complex double *qpms_orbit_action_matrix(
 
        complex double *target,
        const qpms_ss_orbit_type_t *orbit,
        const qpms_vswf_set_spec_t *bspec,
        const struct qpms_finite_group_t *sym,
        const qpms_gmi_t g);
 
 
complex double *qpms_orbit_irrep_projector_matrix(
 
        complex double *target,
        const qpms_ss_orbit_type_t *orbit,
        const qpms_vswf_set_spec_t *bspec,
        const struct qpms_finite_group_t *sym,
        const qpms_iri_t iri);
 
complex double *qpms_orbit_irrep_basis(
        size_t *basis_size,
 
        complex double *target,
        const qpms_ss_orbit_type_t *orbit,
        const qpms_vswf_set_spec_t *bspec,
        const struct qpms_finite_group_t *sym,
        const qpms_iri_t iri);
 
 
 
complex double *qpms_scatsys_incident_field_vector_full(
 
        complex double *target_full,
        const qpms_scatsys_t *ss,
        qpms_incfield_t field_at_point,
        const void *args, 
        bool add 
        );
 
complex double *qpms_scatsysw_apply_Tmatrices_full(
        complex double *target_full, 
        const complex double *inc_full, 
        const qpms_scatsys_at_omega_t *ssw
        );
 
struct beyn_result_t; // See beyn.h for full definition
 
 
struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
        const qpms_scatsys_t *ss, 
        qpms_iri_t iri,
        complex double omega_centre, 
        double omega_rr, 
        double omega_ri, 
        size_t contour_npoints, 
        double rank_tol, 
        size_t rank_sel_min, 
        double res_tol 
        );
 
#if 0
 
struct beyn_result_t *qpms_scatsys_find_eigenmodes(
        const qpms_scatsys_t *ss,
        double eta, 
        const double *beta_, 
        complex double omega_centre, 
        double omega_rr, 
        double omega_ri, 
        size_t contour_npoints, 
        double rank_tol, 
        size_t rank_sel_min, 
        double res_tol 
        );
#endif
 
 
#if 0
 
complex double *qpms_scatsys_incident_field_vector_irrep_packed(
 
        complex double *target_full,
        const qpms_scatsys_t *ss,
        const qpms_iri_t iri, 
        qpms_incfield_t field_at_point,
        const void *args, 
        bool add 
        );
#endif
 
 
ccart3_t qpms_scatsys_scattered_E(
        const qpms_scatsys_t *ss,
        qpms_bessel_t typ, 
        complex double wavenumber, 
        const complex double *scatcoeff_full, 
        cart3_t evalpoint 
        );
 
 
ccart3_t qpms_scatsysw_scattered_E(
        const qpms_scatsys_at_omega_t *ssw,
        qpms_bessel_t typ, 
        const complex double *scatcoeff_full, 
        cart3_t evalpoint 
        );
 
 
ccart3_t qpms_scatsys_scattered_E__alt(
        const qpms_scatsys_t *ss,
        qpms_bessel_t typ, 
        complex double wavenumber, 
        const complex double *scatcoeff_full, 
        cart3_t evalpoint 
        );
 
 
ccart3_t qpms_scatsysw_scattered_E__alt(
        const qpms_scatsys_at_omega_t *ssw,
        qpms_bessel_t typ, 
        const complex double *scatcoeff_full, 
        cart3_t evalpoint 
        );
 
 
 
ccart3_t qpms_scatsyswk_scattered_E(
        const qpms_scatsys_at_omega_k_t *sswk,
        qpms_bessel_t typ, 
        const complex double *scatcoeff_full, 
        cart3_t evalpoint 
        );
 
 
qpms_errno_t qpms_scatsyswk_scattered_field_basis(
        ccart3_t *target, 
        const qpms_scatsys_at_omega_k_t *sswk,
        qpms_bessel_t typ, 
        cart3_t evalpoint 
        );
 
 
// TODO DOC
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]);
 
#if 0
ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
        qpms_iri_t iri, 
        const complex double *coeff_vector, 
        cart3_t where, 
        complex double k 
        );
#endif
 
#endif //QPMS_SCATSYSTEM_H