/*45:*/
#line 3404 "./magbragg.w"

/*46:*/
#line 3437 "./magbragg.w"

#include <math.h> 
#include <time.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include <string.h> 
#include <ctype.h> 

/*:46*/
#line 3405 "./magbragg.w"

/*47:*/
#line 3454 "./magbragg.w"

#define VERSION "1.43"
#define COPYRIGHT "Copyright (C) 2002-2007, Fredrik Jonsson"
#define SUCCESS (0)
#define FAILURE (1)
#define NCHMAX (256)

/*:47*/
#line 3406 "./magbragg.w"

/*48:*/
#line 3466 "./magbragg.w"

typedef struct DCOMPLEX{double r,i;}dcomplex;

/*:48*/
#line 3407 "./magbragg.w"

/*49:*/
#line 3475 "./magbragg.w"

extern char*optarg;
char*progname;

/*:49*/
#line 3408 "./magbragg.w"

/*50:*/
#line 3481 "./magbragg.w"

/*88:*/
#line 5491 "./magbragg.w"

int numeric(char ch){
if((ch=='0')||(ch=='1')||(ch=='2')||(ch=='3')||(ch=='4')||(ch=='5')
||(ch=='6')||(ch=='7')||(ch=='8')||(ch=='9')||(ch=='+')||(ch=='-')){
return(1);
}else{
return(0);
}
}

/*:88*/
#line 3482 "./magbragg.w"

/*89:*/
#line 5503 "./magbragg.w"

void init_cantor_fractal_grating(double*z,int indxmin,int indxmax,
double llmin,double llmax,double n1,double n2){
int indxmid;
double dll;
if(indxmax-indxmin==1){
z[indxmin]= llmin;
z[indxmax]= llmax;
}else if(indxmax-indxmin>=3){
indxmid= (indxmin+indxmax)/2;
dll= (n2/(n1+2*n2))*(llmax-llmin);
init_cantor_fractal_grating(z,indxmin,indxmid,llmin,llmin+dll,n1,n2);
init_cantor_fractal_grating(z,indxmid+1,indxmax,llmax-dll,llmax,n1,n2);
}else{
fprintf(stderr,"%s: Error in indexes detected in routine %s",progname,
"init_cantor_fractal_grating()!\n");
fprintf(stderr,"%s: Please verify that the number of discrete layers\n",
progname);
fprintf(stderr,"%s: in the grating is N-1, where N is an integer"
" power of 2.\n",progname);
exit(FAILURE);
}
}

/*:89*/
#line 3483 "./magbragg.w"

/*90:*/
#line 5536 "./magbragg.w"

/*91:*/
#line 5545 "./magbragg.w"

short pathcharacter(int ch){
return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
(ch=='-')||(ch=='+'));
}

/*:91*/
#line 5537 "./magbragg.w"

/*92:*/
#line 5557 "./magbragg.w"

char*strip_away_path(char filename[]){
int j,k= 0;
while(pathcharacter(filename[k]))k++;
j= (--k);
while(isalnum((int)(filename[j])))j--;
j++;
return(&filename[j]);
}

/*:92*/
#line 5538 "./magbragg.w"


/*:90*/
#line 3484 "./magbragg.w"

/*111:*/
#line 5895 "./magbragg.w"

/*112:*/
#line 5906 "./magbragg.w"

double*dvector(long nl,long nh){
double*v;
v= (double*)malloc((size_t)((nh-nl+2)*sizeof(double)));
if(!v){
fprintf(stderr,"Error: Allocation failure in dvector()\n");
exit(FAILURE);
}
return v-nl+1;
}

/*:112*/
#line 5896 "./magbragg.w"

/*113:*/
#line 5922 "./magbragg.w"

dcomplex*dcvector(long nl,long nh){
dcomplex*v;
v= (dcomplex*)malloc((size_t)((nh-nl+2)*sizeof(dcomplex)));
if(!v){
fprintf(stderr,"Error: Allocation failure in dcvector()\n");
exit(FAILURE);
}
return v-nl+1;
}

/*:113*/
#line 5897 "./magbragg.w"

/*114:*/
#line 5936 "./magbragg.w"

void free_dvector(double*v,long nl,long nh){
free((char*)(v+nl-1));
}

/*:114*/
#line 5898 "./magbragg.w"

/*115:*/
#line 5944 "./magbragg.w"

void free_dcvector(dcomplex*v,long nl,long nh){
free((char*)(v+nl-1));
}

/*:115*/
#line 5899 "./magbragg.w"


/*:111*/
#line 3485 "./magbragg.w"

