/* Balloon Laue-lens optimization. Niels Lund March 2006. */
  /* Three-layer lens version December 2018 - March 2019 */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/*                                                                                */        
/*              Three layer Bragg telescope, Niels Lund, December 2018            */        
/*                                                                                */
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/*                                                                                                                 */        
/* Overall description of program assumptions and logic                                                            */
/*         1) The lens consists of three layers of crystal facets, each layer arranged in a hexagonal pattern      */
/*                                                                                                                 */        
/*         2) The centers of the crystals in each layer are offset from the other two layers so as to obtain       */
/*            the uniform coverage of the three layers combined.                                                   */
/*                                                                                                                 */        
/*         3) The lens may be divided into a number of radial zones. Within each zone the diffraction option in    */
/*            all three layers are kept constant. (The diffraction option specifies the crystal type (Z), the      */
/*            diffracton order and type (Miller index) and the mosaic width). However, the crystal thickness may   */
/*            vary smoothly across the zone to optimize the diffraction yield over the energy bands diffracted     */
/*            from this zone.                                                                                      */
/*            The number of zones and their characteristics are defined in the parameter file ZONES.dat            */
/*                                                                                                                 */
/*         Program logic: (Program parts are separated by the string 'QQQQQ'                                       */
/*                                                                                                                 */
/*         A) Define the (X,Y) coordinates of all the facet positions in the three layers. Assign zones to facets  */
/*            Line QQQ01 (422) - QQQ02 (470).                                                                      */        
/*                                                                                                                 */
/*         B) Establish neighbor relations between facets in layer 1 with 2, 2 with 3 and 3 with 1.                */
/*            Line QQQ03 (570) - QQQ04 (470).                                                                      */        
/*                                                                                                                 */        
/*         C) Assign diffraction options to facets according to the specifications in 'ZONES.dat'.                 */
/*            Note that some facet positions may be left empty.                                                    */
/*            Line QQQ05 (570) - QQQ06 (470).                                                                      */        
/*                                                                                                                 */        
/*         D) Choose the thickness of each facet according to its diffraction efficiency, absorption properties    */
/*            and the possible overlap with neighboring crystals in the other two layers.                          */
/*            Line QQQ07 (422) - QQQ08 (470).                                                                      */        
/*                                                                                                                 */        
/*         E) Calculate the flux collection area as function of energy for the lens                                */
/*            Line QQQ09 (422) - QQQ10 (470).                                         				   */
/*                                                                                                                 */        
/*         F) Terminate by preparing various numerical and graphic output                                          */ 
/*            Line QQQ11 (422) - QQQ12 (470).                                                                      */        
/*                                                                                                                 */        
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "mucal.h"
#include "panel.h"
#include <time.h>

