/*******************************************************************************
*
* Mcstas, neutron ray-tracing package
*         Copyright (C) 1997-2012, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: ESS_moderator
*
* %I
* Written by: P Willendrup and E Klinkby, February 2014, derived from K Lefmann ESS_moderator_long
* Modified by: 
* Version: $Revision: 4592 $
* Origin: DTU
* Release: McStas 2.1
*
* A parametrised pulsed source for modelling ESS long pulses.
*
* %D
* Produces a time-of-flight spectrum, from the ESS parameters
* Chooses evenly in lambda, evenly/exponentially decaying in time
* Adapted from ESS_moderator_long by: K Lefmann, 2001
*
* 2014-updates:
* <ol>
* <li>Selectable between simple flat single-moderator (planar=1) and MCNPX-like dual-moderator (thermal and cold) setup.
* <b>IMPORTANT</b>: In case of the dual-moderator geometry, the origin of the component is inside the cylindrical 
* moderator, i.e. take care when positioning the next components! (E Klinkby)
* <li>The spectrum from the source(s) is defined via the sourcedef string input parameter which allows these values:
* <ul>
* <li>sourcedef="2001", legacy "Mezei moderators" from the original F. Mezei documents 
* "ESS reference moderator characteristics for generic instrument performance evaluation", but rescaled to ESS TDR frequency, pulselength and power.
* <li>sourcedef="TDR", Mezei moderators, with a wavelength-dependent correction term to the cold flux, derived from 
* 2012 MCNPX calculations by ESS neutronics group. Corrections calculated by K Lieutenant (Vitess) and 
* implemented here by E Klinkby. NOTE: uses the 2001 brilliance for the thermal moderator!
* <li>sourcedef="2014", updated brilliance using formulation by Troels Schoenfeldt, including support for the "pancacke", i.e. flat geometry.
*</ul>
* <li>The component can use target_index for focusing to a given beam port. Use an Arm() and ROTATED to position 
* relatively to the moderator.
* <li>Time focusing option: Adjusts neutron departure time to match a 'first chopper' defined by parameters tfocus_dist, tfocus_time, tfocus_width (K Lefmann).
* NOTE: Currently only applies to the Mezei-moderator.
* <li>The component relies on the new ess_source-lib which is expected to become further enriched during design-finaliziation and construciton of the ESS.
* </ol>
*
* Units of flux: n/cm^2/s/AA/ster
* (McStas units are in general neutrons/second)
*
* Derived from ESS_moderator_long which was debugged intensively against Mezei note (4/12 2000) and VitESS @ Rencurel 2006.
*
* %VALIDATION 
* Validated against VitESS and Mezei note (4/12 2000) @ Rencurel 2006
*
* %P
* Input parameters:
*
* yheight_c: (m)    Height of the cylindershaped cold source
* yheight_t: (m)    Height of the rectangular thermal source
* cyl_radius:(m) Radius of the cylindershaped cold source - coupled to width_c
* width_t: (m)    Edge of cube shaped thermal source
* width_c: (m)    Arc-length opening of cylindershaped cold source - coupled to cyl_radius. (Used as width of rectangular geometry if planar=1)
* Lmin:   (AA)   Lower edge of wavelength distribution
* Lmax:   (AA)   Upper edge of wavelength distribution
* dist:   (m)    Distance from source to focusing rectangle; at (0,0,dist)
* focus_xw:(m)   Width of focusing rectangle
* focus_yh:(m)   Height of focusing rectangle
* target_index:(1)  relative index of component to focus at, e.g. next is +1
*                this is used to compute 'dist' automatically.
* n_pulses: (1)  Number of pulses simulated. 0 and 1 creates one pulse. 
*                The integrated intensity is constant 
* cold_frac: (1) Fraction of neutron statistics from cold source. It is implicitely assumed 
*                that supermirror allows each beamline to choose the desired fraction
*                of cold and thermal neutrons (i.e. extreme idealization).
* tfocus_dist: (m) Position of time window
* tfocus_time: (s) Time position of window
* tfocus_width: (s) Time width of window
* beamport_angle: (deg) Direction within the beamport sector (0 < angle < 60) to direct neutrons
* acc_power: (MW) Accelerator power in MW
* sourcedef: (string) ESS source "database", values: "TDR", "2001", "2013", "2014"
* planar: (boolean) If set to 1, source geometry is flat and rectangular. Can only be used if cold_frac is 1 or 0, i.e. for fully thermal or fully cold source. May overwrite other geometry input.
* isleft: (boolean) If set to 1, place thermal wing left of cold cylinder, if set to 0, place thermal wing right of cold cylinder.
*
* %E
*******************************************************************************/

