#include #include #include #include #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 This version calculates the scattering efficiency for a three layer lens. The routine must be called for each facet position in each layer (some of which may be empty). For each filled facet the routine checks for overlapping facets in the two neighboring layers. The area of each facet is split into fractions with different overlap situations. Normally we will have 3 different fractions (we will treat the three layers equally and assume that all three layers will have similar absorbtion properties): 1) No overlap 2) Overlap with one layer 3 Overlap with two layers Niels Lund September 2005 and February 2019 */ float facet_Bcontrib( float E, int iOpt, int iOrder, float panel_new[], facet *FACET, material *MATER, reflect *REFL, int LEN[], int L_order[numLayers][N_CryPos], float *refl_E, int index, float NABSTAT[3][30], FILE* gnufil, int Lens_layers) { int i, j, n, nabo_index, status=0, N_Z, Z, err, inabo, ilayer; int t1, t2, t3, t4, t5, t6, t7, t8, t9, layer=0, thislayer=0; static int first=1, prcnt=0, index_old=0; float sigma, ref_specif, abs_specif, dev_sig, fx, fy, fr, dx, dy, fdx, fdy; float reflect, absorb_fact[6], theta, E_B, theta_B, fwhmt, det_ovlp, rr, RR, Dx, Dy; float eff_area, af, B_contr, dE_B, d_theta, dis, overlap_frac, overlap_frac_h, A_final; float kappa, dens, mu0, A_int, th, Pth[6], Pmu0[6]; float AAA, alpha, beta, AH0, AH1, Frac0, Frac1, Frac2, Frac3, refl_E3[numLayers]; float xA, xB, yA, yB, xAmean, yAmean, xBmean, yBmean, distAB, GaussFactor; static double sqrt2PI; static float arcmin2rad, foc_len, facet_radius, facet_hole_radius, A_basic, A_hole, detector_radius, fstep; static float A_crys, RRlimA, RRlimB, facet_sep, sq3half, old_x=0.0, old_y=0.0; static float rel_facet_sep, theta_ov2, rel_facet_sep2, gamma, theta_ov3, triang3, afskaer3; static float overlap0, overlap1, overlap2, overlap3, rel_hexagon, olap2, olap3, theta_ov2a, olap1; static float useFractions[4][4][3], O2, O3, O2b; char unit = 'b', errmsg[100]; int print_flag = 1, numnabo[2]; double energy[10], xsec[12], fl_yield[6]; static int x0=0, x1=0, x2=0, x3=0, x4=0, nabcnt=0; thislayer = index / N_CryPos; if ((thislayer < 0) || (thislayer > 2)) { printf("thislayer: %6d, N_CryPos: %6d, index: %8d\n", thislayer, N_CryPos, index); exit(0); } B_contr = 0.0; if (FACET[index].Z <= 0) return(B_contr); if (FACET[index].thick <= 0.0) return(B_contr); NABSTAT[thislayer][23+thislayer] += 1.0; if (first) { /* The 'first'-section ends in line 241 */ PI = 3.14159; first = 0; foc_len = panel_new[0]; /* cm */ detector_radius = panel_new[13]; /* cm */ facet_radius = panel_new[2] * 0.5; /* cm */ facet_hole_radius = panel_new[12]; /* cm */ facet_sep = panel_new[14]; /* cm */ sq3half = 0.5 * sqrt(3.0); /* cm */ A_basic = PI * facet_radius*facet_radius; /* Area within crystal perimeter cm2 */ A_hole = PI * facet_hole_radius*facet_hole_radius; /* Area within crystal hole cm2 */ A_crys = A_basic - A_hole; GTORAD = PI / 180.0; arcmin2rad = PI / (180.0 * 60.0); sqrt2PI = sqrt(2.0 * PI); fstep = 0.1; /* integration step 1 mm (in both dimensions) */ RRlimA = facet_sep * sqrt(panel_new[5]*panel_new[5] + panel_new[6]*panel_new[6]) * 0.99; RRlimB = RRlimA * 1.02; /* Calculate overlapping area fractions for neighboring crystals in different layers */ /* Overlap between two neighboring crystals, center to center distance: facet_sep, crystal radius: facet_radius */ /* theta_ov is the half opening angle of the overlap sector of the two facets */ rel_facet_sep = facet_sep / facet_radius; rel_facet_sep2 = rel_facet_sep / (2.0 * cos(30.0 * GTORAD)); /* printf("Facets relative dimensions: facet radius: %5.2f, hole radius: %5.2f, facet separation: %5.3f, facet separation2: %5.3f\n", facet_radius/facet_radius, facet_hole_radius/facet_radius, rel_facet_sep, rel_facet_sep2); */ if (rel_facet_sep <= 2.0) { printf("The facets are too close, facets overlap within the same layer\n"); exit(0); } if (rel_facet_sep2 <= (facet_radius + facet_hole_radius) / facet_radius) { printf("The facets are too close, not sufficient room for facet supports. f_rad: %5.2f, h_rad: %5.2f, f_sep: %5.2f, rf_sep2: %6.3f\n", facet_radius, facet_hole_radius, facet_sep, rel_facet_sep2); exit(0); } overlap1 = PI * (1.0 - (facet_hole_radius / facet_radius) * (facet_hole_radius / facet_radius)); olap1 = overlap1 / 12.0; /* prepare for factor applied for normal configuration */ AAA = rel_facet_sep2; if (AAA >= 2.0) { /* no facets overlap */ alpha = 0.0; beta = 0.0; triang3 = 0.0; afskaer3 = 0.0; overlap2 = 0.0; overlap3 = 0.0; } else { /* facet from different layers overlap */ alpha = acos(AAA / 2.0); overlap2 = (alpha - sin(alpha)*cos(alpha)) * 2.0; olap1 = (overlap1 - overlap2); beta = alpha - 30.0 * GTORAD; if (beta <= 0.0) { /* but there are no triple overlaps */ triang3 = 0.0; afskaer3 = 0.0; overlap3 = 0.0; } else { /* triple overlaps occur */ triang3 = sin(beta)*sin(beta) * tan(60.0 *GTORAD); afskaer3 = beta - sin(beta)*cos(beta); overlap3 = triang3 + 3.0 * afskaer3; gamma = asin((AAA / 2.0 - sin(beta)) * tan(30.0 * GTORAD)); overlap1 = 0.5 * AAA * (AAA / 2.0 - sin(beta)) * tan(30.0 * GTORAD) - 0.5 * gamma; olap1 = 12.0 * overlap1; } } rel_hexagon = 6.0 * (rel_facet_sep / 2.0) * (rel_facet_sep / 2.0) * tan(30.0 * GTORAD) / PI; /* printf("Overlap calculations: relative facet separation: %5.3f, rel_sep2: %5.3f\n", rel_facet_sep, rel_facet_sep2); printf("alpha: %5.3f overlp2 %5.3f beta %5.3f tri3 %5.3f afsk3 %5.3f, overlp3 %5.3f gamma %5.3f overlp1: %5.3f, rel.hexag %5.3f\n\n", alpha, overlap2, beta, triang3, afskaer3, overlap3, gamma, olap1, rel_hexagon); */ /* calculation of total overlap areas within a unit hexagon */ overlap0 = 3.0*A_hole; /* one central hole plus a third of six corner holes */ olap3 = 6.0 * overlap3 * A_basic /PI; /* six overlap3 areas within the central facet in the hexagon */ olap2 = 18.0 * (overlap2 / 2.0 - overlap3) * A_basic /PI; /* the real overlap2 consists of (half the primary overlap2 - overlap3) */ /* in each hexagon there are 18 such contributions */ overlap1 = olap1 * A_basic / PI; /* AH0 = rel_hexagon * A_basic * 0.01; AH1 = (overlap1 +olap2 + olap3); printf("Areas (cm2): cryst: %6.2f, central hole: %6.2f, hexagon: %6.2f, overlap0: %6.2f, overlap1: %6.2f, overlap2: %6.2f, overlap3: %6.2f\n", A_crys, A_hole, rel_hexagon*A_basic, overlap0, overlap1, olap2, olap3); printf("Fractions (%): hexagon: %6.1f, cryst: %6.1f, central hole: %6.1f, overlap0: %6.1f, overlap1: %6.1f, overlap2: %6.1f, overlap3: %6.1f\n", rel_hexagon*A_basic/AH0, A_crys/AH0, A_hole/AH0, overlap0/AH0, overlap1/AH0, olap2/AH0, olap3/AH0); printf("Fractions (%): crystal: %6.1f, overlap1: %6.1f, overlap2: %6.1f, overlap3: %6.1f\n", 100.0, 100.0*overlap1/AH1, 100.0*olap2/AH1, 100.0*olap3/AH1); */ alpha = acos(rel_facet_sep2 / 2.0); overlap2 = (alpha - sin(alpha)*cos(alpha)) * 2.0; O2 = overlap2 / (PI * 0.96); /* Area of overlap between two unit-circles with center-to-center separation of 'rel_facet_sep2' */ beta = alpha - 30.0 * GTORAD; triang3 = sin(beta)*sin(beta) * tan(60.0 *GTORAD); afskaer3 = beta - sin(beta)*cos(beta); overlap3 = triang3 + 3.0 * afskaer3; O3 = overlap3 / (PI * 0.96); /* Area of overlap between three unit-circles with center-to-center separation of 'rel_facet_sep2' */ O2b = (overlap2 - 2.0 * overlap3) / 6.0; /* Area of double overlap between two adjacent triple overlap regions */ /* printf("Overlap2: %6.4f (should be 0.2824 for nominal lens configuration (40 mm facets with 8 mm holes on 42.5 mm pitch) \n", O2); printf("Overlap3: %6.4f (should be 0.1189 for nominal lens configuration (40 mm facets with 8 mm holes on 42.5 mm pitch) \n", O3); printf("Overlap2b: %6.4f (should be 0.0223 for nominal lens configuration (40 mm facets with 8 mm holes on 42.5 mm pitch) \n", O2b); */ useFractions[3][3][0] = overlap1 / AH1; useFractions[3][3][1] = olap2 / AH1; useFractions[3][3][2] = olap3 /AH1; useFractions[3][3][0] = 1.0 - O2 * 3.0 - O2b * 3.0; useFractions[3][3][1] = O2 * 3.0 - O3 * 6.0 + O2b * 3.0; useFractions[3][3][2] = O3 * 6.0; useFractions[3][2][0] = 1.0 - O2 * 3.0 - O2b * 2.0; useFractions[3][2][1] = O2 * 3.0 - O3 * 4.0 + O2b * 2.0; useFractions[3][2][2] = O3 * 4.0; useFractions[3][1][0] = 1.0 - O2 * 3.0 - O2b; useFractions[3][1][1] = O2 * 3.0 - O3 * 2.0 + O2b; useFractions[3][1][2] = O3 * 2.0; useFractions[3][0][0] = 1.0 - O2 * 3.0; useFractions[3][0][1] = O2 * 3.0; useFractions[3][0][2] = 0.0; useFractions[2][3][0] = 1.0 - O2 * 3.0 - O2b * 2.0; useFractions[2][3][1] = O2 * 3.0 - O3 * 4.0 + O2b * 2.0; useFractions[2][3][2] = O3 * 4.0; useFractions[2][2][0] = -1.0; useFractions[2][2][1] = 0.0; useFractions[2][2][2] = 0.0; useFractions[2][1][0] = -2.0; useFractions[2][1][1] = 0.0; useFractions[2][1][2] = 0.0; useFractions[2][0][0] = 1.0 - O2 * 2.0; useFractions[2][0][1] = O2 * 2.0; useFractions[2][0][2] = 0.0; useFractions[1][3][0] = 1.0 - O2 * 3.0 - O2b; useFractions[1][3][1] = O2 * 3.0 - O3 * 2.0 + O2b; useFractions[1][3][2] = O3 * 2.0; useFractions[1][2][0] = -3.0; useFractions[1][2][1] = 0.0; useFractions[1][2][2] = 0.0; useFractions[1][1][0] = -4.0; useFractions[1][1][1] = 0.0; useFractions[1][1][2] = 0.0; useFractions[1][0][0] = 1.0 - O2; useFractions[1][0][1] = O2; useFractions[1][0][2] = 0.0; useFractions[0][3][0] = 1.0 - O2 * 3.0; useFractions[0][3][1] = O2 * 3.0; useFractions[0][3][2] = 0.0; useFractions[0][2][0] = 1.0 - O2 * 2.0; useFractions[0][2][1] = O2 * 2.0; useFractions[0][2][2] = 0.0; useFractions[0][1][0] = 1.0 - O2; useFractions[0][1][1] = O2; useFractions[0][1][2] = 0.0; useFractions[0][0][0] = 1.0; useFractions[0][1][1] = 0.0; useFractions[0][1][2] = 0.0; /* for (i=0; i<4; i++) for (j=0; j<4; j++) printf("useFractions: %1d %1d %6.4f %6.4f %6.4f sum: %6.4f\n", i, j, useFractions[i][j][0], useFractions[i][j][1], useFractions[i][j][2], useFractions[i][j][0]+useFractions[i][j][1]+useFractions[i][j][2]); printf("\nNeighbor situations with two configurations, here with axial symmetry:\n"); useFractions[2][2][0] = 1.0 - 4.0 * (O2 - O3) - 2.0 * O3; useFractions[2][2][1] = 4.0 * (O2 - O3); useFractions[2][2][2] = 2.0 * O3; printf("useFractions: %1d %1d %6.4f %6.4f %6.4f sum: %6.4f\n", 2, 2, useFractions[2][2][0], useFractions[2][2][1], useFractions[2][2][2], useFractions[2][2][0]+useFractions[2][2][1]+useFractions[2][2][2]); useFractions[1][1][0] = 1.0 - 2.0 * O2; useFractions[1][1][1] = 2.0 * O2; useFractions[1][1][2] = 0.0; printf("useFractions: %1d %1d %6.4f %6.4f %6.4f sum: %6.4f\n", 1, 1, useFractions[1][1][0], useFractions[1][1][1], useFractions[1][1][2], useFractions[1][1][0]+useFractions[1][1][1]+useFractions[1][1][2]); useFractions[2][1][0] = useFractions[1][2][0] = 1.0 - 2.0 * (O2 - O3) - 2.0 * O3; useFractions[2][1][1] = useFractions[1][2][1] = 2.0 * (O2 - O3); useFractions[2][1][2] = useFractions[1][2][2] = 2.0 * O3; printf("useFractions: %1d %1d %6.4f %6.4f %6.4f sum: %6.4f\n", 2, 1, useFractions[2][1][0], useFractions[2][1][1], useFractions[2][1][2], useFractions[2][1][0]+useFractions[2][1][1]+useFractions[2][1][2]); printf("\nNeighbor situations with two configurations, here without axial symmetry:\n"); useFractions[2][2][0] = 1.0-2.0*(O2-O3) - 2.0*O2b - 3.0*O3; useFractions[2][2][1] = 2.0*(O2-O3)+2.0*O2b; useFractions[2][2][2] = 3.0 * O3; printf("useFractions: %1d %1d %6.4f %6.4f %6.4f sum: %6.4f\n", 2, 2, useFractions[2][2][0], useFractions[2][2][1], useFractions[2][2][2], useFractions[2][2][0]+useFractions[2][2][1]+useFractions[2][2][2]); useFractions[1][1][0] = 1.0 - 2.0 * (O2 -O3) - O3; useFractions[1][1][1] = 2.0 * (O2 - O3); useFractions[1][1][2] = O3; printf("useFractions: %1d %1d %6.4f %6.4f %6.4f sum: %6.4f\n", 1, 1, useFractions[1][1][0], useFractions[1][1][1], useFractions[1][1][2], useFractions[1][1][0]+useFractions[1][1][1]+useFractions[1][1][2]); useFractions[2][1][0] = useFractions[1][2][0] = 1.0 - O2 - 2.0 * (O2 - O3) - O3; useFractions[2][1][1] = useFractions[1][2][1] = O2 + 2.0 * (O2 - O3); useFractions[2][1][2] = useFractions[1][2][2] = O3; printf("useFractions: %1d %1d %6.4f %6.4f %6.4f sum: %6.4f\n", 2, 1, useFractions[2][1][0], useFractions[2][1][1], useFractions[2][1][2], useFractions[2][1][0]+useFractions[2][1][1]+useFractions[2][1][2]); */ useFractions[2][2][0] = -1.0; useFractions[2][2][1] = 0.0; useFractions[2][2][2] = 0.0; useFractions[2][1][0] = -2.0; useFractions[2][1][1] = 0.0; useFractions[2][1][2] = 0.0; useFractions[1][2][0] = -3.0; useFractions[1][2][1] = 0.0; useFractions[1][2][2] = 0.0; useFractions[1][1][0] = -4.0; useFractions[1][1][1] = 0.0; useFractions[1][1][2] = 0.0; } /* end of 'first'-section */ sigma = REFL->fwhmt * arcmin2rad / 2.37; /* Rocking curve standard deviation rad */ theta_B = FACET[index].thetaB; /* Nominal Bragg angle at center of crystal rad */ E_B = iOrder * FACET[index].E[iOpt]; /* Nominal Bragg energy for this crystal keV */ dE_B = iOrder * FACET[index].dE[iOpt]; /* Nominal energy bandpass for this crystal keV */ theta = theta_B * E_B / E; /* Actual Bragg angle at specified energy rad */ d_theta = theta - theta_B; /* Deviation of actual from nominal angle rad */ dev_sig = fabs(d_theta / sigma); /* Deviation in units of the standard deviation */ dis = foc_len * tan(2.0 * d_theta); /* Displacement of focal circle from nominal position cm */ Z = FACET[index].Z; GaussFactor = exp(-dev_sig * dev_sig * 0.5); detector_ovlap(dis, detector_radius, facet_radius, &overlap_frac); if (overlap_frac <= 0.0) return(B_contr); detector_ovlap(dis, detector_radius, facet_hole_radius, &overlap_frac_h); A_final = A_basic * overlap_frac - A_hole * overlap_frac_h; /* Actual overlap area of crystal image and detector cm2 */ /* printf("A_final: %5.1f, overlap_frac: %5.3f, overlap_hole: %5.3f, Option: %2d, iOrder: %1d, dis: %5.2f, rad: %5.1f %5.1f %5.1f cm\n", A_final, overlap_frac, overlap_frac_h, iOpt, iOrder, dis, detector_radius, facet_radius, facet_hole_radius); */ /* get crystal reflectivity at specific energy */ /* status = reflect3A(E, REFL, &(FACET[index]), MATER, iOrder, refl_E); if (status < 0.0) return(0.0); B_contr = A_final * *refl_E; return(B_contr); */ /* get crystal reflectivity at specific energy (reflectivity calculated for 1, 2 and 3 times the single layer absorbtion */ status = reflect3B(E, REFL, &(FACET[index]), MATER, iOrder, refl_E3); if ((refl_E3[0] < 0.0) || (refl_E3[0] > 1.0)) { printf("R3status: %3d, E: %6.1f, dev_sig: %5.2f, iOrder: %1d, refl_E: %6.4f %6.4f %6.4f\n", status, E, dev_sig, iOrder, refl_E3[0], refl_E3[1], refl_E3[2]); exit(0); } /* if ((refl_E3[0]+refl_E3[1]+refl_E3[2]) > 0.5) { printf("R3status: %3d, E: %6.1f, refl_E: %6.4f, iOrder: %1d, thick: %5.2f, refl_E: %6.4f %6.4f %6.4f\n", status, E, *refl_E, iOrder, FACET[index].thick, refl_E3[0], refl_E3[1], refl_E3[2]); } */ if (status < 0.0) {t9++; return(B_contr);} /* determine neighbor configuration for this crystal */ numnabo[0] = numnabo[1] = 0; for (i=0; i<6; i++) { /* loop over (potential) neighbors in the two secondary layers */ inabo = i % 3; ilayer = (i - inabo) / 3; nabo_index = FACET[index].nabo[ilayer][inabo]; if ((nabo_index < -1) || (nabo_index > numLayers*N_CryPos)) { printf("nAbo: %6d, LEN: %6d \n", nabo_index, LEN[layer]); exit(0); } if (nabo_index >= 0) numnabo[ilayer]++; } if (((7+numnabo[0] + numnabo[1] * 4) < 0) || ((7+numnabo[0] + numnabo[1] * 4) > 29)) { printf("Index problem numnabo[0]: %6d, numnabo[1]: %6d\n", numnabo[0], numnabo[1]); exit(0); } if (useFractions[numnabo[0]][numnabo[1]][0] > 0.0) goto ready; /* find maximum distance between a facet in layer A and one in layer B */ for (i=0; i 1.5 * facet_sep) { /* symmetric configuration */ useFractions[2][2][0] = 1.0 - 4.0 * (O2 - O3) - 2.0 * O3; useFractions[2][2][1] = 4.0 * (O2 - O3); useFractions[2][2][2] = 2.0 * O3; goto ready; } else { /* asymmetric configuration */ useFractions[2][2][0] = 1.0 - 2.0 * (O2 - O3) - 2.0 * O2b - 3.0 * O3; useFractions[2][2][1] = 2.0 * (O2 - O3) + 2.0 * O2b; useFractions[2][2][2] = 3.0 * O3; goto ready; } } if (useFractions[numnabo[0]][numnabo[1]][0] == -4.0) { /* one neighbor in both layers */ if (distAB > 1.5 * facet_sep) { /* symmetric configuration */ useFractions[1][1][0] = 1.0 - 2.0 * O2; useFractions[1][1][1] = 2.0 * O2; useFractions[1][1][2] = 0.0; goto ready; } else { /* asymmetric configuration */ useFractions[1][1][0] = 1.0 - 2.0 * (O2 -O3) - O3; useFractions[1][1][1] = 2.0 * (O2 - O3); useFractions[1][1][2] = O3; goto ready; } } /* two neighbors in one layer, one in the other */ if (distAB > 1.5 * facet_sep) { /* symmetric configuration */ useFractions[2][1][0] = useFractions[1][2][0] = 1.0 - 2.0 * (O2 - O3) - 2.0 * O3; useFractions[2][1][1] = useFractions[1][2][1] = 2.0 * (O2 - O3); useFractions[2][1][2] = useFractions[1][2][2] = 2.0 * O3; goto ready; } else { /* asymmetric configuration */ useFractions[2][1][0] = useFractions[1][2][0] = 1.0 - O2 - 2.0 * (O2 - O3) - O3; useFractions[2][1][1] = useFractions[1][2][1] = O2 + 2.0 * (O2 - O3); useFractions[2][1][2] = useFractions[1][2][2] = O3; goto ready; } ready: /* ready */ NABSTAT[thislayer][1] += refl_E3[0]; NABSTAT[thislayer][2] += refl_E3[1]; NABSTAT[thislayer][3] += refl_E3[2]; if (index != index_old) { index_old = index; NABSTAT[thislayer][0] += 1.0; NABSTAT[thislayer][7+numnabo[0] + numnabo[1] * 4] += 1.0; NABSTAT[thislayer][4] += useFractions[numnabo[0]][numnabo[1]][0]; NABSTAT[thislayer][5] += useFractions[numnabo[0]][numnabo[1]][1]; NABSTAT[thislayer][6] += useFractions[numnabo[0]][numnabo[1]][2]; } B_contr = A_final * useFractions[numnabo[0]][numnabo[1]][0] * refl_E3[0] * GaussFactor; B_contr += A_final * useFractions[numnabo[0]][numnabo[1]][1] * refl_E3[1] * GaussFactor; B_contr += A_final * useFractions[numnabo[0]][numnabo[1]][2] * refl_E3[2] * GaussFactor; *refl_E = refl_E3[0]; /* if ((iOrder == 1) && (iOpt == 5)) printf("R3status: %3d, E: %6.1f, dev_sig: %5.2f, iOrder: %1d, refl_E: %6.4f %6.4f %6.4f, GaussFactor: %6.4f, sigma: %f\n", status, E, dev_sig, iOrder, refl_E3[0], refl_E3[1], refl_E3[2], GaussFactor, sigma); */ eff_area /= 100.0; /* convert to cm2 */ A_int /= 100.0; /* convert to cm2 */ return(B_contr); }