/*93:*/
#line 5571 "./magbragg.w"

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran1(long*idum)
{
int j;
long k;
static long iy= 0;
static long iv[NTAB];
float temp;

if(*idum<=0||!iy){
if(-(*idum)<1)*idum= 1;
else*idum= -(*idum);
for(j= NTAB+7;j>=0;j--){
k= (*idum)/IQ;
*idum= IA*(*idum-k*IQ)-IR*k;
if(*idum<0)*idum+= IM;
if(j<NTAB)iv[j]= *idum;
}
iy= iv[0];
}
k= (*idum)/IQ;
*idum= IA*(*idum-k*IQ)-IR*k;
if(*idum<0)*idum+= IM;
j= iy/NDIV;
iy= iv[j];
iv[j]= *idum;
if((temp= AM*iy)> RNMX)return RNMX;
else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX

/*:93*/
#line 3486 "./magbragg.w"

/*94:*/
#line 5626 "./magbragg.w"

/*95:*/
#line 5643 "./magbragg.w"

dcomplex complex(double re,double im){
dcomplex c;
c.r= re;
c.i= im;
return c;
}

/*:95*/
#line 5627 "./magbragg.w"

/*96:*/
#line 5655 "./magbragg.w"

dcomplex conjg(dcomplex z){
dcomplex c;
c.r= z.r;
c.i= -z.i;
return c;
}

/*:96*/
#line 5628 "./magbragg.w"

/*97:*/
#line 5667 "./magbragg.w"

double cdabs(dcomplex z){
double x,y,c,tmp;
x= fabs(z.r);
y= fabs(z.i);
if(x==0.0)
c= y;
else if(y==0.0)
c= x;
else if(x> y){
tmp= y/x;
c= x*sqrt(1.0+tmp*tmp);
}else{
tmp= x/y;
c= y*sqrt(1.0+tmp*tmp);
}
return c;
}

/*:97*/
#line 5629 "./magbragg.w"

/*98:*/
#line 5691 "./magbragg.w"

double cabs2(dcomplex z){
double x,y,c,tmp;
x= fabs(z.r);
y= fabs(z.i);
if(x==0.0)
c= y*y;
else if(y==0.0)
c= x*x;
else if(x> y){
tmp= y/x;
c= x*x*(1.0+tmp*tmp);
}else{
tmp= x/y;
c= y*y*(1.0+tmp*tmp);
}
return c;
}

/*:98*/
#line 5630 "./magbragg.w"

/*99:*/
#line 5713 "./magbragg.w"

dcomplex cadd(dcomplex a,dcomplex b){
dcomplex c;
c.r= a.r+b.r;
c.i= a.i+b.i;
return c;
}

/*:99*/
#line 5631 "./magbragg.w"

/*100:*/
#line 5724 "./magbragg.w"

dcomplex csub(dcomplex a,dcomplex b){
dcomplex c;
c.r= a.r-b.r;
c.i= a.i-b.i;
return c;
}

/*:100*/
#line 5632 "./magbragg.w"

/*101:*/
#line 5737 "./magbragg.w"

/*102:*/
#line 5744 "./magbragg.w"

dcomplex cmul(dcomplex a,dcomplex b){
dcomplex c;
c.r= a.r*b.r-a.i*b.i;
c.i= a.i*b.r+a.r*b.i;
return c;
}

/*:102*/
#line 5738 "./magbragg.w"

/*103:*/
#line 5756 "./magbragg.w"

dcomplex rcmul(double a,dcomplex b){
dcomplex c;
c.r= a*b.r;
c.i= a*b.i;
return c;
}

/*:103*/
#line 5739 "./magbragg.w"


/*:101*/
#line 5633 "./magbragg.w"

/*104:*/
#line 5769 "./magbragg.w"

/*105:*/
#line 5779 "./magbragg.w"

dcomplex cdiv(dcomplex a,dcomplex b){
dcomplex c;
double r,den;
if(cdabs(b)> 0.0){
if(fabs(b.i)==0.0){
c.r= a.r/b.r;
c.i= a.i/b.r;
}else{
if(fabs(b.r)>=fabs(b.i)){
r= b.i/b.r;
den= b.r+r*b.i;
c.r= (a.r+r*a.i)/den;
c.i= (a.i-r*a.r)/den;
}else{
r= b.r/b.i;
den= b.i+r*b.r;
c.r= (a.r*r+a.i)/den;
c.i= (a.i*r-a.r)/den;
}
}
}else{
fprintf(stderr,"Error in cdiv(): Division by zero!\n");
exit(FAILURE);
}
return c;
}

/*:105*/
#line 5770 "./magbragg.w"

/*106:*/
#line 5811 "./magbragg.w"

dcomplex crdiv(dcomplex a,double b){
dcomplex c;
if(fabs(b)> 0.0){
c.r= a.r/b;
c.i= a.i/b;
}else{
fprintf(stderr,"Error in crdiv(): Division by zero!\n");
exit(FAILURE);
}
return c;
}

/*:106*/
#line 5771 "./magbragg.w"


/*:104*/
#line 5634 "./magbragg.w"

/*107:*/
#line 5827 "./magbragg.w"

dcomplex csqrt(dcomplex z){
dcomplex c;
double x,y,w,r;
if((z.r==0.0)&&(z.i==0.0)){
c.r= 0.0;
c.i= 0.0;
return c;
}else{
x= fabs(z.r);
y= fabs(z.i);
if(x>=y){
r= y/x;
w= sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
}else{
r= x/y;
w= sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
}
if(z.r>=0.0){
c.r= w;
c.i= 0.5*z.i/w;
}else{
c.i= ((z.i>=0)?w:-w);
c.r= 0.5*z.i/c.i;
}
return c;
}
}

/*:107*/
#line 5635 "./magbragg.w"

/*108:*/
#line 5862 "./magbragg.w"

/*109:*/
#line 5870 "./magbragg.w"

dcomplex cexp(dcomplex a){
dcomplex c;
double tmp= exp(a.r);
c.r= tmp*cos(a.i);
c.i= tmp*sin(a.i);
return c;
}

/*:109*/
#line 5863 "./magbragg.w"

/*110:*/
#line 5882 "./magbragg.w"

dcomplex crexpi(double a){
dcomplex c;
c.r= cos(a);
c.i= sin(a);
return c;
}

/*:110*/
#line 5864 "./magbragg.w"


/*:108*/
#line 5636 "./magbragg.w"


/*:94*/
#line 3487 "./magbragg.w"

/*130:*/
#line 7271 "./magbragg.w"

/*131:*/
#line 7282 "./magbragg.w"

void hl(char firststring[],char secondstring[]){
if(strlen(firststring)> 25){
fprintf(stderr,"%s:******* Error in hl() routine! *******\n",progname);
fprintf(stderr,"%s: The first string argument is too long:\n",progname);
fprintf(stderr,"%s:   '%s'\n",progname,firststring);
fprintf(stderr,"%s: String lengths is %d characters\n",
progname,(int)strlen(firststring));
fprintf(stderr,"%s: (Maximum 25 characters for first argument.)\n",
progname);
fprintf(stderr,"%s:******** End of error message *********\n",progname);
exit(FAILURE);
}
if(strlen(secondstring)> 55){
fprintf(stderr,"%s:******* Error in hl() routine! *******\n",progname);
fprintf(stderr,"%s: The second string argument is too long:\n",progname);
fprintf(stderr,"%s:   '%s'\n",progname,secondstring);
fprintf(stderr,"%s: String lengths is %d characters\n",
progname,(int)strlen(secondstring));
fprintf(stderr,"%s: (Maximum 55 characters for second argument.)\n",
progname);
fprintf(stderr,"%s:******** End of error message *********\n",progname);
exit(FAILURE);
}
fprintf(stderr,"%-25.25s%1.55s\n",firststring,secondstring);
}

/*:131*/
#line 7272 "./magbragg.w"

/*132:*/
#line 7313 "./magbragg.w"

void fhl(char linestring[]){
if(strlen(linestring)> 80){
fprintf(stderr,"%s:******* Error in fhl() routine! *******\n",progname);
fprintf(stderr,"%s: The following help line is too long:\n",progname);
fprintf(stderr,"%s:   '%s'\n",progname,linestring);
fprintf(stderr,"%s: String is %d characters\n",
progname,(int)strlen(linestring));
fprintf(stderr,"%s: (Maximum 80 characters per help line is allowed.)\n",
progname);
fprintf(stderr,"%s:******** End of error message *********\n",progname);
exit(FAILURE);
}
fprintf(stderr,"%s\n",linestring);
}

/*:132*/
#line 7273 "./magbragg.w"

/*133:*/
#line 7332 "./magbragg.w"

void showsomehelp(void){
fprintf(stderr," Usage: %s [options]\n",progname);
fhl(" Options:");
hl("  -h, --help","Display this help message and exit clean.");
hl("  -N <int>","");
hl("  -M <int>","");
hl("  -v, --verbose","");
hl("  -o, --outputfile <str>","");
fhl("  --fieldevolution {efield|stoke} <n> <str>");
fhl("  --intensityevolution <lambda> <str>");
fhl("  --normalize_length_to_um");
hl("  --normalize_intensity","");
hl("","When saving the spatial evolution of the intra-");
hl("","grating intensity, normalize the intensity with");
hl("","respect to the intensity at z=0, inside the");
hl("","grating (that is to say, normalize with respect");
hl("","to the initial intra-grating intensity). This");
hl("","option *only* affects the fields saved with the");
hl("","--intensityevolution or --fieldevolution");
hl("","options.");
hl("  -r, --random","");
hl("  -a, --apodize <real>","");
hl("","Apodize the grating structure over geometrical");
hl("","distance <real> at each end of the grating.");
hl("","The option only applies to gratings with sinus-");
hl("","oidal modulation of refractive index and");
hl("","gyration coefficient, in constant or chirped");
hl("","periodic configurations, as specified with the");
hl("","'--grating sinusoidal' or '--grating chirped'");
hl("","options respectively.");
hl("","   The apodization is applied to the refractive");
hl("","index modulation and, in cases where the linear");
hl("","magneto-optical is spatially modulated as well,");
hl("","to the linear gyration coefficient.");
fhl("  -j, --phasejump <r1 (angle)> <r2 (position)>");
hl("","Apply discrete phase jump in the spatial phase");
hl("","of the grating profile. This option adds the");
hl("","real number <r1> to the argument of the sinus-");
hl("","oidal function for the grating profile for all");
hl("","spatial coordinates z >= <r2>.");
hl("","   As in the case of apodization, this option");
hl("","only applies to gratings with sinusoidal modu-");
hl("","lation of refractive index and gyration coeffi-");
hl("","cient, in constant or chirped periodic configu-");
hl("","rations, as specified with the '--grating");
hl("","sinusoidal' or '--grating chirped' options");
hl("","respectively. The discrete phase jump applies");
hl("","to the linear linear refractive index modula-");
hl("","tion and, in cases where the linear magneto-");
hl("","optical interaction is spatially modulated as");
hl("","well, also to the linear gyration coefficient.");
fhl("  -w, --writegratingfile <str>");
hl("  --spectrumfile <str>","");
hl("","Generates the complex reflectance as function");
hl("","of the vacuum wavelength in meters, and save");
hl("","the spectrum in file named according to the");
hl("","supplied character string <str>.");
fhl("  --intensityspectrumfile <str>");
hl("","Generates the intensity reflectance as function");
hl("","of the vacuum wavelength in meters, and save");
hl("","the spectrum in file named according to the");
hl("","supplied character string <str>.");
fhl("  -g, --grating <grating options>");
hl("","Specifies the grating type, where");
hl("","<grating options> = ");
hl("","   [stepwise <stepwise options> |");
hl("","    sinusoidal <sinusoidal options> |");
hl("","    chirped <chirped options>]");
hl("","<stepwise options> = ");
hl("","   twolevel t1  <f> t2 <f>");
hl("","            n1  <f> n2 <f>  g1  <f> g2  <f>");
hl("","            pe1 <f> pe2 <f> pm1 <f> pm2 <f>");
hl("","            qe1 <f> qe2 <f> qm1 <f> qm2 <f>");
hl("","<sinusoidal options> = ");
hl("","            n  <n0>  <dn>  <nper>");
hl("","            g  <g0>  <dg>  <gper>");
hl("","            pe <pe0> <dpe> <peper>");
hl("","            pm <pm0> <dpm> <pmper>");
hl("","            qe <qe0> <dqe> <qeper>");
hl("","            qm <qm0> <dqm> <qmper>");
hl("","<chirped options> = ");
hl("","            n  <n0>  <dn>  <nper>  <ncrp>");
hl("","            g  <g0>  <dg>  <gper>  <gcrp>");
hl("","            pe <pe0> <dpe> <peper> <pecrp>");
hl("","            pm <pm0> <dpm> <pmper> <pmcrp>");
hl("","            qe <qe0> <dqe> <qeper> <qecrp>");
hl("","            qm <qm0> <dqm> <qmper> <qmcrp>");
hl("  -L,--gratinglength <f>","Physical length of grating in meter [m]");
fhl("     --refindsurr <f>");
fhl("     --trmtraject <str>");
fhl("     --trmintensity <istart> <istop> <mmi>");
fhl("   (intensity measured in Watts per square meter)");
fhl("     --trmellipticity <estart> <estop> <mme>");
fhl("     --lambdastart <lambda>");
fhl("   (start vacuum wavelength measured in meter)");
fhl("     --lambdastop <lambda>");
fhl("   (stop vacuum wavelength measured in meter)");
exit(FAILURE);
}

/*:133*/
#line 7274 "./magbragg.w"


/*:130*/
#line 3488 "./magbragg.w"


/*:50*/
#line 3409 "./magbragg.w"


int main(int argc,char*argv[])
{
/*51:*/
#line 3501 "./magbragg.w"

/*52:*/
#line 3565 "./magbragg.w"

static double pi= 3.1415926535897932384626433832795;
static double twopi= 2.0*3.1415926535897932384626433832795;
static double c= 2.997925e8;
static double epsilon0= 8.854e-12;

/*:52*/
#line 3502 "./magbragg.w"

/*53:*/
#line 3580 "./magbragg.w"

time_t initime= time(NULL),now= time(NULL),eta= time(NULL);

/*:53*/
#line 3503 "./magbragg.w"

/*54:*/
#line 3607 "./magbragg.w"

dcomplex*efp,*efm,*ebp,*ebm;

/*:54*/
#line 3504 "./magbragg.w"

/*55:*/
#line 3618 "./magbragg.w"

short verbose;

short randomdistribution;

short writegratingtofile;

short scale_stokesparams;

short normalize_length_to_micrometer;
short normalize_intensity;
short normalize_ellipticity;
short normalize_internally;
short odd_layer;
short chirpflag;

short apodize;
short phasejump;

short fieldevoflag;
short fieldevoflag_efield;

short intensityevoflag;

short fieldevoflag_stoke;

short intensityinfo;


short saveintensityinfologfile;

short trmtraject_specified;
short save_dbspectra;

short stokes_parameter_spectrum;

short display_surrounding_media;

short perturbed_gyration_constant;

/*:55*/
#line 3505 "./magbragg.w"

/*56:*/
#line 3671 "./magbragg.w"

long mm,mme,mmi;

/*:56*/
#line 3506 "./magbragg.w"

/*57:*/
#line 3724 "./magbragg.w"

long nn,modnum;
double ll,*z,*dz,*n,*g,*pe,*pm,*qe,*qm,*etafp,*etafm,*etabp,*etabm;
double*taup,*taum,*taupp,*taupm,*rhop,*rhom,*rhopp,*rhopm;

/*:57*/
#line 3507 "./magbragg.w"

/*58:*/
#line 3731 "./magbragg.w"

FILE*fp_s0,*fp_s1,*fp_s2,*fp_s3;
FILE*fp_v0,*fp_v1,*fp_v2,*fp_v3;
FILE*fp_w0,*fp_w1,*fp_w2,*fp_w3;
FILE*fp_evo= NULL,*fp_ievo= NULL,*fp_spec= NULL,*fp_traject= NULL;
FILE*fp_irspec= NULL,*fp_itspec= NULL,*fp_icspec= NULL,*fp_gr= NULL;
FILE*fp_evo_s0= NULL,*fp_evo_s1= NULL,*fp_evo_s2= NULL,*fp_evo_s3= NULL;
FILE*intensinfologfile= NULL;

/*:58*/
#line 3508 "./magbragg.w"

/*59:*/
#line 3745 "./magbragg.w"

char gratingtype[NCHMAX],gratingsubtype[NCHMAX];
char fieldevofilename[NCHMAX],intensityevofilename[NCHMAX];
char spectrumfilename[NCHMAX],trmtraject_filename[NCHMAX];
char intensity_reflection_spectrumfilename[NCHMAX],
intensity_transmission_spectrumfilename[NCHMAX],
intensity_check_spectrumfilename[NCHMAX];
char fieldevofilename_s0[NCHMAX],fieldevofilename_s1[NCHMAX],
fieldevofilename_s2[NCHMAX],fieldevofilename_s3[NCHMAX];
char intensinfologfilename[NCHMAX];
char outfilename[NCHMAX],gratingfilename[NCHMAX];
char outfilename_s0[NCHMAX],outfilename_s1[NCHMAX],
outfilename_s2[NCHMAX],outfilename_s3[NCHMAX];
char outfilename_v0[NCHMAX],outfilename_v1[NCHMAX],
outfilename_v2[NCHMAX],outfilename_v3[NCHMAX];
char outfilename_w0[NCHMAX],outfilename_w1[NCHMAX],
outfilename_w2[NCHMAX],outfilename_w3[NCHMAX];

/*:59*/
#line 3509 "./magbragg.w"

/*60:*/
#line 3770 "./magbragg.w"

dcomplex tmpfp,tmpfm,tmpbp,tmpbm;
double tmp,zt,s0,s1,s2,s3,v0,v1,v2,v3,w0,w1,w2,w3,stn,
maxintens= 0.0,maxintens_inintens= 0.0,maxintens_inellip= 0.0,
maxintens_trintens= 0.0,maxintens_trellip= 0.0;
int no_arg,status= 0,tmpch;
long j,k,ke,ki;

/*:60*/
#line 3510 "./magbragg.w"

long nne,jje,mmtraject;
long ranseed;
long maxintens_layer= 0;
int fractal_level= 0;
int maximum_fractal_level= 0;

double*w0traj= NULL,*w3traj= NULL;
double lambdastart;

double lambdastop;

double lambda;

double omega;

double ievolambda;
double apolength;

double phasejumpangle;

double phasejumpposition;

double phi;

double gyroperturb_position;


double gyroperturb_amplitude;


double gyroperturb_width;


double nsurr;

double n1,n2,t1,t2,nper,ncrp,g1,g2,gper,gcrp,
pe1,pe2,peper,pecrp,pm1,pm2,pmper,pmcrp,
qe1,qe2,qeper,qecrp,qm1,qm2,qmper,qmcrp;
double modn1,modt1,modg1,modpe1,modpm1,modqe1,modqm1,
aafp2,aafm2,aabp2,aabm2;
double trmintensity,trmintenstart,trmintenstop,
trmellipticity,trmellipstart,trmellipstop,stoke_scalefactor;
static double nndef= 600;
static double mmdef= 1000;
static double mmedef= 1;
static double mmidef= 1;
static double lldef= 0.050;
static double nsurrdef= 1.0;
static double lambdastartdef= 1536.0e-9,lambdastopdef= 1546.0e-9;

/*:51*/
#line 3413 "./magbragg.w"

/*61:*/
#line 3780 "./magbragg.w"

trmtraject_specified= 0;
intensityinfo= 0;
saveintensityinfologfile= 0;
randomdistribution= 0;
writegratingtofile= 0;
scale_stokesparams= 0;
stoke_scalefactor= 1.0;
normalize_length_to_micrometer= 0;
normalize_intensity= 0;
normalize_ellipticity= 0;
ranseed= 1097;
modnum= -1;
verbose= 0;
apodize= 0;
phasejump= 0;
normalize_internally= 0;
odd_layer= 0;
chirpflag= 1;
save_dbspectra= 0;
stokes_parameter_spectrum= 0;
display_surrounding_media= 1;


perturbed_gyration_constant= 0;
fieldevoflag= 0;
fieldevoflag_efield= 0;
fieldevoflag_stoke= 0;
intensityevoflag= 0;
mm= mmdef;
nn= nndef;
mme= mmedef;
mmi= mmidef;
nne= 1;
ll= lldef;
lambdastart= lambdastartdef;
lambdastop= lambdastopdef;
phasejumpangle= 0.0;
phasejumpposition= 0.0;
phi= 0.0;
nsurr= nsurrdef;
strcpy(outfilename,"out.stok");
strcpy(fieldevofilename,"out.fevo.dat");
strcpy(fieldevofilename_s0,"out.fevo.s0.dat");
strcpy(fieldevofilename_s1,"out.fevo.s1.dat");
strcpy(fieldevofilename_s2,"out.fevo.s2.dat");
strcpy(fieldevofilename_s3,"out.fevo.s3.dat");
strcpy(spectrumfilename,"out.rsp.dat");
strcpy(intensity_reflection_spectrumfilename,"out.irsp.dat");
strcpy(intensity_transmission_spectrumfilename,"out.trsp.dat");
strcpy(intensity_check_spectrumfilename,"out.chec.dat");
fp_s0= NULL;fp_s1= NULL;fp_s2= NULL;fp_s3= NULL;
fp_v0= NULL;fp_v1= NULL;fp_v2= NULL;fp_v3= NULL;
fp_w0= NULL;fp_w1= NULL;fp_w2= NULL;fp_w3= NULL;

/*:61*/
#line 3414 "./magbragg.w"

/*116:*/
#line 5955 "./magbragg.w"

{
progname= strip_away_path(argv[0]);
fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
no_arg= argc;
while(--argc){
if(!strcmp(argv[no_arg-argc],"-o")||
!strcmp(argv[no_arg-argc],"--outputfile")){
--argc;
strcpy(outfilename,argv[no_arg-argc]);
}else if(!strcmp(argv[no_arg-argc],"-w")||
!strcmp(argv[no_arg-argc],"--writegratingfile")){
--argc;
strcpy(gratingfilename,argv[no_arg-argc]);
writegratingtofile= 1;
}else if(!strcmp(argv[no_arg-argc],"--displaysurrmedia")){
display_surrounding_media= (display_surrounding_media?0:1);
}else if(!strcmp(argv[no_arg-argc],"-a")||
!strcmp(argv[no_arg-argc],"--apodize")){
/*118:*/
#line 6144 "./magbragg.w"

{
apodize= 1;
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&apolength)){
fprintf(stderr,"%s: Error in '--apodize' option.\n",progname);
exit(FAILURE);
}
}

/*:118*/
#line 5974 "./magbragg.w"
;
}else if(!strcmp(argv[no_arg-argc],"--phasejump")){
/*119:*/
#line 6157 "./magbragg.w"

{
phasejump= 1;
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&phasejumpangle)){
fprintf(stderr,"%s: Error in 1st arg of '--phasejump' option.\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&phasejumpposition)){
fprintf(stderr,"%s: Error in 2nd arg of '--phasejump' option.\n",
progname);
exit(FAILURE);
}
}

/*:119*/
#line 5976 "./magbragg.w"
;
}else if(!strcmp(argv[no_arg-argc],"--fieldevolution")){
/*121:*/
#line 6205 "./magbragg.w"

{
--argc;
if(!strcmp(argv[no_arg-argc],"efield")){
fieldevoflag_efield= 1;
}else if(!strcmp(argv[no_arg-argc],"stoke")){
fieldevoflag_stoke= 1;
}else{
fprintf(stderr,
"%s: Unknown field evolution flag '%s' in second argument of\n"
" --fieldevolution option.\n",progname,argv[no_arg-argc]);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&nne)){
fprintf(stderr,"%s: Error in '--fieldevolution' option.\n",
progname);
exit(FAILURE);
}
--argc;
strcpy(fieldevofilename,argv[no_arg-argc]);
fieldevoflag= 1;
if(fieldevoflag_stoke){
sprintf(fieldevofilename_s0,"%s.s0.dat",fieldevofilename);
sprintf(fieldevofilename_s1,"%s.s1.dat",fieldevofilename);
sprintf(fieldevofilename_s2,"%s.s2.dat",fieldevofilename);
sprintf(fieldevofilename_s3,"%s.s3.dat",fieldevofilename);
}
}

/*:121*/
#line 5978 "./magbragg.w"
;
}else if(!strcmp(argv[no_arg-argc],"--intensityevolution")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&ievolambda)){
fprintf(stderr,"%s: Error in '--intensityevolution' option.\n",
progname);
exit(FAILURE);
}
--argc;
strcpy(intensityevofilename,argv[no_arg-argc]);
intensityevoflag= 1;
}else if(!strcmp(argv[no_arg-argc],"--spectrumfile")){
--argc;
strcpy(spectrumfilename,argv[no_arg-argc]);
}else if(!strcmp(argv[no_arg-argc],"--stokesspectrum")){
stokes_parameter_spectrum= 1;
}else if(!strcmp(argv[no_arg-argc],"--trmtraject")){
--argc;
strcpy(trmtraject_filename,argv[no_arg-argc]);
trmtraject_specified= 1;
}else if(!strcmp(argv[no_arg-argc],"--intensityspectrumfile")){
--argc;
strcpy(intensity_reflection_spectrumfilename,argv[no_arg-argc]);
strcpy(intensity_transmission_spectrumfilename,argv[no_arg-argc]);
strcpy(intensity_check_spectrumfilename,argv[no_arg-argc]);
sprintf(intensity_reflection_spectrumfilename,"%s.irsp.dat",
intensity_reflection_spectrumfilename);
sprintf(intensity_transmission_spectrumfilename,"%s.itsp.dat",
intensity_transmission_spectrumfilename);
sprintf(intensity_check_spectrumfilename,"%s.chec.dat",
intensity_check_spectrumfilename);

}else if(!strcmp(argv[no_arg-argc],"--logarithmicspectrum")){
save_dbspectra= 1;

}else if(!strcmp(argv[no_arg-argc],"-v")||
!strcmp(argv[no_arg-argc],"--verbose")){
verbose= (verbose?0:1);

}else if(!strcmp(argv[no_arg-argc],"--scale_stokesparams")){
scale_stokesparams= 1;
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&stoke_scalefactor)){
fprintf(stderr,"%s: Error in '--scale_stokesparams' option.\n",\
progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"--intensityinfo")){
intensityinfo= (intensityinfo?0:1);
}else if(!strcmp(argv[no_arg-argc],"--intensityinfologfile")){
intensityinfo= 1;
saveintensityinfologfile= 1;
--argc;
strcpy(intensinfologfilename,argv[no_arg-argc]);
}else if(!strcmp(argv[no_arg-argc],"--gyroperturb")){
/*120:*/
#line 6177 "./magbragg.w"

{
perturbed_gyration_constant= 1;
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&gyroperturb_position)){
fprintf(stderr,"%s: Error in '--gyroperturb' option.\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&gyroperturb_amplitude)){
fprintf(stderr,"%s: Error in '--gyroperturb' option.\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&gyroperturb_width)){
fprintf(stderr,"%s: Error in '--gyroperturb' option.\n",progname);
exit(FAILURE);
}
}

/*:120*/
#line 6033 "./magbragg.w"
;
}else if(!strcmp(argv[no_arg-argc],"--normalize_length_to_um")){
normalize_length_to_micrometer= (normalize_length_to_micrometer?0:1);
}else if(!strcmp(argv[no_arg-argc],"--normalize_intensity")){
normalize_intensity= (normalize_intensity?0:1);
}else if(!strcmp(argv[no_arg-argc],"--normalize_ellipticity")){
normalize_ellipticity= (normalize_ellipticity?0:1);
}else if(!strcmp(argv[no_arg-argc],"--normalizedrepresentation")){
normalize_internally= (normalize_internally?0:1);
}else if(!strcmp(argv[no_arg-argc],"-r")||
!strcmp(argv[no_arg-argc],"--random")){
randomdistribution= (randomdistribution?0:1);
}else if(!strcmp(argv[no_arg-argc],"-h")||
!strcmp(argv[no_arg-argc],"--help")){
showsomehelp();
}else if(!strcmp(argv[no_arg-argc],"-g")||
!strcmp(argv[no_arg-argc],"--grating")){
--argc;
if(!strcmp(argv[no_arg-argc],"stepwise")){
/*122:*/
#line 6238 "./magbragg.w"

{
strcpy(gratingtype,argv[no_arg-argc]);
--argc;
if(!strcmp(argv[no_arg-argc],"twolevel")){
strcpy(gratingsubtype,argv[no_arg-argc]);
--argc;
if(strcmp(argv[no_arg-argc],"t1")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 't1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&t1)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for t1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"t2")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 't2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&t2)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for t2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"n1")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'n1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"n2")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'n2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"g1")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'g1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"g2")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'g2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pe1")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pe1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pe2")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pe2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pm1")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pm1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm1).\n",\
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pm2")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pm2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm2).\n",\
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qe1")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qe1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe1).\n",\
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qe2")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qe2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qm1")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qm1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qm2")){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qm2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
fprintf(stderr,"%s: Error in 'twolevel' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm2).\n",progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Error in '-g' or '--grating' option.\n",progname);
fprintf(stderr,"%s: (No valid stepwise grating type found!)\n",progname);
exit(FAILURE);
}
}