DEFINE COMPONENT ESS_moderator
DEFINITION PARAMETERS ()
  SETTING PARAMETERS (cyl_radius=0.085, width_c=0, yheight_c=0.12, yheight_t=0.12, width_t=0.12, beamport_angle=30, planar=0, isleft=1,
		      Lmin, Lmax, cold_frac=1.0,
		      dist=0, focus_xw, focus_yh, int target_index=0, 
		      tfocus_dist=0.1, tfocus_time=0.0, tfocus_width=0.0,
		      int n_pulses=1, string sourcedef="TDR", acc_power=5)
  OUTPUT PARAMETERS (cold_bril, thermal_bril, l_range, w_mult, w_geom, w_geom, w_geom_t, yheight)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ 
SHARE
%{
  %include "ess_source-lib"
%}

DECLARE
%{
  double l_range, w_mult, w_geom, w_geom_c, w_geom_t, yheight;
  double tx,ty,tz;
  double t1x,t1y,t1z,t2x,t2y,t2z;
  
  /* Target station geometry... */
  double r_empty = 2.0; /* two meters from moderator surface and out... */
  double r_optics;
  ess_moderator_struct modextras;
  functype cold_bril;
  functype thermal_bril;
%}

INITIALIZE
%{
  n_pulses=(double)floor(n_pulses);
  if (n_pulses == 0) n_pulses=1;
 
  if (target_index && !dist)
  {
    Coords ToTarget;
    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);
  } else if (target_index && !dist) {
    printf("ESS_moderator: %s: Please choose to set either the dist parameter or specify a target_index.\nExit\n", NAME_CURRENT_COMP);
    exit(-1);
  } else {
    tx=0, ty=0, tz=dist;
  }

  if (focus_xw < 0 || focus_yh < 0)
  {
    printf("ESS_moderator: %s: Please specify both focus_xw and focus_yh as positive numbers.\nExit\n", NAME_CURRENT_COMP);
    exit(-1);
  }

  if (dist < r_empty && dist > 0)
  {
    printf("ESS_moderator: %s WARNING: Provided dist parameter is %g and hence inside the vacated zone of the beam extraction system!\nYou might be placing optics in a restricted area!!!\n", NAME_CURRENT_COMP, dist);
  }
    
  if (beamport_angle < 0 || beamport_angle > 60)
  {
    printf("ESS_moderator: %s: Please select a beamport_angle between 0 and 60 degrees!\nExit\n", NAME_CURRENT_COMP);
    exit(-1);
  }
  
  if (width_c && cyl_radius) {
    printf("ESS_moderator: %s: Please specify EITHER cold-moderator radius (cyl_radius) or length of visible arch (width_c)!\nExit\n", NAME_CURRENT_COMP);
    exit(-1);
  } else if (cyl_radius) {
    width_c = 2*PI*cyl_radius*60/360;
  } else {
    cyl_radius = 360*width_c/(2*PI*60);
  }
  r_optics = 6.0 - r_empty - cyl_radius;

  if (Lmin<=0 || Lmax <=0 || dist == 0)
  {
    printf("ESS_moderator: %s: Check parameters (lead to Math Error).\n Avoid 0 value for {Lmin Lmax dist d tau branch1/2/tail} and 1 value for {n n2 branch1/2/tail}\n", NAME_CURRENT_COMP);
    exit(-1);
  }

  if (planar==1 && !(cold_frac==1 || cold_frac==0))
  {
    printf("ESS_moderator: %s: Planar==1 can only be used with cold_frac==1 or cold_frac==0, i.e. purely cold or thermal moderator\n", NAME_CURRENT_COMP);
    exit(-1);
  }
  
  l_range = Lmax-Lmin;
  w_geom_c  = width_c*yheight_c*1.0e4;     /* source area correction */
  w_geom_t  = width_t*yheight_t*1.0e4;
  w_mult  = l_range;            /* wavelength range correction */
  w_mult *= 1.0/mcget_ncount();   /* Correct for number of rays */

  if (planar==1) {
    if (cold_frac==0) {
      /* Flat thermal source, horz basis vector(s) along x */
      t1x=1;
      t1y=0;
      t1z=0;
      t2x=1;
      t2y=0;
      t2z=0;      
    } /* Flat cold source is handled below */
  } else if (planar==0){
    /* Calculate location of thermal wings wrt beamport_angle (z) direction */
    /* Wing 1 (left) is at -beamport_angle */
    t1z = cyl_radius*cos(-DEG2RAD*beamport_angle);
    t1x = cyl_radius*sin(-DEG2RAD*beamport_angle);
    t1y = 0;
    /* Wing 2 (right) is at 60-beamport_angle */
    t2z = cyl_radius*cos(DEG2RAD*(60-beamport_angle));
    t2x = cyl_radius*sin(DEG2RAD*(60-beamport_angle));
    t2y = 0;
  /* We want unit vectors... */
    NORM(t1x,t1y,t1z);
    NORM(t2x,t2y,t2z);
  } else {
    printf("ESS_moderator: %s: Planar is a boolean input parameter, please give either 0 or 1\n", NAME_CURRENT_COMP);
    exit(-1);
  }
  
  if (strcasestr(sourcedef,"TDR")) {
    w_mult *= ESS_SOURCE_FREQUENCY;               /* Correct for frequency */
    w_mult *= acc_power/5;                        /* Correct for accelerator power */
    printf("Using ESS TDR brilliance\n");
    cold_bril=ESS_2012_Lieutenant_cold;
    thermal_bril=ESS_Mezei_thermal;
  } else if (strcasestr(sourcedef,"2001")) {
    w_mult *= ESS_SOURCE_FREQUENCY;               /* Correct for frequency */
    w_mult *= acc_power/5;                        /* Correct for accelerator power */
    printf("Using ESS 2001 brilliance\n");
    cold_bril=ESS_Mezei_cold;
    thermal_bril=ESS_Mezei_thermal;
  } else if (strcasestr(sourcedef,"2014")) {
    // w_mult *= ESS_SOURCE_FREQUENCY;               /* Is already in per-second units */
    w_mult *= acc_power;                             /* Is already in per-MW units */
    printf("Using ESS 2014 brilliance\n");
    // Fix: transfer in m units rather than cm's...
    modextras.height_c=100.0*yheight_c;
    modextras.height_t=100.0*yheight_t;
    cold_bril=ESS_2013_Schoenfeldt_cold;
    thermal_bril=ESS_2013_Schoenfeldt_thermal;
  } else {
    printf("Soure def %s not implemented!\n",sourcedef);
    exit(-1);
  }

  /* For drawing */ 
  yheight = fmax(yheight_t, yheight_c);

