?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix missing import in finiterectlat-scatter.py
[qpms.git] / qpms / tests / test_planewave_decomposition.c
blobae8f38ec41d02a7c3660f4f128a72fcdd646f86f
1 // c99 -o test_planewave_decomposition -ggdb -I .. test_planewave_decomposition.c ../vswf.c -lgsl -lm -lblas ../legendre.c ../bessel.c
2 #include <gsl/gsl_rng.h>
3 #include <gsl/gsl_math.h>
4 #include <gsl/gsl_randist.h>
5 #include <assert.h>
6 #include "vswf.h"
7 #include "vectors.h"
8 #include "string.h"
9 #include "indexing.h"
10
11 #define _GNU_SOURCE
12 #include <fenv.h>
13
14 char *normstr(qpms_normalisation_t norm) {
15 //int csphase = qpms_normalisation_t_csphase(norm);
16 norm = qpms_normalisation_t_normonly(norm);
17 switch (norm) {
18 case QPMS_NORMALISATION_NONE:
19 return "none";
20 case QPMS_NORMALISATION_SPHARM:
21 return "spharm";
22 case QPMS_NORMALISATION_POWER:
23 return "power";
24 case QPMS_NORMALISATION_XU:
25 return "xu";
26 default:
27 return "!!!undef!!!";
28 }
29 }
30 int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points);
31 int test_planewave_decomposition_silent(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points, double relerrthreshold, double *relerrs);
32
33 int main() {
34 gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxs0);
35 gsl_rng_set(rng, 666);
36 feenableexcept(FE_INVALID | FE_OVERFLOW);
37
38 qpms_l_t lMax = 30;
39 int npoints = 10;
40 double sigma = 1.3;
41 cart3_t points[npoints];
42 double relerrs[npoints];
43 memset(points, 0, npoints * sizeof(cart3_t));
44 points[1].x = points[2].y = points[3].z = 1.;
45 double relerrthreshold = 1e-11;
46 for (unsigned i = 4; i < npoints; ++i) {
47 cart3_t *w = points+i;
48 w->x = gsl_ran_gaussian(rng, sigma);
49 w->y = gsl_ran_gaussian(rng, sigma);
50 w->z = gsl_ran_gaussian(rng, sigma);
51 }
52
53 double ksigma = 15;
54
55 for(int i = 1; i < 2; ++i) {
56 cart3_t k = {0, 0, 2};
57 ccart3_t E = {0., 1.1+I, 0.};
58 if(i){
59 k.x = gsl_ran_gaussian(rng, ksigma);
60 k.y = gsl_ran_gaussian(rng, ksigma);
61 k.z = gsl_ran_gaussian(rng, ksigma);
62 E.x = gsl_ran_gaussian(rng, ksigma);
63 E.x += I* gsl_ran_gaussian(rng, ksigma);
64 E.y = gsl_ran_gaussian(rng, ksigma);
65 E.y += I* gsl_ran_gaussian(rng, ksigma);
66 E.z = gsl_ran_gaussian(rng, ksigma);
67 E.z += I* gsl_ran_gaussian(rng, ksigma);
68 }
69 #if 0
70 printf("Test wave k = %gx̂ + %gŷ + %gẑ", k.x, k.y, k.z);
71 printf(", E_0 = (%g+%gj)x̂ + (%g+%gj)ŷ + (%g+%gj)ẑ\n",
72 creal(E.x),cimag(E.x),
73 creal(E.y),cimag(E.y),
74 creal(E.z),cimag(E.z));
75 printf("%d points, sigma = %g, rel. err. threshold = %g\n", npoints, sigma, relerrthreshold);
76 for(int i = 1; i <= 3; ++i){
77 int res = test_planewave_decomposition_silent(k, E, lMax, i, npoints, points, relerrthreshold, relerrs);
78 printf("\n%s %d/%d %s %s\n", res ? "!!" : "OK", npoints-res, npoints, normstr(i), "");
79 for(int p = 0; p < npoints; ++p) printf("%.3g ", relerrs[p]);
80 res = test_planewave_decomposition_silent(k, E, lMax, i |QPMS_NORMALISATION_T_CSBIT, npoints, points, relerrthreshold, relerrs);
81 printf("\n%s %d/%d %s %s\n", res ? "!!" : "OK", npoints-res, npoints, normstr(i), "with C.-S. phase");
82 for(int p = 0; p < npoints; ++p) printf("%.3g ", relerrs[p]);
83 }
84 #endif
85 for(int i = 1; i <= 4; ++i){
86 test_planewave_decomposition(k, E, lMax, i, npoints, points);
87 test_planewave_decomposition(k, E, lMax, i | QPMS_NORMALISATION_T_CSBIT, npoints, points);
88 }
89
90 }
91 gsl_rng_free(rng);
92 }
93
94 int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points){
95 qpms_y_t nelem = qpms_lMax2nelem(lMax);
96 complex double lc[nelem], mc[nelem], ec[nelem];
97 if (QPMS_SUCCESS !=
98 qpms_planewave2vswf_fill_cart(k, E, lc, mc, ec, lMax, norm)) {
99 printf("Error\n");
100 return -1;
101 }
102 printf("==============================================================\n");
103 printf("Test wave k = %fx̂ + %fŷ + %fẑ", k.x, k.y, k.z);
104 sph_t k_sph = cart2sph(k);
105 printf(" = %fr̂ @ θ = %f, φ = %f\n", k_sph.r, k_sph.theta, k_sph.phi);
106 printf(" E_0 = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
107 creal(E.x),cimag(E.x),
108 creal(E.y),cimag(E.y),
109 creal(E.z),cimag(E.z));
110 csphvec_t E_s = ccart2csphvec(E, k_sph);
111 printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂ @ k\n",
112 creal(E_s.rc), cimag(E_s.rc), creal(E_s.thetac), cimag(E_s.thetac),
113 creal(E_s.phic), cimag(E_s.phic));
114 printf(" lMax = %d, norm: %s, csphase = %d\n",
115 (int)lMax, normstr(norm), qpms_normalisation_t_csphase(norm));
116 printf("a_L: ");
117 for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(lc[y]), cimag(lc[y]));
118 printf("\na_M: ");
119 for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(mc[y]), cimag(mc[y]));
120 printf("\na_N: ");
121 for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(ec[y]), cimag(ec[y]));
122 printf("\n");
123 for (int i = 0; i < npoints; i++) {
124 cart3_t w = points[i];
125 sph_t w_sph = cart2sph(w);
126 printf("Point %d: x = %f, y = %f, z = %f,\n", i, w.x, w.y, w.z);
127 printf(" |r| = %f, θ = %f, φ = %f:\n", w_sph.r, w_sph.theta, w_sph.phi);
128 double phase = cart3_dot(k,w);
129 printf(" k.r = %f\n", phase);
130 complex double phfac = cexp(phase * I);
131 ccart3_t Ew = ccart3_scale(phfac, E);
132 printf(" pw E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
133 creal(Ew.x),cimag(Ew.x),
134 creal(Ew.y),cimag(Ew.y),
135 creal(Ew.z),cimag(Ew.z));
136 csphvec_t Ew_s = ccart2csphvec(Ew, w_sph);
137 printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂ @ r\n",
138 creal(Ew_s.rc), cimag(Ew_s.rc),
139 creal(Ew_s.thetac), cimag(Ew_s.thetac),
140 creal(Ew_s.phic), cimag(Ew_s.phic));
141 w_sph.r *= k_sph.r; /// NEVER FORGET THIS!!!
142 csphvec_t Ew_s_recomp = qpms_eval_vswf(w_sph,
143 lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm);
144 ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph);
145 printf(" rec E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
146 creal(Ew_recomp.x),cimag(Ew_recomp.x),
147 creal(Ew_recomp.y),cimag(Ew_recomp.y),
148 creal(Ew_recomp.z),cimag(Ew_recomp.z));
149 printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂ @ r\n",
150 creal(Ew_s_recomp.rc), cimag(Ew_s_recomp.rc),
151 creal(Ew_s_recomp.thetac), cimag(Ew_s_recomp.thetac),
152 creal(Ew_s_recomp.phic), cimag(Ew_s_recomp.phic));
153 double relerrfac = 2/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc)
154 +cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac)
155 +cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic));
156 printf(" rel. err. magnitude: %g @ r̂, %g @ θ̂, %g @ φ̂\n",
157 cabs(Ew_s_recomp.rc - Ew_s.rc) * relerrfac,
158 cabs(Ew_s_recomp.thetac - Ew_s.thetac) * relerrfac,
159 cabs(Ew_s_recomp.phic - Ew_s.phic) * relerrfac
160 );
161 }
162 return 0;
163 }
164
165 int test_planewave_decomposition_silent(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points, double relerrthreshold, double *relerrs) {
166 qpms_y_t nelem = qpms_lMax2nelem(lMax);
167 int failcount = 0;
168 complex double lc[nelem], mc[nelem], ec[nelem];
169 if (QPMS_SUCCESS !=
170 qpms_planewave2vswf_fill_cart(k, E, lc, mc, ec, lMax, norm)) {
171 printf("Error\n");
172 return -1;
173 }
174 sph_t k_sph = cart2sph(k);
175 csphvec_t E_s = ccart2csphvec(E, k_sph);
176 for (int i = 0; i < npoints; i++) {
177 cart3_t w = points[i];
178 sph_t w_sph = cart2sph(w);
179 w_sph.r *= k_sph.r;
180 double phase = cart3_dot(k,w);
181 complex double phfac = cexp(phase * I);
182 ccart3_t Ew = ccart3_scale(phfac, E);
183 csphvec_t Ew_s = ccart2csphvec(Ew, w_sph);
184 csphvec_t Ew_s_recomp = qpms_eval_vswf(w_sph,
185 lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm);
186 ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph);
187 double relerrfac = 2/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc)
188 +cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac)
189 +cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic));
190
191 double relerr = (cabs(Ew_s_recomp.rc - Ew_s.rc)
192 + cabs(Ew_s_recomp.thetac - Ew_s.thetac)
193 + cabs(Ew_s_recomp.phic - Ew_s.phic)
194 ) * relerrfac;
195 if(relerrs) relerrs[i] = relerr;
196 if(relerr > relerrthreshold) ++failcount;
197 }
198 return failcount;
199 }
QPMS photonic multiple scattering suite