/*******************************************************************************
*
* 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: J. Saroun (saroun@ujf.cas.cz)
* Origin: DTU
*
* 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
*
* Updates and simplified interface:
* <ol>
* <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.
* <li>sourcedef="2015", updated brilliance using formulation by Troels Schoenfeldt, new butterfly baseline.
*</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>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>
*
* <p><b>Note that this component does not implement "engineering reality" and currently uses a coordinate system centered on the moderator assembly. An
* updated moderator component which references the "Moderator focus coordinate system" will be released later during the spring of 2016.</b>
*
* <p>Derived from ESS_moderator_long which was debugged intensively against Mezei note (4/12 2000) and VitESS @ Rencurel 2006.
*
*-----------------------------------------------
* Correction by J. Saroun, NPI Rez:
* 1) version 2015: accepts negative port angles
* 2) version 2015: weight by cosine of the port angle
* Warning: The negative beamport angle is not taken into acccount by mcplot
*
* %VALIDATION
* Mezei-modererators validated against VitESS and Mezei note (4/12 2000) @ Rencurel 2006
* Benchmarked against multiple versions of ESS moderator group simulation data 2013-2015
*
* %P
* Input parameters:
*
* yheight_c: [m]             Height of the cold source
* yheight_t: [m]             Height of the thermal source
* xwidth_t: [m]              Edge of thermal source
* xwidth_c: [m]              Width / arc-length opening of the cold source.
* 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).
* tmax_multiplier: [1]       Defined maximum emission time at moderator, tmax= tmax_multiplier * ESS_PULSE_DURATION. Only in combination with sourcedef="2013", "2014" or "2015"
* beamport_angle: [deg]      Direction within the beamport sector (0 < angle < extraction_opening for 2014, -extraction_opening/2 < angle < extraction_opening/2 for 2015) to direct neutrons. For sourcedef="2015", the only allowed values are 5,15,...,55 degrees measured from the central point.
* extraction_opening: [deg]  Width of extraction-area in degrees (60 or 120 degrees). 120 deg only in combination with sourcedef="2014" and "2015".
* acc_power: [MW]            Accelerator power in MW
* sourcedef: [string]        ESS source "database", values: "TDR", "2001", "2013", "2014", "2015"
* isleft: [1]                Fraction of thermal neutrons generated at the "left" moderator slab in case of "2013" or "2014"
*
* %E
*******************************************************************************/

DEFINE COMPONENT ESS_moderator

SETTING PARAMETERS (isleft=0.9,
Lmin, Lmax, cold_frac=1.0, dist=0, focus_xw, focus_yh, int target_index=0, tmax_multiplier=3,yheight_c=0.12, yheight_t=0.12,
int n_pulses=1, acc_power=5, beamport_angle=-1, string sourcedef="2015", xwidth_c=0.1, xwidth_t=0.18, extraction_opening=120)

SHARE
%{
%include "ess_source-lib"
%}

DECLARE
%{
// Weighting-oriented scalars
double l_range, w_mult, w_geom, w_geom_c, w_geom_t, w_stat;
// Flag to indicate if originating from cold or thermal moderator
int cold;
// Flag to indicate if originating from left or right thermal slab
int wasleft;
// Flag to indicate if we should apply the solid angle correction from randvec
int cosine;
// Target coordinates (focusing)
  double tx,ty,tz;
  // Displacement of cold moderator centre from comp origin, cold moderator radius 2014 only!
  double cx,cy,cz,cyl_radius;
  // Vectors giving direction of the thermal slabs
  double t1x,t1y,t1z,t2x,t2y,t2z;

  /* Struct for transfer to ess-source-lib funcs  */
  ess_moderator_struct modextras;
  /* Cold and thermal function pointers */
  functype cold_bril;
  functype thermal_bril;
  /* Time-focusing parameters - currently unused */
  double tfocus_width,  tfocus_time,  dt;
  /* radius of vacated zone around the moderator */
  double r_empty;
  double lambda,xprime,yprime,zprime;
  double abs_beamport_angle,cos_beamport_angle,sin_beamport_angle;
  /* variables needed to correct for the emission surface angle*/
  double cos_thermal,cos_cold, edge_thermal;
  int sgnport;
%}

