/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Monochromator_2foc * * %Identification * Written by: Peter Link. * Date: Feb. 12,1999 * Origin: Uni. Gottingen (Germany) * Modified by: Peter Link Feb. 12,1999, Added double bent feature by: * Modified by: Peter Link Sep. 24. 1999, corrected bug in rotation of v-coords: * Modified by: EF, Feb 13th 2002: Read reflectivity table * * Double bent monochromator with multiple slabs * * %Description * Double bent monochromator which uses a small-mosaicity approximation as well * as the approximation vy^2 << vz^2 + vx^2. * Second order scattering is neglected. * For an unrotated monochromator component, the crystal plane lies in the y-z * plane (ie. parallel to the beam). * 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); * * Example: Monochromator_2foc(zwidth=0.02, yheight=0.02, gap=0.0005, NH=11, NV=11, * mosaich=30, mosaicv=30, r0=0.7, Q=1.8734) * * Example values for lattice parameters * PG 002 DM=3.355 AA (Highly Oriented Pyrolytic Graphite) * PG 004 DM=1.607 AA * Heusler 111 DM=3.362 AA (Cu2MnAl) * CoFe DM=1.771 AA (Co0.92Fe0.08)b * Ge 311 DM=1.714 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: [columns] number of slabs horizontal * NV: [rows] number of slabs vertical * mosaich: [] Horisontal mosaic FWHM (arc minutes) * mosaicv: [] Vertical mosaic FWHM (arc minutes) * r0: [1] Maximum reflectivity * 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: [Angstrom] monochromator d-spacing instead of Q=2*pi/DM * mosaic: [] sets mosaich=mosaicv (arc minutes) * width: [m] total width of monochromator * height: [m] total height of monochromator * verbose: [0/1] verbosity level * reflect: [k, R] reflectivity file name of text file as 2 columns * * %Link * Additional note from Peter Link. * * %End *******************************************************************************/ DEFINE COMPONENT Monochromator_2foc DEFINITION PARAMETERS () SETTING PARAMETERS (string reflect=0, 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, RV=0, RH=0, DM=0, mosaic=0, width=0, height=0, verbose=0) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "read_table-lib" #ifndef DIV_CUTOFF #define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */ #endif %} DECLARE %{ double mos_y; /* mosaic - in radians */ double mos_z; double mono_Q; double SlabWidth; double SlabHeight; t_Table rTable; %} INITIALIZE %{ if (mosaic != 0) { mos_y = mosaic; mos_z = mos_y; } else { mos_y = mosaich; mos_z = mosaicv; } mono_Q = Q; if (DM != 0) mono_Q = 2*PI/DM; if (mono_Q == 0) { fprintf(stderr,"Monochromator_2foc: %s: Error scattering vector Q = 0\n", NAME_CURRENT_COMP); exit(-1); } if (r0 == 0) { fprintf(stderr,"Monochromator_2foc: %s: Error reflectivity r0 is null\n", NAME_CURRENT_COMP); exit(-1); } if (NH*NV == 0) { fprintf(stderr,"Monochromator_2foc: %s: no slabs ??? (NH or NV=0)\n", NAME_CURRENT_COMP); exit(-1); } if (verbose) { printf("Monochromator_2foc: 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) { if (verbose) fprintf(stdout, "Monochromator_2foc: %s : Reflectivity data (k, R)\n", NAME_CURRENT_COMP); 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); Table_Info(rTable); } else rTable.data = NULL; if (width == 0) SlabWidth = zwidth; else SlabWidth = (width+gap)/NH - gap; if (height == 0) SlabHeight = yheight; else SlabHeight = (height+gap)/NV - gap; %} TRACE %{ double dt; if(vx != 0.0 && (dt = -x/vx) >= 0.0) { double zmin,zmax, ymin,ymax, zp,yp, y1,z1,t1; zmax = ((NH*(SlabWidth+gap))-gap)/2; zmin = -1*zmax; ymax = ((NV*(SlabHeight+gap))-gap)/2; ymin = -1*ymax; y1 = y + vy*dt; /* Propagate to crystal plane */ z1 = z + vz*dt; t1 = t + dt; zp = fmod ( (z1-zmin),(SlabWidth+gap) ); yp = fmod ( (y1-ymin),(SlabHeight+gap) ); /* hit a slab or a gap ? */ if (z1>zmin && z1ymin && y1= 1) { #ifndef OPENACC if (verbose) fprintf(stdout, "Warning: Monochromator_2foc: %s: lowered reflectivity from %f to 0.99 (k=%f)\n", NAME_CURRENT_COMP, my_r0, k); #endif my_r0=0.99; } if (my_r0 < 0) { #ifndef OPENACC if (verbose) fprintf(stdout, "Warning: Monochromator_2foc: %s: raised reflectivity from %f to 0 (k=%f)\n", NAME_CURRENT_COMP, my_r0, k); #endif my_r0=0; } x = 0.0; y = y1; z = z1; t = t1; /* reflectivity */ t1 = fabs(my_r0)*exp(-tmp3*tmp3*4*log(2)); if (t1 <= 0) ABSORB; if (t1 > 1) t1 = 1; p *= t1; /* Use mosaics */ tmp1 = 2*theta; cs = cos(tmp1); sn = sin(tmp1); tmp2 = cs*vxp - sn*vzp; vyp = vyp; /* vz = cs*vz + sn*vx; diese Zeile wurde durch die folgende ersetzt */ tmp4 = vyp/vzp; /* korrigiert den schr�en Einfall aufs Pl�tchen */ vzp = cs*(-vyp*sin(tmp4)+vzp*cos(tmp4)) + sn*vxp; vxp = tmp2; /* Second: scatering out of plane. Approximation is that Debye-Scherrer cone is a plane */ phi = atan2(vyp,vzp); /* out-of plane angle */ dphi = (MIN2RAD*mos_z)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */ /* Vertical angle of the crystallite */ vyp = vzp*tan(phi+2*dphi*sin(theta)); vratio = v/sqrt(vxp*vxp+vyp*vyp+vzp*vzp); vzp = vzp*vratio; vyp = vyp*vratio; /* Renormalize v */ vxp = vxp*vratio; /* rotate v coords back */ vx = vxp*csb*csa-vyp*snb*csa+vzp*sna; vy = vxp*snb+vyp*csb; vz = -vxp*csb*sna+vyp*snb*sna+vzp*csa; /* v=sqrt(vx*vx+vy*vy+vz*vz); */ SCATTER; } /* end if Bragg ok */ } /* End intersect the crystal (if z1) */ } /* End neutron moving towards crystal (if vx)*/ %} MCDISPLAY %{ int ih; for(ih = 0; ih < NH; ih++) { int iv; for(iv = 0; iv < NV; iv++) { double zmin,zmax,ymin,ymax; double xt, xt1, yt, yt1; 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 = zmin*zmin/RH; xt1 = zmax*zmax/RH; } else { xt = 0; xt1 = 0; } if (RV) { yt = ymin*ymin/RV; yt1 = ymax*ymax/RV; } else { yt = 0; yt1 = 0; } multiline(5, xt+yt, (double)ymin, (double)zmin, xt+yt1, (double)ymax, (double)zmin, xt1+yt1, (double)ymax, (double)zmax, xt1+yt, (double)ymin, (double)zmax, xt+yt, (double)ymin, (double)zmin); } } %} END