/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Source_Maxwell_3
*
* %I
* Written by: Kim Lefmann
* Date: March 2001
* Origin: Risoe
*
* Source with up to three Maxwellian distributions
*
* %D
* A parametrised continuous source for modelling a (cubic) source
* with (up to) 3 Maxwellian distributions.
* The source produces a continuous spectrum.
* The sampling of the neutrons is uniform in wavelength.
*
* Units of flux: neutrons/cm^2/second/ster
* (McStas units are in general neutrons/second)
*
* Example:  PSI cold source T1=150.42 K / 2.51 AA     I1 = 3.67 E11
*                           T2=38.74 K / 4.95 AA      I2 = 3.64 E11
*                           T3=14.84 K / 9.5 AA       I3 = 0.95 E11
*
* %P
* Input parameters:
*
* yheight: [m]        Height of rectangular source
* xwidth: [m]         Width of rectangular source
* Lmin: [AA]          Lower edge of lambda distribution
* Lmax: [AA]          Upper edge of lambda distribution
* lambda0: [AA]       Mean wavelength of neutrons.
* dlambda: [AA]       Wavelength spread of neutrons.
* target_index: [1]   relative index of component to focus at, e.g. next is +1 this is used to compute 'dist' automatically.
* focus_xw: [m]       Width of focusing rectangle
* focus_yh: [m]       Height of focusing rectangle
* T1: [K]             1st temperature of thermal distribution
* I1: [1/(cm**2*st)]  flux, 1 (in flux units, see above)
*
* Optional parameters:
* dist: [m]           Distance from source to focusing rectangle; at (0,0,dist)
* T2: [K]             2nd temperature of thermal distribution
* T3: [K]             3nd temperature of  - - -
* I2: [1/(cm**2*st)]  flux, 2 (in flux units, see above)
* I3: [1/(cm**2*st)]  flux, 3  - - -
* size: [m]           Edge of cube shaped source (for backward compatibility)
*
* %E
*******************************************************************************/

DEFINE COMPONENT Source_Maxwell_3



SETTING PARAMETERS (size=0, yheight=0, xwidth=0, Lmin, Lmax, dist,
focus_xw, focus_yh,
T1, T2=300, T3=300, I1, I2=0, I3=0,
int target_index=1,lambda0=0, dlambda=0)


/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

SHARE
%{
/* A normalised Maxwellian distribution : Integral over all l = 1 */
#pragma acc routine seq
double SM3_Maxwell(double l, double temp)
  {
    double a=949.0/temp;
    return 2*a*a*exp(-a/(l*l))/(l*l*l*l*l);
  }
%}

DECLARE
%{
  double l_range;
  double w_mult;
  double w_source;
  double h_source;
%}

INITIALIZE
%{
  if (target_index && !dist)
  {
    Coords ToTarget;
    double tx,ty,tz;
    ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),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);
  }

  if (size>0) {
    w_source = h_source = size;
  } else {
    w_source = xwidth;
    h_source = yheight;
  }
  if (lambda0) {
    Lmin=lambda0-dlambda;
    Lmax=lambda0+dlambda;
  }
  l_range = Lmax-Lmin;
  w_mult = w_source*h_source*1.0e4;     /* source area correction */
  w_mult *= l_range;            /* wavelength range correction */
  w_mult *= 1.0/mcget_ncount();   /* correct for # neutron rays */

  if (w_source <0 || h_source < 0 || Lmin <= 0 || Lmax <= 0 || dist <= 0 || T1 <= 0 || T2 <= 0|| T3 <= 0 || Lmax<=Lmin) {
      printf("Source_Maxwell_3: %s: Error in input parameter values!\n"
             "ERROR          Exiting\n",
           NAME_CURRENT_COMP);
      exit(0);
  }

%}

TRACE
%{
  double v,tau_l,E,lambda,k,r,xf,yf,dx,dy,w_focus;

  t=0;
  z=0;
  x = 0.5*w_source*randpm1();
  y = 0.5*h_source*randpm1();         /* Choose initial position */

  randvec_target_rect_real(&xf, &yf, &r, &w_focus,
		      0, 0, dist, focus_xw, focus_yh, ROT_A_CURRENT_COMP, x, y, z, 2);

  dx = xf-x;
  dy = yf-y;
  r = sqrt(dx*dx+dy*dy+dist*dist);

  lambda = Lmin+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);  */

  p *= w_mult*w_focus;                /* Correct for target focusing etc */
  p *= I1*SM3_Maxwell(lambda,T1)+I2*SM3_Maxwell(lambda,T2)+I3*SM3_Maxwell(lambda,T3);
                                        /* Calculate true intensity */
%}

MCDISPLAY
%{
  
  multiline(5, -(double)focus_xw/2.0, -(double)focus_yh/2.0, 0.0,
                (double)focus_xw/2.0, -(double)focus_yh/2.0, 0.0,
                (double)focus_xw/2.0,  (double)focus_yh/2.0, 0.0,
               -(double)focus_xw/2.0,  (double)focus_yh/2.0, 0.0,
               -(double)focus_xw/2.0, -(double)focus_yh/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