#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 
   sub-peak thicknesses of the crystals
*/

int reflect3A(float E, reflect *REFL, facet *FACET, material *MATER, int iOrder, float *ref_E)
{
        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;

	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_E = 0.5 * exp(-mu0 * tt) * ( 1.0 - exp(-R * tt / (sigmarad * sqrt2pi)));

/*        if ((iOrder == 1) && (((int)(E) % 50) == 0) && (fabs(tmax-tt) < 0.01)) 
           printf("reflect3A: E: %6.1f, inOpt: %2d, topt: %5.2f, tt: %5.2f (mm), refl: %5.2f\n", E, FACET->nOption, tmax*10.0, tt*10.0, *ref_E);
*/
        return(1);
stop:
	return(-1);
}