/* Calculation of Bragg scattering efficiency from mosaic crystals
   of finite thickness. Nomenclature follows "Afstudy of focusing...", nl, 1992
   Niels Lund December 2001
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

double pi = 3.141593, gtorad, mtorad;

void main(int argc, char *argv[])
{
	int i, j, k, n, Z;
	double phi, thetaB, mu, t, sigmaE, E, alpha, nypeak, R, x;
	double lambda, d, a, kappa, B, L, sigmat, fwhmt, f, sqrt2pi;
	double A[100], dens[100], a_dens[100], Debye[100], Cgeo[4][2];
        double dist[100], melt[100], klist[4];
	double tmax, nymax, re, re2, alfa, alfa1;
        int Cstruc, Cmiller;
        char iline[200], elem[100][6], struc[100][10], Miller[10];

        FILE *in, *out;

	gtorad = pi/180.0;
	mtorad = gtorad / 60.0;
	sqrt2pi = sqrt(2.0 * pi);
	re = 2.818e-11;	/* classical electron radius (m) */
	re2 = re * re;

        Cgeo[0][0] = 0.5774;    /* plane separation factor for 'fcc 111' */
        Cgeo[0][1] = 0.5000;    /* plane separation factor for 'fcc 200' */
        Cgeo[1][0] = 0.7071;    /* plane separation factor for 'bcc 110' */
        Cgeo[1][1] = 0.5000;    /* plane separation factor for 'bcc 200' */
        Cgeo[2][0] = 0.5769;    /* plane separation factor for 'diamond 111' */
        Cgeo[2][1] = 0.3534;    /* plane separation factor for 'diamond 200' */
        Cgeo[3][0] = 0.5000;    /* plane separation factor for 'cph 1000' */
        Cgeo[3][1] = 0.5000;    /* plane separation factor for 'cph 1000' */
        Cgeo[4][0] = 1.0000;    /* plane separation factor for 'cub 100' */
        Cgeo[4][1] = 0.7071;    /* plane separation factor for 'fcc 110' */

        klist[0] = 300.0;
        klist[1] = 100.0;
        klist[2] = 100.0;
        klist[3] = 600.0;

        out = fopen("bragg_out.txt", "wt");

        in = fopen("bragg_t.txt", "rt");
        if (in == NULL) {
                fprintf(out, "Could not open 'bragg_t.txt'\n");
                goto stop;
        }

        fgets(iline, 200, in);
        printf("%s", iline);
next:
        fgets(iline, 200, in);
        if (feof(in)) goto got_it;

        printf("%s", iline);
        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], &Debye[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;
        
        printf("%2d      %2s      %4.1lf    %5.2lf    %5.2lf  %6.1lf   %3s  %6.2lf  %6.1lf\n",
                  Z, elem[Z], A[Z], dens[Z], a_dens[Z], Debye[Z], struc[Z], dist[Z], melt[Z]);

        Cstruc = -1;
        if (strstr(struc[Z], "fcc") != NULL) Cstruc = 0;
        if (strstr(struc[Z], "bcc") != NULL) Cstruc = 1;
        if (strstr(struc[Z], "dia") != NULL) Cstruc = 2;
        if (strstr(struc[Z], "cph") != NULL) Cstruc = 3;
        if (strstr(struc[Z], "cub") != NULL) Cstruc = 4;
        if (Cstruc < 0) {
                fprintf(out, "Z: %2d. Unknown crystal structure: %3s\n", Z, struc[Z]);
                printf("Z: %2d. Unknown crystal structure: %3s\n", Z, struc[Z]);
                goto stop;
        }
        goto next;
got_it:
	E = 511.0;      /* default: annihilation photons */
	n = 1;		/* default: first order diffraction */
	Z = 29;		/* default: copper crystals */
	L = 15.0;	/* default: 15 m focal length */
	fwhmt = 1.0;	/* default: 1 arcminute mosaic width */
        Cmiller = 0;    /* default: Miller index corresponding to minimum Bragg angle */

	if (argc > 1) fwhmt = atof(argv[1]);	/* argv in arcminutes */
	if (argc > 2) E = atof(argv[2]);	/* argv in keV */
	if (argc > 3) n = atoi(argv[3]);        /* diffraction order */
	if (argc > 4) Z = atoi(argv[4]);
        if (argc > 5) {
                strcpy(Miller, argv[5]);
                Cmiller = -1;
                if (strstr(Miller, "111") != NULL) Cmiller = 0;
                if (strstr(Miller, "200") != NULL) Cmiller = 1;
                if (strstr(Miller, "220") != NULL) Cmiller = 1;
                if (strstr(Miller, "110") != NULL) Cmiller = 0;
                if (strstr(Miller, "100") != NULL) Cmiller = 0;
                if (strstr(Miller, "101") != NULL) Cmiller = 0;
                if (Cmiller < 0) {
                        fprintf(out, "Unknown Miller index: %s\n", Miller);
                        goto stop;
                }
        }
	d = dist[Z] * Cgeo[Cstruc][Cmiller];    /* separation between crystal planes */
	B = Debye[Z] / 300.0;                   /* thermal factor */
	a = 1.0 / a_dens[Z];                    /* atomic volume */
	f = (double)(Z);                        /* atomic form factor - in first approximation = Z */

                /* The total atomic cross section is expressed as a
                   sum of the photoelectric and pair cross sections */
	if (E < klist[3])
		kappa = klist[0] + (klist[1] - klist[0]) * (E - 100.0) / (klist[3] - 100.0);
	else 	kappa = klist[1] * exp(log(E - klist[3]) * 1.5);  /* kappa is the total atomic cross section */

	sigmat = fwhmt / 2.36;	/* convert from fwhm to sigma */
	lambda = 12.39 / E;	/* lambda in Angstroem */
	thetaB = n * lambda / 2.0 * d;
	sigmaE = E * sigmat / thetaB;
	x = (double)(n) /(2.0 * d);

	mu = kappa / a;
	R = 2.0 * re2 * lambda * lambda * f * exp(B * x * x) * d / n;
	alfa = R / (mu * sigmaE * sqrt2pi);
	alfa1 = 1.0 + alfa;
	tmax = log(alfa1) / (mu * alfa);
	nymax = 0.5 * alfa * exp(log(alfa1) * (-alfa1 / alfa));
        fprintf(out, "Calculated peak reflectivity of %s (%s) for the %s reflection at %5.1lf keV: %4.1lf\%\n",
                elem[Z], struc[Z], Miller, E, nymax);
        fprintf(out, "Optimal crystal thickness: %4.1lf mm for mosaic width (sigma): %4.2lf arcmin\n",
                tmax, sigmat);
stop:
        fclose(out);
        fclose(in);
	return;
}