/******************************************************************************* * * Mcstas, neutron ray-tracing package * Copyright (C) 1997-2012, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: ESS_moderator * * %I * Written by: P Willendrup and E Klinkby, February 2014, derived from K Lefmann ESS_moderator_long * Modified by: J. Saroun (saroun@ujf.cas.cz) * Origin: DTU * * A parametrised pulsed source for modelling ESS long pulses. * * %D * Produces a time-of-flight spectrum, from the ESS parameters * Chooses evenly in lambda, evenly/exponentially decaying in time * Adapted from ESS_moderator_long by: K Lefmann, 2001 * * Updates and simplified interface: *
Note that this component does not implement "engineering reality" and currently uses a coordinate system centered on the moderator assembly. An * updated moderator component which references the "Moderator focus coordinate system" will be released later during the spring of 2016. * *
Derived from ESS_moderator_long which was debugged intensively against Mezei note (4/12 2000) and VitESS @ Rencurel 2006.
*
*-----------------------------------------------
* Correction by J. Saroun, NPI Rez:
* 1) version 2015: accepts negative port angles
* 2) version 2015: weight by cosine of the port angle
* Warning: The negative beamport angle is not taken into acccount by mcplot
*
* %VALIDATION
* Mezei-modererators validated against VitESS and Mezei note (4/12 2000) @ Rencurel 2006
* Benchmarked against multiple versions of ESS moderator group simulation data 2013-2015
*
* %P
* Input parameters:
*
* yheight_c: [m] Height of the cold source
* yheight_t: [m] Height of the thermal source
* xwidth_t: [m] Edge of thermal source
* xwidth_c: [m] Width / arc-length opening of the cold source.
* Lmin: [AA] Lower edge of wavelength distribution
* Lmax: [AA] Upper edge of wavelength distribution
* dist: [m] Distance from source to focusing rectangle; at (0,0,dist)
* focus_xw: [m] Width of focusing rectangle
* focus_yh: [m] Height of focusing rectangle
* target_index: [1] relative index of component to focus at, e.g. next is +1 this is used to compute 'dist' automatically.
* n_pulses: [1] Number of pulses simulated. 0 and 1 creates one pulse. The integrated intensity is constant
* cold_frac: [1] Fraction of neutron statistics from cold source. It is implicitely assumed that supermirror allows each beamline to choose the desired fraction of cold and thermal neutrons (i.e. extreme idealization).
* tmax_multiplier: [1] Defined maximum emission time at moderator, tmax= tmax_multiplier * ESS_PULSE_DURATION. Only in combination with sourcedef="2013", "2014" or "2015"
* beamport_angle: [deg] Direction within the beamport sector (0 < angle < extraction_opening for 2014, -extraction_opening/2 < angle < extraction_opening/2 for 2015) to direct neutrons. For sourcedef="2015", the only allowed values are 5,15,...,55 degrees measured from the central point.
* extraction_opening: [deg] Width of extraction-area in degrees (60 or 120 degrees). 120 deg only in combination with sourcedef="2014" and "2015".
* acc_power: [MW] Accelerator power in MW
* sourcedef: [string] ESS source "database", values: "TDR", "2001", "2013", "2014", "2015"
* isleft: [1] Fraction of thermal neutrons generated at the "left" moderator slab in case of "2013" or "2014"
*
* %E
*******************************************************************************/
DEFINE COMPONENT ESS_moderator
DEFINITION PARAMETERS ()
SETTING PARAMETERS (isleft=0.9,
Lmin, Lmax, cold_frac=1.0, dist=0, focus_xw, focus_yh, int target_index=0, tmax_multiplier=3,yheight_c=0.12, yheight_t=0.12,
int n_pulses=1, acc_power=5, beamport_angle=-1, string sourcedef="2015", xwidth_c=0.1, xwidth_t=0.18, extraction_opening=120)
OUTPUT PARAMETERS ()
SHARE
%{
%include "ess_source-lib"
%}
DECLARE
%{
// Weighting-oriented scalars
double l_range, w_mult, w_geom, w_geom_c, w_geom_t, w_stat;
// Flag to indicate if originating from cold or thermal moderator
int cold;
// Flag to indicate if originating from left or right thermal slab
int wasleft;
// Flag to indicate if we should apply the solid angle correction from randvec
int cosine;
// Target coordinates (focusing)
double tx,ty,tz;
// Displacement of cold moderator centre from comp origin, cold moderator radius 2014 only!
double cx,cy,cz,cyl_radius;
// Vectors giving direction of the thermal slabs
double t1x,t1y,t1z,t2x,t2y,t2z;
/* Struct for transfer to ess-source-lib funcs */
ess_moderator_struct modextras;
/* Cold and thermal function pointers */
functype cold_bril;
functype thermal_bril;
/* Time-focusing parameters - currently unused */
double tfocus_width, tfocus_time, dt;
/* radius of vacated zone around the moderator */
double r_empty;
double lambda,xprime,yprime,zprime;
double abs_beamport_angle,cos_beamport_angle,sin_beamport_angle;
/* variables needed to correct for the emission surface angle*/
double cos_thermal,cos_cold, edge_thermal;
int sgnport;
%}
INITIALIZE
%{
if (Lmin>=Lmax) {
printf("ESS_moderator: %s: Unmeaningful definition of wavelength range!\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(0);
}
sgnport=(beamport_angle>0 ? 1:-1);
if (sgnport<0) {
printf("WARNING, %s: beamport angle < 0, experimental option\n", NAME_CURRENT_COMP);
}
/* ESS 2015 - variables needed to correct for the emission surface angle */
abs_beamport_angle=fabs(beamport_angle);
cos_beamport_angle=cos(beamport_angle*DEG2RAD);
sin_beamport_angle=sin(beamport_angle*DEG2RAD);
/* correction for projection along the beam / projection on the z=0 plane */
cos_thermal=cos_beamport_angle;
cos_cold=cos((abs_beamport_angle-24.24)*DEG2RAD)/cos(24.24*DEG2RAD);
printf("input to cosine %g, output %g\n",abs_beamport_angle, cos_cold);
/* cross-section between the z=const and tilted surfaces [m] */
edge_thermal=-0.0716;
r_empty=2;
/* Weight multipliers */
w_mult=acc_power/5;
w_geom=1;
w_geom_c=1;
w_geom_t=1;
w_stat=1.0/mcget_ncount();
tfocus_width=0; tfocus_time=0; dt=0;
n_pulses=(double)floor(n_pulses);
if (n_pulses == 0) n_pulses=1;
/* Calculation of location of focusing rectangle */
if (target_index && !dist) {
Coords ToTarget;
ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
coords_get(ToTarget, &tx, &ty, &tz);
printf("%s: Focusing on a window centered at [%g, %g, %g]\n", NAME_CURRENT_COMP,tx,ty,tz);
dist=sqrt(tx*tx+ty*ty+tz*tz);
} else if (target_index && !dist) {
printf("ESS_moderator: %s: Please choose to set either the dist parameter or specify a target_index.\nExit\n", NAME_CURRENT_COMP);
exit(-1);
} else {
tx=0, ty=0, tz=dist;
}
if (focus_xw < 0 || focus_yh < 0) {
printf("ESS_moderator: %s: Please specify both focus_xw and focus_yh as positive numbers.\nExit\n", NAME_CURRENT_COMP);
exit(-1);
}
if (Lmin<=0 || Lmax <=0 || dist == 0)
{
printf("ESS_moderator: %s: Check parameters (lead to Math Error).\n Avoid 0 value for {Lmin Lmax dist} \n", NAME_CURRENT_COMP);
exit(-1);
}
/* Calculate cold cylinder radius from extraction_opening */
cyl_radius = 360*xwidth_c/(2*PI*extraction_opening);
/* Is abs_beamport_angle set at a meaningful angle - otherwise recalc */
if ( abs_beamport_angle<0 || abs_beamport_angle>extraction_opening){
beamport_angle = extraction_opening/2.0;
abs_beamport_angle=beamport_angle;
printf("%s: WARNING: Resetting your beamport_angle to %g (%g/2) - central beamslot\n",NAME_CURRENT_COMP,beamport_angle,extraction_opening);
}
/* Calculate positioning of thermal slabs and cold cylinder - NOTE: overwritten in case of ESS 2014 */
t1z = cyl_radius*cos(-DEG2RAD*abs_beamport_angle);
t1x = cyl_radius*sin(-DEG2RAD*abs_beamport_angle);
t1y = 0;
/* Wing 2 (right) is at extraction_width-beamport_angle */
t2z = cyl_radius*cos(DEG2RAD*(extraction_opening-abs_beamport_angle));
t2x = cyl_radius*sin(DEG2RAD*(extraction_opening-abs_beamport_angle));
t2y = 0;
/* We want unit vectors... */
NORM(t1x,t1y,t1z);
NORM(t2x,t2y,t2z);
cx = 0; cy=0; cz=0;
/* Geometry parameters and brilliances specified via the sourcedef string */
/* ESS-TDR definition: */
if (strcasestr(sourcedef,"TDR")) {
w_mult *= ESS_SOURCE_FREQUENCY; /* Correct for frequency */
printf("Using ESS TDR brilliance\n");
cold_bril=ESS_2012_Lieutenant_cold;
thermal_bril=ESS_Mezei_thermal;
/* ESS-2001 definition: */
} else if (strcasestr(sourcedef,"2001")) {
w_mult *= ESS_SOURCE_FREQUENCY; /* Correct for frequency */
printf("Using ESS 2001 brilliance\n");
cold_bril=ESS_Mezei_cold;
thermal_bril=ESS_Mezei_thermal;
/* ESS-2013 post-TDR updated definition, before ESS "pancake": */
} else if (strcasestr(sourcedef,"2013")) {
w_mult *= acc_power; /* Is already in per-MW units */
printf("Using ESS 2013 brilliance\n");
modextras.height_c=yheight_c;
modextras.height_t=yheight_t;
modextras.tmultiplier=tmax_multiplier;
cold_bril=ESS_2013_Schoenfeldt_cold;
thermal_bril=ESS_2013_Schoenfeldt_thermal;
/* ESS-2014 definition: */
} else if (strcasestr(sourcedef,"2014")) {
if (xwidth_c!=0.23 || xwidth_t!=0.12 || extraction_opening!=120) {
fprintf(stderr,"FATAL: sourcedef %s only allows xwidth_c=0.23, xwidth_t=0.12 and extraction_opening=120 !\n",sourcedef);
exit(-1);
}
/* Specify brilliance fct.'s */
cold_bril=ESS_2014_Schoenfeldt_cold;
thermal_bril=ESS_2014_Schoenfeldt_thermal;
/* Cold moderator geometry */
/* In this mode of operation, the emission is from the x-y plane*/
cx = 0; cy=0; cz=-0.12; cyl_radius=0.12;
/* 120 degree beam extraction opening, defining thermal moderator directions */
t1x=sin(DEG2RAD*(120/2-abs_beamport_angle));t1z=cos(DEG2RAD*120/2);
t2x=-sin(DEG2RAD*120/2);t2z=cos(DEG2RAD*120/2);
t1y=0; t2y=0;
/* ESS-2015 definition: */
} else if (strcasestr(sourcedef,"2015")) {
if (xwidth_c!=0.1 || xwidth_t!=0.18 || extraction_opening!=120) {
fprintf(stderr,"FATAL: sourcedef %s only allows xwidth_c=0.1, xwidth_t=0.18 and extraction_opening=120 !\n",sourcedef);
exit(-1);
}
/* Specify brilliance fct.'s */
cold_bril=ESS_2015_Schoenfeldt_cold;
thermal_bril=ESS_2015_Schoenfeldt_thermal;
/* Moderator geometry */
cx = 0; cy=0; cz=0; cyl_radius=0;
/* 120 degree beam extraction opening */
t1x=sin(DEG2RAD*(120/2+abs_beamport_angle));t1z=cos(DEG2RAD*(120/2+abs_beamport_angle));
t2x=-sin(DEG2RAD*(120/2-abs_beamport_angle));t2z=cos(DEG2RAD*(120/2-abs_beamport_angle));
t1y=0; t2y=0;
} else {
fprintf(stderr,"FATAL: sourcedef %s is not defined!\n",sourcedef);
exit(-1);
}
modextras.height_c=yheight_c;
modextras.Width_c=xwidth_c;
modextras.Width_t=xwidth_t;
modextras.height_t=yheight_t;
modextras.tmultiplier=tmax_multiplier;
modextras.extractionangle=extraction_opening;
if ((strcasestr(sourcedef,"2014")) && abs_beamport_angle != modextras.extractionangle/2.0) {
printf("%s: WARNING: beamport_angle is not at central slot. With sourcedef %s this means application of a partially analytical brightness model for B(\\theta) \n",NAME_CURRENT_COMP,sourcedef);
}
if ((strcasestr(sourcedef,"2015")) && (abs_beamport_angle!=5.0 && abs_beamport_angle!=15.0 && abs_beamport_angle!=25.0 && abs_beamport_angle!=35.0 && abs_beamport_angle!=45.0 && abs_beamport_angle!=55.0)) {
printf("%s: FATAL: beamport_angle is not at a defined value. Allowed are: (5, 15, ...) and < 60 \n",NAME_CURRENT_COMP);
exit(-1);
}
modextras.beamportangle=abs_beamport_angle;
l_range = Lmax-Lmin;
w_geom_c = xwidth_c*yheight_c*1.0e4; /* source area correction */
w_geom_t = xwidth_t*yheight_t*1.0e4;
w_mult *= l_range; /* wavelength range correction */
%}
TRACE
%{
p=1;
double v,E,k,r,xf,yf,dx,dy,dz,w_focus,tail_flag,cor;
/* Bispectral source - choice of spectrum and initial position */
cold = ( rand01() < cold_frac );
/* Use standard McStas routine or -not for specifying emission directionality */
cosine = 0;
/* Emission geometry adapted from ESS MCNPX model, mid 2012 */
if (strcasestr(sourcedef,"2014")) {
if (cold) { //case: cold moderator
x = 0.5*randpm1()*xwidth_c;
y = 0.5*randpm1()*yheight_c;
z = 0;
cosine = 0;
w_geom = w_geom_c;
} else { //case: thermal moderator
y = 0.5*randpm1()*yheight_t;
if (rand01()