/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Source_gen
*
* %I
* Written by: Emmanuel Farhi, Kim Lefmann
* Date: Aug 27, 2001
* Version: $Revision: 1.21 $
* Origin: ILL/Risoe
* Release: McStas 1.6
* Modified by: EF, Aug 27, 2001 ; can use Energy/wavelength and I1
* Modified by: EF, Sep 18, 2001 ; corrected illumination bug. Add options
* Modified by: EF, Oct 30, 2001 ; minor changes. mccompcurname replaced
*
* Circular/squared neutron source with flat or Maxwellian energy/wavelength
* spectrum (possibly spatially gaussian)
*
* %D
* This routine is a neutron source (rectangular or circular), 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 energy/wavelength is distributed between Emin=E0-dE and
* Emax=E0+dE or Lmin=Lambda0-dLambda and Lmax=Lambda0+dLambda. The I1 may
* be either arbitrary (I1=0), or specified in neutrons per steradian per
* square cm per Angstrom. A Maxwellian spectra may be selected if you give
* the source temperatures (up to 3).
* The source shape is defined by its radius, or can alternatively be squared
* if you specify non-zero h and w parameters.
* The beam is spatially uniform, but becomes gaussian if one of the source
* dimensions (radius, h or w) is negative, or you set the gaussian flag.
* Divergence profiles are triangular.
* The source may have a thickness, which will broaden the default zero time
* distribution.
* For the Maxwellian spectra, the generated intensity is dPhi/dLambda (n/s/AA)
* For flat spectra, the generated intensity is Phi (n/s).
*
* Usage example:
*   Source_gen(radius=0.1,Lambda0=2.36,dLambda=0.16, T1=20, T2=38, I1=1e13)
*   Source_gen(h=0.1,w=0.1,Emin=1,Emax=3, I1=1e13, verbose=1, gaussian=1)
*   EXTEND
*   %{
*      t = rand0max(1e-3); // set time from 0 to 1 ms for TOF instruments.
*   %}
*
* <b>Some sources parameters:</b>
* PSI cold source T1=150.42,I1=3.67e11 T2=38.74,I2=3.64e11
*                 T3=14.84,I3=0.95e11
* ILL VCS cold source T1=216.8,I1=1.24e+13,T2=33.9,I2=1.02e+13
*        (H1)         T3=16.7 ,I3=3.0423e+12
* ILL HCS cold source T1=413.5,I1=10.22e12,T2=145.8,I2=3.44e13
*        (H5)         T3=40.1 ,I3=2.78e13
* ILL Thermal tube    T1=683.7,I1=0.5874e+13,T2=257.7,I2=2.5099e+13
*        (H12)        T3=16.7 ,I3=1.0343e+12
* ILL Hot source      T1=1695,I1=1.74e13,T2=708,I2=3.9e12
*
* %VALIDATION
* Feb 2005: output cross-checked for 3 Maxwellians against VITESS source
*           I(lambda), I(hor_div), I(vert_div) identical in shape and absolute values
* Validated by: K. Lieutenant
*
* %P
* radius: [m]   Radius of circle in (x,y,0) plane where neutrons
*               are generated. You may also use 'h' and 'w' for a square source
* dist:   [m]   Distance to target along z axis.
* xw:     [m]   Width(x) of target.
*                 If dist=0.0, xw = horz. div  in deg
* yh:     [m]   Height(y) of target.
*                 If dist=0.0, yh = vert. div in deg
* E0:     [meV] Mean energy of neutrons. You may also use Lambda0.
* dE:     [meV] Energy spread of neutrons, half width.
*
* Optional parameters:
* Lambda0:  [AA]  Mean wavelength of neutrons.
* dLambda:  [AA]  Wavelength spread of neutrons,half width
* Emin:     [meV] Minimum energy of neutrons
* Emax:     [meV] Maximum energy of neutrons
* Lmin:     [AA]  Minimum wavelength of neutrons
* Lmax:     [AA]  Maximum wavelength of neutrons
* h:        [m]   Source y-height, then does not use radius parameter
* w:        [m]   Source x-width, then does not use radius parameter
* length:   [m]   Source z-length, not anymore flat
* T1:       [K]   Temperature of the Maxwellian source, 0=none
* I1:       [1/(cm**2*st*AA)] Source flux per solid angle, area and Angstrom
* T2:       [K]   Second Maxwellian source Temperature, 0=none
* I2:       [1/(cm**2*st*AA)] Second Maxwellian Source flux
* T3:       [K]   Third Maxwellian source Temperature, 0=none
* I3:       [1/(cm**2*st*AA)] Third Maxwellian Source flux
* gaussian: [0/1] Use gaussian beam, normal distributions
* verbose:  [0/1] display info about the source. -1 unactivate source.
*
* %L
* <a href="http://www.ill.fr/Yellowbook">The ILL Yellow book</a>
* %L
* P. Ageron, Nucl. Inst. Meth. A 284 (1989) 197
* %E
******************************************************************************/