INITIALIZE
%{

  if (Lmin>=Lmax) {
    printf("ESS_moderator: %s: Unmeaningful definition of wavelength range!\n ERROR - Exiting\n",
           NAME_CURRENT_COMP);
      exit(0);
  }

  sgnport=(beamport_angle>0 ? 1:-1);
  if (sgnport<0) {
	printf("WARNING, %s: beamport angle < 0, experimental option\n", NAME_CURRENT_COMP);
  }

/* ESS 2015 - variables needed to correct for the emission surface angle */
  abs_beamport_angle=fabs(beamport_angle);
  cos_beamport_angle=cos(beamport_angle*DEG2RAD);
  sin_beamport_angle=sin(beamport_angle*DEG2RAD);
  /* correction for projection along the beam / projection on the z=0 plane */
  cos_thermal=cos_beamport_angle;
  cos_cold=cos((abs_beamport_angle-24.24)*DEG2RAD)/cos(24.24*DEG2RAD);
  printf("input to cosine %g, output %g\n",abs_beamport_angle, cos_cold);
  /* cross-section between the z=const and tilted surfaces [m] */
  edge_thermal=-0.0716;

  r_empty=2;

  /* Weight multipliers */
  w_mult=acc_power/5;
  w_geom=1;
  w_geom_c=1;
  w_geom_t=1;
  w_stat=1.0/mcget_ncount();

  tfocus_width=0;  tfocus_time=0;  dt=0;

  n_pulses=(double)floor(n_pulses);
  if (n_pulses == 0) n_pulses=1;

  /* Calculation of location of focusing rectangle */
  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);
	printf("%s: Focusing on a window centered at [%g, %g, %g]\n", NAME_CURRENT_COMP,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 (Lmin<=0 || Lmax <=0 || dist == 0)
  {
    printf("ESS_moderator: %s: Check parameters (lead to Math Error).\n Avoid 0 value for {Lmin Lmax dist} \n", NAME_CURRENT_COMP);
    exit(-1);
  }

  /* Calculate cold cylinder radius from extraction_opening */
  cyl_radius = 360*xwidth_c/(2*PI*extraction_opening);

  /* Is abs_beamport_angle set at a meaningful angle - otherwise recalc */
  if ( abs_beamport_angle<0 || abs_beamport_angle>extraction_opening){
    beamport_angle = extraction_opening/2.0;
	abs_beamport_angle=beamport_angle;
    printf("%s: WARNING: Resetting your beamport_angle to %g (%g/2) - central beamslot\n",NAME_CURRENT_COMP,beamport_angle,extraction_opening);
  }

  /* Calculate positioning of thermal slabs and cold cylinder - NOTE: overwritten in case of ESS 2014 */
  t1z = cyl_radius*cos(-DEG2RAD*abs_beamport_angle);
  t1x = cyl_radius*sin(-DEG2RAD*abs_beamport_angle);
  t1y = 0;
  /* Wing 2 (right) is at extraction_width-beamport_angle */
  t2z = cyl_radius*cos(DEG2RAD*(extraction_opening-abs_beamport_angle));
  t2x = cyl_radius*sin(DEG2RAD*(extraction_opening-abs_beamport_angle));
  t2y = 0;
  /* We want unit vectors... */
  NORM(t1x,t1y,t1z);
  NORM(t2x,t2y,t2z);
  cx = 0; cy=0; cz=0;

  /* Geometry parameters and brilliances specified via the sourcedef string */
  /* ESS-TDR definition: */
  if (strcasestr(sourcedef,"TDR")) {
    w_mult *= ESS_SOURCE_FREQUENCY;               /* Correct for frequency */
    printf("Using ESS TDR brilliance\n");
    cold_bril=ESS_2012_Lieutenant_cold;
    thermal_bril=ESS_Mezei_thermal;
 /* ESS-2001 definition: */
  } else if (strcasestr(sourcedef,"2001")) {
    w_mult *= ESS_SOURCE_FREQUENCY;               /* Correct for frequency */
    printf("Using ESS 2001 brilliance\n");
    cold_bril=ESS_Mezei_cold;
    thermal_bril=ESS_Mezei_thermal;
 /* ESS-2013 post-TDR updated definition, before ESS "pancake": */
  } else if (strcasestr(sourcedef,"2013")) {
    w_mult *= acc_power;                             /* Is already in per-MW units */
    printf("Using ESS 2013 brilliance\n");
    modextras.height_c=yheight_c;
    modextras.height_t=yheight_t;
    modextras.tmultiplier=tmax_multiplier;
    cold_bril=ESS_2013_Schoenfeldt_cold;
    thermal_bril=ESS_2013_Schoenfeldt_thermal;
  /* ESS-2014 definition: */
  } else if (strcasestr(sourcedef,"2014")) {
    if (xwidth_c!=0.23 || xwidth_t!=0.12 || extraction_opening!=120) {
      fprintf(stderr,"FATAL: sourcedef %s only allows xwidth_c=0.23, xwidth_t=0.12 and extraction_opening=120 !\n",sourcedef);
      exit(-1);
    }
    /* Specify brilliance fct.'s */
    cold_bril=ESS_2014_Schoenfeldt_cold;
    thermal_bril=ESS_2014_Schoenfeldt_thermal;
    /* Cold moderator geometry */
    /* In this mode of operation, the emission is from the x-y plane*/
    cx = 0; cy=0; cz=-0.12; cyl_radius=0.12;
    /* 120 degree beam extraction opening, defining thermal moderator directions */
    t1x=sin(DEG2RAD*(120/2-abs_beamport_angle));t1z=cos(DEG2RAD*120/2);
    t2x=-sin(DEG2RAD*120/2);t2z=cos(DEG2RAD*120/2);
    t1y=0; t2y=0;
    /* ESS-2015 definition: */
   } else if (strcasestr(sourcedef,"2015")) {
    if (xwidth_c!=0.1 || xwidth_t!=0.18 || extraction_opening!=120) {
      fprintf(stderr,"FATAL: sourcedef %s only allows xwidth_c=0.1, xwidth_t=0.18 and extraction_opening=120 !\n",sourcedef);
      exit(-1);
    }
    /* Specify brilliance fct.'s */
    cold_bril=ESS_2015_Schoenfeldt_cold;
    thermal_bril=ESS_2015_Schoenfeldt_thermal;
    /* Moderator geometry */
    cx = 0; cy=0; cz=0; cyl_radius=0;
    /* 120 degree beam extraction opening */
    t1x=sin(DEG2RAD*(120/2+abs_beamport_angle));t1z=cos(DEG2RAD*(120/2+abs_beamport_angle));
    t2x=-sin(DEG2RAD*(120/2-abs_beamport_angle));t2z=cos(DEG2RAD*(120/2-abs_beamport_angle));
    t1y=0; t2y=0;
   } else {
    fprintf(stderr,"FATAL: sourcedef %s is not defined!\n",sourcedef);
    exit(-1);
  }
  modextras.height_c=yheight_c;
  modextras.Width_c=xwidth_c;
  modextras.Width_t=xwidth_t;
  modextras.height_t=yheight_t;
  modextras.tmultiplier=tmax_multiplier;
  modextras.extractionangle=extraction_opening;

  if ((strcasestr(sourcedef,"2014")) && abs_beamport_angle != modextras.extractionangle/2.0) {
    printf("%s: WARNING: beamport_angle is not at central slot. With sourcedef %s this means application of a partially analytical brightness model for B(\\theta) \n",NAME_CURRENT_COMP,sourcedef);
  }
  if ((strcasestr(sourcedef,"2015")) && (abs_beamport_angle!=5.0 && abs_beamport_angle!=15.0 && abs_beamport_angle!=25.0 && abs_beamport_angle!=35.0 && abs_beamport_angle!=45.0 && abs_beamport_angle!=55.0)) {
    printf("%s: FATAL: beamport_angle is not at a defined value. Allowed are: (5, 15, ...) and < 60  \n",NAME_CURRENT_COMP);
    exit(-1);
  }
  modextras.beamportangle=abs_beamport_angle;

  l_range = Lmax-Lmin;
  w_geom_c  = xwidth_c*yheight_c*1.0e4;     /* source area correction */
  w_geom_t  = xwidth_t*yheight_t*1.0e4;
  w_mult *= l_range;            /* wavelength range correction */

 %}
