/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Source_div
*
* %I
* Written by: KL
* Date: November 20, 1998
* Modified by: KL, 8 October 2001
* Origin: Risoe
*
* Neutron source with Gaussian or uniform divergence
*
* %D
* The routine is a rectangular neutron source, which has a gaussian or uniform
* divergent output in the forward direction.
* The neutron energy is distributed between lambda0-dlambda and
* lambda0+dlambda or between E0-dE and E0+dE. The flux unit is specified
* in n/cm2/s/st/energy unit (meV or Angs).
* In the case of uniform distribution (gauss=0), angles are uniformly distributed
* between -focus_aw and +focus_aw as well as -focus_ah and +focus_ah.
* For Gaussian distribution (gauss=1), 'focus_aw' and 'focus_ah' define the
* FWHM of a Gaussian distribution. Energy/wavelength distribution is also
* Gaussian.
*
* Example: Source_div(xwidth=0.1, yheight=0.1, focus_aw=2, focus_ah=2, E0=14, dE=2, gauss=0)
*
* %VALIDATION
* Feb 2005: tested by Kim Lefmann    (o.k.)
* Apr 2005: energy distribution used in external tests of Fermi choppers (o.k.)
* Jun 2005: wavelength distribution used in external tests of velocity selectors (o.k.)
* Validated by: K. Lieutenant
*
* %BUGS
* distribution is uniform in (hor. and vert.) angle (relative to moderator normal),
* therefore not suited for large angles
*
* %P
* xwidth: [m]                        Width of source
* yheight: [m]                       Height of source
* focus_aw: [deg]                    FWHM (Gaussian) or maximal (uniform) horz. width divergence
* focus_ah: [deg]                    FWHM (Gaussian) or maximal (uniform) vert. height divergence
* E0: [meV]                          Mean energy of neutrons.
* dE: [meV]                          Energy half spread of neutrons.
* lambda0: [Ang]                     Mean wavelength of neutrons (only relevant for E0=0)
* dlambda: [Ang]                     Wavelength half spread of neutrons.
* gauss: [0|1]                       Criterion: 0: uniform, 1: Gaussian distributions
* flux: [1/(s cm 2 st energy_unit)]  flux per energy unit, Angs or meV
*
* CALCULATED PARAMETERS:
* sigmah: [rad]                      parameter 'sigma' of the Gaussian distribution for horizontal divergence
* sigmav: [rad]                      parameter 'sigma' of the Gaussian distribution for vertical divergence
* p_init: [1]                        normalisation factor 1/'neutron_count'
*
* %E
*******************************************************************************/

DEFINE COMPONENT Source_div

SETTING PARAMETERS (xwidth, yheight, focus_aw, focus_ah, E0=0.0, dE=0.0, lambda0=0.0, dlambda=0.0, gauss=0, flux=1)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
DECLARE
%{
double sigmah;
double sigmav;
double p_init;
double dist;
double focus_xw;
double focus_yh;
%}
INITIALIZE
%{
sigmah = DEG2RAD*focus_aw/(sqrt(8.0*log(2.0)));
  sigmav = DEG2RAD*focus_ah/(sqrt(8.0*log(2.0)));

  if (xwidth < 0 || yheight < 0 || focus_aw < 0 || focus_ah < 0) {
      printf("Source_div: %s: Error in input parameter values!\n"
             "ERROR       Exiting\n",
           NAME_CURRENT_COMP);
      exit(0);
  }
  if ((!lambda0 && !E0 && !dE && !dlambda)) {
    printf("Source_div: %s: You must specify either a wavelength or energy range!\n ERROR - Exiting\n",
           NAME_CURRENT_COMP);
    exit(0);
  }
  if ((!lambda0 && !dlambda && (E0 <= 0 || dE < 0 || E0-dE <= 0))
    || (!E0 && !dE && (lambda0 <= 0 || dlambda < 0 || lambda0-dlambda <= 0))) {
    printf("Source_div: %s: Unmeaningful definition of wavelength or energy range!\n ERROR - Exiting\n",
           NAME_CURRENT_COMP);
      exit(0);
  }
  /* compute distance to next component */
  Coords ToTarget;
  double tx,ty,tz;
  ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+1),POS_A_CURRENT_COMP);
  ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
  coords_get(ToTarget, &tx, &ty, &tz);
  dist=sqrt(tx*tx+ty*ty+tz*tz);
  /* compute target area */
  if (dist) {
    focus_xw=dist*tan(focus_aw*DEG2RAD);
    focus_yh=dist*tan(focus_ah*DEG2RAD);
  }

  p_init  = flux*1e4*xwidth*yheight/mcget_ncount();
  if (!focus_aw || !focus_ah)
    exit(printf("Source_div: %s: Zero divergence defined. \n"
                "ERROR       Use non zero values for focus_aw and focus_ah.\n",
           NAME_CURRENT_COMP));
  p_init *= 2*fabs(DEG2RAD*focus_aw*sin(DEG2RAD*focus_ah/2));  /* solid angle */
  if (dlambda)
    p_init *= 2*dlambda;
  else if (dE)
    p_init *= 2*dE;
%}
TRACE
%{
  double E,lambda,v;
  double tan_h;
  double tan_v;
  double thetah;
  double thetav;

  p=p_init;
  z=0;
  t=0;

  x=randpm1()*xwidth/2.0;
  y=randpm1()*yheight/2.0;
  if(lambda0==0) {
    if (!gauss) {
      E=E0+dE*randpm1();              /*  Choose from uniform distribution */
    } else {
      E=E0+randnorm()*dE;
    }
    v=sqrt(E)*SE2V;
  } else {
    if (!gauss) {
      lambda=lambda0+dlambda*randpm1();
    } else {
      lambda=lambda0+randnorm()*dlambda;
    }
    v = K2V*(2*PI/lambda);
  }

  if (gauss==1) {
    thetah = randnorm()*sigmah;
    thetav = randnorm()*sigmav;
  } else {
    thetah = randpm1()*focus_aw*DEG2RAD/2;
    thetav = randpm1()*focus_ah*DEG2RAD/2;
  }

  tan_h = tan(thetah);
  tan_v = tan(thetav);

  /* Perform the correct treatment - no small angle approx. here! */
  vz = v / sqrt(1 + tan_v*tan_v + tan_h*tan_h);
  vy = tan_v * vz;
  vx = tan_h * vz;
%}

MCDISPLAY
%{
  
  multiline(5, -xwidth/2.0, -yheight/2.0, 0.0,
                xwidth/2.0, -yheight/2.0, 0.0,
                xwidth/2.0,  yheight/2.0, 0.0,
               -xwidth/2.0,  yheight/2.0, 0.0,
               -xwidth/2.0, -yheight/2.0, 0.0);
  if (dist) {
    dashed_line(0,0,0, -focus_xw/2,-focus_yh/2,dist, 4);
    dashed_line(0,0,0,  focus_xw/2,-focus_yh/2,dist, 4);
    dashed_line(0,0,0,  focus_xw/2, focus_yh/2,dist, 4);
    dashed_line(0,0,0, -focus_xw/2, focus_yh/2,dist, 4);
  }
%}

END