/*******************************************************************************
*
* McStas, the neutron ray-tracing package
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Copyright 1997-2000 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: <a href="mailto:farhi@ill.fr">Emmanuel Farhi</a> and KL.
* Date: 22 Jan 2001
* Version: 0.1
* Origin: McStas 1.5/<a href="http://www.ill.fr">ILL (France)</a> (Obsolete)
* Modified by: KN, 1998 from Source_flat.comp
* Modified by: EF, 2000 from Convert_FlatE_2_Maxwell by Thomas hansen
* Modified by: EF, KL, 2001 correct focusing pb with KL method
* Modified by: EF, 2004. Obsoleted as it is not working
*
* Neutron source with Maxwellian wavelength spectrum and
* user-specified flux.
*
* %D
* The routine is a rectangular neutron source, which aims at a square target
* centered at the beam (in order to improve MC-acceptance rate).  The angular
* divergence is then given by the dimensions of the target. The neutron
* wavelength is distributed between lambda_0 - d_lambda and
* lambda_0 + d_lambda with a Maxwellian profile.
* The source flux is specified in neutrons per steradian per square cm for the
* Maxwellian distribution. Source output is in [n/s/cm^2/Angs]
* The flux profile is flat for T=0.
* OBSOLETE: wrong results ! rather use sources/Source_Maxwell_3 or Source_gen
*
* %BUGS
* Produces wrong results
*
* <b>Example:</b>
* COMPONENT Source = Source_Maxwell(
*   h        = 0.220, w        = 0.160,
*   dist     = 2.33, xw       = 0.038, yh       = 0.2,
*   lambda_0 = 2,    d_lambda = 1.5,
*   flux     = 1e13, T        = 315)
* AT (0,0,0) ABSOLUTE
*
* %P
* h:        (m)               Height of rectangle in (x,y,0) plane where
*                             neutrons are generated.
* w:        (m)               Width of rectangle
* dist:     (m)               Distance to target along z axis.
* xw:       (m)               Width(x) of target
* yh:       (m)               Height(y) of target
* lambda_0: (AA)              Mean wavelength of neutrons.
* d_lambda: (AA)              Wavelength spread of neutrons.
* flux:     (1/(cm**2*st))    Source flux
* T:        (K)               The temperature of the source (293 K)
*
* %E
*******************************************************************************/

DEFINE COMPONENT Source_Maxwell
DEFINITION PARAMETERS ()
SETTING PARAMETERS (h=0.1, w=0.1, dist=2.33, xw=0.05, yh=0.1, lambda_0=2.36, d_lambda=0.1627, flux=1e13, T=315)
OUTPUT PARAMETERS (p_in)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
%{
 double p_in;
 double lambda0, lambda02, L2P, Lmin, Lmax;
%}
INITIALIZE
%{
  double factor, delta_lambda, source_area, k;

  source_area = h*w*1e4; /* cm^2 */
  p_in = flux/mcget_ncount()*source_area;

  k  = 1.38066e-23; /* k_B */
  if (T > 0)
        lambda0  = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T);
  else
        lambda0 = 0;
  lambda02 = lambda0*lambda0;
  L2P      = 2*lambda02*lambda02;
  Lmin = lambda_0 - d_lambda; /* AAngstroem */
  Lmax = lambda_0 + d_lambda;
%}

TRACE
%{
  double lambda, v;
  double E, Maxwell, dE_dl, dE, dE1, dE2;
  double lambda2, lambda5;
  double rf,xf,yf,dx,dy;

  z=0;

  x=w*randpm1()/2;                      /* Choose point on source */
  y=h*randpm1()/2;                      /* with uniform distribution. */

  xf = 0.5*xw*randpm1();          /* Choose focusing position uniformly */
  yf = 0.5*yh*randpm1();
  dx = xf-x;
  dy = yf-y;
  rf = sqrt(dx*dx+dy*dy+dist*dist);
  p = p_in*xw*yh*dist/(rf*rf*rf);

  lambda  = lambda_0+d_lambda*randpm1();
  lambda2 = lambda*lambda;
  lambda5 = lambda2*lambda2*lambda;

  if (lambda == 0)
  {
         printf("Source_Maxwell: warning: lambda is 0\n");
         ABSORB;
  }
  v  = K2V*(2*PI/lambda);
  vz = v*dist/rf;
  vy = v*dy/rf;
  vx = v*dx/rf;

  if (lambda0 != 0)
  {
          E  = VS2E * v*v;

          dE1 = 2*PI/(V2K*SE2V*(lambda_0-d_lambda)); dE1 *= dE1;
          dE2 = 2*PI/(V2K*SE2V*(lambda_0+d_lambda)); dE2 *= dE2;

          dE = fabs(dE1-dE2);

          dE_dl  = E/lambda*2.0;
          Maxwell= L2P/lambda5*exp(-lambda02/lambda2);
          p     *= Maxwell*dE/d_lambda/dE_dl;
  }
  SCATTER;
%}

MCDISPLAY
%{
  magnify("xy");
    multiline(5,
            -w/2.0, -h/2.0, 0.0,
             w/2.0, -h/2.0, 0.0,
             w/2.0,  h/2.0, 0.0,
            -w/2.0,  h/2.0, 0.0,
            -w/2.0, -h/2.0, 0.0);
%}

END