/*******************************************************************************
*
* McStas, the neutron ray-tracing package
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Copyright 1997-2000 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: KN, M.Hagen
* Date: August 1998
* Version: $Revision: 1.1.1.1 $
* Origin: McStas release
*
* A simple pulsed source for time-of-flight.
*
* %D
* Produces a simple time-of-flight spectrum, with a flat energy distribution
*
* %P
* Input parameters:
*
* radius: (m)   Radius of source
* E0:     (meV) Lower edge of energy distribution
* E1:     (meV) Upper edge of energy distribution
* dist:   (m)   Distance from source to the focusing rectangle
* xw:     (m)   Width of focusing rectangle
* yh:     (m)   Height of focusing rectangle
* t0:     (mus) decay constant for low-energy neutrons
* Ec:     (meV) Critical energy, below which the flux decay is constant
* gam:    (meV) energy dependence of decay time
*
* %E
*******************************************************************************/

DEFINE COMPONENT Moderator
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius, E0, E1, dist, xw, yh, t0, Ec, gam)
OUTPUT PARAMETERS (hdiv,vdiv,p_in)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
  double hdiv,vdiv;
  double p_in;
%}
INITIALIZE
%{
  hdiv = atan(xw/(2.0*dist));
  vdiv = atan(yh/(2.0*dist));
  p_in = (4*hdiv*vdiv)/(4*PI);
%}
TRACE
%{
  double theta0,phi0,chi,theta,phi,v,r,tauv,E;

  p=p_in;
  z=0;

  chi = 2*PI*rand01();          /* Choose point on source */
  r = sqrt(rand01())*radius;    /* with uniform distribution. */
  x = r*cos(chi);
  y = r*sin(chi);

  theta0 = -atan(x/dist);       /* Angles to aim at target centre */
  phi0 = -atan(y/dist);

  theta = theta0 + hdiv*randpm1(); /* Small angle approx. */ 
  phi = phi0 + vdiv*randpm1();

  E = E0+(E1-E0)*rand01();      /* Assume linear distribution */
  v = SE2V*sqrt(E);

  vz = v*cos(phi)*cos(theta);   /* Small angle approx. */
  vy = v*sin(phi);
  vx = v*cos(phi)*sin(theta);

  if(E < Ec)
  {
    tauv = t0;
  }
  else
  {
    double tmp = ((E - Ec) / gam);
    tauv = t0 / (1 + (tmp*tmp));
  }
  t = -tauv*log(rand01())*1E-6;
%}

MCDISPLAY
%{
  magnify("xy");
  circle("xy",0,0,0,radius);
%}

END