DEFINE COMPONENT Source_gen
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius=0.0, dist=57.296, xw=0, yh=0, E0=14.68, dE=0, Lambda0=2.36, dLambda=0.1627, I1=0, h=0, w=0, gaussian=0, verbose=0, T1=0,Lmin=0,Lmax=0,Emin=0,Emax=0,T2=0,I2=0,T3=0,I3=0,length=0)
OUTPUT PARAMETERS (p_in,lambda0, lambda02, L2P,lambda_0,d_lambda,lambda0b, lambda02b, L2Pb,lambda0c, lambda02c, L2Pc)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
%{
 double p_in;
 double lambda_0=0,d_lambda=0;
 double lambda0, lambda02, L2P; /* first Maxwellian source */
 double lambda0b, lambda02b, L2Pb;  /* second Maxwellian source */
 double lambda0c, lambda02c, L2Pc;  /* third Maxwellian source */
%}
INITIALIZE
%{
  double factor, delta_lambda, source_area, k;

  if (radius == 0 && h == 0 && w == 0)
  {
    fprintf(stderr,"Source_gen: Error: Please precise source geometry (radius, h, w)\n");
    exit(-1);
  }
  if (xw*yh == 0)
  {
    fprintf(stderr,"Source_gen: Error: Please precise source target (xw, yh)\n");
    exit(-1);
  }
  if ((radius < 0) || (h < 0) || (w < 0))
    gaussian = 1;
  else
    gaussian = 0;

  if ((Emin != 0) && (Emax != 0))
  { E0 = (Emin+Emax)/2; dE=fabs(Emin-Emax)/2; }
  if ((E0 != 0) && (dE != 0))
  {
    Lmin = sqrt(81.81/(E0+dE)); /* AAngstroem */
    Lmax = sqrt(81.81/(E0-dE));
  }
  if ((Lmin != 0) && (Lmax != 0))
  { Lambda0 = (Lmin+Lmax)/2; dLambda=fabs(Lmin-Lmax)/2; }
  if ((Lambda0 != 0) || (dLambda != 0))
  { lambda_0 = Lambda0; d_lambda=dLambda; }

  radius = fabs(radius); w=fabs(w); h=fabs(h);  I1=fabs(I1);
  lambda_0=fabs(lambda_0); d_lambda=fabs(d_lambda);
  xw = fabs(xw); yh=fabs(yh); dist=fabs(dist);

  if ((lambda_0-d_lambda <= 0) || (lambda_0+d_lambda <= 0))
  {
    fprintf(stderr,"Source_gen: Warning: Wavelength will reach negative values\n");

  }

  if (dist == 0)
  {
    fprintf(stderr,"Source_gen: warning: focusing distance is null. xw and yh interpreted as divergence in deg\n");
  }

  delta_lambda = 2*d_lambda;
  Lmin = lambda_0 - d_lambda; /* AAngstroem */
  Lmax = lambda_0 + d_lambda;

  if (I1 != 0)
  {

    if ((h == 0) || (w == 0))
      source_area = radius*radius*PI*1e4; /* circular cm^2 */
    else
      source_area = h*w*1e4; /* square cm^2 */
    factor = I1*source_area*delta_lambda;
    p_in = factor;  /* p_in = (4*hdiv*vdiv)*factor; */
  }
  else
    p_in = 1.0/4/PI; /* Small angle approx. */
  p_in /= mcget_ncount();


  k  = 1.38066e-23; /* k_B */
  if (T1 > 0)
  {
    lambda0  = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T1);
    lambda02 = lambda0*lambda0;
    L2P      = 2*lambda02*lambda02;
  }
  else
    { lambda0 = lambda_0; }

  if (T2 > 0)
  {
    lambda0b  = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T2);
    lambda02b = lambda0b*lambda0b;
    L2Pb      = 2*lambda02b*lambda02b;
  }
  else
    { lambda0b = lambda_0; }

  if (T3 > 0)
  {
    lambda0c  = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T3);
    lambda02c = lambda0c*lambda0c;
    L2Pc      = 2*lambda02c*lambda02c;
  }
  else
    { lambda0c = lambda_0; }

  if (verbose == 1)
  {
    printf("Source_gen: component %s ", NAME_CURRENT_COMP);
    if ((h == 0) || (w == 0))
      printf("(disk)");
    else
      printf("(square)");
    printf("\n            spectra ");
    printf("%.2f to %.2f AA (%.2f to %.2f meV)", Lmin, Lmax, 81.81/Lmax/Lmax, 81.81/Lmin/Lmin);
    if (gaussian)
      printf(", gaussian beam");
    printf("\n");
    if (T1 != 0)
      printf("            T1=%.1f K (%.2f AA)", T1, lambda0);
    if (T2*I2 != 0)
      printf(", T2=%.1f K (%.2f AA)", T2, lambda0b);
    if (T3*I3 != 0)
      printf(", T3=%.1f K (%.2f AA)", T3, lambda0c);
    if (T1) printf("\n");
    if (T1 && I1) printf("  Flux is dPhi/dLambda in [n/s/AA].\n");
    else          printf("  Flux is Phi in [n/s].\n");
  }
  else
    if (verbose == -1)
      printf("Source_gen: component %s unactivated", NAME_CURRENT_COMP);