/*:122*/
#line 6052 "./magbragg.w"

}else if(!strcmp(argv[no_arg-argc],"sinusoidal")){
/*123:*/
#line 6425 "./magbragg.w"

{
strcpy(gratingtype,argv[no_arg-argc]);
--argc;
if(strcmp(argv[no_arg-argc],"n")){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'n').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&nper)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"g")){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'g').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&gper)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pe")){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pe').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&peper)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pm")){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pm').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pmper)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qe")){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qe').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qeper)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qm")){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qm').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qmper)){
fprintf(stderr,"%s: Error in 'sinusoidal' option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm [3rd arg]).\n",progname);
exit(FAILURE);
}
}

/*:123*/
#line 6054 "./magbragg.w"

}else if(!strcmp(argv[no_arg-argc],"chirped")){
/*124:*/
#line 6577 "./magbragg.w"

{
strcpy(gratingtype,argv[no_arg-argc]);
--argc;
if(strcmp(argv[no_arg-argc],"n")){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'n').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&nper)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&ncrp)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [4th arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"g")){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'g').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&gper)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&gcrp)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [4th arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pe")){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pe').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&peper)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pecrp)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe [4th arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pm")){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pm').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pmper)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pmcrp)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm [4th arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qe")){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qe').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qeper)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qecrp)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe [4th arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qm")){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qm').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm [1st arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm [2nd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qmper)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm [3rd arg]).\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qmcrp)){
fprintf(stderr,"%s: Error in 'chirped' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm [4th arg]).\n",progname);
exit(FAILURE);
}
}

