-
Notifications
You must be signed in to change notification settings - Fork 4
/
tdms_iterator.cpp
7437 lines (6557 loc) · 266 KB
/
tdms_iterator.cpp
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
/*****************************************************************
*
* Project.....: isotropic FDTD code
* Application.: main FDTD algorithm
* Module......: iterater.cpp
* Description.: Contains the main FDTD loop as well as other functions
* such as phasor extraction etc. Works in both pulsed
* and steady state mode.
* Compiler....: g++
* Written by..: Peter Munro, Imperial College London, 2002-2008
* Environment.: Linux
* Modified....: Numerous times
*
******************************************************************/
/*iterater_OMP_reduced.cpp iterater_OMP.cpp are the same
*iterater_OMP_full.cpp predates iterater_OMP.cpp and doesn't skip any steps in phasor extraction.
*
*/
/*---------------------------------------------------------------*/
// INCLUDE section
/*---------------------------------------------------------------*/
#include <omp.h>
#include "matio.h"
#include <complex>
#include "tdms_iterator.h"
#include <string.h>
#include "interpolate.h"
#include "numeric.h"
#include "mesh_base.h"
#include "numerical_derivative.h"
#include <time.h>
using namespace std;
#include "globals.h"
#include "matlabio.h"
#include "mesh_base.h"
/*---------------------------------------------------------------*/
// DEFINE section
/*---------------------------------------------------------------*/
//whether of not to time execution
#define TIME_EXEC 0
//time main loop
#define TIME_MAIN_LOOP 1
//threshold used to terminate the steady state iterations
#define TOL 1e-6
//parameter to control the with of the ramp when introducing the waveform in steady state mode
#define ramp_width 4.
//different run modes
#define rm_complete 0
#define rm_analyse 1
//different source modes
#define sm_steadystate 0
#define sm_pulsed 1
//dimensionallity of simulation
#define THREE 0
#define TE 1
#define TM 2
//Used for the 2D case when J_tot=0
inline int max ( int a, int b ) { return a > b ? a : b; }
inline int min ( int a, int b ) { return a < b ? a : b; }
/*This mex function will take in the following arguments
fdtdgrid
Cmaterial
Dmaterial
C
D
freespace
disp_params
delta
interface
Isource
Jsource
Ksource
grid_labels
omega_an
to_l
hwhm
Dxl
Dxu
Dyl
Dyu
Dzl
Dzu
Nt
dt
start_tind
sourcemode
runmode
exphasorsvolume
exphasorssurface
intphasorsurface
phasorsurface
phasorinc
dimension
conductive_aux
dispersive_aux
structure
tdfield
fdtdgrid - A structre with the following members, each of which is a 3 dimensional
array:
fdtdgrid.Exy (double)
fdtdgrid.Exz
fdtdgrid.Eyx
fdtdgrid.Eyz
fdtdgrid.Ezx
fdtdgrid.Ezy
fdtdgrid.Hxy
fdtdgrid.Hxz
fdtdgrid.Hyx
fdtdgrid.Hyz
fdtdgrid.Hzx
fdtdgrid.Hzy
fdtdgrid.materials (uint8)
Cmaterial - A structure with the following members, each of which is a 1 dimensional
array indexed by indices od scattering materials:
Cmaterial.Cax (double)
Cmaterial.Cay
Cmaterial.Caz
Cmaterial.Cbx
Cmaterial.Cby
Cmaterial.Cbz
Cmaterial.Ccx
Cmaterial.Ccy
Cmaterial.Ccz
Dmaterial - A structure with the following members, each of which is a 1 dimensional
array indexed by indices od scattering materials:
Dmaterial.Cax (double)
Dmaterial.Cay
Dmaterial.Caz
Dmaterial.Cbx
Dmaterial.Cby
Dmaterial.Cbz
C - A structure with the same elements as Cmaterial whose elements relate to the background
region. Elements are not restricted to 1 dimension.
D - A structure with the same elements as Cmaterial whose elements relate to the background
region. Elements are not restricted to 1 dimension.
freespace - A structure with the following members each of which is a double
freespace.Cbx (double)
freespace.Cby
freespace.Cbz
freespace.Dbx
freespace.Dby
freespace.Dbz
disp_params - A structure with the following members:
disp_params.alpha (double)
disp_params.beta
disp_params.gamma
delta - A structure with the following elements representing Yee cell dimension
delta.x (double)
delta.y
delta.z
interface - A structure with the following members each of which is a double
interface.I0 - [I0 apply]
interface.I1 - [I1 apply]
interface.J0 - [J0 apply]
interface.J1 - [J1 apply]
interface.K0 - [K0 apply]
interface.K1 - [K1 apply]
where, for example, I0 is the position of the interface plane between
total and scattered formulations. apply is a boolean which indicates
whether or not to apply the boundary condition
Isource - A 3d matrix of dimension 8x(J1-J0+1)x(K1-K0+1) which contains
complex amplitude information for the interface equations at the I0 and I1
planes. Should be Isource(1,:,:) - Ey, Isource(2,:,:) - Ez, Isource(3,:,:) - Hy,
Isource(4,:,:) - Hz,Isource(5,:,:) - Ey, Isource(6,:,:) - Ez, Isource(7,:,:) - Hy,
Isource(8,:,:) - Hz
Jsource - A 3d matrix of dimension 8x(I1-I0+1)x(K1-K0+1) which contains
complex amplitude information for the interface equations at the J0 and J1
planes. Should be Jsource(1,:,:) - Ex, Jsource(2,:,:) - Ez, Jsource(3,:,:) - Hx,
Jsource(4,:,:) - Hz,Jsource(5,:,:) - Ex, Jsource(6,:,:) - Ez, Jsource(7,:,:) - Hx,
Jsource(8,:,:) - Hz
Ksource - A 3d matrix of dimension 8x(I1-I0+1)x(J1-J0+1) which contains
complex amplitude information for the interface equations at the K0 and K1
planes. Should be Ksource(1,:,:) - Ex, Ksource(2,:,:) - Ey, Ksource(3,:,:) - Hx,
Ksource(4,:,:) - Hy,Ksource(5,:,:) - Ex, Ksource(6,:,:) - Ey, Ksource(7,:,:) - Hx,
Ksource(8,:,:) - Hy
grid_labels - A structure with 3 elements, represents the axis labels:
x_grid_label - cartesian labels of Yee cell origins (x-coordinates)
y_grid_label - cartesian labels of Yee cell origins (y-coordinates)
z_grid_label - cartesian labels of Yee cell origins (z-coordinates)
omega_an - analytic angular frequency source
to_l - time delay of pulse
hwhm - hwhm of pulse
Dxl - thickness of upper pml in the x direction
Dxu - thickness of upper pml in the x direction
Dyl - thickness of upper pml in the y direction
Dyu - thickness of upper pml in the y direction
Dzl - thickness of upper pml in the z direction
Dzu - thickness of upper pml in the z direction
lower_boundary_update - boolean for applying the lower update equations
Nt - the number of iterations
start_tind - the starting value of the tind index which controls the current time
step of each fdtd iteration. The default for this would be 0;
sourcemode - integer of value 0 for s steady state simulation and value 1 for a pulse simulation
runmode - integer of value 0 for "complete" and 1 for "analyse"
exphasorsvolume - Extract phasors in the FDTD volume if set to true
exphasorssurface - Extract phasors about a specified cuboid surface if set to true
intphasorssurface - Interpolate surface phasors onto a common point if true
phasorsurface - A list of indices defining the cuboid to extract the phasors at
phasorinc - An integer vector of three elements describing the factor by which to reduce the
density of vertices in the enclosing observation mesh
dimension - A string of value "3", "TE" or "TM"
conductive_aux - auxilliary parameters required to model conductive multilayer
dispersive_aux - auxilliary parameters required to model dispersive multilayer
structure - 2 x (I_tot+1) integer array describing the grating stucture, if one is present
f_ex_vec - 1xN or Nx1 vector of frequencies at which to perform the extraction of complex amplitudes
tdfield - structure containing elements exi and eyi which have dimension (I1-I0+1)x(J1-J0+1)xNt
fieldsample.i - indices along the x-direction of locations at which to sample the field
fieldsample.j - indices along the y-direction of locations at which to sample the field
fieldsample.k - indices along the z-direction of locations at which to sample the field
fieldsample.n - vector of the moments of the field to sample
campssample.vertices - N x 3 matrix of indices where the complex amplitudes will be sampled
campssample.components - numerical array of up to six elements which defines which field components
will be sampled, 1 means Ex, 2 Ey etc.
*/
/*void mexFunction(int nlhs, mxArray *plhs[], int nrhs,const mxArray *prhs[])
The principle of structuring the program in this way is to be able to compile both a mex-file
and an executable. This is a monolithic function which performs the entire FDTD simulation. It
does farm out tasks to other functions however much of the core functionallity is found here. This
is certainly not an example of good programming style however it reduces the amount of argument
passing.
*/
void mexfprintf(const char err[]){
fprintf(stderr,(char *)err);
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/*Local variables*/
#ifdef FDFLAG
int useCD = 1; //determines whether FDTD or PSTD is used
#else
int useCD = 0;
#endif
double time_0, time_1, time_ml_1, time_ml_0;
double secs;
complex<double> commonPhase;
complex<double> *E_norm;
complex<double> *H_norm;
complex<double> E_norm_an=0.;
complex<double> H_norm_an=0.;
double fte = 0., fth = 0.;
double commonAmplitude;
double *x_grid_labels, *y_grid_labels, *z_grid_labels, *x_grid_labels_out, *y_grid_labels_out, *z_grid_labels_out;
double ***Exy, ***Exz, ***Eyx, ***Eyz, ***Ezx, ***Ezy, ***Hxy, ***Hxz, ***Hyx, ***Hyz, ***Hzx, ***Hzy;
double ***Exy_nm1, ***Exz_nm1, ***Eyx_nm1, ***Eyz_nm1, ***Ezx_nm1, ***Ezy_nm1;
double ***Jxy, ***Jxz, ***Jyx, ***Jyz, ***Jzx, ***Jzy, ***Jcxy, ***Jcxz, ***Jcyx, ***Jcyz, ***Jczx, ***Jczy;
double ***Jxy_nm1, ***Jxz_nm1, ***Jyx_nm1, ***Jyz_nm1, ***Jzx_nm1, ***Jzy_nm1;
double ***ExR, ***ExI, ***EyR, ***EyI, ***EzR, ***EzI,***HxR, ***HxI, ***HyR, ***HyI, ***HzR, ***HzI;
double ***ExR2, ***ExI2, ***EyR2, ***EyI2, ***EzR2, ***EzI2;
double **ex_tdf;
mxArray *ex_tdf_array;
double Ex_temp, Ey_temp, Ez_temp;
double ***exi, ***eyi;
int exi_present, eyi_present;
double *I0, *I1, *J0, *J1, *K0, *K1;
double ***IsourceI, ***JsourceI, ***KsourceI, ***IsourceR, ***JsourceR, ***KsourceR;
double ***surface_EHr,***surface_EHi;
double *alpha, *beta, *gamma;
double *ml_alpha, *ml_beta, *ml_gamma, *ml_kappa_x, *ml_kappa_y, *ml_kappa_z, *ml_sigma_x, *ml_sigma_y, *ml_sigma_z;
double *rho_x, *rho_y, *rho_z, rho;
double alpha_l, beta_l, gamma_l;
double kappa_l, sigma_l;
double dx, dy, dz;
double t0;
double *Cmaterial_Cax,*Cmaterial_Cay,*Cmaterial_Caz, *Cmaterial_Cbx, *Cmaterial_Cby, *Cmaterial_Cbz, *Cmaterial_Ccx, *Cmaterial_Ccy, *Cmaterial_Ccz, *Dmaterial_Dax, *Dmaterial_Day, *Dmaterial_Daz, *Dmaterial_Dbx, *Dmaterial_Dby, *Dmaterial_Dbz;//non free space material parameters
double *freespace_Cbx, *freespace_Cby __attribute__ ((unused)),*freespace_Cbz __attribute__ ((unused)),*freespace_Dbx __attribute__ ((unused)), *freespace_Dby __attribute__ ((unused)),*freespace_Dbz __attribute__ ((unused));//freespace variables
double Ca, Cb, Cc;//used by interpolation scheme
double *Cax, *Cay, *Caz,*Cbx, *Cby, *Cbz,*Ccx, *Ccy, *Ccz,*Dax, *Day, *Daz,*Dbx, *Dby, *Dbz;
double *f_ex_vec;
int N_f_ex_vec;
//the C and D vars for free space and pml
double Enp1, Jnp1;
//these are used for boot strapping. There is currently no way of exporting this.
double **iwave_lEx_Rbs, **iwave_lEy_Rbs, **iwave_lHx_Rbs, **iwave_lHy_Rbs,**iwave_lEx_Ibs, **iwave_lEy_Ibs, **iwave_lHx_Ibs, **iwave_lHy_Ibs;
double *to_l, *hwhm, *omega_an, *dt;
double maxfield = 0, tempfield;
double *place_holder;
double *array_ptr_dbl;
int intmatprops=1;//means the material properties will be interpolated
int intmethod;//method of interpolating surface field quantities
double *fieldsample_i, *fieldsample_j, *fieldsample_k, *fieldsample_n;
int N_fieldsample_i, N_fieldsample_j, N_fieldsample_k, N_fieldsample_n;
double air_interface;
int air_interface_present;
//refractive index of the first layer of the multilayer, or of the bulk of homogeneous
double refind;
//PSTD storage
double **ca_vec, **cb_vec, **cc_vec;
fftw_complex **eh_vec, *dk_e_x, *dk_e_y, *dk_e_z, *dk_h_x, *dk_h_y, *dk_h_z;
fftw_plan *pb_exy, *pf_exy, *pb_exz, *pf_exz, *pb_eyx, *pf_eyx, *pb_eyz, *pf_eyz, *pb_ezx, *pf_ezx,*pb_ezy, *pf_ezy, *pb_hxy, *pf_hxy, *pb_hxz, *pf_hxz, *pb_hyx, *pf_hyx, *pb_hyz, *pf_hyz, *pb_hzx, *pf_hzx,*pb_hzy, *pf_hzy;
fftw_plan pex_t, pey_t;
int N_e_x, N_e_y, N_e_z, N_h_x, N_h_y, N_h_z;
int exdetintegral;
int Ndetmodes;
complex<double> ***Dx_tilde, ***Dy_tilde;
double phaseTermE;
complex<double> cphaseTermE;
//these are 2d matrices and must be in row-major order, which diffes from Matlab
fftw_complex *Ex_t, *Ey_t;
complex<double> **Ex_t_cm, **Ey_t_cm;
double lambda_an_t;
double ***D_temp_re, ***D_temp_im;
double **Pupil;
int k_det_obs_global;
double *fx_vec, *fy_vec;
int Nfx_vec, Nfy_vec;
double z_obs = 0.;
//end PSTD storage
unsigned char ***materials;
unsigned char *array_ptr_uint8;
// int *lower_boundary_update;
int *Dxl, *Dxu, *Dyl, *Dyu, *Dzl, *Dzu, *Nt;
int i,j,k, material_nlayers, is_disp, is_cond, is_disp_ml = 0;
int k_loc;
int tind;
int input_counter = 0;
int is_multilayer = 0;
int start_tind;
int pind_il, pind_iu, pind_jl, pind_ju, pind_kl, pind_ku;
int cuboid[6];
int **vertices;
int nvertices;
int *components;
int ncomponents;
double ***camplitudesR, ***camplitudesI;
mxArray *mx_camplitudes;
int sourcemode = sm_steadystate;//0 - steadystate, 1 - pulsed
int runmode = rm_complete;//0 - complete, 1 - analyse
int exphasorsvolume, exphasorssurface, intphasorssurface;
int phasorinc[3];
int dimension = 0;
int num_fields = 0;
int ndims;
int **structure, is_structure = 0;
int I_tot, J_tot, K_tot, K, max_IJK;
int Nsteps = 0, dft_counter = 0;
int **surface_vertices, n_surface_vertices;
int poutfile = 0;
int Np=0; //The phasor extraction algorithm will be executed every Np iterations.
int Npe=0; //The number of terms in the algorithm to extract the phasors
double dtp=0.; //The phasor extraction time step
int Ni_tdf=0, Nk_tdf=0;
#ifdef FDFLAG
int skip_tdf = 6;
#else
int skip_tdf = 1;
#endif
const mwSize *dimptr_out;
mwSize *dims;dims = (mwSize *)malloc(3*sizeof(mwSize));
mwSize *label_dims;label_dims = (mwSize *)malloc(2*sizeof(mwSize));
mxArray *dummy_array[3];
mxArray *element;
mxArray *mx_surface_vertices, *mx_surface_facets, *mx_surface_amplitudes;
mxArray *mx_fieldsample, *mx_fieldsample_x, *mx_fieldsample_y, *mx_fieldsample_z;
double ****fieldsample,****fieldsample_x,****fieldsample_y,****fieldsample_z;
mxArray *mx_Idx, *mx_Idy;
double **Idx_re, **Idx_im, **Idy_re, **Idy_im;
complex<double> **Idx, **Idy;
complex<double> Idxt,Idyt,kprop;
char message_buffer[500];
char dimension_str[3];
const char fdtdgrid_elements[][15] = {"Exy","Exz","Eyx","Eyz","Ezx","Ezy","Hxy","Hxz","Hyx","Hyz","Hzx","Hzy","materials"};
const char Cmaterial_elements[][10] = {"Cax","Cay","Caz","Cbx","Cby","Cbz","Ccx","Ccy","Ccz"};
const char Dmaterial_elements[][10] = {"Dax","Day","Daz","Dbx","Dby","Dbz"};
const char C_elements[][10] = {"Cax","Cay","Caz","Cbx","Cby","Cbz"};
const char C_elements_disp_ml[][10] = {"Cax","Cay","Caz","Cbx","Cby","Cbz","Ccx","Ccy","Ccz"};
const char D_elements[][10] = {"Dax","Day","Daz","Dbx","Dby","Dbz"};
const char freespace_elements[][10] = {"Cbx","Cby","Cbz","Dbx","Dby","Dbz"};
const char disp_params_elements[][10] = {"alpha","beta","gamma"};
const char conductive_aux_elements[][10] = {"rho_x","rho_y","rho_z"};
const char dispersive_aux_elements[][10] = {"alpha","beta","gamma","kappa_x","kappa_y","kappa_z","sigma_x","sigma_y","sigma_z"};
const char delta_elements[][10] = {"x","y","z"};
const char interface_fields[][5] = {"I0","I1","J0","J1","K0","K1"};
const char grid_labels_fields[][15] = {"x_grid_labels","y_grid_labels","z_grid_labels"};
const char fieldsample_elements[][2] = {"i","j","k","n"};
const char campssample_elements[][15] = {"vertices","components"};
char *sourcemodestr;
char *tdfdirstr;
FILE *outfile;
if(poutfile){
outfile = fopen("out.1.2.txt","w");
}
// FILE *eyfile;
// FILE *jyfile;
// eyfile = fopen("Eyz.txt","w");
// jyfile = fopen("Jyz.txt","w");
if( nrhs != 49 ){
fprintf(stderr,"%d\n",nrhs);
mexfprintf("Expected 49 inputs.");
}
if (nlhs != 27)
mexfprintf("27 outputs required.");
/*Get fdtdgrid*/
if(mxIsStruct(prhs[input_counter])){
num_fields = mxGetNumberOfFields(prhs[input_counter]);
//check that all fields are present
if(num_fields != 13){
sprintf(message_buffer, "fdtdgrid should have 13 members, it only has %d",num_fields);
mexfprintf(message_buffer);
}
//now loop over the fields
for(int i=0;i<num_fields;i++){
//element = mxGetField(prhs[input_counter], 0, fdtdgrid_elements[i]);
element = mxGetField( (mxArray *)prhs[input_counter], 0, fdtdgrid_elements[i]);
if(mxIsDouble(element))
array_ptr_dbl = mxGetPr(element);
else if(mxIsUint8(element)){
array_ptr_uint8 = (unsigned char *)mxGetPr(element);
}
else{
sprintf(message_buffer, "Incorrect data type in fdtdgrid.%s",fdtdgrid_elements[i]);
mexfprintf(message_buffer);
}
ndims = mxGetNumberOfDimensions(element);
dimptr_out = mxGetDimensions(element);
if( (ndims != 2) && (ndims != 3) ){
sprintf(message_buffer, "field matrix %s should be 2- or 3-dimensional",fdtdgrid_elements[i]);
mexfprintf(message_buffer);
}
//start
if(!strcmp(fdtdgrid_elements[i],"Exy")){
if(ndims==2)
Exy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else{
//fprintf(stderr,"Dims Exy=%d %d %d\n",dimptr_out[0], dimptr_out[1], dimptr_out[2]);
Exy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
}
else if(!strcmp(fdtdgrid_elements[i],"Exz")){
if(ndims==2)
Exz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Exz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Eyx")){
if(ndims==2)
Eyx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Eyx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Eyz")){
if(ndims==2)
Eyz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Eyz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Ezx")){
if(ndims==2)
Ezx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Ezx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Ezy")){
if(ndims==2)
Ezy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Ezy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Hxy")){
if(ndims==2)
Hxy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Hxy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Hxz")){
if(ndims==2)
Hxz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Hxz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Hyx")){
if(ndims==2)
Hyx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Hyx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Hyz")){
if(ndims==2)
Hyz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Hyz = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Hzx")){
if(ndims==2)
Hzx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Hzx = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"Hzy")){
if(ndims==2)
Hzy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], 0);
else
Hzy = castMatlab3DArray(array_ptr_dbl, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
}
else if(!strcmp(fdtdgrid_elements[i],"materials")){
if(ndims==2)
materials = castMatlab3DArrayUint8(array_ptr_uint8, dimptr_out[0], dimptr_out[1], 0);
else
materials = castMatlab3DArrayUint8(array_ptr_uint8, dimptr_out[0], dimptr_out[1], dimptr_out[2]);
//save this for later when freeing memory
material_nlayers = dimptr_out[2];
I_tot = dimptr_out[0]-1;//The _tot variables do NOT include the additional cell at the
J_tot = dimptr_out[1]-1;//edge of the grid which is only partially used
if(ndims==2)
K_tot = 0;
else
K_tot = dimptr_out[2]-1;
}
else{
sprintf(message_buffer, "element fdtdgrid.%s not handled",fdtdgrid_elements[i]);
mexfprintf(message_buffer);
}
}//end
input_counter++;
}
else{
sprintf(message_buffer, "Argument %d was expected to be a structure",input_counter);
mexfprintf(message_buffer);
}
/*Got fdtdgrid*/
//fprintf(stderr,"Got fdtdgrid\n");
/*Get Cmaterials */
if(mxIsStruct(prhs[input_counter])){
num_fields = mxGetNumberOfFields(prhs[input_counter]);
//check that all fields are present
if(num_fields != 9){
sprintf(message_buffer, "Cmaterials should have 9 members, it has %d",num_fields);
mexfprintf(message_buffer);
}
for(int i=0;i<9;i++){
element = mxGetField( (mxArray *)prhs[input_counter], 0, Cmaterial_elements[i]);
ndims = mxGetNumberOfDimensions(element);
if( ndims == 2 ){
dimptr_out = mxGetDimensions(element);
if(dimptr_out[0] != 1){
sprintf(message_buffer, "Incorrect dimension on Cmaterial.%s",Cmaterial_elements[i]);
mexfprintf(message_buffer);
}
if(!strcmp(Cmaterial_elements[i],"Cax")){
Cmaterial_Cax = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Cay")){
Cmaterial_Cay = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Caz")){
Cmaterial_Caz = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Cbx")){
Cmaterial_Cbx = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Cby")){
Cmaterial_Cby = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Cbz")){
Cmaterial_Cbz = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Ccx")){
Cmaterial_Ccx = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Ccy")){
Cmaterial_Ccy = mxGetPr(element);
}
else if(!strcmp(Cmaterial_elements[i],"Ccz")){
Cmaterial_Ccz = mxGetPr(element);
}
else{
sprintf(message_buffer, "element Cmaterial.%s not handled",Cmaterial_elements[i]);
mexfprintf(message_buffer);
}
}
else
mexfprintf("Incorrect dimension on Cmaterial.Ca");
}
input_counter++;
}
else{
sprintf(message_buffer, "Argument %d was expected to be a structure",input_counter);
mexfprintf(message_buffer);
}
/*Got Cmaterials */
//fprintf(stderr,"Got Cmaterials\n");
/*Get Dmaterials */
if(mxIsStruct(prhs[input_counter])){
num_fields = mxGetNumberOfFields(prhs[input_counter]);
//check that all fields are present
if(num_fields != 6){
sprintf(message_buffer, "Dmaterials should have 6 members, it has %d",num_fields);
mexfprintf(message_buffer);
}
for(int i=0;i<6;i++){
element = mxGetField( (mxArray *)prhs[input_counter], 0, Dmaterial_elements[i]);
ndims = mxGetNumberOfDimensions(element);
if( ndims == 2 ){
dimptr_out = mxGetDimensions(element);
if(dimptr_out[0] != 1){
sprintf(message_buffer, "Incorrect dimension on Dmaterial.%s",Dmaterial_elements[i]);
mexfprintf(message_buffer);
}
if(!strcmp(Dmaterial_elements[i],"Dax")){
Dmaterial_Dax = mxGetPr(element);
}
else if(!strcmp(Dmaterial_elements[i],"Day")){
Dmaterial_Day = mxGetPr(element);
}
else if(!strcmp(Dmaterial_elements[i],"Daz")){
Dmaterial_Daz = mxGetPr(element);
}
else if(!strcmp(Dmaterial_elements[i],"Dbx")){
Dmaterial_Dbx = mxGetPr(element);
}
else if(!strcmp(Dmaterial_elements[i],"Dby")){
Dmaterial_Dby = mxGetPr(element);
}
else if(!strcmp(Dmaterial_elements[i],"Dbz")){
Dmaterial_Dbz = mxGetPr(element);
}
else{
sprintf(message_buffer, "element Dmaterial.%s not handled",Dmaterial_elements[i]);
mexfprintf(message_buffer);
}
}
else
mexfprintf("Incorrect dimension on Dmaterial.Da");
}
input_counter++;
}
else{
sprintf(message_buffer, "Argument %d was expected to be a structure",input_counter);
mexfprintf(message_buffer);
}
/*Got Dmaterials */
//fprintf(stderr,"Got Dmaterials\n");
/*Get C */
if(mxIsStruct(prhs[input_counter])){
num_fields = mxGetNumberOfFields(prhs[input_counter]);
//check that all fields are present
if( (num_fields != 6) && (num_fields != 9) ){
sprintf(message_buffer, "C should have 6 or 9 members, it has %d",num_fields);
mexfprintf(message_buffer);
}
if( num_fields==9 ){
is_disp_ml = 1;
}
if( num_fields==6 ){
for(int i=0;i<num_fields;i++){
element = mxGetField( (mxArray *)prhs[input_counter], 0, C_elements[i]);
ndims = mxGetNumberOfDimensions(element);
if( ndims == 2 ){
dimptr_out = mxGetDimensions(element);
if(dimptr_out[0] != 1){
is_multilayer = 1;
}
if(!strcmp(C_elements[i],"Cax")){
Cax = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Cay")){
Cay = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Caz")){
Caz = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Cbx")){
Cbx = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Cby")){
Cby = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Cbz")){
Cbz = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Ccx")){
Ccx = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Ccy")){
Ccy = mxGetPr(element);
}
else if(!strcmp(C_elements[i],"Ccz")){
Ccz = mxGetPr(element);
}
else{
sprintf(message_buffer, "element C.%s not handled",C_elements[i]);
mexfprintf(message_buffer);
}
}
else
mexfprintf("Incorrect dimension on C");
}
}
else{
for(int i=0;i<num_fields;i++){
element = mxGetField( (mxArray *)prhs[input_counter], 0, C_elements_disp_ml[i]);
ndims = mxGetNumberOfDimensions(element);
if( ndims == 2 ){
dimptr_out = mxGetDimensions(element);
if(dimptr_out[0] != 1){
//sprintf(message_buffer, "Incorrect dimension on C.%s",C_elements[i]);
//mexfprintf(message_buffer);
is_multilayer = 1;
}
if(!strcmp(C_elements_disp_ml[i],"Cax")){
Cax = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Cay")){
Cay = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Caz")){
Caz = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Cbx")){
Cbx = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Cby")){
Cby = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Cbz")){
Cbz = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Ccx")){
Ccx = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Ccy")){
Ccy = mxGetPr(element);
}
else if(!strcmp(C_elements_disp_ml[i],"Ccz")){
Ccz = mxGetPr(element);
}
else{
sprintf(message_buffer, "element C.%s not handled",C_elements_disp_ml[i]);
mexfprintf(message_buffer);
}
}
else
mexfprintf("Incorrect dimension on C");
}
}
input_counter++;
}
else{
sprintf(message_buffer, "Argument %d was expected to be a structure",input_counter);
mexfprintf(message_buffer);
}
/*Got C */
//fprintf(stderr,"Got C\n");
/*Get D */
if(mxIsStruct(prhs[input_counter])){
num_fields = mxGetNumberOfFields(prhs[input_counter]);
//check that all fields are present
if(num_fields != 6){
sprintf(message_buffer, "D should have 6 members, it has %d",num_fields);
mexfprintf(message_buffer);
}
for(int i=0;i<6;i++){
element = mxGetField( (mxArray *)prhs[input_counter], 0, D_elements[i]);
ndims = mxGetNumberOfDimensions(element);
if( ndims == 2 ){
dimptr_out = mxGetDimensions(element);
if(dimptr_out[0] != 1){
//sprintf(message_buffer, "Incorrect dimension on D.%s",D_elements[i]);
//mexfprintf(message_buffer);
}
if(!strcmp(D_elements[i],"Dax")){
Dax = mxGetPr(element);
}
else if(!strcmp(D_elements[i],"Day")){
Day = mxGetPr(element);
}
else if(!strcmp(D_elements[i],"Daz")){
Daz = mxGetPr(element);
}
else if(!strcmp(D_elements[i],"Dbx")){
Dbx = mxGetPr(element);
}
else if(!strcmp(D_elements[i],"Dby")){
Dby = mxGetPr(element);
}
else if(!strcmp(D_elements[i],"Dbz")){
Dbz = mxGetPr(element);
}
else{
sprintf(message_buffer, "element D.%s not handled",D_elements[i]);
mexfprintf(message_buffer);
}
}
else
mexfprintf("Incorrect dimension on D");
}
input_counter++;
}
else{
sprintf(message_buffer, "Argument %d was expected to be a structure",input_counter);
mexfprintf(message_buffer);
}
/*Got D */
//fprintf(stderr,"Got D\n");
/*Get freespace*/
if(mxIsStruct(prhs[input_counter])){
num_fields = mxGetNumberOfFields(prhs[input_counter]);
//check that all fields are present
if(num_fields != 6){
sprintf(message_buffer, "freespace should have 6 members, it has %d",num_fields);
mexfprintf(message_buffer);
}
for(int i=0;i<6;i++){
element = mxGetField( (mxArray *)prhs[input_counter], 0, freespace_elements[i]);
ndims = mxGetNumberOfDimensions(element);
if( ndims == 2 ){
dimptr_out = mxGetDimensions(element);
if(dimptr_out[0] != 1){
sprintf(message_buffer, "Incorrect dimension on freespace.%s",freespace_elements[i]);
mexfprintf(message_buffer);
}
if(!strcmp(freespace_elements[i],"Cbx")){
freespace_Cbx = mxGetPr(element);
}
else if(!strcmp(freespace_elements[i],"Cby")){
freespace_Cby = mxGetPr(element);
}
else if(!strcmp(freespace_elements[i],"Cbz")){
freespace_Cbz = mxGetPr(element);
}
else if(!strcmp(freespace_elements[i],"Dbx")){
freespace_Dbx = mxGetPr(element);
}
else if(!strcmp(freespace_elements[i],"Dby")){
freespace_Dby = mxGetPr(element);
}
else if(!strcmp(freespace_elements[i],"Dbz")){
freespace_Dbz = mxGetPr(element);
}
else{
sprintf(message_buffer, "element freespace.%s not handled",freespace_elements[i]);
mexfprintf(message_buffer);
}
}
else
mexfprintf("Incorrect dimension on freespace");
}
input_counter++;
}
else{
sprintf(message_buffer, "Argument %d was expected to be a structure",input_counter);
mexfprintf(message_buffer);
}
/*Got freespace*/
//fprintf(stderr,"Got freespace\n");
/*Get disp_params */
if(mxIsStruct(prhs[input_counter])){
num_fields = mxGetNumberOfFields(prhs[input_counter]);
//check that all fields are present
if(num_fields != 3){
sprintf(message_buffer, "disp_params should have 3 members, it has %d",num_fields);
mexfprintf(message_buffer);
}
for(int i=0;i<3;i++){
element = mxGetField( (mxArray *)prhs[input_counter], 0, disp_params_elements[i]);
ndims = mxGetNumberOfDimensions(element);
if( ndims == 2 ){
dimptr_out = mxGetDimensions(element);
if( !(dimptr_out[0] == 1 ||dimptr_out[0] == 0) ){
sprintf(message_buffer, "Incorrect dimension on disp_params.%s",disp_params_elements[i]);
mexfprintf(message_buffer);
}
if(!strcmp(disp_params_elements[i],"alpha")){
alpha = mxGetPr(element);
}
else if(!strcmp(disp_params_elements[i],"beta")){
beta = mxGetPr(element);
}
else if(!strcmp(disp_params_elements[i],"gamma")){
gamma = mxGetPr(element);
}
else{
sprintf(message_buffer, "element disp_params.%s not handled",disp_params_elements[i]);
mexfprintf(message_buffer);
}
}
else
mexfprintf("Incorrect dimension on disp_params");
}
input_counter++;
}
else{
sprintf(message_buffer, "Argument %d was expected to be a structure",input_counter);
mexfprintf(message_buffer);
}