/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: ESS_moderator_short
*
* %I
* Written by: KL, February 2001
* Modified by: KL, 18 November 2001
* Version: $Revision: 1.21 $
* Origin: Risoe
* Release: McStas 1.6
*
* A parametrised pulsed source for modelling ESS short pulses.
*
* %D
* Produces a time-of-flight spectrum, from the ESS parameters
* Chooses evenly in lambda, exponentially decaying in time .
* Adapted from Moderator by: KN, M.Hagen, August 1998
*
* Units of flux: n/cm^2/s/AA/ster
* (McStas units are in general neutrons/second)
*
* Example general parameters (general):
* size=0.12 l_low=0.1 l_high=10 dist=2 xw=0.06 yh=0.12 freq=50
* branchframe=0.5
*
* Example moderator specific parameters
* (From F. Mezei, "ESS reference moderator characteristics for ...", 4/12/00:
* Defining the normalised Maxwelian
* M(lam,T) = 2 a^2 lam^-5 exp(-a/lam^2); a=949/T; lam in AA; T in K,
* and the normalised pulse shape function
* F(t,tau,n) = ( exp(-t/tau) - exp(-nt/tau) ) n/(n-1)/tau,
* the flux distribution is given as
* Phi(t,lam) = I0 M(lam,T) F(t,tau,n)
* + I2/(1+exp(chi2 lam-2.2))/lam*F(t,tau2*lam,n2) )
*
* a1: Ambient H20, short pulse, decoupled poisoned
* T=325 tau=22e-6 tau1=0 tau2=7e-6 n=5 n2=5 chi2=2.5
* I0=9e10 I2=4.6e10 branch1=0 branch2=0.5
*
* a2: Ambient H20, short pulse, decoupled un-poisoned
* T=325 tau=35e-6 tau1=0 tau2=12e-6 n=5 n2=5 chi2=2.5
* I0=1.8e11 I2=9.2e10 branch1=0 branch2=0.5
*
* a3: Ambient H20, short pulse, coupled
* T=325 tau=80e-6 tau1=400e-6 tau2=12e-6 n=20 n2=5 chi2=2.5
* I0=4.5e11 I2=9.2e10 branch1=0.5 branch2=0.5
*
* b1: Liquid H2, short pulse, decoupled poisoned
* T=50 tau=49e-6 tau1=0 tau2=7e-6 n=5 n2=5 chi2=0.9
* I0=2.7e10 I2=4.6e10 branch1=0 branch2=0.5
*
* b2: Liquid H2, short pulse, decoupled un-poisoned
* T=50 tau=78e-6 tau1=0 tau2=12e-6 n=5 n2=5 chi2=0.9
* I0=5.4e10 I2=9.2e10 branch1=0 branch2=0.5
*
* b3: Liquid H2, short pulse, coupled
* T=50 tau=287e-6 tau1=0 tau2=12e-6 n=20 n2=5 chi2=0.9
* I0=2.3e11 I2=9.2e10 branch1=0 branch2=0.5
*
* %P
* Input parameters:
*
* size: (m) Edge of cube shaped source
* l_low: (AA) Lower edge of lambda distribution
* l_high: (AA) Upper edge of energy distribution
* dist: (m) Distance from source to focusing rectangle; at (0,0,dist)
* xw: (m) Width of focusing rectangle
* yh: (m) Height of focusing rectangle
* freq: (Hz) Frequency of pulses
* T: (K) Temperature of source
* tau: (s) long time decay constant for pulse, 1a
* tau1: (s) long time decay constant for pulse, 1b
* (only for coupled water, else 0)
* tau2: (s/AA) long time decay constant for pulse, 2
* n: (1) pulse shape parameter 1
* n2: (1) pulse shape parameter 2
* chi2: (1/AA) lambda-distribution parameter in pulse 2
* I0: (flux) integrated flux 1 (in flux units, see above) (default 9e10)
* I2: (flux*AA) integrated flux 2 (default 4.6e10)
* branch1: (1) limit for switching between distribution 1 and 2.
* (has effect only for coupled water, tau1>0)
* branch2: (1) limit for switching between distribution 1 and 2.
* (default value 0.5)
* branchframe: (1) limit for switching between 1st and 2nd pulse
* (if only one pulse wanted: 1)
*
* %E
*******************************************************************************/
DEFINE COMPONENT ESS_moderator_short
DEFINITION PARAMETERS ()
SETTING PARAMETERS (size, l_low, l_high, dist, xw, yh,
freq, T, tau, tau1, tau2, n, n2, chi2, I0, I2,
branch1, branch2, branchframe)
OUTPUT PARAMETERS (M, F, l_range, w_mult)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
double l_range, w_mult, Delta_t ;
double M(double l, double temp)
{
double a=949.0/temp;
return 2*a*a*exp(-a/(l*l))/(l*l*l*l*l);
}
double F(double t, double tau, int n)
{
return (exp(-t/tau)-exp(-n*t/tau))*n/(n-1)/tau;
}
%}
INITIALIZE
%{
if (n == 1 || n2 == 1 || l_low<=0 || l_high <=0 || dist == 0
|| branch1 == 1 || branch2 == 1 || tau == 0)
{
printf("ESS_moderator_short: %s: Check parameters (lead to Math Error).\n Avoid 0 value for {l_low l_high dist tau} and 1 value for {n n2 branch1/2/frame}\n", NAME_CURRENT_COMP);
exit(-1);
}
Delta_t = 1.0/(double)freq;
l_range = l_high-l_low;
w_mult = size*size*1.0e4; /* source area correction */
w_mult *= l_range; /* wavelength range correction */
w_mult *= 1/mcget_ncount(); /* Flux value */
%}
TRACE
%{
double v,tau_l,E,lambda,k,r,xf,yf,dx,dy,w_focus;
z=0;
x = 0.5*size*randpm1();
y = 0.5*size*randpm1(); /* Choose initial position */
randvec_target_rect(&xf, &yf, &r, &w_focus,
0, 0, dist, xw, yh, ROT_A_CURRENT_COMP);
dx = xf-x;
dy = yf-y;
r = sqrt(dx*dx+dy*dy+dist*dist);
lambda = l_low+l_range*rand01(); /* Choose from uniform distribution */
k = 2*PI/lambda;
v = K2V*k;
vz = v*dist/r;
vy = v*dy/r;
vx = v*dx/r;
/* printf("pos0 (%g %g %g), pos1 (%g %g %g), r: %g, v (%g %g %g), v %g\n",
x,y,z,xf,yf,dist,r,vx,vy,vz, v);
printf("l %g, w_focus %g \n", lambda, w_focus); */
if (rand01() < branch2)
{
if (tau1>0) /* coupled water */
if (rand01() < branch1)
{ /* FIRST CASE a */
tau_l = tau;
p = 1/(branch1*branch2); /* Correct for switching prob. */
}
else
{ /* FIRST CASE b */
tau_l = tau1;
p = 1/((1-branch1)*branch2); /* Correct for switching prob. */
}
else /* all other moderators */
{
tau_l = tau;
p = 1/branch2; /* Correct for switching prob. */
}
t = -tau_l*log(1e-12+rand01()); /* Sample from long-time tail a */
p *= n/(n-1)*(1-exp(-(n-1)*t/tau_l)); /* Correct for true pulse shape */
p *= w_focus; /* Correct for target focusing */
p *= I0*w_mult*M(lambda,T); /* Calculate true intensity */
}
else
{
/* SECOND CASE */
tau_l = tau2*lambda; /* CHECK THIS */
t = -tau_l*log(1e-12+rand01()); /* Sample from long-time tail */
p = n2/(n2-1)*(1-exp(-(n2-1)*t/tau_l));/* Correct for true pulse shape */
p /= (1-branch2); /* Correct for switching prob. */
p *= w_focus; /* Correct for target focusing */
p *= I2*w_mult/(1+exp(chi2*lambda-2.2))/lambda;
/* Calculate true intensity */
}
if (branchframe > 0 && rand01() <= branchframe)
{
p /= branchframe;
}
else
{
if (rand01() < 0.5)
t += Delta_t;
else
t -= Delta_t;
p *= 2.0/(1-branchframe);
}
%}
MCDISPLAY
%{
magnify("xy");
multiline(5, -(double)dist/2.0, -(double)dist/2.0, 0.0,
(double)dist/2.0, -(double)dist/2.0, 0.0,
(double)dist/2.0, (double)dist/2.0, 0.0,
-(double)dist/2.0, (double)dist/2.0, 0.0,
-(double)dist/2.0, -(double)dist/2.0, 0.0);
%}
END