/*:124*/
#line 6056 "./magbragg.w"

}else if(!strcmp(argv[no_arg-argc],"fractal")){
/*125:*/
#line 6766 "./magbragg.w"

{
strcpy(gratingtype,argv[no_arg-argc]);
--argc;
if(!strcmp(argv[no_arg-argc],"cantor")){
strcpy(gratingsubtype,argv[no_arg-argc]);
--argc;
if(strcmp(argv[no_arg-argc],"fractal_level")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'fractal_level').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%d",&fractal_level)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for fractal_level).\n",
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"maximum_fractal_level")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'maximum_fractal_level').\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%d",&maximum_fractal_level)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for maximum_fractal_level).\n",
progname);
exit(FAILURE);
}
nn= 1;
for(j= 1;j<=fractal_level;j++)nn= 2*nn;
--argc;
if(strcmp(argv[no_arg-argc],"t1")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 't1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&t1)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for t1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"t2")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 't2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&t2)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for t2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"n1")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'n1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n1)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"n2")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'n2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&n2)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for n2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"g1")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'g1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g1)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"g2")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'g2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&g2)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for g2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pe1")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pe1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe1)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pe2")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pe2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pe2)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pm1")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pm1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm1)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm1).\n",\
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pm2")){
fprintf(stderr,"%s: Error.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pm2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&pm2)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm2).\n",\
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qe1")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qe1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe1)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe1).\n",\
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qe2")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qe2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qe2)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe2).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qm1")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qm1').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm1)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm1).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qm2")){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qm2').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&qm2)){
fprintf(stderr,"%s: Error in 'cantor' grating option.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm2).\n",progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Error in 'fractal' grating option.\n",progname);
fprintf(stderr,"%s: (No valid fractal type found!)\n",progname);
fprintf(stderr,"%s: (Currently only Cantor type implemented)\n",progname);
exit(FAILURE);
}
}

/*:125*/
#line 6058 "./magbragg.w"

}else{
fprintf(stderr,"%s: Error in '-g' or '--grating' option.\n",
progname);
fprintf(stderr,"%s: (No valid grating type found!)\n",progname);
exit(FAILURE);
}
}else if((!strcmp(argv[no_arg-argc],"-L"))||
(!strcmp(argv[no_arg-argc],"--gratinglength"))){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&ll)){
fprintf(stderr,"%s: Error in '-L' option.\n",progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"--modifylayer")){
/*126:*/
#line 6984 "./magbragg.w"

{
--argc;
if(strcmp(argv[no_arg-argc],"num")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 'num' after '--modifylayer').\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&modnum)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for num).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"t")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 't' [for layer thickness]).\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&modt1)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for t [layer thickness]).\n",
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"n")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 'n' [for refractive index]).\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&modn1)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for n [refractive index]).\n",
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"g")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 'g' [for gyration constant]).\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&modg1)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for g [gyration constant]).\n",
progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pe")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pe').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&modpe1)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for pe).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"pm")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 'pm').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&modpm1)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for pm).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qe")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qe').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&modqe1)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for qe).\n",progname);
exit(FAILURE);
}
--argc;
if(strcmp(argv[no_arg-argc],"qm")){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Expecting string 'qm').\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&modqm1)){
fprintf(stderr,"%s: Error in '--modifylayer'.\n",progname);
fprintf(stderr,"%s: (Could not read data for qm).\n",progname);
exit(FAILURE);
}
}

/*:126*/
#line 6073 "./magbragg.w"

}else if(!strcmp(argv[no_arg-argc],"--refindsurr")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&nsurr)){
fprintf(stderr,"%s: Error in '--refindsurr' option.\n",progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"--trmintensity")){
/*127:*/
#line 7103 "./magbragg.w"

{
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&trmintenstart)){
fprintf(stderr,"%s: Error in '--trmintensity' option.\n",\
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&trmintenstop)){
fprintf(stderr,"%s: Error in '--trmintensity' option.\n",\
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&mmi)){
fprintf(stderr,"%s: Error in '--trmintensity' option.\n",\
progname);
exit(FAILURE);
}
}

/*:127*/
#line 6081 "./magbragg.w"

}else if(!strcmp(argv[no_arg-argc],"--trmellipticity")){
/*128:*/
#line 7137 "./magbragg.w"

{
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&trmellipstart)){
fprintf(stderr,"%s: Error in '--trmellipticity' option.\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&trmellipstop)){
fprintf(stderr,"%s: Error in '--trmellipticity' option.\n",progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&mme)){
fprintf(stderr,"%s: Error in '--trmellipticity' option.\n",progname);
exit(FAILURE);
}
}

/*:128*/
#line 6083 "./magbragg.w"

}else if(!strcmp(argv[no_arg-argc],"--lambdastart")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&lambdastart)){
fprintf(stderr,"%s: Error in '--lambdastart' option.\n",progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"--lambdastop")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&lambdastop)){
fprintf(stderr,"%s: Error in '--lambdastop' option.\n",progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"--wlenlinz")){
chirpflag= 1;
}else if(!strcmp(argv[no_arg-argc],"--freqlinz")){
chirpflag= 0;
}else if(!strcmp(argv[no_arg-argc],"-N")){
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&nn)){
fprintf(stderr,"%s: Error in '-N' option.\n",progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"-M")){
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&mm)){
fprintf(stderr,"%s: Error in '-M' option.\n",progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Specified option invalid!\n",progname);
showsomehelp();
exit(FAILURE);
}
}
/*117:*/
#line 6125 "./magbragg.w"

{
sprintf(outfilename_s0,"%s.s0.dat",outfilename);
sprintf(outfilename_s1,"%s.s1.dat",outfilename);
sprintf(outfilename_s2,"%s.s2.dat",outfilename);
sprintf(outfilename_s3,"%s.s3.dat",outfilename);
sprintf(outfilename_v0,"%s.v0.dat",outfilename);
sprintf(outfilename_v1,"%s.v1.dat",outfilename);
sprintf(outfilename_v2,"%s.v2.dat",outfilename);
sprintf(outfilename_v3,"%s.v3.dat",outfilename);
sprintf(outfilename_w0,"%s.w0.dat",outfilename);
sprintf(outfilename_w1,"%s.w1.dat",outfilename);
sprintf(outfilename_w2,"%s.w2.dat",outfilename);
sprintf(outfilename_w3,"%s.w3.dat",outfilename);
}

/*:117*/
#line 6118 "./magbragg.w"

/*129:*/
#line 7161 "./magbragg.w"

{
if(verbose){
for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
fprintf(stdout,"Input parameters:\n");
fprintf(stdout,
"Grating type: %s\n",gratingtype);
if(!strcmp(gratingtype,"sinusoidal")){
fprintf(stdout,
"Bias refractive index:          %10.6e\n",n1);
fprintf(stdout,
"Refractive index modulation:    %10.6e\n",n2);
fprintf(stdout,
"Modulation period:              %10.6e [m]\n",nper);
}else if(!strcmp(gratingtype,"chirped")){
}else if(!strcmp(gratingtype,"stepwise")){
}else if(!strcmp(gratingtype,"fractal")){
}else{
fprintf(stdout,"%s: Error: Unknown grating type '%s'\n",
progname,gratingtype);
}
fprintf(stdout,"Geometrical length:            "
"L=%10.6e [m]\n",ll);
fprintf(stdout,"Surrounding refractive index:  "
"nsurr=%10.6e\n",nsurr);
fprintf(stdout,"Begin wavelength of spectrum:  "
"lambda_start=%10.6e [m]\n",lambdastart);
fprintf(stdout,"End wavelength of spectrum:    "
"lambda_stop=%10.6e [m]\n",lambdastop);
fprintf(stdout,"Number of samples in spectrum: "
"M=%-12ld\n",mm);
fprintf(stdout,"Number of discrete layers:     "
"N=%-12ld\n",nn);
if(trmtraject_specified){
fprintf(stdout,"Trajectory for transmitted Stokes parameters: %s\n",
trmtraject_filename);
}else{
fprintf(stdout,"Number of samples in output intensity:   "
"mmi=%-12ld\n",mmi);
fprintf(stdout,"Number of samples in output ellipticity: "
"mme=%-12ld\n",mme);
}
fprintf(stdout,"Stokes parameters will be written to files:\n");
fprintf(stdout,"   %s  [S0 (incident wave)],\n",outfilename_s0);
fprintf(stdout,"   %s  [S1 (incident wave)],\n",outfilename_s1);
fprintf(stdout,"   %s  [S2 (incident wave)],\n",outfilename_s2);
fprintf(stdout,"   %s  [S3 (incident wave)],\n",outfilename_s3);
fprintf(stdout,"   %s  [V0 (reflected wave)],\n",outfilename_v0);
fprintf(stdout,"   %s  [V1 (reflected wave)],\n",outfilename_v1);
fprintf(stdout,"   %s  [V2 (reflected wave)],\n",outfilename_v2);
fprintf(stdout,"   %s  [V3 (reflected wave)],\n",outfilename_v3);
fprintf(stdout,"   %s  [W0 (transmitted wave)],\n",outfilename_w0);
fprintf(stdout,"   %s  [W1 (transmitted wave)],\n",outfilename_w1);
fprintf(stdout,"   %s  [W2 (transmitted wave)],\n",outfilename_w2);
fprintf(stdout,"   %s  [W3 (transmitted wave)],\n",outfilename_w3);
if(fieldevoflag){
if(fieldevoflag_efield){
if(strcmp(fieldevofilename,"")){
fprintf(stdout,"Intra grating optical field evolution will "
"be written to file:\n");
fprintf(stdout,"   %s\n",fieldevofilename);
}else{
fprintf(stderr,
"%s: Error: No file name specified for saving spatial\n"
"field evolution. (efield option)\n",progname);
exit(FAILURE);
}
fprintf(stdout,
"(Intra grating field evolution will be presented in terms of\n"
"the electrical field displacement.)\n");
}else if(fieldevoflag_stoke){
if(strcmp(fieldevofilename_s0,"")
&&strcmp(fieldevofilename_s1,"")
&&strcmp(fieldevofilename_s2,"")
&&strcmp(fieldevofilename_s3,"")){
fprintf(stdout,"Intra grating optical field evolution will "
"be written to files:\n");
fprintf(stdout,"   %s\n   %s\n   %s\n   %s\n",
fieldevofilename_s0,fieldevofilename_s1,
fieldevofilename_s2,fieldevofilename_s3);
}else{
fprintf(stderr,
"%s: Error: No file name specified for saving spatial\n"
"field evolution. (stoke option)\n",progname);
exit(FAILURE);
}
fprintf(stdout,
"(Intra grating field evolution will be presented in terms of\n"
"the Stokes parameters of the forward propagating field.)\n");
}else{
fprintf(stderr,"%s: Unknown field evolution flag.\n",progname);
exit(FAILURE);
}
fprintf(stdout,
"Number of intermediate samples within each layer: %-12ld\n",nne);
}
if(intensityevoflag){
if(strcmp(intensityevofilename,"")){
fprintf(stdout,"Intra grating optical intensity evolution will "
"be written to file:\n");
fprintf(stdout,"   %s\n",intensityevofilename);
}
}
fprintf(stdout,"Program execution started %s",ctime(&initime));
for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
}
}

/*:129*/
#line 6119 "./magbragg.w"

}

/*:116*/
#line 3415 "./magbragg.w"

/*134:*/
#line 7435 "./magbragg.w"

{
mmtraject= 0;
if(trmtraject_specified){
if((fp_traject= fopen(trmtraject_filename,"r"))==NULL){
fprintf(stderr,
"%s: Could not open file %s for reading Stokes parameter"
" trajectory of transmitted wave!\n",
progname,trmtraject_filename);
exit(FAILURE);
}
fseek(fp_traject,0L,SEEK_SET);


while((tmpch= getc(fp_traject))!=EOF){
ungetc(tmpch,fp_traject);
fscanf(fp_traject,"%lf",&tmp);
fscanf(fp_traject,"%lf",&tmp);
mmtraject++;


tmpch= getc(fp_traject);
while((tmpch!=EOF)&&(!numeric(tmpch))){
tmpch= getc(fp_traject);
}
if(tmpch!=EOF)ungetc(tmpch,fp_traject);
}

if(verbose){
fprintf(stdout,
"%s: I have now pre-parsed the specified trajectory of transmitted\n"
"Stokes parameters (W0,W3) in file %s, and I found %-ld points.\n",
progname,trmtraject_filename,mmtraject);
fprintf(stdout,
"%s: Now allocating the vectors for the transmitted trajectory...",
progname);
}

fseek(fp_traject,0L,SEEK_SET);
w0traj= dvector(1,mmtraject);
w3traj= dvector(1,mmtraject);
for(ke= 1;ke<=mmtraject;ke++){
fscanf(fp_traject,"%le",&w0traj[ke]);
fscanf(fp_traject,"%le",&w3traj[ke]);
}

if(0==1){
for(ke= 1;ke<=mmtraject;ke++)
fprintf(stdout,"w0=%e   w3=%e\n",w0traj[ke],w3traj[ke]);
}

fclose(fp_traject);
}
}