TRACE
%{
  p=1;
  double v,E,k,r,xf,yf,dx,dy,dz,w_focus,tail_flag,cor;

  /* Bispectral source - choice of spectrum and initial position */
  cold = ( rand01() < cold_frac );

  /* Use standard McStas routine or -not for specifying emission directionality */
  cosine = 0;

  /* Emission geometry adapted from ESS MCNPX model, mid 2012 */
  if (strcasestr(sourcedef,"2014")) {
    if (cold) { //case: cold moderator
      x = 0.5*randpm1()*xwidth_c;
      y = 0.5*randpm1()*yheight_c;
      z = 0;
      cosine = 0;
      w_geom = w_geom_c;
    }  else  { //case: thermal moderator
      y = 0.5*randpm1()*yheight_t;
      if (rand01()<isleft) {
	wasleft=1;
	x = -xwidth_c/2.0 - rand01()*xwidth_t;
      } else {
	wasleft=0;
	x = xwidth_c/2.0 + rand01()*xwidth_t;
      }
      z = 0;
      cosine = 0;
      w_geom = w_geom_t;
    }
  } else if (strcasestr(sourcedef,"2015")) {
    if (cold) { //case: cold moderator
      xprime = -(0.06 + 0.1*rand01());
      yprime = 0.5*randpm1()*yheight_c;
      zprime = TSC2015_z0_BF3cm(100*xprime)/100.0;
      cosine = 0;
      w_geom = w_geom_c;
    }  else  { //case: thermal moderator
      xprime = -(0.09*randpm1());
      yprime = 0.5*randpm1()*yheight_t;
      zprime = TSC2015_z0_BF3cm(100*xprime)/100.0;
      cosine = 0;
      w_geom = w_geom_t;
    }
    y = yprime;
    x = sgnport*xprime*cos_beamport_angle + zprime*sin_beamport_angle;
    z = cos_beamport_angle*zprime - sgnport*xprime*sin_beamport_angle;
	/* the position was projected in the beamport direction - we must multiply by cosine*/
	/* if (xprime>edge_thermal) { */
	/* 	w_geom *=cos_thermal; */
	/* } else { */
	/* 	w_geom *=cos_cold; */
	/* } */


  } else {
    if (cold) {          //case: cold moderator
      //choose random point on cylinder surface
      double theta_tmp = (rand01()*extraction_opening - abs_beamport_angle)*DEG2RAD;
      x     = cyl_radius * sin(theta_tmp);
      y     = 0.5*randpm1()*yheight_c;
      z     = cyl_radius * cos(theta_tmp);
      cosine = 0;
      w_geom = w_geom_c;
    }  else  {                      //case: thermal moderator
      double poshorz, posvert;
      poshorz = cyl_radius+rand01()*xwidth_t;
      posvert = 0.5*randpm1()*yheight_t;

      if (rand01()<isleft) {
	wasleft=1;
	x = t1x * poshorz;
	z = t1z * poshorz;
      } else {
	wasleft=0;
	x = t2x * poshorz;
	z = t2z * poshorz;
      }
      y = posvert;
      cosine = 1;
      w_geom = w_geom_t;
    }
    // Correct for comp origin
    x = x-cx;
    y = y-cy;
    z = z-cz;
  }

  randvec_target_rect_real(&xf, &yf, &r, &w_focus,
			   tx, ty, tz, focus_xw, focus_yh, ROT_A_CURRENT_COMP, x, y, z, cosine);

  /* In case of 2014 source, solid angle correction is done internally  in the brill fcts */
  if (strcasestr(sourcedef,"2014") || strcasestr(sourcedef,"2015")) {
    w_focus=focus_xw*focus_yh/(tx*tx+ty*ty+tz*tz);
  }

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

  lambda = Lmin+l_range*rand01();    /* Choose from uniform distribution */

  k = 2*PI/lambda;
  v = K2V*k;

  vz = v*dz/r;
  vy = v*dy/r;
  vx = v*dx/r;

  if (strcasestr(sourcedef,"2015")) {
    modextras.X=xprime; modextras.Y=yprime; modextras.Z=zprime; //modextras.Wasleft=1;
  } else {
    modextras.X=x; modextras.Y=y; modextras.Z=z; modextras.Wasleft=wasleft;
  }

  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_stat*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);
    // In the 2015 case we only describe one thermal emission area
    if (!strcasestr(sourcedef,"2015")) {
      if (wasleft) {
	p/=(isleft);
      } else {
	p/=(1-isleft);
      }
    }
  }
  SCATTER;

