Contents of file 'magbragg/magbragg.c':
1 /*45:*/
2 #line 3404 "./magbragg.w"
3
4 /*46:*/
5 #line 3437 "./magbragg.w"
6
7 #include <math.h>
8 #include <time.h>
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <string.h>
12 #include <ctype.h>
13
14 /*:46*/
15 #line 3405 "./magbragg.w"
16
17 /*47:*/
18 #line 3454 "./magbragg.w"
19
20 #define VERSION "1.43"
21 #define COPYRIGHT "Copyright (C) 2002-2007, Fredrik Jonsson"
22 #define SUCCESS (0)
23 #define FAILURE (1)
24 #define NCHMAX (256)
25
26 /*:47*/
27 #line 3406 "./magbragg.w"
28
29 /*48:*/
30 #line 3466 "./magbragg.w"
31
32 typedef struct DCOMPLEX{double r,i;}dcomplex;
33
34 /*:48*/
35 #line 3407 "./magbragg.w"
36
37 /*49:*/
38 #line 3475 "./magbragg.w"
39
40 extern char*optarg;
41 char*progname;
42
43 /*:49*/
44 #line 3408 "./magbragg.w"
45
46 /*50:*/
47 #line 3481 "./magbragg.w"
48
49 /*88:*/
50 #line 5491 "./magbragg.w"
51
52 int numeric(char ch){
53 if((ch=='0')||(ch=='1')||(ch=='2')||(ch=='3')||(ch=='4')||(ch=='5')
54 ||(ch=='6')||(ch=='7')||(ch=='8')||(ch=='9')||(ch=='+')||(ch=='-')){
55 return(1);
56 }else{
57 return(0);
58 }
59 }
60
61 /*:88*/
62 #line 3482 "./magbragg.w"
63
64 /*89:*/
65 #line 5503 "./magbragg.w"
66
67 void init_cantor_fractal_grating(double*z,int indxmin,int indxmax,
68 double llmin,double llmax,double n1,double n2){
69 int indxmid;
70 double dll;
71 if(indxmax-indxmin==1){
72 z[indxmin]= llmin;
73 z[indxmax]= llmax;
74 }else if(indxmax-indxmin>=3){
75 indxmid= (indxmin+indxmax)/2;
76 dll= (n2/(n1+2*n2))*(llmax-llmin);
77 init_cantor_fractal_grating(z,indxmin,indxmid,llmin,llmin+dll,n1,n2);
78 init_cantor_fractal_grating(z,indxmid+1,indxmax,llmax-dll,llmax,n1,n2);
79 }else{
80 fprintf(stderr,"%s: Error in indexes detected in routine %s",progname,
81 "init_cantor_fractal_grating()!\n");
82 fprintf(stderr,"%s: Please verify that the number of discrete layers\n",
83 progname);
84 fprintf(stderr,"%s: in the grating is N-1, where N is an integer"
85 " power of 2.\n",progname);
86 exit(FAILURE);
87 }
88 }
89
90 /*:89*/
91 #line 3483 "./magbragg.w"
92
93 /*90:*/
94 #line 5536 "./magbragg.w"
95
96 /*91:*/
97 #line 5545 "./magbragg.w"
98
99 short pathcharacter(int ch){
100 return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
101 (ch=='-')||(ch=='+'));
102 }
103
104 /*:91*/
105 #line 5537 "./magbragg.w"
106
107 /*92:*/
108 #line 5557 "./magbragg.w"
109
110 char*strip_away_path(char filename[]){
111 int j,k= 0;
112 while(pathcharacter(filename[k]))k++;
113 j= (--k);
114 while(isalnum((int)(filename[j])))j--;
115 j++;
116 return(&filename[j]);
117 }
118
119 /*:92*/
120 #line 5538 "./magbragg.w"
121
122
123 /*:90*/
124 #line 3484 "./magbragg.w"
125
126 /*111:*/
127 #line 5895 "./magbragg.w"
128
129 /*112:*/
130 #line 5906 "./magbragg.w"
131
132 double*dvector(long nl,long nh){
133 double*v;
134 v= (double*)malloc((size_t)((nh-nl+2)*sizeof(double)));
135 if(!v){
136 fprintf(stderr,"Error: Allocation failure in dvector()\n");
137 exit(FAILURE);
138 }
139 return v-nl+1;
140 }
141
142 /*:112*/
143 #line 5896 "./magbragg.w"
144
145 /*113:*/
146 #line 5922 "./magbragg.w"
147
148 dcomplex*dcvector(long nl,long nh){
149 dcomplex*v;
150 v= (dcomplex*)malloc((size_t)((nh-nl+2)*sizeof(dcomplex)));
151 if(!v){
152 fprintf(stderr,"Error: Allocation failure in dcvector()\n");
153 exit(FAILURE);
154 }
155 return v-nl+1;
156 }
157
158 /*:113*/
159 #line 5897 "./magbragg.w"
160
161 /*114:*/
162 #line 5936 "./magbragg.w"
163
164 void free_dvector(double*v,long nl,long nh){
165 free((char*)(v+nl-1));
166 }
167
168 /*:114*/
169 #line 5898 "./magbragg.w"
170
171 /*115:*/
172 #line 5944 "./magbragg.w"
173
174 void free_dcvector(dcomplex*v,long nl,long nh){
175 free((char*)(v+nl-1));
176 }
177
178 /*:115*/
179 #line 5899 "./magbragg.w"
180
181
182 /*:111*/
183 #line 3485 "./magbragg.w"
184
185 /*93:*/
186 #line 5571 "./magbragg.w"
187
188 #define IA 16807
189 #define IM 2147483647
190 #define AM (1.0/IM)
191 #define IQ 127773
192 #define IR 2836
193 #define NTAB 32
194 #define NDIV (1+(IM-1)/NTAB)
195 #define EPS 1.2e-7
196 #define RNMX (1.0-EPS)
197
198 float ran1(long*idum)
199 {
200 int j;
201 long k;
202 static long iy= 0;
203 static long iv[NTAB];
204 float temp;
205
206 if(*idum<=0||!iy){
207 if(-(*idum)<1)*idum= 1;
208 else*idum= -(*idum);
209 for(j= NTAB+7;j>=0;j--){
210 k= (*idum)/IQ;
211 *idum= IA*(*idum-k*IQ)-IR*k;
212 if(*idum<0)*idum+= IM;
213 if(j<NTAB)iv[j]= *idum;
214 }
215 iy= iv[0];
216 }
217 k= (*idum)/IQ;
218 *idum= IA*(*idum-k*IQ)-IR*k;
219 if(*idum<0)*idum+= IM;
220 j= iy/NDIV;
221 iy= iv[j];
222 iv[j]= *idum;
223 if((temp= AM*iy)> RNMX)return RNMX;
224 else return temp;
225 }
226 #undef IA
227 #undef IM
228 #undef AM
229 #undef IQ
230 #undef IR
231 #undef NTAB
232 #undef NDIV
233 #undef EPS
234 #undef RNMX
235
236 /*:93*/
237 #line 3486 "./magbragg.w"
238
239 /*94:*/
240 #line 5626 "./magbragg.w"
241
242 /*95:*/
243 #line 5643 "./magbragg.w"
244
245 dcomplex complex(double re,double im){
246 dcomplex c;
247 c.r= re;
248 c.i= im;
249 return c;
250 }
251
252 /*:95*/
253 #line 5627 "./magbragg.w"
254
255 /*96:*/
256 #line 5655 "./magbragg.w"
257
258 dcomplex conjg(dcomplex z){
259 dcomplex c;
260 c.r= z.r;
261 c.i= -z.i;
262 return c;
263 }
264
265 /*:96*/
266 #line 5628 "./magbragg.w"
267
268 /*97:*/
269 #line 5667 "./magbragg.w"
270
271 double cdabs(dcomplex z){
272 double x,y,c,tmp;
273 x= fabs(z.r);
274 y= fabs(z.i);
275 if(x==0.0)
276 c= y;
277 else if(y==0.0)
278 c= x;
279 else if(x> y){
280 tmp= y/x;
281 c= x*sqrt(1.0+tmp*tmp);
282 }else{
283 tmp= x/y;
284 c= y*sqrt(1.0+tmp*tmp);
285 }
286 return c;
287 }
288
289 /*:97*/
290 #line 5629 "./magbragg.w"
291
292 /*98:*/
293 #line 5691 "./magbragg.w"
294
295 double cabs2(dcomplex z){
296 double x,y,c,tmp;
297 x= fabs(z.r);
298 y= fabs(z.i);
299 if(x==0.0)
300 c= y*y;
301 else if(y==0.0)
302 c= x*x;
303 else if(x> y){
304 tmp= y/x;
305 c= x*x*(1.0+tmp*tmp);
306 }else{
307 tmp= x/y;
308 c= y*y*(1.0+tmp*tmp);
309 }
310 return c;
311 }
312
313 /*:98*/
314 #line 5630 "./magbragg.w"
315
316 /*99:*/
317 #line 5713 "./magbragg.w"
318
319 dcomplex cadd(dcomplex a,dcomplex b){
320 dcomplex c;
321 c.r= a.r+b.r;
322 c.i= a.i+b.i;
323 return c;
324 }
325
326 /*:99*/
327 #line 5631 "./magbragg.w"
328
329 /*100:*/
330 #line 5724 "./magbragg.w"
331
332 dcomplex csub(dcomplex a,dcomplex b){
333 dcomplex c;
334 c.r= a.r-b.r;
335 c.i= a.i-b.i;
336 return c;
337 }
338
339 /*:100*/
340 #line 5632 "./magbragg.w"
341
342 /*101:*/
343 #line 5737 "./magbragg.w"
344
345 /*102:*/
346 #line 5744 "./magbragg.w"
347
348 dcomplex cmul(dcomplex a,dcomplex b){
349 dcomplex c;
350 c.r= a.r*b.r-a.i*b.i;
351 c.i= a.i*b.r+a.r*b.i;
352 return c;
353 }
354
355 /*:102*/
356 #line 5738 "./magbragg.w"
357
358 /*103:*/
359 #line 5756 "./magbragg.w"
360
361 dcomplex rcmul(double a,dcomplex b){
362 dcomplex c;
363 c.r= a*b.r;
364 c.i= a*b.i;
365 return c;
366 }
367
368 /*:103*/
369 #line 5739 "./magbragg.w"
370
371
372 /*:101*/
373 #line 5633 "./magbragg.w"
374
375 /*104:*/
376 #line 5769 "./magbragg.w"
377
378 /*105:*/
379 #line 5779 "./magbragg.w"
380
381 dcomplex cdiv(dcomplex a,dcomplex b){
382 dcomplex c;
383 double r,den;
384 if(cdabs(b)> 0.0){
385 if(fabs(b.i)==0.0){
386 c.r= a.r/b.r;
387 c.i= a.i/b.r;
388 }else{
389 if(fabs(b.r)>=fabs(b.i)){
390 r= b.i/b.r;
391 den= b.r+r*b.i;
392 c.r= (a.r+r*a.i)/den;
393 c.i= (a.i-r*a.r)/den;
394 }else{
395 r= b.r/b.i;
396 den= b.i+r*b.r;
397 c.r= (a.r*r+a.i)/den;
398 c.i= (a.i*r-a.r)/den;
399 }
400 }
401 }else{
402 fprintf(stderr,"Error in cdiv(): Division by zero!\n");
403 exit(FAILURE);
404 }
405 return c;
406 }
407
408 /*:105*/
409 #line 5770 "./magbragg.w"
410
411 /*106:*/
412 #line 5811 "./magbragg.w"
413
414 dcomplex crdiv(dcomplex a,double b){
415 dcomplex c;
416 if(fabs(b)> 0.0){
417 c.r= a.r/b;
418 c.i= a.i/b;
419 }else{
420 fprintf(stderr,"Error in crdiv(): Division by zero!\n");
421 exit(FAILURE);
422 }
423 return c;
424 }
425
426 /*:106*/
427 #line 5771 "./magbragg.w"
428
429
430 /*:104*/
431 #line 5634 "./magbragg.w"
432
433 /*107:*/
434 #line 5827 "./magbragg.w"
435
436 dcomplex csqrt(dcomplex z){
437 dcomplex c;
438 double x,y,w,r;
439 if((z.r==0.0)&&(z.i==0.0)){
440 c.r= 0.0;
441 c.i= 0.0;
442 return c;
443 }else{
444 x= fabs(z.r);
445 y= fabs(z.i);
446 if(x>=y){
447 r= y/x;
448 w= sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
449 }else{
450 r= x/y;
451 w= sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
452 }
453 if(z.r>=0.0){
454 c.r= w;
455 c.i= 0.5*z.i/w;
456 }else{
457 c.i= ((z.i>=0)?w:-w);
458 c.r= 0.5*z.i/c.i;
459 }
460 return c;
461 }
462 }
463
464 /*:107*/
465 #line 5635 "./magbragg.w"
466
467 /*108:*/
468 #line 5862 "./magbragg.w"
469
470 /*109:*/
471 #line 5870 "./magbragg.w"
472
473 dcomplex cexp(dcomplex a){
474 dcomplex c;
475 double tmp= exp(a.r);
476 c.r= tmp*cos(a.i);
477 c.i= tmp*sin(a.i);
478 return c;
479 }
480
481 /*:109*/
482 #line 5863 "./magbragg.w"
483
484 /*110:*/
485 #line 5882 "./magbragg.w"
486
487 dcomplex crexpi(double a){
488 dcomplex c;
489 c.r= cos(a);
490 c.i= sin(a);
491 return c;
492 }
493
494 /*:110*/
495 #line 5864 "./magbragg.w"
496
497
498 /*:108*/
499 #line 5636 "./magbragg.w"
500
501
502 /*:94*/
503 #line 3487 "./magbragg.w"
504
505 /*130:*/
506 #line 7271 "./magbragg.w"
507
508 /*131:*/
509 #line 7282 "./magbragg.w"
510
511 void hl(char firststring[],char secondstring[]){
512 if(strlen(firststring)> 25){
513 fprintf(stderr,"%s:******* Error in hl() routine! *******\n",progname);
514 fprintf(stderr,"%s: The first string argument is too long:\n",progname);
515 fprintf(stderr,"%s: '%s'\n",progname,firststring);
516 fprintf(stderr,"%s: String lengths is %d characters\n",
517 progname,(int)strlen(firststring));
518 fprintf(stderr,"%s: (Maximum 25 characters for first argument.)\n",
519 progname);
520 fprintf(stderr,"%s:******** End of error message *********\n",progname);
521 exit(FAILURE);
522 }
523 if(strlen(secondstring)> 55){
524 fprintf(stderr,"%s:******* Error in hl() routine! *******\n",progname);
525 fprintf(stderr,"%s: The second string argument is too long:\n",progname);
526 fprintf(stderr,"%s: '%s'\n",progname,secondstring);
527 fprintf(stderr,"%s: String lengths is %d characters\n",
528 progname,(int)strlen(secondstring));
529 fprintf(stderr,"%s: (Maximum 55 characters for second argument.)\n",
530 progname);
531 fprintf(stderr,"%s:******** End of error message *********\n",progname);
532 exit(FAILURE);
533 }
534 fprintf(stderr,"%-25.25s%1.55s\n",firststring,secondstring);
535 }
536
537 /*:131*/
538 #line 7272 "./magbragg.w"
539
540 /*132:*/
541 #line 7313 "./magbragg.w"
542
543 void fhl(char linestring[]){
544 if(strlen(linestring)> 80){
545 fprintf(stderr,"%s:******* Error in fhl() routine! *******\n",progname);
546 fprintf(stderr,"%s: The following help line is too long:\n",progname);
547 fprintf(stderr,"%s: '%s'\n",progname,linestring);
548 fprintf(stderr,"%s: String is %d characters\n",
549 progname,(int)strlen(linestring));
550 fprintf(stderr,"%s: (Maximum 80 characters per help line is allowed.)\n",
551 progname);
552 fprintf(stderr,"%s:******** End of error message *********\n",progname);
553 exit(FAILURE);
554 }
555 fprintf(stderr,"%s\n",linestring);
556 }
557
558 /*:132*/
559 #line 7273 "./magbragg.w"
560
561 /*133:*/
562 #line 7332 "./magbragg.w"
563
564 void showsomehelp(void){
565 fprintf(stderr," Usage: %s [options]\n",progname);
566 fhl(" Options:");
567 hl(" -h, --help","Display this help message and exit clean.");
568 hl(" -N <int>","");
569 hl(" -M <int>","");
570 hl(" -v, --verbose","");
571 hl(" -o, --outputfile <str>","");
572 fhl(" --fieldevolution {efield|stoke} <n> <str>");
573 fhl(" --intensityevolution <lambda> <str>");
574 fhl(" --normalize_length_to_um");
575 hl(" --normalize_intensity","");
576 hl("","When saving the spatial evolution of the intra-");
577 hl("","grating intensity, normalize the intensity with");
578 hl("","respect to the intensity at z=0, inside the");
579 hl("","grating (that is to say, normalize with respect");
580 hl("","to the initial intra-grating intensity). This");
581 hl("","option *only* affects the fields saved with the");
582 hl("","--intensityevolution or --fieldevolution");
583 hl("","options.");
584 hl(" -r, --random","");
585 hl(" -a, --apodize <real>","");
586 hl("","Apodize the grating structure over geometrical");
587 hl("","distance <real> at each end of the grating.");
588 hl("","The option only applies to gratings with sinus-");
589 hl("","oidal modulation of refractive index and");
590 hl("","gyration coefficient, in constant or chirped");
591 hl("","periodic configurations, as specified with the");
592 hl("","'--grating sinusoidal' or '--grating chirped'");
593 hl("","options respectively.");
594 hl(""," The apodization is applied to the refractive");
595 hl("","index modulation and, in cases where the linear");
596 hl("","magneto-optical is spatially modulated as well,");
597 hl("","to the linear gyration coefficient.");
598 fhl(" -j, --phasejump <r1 (angle)> <r2 (position)>");
599 hl("","Apply discrete phase jump in the spatial phase");
600 hl("","of the grating profile. This option adds the");
601 hl("","real number <r1> to the argument of the sinus-");
602 hl("","oidal function for the grating profile for all");
603 hl("","spatial coordinates z >= <r2>.");
604 hl(""," As in the case of apodization, this option");
605 hl("","only applies to gratings with sinusoidal modu-");
606 hl("","lation of refractive index and gyration coeffi-");
607 hl("","cient, in constant or chirped periodic configu-");
608 hl("","rations, as specified with the '--grating");
609 hl("","sinusoidal' or '--grating chirped' options");
610 hl("","respectively. The discrete phase jump applies");
611 hl("","to the linear linear refractive index modula-");
612 hl("","tion and, in cases where the linear magneto-");
613 hl("","optical interaction is spatially modulated as");
614 hl("","well, also to the linear gyration coefficient.");
615 fhl(" -w, --writegratingfile <str>");
616 hl(" --spectrumfile <str>","");
617 hl("","Generates the complex reflectance as function");
618 hl("","of the vacuum wavelength in meters, and save");
619 hl("","the spectrum in file named according to the");
620 hl("","supplied character string <str>.");
621 fhl(" --intensityspectrumfile <str>");
622 hl("","Generates the intensity reflectance as function");
623 hl("","of the vacuum wavelength in meters, and save");
624 hl("","the spectrum in file named according to the");
625 hl("","supplied character string <str>.");
626 fhl(" -g, --grating <grating options>");
627 hl("","Specifies the grating type, where");
628 hl("","<grating options> = ");
629 hl(""," [stepwise <stepwise options> |");
630 hl(""," sinusoidal <sinusoidal options> |");
631 hl(""," chirped <chirped options>]");
632 hl("","<stepwise options> = ");
633 hl(""," twolevel t1 <f> t2 <f>");
634 hl(""," n1 <f> n2 <f> g1 <f> g2 <f>");
635 hl(""," pe1 <f> pe2 <f> pm1 <f> pm2 <f>");
636 hl(""," qe1 <f> qe2 <f> qm1 <f> qm2 <f>");
637 hl("","<sinusoidal options> = ");
638 hl(""," n <n0> <dn> <nper>");
639 hl(""," g <g0> <dg> <gper>");
640 hl(""," pe <pe0> <dpe> <peper>");
641 hl(""," pm <pm0> <dpm> <pmper>");
642 hl(""," qe <qe0> <dqe> <qeper>");
643 hl(""," qm <qm0> <dqm> <qmper>");
644 hl("","<chirped options> = ");
645 hl(""," n <n0> <dn> <nper> <ncrp>");
646 hl(""," g <g0> <dg> <gper> <gcrp>");
647 hl(""," pe <pe0> <dpe> <peper> <pecrp>");
648 hl(""," pm <pm0> <dpm> <pmper> <pmcrp>");
649 hl(""," qe <qe0> <dqe> <qeper> <qecrp>");
650 hl(""," qm <qm0> <dqm> <qmper> <qmcrp>");
651 hl(" -L,--gratinglength <f>","Physical length of grating in meter [m]");
652 fhl(" --refindsurr <f>");
653 fhl(" --trmtraject <str>");
654 fhl(" --trmintensity <istart> <istop> <mmi>");
655 fhl(" (intensity measured in Watts per square meter)");
656 fhl(" --trmellipticity <estart> <estop> <mme>");
657 fhl(" --lambdastart <lambda>");
658 fhl(" (start vacuum wavelength measured in meter)");
659 fhl(" --lambdastop <lambda>");
660 fhl(" (stop vacuum wavelength measured in meter)");
661 exit(FAILURE);
662 }
663
664 /*:133*/
665 #line 7274 "./magbragg.w"
666
667
668 /*:130*/
669 #line 3488 "./magbragg.w"
670
671
672 /*:50*/
673 #line 3409 "./magbragg.w"
674
675
676 int main(int argc,char*argv[])
677 {
678 /*51:*/
679 #line 3501 "./magbragg.w"
680
681 /*52:*/
682 #line 3565 "./magbragg.w"
683
684 static double pi= 3.1415926535897932384626433832795;
685 static double twopi= 2.0*3.1415926535897932384626433832795;
686 static double c= 2.997925e8;
687 static double epsilon0= 8.854e-12;
688
689 /*:52*/
690 #line 3502 "./magbragg.w"
691
692 /*53:*/
693 #line 3580 "./magbragg.w"
694
695 time_t initime= time(NULL),now= time(NULL),eta= time(NULL);
696
697 /*:53*/
698 #line 3503 "./magbragg.w"
699
700 /*54:*/
701 #line 3607 "./magbragg.w"
702
703 dcomplex*efp,*efm,*ebp,*ebm;
704
705 /*:54*/
706 #line 3504 "./magbragg.w"
707
708 /*55:*/
709 #line 3618 "./magbragg.w"
710
711 short verbose;
712
713 short randomdistribution;
714
715 short writegratingtofile;
716
717 short scale_stokesparams;
718
719 short normalize_length_to_micrometer;
720 short normalize_intensity;
721 short normalize_ellipticity;
722 short normalize_internally;
723 short odd_layer;
724 short chirpflag;
725
726 short apodize;
727 short phasejump;
728
729 short fieldevoflag;
730 short fieldevoflag_efield;
731
732 short intensityevoflag;
733
734 short fieldevoflag_stoke;
735
736 short intensityinfo;
737
738
739 short saveintensityinfologfile;
740
741 short trmtraject_specified;
742 short save_dbspectra;
743
744 short stokes_parameter_spectrum;
745
746 short display_surrounding_media;
747
748 short perturbed_gyration_constant;
749
750 /*:55*/
751 #line 3505 "./magbragg.w"
752
753 /*56:*/
754 #line 3671 "./magbragg.w"
755
756 long mm,mme,mmi;
757
758 /*:56*/
759 #line 3506 "./magbragg.w"
760
761 /*57:*/
762 #line 3724 "./magbragg.w"
763
764 long nn,modnum;
765 double ll,*z,*dz,*n,*g,*pe,*pm,*qe,*qm,*etafp,*etafm,*etabp,*etabm;
766 double*taup,*taum,*taupp,*taupm,*rhop,*rhom,*rhopp,*rhopm;
767
768 /*:57*/
769 #line 3507 "./magbragg.w"
770
771 /*58:*/
772 #line 3731 "./magbragg.w"
773
774 FILE*fp_s0,*fp_s1,*fp_s2,*fp_s3;
775 FILE*fp_v0,*fp_v1,*fp_v2,*fp_v3;
776 FILE*fp_w0,*fp_w1,*fp_w2,*fp_w3;
777 FILE*fp_evo= NULL,*fp_ievo= NULL,*fp_spec= NULL,*fp_traject= NULL;
778 FILE*fp_irspec= NULL,*fp_itspec= NULL,*fp_icspec= NULL,*fp_gr= NULL;
779 FILE*fp_evo_s0= NULL,*fp_evo_s1= NULL,*fp_evo_s2= NULL,*fp_evo_s3= NULL;
780 FILE*intensinfologfile= NULL;
781
782 /*:58*/
783 #line 3508 "./magbragg.w"
784
785 /*59:*/
786 #line 3745 "./magbragg.w"
787
788 char gratingtype[NCHMAX],gratingsubtype[NCHMAX];
789 char fieldevofilename[NCHMAX],intensityevofilename[NCHMAX];
790 char spectrumfilename[NCHMAX],trmtraject_filename[NCHMAX];
791 char intensity_reflection_spectrumfilename[NCHMAX],
792 intensity_transmission_spectrumfilename[NCHMAX],
793 intensity_check_spectrumfilename[NCHMAX];
794 char fieldevofilename_s0[NCHMAX],fieldevofilename_s1[NCHMAX],
795 fieldevofilename_s2[NCHMAX],fieldevofilename_s3[NCHMAX];
796 char intensinfologfilename[NCHMAX];
797 char outfilename[NCHMAX],gratingfilename[NCHMAX];
798 char outfilename_s0[NCHMAX],outfilename_s1[NCHMAX],
799 outfilename_s2[NCHMAX],outfilename_s3[NCHMAX];
800 char outfilename_v0[NCHMAX],outfilename_v1[NCHMAX],
801 outfilename_v2[NCHMAX],outfilename_v3[NCHMAX];
802 char outfilename_w0[NCHMAX],outfilename_w1[NCHMAX],
803 outfilename_w2[NCHMAX],outfilename_w3[NCHMAX];
804
805 /*:59*/
806 #line 3509 "./magbragg.w"
807
808 /*60:*/
809 #line 3770 "./magbragg.w"
810
811 dcomplex tmpfp,tmpfm,tmpbp,tmpbm;
812 double tmp,zt,s0,s1,s2,s3,v0,v1,v2,v3,w0,w1,w2,w3,stn,
813 maxintens= 0.0,maxintens_inintens= 0.0,maxintens_inellip= 0.0,
814 maxintens_trintens= 0.0,maxintens_trellip= 0.0;
815 int no_arg,status= 0,tmpch;
816 long j,k,ke,ki;
817
818 /*:60*/
819 #line 3510 "./magbragg.w"
820
821 long nne,jje,mmtraject;
822 long ranseed;
823 long maxintens_layer= 0;
824 int fractal_level= 0;
825 int maximum_fractal_level= 0;
826
827 double*w0traj= NULL,*w3traj= NULL;
828 double lambdastart;
829
830 double lambdastop;
831
832 double lambda;
833
834 double omega;
835
836 double ievolambda;
837 double apolength;
838
839 double phasejumpangle;
840
841 double phasejumpposition;
842
843 double phi;
844
845 double gyroperturb_position;
846
847
848 double gyroperturb_amplitude;
849
850
851 double gyroperturb_width;
852
853
854 double nsurr;
855
856 double n1,n2,t1,t2,nper,ncrp,g1,g2,gper,gcrp,
857 pe1,pe2,peper,pecrp,pm1,pm2,pmper,pmcrp,
858 qe1,qe2,qeper,qecrp,qm1,qm2,qmper,qmcrp;
859 double modn1,modt1,modg1,modpe1,modpm1,modqe1,modqm1,
860 aafp2,aafm2,aabp2,aabm2;
861 double trmintensity,trmintenstart,trmintenstop,
862 trmellipticity,trmellipstart,trmellipstop,stoke_scalefactor;
863 static double nndef= 600;
864 static double mmdef= 1000;
865 static double mmedef= 1;
866 static double mmidef= 1;
867 static double lldef= 0.050;
868 static double nsurrdef= 1.0;
869 static double lambdastartdef= 1536.0e-9,lambdastopdef= 1546.0e-9;
870
871 /*:51*/
872 #line 3413 "./magbragg.w"
873
874 /*61:*/
875 #line 3780 "./magbragg.w"
876
877 trmtraject_specified= 0;
878 intensityinfo= 0;
879 saveintensityinfologfile= 0;
880 randomdistribution= 0;
881 writegratingtofile= 0;
882 scale_stokesparams= 0;
883 stoke_scalefactor= 1.0;
884 normalize_length_to_micrometer= 0;
885 normalize_intensity= 0;
886 normalize_ellipticity= 0;
887 ranseed= 1097;
888 modnum= -1;
889 verbose= 0;
890 apodize= 0;
891 phasejump= 0;
892 normalize_internally= 0;
893 odd_layer= 0;
894 chirpflag= 1;
895 save_dbspectra= 0;
896 stokes_parameter_spectrum= 0;
897 display_surrounding_media= 1;
898
899
900 perturbed_gyration_constant= 0;
901 fieldevoflag= 0;
902 fieldevoflag_efield= 0;
903 fieldevoflag_stoke= 0;
904 intensityevoflag= 0;
905 mm= mmdef;
906 nn= nndef;
907 mme= mmedef;
908 mmi= mmidef;
909 nne= 1;
910 ll= lldef;
911 lambdastart= lambdastartdef;
912 lambdastop= lambdastopdef;
913 phasejumpangle= 0.0;
914 phasejumpposition= 0.0;
915 phi= 0.0;
916 nsurr= nsurrdef;
917 strcpy(outfilename,"out.stok");
918 strcpy(fieldevofilename,"out.fevo.dat");
919 strcpy(fieldevofilename_s0,"out.fevo.s0.dat");
920 strcpy(fieldevofilename_s1,"out.fevo.s1.dat");
921 strcpy(fieldevofilename_s2,"out.fevo.s2.dat");
922 strcpy(fieldevofilename_s3,"out.fevo.s3.dat");
923 strcpy(spectrumfilename,"out.rsp.dat");
924 strcpy(intensity_reflection_spectrumfilename,"out.irsp.dat");
925 strcpy(intensity_transmission_spectrumfilename,"out.trsp.dat");
926 strcpy(intensity_check_spectrumfilename,"out.chec.dat");
927 fp_s0= NULL;fp_s1= NULL;fp_s2= NULL;fp_s3= NULL;
928 fp_v0= NULL;fp_v1= NULL;fp_v2= NULL;fp_v3= NULL;
929 fp_w0= NULL;fp_w1= NULL;fp_w2= NULL;fp_w3= NULL;
930
931 /*:61*/
932 #line 3414 "./magbragg.w"
933
934 /*116:*/
935 #line 5955 "./magbragg.w"
936
937 {
938 progname= strip_away_path(argv[0]);
939 fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
940 no_arg= argc;
941 while(--argc){
942 if(!strcmp(argv[no_arg-argc],"-o")||
943 !strcmp(argv[no_arg-argc],"--outputfile")){
944 --argc;
945 strcpy(outfilename,argv[no_arg-argc]);
946 }else if(!strcmp(argv[no_arg-argc],"-w")||
947 !strcmp(argv[no_arg-argc],"--writegratingfile")){
948 --argc;
949 strcpy(gratingfilename,argv[no_arg-argc]);
950 writegratingtofile= 1;
951 }else if(!strcmp(argv[no_arg-argc],"--displaysurrmedia")){
952 display_surrounding_media= (display_surrounding_media?0:1);
953 }else if(!strcmp(argv[no_arg-argc],"-a")||
954 !strcmp(argv[no_arg-argc],"--apodize")){
955 /*118:*/
956 #line 6144 "./magbragg.w"
957
958 {
959 apodize= 1;
960 --argc;
961 if(!sscanf(argv[no_arg-argc],"%lf",&apolength)){
962 fprintf(stderr,"%s: Error in '--apodize' option.\n",progname);
963 exit(FAILURE);
964 }
965 }
966
967 /*:118*/
968 #line 5974 "./magbragg.w"
969 ;
970 }else if(!strcmp(argv[no_arg-argc],"--phasejump")){
971 /*119:*/
972 #line 6157 "./magbragg.w"
973
974 {
975 phasejump= 1;
976 --argc;
977 if(!sscanf(argv[no_arg-argc],"%lf",&phasejumpangle)){
978 fprintf(stderr,"%s: Error in 1st arg of '--phasejump' option.\n",
979 progname);
980 exit(FAILURE);
981 }
982 --argc;
983 if(!sscanf(argv[no_arg-argc],"%lf",&phasejumpposition)){
984 fprintf(stderr,"%s: Error in 2nd arg of '--phasejump' option.\n",
985 progname);
986 exit(FAILURE);
987 }
988 }
989
990 /*:119*/
991 #line 5976 "./magbragg.w"
992 ;
993 }else if(!strcmp(argv[no_arg-argc],"--fieldevolution")){
994 /*121:*/
995 #line 6205 "./magbragg.w"
996
997 {
998 --argc;
999 if(!strcmp(argv[no_arg-argc],"efield")){
1000 fieldevoflag_efield= 1;
1001 }else if(!strcmp(argv[no_arg-argc],"stoke")){
1002 fieldevoflag_stoke= 1;
1003 }else{
1004 fprintf(stderr,
1005 "%s: Unknown field evolution flag '%s' in second argument of\n"
1006 " --fieldevolution option.\n",progname,argv[no_arg-argc]);
1007 exit(FAILURE);
1008 }
1009 --argc;
1010 if(!sscanf(argv[no_arg-argc],"%ld",&nne)){
1011 fprintf(stderr,"%s: Error in '--fieldevolution' option.\n",
1012 progname);
1013 exit(FAILURE);
1014 }
1015 --argc;
1016 strcpy(fieldevofilename,argv[no_arg-argc]);
1017 fieldevoflag= 1;
1018 if(fieldevoflag_stoke){
1019 sprintf(fieldevofilename_s0,"%s.s0.dat",fieldevofilename);
1020 sprintf(fieldevofilename_s1,"%s.s1.dat",fieldevofilename);
1021 sprintf(fieldevofilename_s2,"%s.s2.dat",fieldevofilename);
1022 sprintf(fieldevofilename_s3,"%s.s3.dat",fieldevofilename);
1023 }
1024 }
1025
1026 /*:121*/
1027 #line 5978 "./magbragg.w"
1028 ;
1029 }else if(!strcmp(argv[no_arg-argc],"--intensityevolution")){
1030 --argc;
1031 if(!sscanf(argv[no_arg-argc],"%lf",&ievolambda)){
1032 fprintf(stderr,"%s: Error in '--intensityevolution' option.\n",
1033 progname);
1034 exit(FAILURE);
1035 }
1036 --argc;
1037 strcpy(intensityevofilename,argv[no_arg-argc]);
1038 intensityevoflag= 1;
1039 }else if(!strcmp(argv[no_arg-argc],"--spectrumfile")){
1040 --argc;
1041 strcpy(spectrumfilename,argv[no_arg-argc]);
1042 }else if(!strcmp(argv[no_arg-argc],"--stokesspectrum")){
1043 stokes_parameter_spectrum= 1;
1044 }else if(!strcmp(argv[no_arg-argc],"--trmtraject")){
1045 --argc;
1046 strcpy(trmtraject_filename,argv[no_arg-argc]);
1047 trmtraject_specified= 1;
1048 }else if(!strcmp(argv[no_arg-argc],"--intensityspectrumfile")){
1049 --argc;
1050 strcpy(intensity_reflection_spectrumfilename,argv[no_arg-argc]);
1051 strcpy(intensity_transmission_spectrumfilename,argv[no_arg-argc]);
1052 strcpy(intensity_check_spectrumfilename,argv[no_arg-argc]);
1053 sprintf(intensity_reflection_spectrumfilename,"%s.irsp.dat",
1054 intensity_reflection_spectrumfilename);
1055 sprintf(intensity_transmission_spectrumfilename,"%s.itsp.dat",
1056 intensity_transmission_spectrumfilename);
1057 sprintf(intensity_check_spectrumfilename,"%s.chec.dat",
1058 intensity_check_spectrumfilename);
1059
1060 }else if(!strcmp(argv[no_arg-argc],"--logarithmicspectrum")){
1061 save_dbspectra= 1;
1062
1063 }else if(!strcmp(argv[no_arg-argc],"-v")||
1064 !strcmp(argv[no_arg-argc],"--verbose")){
1065 verbose= (verbose?0:1);
1066
1067 }else if(!strcmp(argv[no_arg-argc],"--scale_stokesparams")){
1068 scale_stokesparams= 1;
1069 --argc;
1070 if(!sscanf(argv[no_arg-argc],"%lf",&stoke_scalefactor)){
1071 fprintf(stderr,"%s: Error in '--scale_stokesparams' option.\n",\
1072 progname);
1073 exit(FAILURE);
1074 }
1075 }else if(!strcmp(argv[no_arg-argc],"--intensityinfo")){
1076 intensityinfo= (intensityinfo?0:1);
1077 }else if(!strcmp(argv[no_arg-argc],"--intensityinfologfile")){
1078 intensityinfo= 1;
1079 saveintensityinfologfile= 1;
1080 --argc;
1081 strcpy(intensinfologfilename,argv[no_arg-argc]);
1082 }else if(!strcmp(argv[no_arg-argc],"--gyroperturb")){
1083 /*120:*/
1084 #line 6177 "./magbragg.w"
1085
1086 {
1087 perturbed_gyration_constant= 1;
1088 --argc;
1089 if(!sscanf(argv[no_arg-argc],"%lf",&gyroperturb_position)){
1090 fprintf(stderr,"%s: Error in '--gyroperturb' option.\n",progname);
1091 exit(FAILURE);
1092 }
1093 --argc;
1094 if(!sscanf(argv[no_arg-argc],"%lf",&gyroperturb_amplitude)){
1095 fprintf(stderr,"%s: Error in '--gyroperturb' option.\n",progname);
1096 exit(FAILURE);
1097 }
1098 --argc;
1099 if(!sscanf(argv[no_arg-argc],"%lf",&gyroperturb_width)){
1100 fprintf(stderr,"%s: Error in '--gyroperturb' option.\n",progname);
1101 exit(FAILURE);
1102 }
1103 }
1104
1105 /*:120*/
1106 #line 6033 "./magbragg.w"
1107 ;
1108 }else if(!strcmp(argv[no_arg-argc],"--normalize_length_to_um")){
1109 normalize_length_to_micrometer= (normalize_length_to_micrometer?0:1);
1110 }else if(!strcmp(argv[no_arg-argc],"--normalize_intensity")){
1111 normalize_intensity= (normalize_intensity?0:1);
1112 }else if(!strcmp(argv[no_arg-argc],"--normalize_ellipticity")){
1113 normalize_ellipticity= (normalize_ellipticity?0:1);
1114 }else if(!strcmp(argv[no_arg-argc],"--normalizedrepresentation")){
1115 normalize_internally= (normalize_internally?0:1);
1116 }else if(!strcmp(argv[no_arg-argc],"-r")||
1117 !strcmp(argv[no_arg-argc],"--random")){
1118 randomdistribution= (randomdistribution?0:1);
1119 }else if(!strcmp(argv[no_arg-argc],"-h")||
1120 !strcmp(argv[no_arg-argc],"--help")){
1121 showsomehelp();
1122 }else if(!strcmp(argv[no_arg-argc],"-g")||
1123 !strcmp(argv[no_arg-argc],"--grating")){
1124 --argc;
1125 if(!strcmp(argv[no_arg-argc],"stepwise")){
1126 /*122:*/
1127 #line 6238 "./magbragg.w"
1128
1129 {
1130 strcpy(gratingtype,argv[no_arg-argc]);
1131 --argc;
1132 if(!strcmp(argv[no_arg-argc],"twolevel")){
1133 strcpy(gratingsubtype,argv[no_arg-argc]);
1134 --argc;
1135 if(strcmp(argv[no_arg-argc],"t1")){
1136 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1137 fprintf(stderr,"%s: (Expecting string 't1').\n",progname);
1138 exit(FAILURE);
1139 }
1140 --argc;
1141 if(!sscanf(argv[no_arg-argc],"%lf",&t1)){
1142 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1143 fprintf(stderr,"%s: (Could not read data for t1).\n",progname);
1144 exit(FAILURE);
1145 }
1146 --argc;
1147 if(strcmp(argv[no_arg-argc],"t2")){
1148 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1149 fprintf(stderr,"%s: (Expecting string 't2').\n",progname);
1150 exit(FAILURE);
1151 }
1152 --argc;
1153 if(!sscanf(argv[no_arg-argc],"%lf",&t2)){
1154 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1155 fprintf(stderr,"%s: (Could not read data for t2).\n",progname);
1156 exit(FAILURE);
1157 }
1158 --argc;
1159 if(strcmp(argv[no_arg-argc],"n1")){
1160 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1161 fprintf(stderr,"%s: (Expecting string 'n1').\n",progname);
1162 exit(FAILURE);
1163 }
1164 --argc;
1165 if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
1166 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1167 fprintf(stderr,"%s: (Could not read data for n1).\n",progname);
1168 exit(FAILURE);
1169 }
1170 --argc;
1171 if(strcmp(argv[no_arg-argc],"n2")){
1172 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1173 fprintf(stderr,"%s: (Expecting string 'n2').\n",progname);
1174 exit(FAILURE);
1175 }
1176 --argc;
1177 if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
1178 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1179 fprintf(stderr,"%s: (Could not read data for n2).\n",progname);
1180 exit(FAILURE);
1181 }
1182 --argc;
1183 if(strcmp(argv[no_arg-argc],"g1")){
1184 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1185 fprintf(stderr,"%s: (Expecting string 'g1').\n",progname);
1186 exit(FAILURE);
1187 }
1188 --argc;
1189 if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
1190 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1191 fprintf(stderr,"%s: (Could not read data for g1).\n",progname);
1192 exit(FAILURE);
1193 }
1194 --argc;
1195 if(strcmp(argv[no_arg-argc],"g2")){
1196 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1197 fprintf(stderr,"%s: (Expecting string 'g2').\n",progname);
1198 exit(FAILURE);
1199 }
1200 --argc;
1201 if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
1202 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1203 fprintf(stderr,"%s: (Could not read data for g2).\n",progname);
1204 exit(FAILURE);
1205 }
1206 --argc;
1207 if(strcmp(argv[no_arg-argc],"pe1")){
1208 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1209 fprintf(stderr,"%s: (Expecting string 'pe1').\n",progname);
1210 exit(FAILURE);
1211 }
1212 --argc;
1213 if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
1214 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1215 fprintf(stderr,"%s: (Could not read data for pe1).\n",progname);
1216 exit(FAILURE);
1217 }
1218 --argc;
1219 if(strcmp(argv[no_arg-argc],"pe2")){
1220 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1221 fprintf(stderr,"%s: (Expecting string 'pe2').\n",progname);
1222 exit(FAILURE);
1223 }
1224 --argc;
1225 if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
1226 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1227 fprintf(stderr,"%s: (Could not read data for pe2).\n",progname);
1228 exit(FAILURE);
1229 }
1230 --argc;
1231 if(strcmp(argv[no_arg-argc],"pm1")){
1232 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1233 fprintf(stderr,"%s: (Expecting string 'pm1').\n",progname);
1234 exit(FAILURE);
1235 }
1236 --argc;
1237 if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
1238 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1239 fprintf(stderr,"%s: (Could not read data for pm1).\n",\
1240 progname);
1241 exit(FAILURE);
1242 }
1243 --argc;
1244 if(strcmp(argv[no_arg-argc],"pm2")){
1245 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1246 fprintf(stderr,"%s: (Expecting string 'pm2').\n",progname);
1247 exit(FAILURE);
1248 }
1249 --argc;
1250 if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
1251 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1252 fprintf(stderr,"%s: (Could not read data for pm2).\n",\
1253 progname);
1254 exit(FAILURE);
1255 }
1256 --argc;
1257 if(strcmp(argv[no_arg-argc],"qe1")){
1258 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1259 fprintf(stderr,"%s: (Expecting string 'qe1').\n",progname);
1260 exit(FAILURE);
1261 }
1262 --argc;
1263 if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
1264 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1265 fprintf(stderr,"%s: (Could not read data for qe1).\n",\
1266 progname);
1267 exit(FAILURE);
1268 }
1269 --argc;
1270 if(strcmp(argv[no_arg-argc],"qe2")){
1271 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1272 fprintf(stderr,"%s: (Expecting string 'qe2').\n",progname);
1273 exit(FAILURE);
1274 }
1275 --argc;
1276 if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
1277 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1278 fprintf(stderr,"%s: (Could not read data for qe2).\n",progname);
1279 exit(FAILURE);
1280 }
1281 --argc;
1282 if(strcmp(argv[no_arg-argc],"qm1")){
1283 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1284 fprintf(stderr,"%s: (Expecting string 'qm1').\n",progname);
1285 exit(FAILURE);
1286 }
1287 --argc;
1288 if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
1289 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1290 fprintf(stderr,"%s: (Could not read data for qm1).\n",progname);
1291 exit(FAILURE);
1292 }
1293 --argc;
1294 if(strcmp(argv[no_arg-argc],"qm2")){
1295 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1296 fprintf(stderr,"%s: (Expecting string 'qm2').\n",progname);
1297 exit(FAILURE);
1298 }
1299 --argc;
1300 if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
1301 fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
1302 fprintf(stderr,"%s: (Could not read data for qm2).\n",progname);
1303 exit(FAILURE);
1304 }
1305 }else{
1306 fprintf(stderr,"%s: Error in '-g' or '--grating' option.\n",progname);
1307 fprintf(stderr,"%s: (No valid stepwise grating type found!)\n",progname);
1308 exit(FAILURE);
1309 }
1310 }
1311
1312 /*:122*/
1313 #line 6052 "./magbragg.w"
1314
1315 }else if(!strcmp(argv[no_arg-argc],"sinusoidal")){
1316 /*123:*/
1317 #line 6425 "./magbragg.w"
1318
1319 {
1320 strcpy(gratingtype,argv[no_arg-argc]);
1321 --argc;
1322 if(strcmp(argv[no_arg-argc],"n")){
1323 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1324 fprintf(stderr,"%s: (Expecting string 'n').\n",progname);
1325 exit(FAILURE);
1326 }
1327 --argc;
1328 if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
1329 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1330 fprintf(stderr,"%s: (Could not read data for n [1st arg]).\n",progname);
1331 exit(FAILURE);
1332 }
1333 --argc;
1334 if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
1335 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1336 fprintf(stderr,"%s: (Could not read data for n [2nd arg]).\n",progname);
1337 exit(FAILURE);
1338 }
1339 --argc;
1340 if(!sscanf(argv[no_arg-argc],"%lf",&nper)){
1341 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1342 fprintf(stderr,"%s: (Could not read data for n [3rd arg]).\n",progname);
1343 exit(FAILURE);
1344 }
1345 --argc;
1346 if(strcmp(argv[no_arg-argc],"g")){
1347 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1348 fprintf(stderr,"%s: (Expecting string 'g').\n",progname);
1349 exit(FAILURE);
1350 }
1351 --argc;
1352 if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
1353 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1354 fprintf(stderr,"%s: (Could not read data for g [1st arg]).\n",progname);
1355 exit(FAILURE);
1356 }
1357 --argc;
1358 if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
1359 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1360 fprintf(stderr,"%s: (Could not read data for g [2nd arg]).\n",progname);
1361 exit(FAILURE);
1362 }
1363 --argc;
1364 if(!sscanf(argv[no_arg-argc],"%lf",&gper)){
1365 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1366 fprintf(stderr,"%s: (Could not read data for g [3rd arg]).\n",progname);
1367 exit(FAILURE);
1368 }
1369 --argc;
1370 if(strcmp(argv[no_arg-argc],"pe")){
1371 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1372 fprintf(stderr,"%s: (Expecting string 'pe').\n",progname);
1373 exit(FAILURE);
1374 }
1375 --argc;
1376 if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
1377 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1378 fprintf(stderr,"%s: (Could not read data for pe [1st arg]).\n",progname);
1379 exit(FAILURE);
1380 }
1381 --argc;
1382 if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
1383 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1384 fprintf(stderr,"%s: (Could not read data for pe [2nd arg]).\n",progname);
1385 exit(FAILURE);
1386 }
1387 --argc;
1388 if(!sscanf(argv[no_arg-argc],"%lf",&peper)){
1389 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1390 fprintf(stderr,"%s: (Could not read data for pe [3rd arg]).\n",progname);
1391 exit(FAILURE);
1392 }
1393 --argc;
1394 if(strcmp(argv[no_arg-argc],"pm")){
1395 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1396 fprintf(stderr,"%s: (Expecting string 'pm').\n",progname);
1397 exit(FAILURE);
1398 }
1399 --argc;
1400 if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
1401 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1402 fprintf(stderr,"%s: (Could not read data for pm [1st arg]).\n",progname);
1403 exit(FAILURE);
1404 }
1405 --argc;
1406 if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
1407 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1408 fprintf(stderr,"%s: (Could not read data for pm [2nd arg]).\n",progname);
1409 exit(FAILURE);
1410 }
1411 --argc;
1412 if(!sscanf(argv[no_arg-argc],"%lf",&pmper)){
1413 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1414 fprintf(stderr,"%s: (Could not read data for pm [3rd arg]).\n",progname);
1415 exit(FAILURE);
1416 }
1417 --argc;
1418 if(strcmp(argv[no_arg-argc],"qe")){
1419 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1420 fprintf(stderr,"%s: (Expecting string 'qe').\n",progname);
1421 exit(FAILURE);
1422 }
1423 --argc;
1424 if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
1425 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1426 fprintf(stderr,"%s: (Could not read data for qe [1st arg]).\n",progname);
1427 exit(FAILURE);
1428 }
1429 --argc;
1430 if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
1431 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1432 fprintf(stderr,"%s: (Could not read data for qe [2nd arg]).\n",progname);
1433 exit(FAILURE);
1434 }
1435 --argc;
1436 if(!sscanf(argv[no_arg-argc],"%lf",&qeper)){
1437 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1438 fprintf(stderr,"%s: (Could not read data for qe [3rd arg]).\n",progname);
1439 exit(FAILURE);
1440 }
1441 --argc;
1442 if(strcmp(argv[no_arg-argc],"qm")){
1443 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1444 fprintf(stderr,"%s: (Expecting string 'qm').\n",progname);
1445 exit(FAILURE);
1446 }
1447 --argc;
1448 if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
1449 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1450 fprintf(stderr,"%s: (Could not read data for qm [1st arg]).\n",progname);
1451 exit(FAILURE);
1452 }
1453 --argc;
1454 if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
1455 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1456 fprintf(stderr,"%s: (Could not read data for qm [2nd arg]).\n",progname);
1457 exit(FAILURE);
1458 }
1459 --argc;
1460 if(!sscanf(argv[no_arg-argc],"%lf",&qmper)){
1461 fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
1462 fprintf(stderr,"%s: (Could not read data for qm [3rd arg]).\n",progname);
1463 exit(FAILURE);
1464 }
1465 }
1466
1467 /*:123*/
1468 #line 6054 "./magbragg.w"
1469
1470 }else if(!strcmp(argv[no_arg-argc],"chirped")){
1471 /*124:*/
1472 #line 6577 "./magbragg.w"
1473
1474 {
1475 strcpy(gratingtype,argv[no_arg-argc]);
1476 --argc;
1477 if(strcmp(argv[no_arg-argc],"n")){
1478 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1479 fprintf(stderr,"%s: (Expecting string 'n').\n",progname);
1480 exit(FAILURE);
1481 }
1482 --argc;
1483 if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
1484 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1485 fprintf(stderr,"%s: (Could not read data for n [1st arg]).\n",progname);
1486 exit(FAILURE);
1487 }
1488 --argc;
1489 if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
1490 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1491 fprintf(stderr,"%s: (Could not read data for n [2nd arg]).\n",progname);
1492 exit(FAILURE);
1493 }
1494 --argc;
1495 if(!sscanf(argv[no_arg-argc],"%lf",&nper)){
1496 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1497 fprintf(stderr,"%s: (Could not read data for n [3rd arg]).\n",progname);
1498 exit(FAILURE);
1499 }
1500 --argc;
1501 if(!sscanf(argv[no_arg-argc],"%lf",&ncrp)){
1502 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1503 fprintf(stderr,"%s: (Could not read data for n [4th arg]).\n",progname);
1504 exit(FAILURE);
1505 }
1506 --argc;
1507 if(strcmp(argv[no_arg-argc],"g")){
1508 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1509 fprintf(stderr,"%s: (Expecting string 'g').\n",progname);
1510 exit(FAILURE);
1511 }
1512 --argc;
1513 if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
1514 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1515 fprintf(stderr,"%s: (Could not read data for g [1st arg]).\n",progname);
1516 exit(FAILURE);
1517 }
1518 --argc;
1519 if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
1520 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1521 fprintf(stderr,"%s: (Could not read data for g [2nd arg]).\n",progname);
1522 exit(FAILURE);
1523 }
1524 --argc;
1525 if(!sscanf(argv[no_arg-argc],"%lf",&gper)){
1526 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1527 fprintf(stderr,"%s: (Could not read data for g [3rd arg]).\n",progname);
1528 exit(FAILURE);
1529 }
1530 --argc;
1531 if(!sscanf(argv[no_arg-argc],"%lf",&gcrp)){
1532 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1533 fprintf(stderr,"%s: (Could not read data for g [4th arg]).\n",progname);
1534 exit(FAILURE);
1535 }
1536 --argc;
1537 if(strcmp(argv[no_arg-argc],"pe")){
1538 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1539 fprintf(stderr,"%s: (Expecting string 'pe').\n",progname);
1540 exit(FAILURE);
1541 }
1542 --argc;
1543 if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
1544 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1545 fprintf(stderr,"%s: (Could not read data for pe [1st arg]).\n",progname);
1546 exit(FAILURE);
1547 }
1548 --argc;
1549 if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
1550 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1551 fprintf(stderr,"%s: (Could not read data for pe [2nd arg]).\n",progname);
1552 exit(FAILURE);
1553 }
1554 --argc;
1555 if(!sscanf(argv[no_arg-argc],"%lf",&peper)){
1556 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1557 fprintf(stderr,"%s: (Could not read data for pe [3rd arg]).\n",progname);
1558 exit(FAILURE);
1559 }
1560 --argc;
1561 if(!sscanf(argv[no_arg-argc],"%lf",&pecrp)){
1562 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1563 fprintf(stderr,"%s: (Could not read data for pe [4th arg]).\n",progname);
1564 exit(FAILURE);
1565 }
1566 --argc;
1567 if(strcmp(argv[no_arg-argc],"pm")){
1568 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1569 fprintf(stderr,"%s: (Expecting string 'pm').\n",progname);
1570 exit(FAILURE);
1571 }
1572 --argc;
1573 if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
1574 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1575 fprintf(stderr,"%s: (Could not read data for pm [1st arg]).\n",progname);
1576 exit(FAILURE);
1577 }
1578 --argc;
1579 if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
1580 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1581 fprintf(stderr,"%s: (Could not read data for pm [2nd arg]).\n",progname);
1582 exit(FAILURE);
1583 }
1584 --argc;
1585 if(!sscanf(argv[no_arg-argc],"%lf",&pmper)){
1586 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1587 fprintf(stderr,"%s: (Could not read data for pm [3rd arg]).\n",progname);
1588 exit(FAILURE);
1589 }
1590 --argc;
1591 if(!sscanf(argv[no_arg-argc],"%lf",&pmcrp)){
1592 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1593 fprintf(stderr,"%s: (Could not read data for pm [4th arg]).\n",progname);
1594 exit(FAILURE);
1595 }
1596 --argc;
1597 if(strcmp(argv[no_arg-argc],"qe")){
1598 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1599 fprintf(stderr,"%s: (Expecting string 'qe').\n",progname);
1600 exit(FAILURE);
1601 }
1602 --argc;
1603 if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
1604 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1605 fprintf(stderr,"%s: (Could not read data for qe [1st arg]).\n",progname);
1606 exit(FAILURE);
1607 }
1608 --argc;
1609 if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
1610 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1611 fprintf(stderr,"%s: (Could not read data for qe [2nd arg]).\n",progname);
1612 exit(FAILURE);
1613 }
1614 --argc;
1615 if(!sscanf(argv[no_arg-argc],"%lf",&qeper)){
1616 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1617 fprintf(stderr,"%s: (Could not read data for qe [3rd arg]).\n",progname);
1618 exit(FAILURE);
1619 }
1620 --argc;
1621 if(!sscanf(argv[no_arg-argc],"%lf",&qecrp)){
1622 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1623 fprintf(stderr,"%s: (Could not read data for qe [4th arg]).\n",progname);
1624 exit(FAILURE);
1625 }
1626 --argc;
1627 if(strcmp(argv[no_arg-argc],"qm")){
1628 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1629 fprintf(stderr,"%s: (Expecting string 'qm').\n",progname);
1630 exit(FAILURE);
1631 }
1632 --argc;
1633 if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
1634 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1635 fprintf(stderr,"%s: (Could not read data for qm [1st arg]).\n",progname);
1636 exit(FAILURE);
1637 }
1638 --argc;
1639 if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
1640 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1641 fprintf(stderr,"%s: (Could not read data for qm [2nd arg]).\n",progname);
1642 exit(FAILURE);
1643 }
1644 --argc;
1645 if(!sscanf(argv[no_arg-argc],"%lf",&qmper)){
1646 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1647 fprintf(stderr,"%s: (Could not read data for qm [3rd arg]).\n",progname);
1648 exit(FAILURE);
1649 }
1650 --argc;
1651 if(!sscanf(argv[no_arg-argc],"%lf",&qmcrp)){
1652 fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
1653 fprintf(stderr,"%s: (Could not read data for qm [4th arg]).\n",progname);
1654 exit(FAILURE);
1655 }
1656 }
1657
1658 /*:124*/
1659 #line 6056 "./magbragg.w"
1660
1661 }else if(!strcmp(argv[no_arg-argc],"fractal")){
1662 /*125:*/
1663 #line 6766 "./magbragg.w"
1664
1665 {
1666 strcpy(gratingtype,argv[no_arg-argc]);
1667 --argc;
1668 if(!strcmp(argv[no_arg-argc],"cantor")){
1669 strcpy(gratingsubtype,argv[no_arg-argc]);
1670 --argc;
1671 if(strcmp(argv[no_arg-argc],"fractal_level")){
1672 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1673 fprintf(stderr,"%s: (Expecting string 'fractal_level').\n",progname);
1674 exit(FAILURE);
1675 }
1676 --argc;
1677 if(!sscanf(argv[no_arg-argc],"%d",&fractal_level)){
1678 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1679 fprintf(stderr,"%s: (Could not read data for fractal_level).\n",
1680 progname);
1681 exit(FAILURE);
1682 }
1683 --argc;
1684 if(strcmp(argv[no_arg-argc],"maximum_fractal_level")){
1685 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1686 fprintf(stderr,"%s: (Expecting string 'maximum_fractal_level').\n",
1687 progname);
1688 exit(FAILURE);
1689 }
1690 --argc;
1691 if(!sscanf(argv[no_arg-argc],"%d",&maximum_fractal_level)){
1692 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1693 fprintf(stderr,"%s: (Could not read data for maximum_fractal_level).\n",
1694 progname);
1695 exit(FAILURE);
1696 }
1697 nn= 1;
1698 for(j= 1;j<=fractal_level;j++)nn= 2*nn;
1699 --argc;
1700 if(strcmp(argv[no_arg-argc],"t1")){
1701 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1702 fprintf(stderr,"%s: (Expecting string 't1').\n",progname);
1703 exit(FAILURE);
1704 }
1705 --argc;
1706 if(!sscanf(argv[no_arg-argc],"%lf",&t1)){
1707 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1708 fprintf(stderr,"%s: (Could not read data for t1).\n",progname);
1709 exit(FAILURE);
1710 }
1711 --argc;
1712 if(strcmp(argv[no_arg-argc],"t2")){
1713 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1714 fprintf(stderr,"%s: (Expecting string 't2').\n",progname);
1715 exit(FAILURE);
1716 }
1717 --argc;
1718 if(!sscanf(argv[no_arg-argc],"%lf",&t2)){
1719 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1720 fprintf(stderr,"%s: (Could not read data for t2).\n",progname);
1721 exit(FAILURE);
1722 }
1723 --argc;
1724 if(strcmp(argv[no_arg-argc],"n1")){
1725 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1726 fprintf(stderr,"%s: (Expecting string 'n1').\n",progname);
1727 exit(FAILURE);
1728 }
1729 --argc;
1730 if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
1731 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1732 fprintf(stderr,"%s: (Could not read data for n1).\n",progname);
1733 exit(FAILURE);
1734 }
1735 --argc;
1736 if(strcmp(argv[no_arg-argc],"n2")){
1737 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1738 fprintf(stderr,"%s: (Expecting string 'n2').\n",progname);
1739 exit(FAILURE);
1740 }
1741 --argc;
1742 if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
1743 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1744 fprintf(stderr,"%s: (Could not read data for n2).\n",progname);
1745 exit(FAILURE);
1746 }
1747 --argc;
1748 if(strcmp(argv[no_arg-argc],"g1")){
1749 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1750 fprintf(stderr,"%s: (Expecting string 'g1').\n",progname);
1751 exit(FAILURE);
1752 }
1753 --argc;
1754 if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
1755 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1756 fprintf(stderr,"%s: (Could not read data for g1).\n",progname);
1757 exit(FAILURE);
1758 }
1759 --argc;
1760 if(strcmp(argv[no_arg-argc],"g2")){
1761 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1762 fprintf(stderr,"%s: (Expecting string 'g2').\n",progname);
1763 exit(FAILURE);
1764 }
1765 --argc;
1766 if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
1767 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1768 fprintf(stderr,"%s: (Could not read data for g2).\n",progname);
1769 exit(FAILURE);
1770 }
1771 --argc;
1772 if(strcmp(argv[no_arg-argc],"pe1")){
1773 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1774 fprintf(stderr,"%s: (Expecting string 'pe1').\n",progname);
1775 exit(FAILURE);
1776 }
1777 --argc;
1778 if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
1779 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1780 fprintf(stderr,"%s: (Could not read data for pe1).\n",progname);
1781 exit(FAILURE);
1782 }
1783 --argc;
1784 if(strcmp(argv[no_arg-argc],"pe2")){
1785 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1786 fprintf(stderr,"%s: (Expecting string 'pe2').\n",progname);
1787 exit(FAILURE);
1788 }
1789 --argc;
1790 if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
1791 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1792 fprintf(stderr,"%s: (Could not read data for pe2).\n",progname);
1793 exit(FAILURE);
1794 }
1795 --argc;
1796 if(strcmp(argv[no_arg-argc],"pm1")){
1797 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1798 fprintf(stderr,"%s: (Expecting string 'pm1').\n",progname);
1799 exit(FAILURE);
1800 }
1801 --argc;
1802 if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
1803 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1804 fprintf(stderr,"%s: (Could not read data for pm1).\n",\
1805 progname);
1806 exit(FAILURE);
1807 }
1808 --argc;
1809 if(strcmp(argv[no_arg-argc],"pm2")){
1810 fprintf(stderr,"%s: Error.\n",progname);
1811 fprintf(stderr,"%s: (Expecting string 'pm2').\n",progname);
1812 exit(FAILURE);
1813 }
1814 --argc;
1815 if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
1816 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1817 fprintf(stderr,"%s: (Could not read data for pm2).\n",\
1818 progname);
1819 exit(FAILURE);
1820 }
1821 --argc;
1822 if(strcmp(argv[no_arg-argc],"qe1")){
1823 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1824 fprintf(stderr,"%s: (Expecting string 'qe1').\n",progname);
1825 exit(FAILURE);
1826 }
1827 --argc;
1828 if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
1829 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1830 fprintf(stderr,"%s: (Could not read data for qe1).\n",\
1831 progname);
1832 exit(FAILURE);
1833 }
1834 --argc;
1835 if(strcmp(argv[no_arg-argc],"qe2")){
1836 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1837 fprintf(stderr,"%s: (Expecting string 'qe2').\n",progname);
1838 exit(FAILURE);
1839 }
1840 --argc;
1841 if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
1842 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1843 fprintf(stderr,"%s: (Could not read data for qe2).\n",progname);
1844 exit(FAILURE);
1845 }
1846 --argc;
1847 if(strcmp(argv[no_arg-argc],"qm1")){
1848 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1849 fprintf(stderr,"%s: (Expecting string 'qm1').\n",progname);
1850 exit(FAILURE);
1851 }
1852 --argc;
1853 if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
1854 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1855 fprintf(stderr,"%s: (Could not read data for qm1).\n",progname);
1856 exit(FAILURE);
1857 }
1858 --argc;
1859 if(strcmp(argv[no_arg-argc],"qm2")){
1860 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1861 fprintf(stderr,"%s: (Expecting string 'qm2').\n",progname);
1862 exit(FAILURE);
1863 }
1864 --argc;
1865 if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
1866 fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
1867 fprintf(stderr,"%s: (Could not read data for qm2).\n",progname);
1868 exit(FAILURE);
1869 }
1870 }else{
1871 fprintf(stderr,"%s: Error in 'fractal' grating option.\n",progname);
1872 fprintf(stderr,"%s: (No valid fractal type found!)\n",progname);
1873 fprintf(stderr,"%s: (Currently only Cantor type implemented)\n",progname);
1874 exit(FAILURE);
1875 }
1876 }
1877
1878 /*:125*/
1879 #line 6058 "./magbragg.w"
1880
1881 }else{
1882 fprintf(stderr,"%s: Error in '-g' or '--grating' option.\n",
1883 progname);
1884 fprintf(stderr,"%s: (No valid grating type found!)\n",progname);
1885 exit(FAILURE);
1886 }
1887 }else if((!strcmp(argv[no_arg-argc],"-L"))||
1888 (!strcmp(argv[no_arg-argc],"--gratinglength"))){
1889 --argc;
1890 if(!sscanf(argv[no_arg-argc],"%lf",&ll)){
1891 fprintf(stderr,"%s: Error in '-L' option.\n",progname);
1892 exit(FAILURE);
1893 }
1894 }else if(!strcmp(argv[no_arg-argc],"--modifylayer")){
1895 /*126:*/
1896 #line 6984 "./magbragg.w"
1897
1898 {
1899 --argc;
1900 if(strcmp(argv[no_arg-argc],"num")){
1901 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1902 fprintf(stderr,"%s: (Expecting string 'num' after '--modifylayer').\n",
1903 progname);
1904 exit(FAILURE);
1905 }
1906 --argc;
1907 if(!sscanf(argv[no_arg-argc],"%ld",&modnum)){
1908 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1909 fprintf(stderr,"%s: (Could not read data for num).\n",progname);
1910 exit(FAILURE);
1911 }
1912 --argc;
1913 if(strcmp(argv[no_arg-argc],"t")){
1914 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1915 fprintf(stderr,"%s: (Expecting string 't' [for layer thickness]).\n",
1916 progname);
1917 exit(FAILURE);
1918 }
1919 --argc;
1920 if(!sscanf(argv[no_arg-argc],"%lf",&modt1)){
1921 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1922 fprintf(stderr,"%s: (Could not read data for t [layer thickness]).\n",
1923 progname);
1924 exit(FAILURE);
1925 }
1926 --argc;
1927 if(strcmp(argv[no_arg-argc],"n")){
1928 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1929 fprintf(stderr,"%s: (Expecting string 'n' [for refractive index]).\n",
1930 progname);
1931 exit(FAILURE);
1932 }
1933 --argc;
1934 if(!sscanf(argv[no_arg-argc],"%lf",&modn1)){
1935 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1936 fprintf(stderr,"%s: (Could not read data for n [refractive index]).\n",
1937 progname);
1938 exit(FAILURE);
1939 }
1940 --argc;
1941 if(strcmp(argv[no_arg-argc],"g")){
1942 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1943 fprintf(stderr,"%s: (Expecting string 'g' [for gyration constant]).\n",
1944 progname);
1945 exit(FAILURE);
1946 }
1947 --argc;
1948 if(!sscanf(argv[no_arg-argc],"%lf",&modg1)){
1949 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1950 fprintf(stderr,"%s: (Could not read data for g [gyration constant]).\n",
1951 progname);
1952 exit(FAILURE);
1953 }
1954 --argc;
1955 if(strcmp(argv[no_arg-argc],"pe")){
1956 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1957 fprintf(stderr,"%s: (Expecting string 'pe').\n",progname);
1958 exit(FAILURE);
1959 }
1960 --argc;
1961 if(!sscanf(argv[no_arg-argc],"%lf",&modpe1)){
1962 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1963 fprintf(stderr,"%s: (Could not read data for pe).\n",progname);
1964 exit(FAILURE);
1965 }
1966 --argc;
1967 if(strcmp(argv[no_arg-argc],"pm")){
1968 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1969 fprintf(stderr,"%s: (Expecting string 'pm').\n",progname);
1970 exit(FAILURE);
1971 }
1972 --argc;
1973 if(!sscanf(argv[no_arg-argc],"%lf",&modpm1)){
1974 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1975 fprintf(stderr,"%s: (Could not read data for pm).\n",progname);
1976 exit(FAILURE);
1977 }
1978 --argc;
1979 if(strcmp(argv[no_arg-argc],"qe")){
1980 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1981 fprintf(stderr,"%s: (Expecting string 'qe').\n",progname);
1982 exit(FAILURE);
1983 }
1984 --argc;
1985 if(!sscanf(argv[no_arg-argc],"%lf",&modqe1)){
1986 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1987 fprintf(stderr,"%s: (Could not read data for qe).\n",progname);
1988 exit(FAILURE);
1989 }
1990 --argc;
1991 if(strcmp(argv[no_arg-argc],"qm")){
1992 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1993 fprintf(stderr,"%s: (Expecting string 'qm').\n",progname);
1994 exit(FAILURE);
1995 }
1996 --argc;
1997 if(!sscanf(argv[no_arg-argc],"%lf",&modqm1)){
1998 fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
1999 fprintf(stderr,"%s: (Could not read data for qm).\n",progname);
2000 exit(FAILURE);
2001 }
2002 }
2003
2004 /*:126*/
2005 #line 6073 "./magbragg.w"
2006
2007 }else if(!strcmp(argv[no_arg-argc],"--refindsurr")){
2008 --argc;
2009 if(!sscanf(argv[no_arg-argc],"%lf",&nsurr)){
2010 fprintf(stderr,"%s: Error in '--refindsurr' option.\n",progname);
2011 exit(FAILURE);
2012 }
2013 }else if(!strcmp(argv[no_arg-argc],"--trmintensity")){
2014 /*127:*/
2015 #line 7103 "./magbragg.w"
2016
2017 {
2018 --argc;
2019 if(!sscanf(argv[no_arg-argc],"%lf",&trmintenstart)){
2020 fprintf(stderr,"%s: Error in '--trmintensity' option.\n",\
2021 progname);
2022 exit(FAILURE);
2023 }
2024 --argc;
2025 if(!sscanf(argv[no_arg-argc],"%lf",&trmintenstop)){
2026 fprintf(stderr,"%s: Error in '--trmintensity' option.\n",\
2027 progname);
2028 exit(FAILURE);
2029 }
2030 --argc;
2031 if(!sscanf(argv[no_arg-argc],"%ld",&mmi)){
2032 fprintf(stderr,"%s: Error in '--trmintensity' option.\n",\
2033 progname);
2034 exit(FAILURE);
2035 }
2036 }
2037
2038 /*:127*/
2039 #line 6081 "./magbragg.w"
2040
2041 }else if(!strcmp(argv[no_arg-argc],"--trmellipticity")){
2042 /*128:*/
2043 #line 7137 "./magbragg.w"
2044
2045 {
2046 --argc;
2047 if(!sscanf(argv[no_arg-argc],"%lf",&trmellipstart)){
2048 fprintf(stderr,"%s: Error in '--trmellipticity' option.\n",progname);
2049 exit(FAILURE);
2050 }
2051 --argc;
2052 if(!sscanf(argv[no_arg-argc],"%lf",&trmellipstop)){
2053 fprintf(stderr,"%s: Error in '--trmellipticity' option.\n",progname);
2054 exit(FAILURE);
2055 }
2056 --argc;
2057 if(!sscanf(argv[no_arg-argc],"%ld",&mme)){
2058 fprintf(stderr,"%s: Error in '--trmellipticity' option.\n",progname);
2059 exit(FAILURE);
2060 }
2061 }
2062
2063 /*:128*/
2064 #line 6083 "./magbragg.w"
2065
2066 }else if(!strcmp(argv[no_arg-argc],"--lambdastart")){
2067 --argc;
2068 if(!sscanf(argv[no_arg-argc],"%lf",&lambdastart)){
2069 fprintf(stderr,"%s: Error in '--lambdastart' option.\n",progname);
2070 exit(FAILURE);
2071 }
2072 }else if(!strcmp(argv[no_arg-argc],"--lambdastop")){
2073 --argc;
2074 if(!sscanf(argv[no_arg-argc],"%lf",&lambdastop)){
2075 fprintf(stderr,"%s: Error in '--lambdastop' option.\n",progname);
2076 exit(FAILURE);
2077 }
2078 }else if(!strcmp(argv[no_arg-argc],"--wlenlinz")){
2079 chirpflag= 1;
2080 }else if(!strcmp(argv[no_arg-argc],"--freqlinz")){
2081 chirpflag= 0;
2082 }else if(!strcmp(argv[no_arg-argc],"-N")){
2083 --argc;
2084 if(!sscanf(argv[no_arg-argc],"%ld",&nn)){
2085 fprintf(stderr,"%s: Error in '-N' option.\n",progname);
2086 exit(FAILURE);
2087 }
2088 }else if(!strcmp(argv[no_arg-argc],"-M")){
2089 --argc;
2090 if(!sscanf(argv[no_arg-argc],"%ld",&mm)){
2091 fprintf(stderr,"%s: Error in '-M' option.\n",progname);
2092 exit(FAILURE);
2093 }
2094 }else{
2095 fprintf(stderr,"%s: Specified option invalid!\n",progname);
2096 showsomehelp();
2097 exit(FAILURE);
2098 }
2099 }
2100 /*117:*/
2101 #line 6125 "./magbragg.w"
2102
2103 {
2104 sprintf(outfilename_s0,"%s.s0.dat",outfilename);
2105 sprintf(outfilename_s1,"%s.s1.dat",outfilename);
2106 sprintf(outfilename_s2,"%s.s2.dat",outfilename);
2107 sprintf(outfilename_s3,"%s.s3.dat",outfilename);
2108 sprintf(outfilename_v0,"%s.v0.dat",outfilename);
2109 sprintf(outfilename_v1,"%s.v1.dat",outfilename);
2110 sprintf(outfilename_v2,"%s.v2.dat",outfilename);
2111 sprintf(outfilename_v3,"%s.v3.dat",outfilename);
2112 sprintf(outfilename_w0,"%s.w0.dat",outfilename);
2113 sprintf(outfilename_w1,"%s.w1.dat",outfilename);
2114 sprintf(outfilename_w2,"%s.w2.dat",outfilename);
2115 sprintf(outfilename_w3,"%s.w3.dat",outfilename);
2116 }
2117
2118 /*:117*/
2119 #line 6118 "./magbragg.w"
2120
2121 /*129:*/
2122 #line 7161 "./magbragg.w"
2123
2124 {
2125 if(verbose){
2126 for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
2127 fprintf(stdout,"Input parameters:\n");
2128 fprintf(stdout,
2129 "Grating type: %s\n",gratingtype);
2130 if(!strcmp(gratingtype,"sinusoidal")){
2131 fprintf(stdout,
2132 "Bias refractive index: %10.6e\n",n1);
2133 fprintf(stdout,
2134 "Refractive index modulation: %10.6e\n",n2);
2135 fprintf(stdout,
2136 "Modulation period: %10.6e [m]\n",nper);
2137 }else if(!strcmp(gratingtype,"chirped")){
2138 }else if(!strcmp(gratingtype,"stepwise")){
2139 }else if(!strcmp(gratingtype,"fractal")){
2140 }else{
2141 fprintf(stdout,"%s: Error: Unknown grating type '%s'\n",
2142 progname,gratingtype);
2143 }
2144 fprintf(stdout,"Geometrical length: "
2145 "L=%10.6e [m]\n",ll);
2146 fprintf(stdout,"Surrounding refractive index: "
2147 "nsurr=%10.6e\n",nsurr);
2148 fprintf(stdout,"Begin wavelength of spectrum: "
2149 "lambda_start=%10.6e [m]\n",lambdastart);
2150 fprintf(stdout,"End wavelength of spectrum: "
2151 "lambda_stop=%10.6e [m]\n",lambdastop);
2152 fprintf(stdout,"Number of samples in spectrum: "
2153 "M=%-12ld\n",mm);
2154 fprintf(stdout,"Number of discrete layers: "
2155 "N=%-12ld\n",nn);
2156 if(trmtraject_specified){
2157 fprintf(stdout,"Trajectory for transmitted Stokes parameters: %s\n",
2158 trmtraject_filename);
2159 }else{
2160 fprintf(stdout,"Number of samples in output intensity: "
2161 "mmi=%-12ld\n",mmi);
2162 fprintf(stdout,"Number of samples in output ellipticity: "
2163 "mme=%-12ld\n",mme);
2164 }
2165 fprintf(stdout,"Stokes parameters will be written to files:\n");
2166 fprintf(stdout," %s [S0 (incident wave)],\n",outfilename_s0);
2167 fprintf(stdout," %s [S1 (incident wave)],\n",outfilename_s1);
2168 fprintf(stdout," %s [S2 (incident wave)],\n",outfilename_s2);
2169 fprintf(stdout," %s [S3 (incident wave)],\n",outfilename_s3);
2170 fprintf(stdout," %s [V0 (reflected wave)],\n",outfilename_v0);
2171 fprintf(stdout," %s [V1 (reflected wave)],\n",outfilename_v1);
2172 fprintf(stdout," %s [V2 (reflected wave)],\n",outfilename_v2);
2173 fprintf(stdout," %s [V3 (reflected wave)],\n",outfilename_v3);
2174 fprintf(stdout," %s [W0 (transmitted wave)],\n",outfilename_w0);
2175 fprintf(stdout," %s [W1 (transmitted wave)],\n",outfilename_w1);
2176 fprintf(stdout," %s [W2 (transmitted wave)],\n",outfilename_w2);
2177 fprintf(stdout," %s [W3 (transmitted wave)],\n",outfilename_w3);
2178 if(fieldevoflag){
2179 if(fieldevoflag_efield){
2180 if(strcmp(fieldevofilename,"")){
2181 fprintf(stdout,"Intra grating optical field evolution will "
2182 "be written to file:\n");
2183 fprintf(stdout," %s\n",fieldevofilename);
2184 }else{
2185 fprintf(stderr,
2186 "%s: Error: No file name specified for saving spatial\n"
2187 "field evolution. (efield option)\n",progname);
2188 exit(FAILURE);
2189 }
2190 fprintf(stdout,
2191 "(Intra grating field evolution will be presented in terms of\n"
2192 "the electrical field displacement.)\n");
2193 }else if(fieldevoflag_stoke){
2194 if(strcmp(fieldevofilename_s0,"")
2195 &&strcmp(fieldevofilename_s1,"")
2196 &&strcmp(fieldevofilename_s2,"")
2197 &&strcmp(fieldevofilename_s3,"")){
2198 fprintf(stdout,"Intra grating optical field evolution will "
2199 "be written to files:\n");
2200 fprintf(stdout," %s\n %s\n %s\n %s\n",
2201 fieldevofilename_s0,fieldevofilename_s1,
2202 fieldevofilename_s2,fieldevofilename_s3);
2203 }else{
2204 fprintf(stderr,
2205 "%s: Error: No file name specified for saving spatial\n"
2206 "field evolution. (stoke option)\n",progname);
2207 exit(FAILURE);
2208 }
2209 fprintf(stdout,
2210 "(Intra grating field evolution will be presented in terms of\n"
2211 "the Stokes parameters of the forward propagating field.)\n");
2212 }else{
2213 fprintf(stderr,"%s: Unknown field evolution flag.\n",progname);
2214 exit(FAILURE);
2215 }
2216 fprintf(stdout,
2217 "Number of intermediate samples within each layer: %-12ld\n",nne);
2218 }
2219 if(intensityevoflag){
2220 if(strcmp(intensityevofilename,"")){
2221 fprintf(stdout,"Intra grating optical intensity evolution will "
2222 "be written to file:\n");
2223 fprintf(stdout," %s\n",intensityevofilename);
2224 }
2225 }
2226 fprintf(stdout,"Program execution started %s",ctime(&initime));
2227 for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
2228 }
2229 }
2230
2231 /*:129*/
2232 #line 6119 "./magbragg.w"
2233
2234 }
2235
2236 /*:116*/
2237 #line 3415 "./magbragg.w"
2238
2239 /*134:*/
2240 #line 7435 "./magbragg.w"
2241
2242 {
2243 mmtraject= 0;
2244 if(trmtraject_specified){
2245 if((fp_traject= fopen(trmtraject_filename,"r"))==NULL){
2246 fprintf(stderr,
2247 "%s: Could not open file %s for reading Stokes parameter"
2248 " trajectory of transmitted wave!\n",
2249 progname,trmtraject_filename);
2250 exit(FAILURE);
2251 }
2252 fseek(fp_traject,0L,SEEK_SET);
2253
2254
2255 while((tmpch= getc(fp_traject))!=EOF){
2256 ungetc(tmpch,fp_traject);
2257 fscanf(fp_traject,"%lf",&tmp);
2258 fscanf(fp_traject,"%lf",&tmp);
2259 mmtraject++;
2260
2261
2262 tmpch= getc(fp_traject);
2263 while((tmpch!=EOF)&&(!numeric(tmpch))){
2264 tmpch= getc(fp_traject);
2265 }
2266 if(tmpch!=EOF)ungetc(tmpch,fp_traject);
2267 }
2268
2269 if(verbose){
2270 fprintf(stdout,
2271 "%s: I have now pre-parsed the specified trajectory of transmitted\n"
2272 "Stokes parameters (W0,W3) in file %s, and I found %-ld points.\n",
2273 progname,trmtraject_filename,mmtraject);
2274 fprintf(stdout,
2275 "%s: Now allocating the vectors for the transmitted trajectory...",
2276 progname);
2277 }
2278
2279 fseek(fp_traject,0L,SEEK_SET);
2280 w0traj= dvector(1,mmtraject);
2281 w3traj= dvector(1,mmtraject);
2282 for(ke= 1;ke<=mmtraject;ke++){
2283 fscanf(fp_traject,"%le",&w0traj[ke]);
2284 fscanf(fp_traject,"%le",&w3traj[ke]);
2285 }
2286
2287 if(0==1){
2288 for(ke= 1;ke<=mmtraject;ke++)
2289 fprintf(stdout,"w0=%e w3=%e\n",w0traj[ke],w3traj[ke]);
2290 }
2291
2292 fclose(fp_traject);
2293 }
2294 }
2295
2296 /*:134*/
2297 #line 3416 "./magbragg.w"
2298
2299 /*135:*/
2300 #line 7540 "./magbragg.w"
2301
2302 {
2303 if((mme> 1)&&(mmi> 1)){
2304 if((fp_s0= fopen(outfilename_s0,"w"))==NULL){
2305 fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
2306 progname,outfilename_s0);
2307 exit(FAILURE);
2308 }
2309 if((fp_s1= fopen(outfilename_s1,"w"))==NULL){
2310 fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
2311 progname,outfilename_s1);
2312 exit(FAILURE);
2313 }
2314 if((fp_s2= fopen(outfilename_s2,"w"))==NULL){
2315 fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
2316 progname,outfilename_s2);
2317 exit(FAILURE);
2318 }
2319 if((fp_s3= fopen(outfilename_s3,"w"))==NULL){
2320 fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
2321 progname,outfilename_s3);
2322 exit(FAILURE);
2323 }
2324 if((fp_v0= fopen(outfilename_v0,"w"))==NULL){
2325 fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
2326 progname,outfilename_v0);
2327 exit(FAILURE);
2328 }
2329 if((fp_v1= fopen(outfilename_v1,"w"))==NULL){
2330 fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
2331 progname,outfilename_v1);
2332 exit(FAILURE);
2333 }
2334 if((fp_v2= fopen(outfilename_v2,"w"))==NULL){
2335 fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
2336 progname,outfilename_v2);
2337 exit(FAILURE);
2338 }
2339 if((fp_v3= fopen(outfilename_v3,"w"))==NULL){
2340 fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
2341 progname,outfilename_v3);
2342 exit(FAILURE);
2343 }
2344 if((fp_w0= fopen(outfilename_w0,"w"))==NULL){
2345 fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
2346 progname,outfilename_w0);
2347 exit(FAILURE);
2348 }
2349 if((fp_w1= fopen(outfilename_w1,"w"))==NULL){
2350 fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
2351 progname,outfilename_w1);
2352 exit(FAILURE);
2353 }
2354 if((fp_w2= fopen(outfilename_w2,"w"))==NULL){
2355 fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
2356 progname,outfilename_w2);
2357 exit(FAILURE);
2358 }
2359 if((fp_w3= fopen(outfilename_w3,"w"))==NULL){
2360 fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
2361 progname,outfilename_w3);
2362 exit(FAILURE);
2363 }
2364 }
2365 if((mme> 1)&&(mmi> 1)){
2366 fseek(fp_s0,0L,SEEK_SET);
2367 fseek(fp_s1,0L,SEEK_SET);
2368 fseek(fp_s2,0L,SEEK_SET);
2369 fseek(fp_s3,0L,SEEK_SET);
2370 fseek(fp_v0,0L,SEEK_SET);
2371 fseek(fp_v1,0L,SEEK_SET);
2372 fseek(fp_v2,0L,SEEK_SET);
2373 fseek(fp_v3,0L,SEEK_SET);
2374 fseek(fp_w0,0L,SEEK_SET);
2375 fseek(fp_w1,0L,SEEK_SET);
2376 fseek(fp_w2,0L,SEEK_SET);
2377 fseek(fp_w3,0L,SEEK_SET);
2378 }
2379 if(fieldevoflag){
2380 if(fieldevoflag_efield){
2381 if(strcmp(fieldevofilename,"")){
2382 if((fp_evo= fopen(fieldevofilename,"w"))==NULL){
2383 fprintf(stderr,
2384 "%s: Could not open file %s for fieldevo output!\n",
2385 progname,fieldevofilename);
2386 exit(FAILURE);
2387 }
2388 }
2389 }else if(fieldevoflag_stoke){
2390 if(strcmp(fieldevofilename_s0,"")){
2391 if((fp_evo_s0= fopen(fieldevofilename_s0,"w"))==NULL){
2392 fprintf(stderr,
2393 "%s: Could not open file %s for saving spatial S0"
2394 " distribution!\n",progname,fieldevofilename);
2395 exit(FAILURE);
2396 }
2397 }else{
2398 fprintf(stderr,
2399 "%s: A name for the file for saving spatial S0 distribution"
2400 " is required!\n",progname);
2401 exit(FAILURE);
2402 }
2403 if(strcmp(fieldevofilename_s1,"")){
2404 if((fp_evo_s1= fopen(fieldevofilename_s1,"w"))==NULL){
2405 fprintf(stderr,
2406 "%s: Could not open file %s for saving spatial S1"
2407 " distribution!\n",progname,fieldevofilename);
2408 exit(FAILURE);
2409 }
2410 }else{
2411 fprintf(stderr,
2412 "%s: A name for the file for saving spatial S1 distribution"
2413 " is required!\n",progname);
2414 exit(FAILURE);
2415 }
2416 if(strcmp(fieldevofilename_s2,"")){
2417 if((fp_evo_s2= fopen(fieldevofilename_s2,"w"))==NULL){
2418 fprintf(stderr,
2419 "%s: Could not open file %s for saving spatial S2"
2420 " distribution!\n",progname,fieldevofilename);
2421 exit(FAILURE);
2422 }
2423 }else{
2424 fprintf(stderr,
2425 "%s: A name for the file for saving spatial S2 distribution"
2426 " is required!\n",progname);
2427 exit(FAILURE);
2428 }
2429 if(strcmp(fieldevofilename_s3,"")){
2430 if((fp_evo_s3= fopen(fieldevofilename_s3,"w"))==NULL){
2431 fprintf(stderr,
2432 "%s: Could not open file %s for saving spatial S3"
2433 " distribution!\n",progname,fieldevofilename);
2434 exit(FAILURE);
2435 }
2436 }else{
2437 fprintf(stderr,
2438 "%s: A name for the file for saving spatial S3 distribution"
2439 " is required!\n",progname);
2440 exit(FAILURE);
2441 }
2442 }else{
2443 fprintf(stderr,"%s: Unknown field evolution flag.\n",progname);
2444 exit(FAILURE);
2445 }
2446 }
2447 if(intensityevoflag){
2448 if(strcmp(intensityevofilename,"")){
2449 if((fp_ievo= fopen(intensityevofilename,"w"))==NULL){
2450 fprintf(stderr,
2451 "%s: Could not open file %s for intensityevo output!\n",
2452 progname,intensityevofilename);
2453 exit(FAILURE);
2454 }
2455 }
2456 }
2457 if(!((mme> 1)&&(mmi> 1))){
2458 if((fp_spec= fopen(spectrumfilename,"w"))==NULL){
2459 fprintf(stderr,"%s: Could not open file %s for spectrum output!\n",
2460 progname,spectrumfilename);
2461 exit(FAILURE);
2462 }
2463 if((fp_irspec= fopen(intensity_reflection_spectrumfilename,"w"))==NULL){
2464 fprintf(stderr,
2465 "%s: Could not open %s for intensity reflection spectrum!\n",
2466 progname,intensity_reflection_spectrumfilename);
2467 exit(FAILURE);
2468 }
2469 if((fp_itspec= fopen(intensity_transmission_spectrumfilename,"w"))==NULL){
2470 fprintf(stderr,
2471 "%s: Could not open %s for intensity transmission spectrum!\n",
2472 progname,intensity_transmission_spectrumfilename);
2473 exit(FAILURE);
2474 }
2475 if((fp_icspec= fopen(intensity_check_spectrumfilename,"w"))==NULL){
2476 fprintf(stderr,
2477 "%s: Could not open %s for checking spectra!\n",
2478 progname,intensity_check_spectrumfilename);
2479 exit(FAILURE);
2480 }
2481 }
2482 }
2483
2484 /*:135*/
2485 #line 3417 "./magbragg.w"
2486
2487 /*62:*/
2488 #line 3911 "./magbragg.w"
2489
2490 {
2491 z= dvector(1,nn);
2492 dz= dvector(1,nn-1);
2493 efp= dcvector(0,nn);
2494 efm= dcvector(0,nn);
2495 ebp= dcvector(0,nn);
2496 ebm= dcvector(0,nn);
2497 taup= dvector(1,nn);
2498 taum= dvector(1,nn);
2499 taupp= dvector(1,nn);
2500 taupm= dvector(1,nn);
2501 rhop= dvector(1,nn);
2502 rhom= dvector(1,nn);
2503 rhopp= dvector(1,nn);
2504 rhopm= dvector(1,nn);
2505 n= dvector(0,nn);
2506 g= dvector(0,nn);
2507 pe= dvector(1,nn-1);
2508 pm= dvector(1,nn-1);
2509 qe= dvector(1,nn-1);
2510 qm= dvector(1,nn-1);
2511 etafp= dvector(1,nn-1);
2512 etabp= dvector(1,nn-1);
2513 etafm= dvector(1,nn-1);
2514 etabm= dvector(1,nn-1);
2515 }
2516
2517 /*:62*/
2518 #line 3418 "./magbragg.w"
2519
2520 /*64:*/
2521 #line 3984 "./magbragg.w"
2522
2523 {
2524 if(!strcmp(gratingtype,"stepwise")){
2525 /*65:*/
2526 #line 4050 "./magbragg.w"
2527
2528 {
2529 if(!strcmp(gratingsubtype,"twolevel")){
2530 if(randomdistribution){
2531 for(j= 1;j<=nn-1;j++){
2532 if(j==modnum){
2533 if(j==1)z[j]= 0.0;
2534 z[j+1]= z[j]+modt1;
2535 dz[j]= modt1;
2536 n[j]= modn1;
2537 g[j]= modg1;
2538 pe[j]= modpe1;
2539 pm[j]= modpm1;
2540 qe[j]= modqe1;
2541 qm[j]= modqm1;
2542 }else{
2543 ranseed= ranseed+j;
2544 if(ran1(&ranseed)> 0.5){
2545 if(verbose)fprintf(stdout,"Random number: 1\n");
2546 if(j==1)z[j]= 0.0;
2547 z[j+1]= z[j]+t1;
2548 dz[j]= t1;
2549 n[j]= n1;
2550 g[j]= g1;
2551 pe[j]= pe1;
2552 pm[j]= pm1;
2553 qe[j]= qe1;
2554 qm[j]= qm1;
2555 }else{
2556 if(verbose)fprintf(stdout,"Random number: 0\n");
2557 if(j==1)z[j]= 0.0;
2558 z[j+1]= z[j]+t2;
2559 dz[j]= t2;
2560 n[j]= n2;
2561 g[j]= g2;
2562 pe[j]= pe2;
2563 pm[j]= pm2;
2564 qe[j]= qe2;
2565 qm[j]= qm2;
2566 }
2567 }
2568 }
2569 }else{
2570 if(modnum<0){
2571 for(j= 1;j<=nn-1;j= j+2){
2572 z[j]= 0.5*((double)(j-1))*(t1+t2);
2573 dz[j]= t1;
2574 n[j]= n1;
2575 g[j]= g1;
2576 pe[j]= pe1;
2577 pm[j]= pm1;
2578 qe[j]= qe1;
2579 qm[j]= qm1;
2580 if(j==nn-1)z[nn]= z[nn-1]+t1;
2581 }
2582 for(j= 2;j<=nn-1;j= j+2){
2583 z[j]= z[j-1]+t1;
2584 dz[j]= t2;
2585 n[j]= n2;
2586 g[j]= g2;
2587 pe[j]= pe2;
2588 pm[j]= pm2;
2589 qe[j]= qe2;
2590 qm[j]= qm2;
2591 if(j==nn-1)z[nn]= z[nn-1]+t2;
2592 }
2593 }else{
2594 for(j= 1;j<=nn-1;j++){
2595 if(j==1)z[j]= 0.0;
2596 if(j==modnum){
2597 z[j+1]= z[j]+modt1;
2598 dz[j]= modt1;
2599 n[j]= modn1;
2600 g[j]= modg1;
2601 pe[j]= modpe1;
2602 pm[j]= modpm1;
2603 qe[j]= modqe1;
2604 qm[j]= modqm1;
2605 }else{
2606 tmp= ((double)j)/((double)2);
2607 if(tmp-floor(tmp)> 0.25){
2608 z[j+1]= z[j]+t1;
2609 dz[j]= t1;
2610 n[j]= n1;
2611 g[j]= g1;
2612 pe[j]= pe1;
2613 pm[j]= pm1;
2614 qe[j]= qe1;
2615 qm[j]= qm1;
2616 }else{
2617 z[j+1]= z[j]+t2;
2618 dz[j]= t2;
2619 n[j]= n2;
2620 g[j]= g2;
2621 pe[j]= pe2;
2622 pm[j]= pm2;
2623 qe[j]= qe2;
2624 qm[j]= qm2;
2625 }
2626 }
2627 }
2628 }
2629 }
2630 }else{
2631 fprintf(stderr,"%s: Error.\n",progname);
2632 fprintf(stderr,"%s: (No valid grating subtype found).\n",progname);
2633 exit(FAILURE);
2634 }
2635 }
2636
2637 /*:65*/
2638 #line 3987 "./magbragg.w"
2639 ;
2640 }else if(!strcmp(gratingtype,"sinusoidal")){
2641 /*66:*/
2642 #line 4196 "./magbragg.w"
2643
2644 {
2645 t1= ll/((double)(nn-1));
2646 for(j= 1;j<=nn-1;j++){
2647 z[j]= ((double)(j-1))*t1;
2648 dz[j]= t1;
2649 if(apodize){
2650 if((0.0<=z[j])&&(z[j]<=apolength)){
2651 tmp= (1.0-cos(pi*z[j]/apolength))/2.0;
2652 }else if((apolength<=z[j])&&(z[j]<=ll-apolength)){
2653 tmp= 1.0;
2654 }else if((ll-apolength<=z[j])&&(z[j]<=ll)){
2655 tmp= (1.0-cos(pi*(z[j]-ll)/apolength))/2.0;
2656 }else{
2657 tmp= 0.0;
2658 fprintf(stderr,
2659 "%s: Impossible apodization event occurred.",progname);
2660 fprintf(stderr,
2661 "%s: (Please check grating initialization.)",progname);
2662 }
2663 }
2664 if(phasejump){
2665 if(z[j]>=phasejumpposition){
2666 phi= phasejumpangle;
2667 }else{
2668 phi= 0.0;
2669 }
2670 }
2671 if(apodize){
2672 n[j]= n1+n2*tmp*sin(twopi*z[j]/nper+phi);
2673 g[j]= g1+g2*tmp*sin(twopi*z[j]/gper+phi);
2674 }else{
2675 n[j]= n1+n2*sin(twopi*z[j]/nper+phi);
2676 g[j]= g1+g2*sin(twopi*z[j]/gper+phi);
2677 }
2678 pe[j]= pe1+pe2*sin(twopi*z[j]/peper);
2679 pm[j]= pm1+pm2*sin(twopi*z[j]/pmper);
2680 qe[j]= qe1+qe2*sin(twopi*z[j]/qeper);
2681 qm[j]= qm1+qm2*sin(twopi*z[j]/qmper);
2682 }
2683 z[nn]= ll;
2684 }
2685
2686 /*:66*/
2687 #line 3989 "./magbragg.w"
2688 ;
2689 }else if(!strcmp(gratingtype,"chirped")){
2690 /*67:*/
2691 #line 4249 "./magbragg.w"
2692
2693 {
2694 t1= ll/((double)(nn-1));
2695 for(j= 1;j<=nn-1;j++){
2696 z[j]= ((double)(j-1))*t1;
2697 dz[j]= t1;
2698 if(apodize){
2699 if((0.0<=z[j])&&(z[j]<=apolength)){
2700 tmp= (1.0-cos(pi*z[j]/apolength))/2.0;
2701 }else if((apolength<=z[j])&&(z[j]<=ll-apolength)){
2702 tmp= 1.0;
2703 }else if((ll-apolength<=z[j])&&(z[j]<=ll)){
2704 tmp= (1.0-cos(pi*(z[j]-ll)/apolength))/2.0;
2705 }else{
2706 tmp= 0.0;
2707 fprintf(stderr,
2708 "%s: Impossible apodization event occurred.",progname);
2709 fprintf(stderr,
2710 "%s: (Please check grating initialization.)",progname);
2711 }
2712 }
2713 if(phasejump){
2714 if(z[j]>=phasejumpposition){
2715 phi= phasejumpangle;
2716 }else{
2717 phi= 0.0;
2718 }
2719 }
2720 if(apodize){
2721 if(fabs(ncrp)> 0.0)
2722 n[j]= n1+n2*tmp*sin((twopi/ncrp)*log(1.0+ncrp*z[j]/nper)+phi);
2723 else
2724 n[j]= n1+n2*tmp*sin(twopi*z[j]/nper+phi);
2725 if(fabs(gcrp)> 0.0)
2726 g[j]= g1+g2*tmp*sin((twopi/gcrp)*log(1.0+gcrp*z[j]/gper)+phi);
2727 else
2728 g[j]= g1+g2*tmp*sin(twopi*z[j]/gper+phi);
2729 }else{
2730 if(fabs(ncrp)> 0.0)
2731 n[j]= n1+n2*sin((twopi/ncrp)*log(1.0+ncrp*z[j]/nper)+phi);
2732 else
2733 n[j]= n1+n2*sin(twopi*z[j]/nper+phi);
2734 if(fabs(gcrp)> 0.0)
2735 g[j]= g1+g2*sin((twopi/gcrp)*log(1.0+gcrp*z[j]/gper)+phi);
2736 else
2737 g[j]= g1+g2*sin(twopi*z[j]/gper+phi);
2738 }
2739 if(fabs(pecrp)> 0.0)
2740 pe[j]= pe1+pe2*sin((twopi/pecrp)*log(1.0+pecrp*z[j]/peper));
2741 else
2742 pe[j]= pe1+pe2*sin(twopi*z[j]/peper);
2743 if(pmcrp*pmcrp> 0.0)
2744 pm[j]= pm1+pm2*sin((twopi/pmcrp)*log(1.0+pmcrp*z[j]/pmper));
2745 else
2746 pm[j]= pm1+pm2*sin(twopi*z[j]/pmper);
2747 if(fabs(qecrp)> 0.0)
2748 qe[j]= qe1+qe2*sin((twopi/qecrp)*log(1.0+qecrp*z[j]/qeper));
2749 else
2750 qe[j]= qe1+qe2*sin(twopi*z[j]/qeper);
2751 if(fabs(qmcrp)> 0.0)
2752 qm[j]= qm1+qm2*sin((twopi/qmcrp)*log(1.0+qmcrp*z[j]/qmper));
2753 else
2754 qm[j]= qm1+qm2*sin(twopi*z[j]/qmper);
2755 }
2756 z[nn]= ll;
2757 }
2758
2759 /*:67*/
2760 #line 3991 "./magbragg.w"
2761 ;
2762 }else if(!strcmp(gratingtype,"fractal")){
2763 /*68:*/
2764 #line 4342 "./magbragg.w"
2765
2766 {
2767 tmp= 1.0;
2768 for(j= 1;j<=fractal_level-1;j++)tmp= 2.0*tmp;
2769 for(j= 1;j<=maximum_fractal_level-fractal_level;j++)tmp= 3.0*tmp;
2770
2771 ll= tmp*t1-tmp*t2;
2772 tmp= 1.0;
2773 for(j= 1;j<=maximum_fractal_level-1;j++)tmp= 3.0*tmp;
2774 ll= ll+tmp*t2;
2775 if(verbose){
2776 fprintf(stdout,"%s: Minimum layer thickness at maximum recursion depth:\n",
2777 progname);
2778 fprintf(stdout,"%s: t2=%f [nm]\n",progname,t1*1.0e9);
2779 fprintf(stdout,"%s: t2=%f [nm]\n",progname,t2*1.0e9);
2780 fprintf(stdout,"%s: Fractal grating length calculated as L=%e [m]\n",
2781 progname,ll);
2782 fprintf(stdout,
2783 "%s: Based on fractal recursion of %d out of a maximum of %d levels.\n",
2784 progname,fractal_level,maximum_fractal_level);
2785 }
2786 init_cantor_fractal_grating(z,1,nn,0.0,ll,n1,n2);
2787 for(j= 1;j<=nn-1;j++){
2788 dz[j]= z[j+1]-z[j];
2789 if(j==1){
2790 odd_layer= 1;
2791 }else{
2792 odd_layer= (odd_layer?0:1);
2793 }
2794 n[j]= (odd_layer?n1:n2);
2795 g[j]= (odd_layer?g1:g2);
2796 pe[j]= (odd_layer?pe1:pe2);
2797 pm[j]= (odd_layer?pm1:pm2);
2798 qe[j]= (odd_layer?qe1:qe2);
2799 qm[j]= (odd_layer?qm1:qm2);
2800 }
2801 }
2802
2803 /*:68*/
2804 #line 3993 "./magbragg.w"
2805 ;
2806 }else{
2807 fprintf(stderr,
2808 "%s: Error: Specified grating type is invalid.\n",progname);
2809 showsomehelp();
2810 exit(FAILURE);
2811 }
2812 if(perturbed_gyration_constant){
2813 /*69:*/
2814 #line 4400 "./magbragg.w"
2815 {
2816 for(j= 1;j<=nn-1;j++){
2817 tmp= 2.0*(z[j]-gyroperturb_position)/gyroperturb_width;
2818 g[j]+= gyroperturb_amplitude/(1.0+tmp*tmp);
2819 }
2820 }
2821
2822 /*:69*/
2823 #line 4001 "./magbragg.w"
2824 ;
2825 }
2826 }
2827
2828 /*:64*/
2829 #line 3419 "./magbragg.w"
2830
2831 /*70:*/
2832 #line 4426 "./magbragg.w"
2833
2834 {
2835 n[0]= nsurr;
2836 n[nn]= nsurr;
2837 g[0]= 0.0;
2838 g[nn]= 0.0;
2839 }
2840
2841 /*:70*/
2842 #line 3420 "./magbragg.w"
2843
2844 /*71:*/
2845 #line 4483 "./magbragg.w"
2846
2847 {
2848 for(j= 1;j<=nn;j++){
2849 taup[j]= 2.0*(n[j-1]+g[j-1])/(n[j-1]+n[j]+g[j-1]+g[j]);
2850 taum[j]= 2.0*(n[j-1]-g[j-1])/(n[j-1]+n[j]-g[j-1]-g[j]);
2851 taupp[j]= 2.0*(n[j]-g[j])/(n[j-1]+n[j]-g[j-1]-g[j]);
2852 taupm[j]= 2.0*(n[j]+g[j])/(n[j-1]+n[j]+g[j-1]+g[j]);
2853 rhop[j]= (n[j-1]-n[j]+g[j-1]-g[j])/(n[j-1]+n[j]+g[j-1]+g[j]);
2854 rhom[j]= (n[j-1]-n[j]-g[j-1]+g[j])/(n[j-1]+n[j]-g[j-1]-g[j]);
2855 rhopp[j]= -rhom[j];
2856 rhopm[j]= -rhop[j];
2857 }
2858 }
2859
2860 /*:71*/
2861 #line 3421 "./magbragg.w"
2862
2863 /*72:*/
2864 #line 4515 "./magbragg.w"
2865
2866 {
2867 for(k= 1;k<=mm;k++){
2868 if(mm> 1){
2869 lambda= lambdastart+(((double)(k-1))/((double)(mm-1)))
2870 *(lambdastop-lambdastart);
2871 }else{
2872 lambda= lambdastart;
2873 }
2874 omega= twopi*c/lambda;
2875 if(trmtraject_specified){
2876 mme= mmtraject;
2877 mmi= 1;
2878 }
2879 /*74:*/
2880 #line 4562 "./magbragg.w"
2881
2882 {
2883 for(ke= 1;ke<=mme;ke++){
2884 if(trmtraject_specified){
2885 trmellipticity= w3traj[ke]/w0traj[ke];
2886 }else{
2887 if(mme> 1){
2888 trmellipticity= trmellipstart+(((double)(ke-1))/((double)(mme-1)))
2889 *(trmellipstop-trmellipstart);
2890 }else{
2891 trmellipticity= trmellipstart;
2892 }
2893 }
2894 /*75:*/
2895 #line 4589 "./magbragg.w"
2896
2897 {
2898 for(ki= 1;ki<=mmi;ki++){
2899 if(trmtraject_specified){
2900 trmintensity= (epsilon0*c/2.0)*w0traj[ke];
2901 }else{
2902 if(mmi> 1){
2903 trmintensity= trmintenstart
2904 +(((double)(ki-1))/((double)(mmi-1)))
2905 *(trmintenstop-trmintenstart);
2906 }else{
2907 trmintensity= trmintenstart;
2908 }
2909 }
2910 /*76:*/
2911 #line 4650 "./magbragg.w"
2912
2913 {
2914 ebp[nn]= complex(0.0,0.0);
2915 ebm[nn]= complex(0.0,0.0);
2916 efp[nn]= complex(sqrt((1.0+trmellipticity)*trmintensity/(c*epsilon0)),0.0);
2917 efm[nn]= complex(sqrt((1.0-trmellipticity)*trmintensity/(c*epsilon0)),0.0);
2918 }
2919
2920 /*:76*/
2921 #line 4603 "./magbragg.w"
2922 ;
2923 /*77:*/
2924 #line 4662 "./magbragg.w"
2925
2926 {
2927 efp[nn-1]= cmul(crdiv(efp[nn],taup[nn]),crexpi(-omega*n[nn-1]*dz[nn-1]/c));
2928 efm[nn-1]= cmul(crdiv(efm[nn],taum[nn]),crexpi(-omega*n[nn-1]*dz[nn-1]/c));
2929 ebp[nn-1]= cmul(rcmul(rhom[nn],efm[nn-1]),
2930 crexpi(2.0*omega*n[nn-1]*dz[nn-1]/c));
2931 ebm[nn-1]= cmul(rcmul(rhop[nn],efp[nn-1]),
2932 crexpi(2.0*omega*n[nn-1]*dz[nn-1]/c));
2933 }
2934
2935 /*:77*/
2936 #line 4604 "./magbragg.w"
2937 ;
2938 /*78:*/
2939 #line 4702 "./magbragg.w"
2940
2941 {
2942 for(j= nn-1;j>=1;j--){
2943 /*79:*/
2944 #line 4762 "./magbragg.w"
2945
2946 {
2947 aafp2= cdabs(efp[j]);
2948 aafm2= cdabs(efm[j]);
2949 aabp2= cdabs(ebp[j]);
2950 aabm2= cdabs(ebm[j]);
2951 aafp2*= aafp2;
2952 aafm2*= aafm2;
2953 aabp2*= aabp2;
2954 aabm2*= aabm2;
2955 tmp= 3.0/(8.0*n[j]);
2956 etafp[j]= tmp*((pe[j]+pm[j])*(aafp2+2.0*aabm2)
2957 +(qe[j]+qm[j])*(aafm2+aabp2));
2958 etafm[j]= tmp*((pe[j]-pm[j])*(aafm2+2.0*aabp2)
2959 +(qe[j]-qm[j])*(aafp2+aabm2));
2960 etabp[j]= tmp*((pe[j]-pm[j])*(aabp2+2.0*aafm2)
2961 +(qe[j]-qm[j])*(aabm2+aafp2));
2962 etabm[j]= tmp*((pe[j]+pm[j])*(aabm2+2.0*aafp2)
2963 +(qe[j]+qm[j])*(aabp2+aafm2));
2964 }
2965
2966 /*:79*/
2967 #line 4705 "./magbragg.w"
2968 ;
2969 /*80:*/
2970 #line 4795 "./magbragg.w"
2971
2972 {
2973 tmpfp= cmul(efp[j],crexpi(-omega*(etafp[j]+g[j])*dz[j]/c));
2974 tmpfm= cmul(efm[j],crexpi(-omega*(etafm[j]-g[j])*dz[j]/c));
2975 tmpbp= cmul(ebp[j],crexpi(omega*(etabp[j]-g[j])*dz[j]/c));
2976 tmpbm= cmul(ebm[j],crexpi(omega*(etabm[j]+g[j])*dz[j]/c));
2977 }
2978
2979 /*:80*/
2980 #line 4706 "./magbragg.w"
2981 ;
2982 /*81:*/
2983 #line 4808 "./magbragg.w"
2984
2985 {
2986 if(j> 1){
2987 efp[j-1]= crdiv(cmul(csub(tmpfp,rcmul(rhopm[j],tmpbm)),
2988 crexpi(-omega*n[j-1]*dz[j-1]/c)),taup[j]);
2989 efm[j-1]= crdiv(cmul(csub(tmpfm,rcmul(rhopp[j],tmpbp)),
2990 crexpi(-omega*n[j-1]*dz[j-1]/c)),taum[j]);
2991 ebp[j-1]= cadd(rcmul(taupp[j],
2992 cmul(tmpbp,crexpi(omega*n[j-1]*dz[j-1]/c))),
2993 rcmul(rhom[j],
2994 cmul(efm[j-1],crexpi(2.0*omega*n[j-1]*dz[j-1]/c))));
2995 ebm[j-1]= cadd(rcmul(taupm[j],
2996 cmul(tmpbm,crexpi(omega*n[j-1]*dz[j-1]/c))),
2997 rcmul(rhop[j],
2998 cmul(efp[j-1],crexpi(2.0*omega*n[j-1]*dz[j-1]/c))));
2999 }else{
3000 efp[0]= crdiv(csub(tmpfp,rcmul(rhopm[1],tmpbm)),taup[1]);
3001 efm[0]= crdiv(csub(tmpfm,rcmul(rhopp[1],tmpbp)),taum[1]);
3002 ebp[0]= cadd(rcmul(taupp[1],tmpbp),rcmul(rhom[1],efm[0]));
3003 ebm[0]= cadd(rcmul(taupm[1],tmpbm),rcmul(rhop[1],efp[0]));
3004 }
3005 }
3006
3007 /*:81*/
3008 #line 4707 "./magbragg.w"
3009 ;
3010 if(verbose){
3011 /*82:*/
3012 #line 4852 "./magbragg.w"
3013
3014 {
3015 modf(100.0*((float)((nn-j-1)+(ki-1)*(nn-1)
3016 +(ke-1)*mmi*(nn-1)+(k-1)*mme*mmi*(nn-1)))
3017 /((float)(mm*mme*mmi*(nn-1))),&stn);
3018 if(stn> status+10){
3019 status= status+10;
3020 now= time(NULL);
3021 eta= initime
3022 +((int)((100.0/((double)status))*difftime(now,initime)));
3023 fprintf(stdout," ...%2d percent finished... ",status);
3024 fprintf(stdout," ETA: %s",ctime(&eta));
3025 }
3026 }
3027
3028 /*:82*/
3029 #line 4709 "./magbragg.w"
3030 ;
3031 }
3032 }
3033 }
3034
3035 /*:78*/
3036 #line 4605 "./magbragg.w"
3037 ;
3038 /*83:*/
3039 #line 5021 "./magbragg.w"
3040
3041 {
3042
3043
3044 if(intensityinfo){
3045 for(j= 0;j<=nn;j++){
3046 tmp= (epsilon0*n[j]*c/2)*(cabs2(efp[j])+cabs2(efm[j]));
3047 if(maxintens<tmp){
3048 maxintens= tmp;
3049 maxintens_layer= j;
3050 maxintens_inintens= (epsilon0*nsurr*c/2)
3051 *(cabs2(efp[0])+cabs2(efm[0]));
3052 maxintens_inellip= (cabs2(efp[0])-cabs2(efm[0]))
3053 /(cabs2(efp[0])+cabs2(efm[0]));
3054 maxintens_trintens= (epsilon0*nsurr*c/2)
3055 *(cabs2(efp[nn])+cabs2(efm[nn]));
3056 maxintens_trellip= (cabs2(efp[nn])-cabs2(efm[nn]))
3057 /(cabs2(efp[nn])+cabs2(efm[nn]));
3058 }
3059 }
3060 }
3061 if((mme> 1)&&(mmi> 1)){
3062
3063 s0= cabs2(efp[0])+cabs2(efm[0]);
3064 s1= 2.0*cmul(conjg(efp[0]),efm[0]).r;
3065 s2= 2.0*cmul(conjg(efp[0]),efm[0]).i;
3066 s3= cabs2(efp[0])-cabs2(efm[0]);
3067 if(scale_stokesparams){
3068 s0= s0*stoke_scalefactor;
3069 s1= s1*stoke_scalefactor;
3070 s2= s2*stoke_scalefactor;
3071 s3= s3*stoke_scalefactor;
3072 }
3073 if(normalize_ellipticity)s3= s3/s0;
3074 fprintf(fp_s0,"%16.12e ",s0);
3075 fprintf(fp_s1,"%16.12e ",s1);
3076 fprintf(fp_s2,"%16.12e ",s2);
3077 fprintf(fp_s3,"%16.12e ",s3);
3078 fflush(fp_s0);
3079 fflush(fp_s1);
3080 fflush(fp_s2);
3081 fflush(fp_s3);
3082
3083
3084 v0= cabs2(ebp[0])+cabs2(ebm[0]);
3085 v1= 2.0*cmul(conjg(ebp[0]),ebm[0]).r;
3086 v2= 2.0*cmul(conjg(ebp[0]),ebm[0]).i;
3087 v3= cabs2(ebp[0])-cabs2(ebm[0]);
3088 if(scale_stokesparams){
3089 v0= v0*stoke_scalefactor;
3090 v1= v1*stoke_scalefactor;
3091 v2= v2*stoke_scalefactor;
3092 v3= v3*stoke_scalefactor;
3093 }
3094 if(normalize_ellipticity)v3= v3/v0;
3095 fprintf(fp_v0,"%16.12e ",v0);
3096 fprintf(fp_v1,"%16.12e ",v1);
3097 fprintf(fp_v2,"%16.12e ",v2);
3098 fprintf(fp_v3,"%16.12e ",v3);
3099 fflush(fp_v0);
3100 fflush(fp_v1);
3101 fflush(fp_v2);
3102 fflush(fp_v3);
3103
3104
3105 w0= cabs2(efp[nn])+cabs2(efm[nn]);
3106 w1= 2.0*cmul(conjg(efp[nn]),efm[nn]).r;
3107 w2= 2.0*cmul(conjg(efp[nn]),efm[nn]).i;
3108 w3= cabs2(efp[nn])-cabs2(efm[nn]);
3109 if(scale_stokesparams){
3110 w0= w0*stoke_scalefactor;
3111 w1= w1*stoke_scalefactor;
3112 w2= w2*stoke_scalefactor;
3113 w3= w3*stoke_scalefactor;
3114 }
3115 if(normalize_ellipticity)w3= w3/w0;
3116 fprintf(fp_w0,"%16.12e ",w0);
3117 fprintf(fp_w1,"%16.12e ",w1);
3118 fprintf(fp_w2,"%16.12e ",w2);
3119 fprintf(fp_w3,"%16.12e ",w3);
3120 fflush(fp_w0);
3121 fflush(fp_w1);
3122 fflush(fp_w2);
3123 fflush(fp_w3);
3124 }else{
3125 if(stokes_parameter_spectrum){
3126
3127 s0= cabs2(efp[0])+cabs2(efm[0]);
3128 s1= 2.0*cmul(conjg(efp[0]),efm[0]).r;
3129 s2= 2.0*cmul(conjg(efp[0]),efm[0]).i;
3130 s3= cabs2(efp[0])-cabs2(efm[0]);
3131 fprintf(fp_spec,"%16.12e %16.12e %16.12e %16.12e %16.12e\n",
3132 lambda,omega,s1/s0,s2/s0,s3/s0);
3133 fflush(fp_spec);
3134
3135
3136 v0= cabs2(ebp[0])+cabs2(ebm[0]);
3137 v1= 2.0*cmul(conjg(ebp[0]),ebm[0]).r;
3138 v2= 2.0*cmul(conjg(ebp[0]),ebm[0]).i;
3139 v3= cabs2(ebp[0])-cabs2(ebm[0]);
3140
3141
3142
3143 w0= cabs2(efp[nn])+cabs2(efm[nn]);
3144 w1= 2.0*cmul(conjg(efp[nn]),efm[nn]).r;
3145 w2= 2.0*cmul(conjg(efp[nn]),efm[nn]).i;
3146 w3= cabs2(efp[nn])-cabs2(efm[nn]);
3147
3148 }
3149 if(save_dbspectra){
3150 tmp= (cabs2(ebp[0])+cabs2(ebm[0]))/(cabs2(efp[0])+cabs2(efm[0]));
3151 fprintf(fp_irspec,"%-10.8f %-10.8f\n",lambda*1.0e9,10.0*log10(tmp));
3152 fprintf(fp_itspec,"%-10.8f %-10.8f\n",lambda*1.0e9,10.0*log10(1.0-tmp));
3153 fflush(fp_irspec);
3154 fflush(fp_itspec);
3155 }else{
3156 fprintf(fp_irspec,"%-10.8f %-10.8f\n",lambda*1.0e9,
3157 (cabs2(ebp[0])+cabs2(ebm[0]))/(cabs2(efp[0])+cabs2(efm[0])));
3158 fprintf(fp_itspec,"%-10.8f %-10.8f\n",lambda*1.0e9,
3159 (cabs2(efp[nn])+cabs2(efm[nn]))/(cabs2(efp[0])+cabs2(efm[0])));
3160 fflush(fp_irspec);
3161 fflush(fp_itspec);
3162 }
3163 fprintf(fp_icspec,"%-10.8f %-10.8f\n",lambda*1.0e9,
3164 (cabs2(ebp[0])+cabs2(ebm[0]))/(cabs2(efp[0])+cabs2(efm[0]))
3165 +(cabs2(efp[nn])+cabs2(efm[nn]))/(cabs2(efp[0])+cabs2(efm[0])));
3166 fflush(fp_icspec);
3167 }
3168
3169 if(ki>=mmi){
3170 if((mme> 1)&&(mmi> 1)){
3171 fprintf(fp_s0,"\n");
3172 fprintf(fp_s1,"\n");
3173 fprintf(fp_s2,"\n");
3174 fprintf(fp_s3,"\n");
3175 fprintf(fp_v0,"\n");
3176 fprintf(fp_v1,"\n");
3177 fprintf(fp_v2,"\n");
3178 fprintf(fp_v3,"\n");
3179 fprintf(fp_w0,"\n");
3180 fprintf(fp_w1,"\n");
3181 fprintf(fp_w2,"\n");
3182 fprintf(fp_w3,"\n");
3183 }
3184 }
3185 }
3186
3187 /*:83*/
3188 #line 4606 "./magbragg.w"
3189 ;
3190 /*84:*/
3191 #line 5179 "./magbragg.w"
3192
3193 {
3194 if(fieldevoflag){
3195 if(fieldevoflag_efield){
3196 if(verbose)
3197 fprintf(stdout,
3198 "Writing spatial field evolution to file.\n");
3199 if(fp_evo!=NULL){
3200 for(j= 1;j<=nn-1;j++){
3201 for(jje= 1;jje<=nne;jje++){
3202 if(nne> 1){
3203 zt= z[j]+((double)(jje-1))*dz[j]/((double)(nne));
3204 }else{
3205 zt= z[j];
3206 }
3207 tmpfp= cmul(efp[j],
3208 crexpi(omega*(etafp[j]+g[j])*(zt-z[j])/c));
3209 tmpfm= cmul(efm[j],
3210 crexpi(omega*(etafm[j]-g[j])*(zt-z[j])/c));
3211 tmpbp= cmul(ebp[j],
3212 crexpi(-omega*(etabp[j]-g[j])*(zt-z[j])/c));
3213 tmpbm= cmul(ebm[j],
3214 crexpi(-omega*(etabm[j]+g[j])*(zt-z[j])/c));
3215 if(normalize_length_to_micrometer){
3216 fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e"
3217 " %16.12e %16.12e %16.12e %16.12e"
3218 " %16.12e\n",zt*1.0e6,
3219 tmpfp.r,tmpfp.i,tmpfm.r,tmpfm.i,
3220 tmpbp.r,tmpbp.i,tmpbm.r,tmpbm.i);
3221 }else{
3222 fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e"
3223 " %16.12e %16.12e %16.12e %16.12e"
3224 " %16.12e\n",zt,
3225 tmpfp.r,tmpfp.i,tmpfm.r,tmpfm.i,
3226 tmpbp.r,tmpbp.i,tmpbm.r,tmpbm.i);
3227 }
3228 }
3229 }
3230 if(normalize_length_to_micrometer){
3231 fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e %16.12e"
3232 " %16.12e %16.12e %16.12e %16.12e\n",
3233 z[nn]*1.0e6,efp[nn].r,efp[nn].i,efm[nn].r,efm[nn].i,
3234 ebp[nn].r,ebp[nn].i,ebm[nn].r,ebm[nn].i);
3235 }else{
3236 fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e %16.12e"
3237 " %16.12e %16.12e %16.12e %16.12e\n",
3238 z[nn],efp[nn].r,efp[nn].i,efm[nn].r,efm[nn].i,
3239 ebp[nn].r,ebp[nn].i,ebm[nn].r,ebm[nn].i);
3240 }
3241 }else{
3242 fprintf(stderr,"%s: Could not write to file %s!\n",
3243 progname,fieldevofilename);
3244 exit(FAILURE);
3245 }
3246 }else if(fieldevoflag_stoke){
3247 if(verbose)
3248 fprintf(stdout,
3249 "Writing spatial evolution of Stokes parameters to file.\n");
3250 if((fp_evo_s0!=NULL)&&(fp_evo_s1!=NULL)
3251 &&(fp_evo_s2!=NULL)&&(fp_evo_s3!=NULL)){
3252 for(j= 1;j<=nn-1;j++){
3253 for(jje= 1;jje<=nne;jje++){
3254 if(nne> 1){
3255 zt= z[j]+((double)(jje-1))*dz[j]/((double)(nne));
3256 }else{
3257 zt= z[j];
3258 }
3259 tmpfp= cmul(efp[j],
3260 crexpi(omega*(etafp[j]+g[j])*(zt-z[j])/c));
3261 tmpfm= cmul(efm[j],
3262 crexpi(omega*(etafm[j]-g[j])*(zt-z[j])/c));
3263 tmpbp= cmul(ebp[j],
3264 crexpi(-omega*(etabp[j]-g[j])*(zt-z[j])/c));
3265 tmpbm= cmul(ebm[j],
3266 crexpi(-omega*(etabm[j]+g[j])*(zt-z[j])/c));
3267 s0= cabs2(tmpfp)+cabs2(tmpfm);
3268 s1= 2.0*cmul(conjg(tmpfp),tmpfm).r;
3269 s2= 2.0*cmul(conjg(tmpfp),tmpfm).i;
3270 s3= cabs2(tmpfp)-cabs2(tmpfm);
3271 if(normalize_intensity)s0= s0/(cabs2(efp[1])+cabs2(efm[1]));
3272 if(normalize_length_to_micrometer){
3273 fprintf(fp_evo_s0,"%16.12e %16.12e\n",zt*1.0e6,s0);
3274 fprintf(fp_evo_s1,"%16.12e %16.12e\n",zt*1.0e6,s1);
3275 fprintf(fp_evo_s2,"%16.12e %16.12e\n",zt*1.0e6,s2);
3276 fprintf(fp_evo_s3,"%16.12e %16.12e\n",zt*1.0e6,s3);
3277 }else{
3278 fprintf(fp_evo_s0,"%16.12e %16.12e\n",zt,s0);
3279 fprintf(fp_evo_s1,"%16.12e %16.12e\n",zt,s1);
3280 fprintf(fp_evo_s2,"%16.12e %16.12e\n",zt,s2);
3281 fprintf(fp_evo_s3,"%16.12e %16.12e\n",zt,s3);
3282 }
3283 }
3284 }
3285 s0= cabs2(efp[nn])+cabs2(efm[nn]);
3286 s1= 2.0*cmul(conjg(efp[nn]),efm[nn]).r;
3287 s2= 2.0*cmul(conjg(efp[nn]),efm[nn]).i;
3288 s3= cabs2(efp[nn])-cabs2(efm[nn]);
3289 if(normalize_intensity)s0= s0/(cabs2(efp[1])+cabs2(efm[1]));
3290 if(normalize_length_to_micrometer){
3291 fprintf(fp_evo_s0,"%16.12e %16.12e\n",z[nn]*1.0e6,s0);
3292 fprintf(fp_evo_s1,"%16.12e %16.12e\n",z[nn]*1.0e6,s1);
3293 fprintf(fp_evo_s2,"%16.12e %16.12e\n",z[nn]*1.0e6,s2);
3294 fprintf(fp_evo_s3,"%16.12e %16.12e\n",z[nn]*1.0e6,s3);
3295 }else{
3296 fprintf(fp_evo_s0,"%16.12e %16.12e\n",z[nn],s0);
3297 fprintf(fp_evo_s1,"%16.12e %16.12e\n",z[nn],s1);
3298 fprintf(fp_evo_s2,"%16.12e %16.12e\n",z[nn],s2);
3299 fprintf(fp_evo_s3,"%16.12e %16.12e\n",z[nn],s3);
3300 }
3301 }else{
3302 fprintf(stderr,"%s: Could not write to file %s, %s, %s, or %s!\n",
3303 progname,fieldevofilename_s0,fieldevofilename_s1,
3304 fieldevofilename_s2,fieldevofilename_s3);
3305 exit(FAILURE);
3306 }
3307 }else{
3308 fprintf(stderr,"%s: Unknown field evolution flag.\n"
3309 "%s: (This cannot happen)\n",progname,progname);
3310 exit(FAILURE);
3311 }
3312 }
3313 }
3314
3315 /*:84*/
3316 #line 4607 "./magbragg.w"
3317 ;
3318 /*85:*/
3319 #line 5306 "./magbragg.w"
3320
3321 {
3322 if(intensityevoflag){
3323 if(fabs(ievolambda-lambda)
3324 <fabs(lambdastop-lambdastart)/((double)(mm))){
3325 if(verbose)
3326 fprintf(stdout,"%s: Saving intensity evolution at lambda=%8.4e\n",
3327 progname,lambda);
3328 if(strcmp(intensityevofilename,"")){
3329 if(fp_ievo!=NULL){
3330 for(j= 1;j<=nn-1;j++){
3331 if(normalize_length_to_micrometer){
3332 fprintf(fp_ievo,"%16.12e %16.12e\n",z[j]*1.0e6,
3333 (cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
3334 /(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
3335 fprintf(fp_ievo,"%16.12e %16.12e\n",z[j+1]*1.0e6,
3336 (cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
3337 /(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
3338 }else{
3339 fprintf(fp_ievo,"%16.12e %16.12e\n",z[j],
3340 (cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
3341 /(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
3342 fprintf(fp_ievo,"%16.12e %16.12e\n",z[j+1],
3343 (cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
3344 /(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
3345 }
3346 }
3347 }else{
3348 fprintf(stderr,"%s: Could not write to file %s!\n",
3349 progname,intensityevofilename);
3350 exit(FAILURE);
3351 }
3352 }
3353 }
3354 }
3355 }
3356
3357 /*:85*/
3358 #line 4608 "./magbragg.w"
3359 ;
3360 /*86:*/
3361 #line 5348 "./magbragg.w"
3362
3363 {
3364 if(writegratingtofile){
3365 if((fp_gr= fopen(gratingfilename,"w"))==NULL){
3366 fprintf(stderr,"%s: Could not open file %s for output!\n",
3367 progname,gratingfilename);
3368 exit(FAILURE);
3369 }
3370 if(normalize_length_to_micrometer){
3371 if(display_surrounding_media){
3372 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3373 " %16.12e %16.12e %16.12e\n",
3374 (z[1]-0.1*(z[nn]-z[1]))*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
3375 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3376 " %16.12e %16.12e %16.12e\n",
3377 z[1]*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
3378 }
3379 for(j= 1;j<=nn-1;j++){
3380 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3381 " %16.12e %16.12e %16.12e\n",
3382 z[j]*1.0e6,n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
3383 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3384 " %16.12e %16.12e %16.12e\n",
3385 z[j+1]*1.0e6,n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
3386 }
3387 if(display_surrounding_media){
3388 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3389 " %16.12e %16.12e %16.12e\n",
3390 z[nn]*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
3391 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3392 " %16.12e %16.12e %16.12e\n",
3393 (z[nn]+0.1*(z[nn]-z[1]))*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
3394 }
3395 }else{
3396 if(display_surrounding_media){
3397 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3398 " %16.12e %16.12e %16.12e\n",
3399 z[1]-0.1*(z[nn]-z[1]),nsurr,0.0,0.0,0.0,0.0,0.0);
3400 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3401 " %16.12e %16.12e %16.12e\n",
3402 z[1],nsurr,0.0,0.0,0.0,0.0,0.0);
3403 }
3404 for(j= 1;j<=nn-1;j++){
3405 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3406 " %16.12e %16.12e %16.12e\n",
3407 z[j],n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
3408 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3409 " %16.12e %16.12e %16.12e\n",
3410 z[j+1],n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
3411 }
3412 if(display_surrounding_media){
3413 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3414 " %16.12e %16.12e %16.12e\n",
3415 z[nn],nsurr,0.0,0.0,0.0,0.0,0.0);
3416 fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
3417 " %16.12e %16.12e %16.12e\n",
3418 z[nn]+0.1*(z[nn]-z[1]),nsurr,0.0,0.0,0.0,0.0,0.0);
3419 }
3420 }
3421 fclose(fp_gr);
3422 }
3423 }
3424
3425 /*:86*/
3426 #line 4609 "./magbragg.w"
3427 ;
3428 }
3429 }
3430
3431 /*:75*/
3432 #line 4575 "./magbragg.w"
3433 ;
3434 }
3435 }
3436
3437 /*:74*/
3438 #line 4529 "./magbragg.w"
3439 ;
3440 }
3441 if(verbose)/*73:*/
3442 #line 4538 "./magbragg.w"
3443
3444 {
3445 fprintf(stdout," ...done. ");
3446 now= time(NULL);
3447 fprintf(stdout,"Elapsed execution time: %d s\n",
3448 ((int)difftime(now,initime)));
3449 for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
3450 fprintf(stdout,"Program execution closed %s",ctime(&now));
3451 }
3452
3453 /*:73*/
3454 #line 4531 "./magbragg.w"
3455 ;
3456 }
3457
3458 /*:72*/
3459 #line 3422 "./magbragg.w"
3460
3461 /*87:*/
3462 #line 5413 "./magbragg.w"
3463
3464 {
3465 if(intensityinfo){
3466 for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
3467 fprintf(stdout,"Summary of intra-grating intensities:\n");
3468 fprintf(stdout,
3469 "The maximum intensity %1.4e [W/m^2] (%1.4f GW/cm^2) was detected\n",
3470 maxintens,maxintens*1.0e-13);
3471 fprintf(stdout,
3472 "in layer %ld. The maximum intensity was detected at a transmitted\n",
3473 maxintens_layer);
3474 fprintf(stdout,
3475 "intensity of %1.4e [W/m^2] (%1.4f GW/cm^2), and at a transmitted\n",
3476 maxintens_trintens,maxintens_trintens*1.0e-13);
3477 fprintf(stdout,
3478 "normalized ellipticity of S3/S0=%1.4f. ",maxintens_trellip);
3479 fprintf(stdout,"(where S3/S0=1 for LCP, and -1\n""for RCP).");
3480 fprintf(stdout,
3481 "At this state, the calculated optical intensity incident to the\n");
3482 fprintf(stdout,
3483 "crystal was %1.4e [W/m^2] (%1.4f GW/cm^2), or %1.1f percent\n",
3484 maxintens_inintens,maxintens_inintens*1.0e-13,
3485 100.0*maxintens_inintens/maxintens);
3486 fprintf(stdout,"of the maximum intra-grating optical intensity.\n");
3487 fprintf(stdout,
3488 "The calculated normalized incident ellipticity was %1.4f.\n",
3489 maxintens_inellip);
3490 fprintf(stdout,
3491 "The intensity transmission was for this state %1.1f percent.\n",
3492 100.0*maxintens_trintens/maxintens_inintens);
3493 for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
3494 if(saveintensityinfologfile){
3495 if((intensinfologfile= fopen(intensinfologfilename,"w"))==NULL){
3496 fprintf(stderr,"%s: Could not open %s for intensity log!\n",
3497 progname,intensinfologfilename);
3498 exit(FAILURE);
3499 }
3500 for(k= 1;k<=32;k++)fprintf(intensinfologfile,(k<32?"-":"\n"));
3501 fprintf(intensinfologfile,"Summary of intra-grating intensities:\n");
3502 fprintf(intensinfologfile,
3503 "Maximum intensity %1.4e [W/sq.m] (%1.4f GW/sq.cm) was detected\n",
3504 maxintens,maxintens*1.0e-13);
3505 fprintf(intensinfologfile,
3506 "in layer %ld. Maximum intensity was detected at a transmitted\n",
3507 maxintens_layer);
3508 fprintf(intensinfologfile,
3509 "intensity of %1.4e [W/sq.m] (%1.4f GW/sq.cm), and at transmitted\n",
3510 maxintens_trintens,maxintens_trintens*1.0e-13);
3511 fprintf(intensinfologfile,
3512 "normalized ellipticity of S3/S0=%1.4f. ",maxintens_trellip);
3513 fprintf(intensinfologfile,
3514 "(where S3/S0=1 for LCP, and -1\n""for RCP). ");
3515 fprintf(intensinfologfile,
3516 "At this state, the calculated optical intensity incident to the\n");
3517 fprintf(intensinfologfile,
3518 "crystal was %1.4e [W/sq.m] (%1.4f GW/sq.cm), or %1.1f percent\n",
3519 maxintens_inintens,maxintens_inintens*1.0e-13,
3520 100.0*maxintens_inintens/maxintens);
3521 fprintf(intensinfologfile,
3522 "of the maximum intra-grating optical intensity.\n");
3523 fprintf(intensinfologfile,
3524 "The calculated normalized incident ellipticity was %1.4f.\n",
3525 maxintens_inellip);
3526 fprintf(intensinfologfile,
3527 "The intensity transmission was for this state %1.1f percent.\n",
3528 100.0*maxintens_trintens/maxintens_inintens);
3529 fclose(intensinfologfile);
3530 for(k= 1;k<=32;k++)fprintf(intensinfologfile,(k<32?"-":"\n"));
3531 }
3532 }
3533 }
3534
3535 /*:87*/
3536 #line 3423 "./magbragg.w"
3537
3538 /*63:*/
3539 #line 3944 "./magbragg.w"
3540
3541 {
3542 free_dcvector(efp,0,nn);
3543 free_dcvector(efm,0,nn);
3544 free_dcvector(ebp,0,nn);
3545 free_dcvector(ebm,0,nn);
3546 free_dvector(taup,1,nn);
3547 free_dvector(taum,1,nn);
3548 free_dvector(taupp,1,nn);
3549 free_dvector(taupm,1,nn);
3550 free_dvector(rhop,1,nn);
3551 free_dvector(rhom,1,nn);
3552 free_dvector(rhopp,1,nn);
3553 free_dvector(rhopm,1,nn);
3554 free_dvector(n,0,nn);
3555 free_dvector(g,0,nn);
3556 free_dvector(pe,1,nn-1);
3557 free_dvector(pm,1,nn-1);
3558 free_dvector(qe,1,nn-1);
3559 free_dvector(qm,1,nn-1);
3560 free_dvector(etafp,1,nn-1);
3561 free_dvector(etabp,1,nn-1);
3562 free_dvector(etafm,1,nn-1);
3563 free_dvector(etabm,1,nn-1);
3564 if(trmtraject_specified){
3565 free_dvector(w0traj,1,mmtraject);
3566 free_dvector(w3traj,1,mmtraject);
3567 }
3568 }
3569
3570 /*:63*/
3571 #line 3424 "./magbragg.w"
3572
3573 /*136:*/
3574 #line 7725 "./magbragg.w"
3575
3576 {
3577 if((mme> 1)&&(mmi> 1)){
3578 fclose(fp_s0);
3579 fclose(fp_s1);
3580 fclose(fp_s2);
3581 fclose(fp_s3);
3582 fclose(fp_v0);
3583 fclose(fp_v1);
3584 fclose(fp_v2);
3585 fclose(fp_v3);
3586 fclose(fp_w0);
3587 fclose(fp_w1);
3588 fclose(fp_w2);
3589 fclose(fp_w3);
3590 }
3591 if(fieldevoflag)if(strcmp(fieldevofilename,""))fclose(fp_evo);
3592 if(intensityevoflag)if(strcmp(intensityevofilename,""))fclose(fp_ievo);
3593 if(!((mme> 1)&&(mmi> 1))){
3594 fclose(fp_spec);
3595 fclose(fp_irspec);
3596 fclose(fp_itspec);
3597 fclose(fp_icspec);
3598 }
3599 }
3600
3601 /*:136*/
3602 #line 3425 "./magbragg.w"
3603
3604 return(SUCCESS);
3605 }
3606
3607 /*:45*/
3608
Generated by ::viewsrc::