#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 diffraction efficiency for for specified thicknesses of the crystals
*/

int reflect3B(float E, reflect *REFL, facet *FACET, material *MATER, int iOrder, float ref_E3[])
{
        char iline[200];

	int i, j, k, err, Z, len;
        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;
        float  f_atom, dist, dist_geo, fwhmt, T, dens;
	double tmax, nymax, re, re2, alfa, alfa1, fwhmrad, sigmarad, tt, rr;
        float SL, f_old, ny_minthick, ny_redu;

        static int numcalls = 0;

	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;

        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;
	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;
        }

        tt = FACET->thick;

        ref_E3[0] = 0.5 * exp(-mu0 * tt) * ( 1.0 - exp(-R * tt / (sigmarad * sqrt2pi)));
        ref_E3[1] = 0.5 * exp(-mu0 * tt * 2.0) * ( 1.0 - exp(-R * tt / (sigmarad * sqrt2pi)));
        ref_E3[2] = 0.5 * exp(-mu0 * tt * 3.0) * ( 1.0 - exp(-R * tt / (sigmarad * sqrt2pi)));
        
/*        if ((numcalls++ % 5000) == 0) */
/*        if ((FACET->nOption == 5) && (iOrder == 1))
           printf("reflect3B: E: %6.1f, inOpt: %2d, iOrder: %1d, topt: %5.2f, tt: %5.2f (mm), refl: %5.2f %5.2f %5.2f\n", 
              E, FACET->nOption, iOrder, tmax*10.0, tt*10.0, ref_E3[0], ref_E3[1], ref_E3[2]);
*/
        return(1);
stop:
	return(-1);
}