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

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

/* Calculation of Bragg scattering efficiency from mosaic crystals
   of finite thickness. Nomenclature follows "A study of focusing...", nl, 1992
   Niels Lund June 2016

   This routine calculates the thickness of each facet according to the supplied reflection type and t_factor
*/

int  reflectT2(reflect *REFL, material *MATER, facet *FACET, float t_factor)
{
        char iline[200];

	int i, j, k, err, Z, iOrder = 1;
        static int first = 1;
        
        char unit = 'b', errmsg[100];
   	int print_flag = 1;
   	double energy[10], xsec[12], fl_yield[6];

	double phi, thetaB, mu0, t, sigmaE, alpha, nypeak, R, x, y, F_cell, RR, V, RRR; 
	double lambda, d, a_vol, kappa, B, L, sigmat, f, sqrt2pi, thermal_f, E;
        float  f_atom, dist, dist_geo, fwhmt, T, dens, dev_sig=0.0;
	double tmax, nymax, re, re2, alfa, alfa1, fwhmrad, sigmarad, tt, rr;
        float SL, f_old, ny_minthick, ny_redu;

	sqrt2pi = sqrt(TWOPI);
	re = 2.818e-13;	/* classical electron radius (cm) */
	re2 = re * re;
        
/*        dist = MATER->dist;
        dist_geo = dist * Cgeo[MATER->Cstruc][REFL->CM]; */
        dist_geo = REFL->dist_geo;
        B      = MATER->ThermalFact;
        f_atom = MATER->formfact_atom[iOrder];
        F_cell = MATER->FormFact[REFL->CM][iOrder];
	a_vol = 1.0 / MATER->atom_dens;			/* atomic volume (cm3) */
        Z = MATER->Z;
        T = 300.0;
        E = FACET->E[FACET->nOption];			/* nominal Bragg energy for this facet */

/*        if ((((first++)%N_CryPos) < 10) && (E > 0.0))
           printf("reflectT2: Z %2d %s dist_geo %6.4f E %6.1f (keV)\n", Z, MATER->ELEM_symb, dist_geo * 1.0e+8, E);
*/
        print_flag = 1;
        err = mucal(MATER->ELEM_symb, Z, E, unit, print_flag, energy, xsec, fl_yield, errmsg);

        V = MATER->dist * MATER->dist * MATER->dist;	/* unit cell volume */
        kappa = xsec[3]; 				/* kappa is the total atomic cross section in barns/atom */
        dens  = xsec[7];
        mu0 = kappa / xsec[4] * dens;			/* convert to 1/cm */
        fwhmt = REFL->fwhmt;
	sigmat = fwhmt / 2.36;				/* convert from fwhm to sigma */
	lambda = 12.39 / (E * (double)(iOrder));	/* INCLUDED iOrder 1/4 2016 NL. lambda in Angstroem, E in KeV */
        lambda  *= 1.0e-8;				/* change from Angströms to cm */
        fwhmrad = fwhmt * 3.14159 / (180.0 * 60.0);
        sigmarad =fwhmrad / 2.36;

	thetaB = (double)(iOrder) * lambda / (2.0 * dist_geo);
	sigmaE = E * sigmarad / thetaB;
        FACET->dE[FACET->nOption] = sigmaE;
	x = (double)(iOrder) /(2.0 * dist_geo);
        y = MATER->DebyeT * MATER->DebyeT / (T * T);
        B *= (1.0 + y / 36.0 - (y * y) / 3600.0);
        thermal_f = exp(-B * x * x);

	R = 2.0 * re2 * lambda * lambda * f_atom * f_atom * thermal_f  * dist_geo / ((double)(iOrder) * a_vol * a_vol);
/*	RR = 2.0 * re2 * lambda * lambda * F_cell * F_cell * thermal_f * dist_geo / ((double)(iOrder) * V * V);
        RRR = 2.0 * re2 * lambda * lambda * lambda * F_cell * F_cell * thermal_f / (V * V * sin(2.0 * thetaB));
*/
	alfa = R / (mu0 * sigmarad * sqrt2pi);	/* HERE MUST BE sigmarad !!! */
	alfa1 = 1.0 + alfa;
        
        if (alfa > 0.0) {
           tmax = log(alfa1) / (mu0 * alfa);
           nymax = 0.5 * alfa * exp(log(alfa1) * (-alfa1 / alfa));
        } else {
           tmax = 0.00001;
           nymax = 0.0;
           thetaB = 0.00000001;
        }

        FACET->thick = tmax * t_factor;
	if (FACET->thick < REFL->minthick) FACET->thick = REFL->minthick; 

        return(1);
stop:
	return(-1);
}