/*:134*/
#line 3416 "./magbragg.w"

/*135:*/
#line 7540 "./magbragg.w"

{
if((mme> 1)&&(mmi> 1)){
if((fp_s0= fopen(outfilename_s0,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
progname,outfilename_s0);
exit(FAILURE);
}
if((fp_s1= fopen(outfilename_s1,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
progname,outfilename_s1);
exit(FAILURE);
}
if((fp_s2= fopen(outfilename_s2,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
progname,outfilename_s2);
exit(FAILURE);
}
if((fp_s3= fopen(outfilename_s3,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving incident wave!\n",
progname,outfilename_s3);
exit(FAILURE);
}
if((fp_v0= fopen(outfilename_v0,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
progname,outfilename_v0);
exit(FAILURE);
}
if((fp_v1= fopen(outfilename_v1,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
progname,outfilename_v1);
exit(FAILURE);
}
if((fp_v2= fopen(outfilename_v2,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
progname,outfilename_v2);
exit(FAILURE);
}
if((fp_v3= fopen(outfilename_v3,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving reflected wave!\n",
progname,outfilename_v3);
exit(FAILURE);
}
if((fp_w0= fopen(outfilename_w0,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
progname,outfilename_w0);
exit(FAILURE);
}
if((fp_w1= fopen(outfilename_w1,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
progname,outfilename_w1);
exit(FAILURE);
}
if((fp_w2= fopen(outfilename_w2,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
progname,outfilename_w2);
exit(FAILURE);
}
if((fp_w3= fopen(outfilename_w3,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for saving transmitted wave!\n",
progname,outfilename_w3);
exit(FAILURE);
}
}
if((mme> 1)&&(mmi> 1)){
fseek(fp_s0,0L,SEEK_SET);
fseek(fp_s1,0L,SEEK_SET);
fseek(fp_s2,0L,SEEK_SET);
fseek(fp_s3,0L,SEEK_SET);
fseek(fp_v0,0L,SEEK_SET);
fseek(fp_v1,0L,SEEK_SET);
fseek(fp_v2,0L,SEEK_SET);
fseek(fp_v3,0L,SEEK_SET);
fseek(fp_w0,0L,SEEK_SET);
fseek(fp_w1,0L,SEEK_SET);
fseek(fp_w2,0L,SEEK_SET);
fseek(fp_w3,0L,SEEK_SET);
}
if(fieldevoflag){
if(fieldevoflag_efield){
if(strcmp(fieldevofilename,"")){
if((fp_evo= fopen(fieldevofilename,"w"))==NULL){
fprintf(stderr,
"%s: Could not open file %s for fieldevo output!\n",
progname,fieldevofilename);
exit(FAILURE);
}
}
}else if(fieldevoflag_stoke){
if(strcmp(fieldevofilename_s0,"")){
if((fp_evo_s0= fopen(fieldevofilename_s0,"w"))==NULL){
fprintf(stderr,
"%s: Could not open file %s for saving spatial S0"
" distribution!\n",progname,fieldevofilename);
exit(FAILURE);
}
}else{
fprintf(stderr,
"%s: A name for the file for saving spatial S0 distribution"
" is required!\n",progname);
exit(FAILURE);
}
if(strcmp(fieldevofilename_s1,"")){
if((fp_evo_s1= fopen(fieldevofilename_s1,"w"))==NULL){
fprintf(stderr,
"%s: Could not open file %s for saving spatial S1"
" distribution!\n",progname,fieldevofilename);
exit(FAILURE);
}
}else{
fprintf(stderr,
"%s: A name for the file for saving spatial S1 distribution"
" is required!\n",progname);
exit(FAILURE);
}
if(strcmp(fieldevofilename_s2,"")){
if((fp_evo_s2= fopen(fieldevofilename_s2,"w"))==NULL){
fprintf(stderr,
"%s: Could not open file %s for saving spatial S2"
" distribution!\n",progname,fieldevofilename);
exit(FAILURE);
}
}else{
fprintf(stderr,
"%s: A name for the file for saving spatial S2 distribution"
" is required!\n",progname);
exit(FAILURE);
}
if(strcmp(fieldevofilename_s3,"")){
if((fp_evo_s3= fopen(fieldevofilename_s3,"w"))==NULL){
fprintf(stderr,
"%s: Could not open file %s for saving spatial S3"
" distribution!\n",progname,fieldevofilename);
exit(FAILURE);
}
}else{
fprintf(stderr,
"%s: A name for the file for saving spatial S3 distribution"
" is required!\n",progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Unknown field evolution flag.\n",progname);
exit(FAILURE);
}
}
if(intensityevoflag){
if(strcmp(intensityevofilename,"")){
if((fp_ievo= fopen(intensityevofilename,"w"))==NULL){
fprintf(stderr,
"%s: Could not open file %s for intensityevo output!\n",
progname,intensityevofilename);
exit(FAILURE);
}
}
}
if(!((mme> 1)&&(mmi> 1))){
if((fp_spec= fopen(spectrumfilename,"w"))==NULL){
fprintf(stderr,"%s: Could not open file %s for spectrum output!\n",
progname,spectrumfilename);
exit(FAILURE);
}
if((fp_irspec= fopen(intensity_reflection_spectrumfilename,"w"))==NULL){
fprintf(stderr,
"%s: Could not open %s for intensity reflection spectrum!\n",
progname,intensity_reflection_spectrumfilename);
exit(FAILURE);
}
if((fp_itspec= fopen(intensity_transmission_spectrumfilename,"w"))==NULL){
fprintf(stderr,
"%s: Could not open %s for intensity transmission spectrum!\n",
progname,intensity_transmission_spectrumfilename);
exit(FAILURE);
}
if((fp_icspec= fopen(intensity_check_spectrumfilename,"w"))==NULL){
fprintf(stderr,
"%s: Could not open %s for checking spectra!\n",
progname,intensity_check_spectrumfilename);
exit(FAILURE);
}
}
}

/*:135*/
#line 3417 "./magbragg.w"

/*62:*/
#line 3911 "./magbragg.w"

{
z= dvector(1,nn);
dz= dvector(1,nn-1);
efp= dcvector(0,nn);
efm= dcvector(0,nn);
ebp= dcvector(0,nn);
ebm= dcvector(0,nn);
taup= dvector(1,nn);
taum= dvector(1,nn);
taupp= dvector(1,nn);
taupm= dvector(1,nn);
rhop= dvector(1,nn);
rhom= dvector(1,nn);
rhopp= dvector(1,nn);
rhopm= dvector(1,nn);
n= dvector(0,nn);
g= dvector(0,nn);
pe= dvector(1,nn-1);
pm= dvector(1,nn-1);
qe= dvector(1,nn-1);
qm= dvector(1,nn-1);
etafp= dvector(1,nn-1);
etabp= dvector(1,nn-1);
etafm= dvector(1,nn-1);
etabm= dvector(1,nn-1);
}

/*:62*/
#line 3418 "./magbragg.w"

/*64:*/
#line 3984 "./magbragg.w"

{
if(!strcmp(gratingtype,"stepwise")){
/*65:*/
#line 4050 "./magbragg.w"

{
if(!strcmp(gratingsubtype,"twolevel")){
if(randomdistribution){
for(j= 1;j<=nn-1;j++){
if(j==modnum){
if(j==1)z[j]= 0.0;
z[j+1]= z[j]+modt1;
dz[j]= modt1;
n[j]= modn1;
g[j]= modg1;
pe[j]= modpe1;
pm[j]= modpm1;
qe[j]= modqe1;
qm[j]= modqm1;
}else{
ranseed= ranseed+j;
if(ran1(&ranseed)> 0.5){
if(verbose)fprintf(stdout,"Random number: 1\n");
if(j==1)z[j]= 0.0;
z[j+1]= z[j]+t1;
dz[j]= t1;
n[j]= n1;
g[j]= g1;
pe[j]= pe1;
pm[j]= pm1;
qe[j]= qe1;
qm[j]= qm1;
}else{
if(verbose)fprintf(stdout,"Random number: 0\n");
if(j==1)z[j]= 0.0;
z[j+1]= z[j]+t2;
dz[j]= t2;
n[j]= n2;
g[j]= g2;
pe[j]= pe2;
pm[j]= pm2;
qe[j]= qe2;
qm[j]= qm2;
}
}
}
}else{
if(modnum<0){
for(j= 1;j<=nn-1;j= j+2){
z[j]= 0.5*((double)(j-1))*(t1+t2);
dz[j]= t1;
n[j]= n1;
g[j]= g1;
pe[j]= pe1;
pm[j]= pm1;
qe[j]= qe1;
qm[j]= qm1;
if(j==nn-1)z[nn]= z[nn-1]+t1;
}
for(j= 2;j<=nn-1;j= j+2){
z[j]= z[j-1]+t1;
dz[j]= t2;
n[j]= n2;
g[j]= g2;
pe[j]= pe2;
pm[j]= pm2;
qe[j]= qe2;
qm[j]= qm2;
if(j==nn-1)z[nn]= z[nn-1]+t2;
}
}else{
for(j= 1;j<=nn-1;j++){
if(j==1)z[j]= 0.0;
if(j==modnum){
z[j+1]= z[j]+modt1;
dz[j]= modt1;
n[j]= modn1;
g[j]= modg1;
pe[j]= modpe1;
pm[j]= modpm1;
qe[j]= modqe1;
qm[j]= modqm1;
}else{
tmp= ((double)j)/((double)2);
if(tmp-floor(tmp)> 0.25){
z[j+1]= z[j]+t1;
dz[j]= t1;
n[j]= n1;
g[j]= g1;
pe[j]= pe1;
pm[j]= pm1;
qe[j]= qe1;
qm[j]= qm1;
}else{
z[j+1]= z[j]+t2;
dz[j]= t2;
n[j]= n2;
g[j]= g2;
pe[j]= pe2;
pm[j]= pm2;
qe[j]= qe2;
qm[j]= qm2;
}
}
}
}
}
}else{
fprintf(stderr,"%s: Error.\n",progname);
fprintf(stderr,"%s: (No valid grating subtype found).\n",progname);
exit(FAILURE);
}
}

/*:65*/
#line 3987 "./magbragg.w"
;
}else if(!strcmp(gratingtype,"sinusoidal")){
/*66:*/
#line 4196 "./magbragg.w"

{
t1= ll/((double)(nn-1));
for(j= 1;j<=nn-1;j++){
z[j]= ((double)(j-1))*t1;
dz[j]= t1;
if(apodize){
if((0.0<=z[j])&&(z[j]<=apolength)){
tmp= (1.0-cos(pi*z[j]/apolength))/2.0;
}else if((apolength<=z[j])&&(z[j]<=ll-apolength)){
tmp= 1.0;
}else if((ll-apolength<=z[j])&&(z[j]<=ll)){
tmp= (1.0-cos(pi*(z[j]-ll)/apolength))/2.0;
}else{
tmp= 0.0;
fprintf(stderr,
"%s: Impossible apodization event occurred.",progname);
fprintf(stderr,
"%s: (Please check grating initialization.)",progname);
}
}
if(phasejump){
if(z[j]>=phasejumpposition){
phi= phasejumpangle;
}else{
phi= 0.0;
}
}
if(apodize){
n[j]= n1+n2*tmp*sin(twopi*z[j]/nper+phi);
g[j]= g1+g2*tmp*sin(twopi*z[j]/gper+phi);
}else{
n[j]= n1+n2*sin(twopi*z[j]/nper+phi);
g[j]= g1+g2*sin(twopi*z[j]/gper+phi);
}
pe[j]= pe1+pe2*sin(twopi*z[j]/peper);
pm[j]= pm1+pm2*sin(twopi*z[j]/pmper);
qe[j]= qe1+qe2*sin(twopi*z[j]/qeper);
qm[j]= qm1+qm2*sin(twopi*z[j]/qmper);
}
z[nn]= ll;
}

/*:66*/
#line 3989 "./magbragg.w"
;
}else if(!strcmp(gratingtype,"chirped")){
/*67:*/
#line 4249 "./magbragg.w"

{
t1= ll/((double)(nn-1));
for(j= 1;j<=nn-1;j++){
z[j]= ((double)(j-1))*t1;
dz[j]= t1;
if(apodize){
if((0.0<=z[j])&&(z[j]<=apolength)){
tmp= (1.0-cos(pi*z[j]/apolength))/2.0;
}else if((apolength<=z[j])&&(z[j]<=ll-apolength)){
tmp= 1.0;
}else if((ll-apolength<=z[j])&&(z[j]<=ll)){
tmp= (1.0-cos(pi*(z[j]-ll)/apolength))/2.0;
}else{
tmp= 0.0;
fprintf(stderr,
"%s: Impossible apodization event occurred.",progname);
fprintf(stderr,
"%s: (Please check grating initialization.)",progname);
}
}
if(phasejump){
if(z[j]>=phasejumpposition){
phi= phasejumpangle;
}else{
phi= 0.0;
}
}
if(apodize){
if(fabs(ncrp)> 0.0)
n[j]= n1+n2*tmp*sin((twopi/ncrp)*log(1.0+ncrp*z[j]/nper)+phi);
else
n[j]= n1+n2*tmp*sin(twopi*z[j]/nper+phi);
if(fabs(gcrp)> 0.0)
g[j]= g1+g2*tmp*sin((twopi/gcrp)*log(1.0+gcrp*z[j]/gper)+phi);
else
g[j]= g1+g2*tmp*sin(twopi*z[j]/gper+phi);
}else{
if(fabs(ncrp)> 0.0)
n[j]= n1+n2*sin((twopi/ncrp)*log(1.0+ncrp*z[j]/nper)+phi);
else
n[j]= n1+n2*sin(twopi*z[j]/nper+phi);
if(fabs(gcrp)> 0.0)
g[j]= g1+g2*sin((twopi/gcrp)*log(1.0+gcrp*z[j]/gper)+phi);
else
g[j]= g1+g2*sin(twopi*z[j]/gper+phi);
}
if(fabs(pecrp)> 0.0)
pe[j]= pe1+pe2*sin((twopi/pecrp)*log(1.0+pecrp*z[j]/peper));
else
pe[j]= pe1+pe2*sin(twopi*z[j]/peper);
if(pmcrp*pmcrp> 0.0)
pm[j]= pm1+pm2*sin((twopi/pmcrp)*log(1.0+pmcrp*z[j]/pmper));
else
pm[j]= pm1+pm2*sin(twopi*z[j]/pmper);
if(fabs(qecrp)> 0.0)
qe[j]= qe1+qe2*sin((twopi/qecrp)*log(1.0+qecrp*z[j]/qeper));
else
qe[j]= qe1+qe2*sin(twopi*z[j]/qeper);
if(fabs(qmcrp)> 0.0)
qm[j]= qm1+qm2*sin((twopi/qmcrp)*log(1.0+qmcrp*z[j]/qmper));
else
qm[j]= qm1+qm2*sin(twopi*z[j]/qmper);
}
z[nn]= ll;
}

/*:67*/
#line 3991 "./magbragg.w"
;
}else if(!strcmp(gratingtype,"fractal")){
/*68:*/
#line 4342 "./magbragg.w"

{
tmp= 1.0;
for(j= 1;j<=fractal_level-1;j++)tmp= 2.0*tmp;
for(j= 1;j<=maximum_fractal_level-fractal_level;j++)tmp= 3.0*tmp;

ll= tmp*t1-tmp*t2;
tmp= 1.0;
for(j= 1;j<=maximum_fractal_level-1;j++)tmp= 3.0*tmp;
ll= ll+tmp*t2;
if(verbose){
fprintf(stdout,"%s: Minimum layer thickness at maximum recursion depth:\n",
progname);
fprintf(stdout,"%s:     t2=%f [nm]\n",progname,t1*1.0e9);
fprintf(stdout,"%s:     t2=%f [nm]\n",progname,t2*1.0e9);
fprintf(stdout,"%s: Fractal grating length calculated as L=%e [m]\n",
progname,ll);
fprintf(stdout,
"%s: Based on fractal recursion of %d out of a maximum of %d levels.\n",
progname,fractal_level,maximum_fractal_level);
}
init_cantor_fractal_grating(z,1,nn,0.0,ll,n1,n2);
for(j= 1;j<=nn-1;j++){
dz[j]= z[j+1]-z[j];
if(j==1){
odd_layer= 1;
}else{
odd_layer= (odd_layer?0:1);
}
n[j]= (odd_layer?n1:n2);
g[j]= (odd_layer?g1:g2);
pe[j]= (odd_layer?pe1:pe2);
pm[j]= (odd_layer?pm1:pm2);
qe[j]= (odd_layer?qe1:qe2);
qm[j]= (odd_layer?qm1:qm2);
}
}

/*:68*/
#line 3993 "./magbragg.w"
;
}else{
fprintf(stderr,
"%s: Error: Specified grating type is invalid.\n",progname);
showsomehelp();
exit(FAILURE);
}
if(perturbed_gyration_constant){
/*69:*/
#line 4400 "./magbragg.w"
{
for(j= 1;j<=nn-1;j++){
tmp= 2.0*(z[j]-gyroperturb_position)/gyroperturb_width;
g[j]+= gyroperturb_amplitude/(1.0+tmp*tmp);
}
}

/*:69*/
#line 4001 "./magbragg.w"
;
}
}

/*:64*/
#line 3419 "./magbragg.w"

/*70:*/
#line 4426 "./magbragg.w"

{
n[0]= nsurr;
n[nn]= nsurr;
g[0]= 0.0;
g[nn]= 0.0;
}

/*:70*/
#line 3420 "./magbragg.w"

/*71:*/
#line 4483 "./magbragg.w"

{
for(j= 1;j<=nn;j++){
taup[j]= 2.0*(n[j-1]+g[j-1])/(n[j-1]+n[j]+g[j-1]+g[j]);
taum[j]= 2.0*(n[j-1]-g[j-1])/(n[j-1]+n[j]-g[j-1]-g[j]);
taupp[j]= 2.0*(n[j]-g[j])/(n[j-1]+n[j]-g[j-1]-g[j]);
taupm[j]= 2.0*(n[j]+g[j])/(n[j-1]+n[j]+g[j-1]+g[j]);
rhop[j]= (n[j-1]-n[j]+g[j-1]-g[j])/(n[j-1]+n[j]+g[j-1]+g[j]);
rhom[j]= (n[j-1]-n[j]-g[j-1]+g[j])/(n[j-1]+n[j]-g[j-1]-g[j]);
rhopp[j]= -rhom[j];
rhopm[j]= -rhop[j];
}
}

/*:71*/
#line 3421 "./magbragg.w"

/*72:*/
#line 4515 "./magbragg.w"

{
for(k= 1;k<=mm;k++){
if(mm> 1){
lambda= lambdastart+(((double)(k-1))/((double)(mm-1)))
*(lambdastop-lambdastart);
}else{
lambda= lambdastart;
}
omega= twopi*c/lambda;
if(trmtraject_specified){
mme= mmtraject;
mmi= 1;
}
/*74:*/
#line 4562 "./magbragg.w"

{
for(ke= 1;ke<=mme;ke++){
if(trmtraject_specified){
trmellipticity= w3traj[ke]/w0traj[ke];
}else{
if(mme> 1){
trmellipticity= trmellipstart+(((double)(ke-1))/((double)(mme-1)))
*(trmellipstop-trmellipstart);
}else{
trmellipticity= trmellipstart;
}
}
/*75:*/
#line 4589 "./magbragg.w"

{
for(ki= 1;ki<=mmi;ki++){
if(trmtraject_specified){
trmintensity= (epsilon0*c/2.0)*w0traj[ke];
}else{
if(mmi> 1){
trmintensity= trmintenstart
+(((double)(ki-1))/((double)(mmi-1)))
*(trmintenstop-trmintenstart);
}else{
trmintensity= trmintenstart;
}
}
/*76:*/
#line 4650 "./magbragg.w"

{
ebp[nn]= complex(0.0,0.0);
ebm[nn]= complex(0.0,0.0);
efp[nn]= complex(sqrt((1.0+trmellipticity)*trmintensity/(c*epsilon0)),0.0);
efm[nn]= complex(sqrt((1.0-trmellipticity)*trmintensity/(c*epsilon0)),0.0);
}

/*:76*/
#line 4603 "./magbragg.w"
;
/*77:*/
#line 4662 "./magbragg.w"

{
efp[nn-1]= cmul(crdiv(efp[nn],taup[nn]),crexpi(-omega*n[nn-1]*dz[nn-1]/c));
efm[nn-1]= cmul(crdiv(efm[nn],taum[nn]),crexpi(-omega*n[nn-1]*dz[nn-1]/c));
ebp[nn-1]= cmul(rcmul(rhom[nn],efm[nn-1]),
crexpi(2.0*omega*n[nn-1]*dz[nn-1]/c));
ebm[nn-1]= cmul(rcmul(rhop[nn],efp[nn-1]),
crexpi(2.0*omega*n[nn-1]*dz[nn-1]/c));
}

/*:77*/
#line 4604 "./magbragg.w"
;
/*78:*/
#line 4702 "./magbragg.w"

{
for(j= nn-1;j>=1;j--){
/*79:*/
#line 4762 "./magbragg.w"

{
aafp2= cdabs(efp[j]);
aafm2= cdabs(efm[j]);
aabp2= cdabs(ebp[j]);
aabm2= cdabs(ebm[j]);
aafp2*= aafp2;
aafm2*= aafm2;
aabp2*= aabp2;
aabm2*= aabm2;
tmp= 3.0/(8.0*n[j]);
etafp[j]= tmp*((pe[j]+pm[j])*(aafp2+2.0*aabm2)
+(qe[j]+qm[j])*(aafm2+aabp2));
etafm[j]= tmp*((pe[j]-pm[j])*(aafm2+2.0*aabp2)
+(qe[j]-qm[j])*(aafp2+aabm2));
etabp[j]= tmp*((pe[j]-pm[j])*(aabp2+2.0*aafm2)
+(qe[j]-qm[j])*(aabm2+aafp2));
etabm[j]= tmp*((pe[j]+pm[j])*(aabm2+2.0*aafp2)
+(qe[j]+qm[j])*(aabp2+aafm2));
}

/*:79*/
#line 4705 "./magbragg.w"
;
/*80:*/
#line 4795 "./magbragg.w"

{
tmpfp= cmul(efp[j],crexpi(-omega*(etafp[j]+g[j])*dz[j]/c));
tmpfm= cmul(efm[j],crexpi(-omega*(etafm[j]-g[j])*dz[j]/c));
tmpbp= cmul(ebp[j],crexpi(omega*(etabp[j]-g[j])*dz[j]/c));
tmpbm= cmul(ebm[j],crexpi(omega*(etabm[j]+g[j])*dz[j]/c));
}

/*:80*/
#line 4706 "./magbragg.w"
;
/*81:*/
#line 4808 "./magbragg.w"

{
if(j> 1){
efp[j-1]= crdiv(cmul(csub(tmpfp,rcmul(rhopm[j],tmpbm)),
crexpi(-omega*n[j-1]*dz[j-1]/c)),taup[j]);
efm[j-1]= crdiv(cmul(csub(tmpfm,rcmul(rhopp[j],tmpbp)),
crexpi(-omega*n[j-1]*dz[j-1]/c)),taum[j]);
ebp[j-1]= cadd(rcmul(taupp[j],
cmul(tmpbp,crexpi(omega*n[j-1]*dz[j-1]/c))),
rcmul(rhom[j],
cmul(efm[j-1],crexpi(2.0*omega*n[j-1]*dz[j-1]/c))));
ebm[j-1]= cadd(rcmul(taupm[j],
cmul(tmpbm,crexpi(omega*n[j-1]*dz[j-1]/c))),
rcmul(rhop[j],
cmul(efp[j-1],crexpi(2.0*omega*n[j-1]*dz[j-1]/c))));
}else{
efp[0]= crdiv(csub(tmpfp,rcmul(rhopm[1],tmpbm)),taup[1]);
efm[0]= crdiv(csub(tmpfm,rcmul(rhopp[1],tmpbp)),taum[1]);
ebp[0]= cadd(rcmul(taupp[1],tmpbp),rcmul(rhom[1],efm[0]));
ebm[0]= cadd(rcmul(taupm[1],tmpbm),rcmul(rhop[1],efp[0]));
}
}

/*:81*/
#line 4707 "./magbragg.w"
;
if(verbose){
/*82:*/
#line 4852 "./magbragg.w"

{
modf(100.0*((float)((nn-j-1)+(ki-1)*(nn-1)
+(ke-1)*mmi*(nn-1)+(k-1)*mme*mmi*(nn-1)))
/((float)(mm*mme*mmi*(nn-1))),&stn);
if(stn> status+10){
status= status+10;
now= time(NULL);
eta= initime
+((int)((100.0/((double)status))*difftime(now,initime)));
fprintf(stdout," ...%2d percent finished...   ",status);
fprintf(stdout," ETA: %s",ctime(&eta));
}
}

/*:82*/
#line 4709 "./magbragg.w"
;
}
}
}

/*:78*/
#line 4605 "./magbragg.w"
;
/*83:*/
#line 5021 "./magbragg.w"

{


if(intensityinfo){
for(j= 0;j<=nn;j++){
tmp= (epsilon0*n[j]*c/2)*(cabs2(efp[j])+cabs2(efm[j]));
if(maxintens<tmp){
maxintens= tmp;
maxintens_layer= j;
maxintens_inintens= (epsilon0*nsurr*c/2)
*(cabs2(efp[0])+cabs2(efm[0]));
maxintens_inellip= (cabs2(efp[0])-cabs2(efm[0]))
/(cabs2(efp[0])+cabs2(efm[0]));
maxintens_trintens= (epsilon0*nsurr*c/2)
*(cabs2(efp[nn])+cabs2(efm[nn]));
maxintens_trellip= (cabs2(efp[nn])-cabs2(efm[nn]))
/(cabs2(efp[nn])+cabs2(efm[nn]));
}
}
}
if((mme> 1)&&(mmi> 1)){

s0= cabs2(efp[0])+cabs2(efm[0]);
s1= 2.0*cmul(conjg(efp[0]),efm[0]).r;
s2= 2.0*cmul(conjg(efp[0]),efm[0]).i;
s3= cabs2(efp[0])-cabs2(efm[0]);
if(scale_stokesparams){
s0= s0*stoke_scalefactor;
s1= s1*stoke_scalefactor;
s2= s2*stoke_scalefactor;
s3= s3*stoke_scalefactor;
}
if(normalize_ellipticity)s3= s3/s0;
fprintf(fp_s0,"%16.12e  ",s0);
fprintf(fp_s1,"%16.12e  ",s1);
fprintf(fp_s2,"%16.12e  ",s2);
fprintf(fp_s3,"%16.12e  ",s3);
fflush(fp_s0);
fflush(fp_s1);
fflush(fp_s2);
fflush(fp_s3);


v0= cabs2(ebp[0])+cabs2(ebm[0]);
v1= 2.0*cmul(conjg(ebp[0]),ebm[0]).r;
v2= 2.0*cmul(conjg(ebp[0]),ebm[0]).i;
v3= cabs2(ebp[0])-cabs2(ebm[0]);
if(scale_stokesparams){
v0= v0*stoke_scalefactor;
v1= v1*stoke_scalefactor;
v2= v2*stoke_scalefactor;
v3= v3*stoke_scalefactor;
}
if(normalize_ellipticity)v3= v3/v0;
fprintf(fp_v0,"%16.12e  ",v0);
fprintf(fp_v1,"%16.12e  ",v1);
fprintf(fp_v2,"%16.12e  ",v2);
fprintf(fp_v3,"%16.12e  ",v3);
fflush(fp_v0);
fflush(fp_v1);
fflush(fp_v2);
fflush(fp_v3);


w0= cabs2(efp[nn])+cabs2(efm[nn]);
w1= 2.0*cmul(conjg(efp[nn]),efm[nn]).r;
w2= 2.0*cmul(conjg(efp[nn]),efm[nn]).i;
w3= cabs2(efp[nn])-cabs2(efm[nn]);
if(scale_stokesparams){
w0= w0*stoke_scalefactor;
w1= w1*stoke_scalefactor;
w2= w2*stoke_scalefactor;
w3= w3*stoke_scalefactor;
}
if(normalize_ellipticity)w3= w3/w0;
fprintf(fp_w0,"%16.12e  ",w0);
fprintf(fp_w1,"%16.12e  ",w1);
fprintf(fp_w2,"%16.12e  ",w2);
fprintf(fp_w3,"%16.12e  ",w3);
fflush(fp_w0);
fflush(fp_w1);
fflush(fp_w2);
fflush(fp_w3);
}else{
if(stokes_parameter_spectrum){

s0= cabs2(efp[0])+cabs2(efm[0]);
s1= 2.0*cmul(conjg(efp[0]),efm[0]).r;
s2= 2.0*cmul(conjg(efp[0]),efm[0]).i;
s3= cabs2(efp[0])-cabs2(efm[0]);
fprintf(fp_spec,"%16.12e %16.12e %16.12e %16.12e %16.12e\n",
lambda,omega,s1/s0,s2/s0,s3/s0);
fflush(fp_spec);


v0= cabs2(ebp[0])+cabs2(ebm[0]);
v1= 2.0*cmul(conjg(ebp[0]),ebm[0]).r;
v2= 2.0*cmul(conjg(ebp[0]),ebm[0]).i;
v3= cabs2(ebp[0])-cabs2(ebm[0]);



w0= cabs2(efp[nn])+cabs2(efm[nn]);
w1= 2.0*cmul(conjg(efp[nn]),efm[nn]).r;
w2= 2.0*cmul(conjg(efp[nn]),efm[nn]).i;
w3= cabs2(efp[nn])-cabs2(efm[nn]);

}
if(save_dbspectra){
tmp= (cabs2(ebp[0])+cabs2(ebm[0]))/(cabs2(efp[0])+cabs2(efm[0]));
fprintf(fp_irspec,"%-10.8f %-10.8f\n",lambda*1.0e9,10.0*log10(tmp));
fprintf(fp_itspec,"%-10.8f %-10.8f\n",lambda*1.0e9,10.0*log10(1.0-tmp));
fflush(fp_irspec);
fflush(fp_itspec);
}else{
fprintf(fp_irspec,"%-10.8f %-10.8f\n",lambda*1.0e9,
(cabs2(ebp[0])+cabs2(ebm[0]))/(cabs2(efp[0])+cabs2(efm[0])));
fprintf(fp_itspec,"%-10.8f %-10.8f\n",lambda*1.0e9,
(cabs2(efp[nn])+cabs2(efm[nn]))/(cabs2(efp[0])+cabs2(efm[0])));
fflush(fp_irspec);
fflush(fp_itspec);
}
fprintf(fp_icspec,"%-10.8f %-10.8f\n",lambda*1.0e9,
(cabs2(ebp[0])+cabs2(ebm[0]))/(cabs2(efp[0])+cabs2(efm[0]))
+(cabs2(efp[nn])+cabs2(efm[nn]))/(cabs2(efp[0])+cabs2(efm[0])));
fflush(fp_icspec);
}

if(ki>=mmi){
if((mme> 1)&&(mmi> 1)){
fprintf(fp_s0,"\n");
fprintf(fp_s1,"\n");
fprintf(fp_s2,"\n");
fprintf(fp_s3,"\n");
fprintf(fp_v0,"\n");
fprintf(fp_v1,"\n");
fprintf(fp_v2,"\n");
fprintf(fp_v3,"\n");
fprintf(fp_w0,"\n");
fprintf(fp_w1,"\n");
fprintf(fp_w2,"\n");
fprintf(fp_w3,"\n");
}
}
}

/*:83*/
#line 4606 "./magbragg.w"
;
/*84:*/
#line 5179 "./magbragg.w"

{
if(fieldevoflag){
if(fieldevoflag_efield){
if(verbose)
fprintf(stdout,
"Writing spatial field evolution to file.\n");
if(fp_evo!=NULL){
for(j= 1;j<=nn-1;j++){
for(jje= 1;jje<=nne;jje++){
if(nne> 1){
zt= z[j]+((double)(jje-1))*dz[j]/((double)(nne));
}else{
zt= z[j];
}
tmpfp= cmul(efp[j],
crexpi(omega*(etafp[j]+g[j])*(zt-z[j])/c));
tmpfm= cmul(efm[j],
crexpi(omega*(etafm[j]-g[j])*(zt-z[j])/c));
tmpbp= cmul(ebp[j],
crexpi(-omega*(etabp[j]-g[j])*(zt-z[j])/c));
tmpbm= cmul(ebm[j],
crexpi(-omega*(etabm[j]+g[j])*(zt-z[j])/c));
if(normalize_length_to_micrometer){
fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e %16.12e"
" %16.12e\n",zt*1.0e6,
tmpfp.r,tmpfp.i,tmpfm.r,tmpfm.i,
tmpbp.r,tmpbp.i,tmpbm.r,tmpbm.i);
}else{
fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e %16.12e"
" %16.12e\n",zt,
tmpfp.r,tmpfp.i,tmpfm.r,tmpfm.i,
tmpbp.r,tmpbp.i,tmpbm.r,tmpbm.i);
}
}
}
if(normalize_length_to_micrometer){
fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e %16.12e\n",
z[nn]*1.0e6,efp[nn].r,efp[nn].i,efm[nn].r,efm[nn].i,
ebp[nn].r,ebp[nn].i,ebm[nn].r,ebm[nn].i);
}else{
fprintf(fp_evo,"%16.12e %16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e %16.12e\n",
z[nn],efp[nn].r,efp[nn].i,efm[nn].r,efm[nn].i,
ebp[nn].r,ebp[nn].i,ebm[nn].r,ebm[nn].i);
}
}else{
fprintf(stderr,"%s: Could not write to file %s!\n",
progname,fieldevofilename);
exit(FAILURE);
}
}else if(fieldevoflag_stoke){
if(verbose)
fprintf(stdout,
"Writing spatial evolution of Stokes parameters to file.\n");
if((fp_evo_s0!=NULL)&&(fp_evo_s1!=NULL)
&&(fp_evo_s2!=NULL)&&(fp_evo_s3!=NULL)){
for(j= 1;j<=nn-1;j++){
for(jje= 1;jje<=nne;jje++){
if(nne> 1){
zt= z[j]+((double)(jje-1))*dz[j]/((double)(nne));
}else{
zt= z[j];
}
tmpfp= cmul(efp[j],
crexpi(omega*(etafp[j]+g[j])*(zt-z[j])/c));
tmpfm= cmul(efm[j],
crexpi(omega*(etafm[j]-g[j])*(zt-z[j])/c));
tmpbp= cmul(ebp[j],
crexpi(-omega*(etabp[j]-g[j])*(zt-z[j])/c));
tmpbm= cmul(ebm[j],
crexpi(-omega*(etabm[j]+g[j])*(zt-z[j])/c));
s0= cabs2(tmpfp)+cabs2(tmpfm);
s1= 2.0*cmul(conjg(tmpfp),tmpfm).r;
s2= 2.0*cmul(conjg(tmpfp),tmpfm).i;
s3= cabs2(tmpfp)-cabs2(tmpfm);
if(normalize_intensity)s0= s0/(cabs2(efp[1])+cabs2(efm[1]));
if(normalize_length_to_micrometer){
fprintf(fp_evo_s0,"%16.12e %16.12e\n",zt*1.0e6,s0);
fprintf(fp_evo_s1,"%16.12e %16.12e\n",zt*1.0e6,s1);
fprintf(fp_evo_s2,"%16.12e %16.12e\n",zt*1.0e6,s2);
fprintf(fp_evo_s3,"%16.12e %16.12e\n",zt*1.0e6,s3);
}else{
fprintf(fp_evo_s0,"%16.12e %16.12e\n",zt,s0);
fprintf(fp_evo_s1,"%16.12e %16.12e\n",zt,s1);
fprintf(fp_evo_s2,"%16.12e %16.12e\n",zt,s2);
fprintf(fp_evo_s3,"%16.12e %16.12e\n",zt,s3);
}
}
}
s0= cabs2(efp[nn])+cabs2(efm[nn]);
s1= 2.0*cmul(conjg(efp[nn]),efm[nn]).r;
s2= 2.0*cmul(conjg(efp[nn]),efm[nn]).i;
s3= cabs2(efp[nn])-cabs2(efm[nn]);
if(normalize_intensity)s0= s0/(cabs2(efp[1])+cabs2(efm[1]));
if(normalize_length_to_micrometer){
fprintf(fp_evo_s0,"%16.12e %16.12e\n",z[nn]*1.0e6,s0);
fprintf(fp_evo_s1,"%16.12e %16.12e\n",z[nn]*1.0e6,s1);
fprintf(fp_evo_s2,"%16.12e %16.12e\n",z[nn]*1.0e6,s2);
fprintf(fp_evo_s3,"%16.12e %16.12e\n",z[nn]*1.0e6,s3);
}else{
fprintf(fp_evo_s0,"%16.12e %16.12e\n",z[nn],s0);
fprintf(fp_evo_s1,"%16.12e %16.12e\n",z[nn],s1);
fprintf(fp_evo_s2,"%16.12e %16.12e\n",z[nn],s2);
fprintf(fp_evo_s3,"%16.12e %16.12e\n",z[nn],s3);
}
}else{
fprintf(stderr,"%s: Could not write to file %s, %s, %s, or %s!\n",
progname,fieldevofilename_s0,fieldevofilename_s1,
fieldevofilename_s2,fieldevofilename_s3);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Unknown field evolution flag.\n"
"%s: (This cannot happen)\n",progname,progname);
exit(FAILURE);
}
}
}

/*:84*/
#line 4607 "./magbragg.w"
;
/*85:*/
#line 5306 "./magbragg.w"

{
if(intensityevoflag){
if(fabs(ievolambda-lambda)
<fabs(lambdastop-lambdastart)/((double)(mm))){
if(verbose)
fprintf(stdout,"%s: Saving intensity evolution at lambda=%8.4e\n",
progname,lambda);
if(strcmp(intensityevofilename,"")){
if(fp_ievo!=NULL){
for(j= 1;j<=nn-1;j++){
if(normalize_length_to_micrometer){
fprintf(fp_ievo,"%16.12e %16.12e\n",z[j]*1.0e6,
(cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
/(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
fprintf(fp_ievo,"%16.12e %16.12e\n",z[j+1]*1.0e6,
(cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
/(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
}else{
fprintf(fp_ievo,"%16.12e %16.12e\n",z[j],
(cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
/(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
fprintf(fp_ievo,"%16.12e %16.12e\n",z[j+1],
(cdabs(efp[j])*cdabs(efp[j])+cdabs(efm[j])*cdabs(efm[j]))
/(cdabs(efp[1])*cdabs(efp[1])+cdabs(efm[1])*cdabs(efm[1])));
}
}
}else{
fprintf(stderr,"%s: Could not write to file %s!\n",
progname,intensityevofilename);
exit(FAILURE);
}
}
}
}
}

/*:85*/
#line 4608 "./magbragg.w"
;
/*86:*/
#line 5348 "./magbragg.w"

{
if(writegratingtofile){
if((fp_gr= fopen(gratingfilename,"w"))==NULL){
fprintf(stderr,"%s: Could not open file %s for output!\n",
progname,gratingfilename);
exit(FAILURE);
}
if(normalize_length_to_micrometer){
if(display_surrounding_media){
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
(z[1]-0.1*(z[nn]-z[1]))*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[1]*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
}
for(j= 1;j<=nn-1;j++){
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[j]*1.0e6,n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[j+1]*1.0e6,n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
}
if(display_surrounding_media){
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[nn]*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
(z[nn]+0.1*(z[nn]-z[1]))*1.0e6,nsurr,0.0,0.0,0.0,0.0,0.0);
}
}else{
if(display_surrounding_media){
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[1]-0.1*(z[nn]-z[1]),nsurr,0.0,0.0,0.0,0.0,0.0);
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[1],nsurr,0.0,0.0,0.0,0.0,0.0);
}
for(j= 1;j<=nn-1;j++){
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[j],n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[j+1],n[j],g[j],pe[j],pm[j],qe[j],qm[j]);
}
if(display_surrounding_media){
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[nn],nsurr,0.0,0.0,0.0,0.0,0.0);
fprintf(fp_gr,"%16.12e %16.12e %16.12e %16.12e"
" %16.12e %16.12e %16.12e\n",
z[nn]+0.1*(z[nn]-z[1]),nsurr,0.0,0.0,0.0,0.0,0.0);
}
}
fclose(fp_gr);
}
}

/*:86*/
#line 4609 "./magbragg.w"
;
}
}

/*:75*/
#line 4575 "./magbragg.w"
;
}
}

/*:74*/
#line 4529 "./magbragg.w"
;
}
if(verbose)/*73:*/
#line 4538 "./magbragg.w"

{
fprintf(stdout," ...done.                     ");
now= time(NULL);
fprintf(stdout,"Elapsed execution time: %d s\n",
((int)difftime(now,initime)));
for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
fprintf(stdout,"Program execution closed %s",ctime(&now));
}

/*:73*/
#line 4531 "./magbragg.w"
;
}

/*:72*/
#line 3422 "./magbragg.w"

/*87:*/
#line 5413 "./magbragg.w"

{
if(intensityinfo){
for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
fprintf(stdout,"Summary of intra-grating intensities:\n");
fprintf(stdout,
"The maximum intensity %1.4e [W/m^2] (%1.4f GW/cm^2) was detected\n",
maxintens,maxintens*1.0e-13);
fprintf(stdout,
"in layer %ld. The maximum intensity was detected at a transmitted\n",
maxintens_layer);
fprintf(stdout,
"intensity of %1.4e [W/m^2] (%1.4f GW/cm^2), and at a transmitted\n",
maxintens_trintens,maxintens_trintens*1.0e-13);
fprintf(stdout,
"normalized ellipticity of S3/S0=%1.4f. ",maxintens_trellip);
fprintf(stdout,"(where S3/S0=1 for LCP, and -1\n""for RCP).");
fprintf(stdout,
"At this state, the calculated optical intensity incident to the\n");
fprintf(stdout,
"crystal was %1.4e [W/m^2] (%1.4f GW/cm^2), or %1.1f percent\n",
maxintens_inintens,maxintens_inintens*1.0e-13,
100.0*maxintens_inintens/maxintens);
fprintf(stdout,"of the maximum intra-grating optical intensity.\n");
fprintf(stdout,
"The calculated normalized incident ellipticity was %1.4f.\n",
maxintens_inellip);
fprintf(stdout,
"The intensity transmission was for this state %1.1f percent.\n",
100.0*maxintens_trintens/maxintens_inintens);
for(k= 1;k<=64;k++)fprintf(stdout,(k<64?"-":"\n"));
if(saveintensityinfologfile){
if((intensinfologfile= fopen(intensinfologfilename,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for intensity log!\n",
progname,intensinfologfilename);
exit(FAILURE);
}
for(k= 1;k<=32;k++)fprintf(intensinfologfile,(k<32?"-":"\n"));
fprintf(intensinfologfile,"Summary of intra-grating intensities:\n");
fprintf(intensinfologfile,
"Maximum intensity %1.4e [W/sq.m] (%1.4f GW/sq.cm) was detected\n",
maxintens,maxintens*1.0e-13);
fprintf(intensinfologfile,
"in layer %ld. Maximum intensity was detected at a transmitted\n",
maxintens_layer);
fprintf(intensinfologfile,
"intensity of %1.4e [W/sq.m] (%1.4f GW/sq.cm), and at transmitted\n",
maxintens_trintens,maxintens_trintens*1.0e-13);
fprintf(intensinfologfile,
"normalized ellipticity of S3/S0=%1.4f. ",maxintens_trellip);
fprintf(intensinfologfile,
"(where S3/S0=1 for LCP, and -1\n""for RCP). ");
fprintf(intensinfologfile,
"At this state, the calculated optical intensity incident to the\n");
fprintf(intensinfologfile,
"crystal was %1.4e [W/sq.m] (%1.4f GW/sq.cm), or %1.1f percent\n",
maxintens_inintens,maxintens_inintens*1.0e-13,
100.0*maxintens_inintens/maxintens);
fprintf(intensinfologfile,
"of the maximum intra-grating optical intensity.\n");
fprintf(intensinfologfile,
"The calculated normalized incident ellipticity was %1.4f.\n",
maxintens_inellip);
fprintf(intensinfologfile,
"The intensity transmission was for this state %1.1f percent.\n",
100.0*maxintens_trintens/maxintens_inintens);
fclose(intensinfologfile);
for(k= 1;k<=32;k++)fprintf(intensinfologfile,(k<32?"-":"\n"));
}
}
}

/*:87*/
#line 3423 "./magbragg.w"

/*63:*/
#line 3944 "./magbragg.w"

{
free_dcvector(efp,0,nn);
free_dcvector(efm,0,nn);
free_dcvector(ebp,0,nn);
free_dcvector(ebm,0,nn);
free_dvector(taup,1,nn);
free_dvector(taum,1,nn);
free_dvector(taupp,1,nn);
free_dvector(taupm,1,nn);
free_dvector(rhop,1,nn);
free_dvector(rhom,1,nn);
free_dvector(rhopp,1,nn);
free_dvector(rhopm,1,nn);
free_dvector(n,0,nn);
free_dvector(g,0,nn);
free_dvector(pe,1,nn-1);
free_dvector(pm,1,nn-1);
free_dvector(qe,1,nn-1);
free_dvector(qm,1,nn-1);
free_dvector(etafp,1,nn-1);
free_dvector(etabp,1,nn-1);
free_dvector(etafm,1,nn-1);
free_dvector(etabm,1,nn-1);
if(trmtraject_specified){
free_dvector(w0traj,1,mmtraject);
free_dvector(w3traj,1,mmtraject);
}
}

/*:63*/
#line 3424 "./magbragg.w"

/*136:*/
#line 7725 "./magbragg.w"

{
if((mme> 1)&&(mmi> 1)){
fclose(fp_s0);
fclose(fp_s1);
fclose(fp_s2);
fclose(fp_s3);
fclose(fp_v0);
fclose(fp_v1);
fclose(fp_v2);
fclose(fp_v3);
fclose(fp_w0);
fclose(fp_w1);
fclose(fp_w2);
fclose(fp_w3);
}
if(fieldevoflag)if(strcmp(fieldevofilename,""))fclose(fp_evo);
if(intensityevoflag)if(strcmp(intensityevofilename,""))fclose(fp_ievo);
if(!((mme> 1)&&(mmi> 1))){
fclose(fp_spec);
fclose(fp_irspec);
fclose(fp_itspec);
fclose(fp_icspec);
}
}

/*:136*/
#line 3425 "./magbragg.w"

return(SUCCESS);
}

/*:45*/