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

/* **************************************************************************************** */

/* structure characterizing a reflection 
struct reflect {
   int   Z;			 atomic number 
   int   CM;			 Miller index indicator 
   char  Miller_txt[10];	 Miller index text 
   float fwhmt;			 mosaic spread (arcmin) 
   float minthick;		 minimum thickness acceptable to use (0.1) 
   float dist_geo;		 crystal plane separation for this reflection (AAngstroem)
*/   
int initREFL(int ZZ, material *MATER, float Cgeo[5][5], float atoms_eff[5][5], int atoms_cell[],
             reflect REFL[], zones ZONES[], int n_zones, int nOptions, float fwhmt, float minthick, int OptZ[])
{
   static FILE *in;

   char iline[200];
   static char elem[100][6], Struc[100][10];
   static char Miller[5][4] = {"100", "110", "111", "200", "220"};
   char unit = 'b', errmsg[100];

   int i, j, k, err, type, len, Cmiller, CMILLER, n_order, n_z;
   int print_flag = 1, CM, Z, azimBin;;
   static int Cstruc[100];
   static int first = 1, type_use[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
   static int Optfirst[20];   
   double energy[10], T, E; 
   double phi, thetaB, mu, t, sigmaE, alpha, nypeak, R, x, y, F_cell, RR, V;
   double lambda, d, a_vol, kappa, B, L, sigmat, f_atom, sqrt2pi, thermal_f, SL;
   static double A[100], dens[100], a_dens[100], DebyeT[100];
   static double dist[100], melt[100];
   double tmax, nymax, re, re2, alfa, alfa1, fwhmrad, sigmarad, tt, rr;

   sqrt2pi = sqrt(TWOPI);
   re = 2.818e-13;	/* classical electron radius (cm) */
   re2 = re * re;
   
   T = 273.2;

   if (first) {		/* initialize element data and read from 'bragg_t.dat' file */
      first = 0;
      for (i=0; i<20; i++) Optfirst[i] = 1;
      for (i=0; i<100; i++) {
         sprintf(Struc[i], "%s", "   ");  
         sprintf(elem[i], "%s", "   ");  
         A[i] = dens[i] = a_dens[i] = DebyeT[i] = dist[i] = melt[i] = 0.0;
      }

      for (i=0; i<5; i++) {
         Cgeo[i][0] = Cgeo[i][1] = Cgeo[i][2] = Cgeo[i][3] = Cgeo[i][4] = -1.0;
         atoms_eff[i][0] = atoms_eff[i][1] = atoms_eff[i][2] = atoms_eff[i][3] = atoms_eff[i][4] = 0.0;
      }
      Cgeo[0][2] = 0.57735;    /* plane separation factor for 'fcc 111' */
      Cgeo[0][3] = 0.50000;    /* plane separation factor for 'fcc 200' */
      Cgeo[0][4] = 0.35355;    /* plane separation factor for 'fcc 220' */
      Cgeo[1][1] = 0.70711;    /* plane separation factor for 'bcc 110' */
      Cgeo[1][3] = 0.50000;    /* plane separation factor for 'bcc 200' */
      Cgeo[2][2] = 0.57735;    /* plane separation factor for 'diamond 111' */
      Cgeo[2][4] = 0.35355;    /* plane separation factor for 'diamond 220' */
      Cgeo[3][3] = 0.50000;    /* plane separation factor for 'cph 200' */
      Cgeo[4][0] = 1.00000;    /* plane separation factor for 'cub 100' */
      Cgeo[4][1] = 0.70711;    /* plane separation factor for 'cub 110' */
      Cgeo[4][2] = 0.57735;    /* plane separation factor for 'cub 111' */
      Cgeo[4][3] = 0.50000;    /* plane separation factor for 'cub 200' */
      Cgeo[4][4] = 0.35355;    /* plane separation factor for 'cub 220' */
   
      for (i=0; i<5; i++) {
         for (j=0; j<5; j++) {
             if (Cgeo[i][j] > 0.0) atoms_eff[i][j] = 1.0;
         }
      }
      atoms_eff[2][2] = 0.5;	/* special case for diamond */ 

      atoms_cell[0] = 4;	/* atoms per cell for 'fcc' */
      atoms_cell[1] = 2;	/* atoms per cell for 'bcc' */
      atoms_cell[2] = 8;	/* atoms per cell for 'diamond' */
      atoms_cell[3] = 4;	/* atoms per cell for 'cph' */
      atoms_cell[4] = 1;	/* atoms per cell for 'cub' */

      in = NULL;
      in = fopen("bragg_t.dat", "rt");
      if (in == NULL) {
         printf("Could not open 'bragg_t.dat'\n");
         return(-1);
      }

      fgets(iline, 200, in);	/* get header lines */
      fgets(iline, 200, in);
      fgets(iline, 200, in);
   
      Cstruc[ZZ] = -1;
next:
      fgets(iline, 200, in);
      if (feof(in)) goto completed;

      i = j = 0;
      sscanf(iline, "%d", &Z);

      while ((i<200) && (j<1)) if (iline[i++] == 9) j++;
      strncpy(elem[Z], &iline[i], 2);
      elem[Z][2] = 0;
      sscanf(&iline[5], "%lf %lf %lf %lf", &A[Z], &dens[Z], &a_dens[Z], &DebyeT[Z]);

      while ((i<200) && (j<6)) if (iline[i++] == 9) j++;
      strncpy(Struc[Z], &iline[i], 3);
      Struc[Z][3] = 0;

      sscanf(&iline[i+4], "%lf %lf", &dist[Z], &melt[Z]);
      while ((i<200) && (j<8)) if (iline[i++] == 9) j++;
      if (j<7) melt[Z] = 0.00;
   
      Cstruc[Z] = -1;
      if (strstr(Struc[Z], "fcc") != NULL) Cstruc[Z] = 0;
      if (strstr(Struc[Z], "bcc") != NULL) Cstruc[Z] = 1;
      if (strstr(Struc[Z], "dia") != NULL) Cstruc[Z] = 2;
      if (strstr(Struc[Z], "cph") != NULL) Cstruc[Z] = 3;
      if (strstr(Struc[Z], "cub") != NULL) Cstruc[Z] = 4;
      goto next;
   }

/* structure characterizing a chemical element   
struct material {
   char  ELEM_name[10];		 Element name 
   char  ELEM_symb[3];		 Element symbol 
   char  Struc_txt[4];		 Crystal structure text 
   int   Cstruc;		 Crystal structure indicator (0-5) 
      				 0: "fcc", 1: "bcc", 2: "dia", 3: "cph", 4: "cub" 
   float A;			 mean atomic mass 
   float DebyeT;		 Debye temperature (K) 
   float atom_dens;		 atomic density (1/cm3)
   float dens;			 density (g/cm3) 
   float dist;			 unit cell atomic distance (cm)
   float melt;			 melting point 
   float ThermalFact;		 thermal factor 
   float formfact_atom[4];	 atomic form factor (n_order [1-3])
   float FormFact[5][4];	 Unit cell form factor for reflection type [0-4] and n_order [1-3] 
};

int initREFL(int ZZ, material *MATER, float Cgeo[5][5], float atoms_eff[5][5], int atoms_cell[],
             reflect REFL[], zones ZONES[], int n_zones, int nOptions, float fwhmt, float minthick, int OptZ[])
*/
completed: ;
   if (in != NULL) fclose(in);
   in = NULL;

   if (Cstruc[ZZ] < 0) {
      printf("Z: %2d. Unknown crystal structure: %3s\n", ZZ, Struc[ZZ]);
      return(-1);
   } 

   MATER->Z		= ZZ;
   MATER->A		= A[ZZ];
   MATER->Cstruc	= Cstruc[ZZ];
   MATER->DebyeT	= DebyeT[ZZ];
   MATER->atom_dens	= a_dens[ZZ] * 1.0e22;		/* change to atoms/cm3 */
   MATER->dens		= dens[ZZ];
   MATER->dist		= dist[ZZ] * 1.0e-8;		/* change from Angströms to cm */
   MATER->melt		= melt[ZZ];
   strncpy(MATER->ELEM_symb, elem[ZZ], 2);
   strncpy(MATER->Struc_txt, Struc[ZZ], 3);
      
   B = 2.29e-12 * T / (A[ZZ] * DebyeT[ZZ] * DebyeT[ZZ]); /* thermal factor (note units: cm for all dimensions!) */
   MATER->ThermalFact = B;

   for (CM=0; CM<5; CM++) {
      MATER->FormFact[CM][0] = MATER->FormFact[CM][1] = MATER->FormFact[CM][2] = MATER->FormFact[CM][3] = 0.0;

      if (Cgeo[Cstruc[ZZ]][CM] < 0.0) continue;
          
      d = MATER->dist * Cgeo[Cstruc[ZZ]][CM];  	/* separation between crystal planes (cm) */
      for (n_order=1; n_order<4; n_order++) {
         SL = (double)(n_order) / (2.0e8 * d);
         f_atom = CromerWaber(ZZ, SL); 		/* atomic form factor */
         F_cell = f_atom * (double)(atoms_cell[Cstruc[ZZ]]) * atoms_eff[Cstruc[ZZ]][CM]; /* unit cell structure factor */
         MATER->FormFact[CM][n_order] = F_cell;
         MATER->formfact_atom[n_order] = f_atom;
         /* special case for higher orders in 'diamond structure 111' */
         if ((Cstruc[ZZ] == 2) && (CM == 2) && ((n_order%2) == 0)) MATER->FormFact[CM][n_order] = 0.0;
      }

      REFL[nOptions].dist_geo = d;
      REFL[nOptions].Z = ZZ;  
      REFL[nOptions].fwhmt = fwhmt; 
      REFL[nOptions].CM = CM;
      REFL[nOptions].minthick = minthick;
      REFL[nOptions].atoms_eff = atoms_eff[Cstruc[ZZ]][CM];
      sprintf(REFL[nOptions].Miller_txt, "%s", Miller[CM]);

      for (k=0; k<numLayers; k++) {
         for (n_z=0; n_z<n_zones; n_z++) {
            for (azimBin=0; azimBin<3; azimBin++) {
               if (ZONES[n_z].Z_lay[k][azimBin] == -1) {ZONES[n_z].Opt_lay[k][azimBin] = -1; continue;}
               if ((ZONES[n_z].Z_lay[k][azimBin] != ZZ) || (ZONES[n_z].Millay[k][azimBin] != atoi(Miller[CM]))) continue;
               ZONES[n_z].Opt_lay[k][azimBin] = nOptions;
            }
            for (azimBin=0; azimBin<3; azimBin++) {
               if (ZONES[n_z].secZ_lay[k][azimBin] == -1) {ZONES[n_z].secOpt_lay[k][azimBin] = -1; continue;}
               if ((ZONES[n_z].secZ_lay[k][azimBin] != ZZ) || (ZONES[n_z].secMillay[k][azimBin] != atoi(Miller[CM]))) continue;
               ZONES[n_z].secOpt_lay[k][azimBin] = nOptions;
            }
         }
      }
      OptZ[nOptions] = ZZ;

      nOptions++;
   }
   return(nOptions);
}
/*
** structure characterizing a zone with particular reflections types **
typedef struct {
   float fR_outer;	** Inner radius of zone (as fraction of the lens width) **
   float fR_inner;	** Outer radius of zone (as fraction of the lens width) **
   float R_outer;	** Inner radius of zone **	** DERIVED value **
   float R_inner;	** Outer radius of zone **	** DERIVED value **
   int   Z_lay[2][3];
   int   Millay[2][3];
   float ftlay[2][3];	** Thickness factor relative to 'optimal' thickness **
   float fracA[2][3];		** start fraction of potential redefined facets **
   float fracB[2][3];		** end fraction of potential redefined facets **
   int   Opt_lay[2][3];	** Option value corresponding to Z and Miller-index ** ** DERIVED value **
   int   secZ_lay[2][3];
   int   secMillay[2][3];
   float secftlay[2][3];	** Thickness factor relative to 'optimal' thickness **
   float secfracA[2][3];	** start fraction of potential redefined facets **
   float secfracB[2][3];	** end fraction of potential redefined facets **
   int   secOpt_lay[2][3];	** Option value corresponding to Z and Miller-index ** ** DERIVED value **
} zones;

   sscanf(iline, "%f %f  %d %d %f  %d %d %f  %d %d %f", &(ZONES[nzon].fR_outer), &(ZONES[nzon].fR_inner),
      &ZONES[nzon].Z_lay[0][0], &ZONES[nzon].Millay[0][0], &ZONES[nzon].ftlay[0][0],
      &ZONES[nzon].Z_lay[0][1], &ZONES[nzon].Millay[0][1], &ZONES[nzon].ftlay[0][1],
      &ZONES[nzon].Z_lay[0][2], &ZONES[nzon].Millay[0][2], &ZONES[nzon].ftlay[0][2]);
   sscanf(iline, "%d %d %f %f %f   %d %d %f %f %f   %d %d %f %f %f",
      &ZONES[nzon].secZ_lay[0][0], &ZONES[nzon].secMillay[0][0], &ZONES[nzon].secftlay[0][0],
      &ZONES[nzon].secfracA[0][0], &ZONES[nzon].secfracB[0][0], 
      &ZONES[nzon].secZ_lay[0][1], &ZONES[nzon].secMillay[0][1], &ZONES[nzon].secftlay[0][1],
      &ZONES[nzon].secfracA[0][1], &ZONES[nzon].secfracB[0][1], 
      &ZONES[nzon].secZ_lay[0][2], &ZONES[nzon].secMillay[0][2], &ZONES[nzon].secftlay[0][2],
      &ZONES[nzon].secfracA[0][2], &ZONES[nzon].secfracB[0][2]);
*/