Move T-matrix -related C code to scatsystem.[ch]
Former-commit-id: e9a162c14fe5d91281d79d9af08014787cd7ed13
This commit is contained in:
Marek Nečada
parent
791a46b446
commit
2af9b83e35
|
|
@ -75,6 +75,10 @@ cdef extern from "qpms_types.h":
|
|||
ctypedef int qpms_gmi_t
|
||||
ctypedef int qpms_iri_t
|
||||
ctypedef const char * qpms_permutation_t
|
||||
struct qpms_tmatrix_t:
|
||||
qpms_vswf_set_spec_t *spec
|
||||
cdouble *m
|
||||
bint owns_m # FIXME in fact bool
|
||||
# maybe more if needed
|
||||
|
||||
cdef extern from "indexing.h":
|
||||
|
|
@ -221,17 +225,7 @@ cdef extern from "gsl/gsl_interp.h":
|
|||
const gsl_interp_type *gsl_interp_cspline
|
||||
# ^^^ These are probably the only relevant ones.
|
||||
|
||||
cdef extern from "scatsystem.h":
|
||||
struct qpms_tmatrix_t:
|
||||
qpms_vswf_set_spec_t *spec
|
||||
cdouble *m
|
||||
bint owns_m # FIXME in fact bool
|
||||
struct qpms_particle_t:
|
||||
cart3_t pos
|
||||
const qpms_tmatrix_t *tmatrix
|
||||
struct qpms_particle_tid_t:
|
||||
cart3_t pos
|
||||
qpms_ss_tmi_t tmatrix_id
|
||||
cdef extern from "tmatrices.h":
|
||||
struct qpms_tmatrix_interpolator_t:
|
||||
const qpms_vswf_set_spec_t *bspec
|
||||
void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp)
|
||||
|
|
@ -239,19 +233,6 @@ cdef extern from "scatsystem.h":
|
|||
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, double *freqs,
|
||||
const qpms_tmatrix_t *tmatrices_array, const gsl_interp_type *iptype)
|
||||
void qpms_tmatrix_free(qpms_tmatrix_t *tmatrix)
|
||||
struct qpms_scatsys_t:
|
||||
qpms_tmatrix_t **tm
|
||||
qpms_ss_tmi_t tm_count
|
||||
qpms_particle_tid_t *p
|
||||
qpms_ss_pi_t p_count
|
||||
# We shouldn't need more to construct a symmetric scatsystem ^^^
|
||||
size_t fecv_size
|
||||
size_t *saecv_sizes
|
||||
const qpms_finite_group_t *sym
|
||||
qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym)
|
||||
void qpms_scatsys_free(qpms_scatsys_t *s)
|
||||
qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI
|
||||
qpms_scatsys_t *qpms_scatsys_load(char *path) #NI
|
||||
qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
|
||||
const double rtol, const double atol)
|
||||
qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
|
||||
|
|
@ -277,6 +258,27 @@ cdef extern from "scatsystem.h":
|
|||
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,
|
||||
cdouble **tmdata)
|
||||
|
||||
cdef extern from "scatsystem.h":
|
||||
struct qpms_particle_t:
|
||||
cart3_t pos
|
||||
const qpms_tmatrix_t *tmatrix
|
||||
struct qpms_particle_tid_t:
|
||||
cart3_t pos
|
||||
qpms_ss_tmi_t tmatrix_id
|
||||
struct qpms_scatsys_t:
|
||||
qpms_tmatrix_t **tm
|
||||
qpms_ss_tmi_t tm_count
|
||||
qpms_particle_tid_t *p
|
||||
qpms_ss_pi_t p_count
|
||||
# We shouldn't need more to construct a symmetric scatsystem ^^^
|
||||
size_t fecv_size
|
||||
size_t *saecv_sizes
|
||||
const qpms_finite_group_t *sym
|
||||
qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym)
|
||||
void qpms_scatsys_free(qpms_scatsys_t *s)
|
||||
qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI
|
||||
qpms_scatsys_t *qpms_scatsys_load(char *path) #NI
|
||||
cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed,
|
||||
const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
|
||||
cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full,
|
||||
|
|
|
|||
|
|
@ -330,6 +330,24 @@ typedef int qpms_iri_t;
|
|||
*/
|
||||
typedef const char * qpms_permutation_t;
|
||||
|
||||
/// A T-matrix.
|
||||
/** In the future, I might rather use a more abstract approach in which T-matrix
|
||||
* is a mapping (function) of the field expansion coefficients.
|
||||
* So the interface might change.
|
||||
* For now, let me stick to the square dense matrix representation.
|
||||
*/
|
||||
typedef struct qpms_tmatrix_t {
|
||||
/** \brief VSWF basis specification, NOT owned by qpms_tmatrix_t by default.
|
||||
*
|
||||
* Usually not checked for meaningfulness by the functions (methods),
|
||||
* so the caller should take care that \a spec->ilist does not
|
||||
* contain any duplicities and that for each wave with order \a m
|
||||
* there is also one with order \a −m.
|
||||
*/
|
||||
const qpms_vswf_set_spec_t *spec;
|
||||
complex double *m; ///< Matrix elements in row-major order.
|
||||
bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
|
||||
} qpms_tmatrix_t;
|
||||
|
||||
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
|
||||
#endif // QPMS_TYPES
|
||||
|
|
|
|||
|
|
@ -6,7 +6,6 @@
|
|||
#include "vswf.h"
|
||||
#include "groups.h"
|
||||
#include "symmetries.h"
|
||||
#include <gsl/gsl_spline.h>
|
||||
#include <assert.h>
|
||||
#include <unistd.h>
|
||||
#include "vectors.h"
|
||||
|
|
@ -14,314 +13,13 @@
|
|||
#include <string.h>
|
||||
#include "qpms_error.h"
|
||||
#include "translations.h"
|
||||
#include "tmatrices.h"
|
||||
#include <pthread.h>
|
||||
|
||||
#define SQ(x) ((x)*(x))
|
||||
#define QPMS_SCATSYS_LEN_RTOL 1e-13
|
||||
#define QPMS_SCATSYS_TMATRIX_ATOL 1e-14
|
||||
#define QPMS_SCATSYS_TMATRIX_RTOL 1e-12
|
||||
qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
|
||||
qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t));
|
||||
if (!t) abort();
|
||||
else {
|
||||
t->spec = bspec;
|
||||
size_t n = bspec->n;
|
||||
t->m = calloc(n*n, sizeof(complex double));
|
||||
if (!t->m) abort();
|
||||
t->owns_m = true;
|
||||
}
|
||||
return t;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
size_t n = T->spec->n;
|
||||
for(size_t i = 0; i < n*n; ++i)
|
||||
t->m = T->m;
|
||||
return t;
|
||||
}
|
||||
|
||||
void qpms_tmatrix_free(qpms_tmatrix_t *t){
|
||||
if(t && t->owns_m) free(t->m);
|
||||
free(t);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
//qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
const size_t n = T->spec->n;
|
||||
complex double tmp[n][n];
|
||||
// tmp = M T
|
||||
const complex double one = 1, zero = 0;
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasNoTrans,
|
||||
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
||||
// t->m = tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &zero, T->m, n);
|
||||
return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
|
||||
const qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
const size_t n = T->spec->n;
|
||||
complex double tmp[n][n];
|
||||
// tmp = M T
|
||||
const complex double one = 1, zero = 0;
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasNoTrans,
|
||||
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
||||
// t->m = tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
|
||||
return t;
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
|
||||
complex double *tmdata, const size_t tmcount,
|
||||
const qpms_vswf_set_spec_t *bspec,
|
||||
const size_t n_symops, const qpms_irot3_t *symops) {
|
||||
const size_t n = bspec->n;
|
||||
qpms_tmatrix_t *tmcopy = qpms_tmatrix_init(bspec);
|
||||
complex double *symop_matrices = malloc(n*n*sizeof(complex double) * n_symops);
|
||||
if(!symop_matrices) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"malloc() failed.");
|
||||
for (size_t i = 0; i < n_symops; ++i)
|
||||
qpms_irot3_uvswfi_dense(symop_matrices + i*n*n, bspec, symops[i]);
|
||||
complex double tmp[n][n];
|
||||
const complex double one = 1, zero = 0;
|
||||
for (size_t tmi = 0; tmi < tmcount; ++tmi) {
|
||||
// Move the data in tmcopy; we will then write the sum directly into tmdata.
|
||||
memcpy(tmcopy->m, tmdata+n*n*tmi, n*n*sizeof(complex double));
|
||||
memset(tmdata+n*n*tmi, 0, n*n*sizeof(complex double));
|
||||
for (size_t i = 0; i < n_symops; ++i) {
|
||||
const complex double *const M = symop_matrices + i*n*n;
|
||||
// tmp = M T
|
||||
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
|
||||
n, n, n, &one, M, n, tmcopy->m, n, &zero, tmp, n);
|
||||
// tmdata[...] += tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &one, tmdata + tmi*n*n, n);
|
||||
}
|
||||
for (size_t ii = 0; ii < n*n; ++ii)
|
||||
tmdata[n*n*tmi+ii] /= n_symops;
|
||||
}
|
||||
free(symop_matrices);
|
||||
qpms_tmatrix_free(tmcopy);
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
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) {
|
||||
if (!(pointgroup->rep3d)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"This function requires pointgroup->rep3d to be set correctly!");
|
||||
return qpms_symmetrise_tmdata_irot3arr(tmdata, tmcount, bspec,
|
||||
pointgroup->order, pointgroup->rep3d);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
size_t n_symops,
|
||||
const qpms_irot3_t *symops
|
||||
) {
|
||||
if(qpms_symmetrise_tmdata_irot3arr(T->m, 1,
|
||||
T->spec, n_symops, symops) != QPMS_SUCCESS)
|
||||
return NULL;
|
||||
else return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const qpms_finite_group_t *pointgroup
|
||||
) {
|
||||
if(qpms_symmetrise_tmdata_finite_group(T->m, 1,
|
||||
T->spec, pointgroup) != QPMS_SUCCESS)
|
||||
return NULL;
|
||||
else return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_apply_symop(T, M);
|
||||
const size_t n = T->spec->n;
|
||||
for(size_t i = 0; i < n*n; ++i)
|
||||
T->m[i] = 0.5 * (t->m[i] + T->m[i]);
|
||||
qpms_tmatrix_free(t);
|
||||
return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
|
||||
const qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
const size_t n = T->spec->n;
|
||||
complex double tmp[n][n];
|
||||
// tmp = M T
|
||||
const complex double one = 1, zero = 0;
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasNoTrans,
|
||||
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
||||
// t->m = tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
|
||||
for(size_t i = 0; i < n*n; ++i)
|
||||
t->m[i] = 0.5 * (t->m[i] + T->m[i]);
|
||||
return t;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(const qpms_tmatrix_t *T) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
|
||||
return qpms_tmatrix_symmetrise_C_inf_inplace(t);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(qpms_tmatrix_t *T) {
|
||||
const size_t n = T->spec->n;
|
||||
for (size_t row = 0; row < n; row++) {
|
||||
qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
|
||||
for (size_t col = 0; col < n; col++) {
|
||||
qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
|
||||
if (rm == cm)
|
||||
;// No-op // t->m[n*row + col] = T->m[n*row + col];
|
||||
else
|
||||
T->m[n*row + col] = 0;
|
||||
}
|
||||
}
|
||||
return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(const qpms_tmatrix_t *T, int N) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
|
||||
return qpms_tmatrix_symmetrise_C_N_inplace(t, N);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) {
|
||||
const size_t n = T->spec->n;
|
||||
for (size_t row = 0; row < n; row++) {
|
||||
qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
|
||||
for (size_t col = 0; col < n; col++) {
|
||||
qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
|
||||
if (((rm - cm) % N) == 0)
|
||||
; // T->m[n*row + col] = T->m[n*row + col];
|
||||
else
|
||||
T->m[n*row + col] = 0;
|
||||
}
|
||||
}
|
||||
return T;
|
||||
}
|
||||
|
||||
bool qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
|
||||
const double rtol, const double atol)
|
||||
{
|
||||
if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec))
|
||||
return false;
|
||||
if (A->m == B->m)
|
||||
return true;
|
||||
const size_t n = A->spec->n;
|
||||
for (size_t i = 0; i < n*n; ++i) {
|
||||
const double tol = atol + rtol * (cabs(B->m[i]));
|
||||
if ( cabs(B->m[i] - A->m[i]) > tol )
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incount,
|
||||
const double *freqs, const qpms_tmatrix_t *ta, const gsl_interp_type *iptype//, const bool copy_bspec
|
||||
) {
|
||||
if (incount <= 0) return NULL;
|
||||
qpms_tmatrix_interpolator_t *ip = malloc(sizeof(qpms_tmatrix_interpolator_t));
|
||||
/*
|
||||
if (copy_bspec) {
|
||||
ip->bspec = qpms_vswf_set_spec_copy(ta[0].spec);
|
||||
ip->owns_bspec = true;
|
||||
}
|
||||
else {
|
||||
*/
|
||||
ip->bspec = ta[0].spec;
|
||||
// ip->owns_bspec = false;
|
||||
//}
|
||||
const size_t n = ip->bspec->n;
|
||||
|
||||
// check if all matrices have the same bspec
|
||||
for (size_t i = 0; i < incount; ++i)
|
||||
if (!qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec))
|
||||
abort();
|
||||
|
||||
if (!(ip->splines_real = calloc(n*n,sizeof(gsl_spline *)))) abort();
|
||||
if (!(ip->splines_imag = calloc(n*n,sizeof(gsl_spline *)))) abort();
|
||||
for (size_t row = 0; row < n; ++row)
|
||||
for (size_t col = 0; col < n; ++col) {
|
||||
double y_real[incount], y_imag[incount];
|
||||
bool n0_real = false, n0_imag = false;
|
||||
for (size_t i = 0; i < incount; ++i) {
|
||||
complex double telem = ta[i].m[n * row + col];
|
||||
if ((y_real[i] = creal(telem))) n0_real = true;
|
||||
if ((y_imag[i] = cimag(telem))) n0_imag = true;
|
||||
}
|
||||
if (n0_real) {
|
||||
gsl_spline *s =
|
||||
ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount);
|
||||
if (gsl_spline_init(s, freqs, y_real, incount) != 0 /*GSL_SUCCESS*/) abort();
|
||||
}
|
||||
else ip->splines_real[n * row + col] = NULL;
|
||||
if (n0_imag) {
|
||||
gsl_spline *s =
|
||||
ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount);
|
||||
if (gsl_spline_init(s, freqs, y_imag, incount) != 0 /*GSL_SUCCESS*/) abort();
|
||||
}
|
||||
else ip->splines_imag[n * row + col] = NULL;
|
||||
}
|
||||
return ip;
|
||||
}
|
||||
|
||||
void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *ip) {
|
||||
if (ip) {
|
||||
const size_t n = ip->bspec->n;
|
||||
for (size_t i = 0; i < n*n; ++i) {
|
||||
if (ip->splines_real[i]) gsl_spline_free(ip->splines_real[i]);
|
||||
if (ip->splines_imag[i]) gsl_spline_free(ip->splines_imag[i]);
|
||||
}
|
||||
//if (ip->owns_bspec)
|
||||
// qpms_vswf_set_spec_free(ip->bspec);
|
||||
free(ip);
|
||||
}
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *ip, double freq) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(ip->bspec);
|
||||
const size_t n = ip->bspec->n;
|
||||
for (size_t i = 0; i < n*n; ++i){
|
||||
if (ip->splines_real[i]) t->m[i] = gsl_spline_eval(ip->splines_real[i], freq, NULL /*does this work?*/);
|
||||
if (ip->splines_imag[i]) t->m[i] += I* gsl_spline_eval(ip->splines_imag[i], freq, NULL /*does this work?*/);
|
||||
}
|
||||
return t;
|
||||
}
|
||||
|
||||
|
||||
|
||||
// ------------ Stupid implementation of qpms_scatsys_apply_symmetry() -------------
|
||||
|
||||
|
|
|
|||
|
|
@ -12,44 +12,6 @@
|
|||
#define QPMS_SCATSYSTEM_H
|
||||
#include "qpms_types.h"
|
||||
#include <stdbool.h>
|
||||
#include <gsl/gsl_spline.h>
|
||||
#include <stdio.h> // only because of qpms_read_scuff_tmatrix()
|
||||
|
||||
/// A T-matrix.
|
||||
/** In the future, I might rather use a more abstract approach in which T-matrix
|
||||
* is a mapping (function) of the field expansion coefficients.
|
||||
* So the interface might change.
|
||||
* For now, let me stick to the square dense matrix representation.
|
||||
*/
|
||||
typedef struct qpms_tmatrix_t {
|
||||
/** \brief VSWF basis specification, NOT owned by qpms_tmatrix_t by default.
|
||||
*
|
||||
* Usually not checked for meaningfulness by the functions (methods),
|
||||
* so the caller should take care that \a spec->ilist does not
|
||||
* contain any duplicities and that for each wave with order \a m
|
||||
* there is also one with order \a −m.
|
||||
*/
|
||||
const qpms_vswf_set_spec_t *spec;
|
||||
complex double *m; ///< Matrix elements in row-major order.
|
||||
bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
|
||||
} qpms_tmatrix_t;
|
||||
|
||||
struct qpms_finite_group_t;
|
||||
typedef struct qpms_finite_group_t qpms_finite_group_t;
|
||||
|
||||
/// Returns a pointer to the beginning of the T-matrix row \a rowno.
|
||||
static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
|
||||
return t->m + (t->spec->n * rowno);
|
||||
}
|
||||
|
||||
/// Initialises a zero T-matrix.
|
||||
qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec);
|
||||
|
||||
/// Copies a T-matrix, allocating new array for the T-matrix data.
|
||||
qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T);
|
||||
|
||||
/// Destroys a T-matrix.
|
||||
void qpms_tmatrix_free(qpms_tmatrix_t *t);
|
||||
|
||||
/// A particle, defined by its T-matrix and position.
|
||||
typedef struct qpms_particle_t {
|
||||
|
|
@ -58,227 +20,8 @@ typedef struct qpms_particle_t {
|
|||
const qpms_tmatrix_t *tmatrix; ///< T-matrix; not owned by qpms_particle_t.
|
||||
} qpms_particle_t;
|
||||
|
||||
/// Check T-matrix equality/similarity.
|
||||
/**
|
||||
* This function actually checks for identical vswf specs.
|
||||
* TODO define constants with "default" atol, rtol for this function.
|
||||
*/
|
||||
bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
|
||||
const double rtol, const double atol);
|
||||
|
||||
/// Creates a T-matrix from another matrix and a symmetry operation.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = MTM^\dagger \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
|
||||
const qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = MTM^\dagger \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
|
||||
qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Symmetrizes a T-matrix with an involution symmetry operation.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = \frac{T + MTM^\dagger}{2} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
|
||||
const qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(
|
||||
const qpms_tmatrix_t *T ///< the original T-matrix
|
||||
);
|
||||
/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = \begin{cases}
|
||||
* T{}_{lm}^{\lambda\mu} & (m-\mu)\mod N=0\\
|
||||
* 0 & (m-\mu)\mod N\ne0
|
||||
* \end{cases} . \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
|
||||
const qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
int N ///< number of z-axis rotations in the group
|
||||
);
|
||||
|
||||
/// Symmetrizes a T-matrix with an involution symmetry operation, rewriting the original one.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = \frac{T + MTM^\dagger}{2} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
|
||||
qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(
|
||||
qpms_tmatrix_t *T ///< the original T-matrix
|
||||
);
|
||||
/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix, rewriting the original one.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = \begin{cases}
|
||||
* T{}_{lm}^{\lambda\mu} & (m-\mu)\mod N=0\\
|
||||
* 0 & (m-\mu)\mod N\ne0
|
||||
* \end{cases} . \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(
|
||||
qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
int N ///< number of z-axis rotations in the group
|
||||
);
|
||||
|
||||
/// Reads an open scuff-tmatrix generated file.
|
||||
/**
|
||||
* \a *freqs, \a *freqs_su, \a *tmatrices_array and \a *tmdata
|
||||
* arrays are allocated by this function
|
||||
* and have to be freed by the caller after use.
|
||||
* \a freqs_su and \a tmatrices_array can be NULL, in that case
|
||||
* the respective arrays are not filled nor allocated.
|
||||
*
|
||||
* The contents of tmatrices_array is NOT
|
||||
* supposed to be freed element per element.
|
||||
*
|
||||
* TODO more checks and options regarding NANs etc.
|
||||
*
|
||||
*/
|
||||
qpms_errno_t qpms_load_scuff_tmatrix(
|
||||
const char *path, ///< Path to the TMatrix file
|
||||
const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec
|
||||
size_t *n, ///< Number of successfully loaded t-matrices
|
||||
double **freqs, ///< Frequencies in SI units..
|
||||
double **freqs_su, ///< Frequencies in SCUFF units (optional).
|
||||
/// The resulting T-matrices (optional).
|
||||
qpms_tmatrix_t **tmatrices_array,
|
||||
complex double **tmdata ///< The T-matrices raw contents
|
||||
);
|
||||
/// Loads a scuff-tmatrix generated file.
|
||||
/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a
|
||||
* path instead of open FILE.
|
||||
*/
|
||||
qpms_errno_t qpms_read_scuff_tmatrix(
|
||||
FILE *f, ///< An open stream with the T-matrix data.
|
||||
const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec
|
||||
size_t *n, ///< Number of successfully loaded t-matrices
|
||||
double **freqs, ///< Frequencies in SI units.
|
||||
double **freqs_su, ///< Frequencies in SCUFF units (optional).
|
||||
/// The resulting T-matrices (optional).
|
||||
qpms_tmatrix_t **tmatrices_array,
|
||||
/// The T-matrices raw contents.
|
||||
/** The coefficient of outgoing wave defined by
|
||||
* \a bspec->ilist[desti] as a result of incoming wave
|
||||
* \a bspec->ilist[srci] at frequency \a (*freqs)[fi]
|
||||
* is accessed via
|
||||
* (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci].
|
||||
*/
|
||||
complex double ** tmdata
|
||||
);
|
||||
|
||||
/// In-place application of point group elements on raw T-matrix data.
|
||||
/** \a tmdata can be e.g. obtained by qpms_load_scuff_tmatrix().
|
||||
* The \a symops array should always contain all elements of a finite
|
||||
* point (sub)group, including the identity operation.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
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
|
||||
);
|
||||
|
||||
/// In-place application of a point group on raw T-matrix data.
|
||||
/** This does the same as qpms_symmetrise_tmdata_irot3arr(),
|
||||
* but takes a valid finite point group as an argument.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
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
|
||||
);
|
||||
|
||||
/// In-place application of point group elements on a T-matrix.
|
||||
/** The \a symops array should always contain all elements of a finite
|
||||
* point (sub)group, including the identity operation.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
size_t n_symops,
|
||||
const qpms_irot3_t *symops
|
||||
);
|
||||
|
||||
/// In-place application of point group elements on a T-matrix.
|
||||
/** This does the same as qpms_tmatrix_symmetrise_irot3arr(),
|
||||
* but takes a valid finite point group as an argument.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const qpms_finite_group_t *pointgroup
|
||||
);
|
||||
|
||||
/* 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; ///< There will be a spline object for each nonzero element
|
||||
gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element
|
||||
// gsl_interp_accel **accel_real;
|
||||
// gsl_interp_accel **accel_imag;
|
||||
} qpms_tmatrix_interpolator_t;
|
||||
|
||||
/// Free a T-matrix interpolator.
|
||||
void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp);
|
||||
|
||||
/// Evaluate a T-matrix interpolated value.
|
||||
/** The result is to be freed using qpms_tmatrix_free().*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq);
|
||||
|
||||
/// Create a T-matrix interpolator from frequency and T-matrix arrays.
|
||||
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided.
|
||||
const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers!
|
||||
const gsl_interp_type *iptype
|
||||
//, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
|
||||
/*, ...? */);
|
||||
|
||||
struct qpms_finite_group_t;
|
||||
typedef struct qpms_finite_group_t qpms_finite_group_t;
|
||||
|
||||
/// A particle, defined by its T-matrix INDEX and position, to be used in qpms_scatsys_t.
|
||||
typedef struct qpms_particle_tid_t {
|
||||
|
|
@ -571,41 +314,4 @@ complex double *qpms_orbit_irrep_basis(
|
|||
/// The index of the irreducible representation of sym.
|
||||
const qpms_iri_t iri);
|
||||
|
||||
|
||||
|
||||
#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; ///< Maybe I want some usual integer type instead.
|
||||
|
||||
///Abstract T-matrix type draft.
|
||||
/**
|
||||
* TODO.
|
||||
*/
|
||||
typedef struct qpms_abstract_tmatrix_t{
|
||||
qpms_tmatrix_id_t id;
|
||||
/// Generators of the discrete point group under which T-matrix is invariant.
|
||||
qpms_irot3_t *invar_gens;
|
||||
/// Length of invar_gens.
|
||||
qpms_gmi_t invar_gens_size;
|
||||
|
||||
} qpms_abstract_tmatrix_t;
|
||||
|
||||
|
||||
typedef struct qpms_abstract_particle_t{
|
||||
} qpms_abstract_particle_t;
|
||||
|
||||
/// An abstract particle, defined by its position and abstract T-matrix.
|
||||
typedef struct qpms_abstract_particle_t {
|
||||
cart3_t pos; ///< Particle position in cartesian coordinates.
|
||||
const qpms_abstract_tmatrix_t *tmatrix; ///< T-matrix; not owned by this.
|
||||
} qpms_abstract_particle_t;
|
||||
|
||||
|
||||
/** This is just an alias, as the same index can be used for
|
||||
* abstract T-matrices as well.
|
||||
*/
|
||||
typedef qpms_particle_tid_t qpms_abstract_particle_tid_t;
|
||||
#endif // 0
|
||||
#endif //QPMS_SCATSYSTEM_H
|
||||
|
|
|
|||
|
|
@ -0,0 +1,456 @@
|
|||
#define _POSIX_C_SOURCE 200809L // for getline()
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <cblas.h>
|
||||
#include <unistd.h>
|
||||
#include "scatsystem.h"
|
||||
#include "indexing.h"
|
||||
#include "vswf.h"
|
||||
#include "groups.h"
|
||||
#include "symmetries.h"
|
||||
#include <gsl/gsl_spline.h>
|
||||
#include <assert.h>
|
||||
#include <unistd.h>
|
||||
#include "vectors.h"
|
||||
#include "wigner.h"
|
||||
#include <string.h>
|
||||
#include "qpms_error.h"
|
||||
#include "tmatrices.h"
|
||||
|
||||
#define HBAR (1.05457162825e-34)
|
||||
#define ELECTRONVOLT (1.602176487e-19)
|
||||
#define SCUFF_OMEGAUNIT (3e14)
|
||||
|
||||
#define SQ(x) ((x)*(x))
|
||||
qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
|
||||
qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t));
|
||||
if (!t) abort();
|
||||
else {
|
||||
t->spec = bspec;
|
||||
size_t n = bspec->n;
|
||||
t->m = calloc(n*n, sizeof(complex double));
|
||||
if (!t->m) abort();
|
||||
t->owns_m = true;
|
||||
}
|
||||
return t;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
size_t n = T->spec->n;
|
||||
for(size_t i = 0; i < n*n; ++i)
|
||||
t->m = T->m;
|
||||
return t;
|
||||
}
|
||||
|
||||
void qpms_tmatrix_free(qpms_tmatrix_t *t){
|
||||
if(t && t->owns_m) free(t->m);
|
||||
free(t);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
//qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
const size_t n = T->spec->n;
|
||||
complex double tmp[n][n];
|
||||
// tmp = M T
|
||||
const complex double one = 1, zero = 0;
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasNoTrans,
|
||||
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
||||
// t->m = tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &zero, T->m, n);
|
||||
return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
|
||||
const qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
const size_t n = T->spec->n;
|
||||
complex double tmp[n][n];
|
||||
// tmp = M T
|
||||
const complex double one = 1, zero = 0;
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasNoTrans,
|
||||
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
||||
// t->m = tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
|
||||
return t;
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
|
||||
complex double *tmdata, const size_t tmcount,
|
||||
const qpms_vswf_set_spec_t *bspec,
|
||||
const size_t n_symops, const qpms_irot3_t *symops) {
|
||||
const size_t n = bspec->n;
|
||||
qpms_tmatrix_t *tmcopy = qpms_tmatrix_init(bspec);
|
||||
complex double *symop_matrices = malloc(n*n*sizeof(complex double) * n_symops);
|
||||
if(!symop_matrices) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"malloc() failed.");
|
||||
for (size_t i = 0; i < n_symops; ++i)
|
||||
qpms_irot3_uvswfi_dense(symop_matrices + i*n*n, bspec, symops[i]);
|
||||
complex double tmp[n][n];
|
||||
const complex double one = 1, zero = 0;
|
||||
for (size_t tmi = 0; tmi < tmcount; ++tmi) {
|
||||
// Move the data in tmcopy; we will then write the sum directly into tmdata.
|
||||
memcpy(tmcopy->m, tmdata+n*n*tmi, n*n*sizeof(complex double));
|
||||
memset(tmdata+n*n*tmi, 0, n*n*sizeof(complex double));
|
||||
for (size_t i = 0; i < n_symops; ++i) {
|
||||
const complex double *const M = symop_matrices + i*n*n;
|
||||
// tmp = M T
|
||||
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
|
||||
n, n, n, &one, M, n, tmcopy->m, n, &zero, tmp, n);
|
||||
// tmdata[...] += tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &one, tmdata + tmi*n*n, n);
|
||||
}
|
||||
for (size_t ii = 0; ii < n*n; ++ii)
|
||||
tmdata[n*n*tmi+ii] /= n_symops;
|
||||
}
|
||||
free(symop_matrices);
|
||||
qpms_tmatrix_free(tmcopy);
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
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) {
|
||||
if (!(pointgroup->rep3d)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"This function requires pointgroup->rep3d to be set correctly!");
|
||||
return qpms_symmetrise_tmdata_irot3arr(tmdata, tmcount, bspec,
|
||||
pointgroup->order, pointgroup->rep3d);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
size_t n_symops,
|
||||
const qpms_irot3_t *symops
|
||||
) {
|
||||
if(qpms_symmetrise_tmdata_irot3arr(T->m, 1,
|
||||
T->spec, n_symops, symops) != QPMS_SUCCESS)
|
||||
return NULL;
|
||||
else return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const qpms_finite_group_t *pointgroup
|
||||
) {
|
||||
if(qpms_symmetrise_tmdata_finite_group(T->m, 1,
|
||||
T->spec, pointgroup) != QPMS_SUCCESS)
|
||||
return NULL;
|
||||
else return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_apply_symop(T, M);
|
||||
const size_t n = T->spec->n;
|
||||
for(size_t i = 0; i < n*n; ++i)
|
||||
T->m[i] = 0.5 * (t->m[i] + T->m[i]);
|
||||
qpms_tmatrix_free(t);
|
||||
return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
|
||||
const qpms_tmatrix_t *T,
|
||||
const complex double *M
|
||||
)
|
||||
{
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
||||
const size_t n = T->spec->n;
|
||||
complex double tmp[n][n];
|
||||
// tmp = M T
|
||||
const complex double one = 1, zero = 0;
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasNoTrans,
|
||||
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
||||
// t->m = tmp M* = M T M*
|
||||
cblas_zgemm(CblasRowMajor,
|
||||
CblasNoTrans,
|
||||
CblasConjTrans,
|
||||
n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
|
||||
for(size_t i = 0; i < n*n; ++i)
|
||||
t->m[i] = 0.5 * (t->m[i] + T->m[i]);
|
||||
return t;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(const qpms_tmatrix_t *T) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
|
||||
return qpms_tmatrix_symmetrise_C_inf_inplace(t);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(qpms_tmatrix_t *T) {
|
||||
const size_t n = T->spec->n;
|
||||
for (size_t row = 0; row < n; row++) {
|
||||
qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
|
||||
for (size_t col = 0; col < n; col++) {
|
||||
qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
|
||||
if (rm == cm)
|
||||
;// No-op // t->m[n*row + col] = T->m[n*row + col];
|
||||
else
|
||||
T->m[n*row + col] = 0;
|
||||
}
|
||||
}
|
||||
return T;
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(const qpms_tmatrix_t *T, int N) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
|
||||
return qpms_tmatrix_symmetrise_C_N_inplace(t, N);
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) {
|
||||
const size_t n = T->spec->n;
|
||||
for (size_t row = 0; row < n; row++) {
|
||||
qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
|
||||
for (size_t col = 0; col < n; col++) {
|
||||
qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
|
||||
if (((rm - cm) % N) == 0)
|
||||
; // T->m[n*row + col] = T->m[n*row + col];
|
||||
else
|
||||
T->m[n*row + col] = 0;
|
||||
}
|
||||
}
|
||||
return T;
|
||||
}
|
||||
|
||||
bool qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
|
||||
const double rtol, const double atol)
|
||||
{
|
||||
if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec))
|
||||
return false;
|
||||
if (A->m == B->m)
|
||||
return true;
|
||||
const size_t n = A->spec->n;
|
||||
for (size_t i = 0; i < n*n; ++i) {
|
||||
const double tol = atol + rtol * (cabs(B->m[i]));
|
||||
if ( cabs(B->m[i] - A->m[i]) > tol )
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incount,
|
||||
const double *freqs, const qpms_tmatrix_t *ta, const gsl_interp_type *iptype//, const bool copy_bspec
|
||||
) {
|
||||
if (incount <= 0) return NULL;
|
||||
qpms_tmatrix_interpolator_t *ip = malloc(sizeof(qpms_tmatrix_interpolator_t));
|
||||
/*
|
||||
if (copy_bspec) {
|
||||
ip->bspec = qpms_vswf_set_spec_copy(ta[0].spec);
|
||||
ip->owns_bspec = true;
|
||||
}
|
||||
else {
|
||||
*/
|
||||
ip->bspec = ta[0].spec;
|
||||
// ip->owns_bspec = false;
|
||||
//}
|
||||
const size_t n = ip->bspec->n;
|
||||
|
||||
// check if all matrices have the same bspec
|
||||
for (size_t i = 0; i < incount; ++i)
|
||||
if (!qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec))
|
||||
abort();
|
||||
|
||||
if (!(ip->splines_real = calloc(n*n,sizeof(gsl_spline *)))) abort();
|
||||
if (!(ip->splines_imag = calloc(n*n,sizeof(gsl_spline *)))) abort();
|
||||
for (size_t row = 0; row < n; ++row)
|
||||
for (size_t col = 0; col < n; ++col) {
|
||||
double y_real[incount], y_imag[incount];
|
||||
bool n0_real = false, n0_imag = false;
|
||||
for (size_t i = 0; i < incount; ++i) {
|
||||
complex double telem = ta[i].m[n * row + col];
|
||||
if ((y_real[i] = creal(telem))) n0_real = true;
|
||||
if ((y_imag[i] = cimag(telem))) n0_imag = true;
|
||||
}
|
||||
if (n0_real) {
|
||||
gsl_spline *s =
|
||||
ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount);
|
||||
if (gsl_spline_init(s, freqs, y_real, incount) != 0 /*GSL_SUCCESS*/) abort();
|
||||
}
|
||||
else ip->splines_real[n * row + col] = NULL;
|
||||
if (n0_imag) {
|
||||
gsl_spline *s =
|
||||
ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount);
|
||||
if (gsl_spline_init(s, freqs, y_imag, incount) != 0 /*GSL_SUCCESS*/) abort();
|
||||
}
|
||||
else ip->splines_imag[n * row + col] = NULL;
|
||||
}
|
||||
return ip;
|
||||
}
|
||||
|
||||
void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *ip) {
|
||||
if (ip) {
|
||||
const size_t n = ip->bspec->n;
|
||||
for (size_t i = 0; i < n*n; ++i) {
|
||||
if (ip->splines_real[i]) gsl_spline_free(ip->splines_real[i]);
|
||||
if (ip->splines_imag[i]) gsl_spline_free(ip->splines_imag[i]);
|
||||
}
|
||||
//if (ip->owns_bspec)
|
||||
// qpms_vswf_set_spec_free(ip->bspec);
|
||||
free(ip);
|
||||
}
|
||||
}
|
||||
|
||||
qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *ip, double freq) {
|
||||
qpms_tmatrix_t *t = qpms_tmatrix_init(ip->bspec);
|
||||
const size_t n = ip->bspec->n;
|
||||
for (size_t i = 0; i < n*n; ++i){
|
||||
if (ip->splines_real[i]) t->m[i] = gsl_spline_eval(ip->splines_real[i], freq, NULL /*does this work?*/);
|
||||
if (ip->splines_imag[i]) t->m[i] += I* gsl_spline_eval(ip->splines_imag[i], freq, NULL /*does this work?*/);
|
||||
}
|
||||
return t;
|
||||
}
|
||||
|
||||
|
||||
|
||||
double qpms_SU2eV(double e_SU) {
|
||||
return e_SU * SCUFF_OMEGAUNIT / (ELECTRONVOLT / HBAR);
|
||||
}
|
||||
|
||||
double qpms_SU2SI(double e_SU) {
|
||||
return e_SU * SCUFF_OMEGAUNIT;
|
||||
}
|
||||
|
||||
/// TODO doc and more checks
|
||||
qpms_errno_t qpms_read_scuff_tmatrix(
|
||||
FILE *f, ///< file handle
|
||||
const qpms_vswf_set_spec_t * bs, ///< VSWF set spec
|
||||
size_t *const n, ///< Number of successfully loaded t-matrices
|
||||
double* *const freqs, ///< Frequencies in SI units
|
||||
double* *const freqs_su, ///< Frequencies in SCUFF units (optional)
|
||||
qpms_tmatrix_t* *const tmatrices_array, ///< The resulting T-matrices (optional).
|
||||
complex double* *const tmdata
|
||||
) {
|
||||
if (!(freqs && n && tmdata))
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"freqs, n, and tmdata are mandatory arguments and must not be NULL.");
|
||||
if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT?
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)"
|
||||
" norm supported right now.");
|
||||
int n_alloc = 128; // First chunk to allocate
|
||||
*n = 0;
|
||||
*freqs = malloc(n_alloc * sizeof(double));
|
||||
if (freqs_su) *freqs_su = malloc(n_alloc * sizeof(double));
|
||||
*tmdata = malloc(sizeof(complex double) * bs->n * bs->n * n_alloc);
|
||||
if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"malloc() failed.");
|
||||
size_t linebufsz = 256;
|
||||
char *linebuf = malloc(linebufsz);
|
||||
ssize_t readchars;
|
||||
double lastfreq_su = NAN;
|
||||
while((readchars = getline(&linebuf, &linebufsz, f)) != -1) {
|
||||
if (linebuf[0] == '#') continue;
|
||||
int Alpha, LAlpha, MAlpha, PAlpha, Beta, LBeta, MBeta, PBeta;
|
||||
double currentfreq_su, tr, ti;
|
||||
if (11 != sscanf(linebuf, "%lf %d %d %d %d %d %d %d %d %lf %lf",
|
||||
¤tfreq_su, &Alpha, &LAlpha, &MAlpha, &PAlpha,
|
||||
&Beta, &LBeta, &MBeta, &PBeta, &tr, &ti))
|
||||
abort(); // Malformed T-matrix file
|
||||
if (currentfreq_su != lastfreq_su) { // New frequency -> new T-matrix
|
||||
++*n;
|
||||
lastfreq_su = currentfreq_su;
|
||||
if(*n > n_alloc) {
|
||||
n_alloc *= 2;
|
||||
*freqs = realloc(*freqs, n_alloc * sizeof(double));
|
||||
if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double));
|
||||
*tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc);
|
||||
if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"realloc() failed.");
|
||||
}
|
||||
if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su;
|
||||
(*freqs)[*n-1] = qpms_SU2SI(currentfreq_su);
|
||||
|
||||
for(size_t i = 0; i < bs->n * bs->n; ++i)
|
||||
(*tmdata)[(*n-1)*bs->n*bs->n + i] = NAN + I*NAN;
|
||||
}
|
||||
qpms_vswf_type_t TAlpha, TBeta;
|
||||
switch(PAlpha) {
|
||||
case 0: TAlpha = QPMS_VSWF_MAGNETIC; break;
|
||||
case 1: TAlpha = QPMS_VSWF_ELECTRIC; break;
|
||||
default: assert(0);
|
||||
}
|
||||
switch(PBeta) {
|
||||
case 0: TBeta = QPMS_VSWF_MAGNETIC; break;
|
||||
case 1: TBeta = QPMS_VSWF_ELECTRIC; break;
|
||||
default: assert(0);
|
||||
}
|
||||
qpms_uvswfi_t srcui = qpms_tmn2uvswfi(TAlpha, MAlpha, LAlpha),
|
||||
destui = qpms_tmn2uvswfi(TBeta, MBeta, LBeta);
|
||||
ssize_t srci = qpms_vswf_set_spec_find_uvswfi(bs, srcui),
|
||||
desti = qpms_vswf_set_spec_find_uvswfi(bs, destui);
|
||||
if (srci == -1 || desti == -1)
|
||||
/* This element has not been requested in bs->ilist. */
|
||||
continue;
|
||||
else
|
||||
(*tmdata)[(*n-1)*bs->n*bs->n + desti*bs->n + srci] = tr + I*ti;
|
||||
}
|
||||
free(linebuf);
|
||||
// free some more memory
|
||||
n_alloc = *n;
|
||||
*freqs = realloc(*freqs, n_alloc * sizeof(double));
|
||||
if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double));
|
||||
if (tmatrices_array) *tmatrices_array = realloc(*tmatrices_array, n_alloc * sizeof(qpms_tmatrix_t));
|
||||
*tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc);
|
||||
if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"realloc() failed.");
|
||||
if (tmatrices_array) {
|
||||
*tmatrices_array = malloc(n_alloc * sizeof(qpms_tmatrix_t));
|
||||
if (!*tmatrices_array)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"malloc() failed.");
|
||||
for (size_t i = 0; i < *n; ++i) {
|
||||
(*tmatrices_array)[i].spec = bs;
|
||||
(*tmatrices_array)[i].m = (*tmdata) + i * bs->n * bs->n;
|
||||
(*tmatrices_array)[i].owns_m = false;
|
||||
}
|
||||
}
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_load_scuff_tmatrix(
|
||||
const char *path, ///< file path
|
||||
const qpms_vswf_set_spec_t * bs, ///< VSWF set spec
|
||||
size_t *const n, ///< Number of successfully loaded t-matrices
|
||||
double **const freqs, ///< Frequencies in SI units
|
||||
double ** const freqs_su, ///< Frequencies in SCUFF units (optional)
|
||||
qpms_tmatrix_t ** const tmatrices_array, ///< The resulting T-matrices (optional).
|
||||
complex double ** const tmdata
|
||||
) {
|
||||
FILE *f = fopen(path, "r");
|
||||
if (!f)
|
||||
qpms_pr_error_at_line(__FILE__, __LINE__, __func__,
|
||||
"Could not open T-matrix file %s", path);
|
||||
qpms_errno_t retval =
|
||||
qpms_read_scuff_tmatrix(f, bs, n, freqs, freqs_su, tmatrices_array, tmdata);
|
||||
if(fclose(f)) qpms_pr_error_at_line(__FILE__, __LINE__, __func__,
|
||||
"Could not close the T-matrix file %s (well, that's weird, "
|
||||
"since it's read only).", path);
|
||||
|
||||
return retval;
|
||||
}
|
||||
|
||||
286
qpms/tmatrices.h
286
qpms/tmatrices.h
|
|
@ -1,12 +1,284 @@
|
|||
/* \file tmatrices.h
|
||||
* \brief T-matrices for scattering systems.
|
||||
*/
|
||||
#ifndef TMATRICES_H
|
||||
#define TMATRICES_H
|
||||
/* TODO
|
||||
* This file will contain declarations of functions providing
|
||||
* a) Mie T-matrix for spherical particle
|
||||
* i) using Drude model
|
||||
* ii) using interpolated material data
|
||||
* b) T-matrix from scuff-tmatrix output, using interpolation
|
||||
*/
|
||||
#include "qpms_types.h"
|
||||
#include <gsl/gsl_spline.h>
|
||||
|
||||
struct qpms_finite_group_t;
|
||||
typedef struct qpms_finite_group_t qpms_finite_group_t;
|
||||
|
||||
/// Returns a pointer to the beginning of the T-matrix row \a rowno.
|
||||
static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
|
||||
return t->m + (t->spec->n * rowno);
|
||||
}
|
||||
|
||||
/// Initialises a zero T-matrix.
|
||||
qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec);
|
||||
|
||||
/// Copies a T-matrix, allocating new array for the T-matrix data.
|
||||
qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T);
|
||||
|
||||
/// Destroys a T-matrix.
|
||||
void qpms_tmatrix_free(qpms_tmatrix_t *t);
|
||||
|
||||
/// Check T-matrix equality/similarity.
|
||||
/**
|
||||
* This function actually checks for identical vswf specs.
|
||||
* TODO define constants with "default" atol, rtol for this function.
|
||||
*/
|
||||
bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
|
||||
const double rtol, const double atol);
|
||||
|
||||
/// Creates a T-matrix from another matrix and a symmetry operation.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = MTM^\dagger \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
|
||||
const qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = MTM^\dagger \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
|
||||
qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Symmetrizes a T-matrix with an involution symmetry operation.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = \frac{T + MTM^\dagger}{2} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
|
||||
const qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(
|
||||
const qpms_tmatrix_t *T ///< the original T-matrix
|
||||
);
|
||||
/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = \begin{cases}
|
||||
* T{}_{lm}^{\lambda\mu} & (m-\mu)\mod N=0\\
|
||||
* 0 & (m-\mu)\mod N\ne0
|
||||
* \end{cases} . \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
|
||||
const qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
int N ///< number of z-axis rotations in the group
|
||||
);
|
||||
|
||||
/// Symmetrizes a T-matrix with an involution symmetry operation, rewriting the original one.
|
||||
/** The symmetry operation is expected to be a unitary (square)
|
||||
* matrix \a M with the same
|
||||
* VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
|
||||
* correspond to CHECKME \f[ T' = \frac{T + MTM^\dagger}{2} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
|
||||
qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
const complex double *M ///< the symmetry op matrix in row-major format
|
||||
);
|
||||
|
||||
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(
|
||||
qpms_tmatrix_t *T ///< the original T-matrix
|
||||
);
|
||||
/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix, rewriting the original one.
|
||||
/**
|
||||
* \f[ {T'}_{tlm}^{\tau\lambda\mu} = \begin{cases}
|
||||
* T{}_{lm}^{\lambda\mu} & (m-\mu)\mod N=0\\
|
||||
* 0 & (m-\mu)\mod N\ne0
|
||||
* \end{cases} . \f]
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(
|
||||
qpms_tmatrix_t *T, ///< the original T-matrix
|
||||
int N ///< number of z-axis rotations in the group
|
||||
);
|
||||
|
||||
/// Reads an open scuff-tmatrix generated file.
|
||||
/**
|
||||
* \a *freqs, \a *freqs_su, \a *tmatrices_array and \a *tmdata
|
||||
* arrays are allocated by this function
|
||||
* and have to be freed by the caller after use.
|
||||
* \a freqs_su and \a tmatrices_array can be NULL, in that case
|
||||
* the respective arrays are not filled nor allocated.
|
||||
*
|
||||
* The contents of tmatrices_array is NOT
|
||||
* supposed to be freed element per element.
|
||||
*
|
||||
* TODO more checks and options regarding NANs etc.
|
||||
*
|
||||
*/
|
||||
qpms_errno_t qpms_load_scuff_tmatrix(
|
||||
const char *path, ///< Path to the TMatrix file
|
||||
const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec
|
||||
size_t *n, ///< Number of successfully loaded t-matrices
|
||||
double **freqs, ///< Frequencies in SI units..
|
||||
double **freqs_su, ///< Frequencies in SCUFF units (optional).
|
||||
/// The resulting T-matrices (optional).
|
||||
qpms_tmatrix_t **tmatrices_array,
|
||||
complex double **tmdata ///< The T-matrices raw contents
|
||||
);
|
||||
/// Loads a scuff-tmatrix generated file.
|
||||
/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a
|
||||
* path instead of open FILE.
|
||||
*/
|
||||
qpms_errno_t qpms_read_scuff_tmatrix(
|
||||
FILE *f, ///< An open stream with the T-matrix data.
|
||||
const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec
|
||||
size_t *n, ///< Number of successfully loaded t-matrices
|
||||
double **freqs, ///< Frequencies in SI units.
|
||||
double **freqs_su, ///< Frequencies in SCUFF units (optional).
|
||||
/// The resulting T-matrices (optional).
|
||||
qpms_tmatrix_t **tmatrices_array,
|
||||
/// The T-matrices raw contents.
|
||||
/** The coefficient of outgoing wave defined by
|
||||
* \a bspec->ilist[desti] as a result of incoming wave
|
||||
* \a bspec->ilist[srci] at frequency \a (*freqs)[fi]
|
||||
* is accessed via
|
||||
* (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci].
|
||||
*/
|
||||
complex double ** tmdata
|
||||
);
|
||||
|
||||
/// In-place application of point group elements on raw T-matrix data.
|
||||
/** \a tmdata can be e.g. obtained by qpms_load_scuff_tmatrix().
|
||||
* The \a symops array should always contain all elements of a finite
|
||||
* point (sub)group, including the identity operation.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
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
|
||||
);
|
||||
|
||||
/// In-place application of a point group on raw T-matrix data.
|
||||
/** This does the same as qpms_symmetrise_tmdata_irot3arr(),
|
||||
* but takes a valid finite point group as an argument.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
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
|
||||
);
|
||||
|
||||
/// In-place application of point group elements on a T-matrix.
|
||||
/** The \a symops array should always contain all elements of a finite
|
||||
* point (sub)group, including the identity operation.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
size_t n_symops,
|
||||
const qpms_irot3_t *symops
|
||||
);
|
||||
|
||||
/// In-place application of point group elements on a T-matrix.
|
||||
/** This does the same as qpms_tmatrix_symmetrise_irot3arr(),
|
||||
* but takes a valid finite point group as an argument.
|
||||
*
|
||||
* TODO more doc.
|
||||
*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
|
||||
qpms_tmatrix_t *T,
|
||||
const qpms_finite_group_t *pointgroup
|
||||
);
|
||||
|
||||
/* 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; ///< There will be a spline object for each nonzero element
|
||||
gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element
|
||||
// gsl_interp_accel **accel_real;
|
||||
// gsl_interp_accel **accel_imag;
|
||||
} qpms_tmatrix_interpolator_t;
|
||||
|
||||
/// Free a T-matrix interpolator.
|
||||
void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp);
|
||||
|
||||
/// Evaluate a T-matrix interpolated value.
|
||||
/** The result is to be freed using qpms_tmatrix_free().*/
|
||||
qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq);
|
||||
|
||||
/// Create a T-matrix interpolator from frequency and T-matrix arrays.
|
||||
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided.
|
||||
const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers!
|
||||
const gsl_interp_type *iptype
|
||||
//, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
|
||||
/*, ...? */);
|
||||
|
||||
|
||||
#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; ///< Maybe I want some usual integer type instead.
|
||||
|
||||
///Abstract T-matrix type draft.
|
||||
/**
|
||||
* TODO.
|
||||
*/
|
||||
typedef struct qpms_abstract_tmatrix_t{
|
||||
qpms_tmatrix_id_t id;
|
||||
/// Generators of the discrete point group under which T-matrix is invariant.
|
||||
qpms_irot3_t *invar_gens;
|
||||
/// Length of invar_gens.
|
||||
qpms_gmi_t invar_gens_size;
|
||||
|
||||
} qpms_abstract_tmatrix_t;
|
||||
|
||||
|
||||
typedef struct qpms_abstract_particle_t{
|
||||
} qpms_abstract_particle_t;
|
||||
|
||||
/// An abstract particle, defined by its position and abstract T-matrix.
|
||||
typedef struct qpms_abstract_particle_t {
|
||||
cart3_t pos; ///< Particle position in cartesian coordinates.
|
||||
const qpms_abstract_tmatrix_t *tmatrix; ///< T-matrix; not owned by this.
|
||||
} qpms_abstract_particle_t;
|
||||
|
||||
|
||||
/** This is just an alias, as the same index can be used for
|
||||
* abstract T-matrices as well.
|
||||
*/
|
||||
typedef qpms_particle_tid_t qpms_abstract_particle_tid_t;
|
||||
#endif // 0
|
||||
|
||||
#endif //TMATRICES_H
|
||||
|
|
|
|||
|
|
@ -1,145 +0,0 @@
|
|||
#define _POSIX_C_SOURCE 200809L // for getline()
|
||||
#include "scatsystem.h"
|
||||
#include <unistd.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include "indexing.h"
|
||||
#include "vswf.h" // qpms_vswf_set_spec_find_uvswfi()
|
||||
#include "qpms_error.h"
|
||||
|
||||
|
||||
#define HBAR (1.05457162825e-34)
|
||||
#define ELECTRONVOLT (1.602176487e-19)
|
||||
#define SCUFF_OMEGAUNIT (3e14)
|
||||
|
||||
double qpms_SU2eV(double e_SU) {
|
||||
return e_SU * SCUFF_OMEGAUNIT / (ELECTRONVOLT / HBAR);
|
||||
}
|
||||
|
||||
double qpms_SU2SI(double e_SU) {
|
||||
return e_SU * SCUFF_OMEGAUNIT;
|
||||
}
|
||||
|
||||
/// TODO doc and more checks
|
||||
qpms_errno_t qpms_read_scuff_tmatrix(
|
||||
FILE *f, ///< file handle
|
||||
const qpms_vswf_set_spec_t * bs, ///< VSWF set spec
|
||||
size_t *const n, ///< Number of successfully loaded t-matrices
|
||||
double* *const freqs, ///< Frequencies in SI units
|
||||
double* *const freqs_su, ///< Frequencies in SCUFF units (optional)
|
||||
qpms_tmatrix_t* *const tmatrices_array, ///< The resulting T-matrices (optional).
|
||||
complex double* *const tmdata
|
||||
) {
|
||||
if (!(freqs && n && tmdata))
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"freqs, n, and tmdata are mandatory arguments and must not be NULL.");
|
||||
if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT?
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)"
|
||||
" norm supported right now.");
|
||||
int n_alloc = 128; // First chunk to allocate
|
||||
*n = 0;
|
||||
*freqs = malloc(n_alloc * sizeof(double));
|
||||
if (freqs_su) *freqs_su = malloc(n_alloc * sizeof(double));
|
||||
*tmdata = malloc(sizeof(complex double) * bs->n * bs->n * n_alloc);
|
||||
if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"malloc() failed.");
|
||||
size_t linebufsz = 256;
|
||||
char *linebuf = malloc(linebufsz);
|
||||
ssize_t readchars;
|
||||
double lastfreq_su = NAN;
|
||||
while((readchars = getline(&linebuf, &linebufsz, f)) != -1) {
|
||||
if (linebuf[0] == '#') continue;
|
||||
int Alpha, LAlpha, MAlpha, PAlpha, Beta, LBeta, MBeta, PBeta;
|
||||
double currentfreq_su, tr, ti;
|
||||
if (11 != sscanf(linebuf, "%lf %d %d %d %d %d %d %d %d %lf %lf",
|
||||
¤tfreq_su, &Alpha, &LAlpha, &MAlpha, &PAlpha,
|
||||
&Beta, &LBeta, &MBeta, &PBeta, &tr, &ti))
|
||||
abort(); // Malformed T-matrix file
|
||||
if (currentfreq_su != lastfreq_su) { // New frequency -> new T-matrix
|
||||
++*n;
|
||||
lastfreq_su = currentfreq_su;
|
||||
if(*n > n_alloc) {
|
||||
n_alloc *= 2;
|
||||
*freqs = realloc(*freqs, n_alloc * sizeof(double));
|
||||
if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double));
|
||||
*tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc);
|
||||
if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"realloc() failed.");
|
||||
}
|
||||
if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su;
|
||||
(*freqs)[*n-1] = qpms_SU2SI(currentfreq_su);
|
||||
|
||||
for(size_t i = 0; i < bs->n * bs->n; ++i)
|
||||
(*tmdata)[(*n-1)*bs->n*bs->n + i] = NAN + I*NAN;
|
||||
}
|
||||
qpms_vswf_type_t TAlpha, TBeta;
|
||||
switch(PAlpha) {
|
||||
case 0: TAlpha = QPMS_VSWF_MAGNETIC; break;
|
||||
case 1: TAlpha = QPMS_VSWF_ELECTRIC; break;
|
||||
default: assert(0);
|
||||
}
|
||||
switch(PBeta) {
|
||||
case 0: TBeta = QPMS_VSWF_MAGNETIC; break;
|
||||
case 1: TBeta = QPMS_VSWF_ELECTRIC; break;
|
||||
default: assert(0);
|
||||
}
|
||||
qpms_uvswfi_t srcui = qpms_tmn2uvswfi(TAlpha, MAlpha, LAlpha),
|
||||
destui = qpms_tmn2uvswfi(TBeta, MBeta, LBeta);
|
||||
ssize_t srci = qpms_vswf_set_spec_find_uvswfi(bs, srcui),
|
||||
desti = qpms_vswf_set_spec_find_uvswfi(bs, destui);
|
||||
if (srci == -1 || desti == -1)
|
||||
/* This element has not been requested in bs->ilist. */
|
||||
continue;
|
||||
else
|
||||
(*tmdata)[(*n-1)*bs->n*bs->n + desti*bs->n + srci] = tr + I*ti;
|
||||
}
|
||||
free(linebuf);
|
||||
// free some more memory
|
||||
n_alloc = *n;
|
||||
*freqs = realloc(*freqs, n_alloc * sizeof(double));
|
||||
if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double));
|
||||
if (tmatrices_array) *tmatrices_array = realloc(*tmatrices_array, n_alloc * sizeof(qpms_tmatrix_t));
|
||||
*tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc);
|
||||
if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"realloc() failed.");
|
||||
if (tmatrices_array) {
|
||||
*tmatrices_array = malloc(n_alloc * sizeof(qpms_tmatrix_t));
|
||||
if (!*tmatrices_array)
|
||||
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
|
||||
"malloc() failed.");
|
||||
for (size_t i = 0; i < *n; ++i) {
|
||||
(*tmatrices_array)[i].spec = bs;
|
||||
(*tmatrices_array)[i].m = (*tmdata) + i * bs->n * bs->n;
|
||||
(*tmatrices_array)[i].owns_m = false;
|
||||
}
|
||||
}
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_load_scuff_tmatrix(
|
||||
const char *path, ///< file path
|
||||
const qpms_vswf_set_spec_t * bs, ///< VSWF set spec
|
||||
size_t *const n, ///< Number of successfully loaded t-matrices
|
||||
double **const freqs, ///< Frequencies in SI units
|
||||
double ** const freqs_su, ///< Frequencies in SCUFF units (optional)
|
||||
qpms_tmatrix_t ** const tmatrices_array, ///< The resulting T-matrices (optional).
|
||||
complex double ** const tmdata
|
||||
) {
|
||||
FILE *f = fopen(path, "r");
|
||||
if (!f)
|
||||
qpms_pr_error_at_line(__FILE__, __LINE__, __func__,
|
||||
"Could not open T-matrix file %s", path);
|
||||
qpms_errno_t retval =
|
||||
qpms_read_scuff_tmatrix(f, bs, n, freqs, freqs_su, tmatrices_array, tmdata);
|
||||
if(fclose(f)) qpms_pr_error_at_line(__FILE__, __LINE__, __func__,
|
||||
"Could not close the T-matrix file %s (well, that's weird, "
|
||||
"since it's read only).", path);
|
||||
|
||||
return retval;
|
||||
}
|
||||
|
||||
|
||||
2
setup.py
2
setup.py
|
|
@ -27,7 +27,7 @@ qpms_c = Extension('qpms_c',
|
|||
'qpms/scatsystem.c',
|
||||
'qpms/vswf.c', # FIXME many things from vswf.c are not required by this module, but they have many dependencies (following in this list); maybe I want to move all the "basespec stuff"
|
||||
'qpms/legendre.c',
|
||||
'qpms/tmatrix_io.c',
|
||||
'qpms/tmatrices.c',
|
||||
'qpms/error.c'],
|
||||
extra_compile_args=['-std=c99','-ggdb', '-O3',
|
||||
'-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required
|
||||
|
|
|
|||
Loading…
Reference in New Issue