/
crosssec.c
1352 lines (1239 loc) · 55.9 KB
/
crosssec.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* All the functions to calculate scattering quantities (except Mueller matrix), to read different parameters from
* files, and to initialize orientation of the particle
*
* Copyright (C) ADDA contributors
* This file is part of ADDA.
*
* ADDA is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
*
* ADDA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along with ADDA. If not, see
* <http://www.gnu.org/licenses/>.
*/
#include "const.h" // keep this first
#include "crosssec.h" // corresponding header
// project headers
#include "cmplx.h"
#include "comm.h"
#include "debug.h"
#include "io.h"
#include "memory.h"
#include "Romberg.h"
#include "timing.h"
#include "vars.h"
// system headers
#include <limits.h>
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
// SEMI-GLOBAL VARIABLES
// defined and initialized in calculator.c
extern doublecomplex * restrict E_ad;
extern double * restrict E2_alldir;
extern const doublecomplex cc[][3];
#ifndef SPARSE
extern doublecomplex * restrict expsX,* restrict expsY,* restrict expsZ;
#endif
// defined and initialized in param.c
extern const double incPolX_0[3],incPolY_0[3];
extern const enum scat ScatRelation;
// defined and initialized in timing.c
extern TIME_TYPE Timing_EFieldAD,Timing_EFieldADComm,Timing_EFieldSG,Timing_EFieldSGComm,
Timing_ScatQuanComm;
// used in CalculateE.c
Parms_1D phi_sg;
double ezLab[3]; // basis vector ez of laboratory RF transformed into the RF of particle
double exSP[3]; // second vector, which determines the scattering plane passing through prop and ez (in particle RF)
// used in calculator.c
Parms_1D parms_alpha; // parameters of integration over alpha
Parms_1D parms[2]; // parameters for integration over theta,phi or beta,gamma
angle_set beta_int,gamma_int,theta_int,phi_int; // sets of angles
// used in param.c
const char *avg_string; // string for output of function that reads averaging parameters
// used in Romberg.c
bool full_al_range; // whether full range of alpha angle is used
// LOCAL VARIABLES
static double exLab[3],eyLab[3]; // basis vectors of laboratory RF transformed into the RF of particle
//======================================================================================================================
static inline int AlldirIndex(const int theta,const int phi)
// Convert the (theta,phi) couple into a linear array index
{
return (theta*phi_int.N + phi);
}
//======================================================================================================================
void InitRotation (void)
/* initialize matrices used for reference frame transformation; based on Mishchenko M.I. "Calculation of the amplitude
* matrix for a nonspherical particle in a fixed orientation", Applied Optics 39(6):1026-1031. This is so-called
* zyz-notation or y-convention. Instead of rotating the particle (active rotation) we rotate everything else (passive).
* Then the beta_matr is the transpose of the active one.
*/
{
double ca,sa,cb,sb,cg,sg;
double beta_matr[3][3];
double alph,bet,gam; // in radians
double ex[3]; // second vector of scattering plane in the laboratory RF
// initialization of angle values in radians
alph=Deg2Rad(alph_deg);
bet=Deg2Rad(bet_deg);
gam=Deg2Rad(gam_deg);
// calculation of rotation matrix
ca=cos(alph);
sa=sin(alph);
cb=cos(bet);
sb=sin(bet);
cg=cos(gam);
sg=sin(gam);
beta_matr[0][0]=ca*cb*cg-sa*sg;
beta_matr[0][1]=sa*cb*cg+ca*sg;
beta_matr[0][2]=-sb*cg;
beta_matr[1][0]=-ca*cb*sg-sa*cg;
beta_matr[1][1]=-sa*cb*sg+ca*cg;
beta_matr[1][2]=sb*sg;
beta_matr[2][0]=ca*sb;
beta_matr[2][1]=sa*sb;
beta_matr[2][2]=cb;
// rotation of incident field and basis vectors of laboratory RF
MatrVec(beta_matr,prop_0,prop);
MatrVec(beta_matr,incPolY_0,incPolY);
MatrVec(beta_matr,incPolX_0,incPolX);
MatrColumn(beta_matr,0,exLab);
MatrColumn(beta_matr,1,eyLab);
MatrColumn(beta_matr,2,ezLab);
/* define the second vector (x) of the scattering plane (through prop and ez)
* if prop is along ez, we choose the original ex (that is - we generally assume that beam propagates from the left
* side and positive x-direction is in the right side (left and right is with respect to the possible rotation of
* ex,ey over the z-axis
*/
if (prop_0[2]==1) vCopy(incPolX_0,ex);
else if (prop_0[2]==-1) vMultScal(-1,incPolX_0,ex);
else {
ex[0]=prop_0[0];
ex[1]=prop_0[1];
ex[2]=0;
vNormalize(ex);
}
MatrVec(beta_matr,ex,exSP);
// if needed rotate beam center
if (beam_asym) MatrVec(beta_matr,beam_center_0,beam_center);
}
//======================================================================================================================
static void ReadLineStart(FILE * restrict file, // opened file
const char * restrict fname, // ... its filename
char * restrict buf,const int buf_size, // buffer for line and its size
const char * restrict start) // beginning of the line to search
// reads the first line that starts with 'start'
{
while (!feof(file)) {
fgets(buf,buf_size,file);
if (strstr(buf,start)==buf) { // if correct beginning
if (strstr(buf,"\n")==NULL && !feof(file))
LogError(ONE_POS,"Buffer overflow while reading '%s' (size of essential line > %d)",fname,buf_size-1);
else return; // line found and fits into buffer
} // finish reading unmatched line
else while (strstr(buf,"\n")==NULL && !feof(file)) fgets(buf,buf_size,file);
}
LogError(ONE_POS,"String '%s' is not found (in correct place) in file '%s'",start,fname);
}
//======================================================================================================================
static inline void ScanDouble(FILE * restrict file,const char * restrict fname,char * restrict buf,const int buf_size,
const char * restrict start,double *res)
/* scans double value from a line starting with exactly 'start'; contains the same arguments as ReadLineStart function,
* plus pointer to where the result should be placed
*/
{
ReadLineStart(file,fname,buf,buf_size,start);
if (sscanf(buf+strlen(start),"%lf",res)!=1)
LogError(ONE_POS,"Error reading value after '%s' in file '%s'",start,fname);
}
//======================================================================================================================
static inline void ScanInt(FILE * restrict file,const char * restrict fname,char * restrict buf,const int buf_size,
const char * restrict start,int *res)
/* scans integer value from a line starting with exactly 'start'; contains the same arguments as ReadLineStart function,
* plus pointer to where the result should be placed
*/
{
double tmp;
ReadLineStart(file,fname,buf,buf_size,start);
if (sscanf(buf+strlen(start),"%lf",&tmp)!=1)
LogError(ONE_POS,"Error reading value after '%s' in file '%s'",start,fname);
if (tmp<INT_MIN || tmp>INT_MAX)
LogError(ONE_POS,"Value after '%s' in file '%s' is out of integer bounds",start,fname);
if (sscanf(buf+strlen(start),"%d",res)!=1)
LogError(ONE_POS,"Error reading value after '%s' in file '%s'",start,fname);
}
//======================================================================================================================
static inline void ScanSizet(FILE * restrict file,const char * restrict fname,char * restrict buf,const int buf_size,
const char * restrict start,size_t *res)
/* scans large integer value from a line starting with exactly 'start'; contains the same arguments as ReadLineStart
* function, plus pointer to where the result should be placed. MinGW already provides C99-compliant printf-style
* functions, but not yet scanf-style. So we have to use workarounds instead of straightforward "%zu" format specifier.
* TODO: change to %zu, when libmingwex will include scanf.
*/
{
double tmp;
unsigned long res_tmp;
ReadLineStart(file,fname,buf,buf_size,start);
if (sscanf(buf+strlen(start),"%lf",&tmp)!=1)
LogError(ONE_POS,"Error reading value after '%s' in file '%s'",start,fname);
if (tmp<0 || tmp>SIZE_MAX) LogError(ONE_POS,"Value after '%s' in file '%s' is out of size_t bounds",start,fname);
if (sscanf(buf+strlen(start),"%lu",&res_tmp)!=1)
LogError(ONE_POS,"Error reading value after '%s' in file '%s'",start,fname);
*res=(size_t)res_tmp;
}
//======================================================================================================================
static inline void ScanString(FILE * restrict file,const char * restrict fname,char * restrict buf,const int buf_size,
const char * restrict start,char * restrict res)
/* scans string value from a line starting with exactly 'start'; contains the same arguments as ReadLineStart function,
* plus pointer to where the result should be placed; the memory allocated to 'res' should be at least buf_size and
* independent of buf.
*/
{
ReadLineStart(file,fname,buf,buf_size,start);
if (sscanf(buf+strlen(start),"%s",res)!=1)
LogError(ONE_POS,"Error reading value after '%s' in file '%s'",start,fname);
/* More secure would be to put field width in format string above (like "%.Ns"), however this field width is
* defined by the variable buf_size. The latter can only be implemented by a preliminary printf to get a format
* string.
*/
}
//======================================================================================================================
static void ScanIntegrParms(
FILE * restrict file,const char * restrict fname, // opened file and filename
angle_set *a, // pointer to angle set
Parms_1D *b, // pointer to parameters of integration
const bool ifcos, // if space angles equally in cos
char * restrict buf,char* restrict temp, // 2 independent buffers
const int buf_size) // and their size
// scan integration parameters for angles from file
{
size_t i;
double unit;
// scan file
ScanDouble(file,fname,buf,buf_size,"min=",&(a->min));
ScanDouble(file,fname,buf,buf_size,"max=",&(a->max));
ScanInt(file,fname,buf,buf_size,"Jmin=",&(b->Jmin));
ScanInt(file,fname,buf,buf_size,"Jmax=",&(b->Jmax));
ScanDouble(file,fname,buf,buf_size,"eps=",&(b->eps));
ScanString(file,fname,buf,buf_size,"equiv=",temp);
if (strcmp(temp,"true")==0) b->equival=true;
else if (strcmp(temp,"false")==0) b->equival=false;
else LogError(ONE_POS,"Wrong argument of 'equiv' option in file %s",fname);
ScanString(file,fname,buf,buf_size,"periodic=",temp);
if (strcmp(temp,"true")==0) b->periodic=true;
else if (strcmp(temp,"false")==0) b->periodic=false;
else LogError(ONE_POS,"Wrong argument of 'periodic' option in file %s",fname);
// fill all parameters
if (a->min==a->max) {
a->N=b->Grid_size=1;
b->Jmax=1;
}
else {
// consistency check
if (a->min>a->max) LogError(ONE_POS,
"Wrong range (min="GFORMDEF", max="GFORMDEF") in file %s (max must be >= min)",a->min,a->max,fname);
if (b->Jmax<b->Jmin)
LogError(ONE_POS,"Wrong Jmax (%d) in file %s; it must be >= Jmin (%d)",b->Jmax,fname,b->Jmin);
if (b->Jmin<1) LogError(ONE_POS,"Wrong Jmin (%d) in file %s (must be >=1)",b->Jmin,fname);
if (b->eps<0) LogError(ONE_POS,"Wrong eps ("GFORMDEF") in file %s (must be >=0)",b->eps,fname);
if (b->Jmax >= (int)(8*sizeof(int)))
LogError(ONE_POS,"Too large Jmax(%d) in file %s, it will cause integer overflow",b->Jmax,fname);
a->N=b->Grid_size=(1 << b->Jmax) + 1;
if (b->equival && a->N>1) (a->N)--;
}
// initialize points of integration
MALLOC_VECTOR(a->val,double,a->N,ALL);
memory += a->N*sizeof(double);
if (ifcos) { // make equal intervals in cos(angle)
// consistency check
if (a->min<0) LogError(ONE_POS,"Wrong min ("GFORMDEF") in file %s (must be >=0 for this angle)",a->min,fname);
if (a->max>180)
LogError(ONE_POS,"Wrong max ("GFORMDEF") in file %s (must be <=180 for this angle)",a->max,fname);
b->min=cos(Deg2Rad(a->max));
b->max=cos(Deg2Rad(a->min));
if (fabs(b->min)<ROUND_ERR) b->min=0; // just for convenience of display in log file
if (fabs(b->max)<ROUND_ERR) b->max=0;
if (b->Grid_size==1) a->val[0]=a->min;
else {
unit = (b->max - b->min)/(b->Grid_size-1);
for (i=0;i<a->N;i++) a->val[i] = Rad2Deg(acos(b->min+unit*i));
}
}
else { // make equal intervals in angle
b->min=Deg2Rad(a->min);
b->max=Deg2Rad(a->max);
if (b->Grid_size==1) a->val[0]=a->min;
else {
unit = (a->max - a->min)/(b->Grid_size-1);
for (i=0;i<a->N;i++) a->val[i] = a->min + unit*i;
}
}
}
//======================================================================================================================
static enum angleset ScanAngleSet(
FILE * restrict file,const char * restrict fname, // opened file and filename
angle_set *a, // pointers to angle set
char * restrict buf,char * restrict temp, // 2 independent buffers
const int buf_size) // and their size
// scan range or set of angles (theta or phi) from file (used for scat_grid)
{
size_t i;
double unit;
enum angleset out;
ScanString(file,fname,buf,buf_size,"type=",temp);
ScanSizet(file,fname,buf,buf_size,"N=",&(a->N));
// initialize angle array
MALLOC_VECTOR(a->val,double,a->N,ALL);
memory += a->N*sizeof(double);
if (strcmp(temp,"range")==0) {
ScanDouble(file,fname,buf,buf_size,"min=",&(a->min));
ScanDouble(file,fname,buf,buf_size,"max=",&(a->max));
if (a->min>a->max) LogError(ONE_POS,
"Wrong range (min="GFORMDEF", max="GFORMDEF") in file %s (max must be >= min)",a->min,a->max,fname);
if (a->N==1) a->val[0]=(a->max + a->min)/2;
else {
unit = (a->max - a->min)/(a->N - 1);
for (i=0;i<a->N;i++) a->val[i] = a->min + unit*i;
}
out=AS_RANGE;
}
else if (strcmp(temp,"values")==0) {
ReadLineStart(file,fname,buf,buf_size,"values=");
for (i=0;i<a->N;i++) {
fgets(buf,buf_size,file);
if (strstr(buf,"\n")==NULL && !feof(file))
LogError(ONE_POS,"Buffer overflow while scanning lines in file '%s' (line size > %d)",fname,buf_size-1);
if (sscanf(buf,"%lf\n",a->val+i)!=1)
LogError(ONE_POS,"Failed scanning values from line '%s' in file '%s'",buf,fname);
}
out=AS_VALUES;
}
else LogError(ONE_POS,"Unknown type '%s' in file '%s'",temp,fname);
return out;
}
//======================================================================================================================
void ReadAvgParms(const char * restrict fname)
// read parameters of orientation averaging from a file
{
FILE * restrict input;
char buf[BUF_LINE],temp[BUF_LINE];
TIME_TYPE tstart=GET_TIME();
// open file
input=FOpenErr(fname,"r",ALL_POS);
//scan file
ReadLineStart(input,fname,buf,BUF_LINE,"alpha:");
ScanIntegrParms(input,fname,&alpha_int,&parms_alpha,false,buf,temp,BUF_LINE);
full_al_range=fabs(alpha_int.max-alpha_int.min-FULL_ANGLE)<FULL_ANGLE*ROUND_ERR;
ReadLineStart(input,fname,buf,BUF_LINE,"beta:");
ScanIntegrParms(input,fname,&beta_int,&parms[THETA],true,buf,temp,BUF_LINE);
ReadLineStart(input,fname,buf,BUF_LINE,"gamma:");
ScanIntegrParms(input,fname,&gamma_int,&parms[PHI],false,buf,temp,BUF_LINE);
// close file
FCloseErr(input,fname,ALL_POS);
// print info to string
if (IFROOT) avg_string=dyn_sprintf(
"alpha: from "GFORMDEF" to "GFORMDEF" in %zu steps\n"
"beta: from "GFORMDEF" to "GFORMDEF" in (up to) %zu steps (equally spaced in cosine "
"values)\n"
"gamma: from "GFORMDEF" to "GFORMDEF" in (up to) %zu steps\n"
"see file 'log_orient_avg' for details\n",
alpha_int.min,alpha_int.max,alpha_int.N,beta_int.min,beta_int.max,beta_int.N,gamma_int.min,
gamma_int.max,gamma_int.N);
D("ReadAvgParms finished");
Timing_FileIO+=GET_TIME()-tstart;
}
//======================================================================================================================
void ReadAlldirParms(const char * restrict fname)
/* read integration parameters for asymmetry-parameter & C_sca; should not be used together with orientation averaging
* because they use the same storage space - parms
*/
{
FILE * restrict input;
char buf[BUF_LINE],temp[BUF_LINE];
TIME_TYPE tstart=GET_TIME();
// open file
input=FOpenErr(fname,"r",ALL_POS);
//scan file
ReadLineStart(input,fname,buf,BUF_LINE,"theta:");
ScanIntegrParms(input,fname,&theta_int,&parms[THETA],true,buf,temp,BUF_LINE);
ReadLineStart(input,fname,buf,BUF_LINE,"phi:");
ScanIntegrParms(input,fname,&phi_int,&parms[PHI],false,buf,temp,BUF_LINE);
// close file
FCloseErr(input,fname,ALL_POS);
// print info
if (IFROOT) fprintf(logfile,"\n"
"Scattered field is calculated for all directions (for integrated scattering quantities)\n"
"theta: from "GFORMDEF" to "GFORMDEF" in (up to) %zu steps (equally spaced in cosine values)\n"
"phi: from "GFORMDEF" to "GFORMDEF" in (up to) %zu steps\n"
"see files 'log_int_***' for details\n\n",
theta_int.min,theta_int.max,theta_int.N,phi_int.min,phi_int.max,phi_int.N);
D("ReadAlldirParms finished");
Timing_FileIO+=GET_TIME()-tstart;
}
//======================================================================================================================
void ReadScatGridParms(const char * restrict fname)
// read parameters of the grid on which to calculate scattered field
{
FILE * restrict input;
char buf[BUF_LINE],temp[BUF_LINE];
enum angleset theta_type,phi_type;
size_t i;
TIME_TYPE tstart=GET_TIME();
// redundant initialization to remove warnings
theta_type=phi_type=AS_RANGE;
// open file
input=FOpenErr(fname,"r",ALL_POS);
// scan file
ScanString(input,fname,buf,BUF_LINE,"global_type=",temp);
if (strcmp(temp,"grid")==0) {
angles.type = SG_GRID;
ReadLineStart(input,fname,buf,BUF_LINE,"theta:");
theta_type=ScanAngleSet(input,fname,&(angles.theta),buf,temp,BUF_LINE);
if (phi_integr) {
ReadLineStart(input,fname,buf,BUF_LINE,"phi_integr:");
ScanIntegrParms(input,fname,&(angles.phi),&phi_sg,false,buf,temp,BUF_LINE);
phi_type = AS_RANGE;
}
else {
ReadLineStart(input,fname,buf,BUF_LINE,"phi:");
phi_type=ScanAngleSet(input,fname,&(angles.phi),buf,temp,BUF_LINE);
}
angles.N=MultOverflow(angles.theta.N,angles.phi.N,ONE_POS,"angles.N");;
}
else if (strcmp(temp,"pairs")==0) {
if (phi_integr) LogError(ONE_POS,"Integration over phi can't be done with 'global_type=pairs'");
angles.type = SG_PAIRS;
ScanSizet(input,fname,buf,BUF_LINE,"N=",&(angles.N));
angles.theta.N=angles.phi.N=angles.N;
// malloc angle arrays
MALLOC_VECTOR(angles.theta.val,double,angles.N,ALL);
MALLOC_VECTOR(angles.phi.val,double,angles.N,ALL);
memory += 2*angles.N*sizeof(double);
ReadLineStart(input,fname,buf,BUF_LINE,"pairs=");
for (i=0;i<angles.N;i++) {
fgets(buf,BUF_LINE,input);
if (strstr(buf,"\n")==NULL && !feof(input))
LogError(ONE_POS,"Buffer overflow while scanning lines in file '%s' (line size > %d)",fname,BUF_LINE-1);
if (sscanf(buf,"%lf %lf\n",angles.theta.val+i,angles.phi.val+i)!=2)
LogError(ONE_POS,"Failed scanning values from line '%s' in file '%s'",buf,fname);
}
}
else LogError(ONE_POS,"Unknown global_type '%s' in file '%s'",temp,fname);
// close file
FCloseErr(input,fname,ALL_POS);
// print info
if (IFROOT) {
fprintf(logfile,"\nScattered field is calculated for multiple directions\n");
switch (angles.type) {
case SG_GRID:
switch (theta_type) {
case AS_RANGE:
fprintf(logfile,"theta: from "GFORMDEF" to "GFORMDEF" in %zu steps\n",angles.theta.min,
angles.theta.max,angles.theta.N);
break;
case AS_VALUES: fprintf(logfile,"theta: %zu given values\n",angles.theta.N); break;
}
switch (phi_type) {
case AS_RANGE:
fprintf(logfile,"phi: from "GFORMDEF" to "GFORMDEF" in %zu steps\n",angles.phi.min,
angles.phi.max,angles.phi.N);
if (phi_integr) fprintf(logfile,"(Mueller matrix is integrated over phi)\n");
break;
case AS_VALUES: fprintf(logfile,"phi: %zu given values\n",angles.phi.N); break;
}
break;
case SG_PAIRS: fprintf(logfile,"Total %zu given (theta,phi) pairs\n",angles.N); break;
}
fprintf(logfile,"\n");
}
D("ReadScatGridParms finished");
Timing_FileIO+=GET_TIME()-tstart;
}
//======================================================================================================================
static inline double eta2(const double n[static restrict 3])
/* calculates IGT_SO correction for scattering at direction n. Exact formula is based on integration of exp(ikn.r) over
* the dipole volume, resulting in Product(sinc(kd[mu]*n[mu]/2),mu). But here we use a second-order approximation.
* Does not depend on n for cubical dipoles.
*/
{
return 1-(kdX*kdX*n[0]*n[0]+kdY*kdY*n[1]*n[1]+kdZ*kdZ*n[2]*n[2])/24;
}
//======================================================================================================================
static inline doublecomplex eta2cmplx(const doublecomplex n[static restrict 3])
// same as eta2, but for complex input vector. Does not depend on n (and is real) for cubical dipoles if n.n=1
{
return 1-(kdX*kdX*n[0]*n[0]+kdY*kdY*n[1]*n[1]+kdZ*kdZ*n[2]*n[2])/24;
}
//======================================================================================================================
static void CalcFieldFree(doublecomplex ebuff[static restrict 3], // where to write calculated scattering amplitude
const double n[static restrict 3]) // scattering direction
/* Near-optimal routine to compute the scattered fields at one specific angle (more exactly - scattering amplitude);
* Specific optimization are possible when e.g. n[0]=0 for scattering in yz-plane, however in this case it is very
* improbable that the routine will become a bottleneck. The latter happens mostly for cases, when grid of scattering
* angles is used with only small fraction of n, allowing simplifications.
*/
{
double kkk;
doublecomplex a,dpr;
doublecomplex sum[3],tbuff[3],tmp=0; // redundant initialization to remove warnings
int i;
unsigned short ix,iy1,iy2,iz1,iz2;
size_t j,jjj;
#ifdef SPARSE
doublecomplex expX, expY, expZ;
#endif
cvInit(sum);
#ifndef SPARSE
// prepare values of exponents, along each of the coordinates
imExp_arr(-kdX*n[0],boxX,expsX);
imExp_arr(-kdY*n[1],boxY,expsY);
imExp_arr(-kdZ*n[2],local_Nz_unif,expsZ);
#endif // !SPARSE
/* this piece of code tries to use that usually only x position changes from dipole to dipole, saving a complex
* multiplication seems to be beneficial, even considering bookkeeping overhead; it may not be as good for very
* porous particles though, but for them this part of code is anyway fast relative to the FFT on a large grid;
* Further optimization is possible using some kind of plans, i.e. by preliminary analyzing the position of the
* real dipoles on the grid.
*/
iy1=iz1=UNDEF;
for (j=0;j<local_nvoid_Ndip;++j) {
jjj=3*j;
// a=exp(-ikr.n), but r is taken relative to the first dipole of the local box
ix=position[jjj];
iy2=position[jjj+1];
iz2=position[jjj+2];
// the second part is very improbable, but needed for robustness
if (iy2!=iy1 || iz2!=iz1) {
iy1=iy2;
iz1=iz2;
#ifndef SPARSE // FFT mode
tmp=expsY[iy2]*expsZ[iz2];
}
a=tmp*expsX[ix];
#else // sparse mode - the difference is that exponents are not precomputed
expY=imExp(-kdY*n[1]*iy2);
expZ=imExp(-kdZ*n[2]*iz2);
tmp=expY*expZ;
}
expX=imExp(-kdX*n[0]*ix);
a=tmp*expX;
#endif // SPARSE
// sum(P*exp(-ik*r.n))
for(i=0;i<3;i++) sum[i]+=pvec[jjj+i]*a;
} /* end for j */
// tbuff=(I-nxn).sum=sum-n*(n.sum)
dpr=crDotProd(sum,n);
cvMultScal_RVec(dpr,n,tbuff);
cvSubtr(sum,tbuff,tbuff);
// ebuff=(-i*k^3)*exp(-ikr0.n)*tbuff, where r0=box_origin_unif
a=imExp(-WaveNum*DotProd(box_origin_unif,n)); // a=exp(-ikr0.n)
kkk=WaveNum*WaveNum*WaveNum;
// the following additional multiplier implements IGT_SO
if (ScatRelation==SQ_IGT_SO) kkk*=eta2(n);
tmp=-I*a*kkk; // tmp=(-i*k^3)*exp(-ikr0.n)
cvMultScal_cmplx(tmp,tbuff,ebuff);
}
//======================================================================================================================
static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to write calculated scattering amplitude
const double nF[static restrict 3]) // scattering direction (at infinity)
/* Same as CalcFieldFree but for particle near surface.
* For scattering into the substrate we employ the reciprocity principle. The scattered field is obtained from field of
* the plane wave incoming from the scattered direction at the dipole position. In particular,
* E_sca(s,p) = eF(s,p)*(k_0^2/r)*exp(ikr)*t'(s,p) * Sum[P_j.eN_(s,P)*exp(-i*k_0*nN.r_j)],
* where eF,eN are unit [e.e=1] vectors at far and near-field, nN is the normalized transmitted k-vector (also nN.nN=1).
* t' is transmittance coefficient from substrate into the vacuum, k is wavevector in the substrate.
* Total scattered field is obtained by summing s and p components. The actual computed quantity is scattering
* amplitude F, defined as E_sca = F*exp(ikr)/(-ikr).
* Reciprocity should be valid even for absorbing substrate (with any symmetric tensor), not depending on the symmetry
* of the refractive index of the particle itself. In principle, the same formula can be obtained by transmitting
* cylindrical waves emitted by dipole, but that needs additional coefficient due to stretching of wavefront during
* transmission (similar to the difference between amplitude and intensity transmission coefficients).
* Moreover, the consideration of specific scattering angle (and wavenumber) implies that we completely ignore the
* surface plasmon polaritons (SPP), which can, in principle, be considered as scattering at 90 degrees. These SPPs may
* be important for energy balance for metallic substrates. However, for such substrates energy balance is not perfect
* anyway due to absorption.
*/
{
doublecomplex aF,aN,phSh;
doublecomplex cs,cp; // coefficients (reflectance or transmittance) for s- and p-polarizations
doublecomplex ki,kt; // normal component of wavevector above and below the surface
doublecomplex sumF[3],sumN[3],t3[3],tmpF=0,tmpN=0; // redundant initialization to remove warnings
doublecomplex nN[3]; // scattering direction (n.n=1) at near field (corresponds to ktVec in GenerateB.c)
double epF[3],es[3]; // unit vectors of s- and p-polarization, ep differs for near- and far-field
doublecomplex epN[3]; // ep at near-field can be complex
int i;
unsigned short ix,iy1,iy2,iz1,iz2;
size_t j,jjj;
#ifdef SPARSE
doublecomplex expX, expY, expZ;
#endif
const bool above=(nF[2]>-ROUND_ERR); // we assume above-the-surface scattering for all boundary cases (like 90 deg)
cvInit(sumN);
if (above) cvInit(sumF); //additional storage for directly propagated scattering
/* There is an inherent discontinuity for msub approaching 1 and scattering angle 90 degrees (nF[2]=0). The problem
* is that for m=1+-0, but |m-1|>>(nF[2])^2, ki<<kt<<1 => rs=rp=-1
* while for m=1 (exactly) the limit of nF[2]->0 results in kt=ki => rs=rp=0
* Therefore, below is a certain logic, which behaves in an intuitively expected way, for common special cases.
* However, it is still not expected to be continuous for fine-changing parameters (like msub approaching 1).
* In particular, the jump occurs when msub crosses 1+-ROUND_ERR boundary.
* Still, the discontinuity should apply only to scattering at exactly 90 degrees, but not to, e.g., integral
* quantities, like Csca (if sufficient large number of integration points is chosen).
*/
// calculate nN, ki, kt, cs, cp, and phSh
if (above) { // simple reflection
/* No scattering at exactly 90 degrees for non-trivial surface (to avoid randomness for this case).
* See A. Small, J. Fung, and V.N. Manoharan, "Generalization of the optical theorem for light scattering from
* a particle at a planar interface," J. Opt. Soc. Am. A 30, 2519-2525 (2013) for theoretical discussion of
* this fact.
*/
if (fabs(nF[2])<ROUND_ERR && cabs(msub-1)>ROUND_ERR) {
cvInit(ebuff);
return;
}
cvBuildRe(nF,nN);
nN[2]*=-1;
ki=nF[2]; // real
if (msubInf) {
cs=-1;
cp=1;
}
// since kt is not further needed, we directly calculate cs and cp (equivalent to kt=ki)
else if (cabs(msub-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) cs=cp=0;
else { // no special treatment here, since other cases, including 90deg-scattering, are taken care above.
kt=cSqrtCut(msub*msub - (nN[0]*nN[0]+nN[1]*nN[1]));
cs=FresnelRS(ki,kt);
cp=FresnelRP(ki,kt,msub);
}
phSh=imExp(2*WaveNum*hsub*creal(ki)); // assumes real ki
}
else { // transmission; here nF[2] is negative
// formulae correspond to plane wave incoming from below, but with change ki<->kt
if (msubInf) { // no transmission for perfectly reflecting substrate => zero result
cvInit(ebuff);
return;
}
kt=-msub*nF[2];
if (cabs(msub-1)<ROUND_ERR && cabs(kt)<SQRT_RND_ERR) ki=kt;
else ki=cSqrtCut(1 - msub*msub*(nF[0]*nF[0]+nF[1]*nF[1]));
// here nN may be complex, but normalized to n.n=1
nN[0]=msub*nF[0];
nN[1]=msub*nF[1];
nN[2]=-ki;
// these formulae works fine for ki=kt (even very small), and ki=kt=0 is impossible here
cs=FresnelTS(kt,ki);
cp=FresnelTP(kt,ki,1/msub);
// coefficient comes from k0->k in definition of F(n) (in denominator)
phSh=msub*cexp(I*WaveNum*hsub*(ki-kt));
}
#ifndef SPARSE
// prepare values of exponents, along each of the coordinates
imExp_arr(-kdX*nN[0],boxX,expsX);
imExp_arr(-kdY*nN[1],boxY,expsY);
imExp_arr(-kdZ*nN[2],local_Nz_unif,expsZ);
#endif // !SPARSE
/* this piece of code tries to use that usually only x position changes from dipole to dipole, saving a complex
* multiplication seems to be beneficial, even considering bookkeeping overhead; it may not be as good for very
* porous particles though, but for them this part of code is anyway fast relative to the FFT on a large grid;
* Further optimization is possible using some kind of plans, i.e. by preliminary analyzing the position of the
* real dipoles on the grid.
*/
iy1=iz1=UNDEF;
if (above) for (j=0;j<local_nvoid_Ndip;++j) { // two sums need to be calculated
jjj=3*j;
// a=exp(-ikr.n), but r is taken relative to the first dipole of the local box
ix=position[jjj];
iy2=position[jjj+1];
iz2=position[jjj+2];
// the second part is very improbable, but needed for robustness
if (iy2!=iy1 || iz2!=iz1) {
iy1=iy2;
iz1=iz2;
#ifndef SPARSE // FFT mode
tmpN=expsY[iy2]*expsZ[iz2];
tmpF=expsY[iy2]*conj(expsZ[iz2]);
}
aN=tmpN*expsX[ix];
aF=tmpF*expsX[ix];
#else // sparse mode - the difference is that exponents are not precomputed; cexp is used since argument can be complex
expY=cexp(-I*kdY*nN[1]*iy2);
expZ=cexp(-I*kdZ*nN[2]*iz2);
tmpN=expY*expZ;
tmpF=expY*conj(expZ);
}
expX=cexp(-I*kdX*nN[0]*ix);
aN=tmpN*expX;
aF=tmpF*expX;
#endif // SPARSE
// sum(P*exp(-ik*r.nN,F))
for(i=0;i<3;i++) {
sumN[i]+=pvec[jjj+i]*aN;
sumF[i]+=pvec[jjj+i]*aF;
}
} /* end for j above surface */
else for (j=0;j<local_nvoid_Ndip;++j) { // below surface, single sum - similar to free-space scattering
jjj=3*j;
// a=exp(-ikr.n), but r is taken relative to the first dipole of the local box
ix=position[jjj];
iy2=position[jjj+1];
iz2=position[jjj+2];
// the second part is very improbable, but needed for robustness
if (iy2!=iy1 || iz2!=iz1) {
iy1=iy2;
iz1=iz2;
#ifndef SPARSE // FFT mode
tmpN=expsY[iy2]*expsZ[iz2];
}
aN=tmpN*expsX[ix];
#else // sparse mode - the difference is that exponents are not precomputed; cexp is used since argument can be complex
expY=cexp(-I*kdY*nN[1]*iy2);
expZ=cexp(-I*kdZ*nN[2]*iz2);
tmpN=expY*expZ;
}
expX=cexp(-I*kdX*nN[0]*ix);
aN=tmpN*expX;
#endif // SPARSE
// sum(P*exp(-ik*r.nN))
for(i=0;i<3;i++) sumN[i]+=pvec[jjj+i]*aN;
} /* end for j below surface */
// Reflected or transmitted light phSh*(Rs*es(es.sumN) + Rp*epF(epN.sumN)), [dot product w/o conjugation]
/* If reciprocal configuration is rigorously considered signs of vectors nF and nN should be changed, along with
* either es or ep. However, such sign change would not change the final result.
*/
// set unit vectors for s- and p-polarizations; es is the same for all cases
if (vAlongZ(nF)) { // special case - es=ey
es[0]=0;
es[1]=1;
es[2]=0;
}
else { // general case: es = ez x nF /||...||; here we implicitly use that lab RF coincides with particle RF
CrossProd(ezLab,nF,es);
vNormalize(es);
}
CrossProd(es,nF,epF); // epF = es x nF
crCrossProd(es,nN,epN); // epN = es x nN (complex)
// finalize fields
cvMultScal_RVec(phSh*cs*crDotProd(sumN,es),es,ebuff);
cvMultScal_RVec(phSh*cp*cDotProd_conj(sumN,epN),epF,t3);
cvAdd(t3,ebuff,ebuff);
// add directly scattered light, when above the surface (phase shift due to main direction being nN)
if (above) { // ebuff+= [(I-nxn).sum=sum-nF*(nF.sum)] * exp(-2ik*r0*nz), where r0=box_origin_unif
cvMultScal_RVec(crDotProd(sumF,nF),nF,t3);
cvSubtr(sumF,t3,t3);
cvMultScal_cmplx(imExp(-2*WaveNum*creal(ki)*box_origin_unif[2]),t3,t3); // assumes real ki
cvAdd(t3,ebuff,ebuff);
}
// ebuff=(-i*k^3)*exp(-ikr0.n)*tbuff, where r0=box_origin_unif
// All m-scaling for substrate has been accounted in phSh above
doublecomplex sc=-I*WaveNum*WaveNum*WaveNum*cexp(-I*WaveNum*crDotProd(nN,box_origin_unif));
// the following additional multiplier implements IGT_SO; when above, it is the same for nF and nN
if (ScatRelation==SQ_IGT_SO) sc*=eta2cmplx(nN);
cvMultScal_cmplx(sc,ebuff,ebuff);
}
//======================================================================================================================*/
void CalcField(doublecomplex ebuff[static restrict 3], // where to write calculated scattering amplitude
const double n[static restrict 3]) // scattering direction
// wrapper, which redirects the calculation of the field to one of two functions
{
if (surface) CalcFieldSurf(ebuff,n);
else CalcFieldFree(ebuff,n);
}
//======================================================================================================================
double ExtCross(const double * restrict incPol)
// Calculate the Extinction cross-section
{
doublecomplex ebuff[3];
double sum;
size_t i;
// this can be considered a legacy case, which works only for the simplest plane way centered at the particle
if (beamtype==B_PLANE && !surface && !beam_asym) {
CalcField (ebuff,prop);
sum=crDotProd_Re(ebuff,incPol); // incPol is real, so no conjugate is needed
MyInnerProduct(&sum,double_type,1,&Timing_ScatQuanComm);
sum*=FOUR_PI/(WaveNum*WaveNum);
}
/* more general formula; normalization is done assuming the unity amplitude of the electric field in the focal point
* of the beam
*/
else {
sum=0;
for (i=0;i<local_nvoid_Ndip;++i) sum+=cDotProd_Im(pvec+3*i,Einc+3*i); // sum{Im(P.E_inc*)}
MyInnerProduct(&sum,double_type,1,&Timing_ScatQuanComm);
sum*=FOUR_PI*WaveNum;
/* For cubical dipoles the following satisfies IGT_SO, because this factor is applied in CalcField() and is
* independent of propagation or scattering direction. For rectangular dipoles, it is only approximate but
* expected to be accurate for not very elongated dipoles and/or not very spread out incident field.
*
* In principle, the situation is similar for full IGT, but there the correction factor depends on the
* propagation direction even for cubical dipoles
*/
if (ScatRelation==SQ_IGT_SO) sum*=eta2(prop);
}
/* TO ADD NEW BEAM
* The formulae above works only if the amplitude of the beam is unity at the focal point. Either make sure that new
* beam satisfies this condition or add another case here with different formulae.
*/
if (surface) sum*=inc_scale;
return sum;
}
//======================================================================================================================
double AbsCross(void)
// Calculate the Absorption cross-section for process 0
{
size_t dip,index;
int i,j;
unsigned char mat;
double sum,temp1;
double mult[MAX_NMAT][3]; // multiplier (possibly anisotropic)
// Cabs = 4*pi*sum
/* In this function IGT_SO is equivalent to DRAINE. It may seem more logical to make IGT_SO same as FINDIP. However,
* the result is different only for LDR (and similar), for which using IGT does not make a lot of sense anyway.
* Overall, peculiar details related to optical theorem warrant a further study.
*/
switch (ScatRelation) {
/* code below is applicable only for diagonal (for some cases - possibly anisotropic) polarizability and should
* be rewritten otherwise
*/
case SQ_IGT_SO:
case SQ_DRAINE:
/* based on Eq.(35) from Yurkin and Hoekstra, "The discrete dipole approximation: an overview and recent
* developments," JQSRT 106:558-589 (2007).
* summand: Im(P.Eexc(*))-(2/3)k^3*|P|^2=|P|^2*(-Im(1/cc)-(2/3)k^3)
*/
temp1 = 2*WaveNum*WaveNum*WaveNum/3;
for (i=0;i<Nmat;i++) for (j=0;j<3;j++) mult[i][j]=-cimag(1/cc[i][j])-temp1;
for (dip=0,sum=0;dip<local_nvoid_Ndip;++dip) {
mat=material[dip];
index=3*dip;
for(i=0;i<3;i++) sum+=mult[mat][i]*cAbs2(pvec[index+i]);
}
break;
case SQ_FINDIP:
/* based on Eq.(31) or equivalently Eq.(58) from the same paper (ref. above)
* summand: Im(P.E(*))=-|P|^2*Im(chi_inv), chi_inv=1/(V*chi)
*/
temp1 = 2*WaveNum*WaveNum*WaveNum/3;
for (i=0;i<Nmat;i++) for (j=0;j<3;j++) mult[i][j]=-cimag(chi_inv[i][j]);
for (dip=0,sum=0;dip<local_nvoid_Ndip;++dip) {
mat=material[dip];
index=3*dip;
for(i=0;i<3;i++) sum+=mult[mat][i]*cAbs2(pvec[index+i]);
}
break;
}
MyInnerProduct(&sum,double_type,1,&Timing_ScatQuanComm);
if (surface) sum*=inc_scale;
return FOUR_PI*WaveNum*sum;
}
//======================================================================================================================
double DecayCross(void)
// computes total cross section for the dipole incident field; similar to Cext
// 4pi*k*Im[p0(*).Escat(r0)]
{
double sum;
size_t i;
/* This is a correct expression only _if_ exciting p0 is real, then
* (using G(r1,r2) = G(r2,r1)^T, valid also for surface)
* p0(*).Escat_i(r0) = p0(*).G_0i.p_i = p_i.G_i0.p0(*) = p_i.G_i0.p0 = p_i.Einc_i
* => Im(p0(*).Escat(r0)) = sum{Im(P.E_inc)}
*
* For complex p0 an efficient calculation strategy (not to waste evaluations of interaction) is to compute an array
* of G_i0.p0(*) together with Einc and use it here afterwards.
*/
sum=0;
for (i=0;i<local_nvoid_Ndip;++i) sum+=cimag(cDotProd_conj(pvec+3*i,Einc+3*i)); // sum{Im(P.E_inc)}
MyInnerProduct(&sum,double_type,1,&Timing_ScatQuanComm);
return FOUR_PI*WaveNum*sum;
}
//======================================================================================================================
void SetScatPlane(const double ct,const double st,const double phi,double robs[static restrict 3],
double polPer[static restrict 3])
/* Given theta (cos,sin) and phi, calculates scattering direction and unit vector perpendicular to the scattering plane.
* Currently, two alternative definitions are used: theta and phi are either angles in the laboratory reference frame,
* or in the beam reference frame. Generally, polPer = robs x prop (normalized to unit vector), but cases when the two
* latter vectors are aligned requires special treatment.
*
* For special cases with prop||ez we assume that the result is continuous for theta changing from 0 to pi. In other
* words th=0 and pi is effective replaced by th=+0 and pi-0 respectively. Then phi determines the scattering plane.
*
* For other special cases (prop not aligned with ez) the scattering plane is taken to include ez (and hence IncPolX).
*
* Overall, we spend a lot of effort to make consistent definition of polPer. However, its sign is irrelevant because
* change of sign (simultaneously for polPer and both incident and scattering polPar) doesn't change the amplitude
* (and Mueller) matrix.
*/
{
double cphi,sphi;
cphi=cos(phi);
sphi=sin(phi);
if (surface) { // definitions are in the laboratory reference frame
// robs = cos(theta)*ezLab + sin(theta)*[cos(phi)*exLab + sin(phi)*eyLab];
LinComb(exLab,eyLab,cphi,sphi,robs);
LinComb(ezLab,robs,ct,st,robs);
}
else { // definitions are with respect to the incident beam
// robs = cos(theta)*prop + sin(theta)*[cos(phi)*incPolX + sin(phi)*incPolY];
LinComb(incPolX,incPolY,cphi,sphi,robs);
LinComb(prop,robs,ct,st,robs);
}
// set unit vectors which defines Eper
// two special cases; testing ct is more accurate than st==0
// 1) robs has +0 component along exLab (for phi=0): incPolper = (+-)(sin(phi)*exLab - cos(phi)*eyLab)
if (surface && propAlongZ && fabs(st)<ROUND_ERR) LinComb(exLab,eyLab,sphi*propAlongZ,-cphi*propAlongZ,polPer);
// 2) robs has +0 component along incPolX (for phi=0): incPolper = sin(phi)*incPolX - cos(phi)*incPolY
else if (!surface && fabs(st)<ROUND_ERR) LinComb(incPolX,incPolY,sphi,-cphi,polPer);
else { // general case: incPolPer = robserver x prop /||...||
CrossProd(robs,prop,polPer);
// when robs || prop scattering plane contains ez, then polPer = prop x ez/||...|| = -incPolY
// currently this may only occur for surface, since otherwise it is handled by a special case above
if (DotProd(polPer,polPer)==0) vMultScal(-1,incPolY,polPer);
else vNormalize(polPer);
}
}
//======================================================================================================================
void CalcAlldir(void)
// calculate scattered field in many directions
{
int index,npoints,point;
size_t i,j;
TIME_TYPE tstart;
double robserver[3],incPolpar[3],incPolper[3],cthet,sthet,th,ph;
doublecomplex ebuff[3];
// Calculate field
tstart = GET_TIME();
npoints = theta_int.N*phi_int.N;
if (IFROOT) PRINTFB("Calculating scattered field for the whole solid angle:\n");
for (i=0,point=0;i<theta_int.N;++i) {
th=Deg2Rad(theta_int.val[i]);
cthet=cos(th);
sthet=sin(th);
for (j=0;j<phi_int.N;++j) {
ph=Deg2Rad(phi_int.val[j]);
/* set robserver and unit vector for Eper (determines scattering plane); actually Eper and Epar are
* irrelevant, since we are interested only in |E|^2. But projecting the vector on two axes helps to
* somewhat decrease communication time. We also may need these components in the future.
*/
SetScatPlane(cthet,sthet,ph,robserver,incPolper);
// set unit vector for Epar
CrossProd(robserver,incPolper,incPolpar);
// calculate scattered field - main bottleneck
CalcField(ebuff,robserver);
/* Set Epar and Eper - use separate E_ad array to store them (to decrease communications in 1.5 times).
* Writing a special case for sequential mode can eliminate the need of E_ad altogether. Moreover,
* E2_alldir can be stored in 1/4 of memory allocated for E_ad. However, we do not do it, because it doesn't
* seem so significant. And, more importantly, complex fields may also be useful in the future, e.g.
* for radiation force calculation through integration of the far-field
*/
index=2*point;