/**************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2006, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Pol_Bfield
*
* %I
* Written by: Erik B Knudsen, Peter Christiansen and Peter Willendrup
* Date: July 2011
* Origin: RISOE
*
* Magnetic field component.
*
* %D
*
* Region with a definable magnetic field.
*
* The component is nestable if the concentric flag is set (the default). This means that it may have a
* construction like the following:
* // START MAGNETIC FIELD
* COMPONENT msf =
* Pol_Bfield(xwidth=0.08, yheight=0.08, zdepth=0.2, Bx=0, By=-0.678332e-4, Bz=0)
*      AT (0, 0, 0) RELATIVE armMSF
*
* // OTHER COMPONENTS INSIDE THE MAGNETIC FIELD CAN BE PLACED HERE
*
* // STOP MAGNETIC FIELD
* COMPONENT msf_stop = Pol_simpleBfield_stop(magnet_comp_stop=msf)
*      AT (0, 0, 0) RELATIVE armMSF
*
* Note that if you have objects within the magnetic field that extend outside it you may get
* wrong results, because propagation within the field will then possibly extend outside, e.g.
* when using a tabled field. The evaluated field will then use the nearest defined field point
* _also_ outside the defintion area. If these outside-points have a non-zero field precession will
* continue - even after the neutron has left the field region.
*
* In between the two component instances the propagation routine
* PROP_DT also handles spin propagation.
* The current algorithm used for spin propagation is:
* SimpleNumMagnetPrecession
* in pol-lib.
*
* Example: Pol_Bfield(xwidth=0.1, yheight=0.1, zdepth=0.2, Bx=0, By=1, Bz=0)
*          Pol_Bfield(xwidth=0.1, yheight=0.1, zdepth=0.2, field_type=1)
*
* Functions supplied by the system are (the number is the ID of the function to be given as the field_type parameter:
* 1. Constant magnetic field: Constant field (Bx,By,Bz) within the region
* 2. Rotating magnetic_field: Field is initially (0,By,0) but after a length of zdepth
*      has rotated to (By,0,0)
* 3. Majorana magnetic_field: Field is initially (Bx,By,0) with By<<Bx, then linearly transforms to
*      (-Bx,By,0) at z=zdepth.
* 4. MSF field:
* 5. RF field: A radio frequency field is modeled by an implcit use of a roating frame.
* 6. Gradient field: Similar to Majorana by without an x-component to the field. I.e. it (0,By,0) at z=0, and
*    becomes (0,-By,0) AT z=zdepth.
*
* If the concentric parameter is set to 1 the magnetic field is turned on by an
* instance of this component, and turned off by an instance of Pol_simpleBfield_stop.
* Anything in between is considered inside the field.
* If concentric is zero the field is considered a closed component - neutrons are propagated
* through the field region, and no other components are permitted inside the field.
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m]               Width of opening.
* yheight: [m]              Height of opening.
* zdepth: [m]               Length of field.
* radius: [m]               Radius of field if it is cylindrical or spherical.
* Bx: [T]                   Parameter used for x composant of field.
* By: [T]                   Parameter used for y composant of field.
* Bz: [T]                   Parameter used for z composant of field.
* field_type: []            ID of function describing the magnetic field. 1:constant field, 2: rotating fiel, 3: majorana type field, 4: MSF, 5: RF-field, 6: gradient field.
* concentric: [ ]           Allow components and/or other fields inside field. Requires a subsequent Pol_simpleBfield_stop component.
*
* %E
****************************************************************************/

DEFINE COMPONENT Pol_Bfield

SETTING PARAMETERS (int field_type=0, xwidth=0, yheight=0,zdepth=0, radius=0, Bx=0, By=0, Bz=0, concentric=1)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

SHARE
%{
  %include "pol-lib"

  double fmax (double, double);
  double fmin (double, double);
%}


DECLARE
%{
  /* Larmor frequency and scalar product threshold*/
  double Bprms[8];
  int shape;
%}