%}

MCDISPLAY
%{
  /* if (planar==0) { */
  /* Draw cold moderator as cylinder */
  
  if (!strcasestr(sourcedef,"2015")) {
    circle("xz", cx,  cy+yheight_c/2.0, cz, cyl_radius);
    circle("xz", cx,  cy-yheight_c/2.0, cz, cyl_radius);
    line(cx, cy-yheight_c/2.0, cz+cyl_radius, cx, cy+yheight_c/2.0, cz+cyl_radius);
    line(cx, cy-yheight_c/2.0, cz-cyl_radius, cx, cy+yheight_c/2.0, cz-cyl_radius);
    line(cx+cyl_radius, cy-yheight_c/2.0, cz, cx+cyl_radius, cy+yheight_c/2.0, cz);
    line(cx-cyl_radius, cy-yheight_c/2.0, cz, cx-cyl_radius, cy+yheight_c/2.0, cz);

    /* Draw thermal moderators as a couple of squares + some lines */

    // Left
    multiline(5, cx+t1x*cyl_radius, cy-yheight_t/2.0, cz+t1z*cyl_radius,
  	      cx+t1x*cyl_radius + t1x*xwidth_t, cy-yheight_t/2.0, cz+t1z*cyl_radius + t1z*xwidth_t,
  	      cx+t1x*cyl_radius + t1x*xwidth_t, cy+yheight_t/2.0, cz+t1z*cyl_radius + t1z*xwidth_t,
  	      cx+t1x*cyl_radius, cy+yheight_t/2.0, cz+t1z*cyl_radius, cx+t1x*cyl_radius, cy-yheight_t/2.0, cz+t1z*cyl_radius);
    // Right
    multiline(5, cx+t2x*cyl_radius, cy-yheight_t/2.0, cz+t2z*cyl_radius,
	      cx+t2x*cyl_radius + t2x*xwidth_t, cy-yheight_t/2.0, cz+t2z*cyl_radius + t2z*xwidth_t,
	      cx+t2x*cyl_radius + t2x*xwidth_t, cy+yheight_t/2.0, cz+t2z*cyl_radius + t2z*xwidth_t,
	      cx+t2x*cyl_radius, cy+yheight_t/2.0, cz+t2z*cyl_radius, cx+t2x*cyl_radius, cy-yheight_t/2.0, cz+t2z*cyl_radius);

    /* Dashed lines for indicating "beam extraction" area... */
    dashed_line(cx+t1x*cyl_radius,cy-yheight_c/2.0, cz+t1z*cyl_radius, cx+t1x*r_empty, cy-yheight_c/2.0, cz+t1z*r_empty,10);
    dashed_line(cx+t1x*cyl_radius, cy+yheight_c/2.0, cz+t1z*cyl_radius, cx+t1x*r_empty, cy+yheight_c/2.0, cz+t1z*r_empty,10);
    dashed_line(cx+t2x*cyl_radius, cy-yheight_c/2.0, cz+t2z*cyl_radius, cx+t2x*r_empty, cy-yheight_c/2.0, cz+t2z*r_empty,5);
    dashed_line(cx+t2x*cyl_radius, cy+yheight_c/2.0, cz+t2z*cyl_radius, cx+t2x*r_empty, cy+yheight_c/2.0, cz+t2z*r_empty,5);


    /* Circles indicating extent of the "empty" zone where optics is not allowed */
    circle("xz", cx,  cy+yheight_c/2.0, cz, r_empty);
    circle("xz", cx,  cy-yheight_c/2.0, cz, r_empty);

    /* Circles indicating the builk shielding of the target monolith at 6 m */
    circle("xz", cx,  cy+focus_yh/2.0 , cz, 6);
    circle("xz", cx, cy-focus_yh/2.0 , cz, 6);
    circle("xz", cx,  cy+2, cz, 6);
    circle("xz", cx, cy-2, cz, 6);
  } else { // 2015
    double dx = (0.09 + 0.16)/25;
    double xp0, yp0, zp0, xp1, yp1, zp1;
    double xr0, yr0, zr0, xr1, yr1, zr1;
    double dxp0, dzp0, dxr0, dzr0, dxr1, dzr1;
    double ax,az,bbx,bbz,ccx,ccz;
    /* Drawing the "emission planes" of the cold and thermal moderators */
    xp0 = -0.16;
    yp0 = -yheight_c/2.0;
    int k,j;
    for (k=0; k<2; k++) {
      for (j=0; j<26; j++) {
	xp1 = xp0 + dx;
	yp1 = yp0;
	zp0 = TSC2015_z0_BF3cm(100*xp0)/100.0;
	zp1 = TSC2015_z0_BF3cm(100*xp1)/100.0;
	xr0 = cos(-abs_beamport_angle*DEG2RAD)*xp0 - sin(-abs_beamport_angle*DEG2RAD)*zp0;
	zr0 = cos(-abs_beamport_angle*DEG2RAD)*zp0 + sin(-abs_beamport_angle*DEG2RAD)*xp0;
	xr1 = cos(-abs_beamport_angle*DEG2RAD)*xp1 - sin(-abs_beamport_angle*DEG2RAD)*zp1;
	zr1 = cos(-abs_beamport_angle*DEG2RAD)*zp1 + sin(-abs_beamport_angle*DEG2RAD)*xp1;
	if (j==0 || j==8 || j==11) {
	  line(xr0, -yheight_c/2.0, zr0, xr0, yheight_c/2.0, zr0);
	}
	if (j==25) {
	  line(xr1, -yheight_c/2.0, zr1, xr1, yheight_c/2.0, zr1);
	}
	if (j%2==0) { // Every other element skipped to make dashed lines...
	  line(xr0, yp0, zr0, xr1, yp1, zr1);
	}
	xp0 = xp1;
      }
      /* Drawing a sketch-version of the butterfly "cross" and wing(s).... */
      /* Thermal cross: */
      xp0 = 0.0827; xp1 = -xp0;
      zp0 = TSC2015_z0_BF3cm(100*xp0)/100.0;
      zp1 = zp0;
      dxp0 = 0.01; dzp0=0.01;
      xr0 = cos(-abs_beamport_angle*DEG2RAD)*xp0 - sin(-abs_beamport_angle*DEG2RAD)*zp0;
      zr0 = cos(-abs_beamport_angle*DEG2RAD)*zp0 + sin(-abs_beamport_angle*DEG2RAD)*xp0;
      xr1 = cos(-abs_beamport_angle*DEG2RAD)*xp1 - sin(-abs_beamport_angle*DEG2RAD)*zp1;
      zr1 = cos(-abs_beamport_angle*DEG2RAD)*zp1 + sin(-abs_beamport_angle*DEG2RAD)*xp1;
      dxr0 = cos(-abs_beamport_angle*DEG2RAD)*dxp0;
      dzr0 = sin(-abs_beamport_angle*DEG2RAD)*dxp0;
      dxr1 = - sin(-abs_beamport_angle*DEG2RAD)*dzp0;
      dzr1 = cos(-abs_beamport_angle*DEG2RAD)*dzp0;

      line(xr0+dxr0, yp0, zr0+dzr0, -xr0+dxr0, yp0, -zr0+dzr0);
      line(xr1+dxr0, yp0, zr1+dzr0, -xr1+dxr0, yp0, -zr1+dzr0);
      line(xr0-dxr0, yp0, zr0-dzr0, -xr0-dxr0, yp0, -zr0-dzr0);
      line(xr1-dxr0, yp0, zr1-dzr0, -xr1-dxr0, yp0, -zr1-dzr0);
      line(xr0+dxr0, yp0, zr0+dzr0, xr0-dxr0, yp0, zr0-dzr0);
      line(xr1+dxr0, yp0, zr1+dzr0, xr1-dxr0, yp0, zr1-dzr0);
      line(-xr0+dxr0, yp0, -zr0+dzr0, -xr0-dxr0, yp0, -zr0-dzr0);
      line(-xr1+dxr0, yp0, -zr1+dzr0, -xr1-dxr0, yp0, -zr1-dzr0);

      /* Wings: */
      line(xr0+2*dxr0, yp0, zr0+2*dzr0, 2*dxr0, yp0, 2*dzr0);
      line(-xr1+2*dxr0, yp0, -zr1+2*dzr0, 2*dxr0, yp0, 2*dzr0);
      line(xr0+2*dxr0, yp0, zr0+2*dzr0,-xr1+2*dxr0, yp0, -zr1+2*dzr0);

      line(-xr0-2*dxr0, yp0, -zr0-2*dzr0, -2*dxr0, yp0, -2*dzr0);
      line(xr1-2*dxr0, yp0, zr1-2*dzr0, -2*dxr0, yp0, -2*dzr0);
      line(xr1-2*dxr0, yp0, zr1-2*dzr0,-xr0-2*dxr0, yp0, -zr0-2*dzr0);
      /* Back to drawing the emission planes... */
      yp1 = yp0;
      yp0 = yp0 + yheight_c;
      xp0 = -0.16;

    }

    /* Dashed lines for indicating "beam extraction" area... */
    dashed_line(cx+t1x*0.2,cy-yheight_c/2.0, cz+t1z*0.2, cx+t1x*r_empty, cy-yheight_c/2.0, cz+t1z*r_empty,10);
    dashed_line(cx+t1x*0.2, cy+yheight_c/2.0, cz+t1z*0.2, cx+t1x*r_empty, cy+yheight_c/2.0, cz+t1z*r_empty,10);
    dashed_line(cx+t2x*0.2, cy-yheight_c/2.0, cz+t2z*0.2, cx+t2x*r_empty, cy-yheight_c/2.0, cz+t2z*r_empty,5);
    dashed_line(cx+t2x*0.2, cy+yheight_c/2.0, cz+t2z*0.2, cx+t2x*r_empty, cy+yheight_c/2.0, cz+t2z*r_empty,5);

    /* Rectangle indicating the chosen focus rectangle - where the optics starts... */
    rectangle("xy",tx,ty,tz,focus_xw,focus_yh);

    /* Circles indicating extent of the "empty" zone where optics is not allowed */
    circle("xz", cx,  cy+yheight_c/2.0, cz, r_empty);
    circle("xz", cx,  cy-yheight_c/2.0, cz, r_empty);

    /* Circles indicating the builk shielding of the target monolith at 6 m */
    circle("xz", cx,  cy+focus_yh/2.0 , cz, 6);
    circle("xz", cx, cy-focus_yh/2.0 , cz, 6);
    circle("xz", cx,  cy+2, cz, 6);
    circle("xz", cx, cy-2, cz, 6);

    /* Arrow indicating proton beam direction */
    ax=sin(DEG2RAD*(90+abs_beamport_angle));az=cos(DEG2RAD*(90+abs_beamport_angle));
    bbx=sin(DEG2RAD*(110+abs_beamport_angle));bbz=cos(DEG2RAD*(110+abs_beamport_angle));
    ccx=sin(DEG2RAD*(70+abs_beamport_angle));ccz=cos(DEG2RAD*(70+abs_beamport_angle));
    line(-0.15*ax,0,-0.15*az,-5.9*ax,0,-5.9*az);
    line(-0.15*ax,0,-0.15*az,-0.15*ax-0.1*bbx,0,-0.15*az-0.1*bbz);
    line(-0.15*ax,0,-0.15*az,-0.15*ax-0.1*ccx,0,-0.15*az-0.1*ccz);

    /* Rectangle indicating the chosen focus rectangle - where the optics starts... */
    rectangle("xy",tx,ty,tz,focus_xw,focus_yh);
  }

  /* This last bit is relevant only in connection with 2014 source definition */
  if (strcasestr(sourcedef,"2014")) {

    line(cx-5.9,cy,cz,cx-cyl_radius-0.1,cy,cz);
    line(cx-cyl_radius-0.2*cos(10*DEG2RAD),cy,cz+0.2*sin(10*DEG2RAD),cx-cyl_radius-0.1,cy,cz);
    line(cx-cyl_radius-0.2*cos(10*DEG2RAD),cy,cz-0.2*sin(10*DEG2RAD),cx-cyl_radius-0.1,cy,cz);

    /* Rectangle indicating the chosen focus rectangle - where the optics starts... */
    rectangle("xy",tx,ty,tz,focus_xw,focus_yh);
    rectangle("xy",0,0,0,xwidth_c,yheight_c);
    rectangle("xy",-xwidth_c/2-xwidth_t/2,0,0,xwidth_t,yheight_t);
    rectangle("xy",xwidth_c/2+xwidth_t/2,0,0,xwidth_t,yheight_t);
  }

%}

END