oldtests/ss_syms_packs.c

   1 //  c99 -g -I.. ss_syms_packs.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c  ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c  ../qpms/wigner.c -lm -lgsl -lblas -llapacke
   2 typedef int qpms_gmi_t;// There is something wrong in the includes, apparently.
   3 #include <qpms/qpms_types.h>
   4 #include <qpms/scatsystem.h>
   5 #include <stdlib.h>
   6 #include <qpms/vswf.h>
   7 #include <qpms/indexing.h>
   8 #include <stdio.h>
   9 #include "staticgroups.h"
  10
  11 const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h;
  12 const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v;
  13 const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g;
  14 const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v;
  15 const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h;
  16 const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h;
  17
  18 double uniform_random(double min, double max) {
  19   double random_value = min + (max-min)*(double)rand()/RAND_MAX;
  20   return random_value;
  21 }
  22
  23 int main()
  24 {
  25   srand(666);
  26 #if 0
  27   qpms_vswf_set_spec_t
  28     *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS),
  29     *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS);
  30 #else
  31   // Only electric waves
  32   qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(),
  33                        *b2 = qpms_vswf_set_spec_init();
  34   b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS;
  35   for(qpms_l_t l = 1; l <= 1; ++l)
  36     for (qpms_m_t m = -l; m <= l; ++m)
  37       qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
  38   for(qpms_l_t l = 1; l <= 1; ++l)
  39     for (qpms_m_t m = -l; m <= l; ++m)
  40       qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
  41 #endif
  42   qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1);
  43   qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2);
  44
  45 #if 0
  46   // Random diagonal T-matrices
  47   for(size_t i = 0; i < b1->n; ++i)
  48     t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
  49   for(size_t i = 0; i < b2->n; ++i)
  50     t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
  51 #else
  52   for(size_t i = 0; i < b1->n; ++i)
  53     t1->m[i + i*b1->n] = 1;
  54   for(size_t i = 0; i < b2->n; ++i)
  55     t2->m[i + i*b2->n] = 1;
  56 #endif
  57
  58   const cart3_t pp1 = {0, 0, 1}, pp2 = {0,0, 2}, pp3 = {0,0 , 0};
  59   qpms_tmatrix_t * tmlist[] = {t1, t2};
  60   qpms_particle_tid_t plist[] = {{pp1, 0}, {pp2, 0}, {pp3, 1}};
  61
  62   qpms_scatsys_t protoss;
  63   protoss.tm = tmlist;
  64   protoss.tm_count=2;
  65   protoss.p = plist;
  66   protoss.p_count=3;
  67
  68   qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D3h);
  69
  70   printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n",
  71       (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps,
  72       (int)ss->orbit_type_count);
  73
  74   const double k = 1.7;
  75
  76   complex double *S_full = qpms_scatsys_build_translation_matrix_full(
  77       NULL, ss, k);
  78   complex double *S_packed[ss->sym->nirreps];
  79   for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri)
  80     S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL,
  81         S_full, ss, iri);
  82
  83   complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL,
  84       S_packed[0], ss, 0, false);
  85   for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri)
  86     qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri],
  87         ss, iri, true);
  88
  89   double maxerr = 0;
  90   for (size_t i = 0; i < ss->fecv_size; ++i) {
  91     double err = cabs(S_full[i] - S_recfull[i]);
  92     maxerr = (err > maxerr) ? err : maxerr;
  93   }
  94
  95   printf("maxerr: %lg\n", maxerr);
  96
  97   for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]);
  98   free(S_full);
  99   qpms_scatsys_free(ss);
 100   qpms_tmatrix_free(t1);
 101   qpms_tmatrix_free(t2);
 102   qpms_vswf_set_spec_free(b1);
 103   qpms_vswf_set_spec_free(b2);
 104   return 0;
 105 }