INITIALIZE
%{
  enum shapes { NONE = 0, BOX, WINDOW, CYLINDER, SPHERE, ANY };
  enum field_functions ff = field_type;
  /* these default field functions are part of pol-lib */
  if (ff == constant) {
    double* t = Bprms;
    t[0] = Bx;
    t[1] = By;
    t[2] = Bz;
  } else if (ff == rotating) {
    double* t = Bprms;
    t[0] = By;
    t[1] = zdepth;
  } else if (ff == majorana) {
    double* t = Bprms;
    t[0] = Bx;
    t[1] = By;
    t[2] = zdepth;
  }

  if (xwidth && yheight && zdepth) {
    shape = BOX;
  } else if (xwidth && yheight && !zdepth) {
    shape = WINDOW;
  } else if (radius && yheight) {
    shape = CYLINDER;
  } else if (radius) {
    shape = SPHERE;
  } else {
    shape = NONE;
  }

  if (shape == WINDOW && !concentric) {
    printf ("Warning (%s): Magnetic field has window shape and concentric==0 => Field volume == 0.\n", NAME_CURRENT_COMP);
  }
  if (shape == NONE) {
    printf ("Warning (%s): Magnetic field has no geometry. Full Z=0-plane is used as boundary.\n", NAME_CURRENT_COMP);
  }
%}

TRACE
%{
  double t0, t1;
  int hit;
  enum shapes { NONE = 0, BOX, WINDOW, CYLINDER, SPHERE, ANY };

  /*enter through whatever object we are*/
  switch (shape) {
  case BOX:
    hit = box_intersect (&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
    /*terminate neutrons which miss the component*/
    if (!hit)
      ABSORB;
    /*If we do hit - propagate to the start of the field unless the neutron is already there.*/
    if (t0 > FLT_EPSILON) {
      PROP_DT (t0);
      t1 -= t0;
    } else if (t0 < -FLT_EPSILON && concentric) {
      PROP_DT (t1);
    }
    break;
  case CYLINDER:
    hit = cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, radius, yheight);
    /*terminate neutrons which miss the component*/
    if (!hit)
      ABSORB;
    /*If we do hit - propagate to the start of the field unless the neutron is already there.*/
    if (t0 > FLT_EPSILON) {
      PROP_DT (t0);
      t1 -= t0;
    } else if (t0 < -FLT_EPSILON && concentric) {
      PROP_DT (t1);
    }
    break;
  case WINDOW:
    PROP_Z0;
    if (2 * x > xwidth || 2 * x < -xwidth || 2 * y > yheight || 2 * y < -yheight) {
      /*terminate neutrons which miss the component*/
      ABSORB;
    }
    break;
  default:
    PROP_Z0;
  }
  mcmagnet_push (_particle, field_type, &(ROT_A_CURRENT_COMP), &(POS_A_CURRENT_COMP), 0, Bprms);
  #ifdef MCDEBUG
  mcmagnet_print_stack ();
  #endif

  /*does the field "close/stop" itself*/
  if (!concentric) {
    switch (shape) {
    case BOX:
      PROP_DT (t1);
      break;
    case CYLINDER:
      PROP_DT (t1);
      break;
    case WINDOW:
      PROP_Z0;
      /*terminate neutrons which miss the exit window*/
      if (2 * x > xwidth || 2 * x < -xwidth || 2 * y > yheight || 2 * y < -yheight) {
        ABSORB;
      }
      break;
    default:
      PROP_Z0;
    }
    mcmagnet_pop (_particle);
  }
%}

MCDISPLAY
%{
  enum shapes { NONE = 0, BOX, WINDOW, CYLINDER, SPHERE, ANY };
  switch (shape) {
  case BOX:
    box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0);
    break;
  case CYLINDER:
    circle ("xz", 0, yheight / 2.0, 0, radius);
    circle ("xz", 0, -yheight / 2.0, 0, radius);
    line (-radius, -yheight / 2.0, 0, -radius, yheight / 2.0, 0);
    line (radius, -yheight / 2.0, 0, radius, yheight / 2.0, 0);
    line (0, -yheight / 2.0, -radius, 0, yheight / 2.0, -radius);
    line (0, -yheight / 2.0, radius, 0, yheight / 2.0, radius);
    break;
  case SPHERE:
    circle ("xz", 0, 0, 0, radius);
    circle ("xy", 0, 0, 0, radius);
    circle ("yz", 0, 0, 0, radius);
    break;
  case WINDOW:
    rectangle ("xy", 0, 0, 0, xwidth, yheight);
    break;
  }
%}

END