int main(int argc, char* argv[])
{
   FILE *in, *out, *out2, *listfil, *preget, *gnufil, *pltfil, *chkfil, *gnufil2, *gnuFIL[2][20], *zone_file, *thkfil;

   static char Miller[6], iline[300], textlines[20][128], Titl[32], gnutitl[64], Lauetitl[64], Outtitl[64];
   static char ELEMENT[20][4], IDtext[16], gnuFILname[32], gnutitl1[64], thktitl[64];
   char Millerlist[20][5], Chktitl[64], zone_name[64], systemtxt[64], PStitl[64];
   char outline2[64], outline3[64], outline4[64], triptxt[16];

   static int L_order[numLayers][N_CryPos];	/* Sequence numbering of crystals in the three crystal layers [0/1/2]
		                                        - ordered according to distance from the axis */

   int i, ii, i_elem, j, k, m, n_order, tst, newcry, r_imax, idist, loi, E_indx0, nl;
   int FF=1, iR, n_z, num_zones=0, iE, iE_low, iE_high, nOpt, iOrder, mE, gnunum, RING_NUM=0;
   int Z, Cmiller, control, Len, r_index, E1k, E2k, E3k, E4k, E5k, Ek, d_index, EE_index;
   int status = 0, Au_111, Au_200, Cu_111, Cu_200, Ag_111, Ag_200, Ge111, Ge220, r_imin;
   int Cu220, Ag220, Au220, chan[100], L_start[600], gnufnum=0, new_gnufnum, dot_gnufnum;
   int iC, iiC, ring_i, elements, g_index, IE, elem_save, imax, IntE, Efound[20][100];
   int basic_elem, basic_order, j_old, r_num[600], updated[20], update_gate, k_dist[8][20];
   int elem, E_indx, E_imax, num_crys[20], cry_ptr_lo[20], cry_ptr_hi[20], selected[20];    
   int Nrand0=0, Nrand1=0, locnum[20], kk[20], listptr, rest_elem, combnum, n_layer=0;
   int DetArea, y2min, y2max, Cu111, Cu200, Ag111, Ag200, Au111, Au200, Zsave, opt2;
   int maxRing, i_radius, imin, preset_i[20][100], chan_max, chan_index, lim_chan1, lim_chan2;;
   int Lens_layers=2, Z2 = 6, Cmiller2 = 99, n_order2 = 1, nOpt0[2] = { 0, 1};
   int IZ[6] = {0, 0, 0, 0, 0, 0}, znum, kmax, old_numOp, Ngnu[20], OptionZ[20];
   int facet_no[3], facet_init[3], last_init[3], uselen[3], NNopt=0;  
   int npr=0, layer_no, ring_no, good_facets, pr_ok=0;
   int last_ring, first_ring=0, last_good=0, allnabs;
   int azimnum=0, modulLim=0, modulLim0=0, Zside[2], CMside[2], numOptions=0, atoms_cell[5];
   int azim_index, N_z, NN, prcnt=0, GG=0, GGOK=0;
   int iaz, ik, ik_min, mm, KK, ii_ref, iii, jjj, ring_num, Xrest[6], option2_num, z_ring;
   int tlnum, iseq, basic_Z=0, MAXcm=0, INITfacets=0, jold=0, N_sec[2][3]={0, 0, 0, 0, 0, 0};

   unsigned LEN[4];

   double E, E2, E3, fwhmt, T, tant2B, B_radius, radius_old, thetaB, lambda, dE;
   double gain, fwhm_cm[20], t_save;
   double ny_use, t_use, E0, minthick, mufact0;
   double x, y, alfa, fi_lo, fi_hi, fi, fwhmt2 = 1.0, T_FACTOR = 0.35;
   
   static int distr_layer[3][400];
   static float r_area[600], nysum[N_CryPos*3], OAz[500];
   static float  EAALL[3+2][N_CryPos], AZIM_REF[6]={0.0, 180.0, 60.0, 240.0, 120.0, 300.0};
   static float T_Area[500][5], Eval_pos[20][600];
   static float EffArea[3*2][20][4][N_CryPos];   /* EffArea [Layer contribution] [nOption] [nOrder] [energy in 2 keV steps] */
   static float EffAOpt[3*2][20][4][N_CryPos];
   static float Agnu[20][3000];
   static int   Rgnu[3][20][100];

   float RR, ms, sum0, sum1, sum2, sum3, sum4, sum5, sum6, disrad, val, valsum, fwhmth[20], t_Factor[20];
   float m_tot[20], fwhm_fac, sigma_fac, sigma_cm, offpeak, arcmin2rad, Rf, facet_R, dmax;
   float ring_width, radius, panel3[25], LensTopRadius, lastR, thick, thickx;
   float x_focus0, y_focus0, facet_diam, lastBragg, MirrorWidth, MirrorDispl, az_min, az_dif, azim_ref;
   float x_focus2, y_focus2, T_fctop, T_fcbot, ModE;
   float x_focus1, y_focus1, facet_radius, facet_hole_radius, detector_diam, Rdf, dis, overlap_frac;
   float lens_diam, lens_hole_diam, cameraH, cameraD, det_area, EAall, EAall2, E_ref0=511.0, EAL = 1.0;
   float sum_area[600], sumarea_E[20], A_rstloc, E_ref[20], iElo[6];
   float rhoval[20], A_chk[3][5], tri_wgt[100], xmean[100], ymean[100], xsum, ysum;
   float FRAC, Agnusum[20], A_wanted, separation, M_tot=0.0;
   float sigmarad, mu, R, nyuse, wght1, sqrt2pi, hexradius, dist, ov;
   float delta_r, old_radius, beta1, beta2, beta3, beta4, divisor;
   float gamma1, gamma2, gamma3, gamma4, gamma5, alfa1, wg, sumfi;
   float delta1, delta2, delta3, delta4, delta5, delta6;
   float Foc_len, H, H_focus, H_break1, H_break2, dH1, dH2, focus_sep, R_Tel;
   float gamuse1, gamuse2, gamuse3, x_focus, y_focus;
   float tB, tB0, Weight[20], Area[20][500], radius_lo, radius_hi;
   float radi[20], radi_lo[20], radi_hi[20], cry_frac[20][500], sum_Area, sumA;
   float distance, Brg_C, Brg_E, delta_Br, Brg_plus, Brg_minus;
   float E_C, refl_E, ref_plus, ref_minus, ref_C, cry_refl_E, WEIGHT=0.0, WEIGHT2=0.0;
   float Hsat, E_factor[20], radiusE0, accum_A, Evalmax;
   float unit_mass, weight[20], Rand_max, used_fi[20][600], fracfi[20][600];
   float theta1, theta2, panel_width, dw, T_fclayer;
   float D_switch, sep_min, thetaref, A_factor, facet_sep;
   float deltachan, chanstep, chan_val1, chan_val2, preset[21][600];
   float Cgeo[5][5], atoms_eff[5][5], A_basic, A_hole, crystal_mass[20];
   float E_low, E_high, A_contr, B_contr, t_global=1.0, lens_inner_diam;
   float Bragg_min, Bragg_max, NABSTAT[3][30], rLim=500.0;
   float thkarray[2][numFACETS], thk_sum=0.0, use_frac, lo_azim_lim, hi_azim_lim;

   float PassBandLow=300.0, PassBandHigh= 800.0;   /* Energy pass band for the lens */
   
   int CameHere=0;

/* The only legal actions on a structure are:
 * 	copying it
 * 	assigning to it as a unit
 * 	taking its address with &
 * 	accessing its members
 * 	(from Kernigan and Ritchey)
*/
   reflect REFL[20], *PREFL;
   reflect REFL2[20], *PREFL2;
   zones ZONES[20], *PZONES;
   material MATER[100], *PMATER;
   static facet FACET[numFACETS];	/* numFACETS is defined in panel.h */
   facet *PFACET;

   char navn[150], tmptxt[128], tmptxt2[150], dtext[130], dtext2[30];
   char dtext3[70], dtext4[70], Xlabel[256], xlabel[64], date_txt[64];
   
   time_t mytime;

/* QQQQQ */

   PREFL  = &REFL[0];
   PREFL2 = &REFL2[0];
   PZONES = &ZONES[0];
   PMATER = &MATER[0];
   PFACET = &FACET[0];

   PI    = 3.14159;
   TWOPI = PI * 2.0;
   GTORAD = PI / 180.0;
   sqrt2pi = sqrt(TWOPI);
   arcmin2rad = GTORAD / 60.0;
   Rand_max = (float)(RAND_MAX);
   
   for (i=0; i<N_CryPos*numLayers; i++) {
      EAALL[0][i] = EAALL[1][i] = EAALL[2][i] = EAALL[3][i] = 0.0;
      EAALL[4][i] = EAALL[5][i] = 0.0;
      nysum[i] = 0.0;
   }
   for (i=0; i<N_CryPos*numLayers; i++) {
      FACET[i].nOption = -1;
      FACET[i].nabo[0][0] = FACET[i].nabo[0][1] = FACET[i].nabo[0][2] = -1;
      FACET[i].nabo[1][0] = FACET[i].nabo[1][1] = FACET[i].nabo[1][2] = -1;
   }

   iElo[0] = 150;
   iElo[1] = 250;
   iElo[2] = 350;
   iElo[3] = 700;
   iElo[4] = 1000;
   iElo[5] = 2200;
   
   time(&mytime);
   sprintf(dtext2, "%s", ctime(&mytime));
   i = strlen(dtext2);
   dtext2[i-1] = 0;
   sprintf(date_txt, &dtext2[4]);

   for (i=0; i<Lens_layers; i++) for (j=0; j<30; j++) NABSTAT[i][j] = 0.0;

   for (j=0; j<20; j++) {
      for (i=0; i<N_CryPos; i++) {
         for (k=0; k<numLayers; k++) {
            EffArea[k*2+0][j][0][i] = EffArea[k*2+0][j][1][i] = EffArea[k*2+0][j][2][i] = EffArea[k*2+0][j][3][i] = 0.0;
            EffArea[k*2+1][j][0][i] = EffArea[k*2+1][j][1][i] = EffArea[k*2+1][j][2][i] = EffArea[k*2+1][j][3][i] = 0.0;
            L_order[k][i] = -1;
         }
      }
   }
   for (i=0; i<3000; i++) {
      for (k=0; k<20; k++) Agnu[k][i] = 0.0;
   }
   for (i=0; i<100; i++) {
     for (j=0; j<20; j++) {
        Efound[j][i] = 0;
        for (k=0; k<numLayers; k++) Rgnu[k][j][i] = 0;
      }
   }

   for (i=0; i<20; i++) {
       Ngnu[i]             = 0;
       Agnusum[i]          = 0.0;
       m_tot[i]            = 0.0;
       crystal_mass[i]     = 0.0;
       OptionZ[i]          = -1;
       REFL2[i].Z          = REFL[i].Z = -1;
       REFL2[i].CM         = REFL[i].CM = 0;
       REFL2[i].dist_geo   = REFL[i].dist_geo = 0.000000001;
       REFL2[i].fwhmt      = REFL[i].fwhmt    = 0.0;
       REFL2[i].minthick   = REFL[i].minthick = 0.0;
   }


   for (i=0; i<500; i++) {
      for (j=0; j<5; j++) T_Area[i][j] = 0.0;
   }
   for (i=0; i<3; i++) for (j=0; j<5; j++) A_chk[i][j] = 0.0;
   for (i=0; i<8; i++) for (j=0; j<20; j++) k_dist[i][j] = 0;
   for (i=0; i<100; i++) tri_wgt[i] = MATER[i].mass = 0.0;
   for (i=0; i<20; i++) { 	/* loop over crystal types (elem) */
      for (j=0; j<600; j++) {
         fracfi[i][j] = 0.0;   
         preset[i][j] = 0.0;
         preset[20][j] = 0.0;
         Eval_pos[i][j] = 0.0;
      }
      for (j=0; j<100; j++) {
         preset_i[i][j] = 0;
         chan[j] = 0;
      }
   }

   for (i=0; i<25; i++) panel3[i] = 0.0;

   znum       = 0;
   numOptions = 0;
   facet_diam = 4.0;

/* QQQQQ */

   printf("Usage: ../P2level Foc_len_factor fwhmt(arcmin) t_global detector_diam.(cm) [IDtext] Lens_layers\n");

   if (argc < 5) {
badargv: ;
      printf("Usage: ../P2level Foc_len_factor fwhmt(arcmin) t_global detector_diam.(cm) [IDtext] Lens_layers\n");
      exit(0);
   }
   system("\rm *.txt\n");
   
   Foc_len = 10000.0; /* basic focal length of 100 m */
   fwhmt = 1.0;
   minthick = 0.0001;
   detector_diam = 5.0;
   Lens_layers = 2;

/*   if (argc > 1) Foc_len = 100.0 * atof(argv[1]); */  /* convert to cm */
   if (argc > 1) {
      Foc_len *= atof(argv[1]);			/* scale Focal length */
      PassBandLow *= atof(argv[1]);		/* scale pass band */
      PassBandHigh *= atof(argv[1]);
      if (atof(argv[1]) > 5.0) goto badargv;
      GG = (int)(atof(argv[1]));
   }
   if (argc > 2) fwhmt = atof(argv[2]);
   if (argc > 3) t_global = atof(argv[3]);
   if (argc > 4) detector_diam = atof(argv[4]);
   if (argc > 5) {
      if (argc > 6) Lens_layers = atoi(argv[6]); 
      sprintf(IDtext, "_%s", argv[5]); 
      sprintf(zone_name, "%s_ZONES%1d.dat", argv[5], Lens_layers);
   } else {
      sprintf(IDtext, "");
      sprintf(zone_name, "ZONES.dat");
   }
   if (Lens_layers > 2) Lens_layers = 2;

   if (detector_diam < facet_diam) {
      printf("The detector diameter is not allowed to be less than the facet diameter! (%3.1f cm)\n", facet_diam); 
      exit(0);
   }
   if ((zone_file = fopen(zone_name, "rt")) == NULL) {
      printf("Could not open %s for reading!\n", zone_name);
      exit(0);
   }

   /*   Initialize list of reflection options */

   sprintf(Outtitl, "LaueW%03dT%03dR%02d%s_G%1d.out", 
      (int)(fwhmt*100.0+0.5), (int)(100.0*t_global+0.5), (int)(detector_diam*10.0+0.5), IDtext, GG);

   if ((out = fopen(Outtitl, "wt")) == NULL)  {
      printf("Could not open %s for writing!\n", Outtitl);
      exit(0);
   } else {
      printf("Opened %s for writing!\n", Outtitl);
   }

   fwhmt2 = fwhmt;
   T = 300.0;     /* assumed crystal temperature (K) (for thermal factor calculation) */ 
   
   num_zones = get_ZONES4(PZONES, zone_file);

   for (i=0; i<num_zones; i++) {
      for (k=0; k<numLayers; k++) {
         ZONES[i].ftlay[k][0] *= t_global; 
         ZONES[i].ftlay[k][1] *= t_global; 
         ZONES[i].ftlay[k][2] *= t_global; 
         ZONES[i].secftlay[k][0] *= t_global; 
         ZONES[i].secftlay[k][1] *= t_global; 
         ZONES[i].secftlay[k][2] *= t_global; 
      }
   }

   basic_Z = 0; 	/* the lens diameter is defined from the focal length and the 'basic_Z' material properties */
   znum = 0;
   for (i=0; i<num_zones; i++) {
      for (k=0; k<numLayers; k++) {
         for (iaz=0; iaz<3; iaz++) {
            if (ZONES[i].Z_lay[k][iaz] > 0) {
               for (j=0; j<znum; j++) if (IZ[j] == ZONES[i].Z_lay[k][iaz]) goto next_k;   
               IZ[znum++] = ZONES[i].Z_lay[k][iaz];
               if (basic_Z == 0) basic_Z = ZONES[i].Z_lay[k][iaz];	/* the first non-zero Z is taken as 'basic_Z' */
            }
         }
next_k: ;
      }
   }
   
   numOptions = 0;
   for (i=0; i<znum; i++) {
      Z = IZ[i];
      if (Z == 0) continue;
      old_numOp = numOptions;
      
      numOptions = initREFL(Z, &(MATER[Z]), Cgeo, atoms_eff, atoms_cell, PREFL, PZONES, num_zones, numOptions, fwhmt, minthick, OptionZ); 
      if (numOptions <= 0) {
         numOptions = old_numOp;
         printf("Found no data for Z: %2d\n", Z);
      }
/*
      for (j=old_numOp; j<numOptions; j++) {
         printf("After initREFL: znum: %1d, Z: %2d, numOptions: %2d, nOpt: %2d, REFL.Z: %2d, Miller: %s, F_cell: %5.1f, d-spacing: %5.3f\n", 
            znum, Z, numOptions, j, REFL[j].Z, REFL[j].Miller_txt, 
            MATER[REFL[j].Z].FormFact[REFL[j].CM][1], REFL[j].dist_geo*1.0e+8); 
      }
*/
      fprintf(out, "After initREFL: znum: %1d, Z: %2d, numOptions: %2d, nOpt: %2d, REFL.Z: %2d, Miller: %s, F_cell: %5.1f\n", 
         znum, Z, numOptions, old_numOp, REFL[old_numOp].Z, REFL[old_numOp].Miller_txt, 
         MATER[REFL[old_numOp].Z].FormFact[REFL[old_numOp].CM][1]); 
   }

   if (numOptions <= 0) {
      printf("Found no valid crystal elements\n");
      return(0);
   }

   sprintf(Titl, "F%03dW%03dT%03dR%02d%s", (int)(Foc_len*0.01), 
      (int)(fwhmt*100.0+0.5), (int)(100.0*t_global+0.5), (int)(detector_diam*10.0+0.5), IDtext);

   printf( "\n R_out R_in   Z10 Mil10 Thk10   Z11 Mil11 Thk11   Z12 Mil12 Thk12   Z10 Mil10 Thk10   Z11 Mil11 Thk11   Z12 Mil12 Thk12 \n");   
   for (j=0; j<num_zones; j++) {
      printf("L0: %5.1f %5.1f  %2d %3d %4.2f   %2d %3d %4.2f   %2d %3d %4.2f \n",
         ZONES[j].fR_outer, ZONES[j].fR_inner,
         ZONES[j].Z_lay[0][0], ZONES[j].Millay[0][0], ZONES[j].ftlay[0][0], 
         ZONES[j].Z_lay[0][1], ZONES[j].Millay[0][1], ZONES[j].ftlay[0][1], 
         ZONES[j].Z_lay[0][2], ZONES[j].Millay[0][2], ZONES[j].ftlay[0][2]);
      printf( "secL0 Z00 Mil00 Thk00 FrA00 FrB00  Z01 Mil11 Thk01 FrA01 FrB01  Z02 Mil02 Thk02 FrA02 FrB02\n");   
      printf(" %2d %3d %4.2f %4.2f %4.2f   %2d %3d %4.2f %4.2f %4.2f   %2d %3d %4.2f %4.2f %4.2f \n", 
         ZONES[j].secZ_lay[0][0], ZONES[j].secMillay[0][0], ZONES[j].secftlay[0][0], ZONES[j].secfracA[0][0], ZONES[j].secfracB[0][0],
         ZONES[j].secZ_lay[0][1], ZONES[j].secMillay[0][1], ZONES[j].secftlay[0][1], ZONES[j].secfracA[0][1], ZONES[j].secfracB[0][1],
         ZONES[j].secZ_lay[0][2], ZONES[j].secMillay[0][2], ZONES[j].secftlay[0][2], ZONES[j].secfracA[0][2], ZONES[j].secfracB[0][2]);
      printf("L1: %5.1f %5.1f  %2d %3d %4.2f   %2d %3d %4.2f   %2d %3d %4.2f \n",
         ZONES[j].Z_lay[1][0], ZONES[j].Millay[1][1], ZONES[j].ftlay[1][2],
         ZONES[j].Z_lay[1][0], ZONES[j].Millay[1][1], ZONES[j].ftlay[1][2],
         ZONES[j].Z_lay[1][0], ZONES[j].Millay[1][1], ZONES[j].ftlay[1][2]);
      printf( "secL1 Z10 Mil10 Thk10 FrA10 FrB10  Z11 Mil11 Thk11 FrA11 FrB11  Z12 Mil12 Thk12 FrA12 FrB12\n");   
      printf(" %2d %3d %4.2f %4.2f %4.2f   %2d %3d %4.2f %4.2f %4.2f   %2d %3d %4.2f %4.2f %4.2f \n", 
         ZONES[j].secZ_lay[1][0], ZONES[j].secMillay[1][0], ZONES[j].secftlay[1][0], ZONES[j].secfracA[1][0], ZONES[j].secfracB[1][0],
         ZONES[j].secZ_lay[1][1], ZONES[j].secMillay[1][1], ZONES[j].secftlay[1][1], ZONES[j].secfracA[1][1], ZONES[j].secfracB[1][1],
         ZONES[j].secZ_lay[1][2], ZONES[j].secMillay[1][2], ZONES[j].secftlay[1][2], ZONES[j].secfracA[1][2], ZONES[j].secfracB[1][2]);
   }
   printf("\n");

   panel3[0] = Foc_len;

   Bragg_max = atan((12.39 / PassBandLow) / (2.0 * REFL[0].dist_geo * 1.0e+8)) / GTORAD;
   Bragg_min = atan((12.39 / PassBandHigh) / (2.0 * REFL[0].dist_geo * 1.0e+8)) / GTORAD;
   lens_diam = 2.0 * 2.0 * Foc_len * tan(Bragg_max * GTORAD);
   lens_inner_diam = 2.0 * 2.0 * Foc_len * tan(Bragg_min * GTORAD);
   printf("Lens diameter for a Foc_len of %6.1f m and Z = %2d: outer: %6.1f cm, inner: %6.1f cm, Bragg angles: %5.3f %5.3f\n", 
      Foc_len / 100.0, basic_Z, lens_diam, lens_inner_diam, Bragg_max, Bragg_min);  

   panel3[1] = lens_diam;  

   for (j=0; j<num_zones; j++) { 	/* convert fractional radii to true cm */
      ZONES[j].R_outer = 0.5 * lens_diam * ZONES[j].fR_outer / 100.0;
      ZONES[j].R_inner = 0.5 * lens_diam * ZONES[j].fR_inner / 100.0;
   }

   facet_diam    = 4.0;
   facet_radius = 0.5 * facet_diam;
   panel3[2] = facet_diam;  

   MAXcm = (int)((lens_diam + facet_diam + 1.0) / 2.0);
   det_area = PI * detector_diam * detector_diam / 4.0;

   for (i=0; i<numOptions; i++) {
       REFL2[i].Z          = REFL[i].Z;
       REFL2[i].CM         = REFL[i].CM;
       REFL2[i].dist_geo   = REFL[i].dist_geo;
       REFL2[i].fwhmt      = REFL[i].fwhmt  ;
       REFL2[i].minthick   = REFL[i].minthick;
       sprintf(REFL2[i].Miller_txt,"%s", REFL[i].Miller_txt);
       
       REFL[i].dist_geo = 0.000000001;
/*
       sprintf(gnuFILname, "option%02d_00.gnu", i);
       gnuFIL[0][i] = fopen(gnuFILname, "wt");
       sprintf(gnuFILname, "option%02d_01.gnu", i);
       gnuFIL[1][i] = fopen(gnuFILname, "wt");
       sprintf(gnuFILname, "option%02d_02.gnu", i);
       gnuFIL[2][i] = fopen(gnuFILname, "wt");
*/
   }

   for (i=0; i<numOptions; i++) {
      thetaB = 0.5 * 0.5 * lens_diam / Foc_len;	/* Bragg angle at edge of lens for chosen focal length */
      lambda = thetaB * 2.0 * REFL2[i].dist_geo * 1.0e+8;
      E = 12.39 / lambda;				/* lambda in Angstroem, E in KeV */
      dE = E * REFL2[i].fwhmt * arcmin2rad / thetaB;
      printf("Option: %2d, Z: %2d, Element: %s, Miller Index: %3s, d_spacing: %8.6lf (AA), Eedge: %6.1f, dE: %5.2f, FormFact: %4.2f\n",
         i, REFL2[i].Z, MATER[REFL2[i].Z].ELEM_name, REFL2[i].Miller_txt, REFL2[i].dist_geo*1.0e+8, E, dE, 
         MATER[REFL2[i].Z].FormFact[REFL2[i].CM][1]);
      fprintf(out, 
         "Option: %2d, Z: %2d, Element: %s, Miller Index: %3s, d_spacing: %8.6lf (AA), Eedge: %6.1f, dE: %5.2f, FormFact: %4.2f\n",
         i, REFL2[i].Z, MATER[REFL2[i].Z].ELEM_name, REFL2[i].Miller_txt, REFL2[i].dist_geo*1.0e+8, E, dE, 
         MATER[REFL2[i].Z].FormFact[REFL2[i].CM][1]);
      E_ref[i] = E;
   }

/* The facets in layer 1 of the panel is offset from those in layer 0 by the following displacements:
 *       delta-x = 0.5 * facet_diam;
 *       delta-y = delta-x * tan(30deg);
 * The facets in layer 2 of the panel is offset from those in layer 0 by the following displacements:
 *       delta-x = 0.5 * facet_diam;
 *       delta-y = -delta-x * tan(30deg);
*/

   x_focus0 = 1.0 / 7.0;
   y_focus0 = 0.3333;         	/* Displacement of the first facet pattern with respect to the telescope axis */

   x_focus1 = x_focus0 - 0.5;   /* displacement of first additional layer facet pattern rel. to telescope axis */
   y_focus1 = y_focus0 - 0.5 * tan(30.0 * GTORAD);   
   x_focus2 = x_focus0 - 0.5;   /* displacement of second additional layer facet pattern rel. to telescope axis */
   y_focus2 = y_focus0 + 0.5 * tan(30.0 * GTORAD);   
   
   panel3[3] = x_focus0;	/* displacement of basic layer facet pattern rel. to telescope axis */
   panel3[4] = y_focus0;
      
   panel3[5] = x_focus1;	/* displacement of first additional layer facet pattern rel. to telescope axis */
   panel3[6] = y_focus1;
      
   cameraH = 11.0;		/* Camera optical axis, height above facet plane */
   panel3[7] = cameraH;
      
   cameraD =  8.0;		/* Diameter of camera optic */
   panel3[8] = cameraD;

   MirrorDispl = 0.8;           /* Horizontal displacement of alignment mirrors */
   panel3[9] = MirrorDispl;

   MirrorWidth = 0.5;           /* width of alignment mirrors */
   panel3[10] = MirrorWidth;

   LensTopRadius = 296.0;	/* Max radius for upward extension of lens (relevant for vertical launch configuration) */
   panel3[11] = LensTopRadius;

   facet_hole_radius = 0.4;
   panel3[12] = facet_hole_radius;

   panel3[13] = detector_diam / 2.0;

   facet_sep     = facet_diam + 0.25;	/* facts are separated 2.5 mm from their neighbors in the hexagonal pattern */
   panel3[14] = facet_sep;  		/* this is the minimum separation which leaves enough room for the 8 mm diam. pedestals */
   ring_width = facet_sep;

   panel3[15] = lens_inner_diam - facet_sep;  
   panel3[1]  = panel3[1] + facet_sep;

   panel3[16] = x_focus2;	/* displacement of second additional layer facet pattern rel. to telescope axis */
   panel3[17] = y_focus2;
      
   /* **************************************************************************** */
   /*                                                                              */
   /*       QQQ1 Generate LIST with crystal positions on each panel layer          */
   /*            Note: Not all crystal positions are used !                              */
   /*                                                                              */
   /* **************************************************************************** */
   
   good_facets = 0;

   for (layer_no=0; layer_no<numLayers; layer_no++) {
      facet_no[layer_no] = 1;
      lastBragg = 0.0;
      first_ring = last_good = 0;
      maxRing = lens_diam / (2.0 * facet_sep * tan(30.0 * GTORAD)) + 1.0;

      INITfacets = (1+maxRing)*maxRing*3+1; 	/* number of facets in hexagnal field surrounding the lens */

      if (maxRing > maxRings) {
         printf("The required maxRing-value, %3d,  exceeds the allowed value: %3d.  Increase maxRings in panel.h \n", maxRing, maxRings);
         exit(0);
      } else {
/*         printf("Required maxRing-value, %3d,  Allowed value: %3d.  INITfacets: %6d, N_CryPos: %6d, numLayers: %1d, facets: %6d\n", 
            maxRing, maxRings, INITfacets, N_CryPos, numLayers, numFACETS);
*/
      }

      jold = 0;
      for (i=1; i<maxRing; i++) { 
         facet_init[layer_no] = m = facet_no[layer_no];
         FACET[0 + layer_no * N_CryPos].nOption = -1;
         facet_no[layer_no] = hexring(i, layer_no, facet_no[layer_no], panel3, FACET);
         j = 0;
         for (k=0; k<facet_no[layer_no]; k++) if (FACET[k + layer_no * N_CryPos].nOption >= 0) j++;
/*
         printf("p2level: FACET: ring: %3d, INITfacets: %8d, index: %4d  %4d, useful: %3d  %3d\n", 
            FACET[facet_no[layer_no]-1 + layer_no * N_CryPos].ring_num, INITfacets, facet_no[layer_no], 
            facet_no[layer_no] - m, j, j-jold);
*/
         jold = j;
      }

      for (j=0; j<facet_no[layer_no]; j++) {
         if (FACET[j + layer_no * N_CryPos].nOption >= 0) {
            good_facets++; 
            lastBragg = FACET[j + layer_no * N_CryPos].thetaB;
            lastR     = FACET[j + layer_no * N_CryPos].r;
            if (j >= facet_init[layer_no]) {
               last_ring = i; 
               last_init[layer_no]=facet_init[layer_no]; 
            } 
            if (j < last_init[layer_no]) last_good = good_facets;
         }
         if (good_facets == 0) i = 5000; else first_ring = i;
      }

      jjj = 0;
      for (iii=0; iii<facet_no[layer_no]; iii++) {
         if (FACET[iii + layer_no * N_CryPos].nOption >= 0) {jjj++;}
      }

      printf("p2level: After hexring: layer %1d, Max ring %2d facets %5d %4d, valid facets (nOpt>=0): %4d, total good %6d\n", 
         layer_no, maxRing, INITfacets, facet_no[layer_no], jjj, good_facets);
 
      pr_ok = 0;
      LEN[layer_no] = facet_no[layer_no];
   }		/* end of loop over 'layer_no' (start 452) */
   /* QQQQQ2 */

   LEN[numLayers] = 0;	/* LEN[numLayers] is reserved for 'numOptions', to be used after facet initialization */

   for (i=0; i<MAXcm; i++) distr_layer[0][i] = distr_layer[1][i] = distr_layer[2][i] = 0;

   /* ***************************************************************** */
   /*                                                                   */
   /*  Order facets according to their distance from telescope axis     */
   /*                                                                   */
   /* ***************************************************************** */

   for (j=0; j<numLayers; j++) {
      uselen[j] = facet_order(FACET, L_order[j], distr_layer[j], facet_no[j], ring_width, j); 
   }

   printf("Ordered %5d %5d facets on the two crystal layers\n ", uselen[0], uselen[1]);

   printf("Distribution of facets as function of distance (cm) from lens axis\n");
   for (i=(int)(lens_inner_diam * 0.5 - (int)(lens_inner_diam * 0.5)%2 - 5); i< (int)(lens_diam * 0.5 + 5); i+=10) {
      sprintf(outline2, "%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d          ", 
         distr_layer[0][i+0], distr_layer[0][i+1], distr_layer[0][i+2], distr_layer[0][i+3], distr_layer[0][i+4],
         distr_layer[0][i+5], distr_layer[0][i+6], distr_layer[0][i+7], distr_layer[0][i+8], distr_layer[0][i+9]);
      sprintf(outline3, "%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d", 
         distr_layer[1][i+0], distr_layer[1][i+1], distr_layer[1][i+2], distr_layer[1][i+3], distr_layer[1][i+4],
         distr_layer[1][i+5], distr_layer[1][i+6], distr_layer[1][i+7], distr_layer[1][i+8], distr_layer[1][i+9]);
      printf("%2d  %s %s\n", i, outline2, outline3);
   }

   for (j=0; j<numLayers; j++) {
      uselen[j] = 0;
      for (i=1; i< MAXcm; i++) uselen[j] += distr_layer[j][i];
   }

   /* ************************************************************ */
   /*                                                              */
   /*            Assign zone and fzone values to facets            */
   /*          Count number of facets per radius interval          */
   /*                                                              */
   /* ************************************************************ */

   for (nl=0; nl<numLayers; nl++) {	/* Loop over layers */
      npr = 0;
      for (n_z=0; n_z<num_zones; n_z++) {
         printf("layer %1d, num_zones: %2d, zone[%1d]: Rin/out: %6.1f %6.1f (cm), Opt_lay: prim: %2d %2d %2d, sec: %2d %2d %2d\n", 
            nl, num_zones, n_z, ZONES[n_z].R_inner, ZONES[n_z].R_outer, 
            ZONES[n_z].Opt_lay[nl][0], ZONES[n_z].Opt_lay[nl][1], ZONES[n_z].Opt_lay[nl][2],
            ZONES[n_z].secOpt_lay[nl][0], ZONES[n_z].secOpt_lay[nl][1], ZONES[n_z].secOpt_lay[nl][2]);
      }

      for (ii=0; ii<LEN[nl]; ii++) {	/* Loop over facets */
         if (L_order[nl][LEN[nl] - 1 - ii] < 0) continue;
         iseq = L_order[nl][LEN[nl] - 1 - ii] + nl * N_CryPos;

         if (FACET[iseq].Z < 0) continue;	/* this facet is outside the physical panel area */

         FACET[iseq].nOption = -1;	/* initialize nOption for this facet */
         for (n_z=0; n_z<num_zones; n_z++) {
            if ((FACET[iseq].r >= ZONES[n_z].R_inner) && (FACET[iseq].r < ZONES[n_z].R_outer)) {
               FACET[iseq].zone  = n_z;
               FACET[iseq].fzone = (ZONES[n_z].R_outer - FACET[iseq].r) / (ZONES[n_z].R_outer - ZONES[n_z].R_inner);
               FACET[iseq].azim  = atan2(FACET[iseq].y, FACET[iseq].x) / GTORAD;
               if (FACET[iseq].azim < 0.0) FACET[iseq].azim += 360.0;
               azim_index = (int)(FACET[iseq].azim / 120.0);
               if (ZONES[n_z].Z_lay[nl][azim_index] == -1) ZONES[n_z].Opt_lay[nl][azim_index] = -1;
               if (ZONES[n_z].secZ_lay[nl][azim_index] == -1) ZONES[n_z].secOpt_lay[nl][azim_index] = -1;

               FACET[iseq].nOption = ZONES[n_z].Opt_lay[nl][azim_index];	/* define nOption for this facet */
               FACET[iseq].FTlay = ZONES[n_z].ftlay[nl][azim_index];

               use_frac = ZONES[n_z].secfracB[nl][azim_index] * (1.0 - FACET[iseq].fzone)
                        + ZONES[n_z].secfracA[nl][azim_index] * FACET[iseq].fzone;
               lo_azim_lim = azim_index * 120.0 - use_frac * 60.0 + 60.0;
               hi_azim_lim = azim_index * 120.0 + use_frac * 60.0 + 60.0;

               if (npr++ < 30) {
                  printf(
   "%1d %3d FAC %6d Opt %2d scOpt %2d az %5.1f idx %1d r %5.1f fzone %5.3f FrA/B %4.2f %4.2f useF %5.3f lo/hilim %5.1f %5.1f fnOpt ",
                     nl, ii, iseq, FACET[iseq].nOption, ZONES[n_z].secOpt_lay[nl][azim_index], FACET[iseq].azim, azim_index,
                     FACET[iseq].r, FACET[iseq].fzone, ZONES[n_z].secfracA[nl][azim_index], ZONES[n_z].secfracB[nl][azim_index],
                     use_frac, lo_azim_lim, hi_azim_lim);
               }
               if ((lo_azim_lim <= FACET[iseq].azim) && (FACET[iseq].azim < hi_azim_lim)) {
                  FACET[iseq].nOption = ZONES[n_z].secOpt_lay[nl][azim_index];	/* redefine nOption for this facet */
                  FACET[iseq].FTlay = ZONES[n_z].secftlay[nl][azim_index];
                  N_sec[nl][azim_index]++;
               }
               if (npr < 31) printf("%2d %3d %3d %3d\n", FACET[iseq].nOption, N_sec[nl][0], N_sec[nl][1], N_sec[nl][2]);
            }
         }
      }
   }

   /* ***************************************************************** */
   /*                                                                   */
   /*  Establish neighbor connections between the crystal layers        */
   /*                                                                   */
   /* ***************************************************************** */

   facet_R = facet_radius * 1.5;
   allnabs = facet_connect3(PFACET, uselen, L_order, facet_R, Lens_layers);

   printf("After neighbor check, allnabs: %5d\n", allnabs);

   for (i=0; i<numOptions; i++) {
      thetaB = 0.5 * 0.5 * lens_diam / Foc_len;	/* Bragg angle at edge of lens for chosen focal length */
      lambda = thetaB * 2.0 * REFL2[i].dist_geo * 1.0e+8;
      E = 12.39 / lambda;				/* lambda in Angstroem, E in KeV */
      dE = E * REFL2[i].fwhmt * arcmin2rad / thetaB;
/*    printf("p3lA Option: %2d, Z: %2d, Element: %s, Miller Index: %3s, d_spacing: %8.6lf (AA), Eedge: %6.1f, dE: %5.2f, FormFact: %4.2f\n",
         i, REFL2[i].Z, MATER[REFL2[i].Z].ELEM_name, REFL2[i].Miller_txt, REFL2[i].dist_geo*1.0e+8, E, dE, 
         MATER[REFL2[i].Z].FormFact[REFL2[i].CM][1]);
*/
   }

   /* **************************************************************************************** */
   /*                                                                                          */
   /*          QQQ3  Initiate all possible facet reflection options for all facets             */
   /*                                                                                          */
   /* **************************************************************************************** */

newFF1:
   printf("before init_facet_REFL: LEN %6d %6d numOpts %3d 1st facet opts, first facet numbers: %6d %6d  options: %2d %2d\n",
      LEN[0], LEN[1], numOptions, L_order[0][0], L_order[1][0], FACET[L_order[0][0]].nOption, FACET[(L_order[1][0]+N_CryPos)].nOption);
   printf("before init_facet_REFL: LEN %6d %6d numOpts %3d 1st facet opts, layer0/1 %2d %2d, Z prim %2d %2d sec %2d %2d r0/1: %5.1f %5.1f\n", 
      LEN[0], LEN[1], numOptions, FACET[L_order[0][0]].nOption, FACET[(L_order[1][0]+N_CryPos)].nOption,
      ZONES[0].Z_lay[0][0], ZONES[0].Z_lay[1][0], ZONES[0].secZ_lay[0][0], ZONES[0].secZ_lay[1][0], 
      FACET[L_order[0][0]].r, FACET[L_order[1][0]+N_CryPos].r);

               if (ZONES[n_z].secZ_lay[nl][azim_index] == -1) ZONES[n_z].secOpt_lay[nl][azim_index] = -1;
   init_facet_REFL(PREFL2, PMATER, FACET, uselen, L_order, numOptions); 

   LEN[numLayers] = numOptions;

   /* *******************************************************************''''''''''''''''''''''''''''''''*** */
   /*                                                                                                        */
   /*                 POPULATE facets with selected crystals according to Zone specifications                */
   /*                                  Start from lowest diffracted energy                                   */
   /*                                                                                                        */
   /* *******************************************************************''''''''''''''''''''''''''''''''*** */

   A_basic = PI * facet_radius*facet_radius;		/* Area within crystal perimeter cm2 */
   A_hole  = PI * facet_hole_radius*facet_hole_radius;	/* Area within crystal hole cm2 */
   A_basic -= A_hole;

/* QQQQQ begin1 */

   sprintf(Titl, "F%03dW%03dT%03dR%02d%s", (int)(Foc_len*0.01), 
      (int)(fwhmt*100.0+0.5), (int)(100.0*t_global+0.5), (int)(detector_diam*10.0+0.5), IDtext);

   for (nl=0; nl<numLayers; nl++) {		/* begin loop over layers QQQQQ (end in 907) */
      printf("Begin layer %1d\n", nl);

      sprintf(thktitl, "thk%s_L%1d.txt", IDtext, nl);
      if (GG == 1) {
         if ((thkfil = fopen(thktitl, "wt")) != NULL) {
            for (ii=0; ii<LEN[nl]; ii++) {		/* begin short loop over facets to initialize thkarray */
               iseq = L_order[nl][LEN[nl] - 1 - ii] + nl * N_CryPos;
               if (L_order[nl][LEN[nl] - 1 - ii] < 0) continue;
               thkarray[nl][iseq] = 0.0;
            }
            GGOK = 1; 
            printf("Layer%1d, open %s for write, set %6d elements of 'thkarray' to zero, set GGOK to %1d\n", 
               nl, thktitl, LEN[nl], GGOK);
            fprintf(out, "Layer%1d, open %s for write, set %6d elements of 'thkarray' to zero, set GGOK to %1d\n", 
               nl, thktitl, LEN[nl], GGOK);
         } else {
            printf("Could not open  %s for writing\n", thktitl);
            exit(0);
         }
      } else {
         if ((thkfil = fopen(thktitl, "rt")) != NULL) {
            GGOK = 1;	/* stored thk-data found */
            printf("Opened %s for reading of thkarray\n", thktitl);
            thk_sum = 0.0;
            i = 0;
            for (ii=0; ii<LEN[nl]; ii++) {		/* begin short loop over facets to read thkarray from thk-file */
               iseq = L_order[nl][LEN[nl] - 1 - ii] + nl * N_CryPos;
               if (L_order[nl][LEN[nl] - 1 - ii] < 0) continue;
               i++;
               fgets(iline, 16, thkfil);
               if (feof(thkfil)) {
                  printf("Premature eof on %s - after reading %5d elements\n", thktitl, i);
                  exit(0);
               }
               sscanf(iline, "%f", &thkarray[nl][iseq]);
               thk_sum += thkarray[nl][iseq];
            }
            printf("Layer%1d, open %s for read, got %5d %5d elements of 'thkarray', thk_sum: %6.1f cm, ,set GGOK to %1d\n", 
               nl, thktitl, i, LEN[nl], thk_sum, GGOK);
            fprintf(out, "Layer%1d, open %s for read, got %5d %5d elements of 'thkarray', thk_sum: %6.1f cm, set GGOK to %1d\n", 
               nl, thktitl, i, LEN[nl], thk_sum, GGOK);
         } else {
            printf("Could not open  %s for reading\n", thktitl);
            GGOK = 0;   /* No previous thk-data available */
         }
      }

      j = 0;
      azimnum = 0;
      for (ii=0; ii<LEN[nl]; ii++) {		/* begin loop over facets QQQQQ (end in 817) */
         iseq = L_order[nl][LEN[nl] - 1 - ii] + nl * N_CryPos;
         if (L_order[nl][LEN[nl] - 1 - ii] < 0) continue;

         if (FACET[iseq].Z < 0) continue;	/* this facet is outside the physical panel area */

         z_ring = (int)(FACET[iseq].r / ring_width);

         N_z = FACET[iseq].zone;		/* Use the zone-specification assigned to each facet */ 
         nOpt = FACET[iseq].nOption;		/* POPULATE FACETS: Reflection option */

         if (nOpt < 0) continue;		/* this facet shall not be filled */

         if (FF > 1) goto newFF2;		/* on subsequent calculations with modified Foc_len, skip all FACET modifications */
						/* FACET[].thetaB is the only parameter which need to be modified */ 
                                                /* specifically the facet thickness should be kept constant for all Foc_len-values */

         T_fclayer = FACET[iseq].FTlay;

         FACET[iseq].Z       = REFL2[nOpt].Z;	/* POPULATE FACETS: Material */

         Rgnu[nl][nOpt][(int)(FACET[iseq].r/facet_sep + 0.5)]++;

/*         fprintf(gnuFIL[0][nOpt], "%6.2f  %6.2f\n", FACET[iseq].x, FACET[iseq].y);
*/
	 /* get crystal thickness as required by nOpt and thickness factor specific for this layer and option */

         status = reflectT2(&(REFL2[nOpt]), &(MATER[FACET[iseq].Z]), &(FACET[iseq]), T_fclayer);

         if (status < 0.0) return(0.0);

         if (GGOK) {
            if (GG == 1) {
               thkarray[nl][iseq] = FACET[iseq].thick;	/* save this thickness if we are analyzing primary Foc_len choice */
            } else {
               FACET[iseq].thick = thkarray[nl][iseq];	/* use stored thickness values when analyzing subsequent Foc_len values */
            }
         }

         E = FACET[iseq].E[nOpt];	/* nominal Bragg energy for this facet */

         Efound[nOpt][0] = 1;
         IntE = (int)(E / 50.0);
         ModE = E - (float)(IntE * 50.0); 
         if ((ModE < 5.0) && (Efound[nOpt][IntE] == 0)) {
            Efound[nOpt][IntE]++;

            status = reflect3A(IntE*50.0, &(REFL2[nOpt]), &(FACET[iseq]), &(MATER[REFL2[nOpt].Z]), 1, &refl_E);

            if (status < 0.0) {printf("Bad call to reflect3A from p2level\n"); exit(0);}
            printf("Peak reflectivity for %s %s at %5d keV is %5.1f (%%) at a crystal thickness of %5.1f (mm)\n",
               MATER[FACET[iseq].Z].ELEM_symb, REFL2[nOpt].Miller_txt, IntE*50, refl_E*100.0, FACET[iseq].thick*10.0);
         }
         if ((ModE > 45.0) && (Efound[nOpt][IntE+1] == 0)) {
            Efound[nOpt][IntE+1]++;

            status = reflect3A((IntE+1)*50.0, &(REFL2[nOpt]), &(FACET[iseq]), &(MATER[REFL2[nOpt].Z]), 1, &refl_E);

            if (status < 0.0) {printf("Bad call to reflect3A from p2level\n"); exit(0);}
            printf("Peak reflectivity for %s %s at %5d keV is %5.1f (%%) at a crystal thickness of %5.1f (mm)\n",
               MATER[FACET[iseq].Z].ELEM_symb, REFL2[nOpt].Miller_txt, (IntE+1)*50, refl_E*100.0, FACET[iseq].thick*10.0);
         }

         j++;
newFF2:
   /* *************************************************************''''''''''''''''''''''''''''''''*** */
   /*                                                                                                  */
   /*     QQQQQ               Begin integration of effective area of lens                              */
   /*                                                                                                  */
   /* *******************************************************************''''''''''''''''''''''''''*** */

         FACET[iseq].PEAKrefl = 0.0;
         FACET[iseq].PEAK_E = 0.0;

   /* *************************************************************''''''''''''''''''''''''''''''''*** */
   /*                                                                                                  */
   /*     QQQQQ   Loop over reflection orders (1-4) and over the energy interval relevant for the      */
   /*                    current facet to integrate the effective area of the lens                     */
   /*                                                                                                  */
   /* *******************************************************************''''''''''''''''''''''''''*** */

         for (iOrder = 1; iOrder < 4; iOrder++) {      /* loop over reflection orders, ends in 799 */
            E_low  = (FACET[iseq].E[nOpt] - FACET[iseq].dE[nOpt]) * (float)(iOrder);
            E_high = (FACET[iseq].E[nOpt] + FACET[iseq].dE[nOpt]) * (float)(iOrder);
            iE_low  = (int)((E_low - 1.0) / 2.0);
            if (iE_low > 2999) iE_low = 2999;
            if (iE_low < 10  ) iE_low = 10  ;
            iE_high = (int)((E_high + 1.0) / 2.0);
            if (iE_high > 2999) iE_high = 2999;

            for (iE = iE_low; iE <= iE_high; iE++) {	/* Loop over E in 2 keV steps (including higher order E-ranges) */
               E = 2.0 * (float)(iE);

               B_contr = facet_Bcontrib(E, nOpt, iOrder, panel3, &(FACET[0]),
                                     &(MATER[REFL2[nOpt].Z]), &(REFL2[nOpt]), LEN, L_order, &refl_E, iseq, NABSTAT, gnufil2, Lens_layers);

               EffArea[nl][N_z][iOrder][iE] += B_contr;
               EffAOpt[nl][nOpt][iOrder][iE] += B_contr;
               Agnu[nOpt][iE] += B_contr;
               if ((iOrder == 0) && (fabs(FACET[iseq].E[0] - E) < 2.0)) FACET[iseq].PEAKrefl = refl_E; 
            }
            if (FF > 1) continue;

            if ((iOrder == 1) && (B_contr > 0.20)) {
               mE = (int)(FACET[iseq].E[nOpt] / 2.0 + 0.5);
               EffArea[nl][N_z][0][mE] = FACET[iseq].thick; 	/* EffArea[iOrder 0] stores the facet thickness at this energy */
               EffAOpt[nl][nOpt][0][mE] = FACET[iseq].thick; 	/* EffAOpt[iOrder 0] stores the facet thickness at this energy */
            }
         }		/* end of loop over reflection orders */
         if (FF > 1) continue;

         crystal_mass[nOpt] += A_basic * MATER[FACET[iseq].Z].dens * FACET[iseq].thick * 0.001;
         Ngnu[nOpt]++; 
         
         if ((nl == 0) && ((iseq%20000) == 0)) 
            printf("nl: %1d, N_z: %1d, ii: %5d, iseq: %5d, nOpt: %2d, Z: %2d, thick: %5.3f, Mill: %s, mass: %5.1f, Ngnu: %4d %4d %4d %4d %4d %4d\n",
                         nl, N_z, ii, iseq, FACET[iseq].nOption, FACET[iseq].Z,  FACET[iseq].thick, REFL2[FACET[iseq].nOption].Miller_txt,
                         A_basic * MATER[FACET[iseq].Z].dens * FACET[iseq].thick, Ngnu[0], Ngnu[1], Ngnu[2], Ngnu[3], Ngnu[4], Ngnu[5]);
      }  /* End of loop over facets QQQQQ */

      NNopt = 0;
      for (nOpt=0; nOpt<6; nOpt++) {
         if (FF == 1) {
            m_tot[nOpt]  = crystal_mass[nOpt];
            M_tot       += crystal_mass[nOpt];
         }
         if (Ngnu[nOpt] > 0) NNopt++;
      }
      for (n_z=0; n_z<num_zones; n_z++) {
         for (iaz=0; iaz<3; iaz++) {
            nOpt       = ZONES[n_z].Opt_lay[nl][iaz];

            if (nOpt >= 0)
               printf("%2d nOpt: %2d, iaz: %1d, Z: %2d Reference energy %5.1f keV, for %s(%s) t_fact: %4.2f, Mass: %5.1f  %6.1f kg\n", 
                  n_z, nOpt, iaz, REFL2[nOpt].Z, E_ref[nOpt], MATER[REFL2[nOpt].Z].ELEM_symb, 
                  REFL2[nOpt].Miller_txt, ZONES[n_z].ftlay[nl][iaz], crystal_mass[nOpt], M_tot);
            }
      }

      if (GGOK) {
         if (GG == 1) {
            thk_sum = 0.0;
            for (ii=0; ii<LEN[nl]; ii++) {		/* begin short loop over facets QQQQQ */
               iseq = L_order[nl][LEN[nl] - 1 - ii] + nl * N_CryPos;
               if (L_order[nl][LEN[nl] - 1 - ii] < 0) continue;
               fprintf(thkfil, "%6.3f\n", thkarray[nl][iseq]);
               thk_sum += thkarray[nl][iseq];
            }
            printf("Layer%1d, %s, wrote %6d elements of 'thkarray', thk_sum: %6.1f cm\n", 
               nl, thktitl, LEN[nl], thk_sum);
            fprintf(out, "Layer%1d, %s, wrote %6d elements of 'thkarray', thk_sum: %6.1f cm\n", 
               nl, thktitl, LEN[nl], thk_sum);
         }
      }
      if (thkfil != NULL) {
         fclose(thkfil);
         thkfil = NULL;
      }
      printf("Completed layer %1d\n", nl);

   }     /* End of layer loop */

    for (i=0; i<100; i++) if (MATER[i].mass > 0.0) printf("%s  %6.1f kg\n", MATER[i].ELEM_symb, MATER[i].mass);
   printf("\n");
/* end4   */

   printf( "\n %s  Focal len: %3.0f m, Mos wdth: %3.0f arcs, Det. diam: %3.1f cm, T_glob: %4.2f, NZ: %2d\n", 
      date_txt, Foc_len * 0.01, fwhmt * 60.0, detector_diam, t_global, num_zones);

   printf("Ngnu: %5d %5d %5d %5d %5d %5d\n", Ngnu[0], Ngnu[1], Ngnu[2], Ngnu[3], Ngnu[4], Ngnu[5]);
   out2 = fopen("Laue_summary.dat", "at");

   fprintf(out2,
      "\n %s  Focal length: %3.0f m, Mosaic width: %3.0f arcsec, Detector diameter: %3.1f cm, T_global: %4.2f\n", 
      date_txt, Foc_len * 0.01, fwhmt * 60.0, detector_diam, t_global);

   fprintf(out2, "NABSTAT      refl.sum              area sum              neighbor distribution  \n");
   for (i=0; i<Lens_layers; i++) {
      fprintf(out2, "%3.0f   %6.1f %6.1f %6.1f   %6.1f %6.1f %6.1f  ", NABSTAT[i][0],
         NABSTAT[i][1], NABSTAT[i][2], NABSTAT[i][3], NABSTAT[i][4], NABSTAT[i][5], NABSTAT[i][6]);
      fprintf(out2, "%4.0f %4.0f %4.0f %4.0f    %4.0f %4.0f %4.0f %4.0f    %4.0f %4.0f %4.0f %4.0f     %4.0f %4.0f %4.0f %4.0f\n", 
         NABSTAT[i][ 7], NABSTAT[i][ 8], NABSTAT[i][ 9], NABSTAT[i][10],
         NABSTAT[i][11], NABSTAT[i][12], NABSTAT[i][13], NABSTAT[i][14],
         NABSTAT[i][15], NABSTAT[i][16], NABSTAT[i][17], NABSTAT[i][18],
         NABSTAT[i][19], NABSTAT[i][20], NABSTAT[i][21], NABSTAT[i][22]);
   }
   
   fclose(out2);

   sprintf(xlabel,"Energy (keV)");
   i = strlen(xlabel);
   for (j=0; j<60; j++) dtext3[j] = dtext4[j] = ' ';
   dtext3[60 - i/2] = 0;
   dtext4[40 - i/2] = 0;
   sprintf(Xlabel, "%s%s%s%s", dtext3, xlabel, dtext4, date_txt);

   tlnum = 0;
   sprintf(textlines[tlnum++],"Sum of all layers ");
   for (j=0; j<6; j++) sprintf(textlines[j+1],"%s %s ", MATER[OptionZ[j]].ELEM_symb,  REFL[j].Miller_txt);

   sprintf(gnutitl, "gnu_Laue%s.plt", Titl);
   if (FF == 1) {
      sprintf(gnutitl1, "gnu_Laue%s.plt", Titl);
      sprintf(PStitl, "Laue%s", Titl);
   }
   pltfil = fopen(gnutitl, "wt");

   if (Lens_layers == 3) sprintf(triptxt, "triple");
   else if (Lens_layers == 2) sprintf(triptxt, "double"); else sprintf(triptxt, "single");

   fprintf(pltfil,"set terminal postscript color\n");
   fprintf(pltfil,"set output \"Laue%s.ps\"\n", Titl);

   fprintf(pltfil,
      "set title \"Laue lens eff. area: F %3.0f m, Lens diam: %3.0f cm, Detect.diam. %3.1f cm, mosaic: %3.1f arcmin \"\n",
      Foc_len/100.0, lens_diam, detector_diam, fwhmt);
   fprintf(pltfil,"set xlabel \"%s\"\n", Xlabel);
   fprintf(pltfil,"set ylabel \"Effective Area (cm2) \"\n");
   fprintf(pltfil,"set y2label \"Gain \"\n");
   fprintf(pltfil,"set xrange [%4d: %4d]\n", 250*FF, 2800*FF);
   DetArea = PI * detector_diam * detector_diam * 0.25;
   y2min = 0;
   y2max = 800;
   fprintf(pltfil,"set yrange [%4d: %4d]\n", y2min, y2max);
   fprintf(pltfil,"set y2range [%2d:%2d]\n", (int)(y2min / DetArea)+1, (int)(y2max / DetArea)+1);
   fprintf(pltfil,"set y2tics\n");

   sprintf(Lauetitl, "Laue%s.gnu", Titl);
   gnufil = fopen(Lauetitl, "wt");
   for (j=0; j<numLayers; j++) Agnusum[j] = 0.0;

   if (NNopt == 1) strcpy(textlines[1], textlines[0]);

			 /* gnunum 0 */
   gnunum = 0;
   for (i=0; i<3000*FF; i++) fprintf(gnufil, "%6d   %5.1f\n", 
      i*2, Agnu[0][i]+Agnu[1][i]+Agnu[2][i]+Agnu[3][i]+Agnu[4][i]+Agnu[5][i]);
   fprintf(pltfil,
      "plot \"%s\" index %d axes x1y1 title \"%s #crystals %4d  M %5.1f kg\" with lines lw 2\n", 
       Lauetitl, gnunum, textlines[0], Ngnu[0]+Ngnu[1]+Ngnu[2]+Ngnu[3]+Ngnu[4]+Ngnu[5], 
       m_tot[0]+m_tot[1]+m_tot[2]+m_tot[3]+m_tot[4]+m_tot[5]);
   fprintf(gnufil, "\n\n"); gnunum++;

			 /* gnunum 1 2 3 4 5 6 */
   for (j=0; j<6; j++) {
      if ((Ngnu[j] == 0) || (NNopt == 1)) continue;
      for (i=0; i<1500*FF; i++) {
         fprintf(gnufil, "%6d   %5.1f\n", i*2, Agnu[j][i]);
         Agnusum[j] += Agnu[j][i];
      }
      fprintf(pltfil,
         "replot \"%s\" index %d axes x1y1 title \"%s #crystals %3d, Tf %3.2f %3.2f %3.2f, M %3.0f %3.0f %3.0f kg\" with lines lw 2\n", 
          Lauetitl, gnunum, textlines[j+1], Ngnu[j], ZONES[0].ftlay[j][0], ZONES[0].ftlay[j][1], ZONES[0].ftlay[j][2] , m_tot[j]);
      fprintf(gnufil, "\n\n"); gnunum++;
   }

   printf("Integrated reflectivity, FF: %1d, Options 0-5: %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f, total: %7.1f  (keVcm2) \n",
      FF, Agnusum[0]*2.0, Agnusum[1]*2.0, Agnusum[2]*2.0, Agnusum[3]*2.0, Agnusum[4]*2.0, Agnusum[5]*2.0, 
      (Agnusum[0] + Agnusum[1] + Agnusum[2] + Agnusum[3] + Agnusum[4] + Agnusum[5])*2.0); 

   for (i=0; i<6; i++) {
      Agnusum[i] = 0.0;
      for (j=0; j<3000; j++) Agnu[i][j] = 0.0;
   }

/*
   for (j=0; j<3; j++) {
      valsum = 0.0;
      for (i=0; i<500; i++) {
         val = EffAOpt[2][j][1][i]+EffAOpt[2][j][2][i]+EffAOpt[2][j][3][i] +  
               EffAOpt[3][j][1][i]+EffAOpt[3][j][2][i]+EffAOpt[3][j][3][i];  
         fprintf(gnufil, "%6d   %5.1f\n", i*2, val);
         valsum += val;
      }
      if (valsum > 0.0) fprintf(pltfil, "replot \"%s\" index %d axes x1y1 title \"%s \" with lines \n", 
          Lauetitl, gnunum, textlines[gnunum]);
      fprintf(gnufil, "\n\n"); gnunum++;
   }
*/
			 /* gnunum 4 5 6 */
/*
   for (nl=0; nl<numLayers; nl++) {
      for (j=0; j<3; j++) {
         for (i=0; i<100; i++) {
            fprintf(gnufil, "%5.0f   %5d\n", (float)(i)*facet_sep, Rgnu[nl][j][i]);
            Rgnu[nl][3][i] += Rgnu[nl][j][i];
         }
         fprintf(pltfil,
            "#replot \"%s\" index %d axes x1y1 title \"%s \" with lines \n", 
             Lauetitl, gnunum, textlines[gnunum]);
         fprintf(gnufil, "\n\n"); gnunum++;
      }
   }
*/
    			 /* gnunum 10 11 */
/*   for (j=0; j<2; j++) {
      for (i=0; i<100; i++) {
         fprintf(gnufil, "%5.0f   %5d\n", (float)(i)*facet_sep-1.0+2.0*j, Rgnu[j][3][i]);
      }
      fprintf(pltfil,
         "#replot \"%s\" index %d axes x1y1 title \"%s \" with points\n", 
          Lauetitl, gnunum, textlines[gnunum]);
      fprintf(gnufil, "\n\n"); gnunum++;
   }
*/
			 /* gnunum 12 */
/*
   for (i=1; i<3000; i++) fprintf(gnufil, "%6d   %5.1f\n", i*2, 180.0 * E_ref[0] / ((float)(i*2)));
   fprintf(pltfil,
      "#replot \"%s\" index %d axes x1y2 title \"%s \" with lines \n", 
       Lauetitl, gnunum, textlines[gnunum], t_global, fwhmt, m_tot[0]+m_tot[1]);
   fprintf(gnufil, "\n\n"); gnunum++;
*/
   goto finis_gnu;

   for (i=0; i<3000; i++) 
      fprintf(gnufil, "%6d   %5.1f\n", i*2, 
         EffArea[2][2][1][i]+EffArea[2][2][2][i]+EffArea[2][2][3][i] +  
         EffArea[2][3][1][i]+EffArea[2][3][2][i]+EffArea[2][3][3][i] +  
         EffArea[2][4][1][i]+EffArea[2][4][2][i]+EffArea[2][4][3][i]);  
   fprintf(gnufil, "\n\n"); gnunum++;

   for (i=0; i<3000; i++) 
      fprintf(gnufil, "%6d   %5.1f\n", i*2, EffArea[2][5][1][i]+EffArea[2][5][2][i]+EffArea[2][5][3][i]);  
   fprintf(gnufil, "\n\n"); gnunum++;

   for (i=0; i<3000; i++) 
      fprintf(gnufil, "%6d   %5.1f\n", i*2, EffArea[3][0][1][i]+EffArea[3][0][2][i]+EffArea[3][0][3][i]);  
   fprintf(gnufil, "\n\n"); gnunum++;

   for (i=0; i<3000; i++) 
      fprintf(gnufil, "%6d   %5.1f\n", i*2, 
         EffArea[3][1][1][i]+EffArea[3][1][2][i]+EffArea[3][1][3][i] +  
         EffArea[3][2][1][i]+EffArea[3][2][2][i]+EffArea[3][2][3][i] +  
/*   fprintf(gnufil, "\n\n"); gnunum++;
   for (i=0; i<3000; i++) 
      fprintf(gnufil, "%6d   %5.1f\n", i*2, 
*/
         EffArea[3][3][1][i]+EffArea[3][3][2][i]+EffArea[3][3][3][i]);  
   fprintf(gnufil, "\n\n"); gnunum++;

   for (i=0; i<3000; i++) 
      fprintf(gnufil, "%6d   %5.1f\n", i*2, EffArea[3][4][1][i]+EffArea[3][4][2][i]+EffArea[3][4][3][i]);  
   fprintf(gnufil, "\n\n"); gnunum++;

finis_gnu:
   fclose(gnufil);

   fclose(pltfil);

   printf("NABSTAT      refl.sum              area sum              neighbor distribution  \n");
   for (i=0; i<Lens_layers; i++) {
      printf("%3.0f   %6.1f %6.1f %6.1f   %6.1f %6.1f %6.1f  ", NABSTAT[i][0],
         NABSTAT[i][1], NABSTAT[i][2], NABSTAT[i][3], NABSTAT[i][4], NABSTAT[i][5], NABSTAT[i][6]);
      printf("%4.0f %4.0f %4.0f %4.0f    %4.0f %4.0f %4.0f %4.0f    %4.0f %4.0f %4.0f %4.0f     %4.0f %4.0f %4.0f %4.0f\n", 
         NABSTAT[i][ 7], NABSTAT[i][ 8], NABSTAT[i][ 9], NABSTAT[i][10],
         NABSTAT[i][11], NABSTAT[i][12], NABSTAT[i][13], NABSTAT[i][14],
         NABSTAT[i][15], NABSTAT[i][16], NABSTAT[i][17], NABSTAT[i][18],
         NABSTAT[i][19], NABSTAT[i][20], NABSTAT[i][21], NABSTAT[i][22]);
   }
   
   if (FF < 0) {
      FF++;
      if (argc > 1) Foc_len = 100.0 * FF * atof(argv[1]);  /* convert to cm */
      panel3[0] = Foc_len;

      for (nl=0; nl<numLayers; nl++) {		/* begin loop over layers QQQQQ (end in 816) */
         printf("Update thetaB for new focal length for layer %1d\n", nl);
         for (ii=0; ii<LEN[nl]; ii++) {		/* begin loop over facets QQQQQ (end in 800) */
            iseq = L_order[nl][LEN[nl] - 1 - ii] + nl * N_CryPos;
            if (L_order[nl][LEN[nl] - 1 - ii] < 0) continue;

            FACET[iseq].thetaB = 0.5 * (float)(atan(FACET[iseq].r / panel3[0]));	/* Define new thetaB in Radians */
         }
      }
               
      for (j=0; j<20; j++) {
         Agnusum[j] = 0.0;
         for (i=0; i<N_CryPos; i++) {
            for (k=0; k<numLayers; k++) {
               EffArea[k*2+0][j][0][i] = EffArea[k*2+0][j][1][i] = EffArea[k*2+0][j][2][i] = EffArea[k*2+0][j][3][i] = 0.0;
               EffArea[k*2+1][j][0][i] = EffArea[k*2+1][j][1][i] = EffArea[k*2+1][j][2][i] = EffArea[k*2+1][j][3][i] = 0.0;
            }
         }
      }
      for (i=0; i<3000; i++) {
         for (k=0; k<6; k++) Agnu[k][i] = 0.0;
      }
      for (i=0; i<Lens_layers; i++) for (j=0; j<30; j++) NABSTAT[i][j] = 0.0;
 
      goto newFF1;

   }
   printf("Nrand0/1: %6d %6d\n", Nrand0, Nrand1);
   printf("N_sec:  %6d %6d %6d   %6d %6d %6d \n", N_sec[0][0], N_sec[0][1], N_sec[0][2], N_sec[1][0], N_sec[1][0], N_sec[1][2]); 
   
   sprintf(systemtxt, "gnuplot %s\n", gnutitl1);
   printf("%s", systemtxt);
   system(systemtxt); 
   sprintf(systemtxt, "ps2pdf %s.ps\n", PStitl);
   printf("%s", systemtxt);
   system(systemtxt); 
   sprintf(systemtxt, "\\rm %s.ps\n", PStitl);
   printf("%s", systemtxt);
   system(systemtxt); 
   sprintf(systemtxt, "mv *.pdf /ftp/pub/nl/Laue\n" , PStitl);
   printf("%s", systemtxt);
   system(systemtxt); 
  
   return(0);
}