%}
TRACE
%{
  double theta0,phi0,theta1,phi1,chi,theta,phi,v,r, lambda;
  double tan_h, tan_v, Maxwell, lambda2, lambda5;

  if (verbose >= 0)
  {

    z=0;

    if ((h == 0) || (w == 0))
    {
      chi=2*PI*rand01();                          /* Choose point on source */
      r=sqrt(rand01())*radius;                    /* with uniform distribution. */
      x=r*cos(chi);
      y=r*sin(chi);
    }
    else
    {
      x = w*randpm1()/2;   /* select point on source (uniform) */
      y = h*randpm1()/2;
    }
    if (length != 0)
      z = length*randpm1()/2;

    if (dist == 0)
    { 
      theta0 =  DEG2RAD*xw; 
      phi0   =  DEG2RAD*yh; 
      theta1 = -DEG2RAD*xw; 
      phi1   = -DEG2RAD*yh; 
    }
    else
    { 
      theta0= -atan((x-xw/2.0)/dist);              /* Angles to aim at target */
      phi0  = -atan((y-yh/2.0)/dist);
      theta1= -atan((x+xw/2.0)/dist);
      phi1  = -atan((y+yh/2.0)/dist);
    }

    /* shot towards target : flat distribution */

    if (gaussian)
    {
      theta= theta0+(theta1- theta0)*(randnorm()*FWHM2RMS+0.5);
      phi  = phi0  +(phi1  - phi0)  *(randnorm()*FWHM2RMS+0.5);
    }
    else
    {
      theta= theta0+(theta1- theta0)*rand01();
      phi  = phi0  +(phi1  - phi0)  *rand01();
    }
    /* Assume linear distribution */
    lambda = lambda_0+d_lambda*randpm1();
    if (lambda <= 0) ABSORB;
    lambda2 = lambda*lambda;
    lambda5 = lambda2*lambda2*lambda;

    v = K2V*(2*PI/lambda);

    p *= p_in*fabs((theta1- theta0)*(phi1  - phi0));

    if (T1 > 0)
    {
            Maxwell= L2P/lambda5*exp(-lambda02/lambda2);  /* 1/AA */
            if ((T2 > 0) && (I1 != 0))
            {
              Maxwell += (I2/I1)*L2Pb/lambda5*exp(-lambda02b/lambda2);
            }
            if ((T3 > 0) && (I1 != 0))
            {
              Maxwell += (I3/I1)*L2Pc/lambda5*exp(-lambda02c/lambda2);
            }
            p *= Maxwell;
    }

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

MCDISPLAY
%{
  double xmin;
  double xmax;
  double ymin;
  double ymax;

  if ((h == 0) || (w == 0))
  {
    magnify("xy");
    circle("xy",0,0,0,radius);
    if (gaussian)
      circle("xy",0,0,0,radius/2);
  }
  else
  {
    xmin = -w/2; xmax = w/2;
    ymin = -h/2; ymax = h/2;

    magnify("xy");
    multiline(5, (double)xmin, (double)ymin, 0.0,
             (double)xmax, (double)ymin, 0.0,
             (double)xmax, (double)ymax, 0.0,
             (double)xmin, (double)ymax, 0.0,
             (double)xmin, (double)ymin, 0.0);
    if (gaussian)
      circle("xy",0,0,0,sqrt(w*w+h*h)/4);
  }
%}

END