/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2007, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Monochromator_curved * * %I * * Written by: Emmanuel Farhi, Kim, Lefmann, Peter Link * Date: Aug. 24th 2001 * Origin: ILL * Modified by: EF, Aug. 24th 2001: From Mosaic_anisotropic and Mon_2foc * Modified by: EF, Feb 13th 2002: Read reflectivity table * Modified by: EF, Oct 24th 2002: Read transmission table * * Double bent multiple crystal slabs with anisotropic gaussian mosaic. * * %Description * Double bent infinitely thin mosaic crystal, useful as a monochromator or * analyzer. which uses a small-mosaicity approximation and taking into account * higher order scattering. The mosaic is anisotropic gaussian, with different * FWHMs in the Y and Z directions. The scattering vector is perpendicular to the * surface. For an unrotated monochromator component, the crystal plane lies in * the y-z plane (ie. parallel to the beam). The component works in reflection, but * also transmits the non-diffracted beam. Reflectivity and transmission files may * be used. The slabs are positioned in the vertical plane (not on a * cylinder/sphere), and are rotated according to the curvature radius. * When curvatures are set to 0, the monochromator is flat. * The curvatures approximation for parallel beam focusing to distance L, with * monochromator rotation angle A1 are: * RV = 2*L*sin(DEG2RAD*A1); * RH = 2*L/sin(DEG2RAD*A1); * * When you rotate the component by A1 = asin(Q/2/Ki)*RAD2DEG, do not forget to * rotate the following components by A2=2*A1 (for 1st order) ! * * Example: Monochromator_curved(zwidth=0.01, yheight=0.01, gap=0.0005, NH=11, NV=11, mosaich=30.0, mosaicv=30.0, r0=0.7, Q=1.8734) * * Monochromator lattice parameter * PG 002 DM=3.355 AA (Highly Oriented Pyrolythic Graphite) * PG 004 DM=1.677 AA * Heusler 111 DM=3.362 AA (Cu2MnAl) * CoFe DM=1.771 AA (Co0.92Fe0.08) * Ge 111 DM=3.266 AA * Ge 311 DM=1.714 AA * Ge 511 DM=1.089 AA * Ge 533 DM=0.863 AA * Si 111 DM=3.135 AA * Cu 111 DM=2.087 AA * Cu 002 DM=1.807 AA * Cu 220 DM=1.278 AA * Cu 111 DM=2.095 AA * * %Parameters * INPUT PARAMETERS: * * zwidth: [m] horizontal width of an individual slab * yheight: [m] vertical height of an individual slab * gap: [m] typical gap between adjacent slabs * NH: [int] number of slabs horizontal * NV: [int] number of slabs vertical * mosaich: [arc minutes] Horisontal mosaic FWHM * mosaicv: [arc minutes] Vertical mosaic FWHM * r0: [1] Maximum reflectivity. O unactivates component * Q: [AA-1] Scattering vector * RV: [m] radius of vertical focussing. flat for 0 * RH: [m] radius of horizontal focussing. flat for 0 * * optional parameters * DM: [AA] monochromator d-spacing instead of Q=2*pi/DM * mosaic: [arc minutes] sets mosaich=mosaicv * width: [m] total width of monochromator, along Z * height: [m] total height of monochromator, along Y * verbose: [0/1] verbosity level * reflect: [str] reflectivity file name of text file as 2 columns [k, R] * transmit: [str] transmission file name of text file as 2 columns [k, T] * t0: [1] transmission efficiency * order: [1] specify the diffraction order, 1 is usually prefered. Use 0 for all * * %Link * Additional note from Peter Link. * %L * Obsolete Mosaic_anisotropic by Kristian Nielsen * %L * Contributed Monochromator_2foc by Peter Link * * %End *******************************************************************************/ DEFINE COMPONENT Monochromator_curved DEFINITION PARAMETERS () SETTING PARAMETERS (string reflect="NULL", string transmit="NULL", zwidth=0.01, yheight=0.01,gap=0.0005,int NH=11, int NV=11, mosaich=30.0, mosaicv=30.0, r0=0.7, t0=1.0,Q=1.8734, RV=0, RH=0, DM=0, mosaic=0, width=0, height=0, verbose=0, order=0) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ #pragma acc routine double GAUSS_monocurved(double x, double mean, double rms){ return (exp(-((x)-(mean))*((x)-(mean))/(2*(rms)*(rms)))/(sqrt(2*PI)*(rms))); } %include "read_table-lib" %} DECLARE %{ double mos_rms_y; /* root-mean-square of mosaic, in radians */ double mos_rms_z; double mos_rms_max; double mono_Q; double SlabWidth; double SlabHeight; t_Table rTable; t_Table tTable; int rTableFlag; int tTableFlag; double *tiltH; double *tiltV; %} INITIALIZE %{ int i; if (mosaic != 0) { mos_rms_y = MIN2RAD*mosaic/sqrt(8*log(2)); mos_rms_z = mos_rms_y; } else { mos_rms_y = MIN2RAD*mosaich/sqrt(8*log(2)); mos_rms_z = MIN2RAD*mosaicv/sqrt(8*log(2)); } mos_rms_max = mos_rms_y > mos_rms_z ? mos_rms_y : mos_rms_z; mono_Q = Q; if (DM != 0) mono_Q = 2*PI/DM; if (mono_Q <= 0) { fprintf(stderr,"Monochromator_curved: %s: Error scattering vector Q = 0\n", NAME_CURRENT_COMP); exit(-1); } if (r0 < 0) { fprintf(stderr,"Monochromator_curved: %s: Error reflectivity r0 is negative\n", NAME_CURRENT_COMP); exit(-1); } if (r0 == 0) { fprintf(stderr,"Monochromator_curved: %s: Reflectivity r0 is null. Ignoring component.\n", NAME_CURRENT_COMP); } if (NH*NV == 0) { fprintf(stderr,"Monochromator_curved: %s: no slabs ??? (NH or NV=0)\n", NAME_CURRENT_COMP); exit(-1); } if (verbose && r0) { printf("Monochromator_curved: component %s Q=%.3g Angs-1 (DM=%.4g Angs)\n", NAME_CURRENT_COMP, mono_Q, 2*PI/mono_Q); if (NH*NV == 1) printf(" flat.\n"); else { if (NH > 1) { printf(" horizontal: %i blades", (int)NH); if (RH != 0) printf(" focusing with RH=%.3g [m]", RH); printf("\n"); } if (NV > 1) { printf(" vertical: %i blades", (int)NV); if (RV != 0) printf(" focusing with RV=%.3g [m]", RV); printf("\n"); } } } if (reflect != NULL && r0 && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) { if (verbose) fprintf(stdout, "Monochromator_curved: %s: Reflectivity data (k, R) from %s\n", NAME_CURRENT_COMP, reflect); Table_Read(&rTable, reflect, 1); /* read 1st block data from file into rTable */ Table_Rebin(&rTable); /* rebin as evenly, increasing array */ if (rTable.rows < 2) Table_Free(&rTable); if (verbose) Table_Info(rTable); rTableFlag = 1; } else { rTableFlag = 0; } if (transmit != NULL && strlen(transmit) && strcmp(transmit,"NULL") && strcmp(transmit,"0")) { if (verbose) fprintf(stdout, "Monochromator_curved: %s: Transmission data (k, T) from %s\n", NAME_CURRENT_COMP, transmit); Table_Read(&tTable, transmit, 1); /* read 1st block data from file into rTable */ Table_Rebin(&tTable); /* rebin as evenly, increasing array */ if (tTable.rows < 2) Table_Free(&tTable); if (verbose) Table_Info(tTable); tTableFlag = 1; } else { tTableFlag = 0; } if (width == 0) SlabWidth = zwidth; else SlabWidth = (width+gap)/NH - gap; if (height == 0) SlabHeight = yheight; else SlabHeight = (height+gap)/NV - gap; tiltH=calloc((int)2*(NH+1),sizeof(double)); tiltV=calloc((int)2*(NV+1),sizeof(double)); if (!tiltH) printf("Monochromator_curved: %s: Warning: not enough memory to allocate tilts (NH=%g).\n", NAME_CURRENT_COMP, NH); else if (RH) { /* pre-compute tilts */ for (i=0;i<=NH;i++){ tiltH[i]=asin((i-(NH+1)/2.0)*(SlabWidth+gap)/RH); } } if (!tiltV) printf("Monochromator_curved: %s: Warning: not enough memory to allocate tilts (NV=%g).\n", NAME_CURRENT_COMP, NV); else if (RV) { for (i=0;i<=NV;i++){ tiltV[i]=-asin((i-(NV+1)/2.0)*(SlabHeight+gap)/RV); } } %} TRACE %{ double dt; double Gauss_X[] = {-0.987992518020485, -0.937273392400706, -0.848206583410427, -0.724417731360170, -0.570972172608539, -0.394151347077563, -0.201194093997435, 0, 0.201194093997435, 0.394151347077563, 0.570972172608539, 0.724417731360170, 0.848206583410427, 0.937273392400706, 0.987992518020485}; double Gauss_W[] = {0.030753241996117, 0.070366047488108, 0.107159220467172, 0.139570677926154, 0.166269205816994, 0.186161000115562, 0.198431485327111, 0.202578241925561, 0.198431485327111, 0.186161000115562, 0.166269205816994, 0.139570677926154, 0.107159220467172, 0.070366047488108, 0.030753241996117}; if(vx != 0.0 && (dt = -x/vx) >= 0.0 && r0) { /* Moving towards crystal? */ double zmin,zmax, ymin,ymax; double yy,zz; zmax = ((NH*(SlabWidth+gap))-gap)/2; zmin = -zmax; ymax = ((NV*(SlabHeight+gap))-gap)/2; ymin = -ymax; /* Test-propagate to crystal plane */ zz=z+vz*dt; yy=y+vy*dt; if (zz>zmin && zzymin && yy 2*k/mono_Q) Q_order--; if((!order && Q_order > 0) || (Q_order == fabs(order) && order)) { /* Bragg scattering possible? */ double q0, q0x, theta, delta, p_reflect, my_r0; q0 = Q_order*mono_Q; q0x = ratio < 0 ? -q0 : q0; theta = asin(q0/(2*k)); /* Actual bragg angle */ /* Make MC choice: reflect or transmit? */ delta = asin(fabs(kux)) - theta; if (rTableFlag) { my_r0 = r0*Table_Value(rTable, k, 1); /* 2nd column */ } else my_r0 = r0; if (my_r0 > 1) { if (my_r0 > 1.01 && verbose) fprintf(stdout, "Warning: Monochromator_curved : lowered reflectivity from %f to 1 (k=%f)\n", my_r0, k); my_r0=0.999; } if (my_r0 < 0) { if (verbose) fprintf(stdout, "Warning: Monochromator_curved : raised reflectivity from %f to 0 (k=%f)\n", my_r0, k); my_r0=0; } p_reflect = fabs(my_r0)*exp(-kiz*kiz/(kiy*kiy + kiz*kiz)*(delta*delta)/ (2*mos_rms_y*mos_rms_y))* exp(-kiy*kiy/(kiy*kiy + kiz*kiz)*(delta*delta)/ (2*mos_rms_z*mos_rms_z)); double rr=rand01(); if(rr <= p_reflect) { /* Reflect */ double bx,by,bz,ax,ay,az,phi; double cos_2theta,k_sin_2theta,cos_phi,sin_phi,q_x,q_y,q_z; double total,c1x,c1y,c1z,w,mos_sample; int i=0; cos_2theta = cos(2*theta); k_sin_2theta = k*sin(2*theta); /* Get unit normal to plane containing ki and most probable kf */ vec_prod(bx, by, bz, kix, kiy, kiz, q0x, 0, 0); NORM(bx,by,bz); bx = bx * k_sin_2theta; by = by * k_sin_2theta; bz = bz * k_sin_2theta; /* Get unit vector normal to ki and b */ vec_prod(ax, ay, az, bx, by, bz, kux, kuy, kuz); /* Compute the total scattering probability at this ki */ total = 0; /* Choose width of Gaussian distribution to sample the angle * phi on the Debye-Scherrer cone for the scattered neutron. * The radius of the Debye-Scherrer cone is smaller by a * factor 1/cos(theta) than the radius of the (partial) sphere * describing the possible orientations of Q due to mosaicity, so we * start with a width 1/cos(theta) greater than the largest of * the two mosaics. */ mos_sample = mos_rms_max/cos(theta); c1x = kix*(cos_2theta-1); c1y = kiy*(cos_2theta-1); c1z = kiz*(cos_2theta-1); /* Loop, repeatedly reducing the sample width until it is small * enough to avoid sampling scattering directions with * ridiculously low scattering probability. * Use a cut-off at 5 times the gauss width for considering * scattering probability as well as for integration limits * when integrating the sampled distribution below. */ for(i=0; i<100; i++) { w = 5*mos_sample; cos_phi = cos(w); sin_phi = sin(w); q_x = c1x + cos_phi*ax + sin_phi*bx; q_y = (c1y + cos_phi*ay + sin_phi*by)/mos_rms_z; q_z = (c1z + cos_phi*az + sin_phi*bz)/mos_rms_y; /* Stop when we get near a factor of 25=5^2. */ if(q_z*q_z + q_y*q_y < (25/(2.0/3.0))*(q_x*q_x)) break; mos_sample *= (2.0/3.0); } /* Now integrate the chosen sampling distribution, using a * cut-off at five times sigma. */ for(i = 0; i < (sizeof(Gauss_X)/sizeof(double)); i++) { phi = w*Gauss_X[i]; cos_phi = cos(phi); sin_phi = sin(phi); q_x = c1x + cos_phi*ax + sin_phi*bx; q_y = c1y + cos_phi*ay + sin_phi*by; q_z = c1z + cos_phi*az + sin_phi*bz; p_reflect = GAUSS_monocurved((q_z/q_x),0,mos_rms_y)* GAUSS_monocurved((q_y/q_x),0,mos_rms_z); total += Gauss_W[i]*p_reflect; } total *= w; /* Choose point on Debye-Scherrer cone. Sample from a Gaussian of * width 1/cos(theta) greater than the mosaic and correct for any * error by adjusting the neutron weight later. */ phi = mos_sample*randnorm(); /* Compute final wave vector kf and scattering vector q = ki - kf */ cos_phi = cos(phi); sin_phi = sin(phi); q_x = c1x + cos_phi*ax + sin_phi*bx; q_y = c1y + cos_phi*ay + sin_phi*by; q_z = c1z + cos_phi*az + sin_phi*bz; p_reflect = GAUSS_monocurved((q_z/q_x),0,mos_rms_y)* GAUSS_monocurved((q_y/q_x),0,mos_rms_z); vx = K2V*(kix+q_x); vy = K2V*(kiy+q_y); vz = K2V*(kiz+q_z); p_reflect /= total*GAUSS_monocurved(phi,0,mos_sample); if (p_reflect <= 0) ABSORB; if (p_reflect > 1) p_reflect = 1; p = p * p_reflect; } /* End MC choice to reflect or transmit neutron (if tmp 1) { if (my_t0 > 1.01 && verbose) fprintf(stdout, "Warning: Monochromator_curved : lowered transmission from %f to 1 (k=%f)\n", my_t0, k); my_t0=0.999; } if (my_t0 > 0) p = p * my_t0; else ABSORB; } } /* end if not in gap */ /* rotate back in component frame */ Rotation TT; rot_transpose(T, TT); /* now make the coordinate system change */ mccoordschange_polarisation(TT, &vx, &vy, &vz); coords_get(rot_apply(TT,coords_set(x,y,z)),&x,&y,&z); y= y+center_y; z= z+center_z; /* Visualise scattering point in proper, component frame - but only if the neutron is reflected, that is none of: * transmitted * falling outside the slab material */ if(!do_transmit) SCATTER; /* mccoordschange_polarisation(tt, &sx, &sy, &sz); */ } /* End intersect the crystal (if z) */ else { /* restore neutron state when no interaction */ RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } } /* End neutron moving towards crystal (if vx)*/ %} FINALLY %{ if(rTableFlag){ Table_Free(&rTable); } if(tTableFlag){ Table_Free(&tTable); } if (tiltH) free(tiltH); if (tiltV) free(tiltV); %} MCDISPLAY %{ int ih; for(ih = 0; ih < NH; ih++) { int iv; for(iv = 0; iv < NV; iv++) { double zmin,zmax,ymin,ymax; double xt, yt; zmin = (SlabWidth+gap)*(ih-NH/2.0)+gap/2; zmax = zmin+SlabWidth; ymin = (SlabHeight+gap)*(iv-NV/2.0)+gap/2; ymax = ymin+SlabHeight; if (RH) xt = -(zmax*zmax - zmin*zmin)/RH/2; else xt = 0; if (RV) yt = -(ymax*ymax - ymin*ymin)/RV/2; else yt = 0; multiline(5, xt+yt, (double)ymin, (double)zmin, xt-yt, (double)ymax, (double)zmin, -xt-yt, (double)ymax, (double)zmax, -xt+yt, (double)ymin, (double)zmax, xt+yt, (double)ymin, (double)zmin); } } %} END