%}
TRACE
%{
  double v,tau_l,E,lambda,k,r,xf,yf,dx,dy,w_focus,tail_flag,cor,dt,xprime,yprime,zprime;

  /* Bispectral source - choice of spectrum and initial position */
  int cold = ( rand01() < cold_frac ) ? 1 : 0;

  /* Geometry adapted from ESS MCNPX model, mid 2012 */
  if (cold) {          //case: cold moderator
    if (planar==1) {
      x = 0.5*randpm1()*width_c;
      y = 0.5*randpm1()*yheight_c;
      z = 0;
    } else {
      double theta_tmp;
      
      //choose random point on cylinder surface
      theta_tmp = randpm1()*PI/6 + PI/2 + (- 30 + beamport_angle)*DEG2RAD;
      x     = cyl_radius * cos(theta_tmp);
      y     = 0.5*randpm1()*yheight;
      z     = cyl_radius * sin(theta_tmp);
    }
    
    w_geom = w_geom_c;
  }  else  {                      //case: thermal moderator
    double poshorz, posvert;
    
    if (planar==1) {
      x = 0.5*randpm1()*width_t;
      y = 0.5*randpm1()*yheight_t;
      z = 0;
    } else {
      poshorz = cyl_radius+rand01()*width_t;
      posvert = 0.5*randpm1()*yheight;
      
      if (isleft) {
	x = t1x * poshorz;
	z = t1z * poshorz;
      } else {
	x = t2x * poshorz;
      z = t2z * poshorz;
      }
      y = posvert;
    }
    w_geom = w_geom_t;
  }

  randvec_target_rect_real(&xf, &yf, &r, &w_focus,
			   tx, ty, tz, 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;

  /* Determine delta-t needed to reach first chopper */
  if (tfocus_width>0) {
    dt = tfocus_dist/vz;			/* Flight time to time window (chopper) */
  }

  if (cold) {          //case: cold moderator
    cold_bril( &t,  &p,  lambda,  tfocus_width,  tfocus_time,  dt, modextras);
  }  else  {                      //case: thermal moderator
    thermal_bril( &t,  &p,  lambda,  tfocus_width,  tfocus_time,  dt, modextras);
  }
  p*=w_focus*w_geom*w_mult;
  t+=(double)floor((n_pulses)*rand01())/ESS_SOURCE_FREQUENCY;   /* Select a random pulse */

  /* Correct weight for sampling of cold vs. thermal events. */
  if (cold) {
    p /=cold_frac;
  } else {
    p/=(1-cold_frac);
  }
  
%}

MCDISPLAY
%{
  if (planar==0) {
    /* Draw cold moderator as cylinder */
    
    circle("xz", 0,  yheight_c/2.0, 0, cyl_radius);
    circle("xz", 0,  -yheight_c/2.0, 0, cyl_radius);
    line(0, -yheight_c/2.0, cyl_radius, 0, yheight_c/2.0, cyl_radius);
    line(0, -yheight_c/2.0, -cyl_radius, 0, yheight_c/2.0, -cyl_radius);
    line(cyl_radius, yheight_c/2.0, 0, cyl_radius, yheight_c/2.0, 0);  
    line(-cyl_radius, -yheight_c/2.0, 0, -cyl_radius, yheight_c/2.0, 0);
    /* Draw thermal moderators as a couple of squares + some lines */
    // Left
    multiline(4, t1x*cyl_radius, -yheight_t/2.0, t1z*cyl_radius,
	      t1x*(cyl_radius + width_t), -yheight_t/2.0, t1z*(cyl_radius + width_t),
	      t1x*(cyl_radius + width_t), yheight_t/2.0,  t1z*(cyl_radius + width_t),
	      t1x*cyl_radius, yheight_t/2.0, t1z*cyl_radius);
    // Right
    multiline(4, t2x*cyl_radius, -yheight_t/2.0, t2z*cyl_radius,
	      t2x*(cyl_radius + width_t), -yheight_t/2.0, t2z*(cyl_radius + width_t),
	      t2x*(cyl_radius + width_t), yheight_t/2.0,  t2z*(cyl_radius + width_t),
	      t2x*cyl_radius, yheight_t/2.0, t2z*cyl_radius);
    
    /* Dashed lines for indicating "beam extraction" area... */
    dashed_line(t1x*cyl_radius, -yheight/2.0, t1z*cyl_radius, t1x*r_empty, -yheight/2.0, t1z*r_empty,10);
    dashed_line(t1x*cyl_radius, yheight/2.0, t1z*cyl_radius, t1x*r_empty, yheight/2.0, t1z*r_empty,10);
    dashed_line(t2x*cyl_radius, -yheight/2.0, t2z*cyl_radius, t2x*r_empty, -yheight/2.0, t2z*r_empty,5);
    dashed_line(t2x*cyl_radius, yheight/2.0, t2z*cyl_radius, t2x*r_empty, yheight/2.0, t2z*r_empty,5);
    
    /* Circles indicating extent of the "empty" zone where optics is not allowed */
    circle("xz", 0,  yheight/2.0, 0, r_empty);
    circle("xz", 0,  -yheight/2.0, 0, r_empty);
    
    /* Circles indicating the builk shielding of the target monolith at 6 m */
    circle("xz", 0,  focus_yh/2.0 , 0, 6);
    circle("xz", 0, -focus_yh/2.0 , 0, 6);
    circle("xz", 0,  2, 0, 6);
    circle("xz", 0, -2, 0, 6);
    
    /* Rectangle indicating the chosen focus rectangle - where the optics starts... */
    rectangle("xy",tx,ty,tz,focus_xw,focus_yh);
  } else {
    /* Simple 'flat', single-spectrum geometry */
    double tmpheight,tmpwidth;
    if (cold_frac==0) {
      tmpheight = yheight_t;
      tmpwidth = width_t;
    } else {
      tmpheight = yheight_c;
      tmpwidth = width_c;
    }
    rectangle("xy",0,0,0,tmpwidth,tmpheight);
    rectangle("xy",tx,ty,tz,focus_xw,focus_yh);
  }  
%}

END