Search:

Return to previous page

Contents of file 'magbragg/magbragg.c':




Deprecated: Function split() is deprecated in /storage/content/45/2011745/jonsson.eu/public_html/php/htmlicise.php on line 46
    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   

Return to previous page

Generated by ::viewsrc::

Last modified Saturday 24 Dec 2011