/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Guide_four_side.comp
*
* Component: Guide_four_side
*
* %I
* Written by: Tobias Panzner
* Date: 02/12/2009
* Version: $Revision: 1.0 $
* Release: McStas 1.12b
* Origin: PSI
*
* %D
*  This component models a guide with elliptic, parabolic or straight shaped side walls.
*  All four sides of the guide are independent from each other.
*  In the elliptic case the side wall shape follows the equation x^2/b^2+(z+z0)^2/a^2=1
*  (the center of the ellipse is located at (0,-z0)).
*  In the parabolic case the side wall shape follows the equation z=b-ax^2;mc
*  In the straight case the side wall shape follows the equation z=l/(w2-w1)*x-w1.
*
*  The shape selection is done by the focal points. The focal points are located at the 
*  z-axis and are defined by their distance to the entrance or exit window of the guide
*  (in the following called 'focal length').
*
*  If both focal lengths for one wall are zero it will be a straight wall (entrance and 
*  exit width have to be given in the beginning).
*
*  If one of the focal lengths is not zero the shape will be parabolic (only the entrance width 
*  given in the beginning is recognized; exit width will be calculated). If the the entrance 
*  focal length is zero the guide will be a focusing devise.
*  If the exit focal length is zero it will be defocusing devise.
*
*  If both focals are non zero the shape of the wall will be elliptic (only the entrance width 
*  given in the beginning is recognized; exit width will be calculated). 
*
*  Notice: 1.)The focal points are in general located outside the guide (positive focal lengths).
*             Focal points inside the guide need to have negative focal lengths.
*          2.)The exit width parameters (w2r, w2l, h2u,h2d) are only taken into account if the 
*             walls have a linear shape. In the ellitic or parabolic case they will be ignored.
*
*  The outer walls of each side is calculated by the component in depentence of the wallthickness 
*  and the shape of the inner walls.
*
*  Each of the 8 walls (one INNER and one OUTER for every side) can have a own indepenting 
*  reflecting layer or it can be a absorber or it can be transparent.
*
*  The reflectivity properties can be given by an input file (Format [q(Angs-1) R(0-1)]) or by 
*  parameters (Qc, alpha, m, W).
*
* %BUGS
* This component does not work with gravitation on.
*
* This component does not work correctly in GROUP-modus.
*
* %P
* INPUT PARAMETERS (62) :
*
* w1r:     [m]    Width at the right guide entry   (negative x-axis)
*                 (DEFAULT = 0)
*
* w2r:     [m]    Width at the right guide exit    (negative x-axis)
*                 (DEFAULT = 0)
*
* w1l:     [m]    Width at the left guide entry    (positive x-axis)
*                 (DEFAULT = 0)
*
* w2l:     [m]    Width at the left guide exit     (positive x-axis)
*                 (DEFAULT = 0)
*
* h1d:     [m]    Height at the bottom guide entry (negative y-axis)
*                 (DEFAULT = 0)
*
* h2d:     [m]    Height at the bottom guide entry (negative y-axis)
*                 (DEFAULT = 0)
*
* h1u:     [m]    Height at the top guide entry    (positive y-axis)
*                 (DEFAULT = 0)
*
* h2u:     [m]    Height at the top guide entry    (positive y-axis)
*                 (DEFAULT = 0)
*
* linwr     [m]   right horizontal wall: distance of 1st focal point
*                 and guide entry (DEFAULT = 0)
*
* loutwr   [m]    right horizontal wall: distance of 2nd focal point
*                 and guide exit (DEFAULT = 0)
*
* linwl    [m]    left horizontal wall: distance of 1st focal point
*                 and guide entry (DEFAULT = 0)
*
* loutwl   [m]    left horizontal wall: distance of 2nd focal point
*                 and guide exit  (DEFAULT = 0)
*
* linhd    [m]    lower vertical wall: distance of 1st focal point
*                 and guide entry (DEFAULT = 0)
*
* louthd   [m]    lower vertical wall: distance of 2nd focal point
*                 and guide exit (DEFAULT = 0)
*
* linhu    [m]    upper vertical wall: distance of 1st focal point
*                 and guide entry (DEFAULT = 0) 
*
* louthu   [m]    upper vertical wall: distance of 2nd focal point
*                 and guide exit (DEFAULT = 0)
*
*
* l:       [m]    length of guide (DEFAULT = 0)
*
* R0:      [1]    Low-angle reflectivity (DEFAULT = 0.99)
*
* RIreflect: (str)  Name of relfectivity file for the right inner wall.
*                   Format [q(Angs-1) R(0-1)]  (DEFAULT : no file)
*
* LIreflect: (str)  Name of relfectivity file for the left inner wall.
*                   Format [q(Angs-1) R(0-1)]   (DEFAULT : no file)
* 
* DIreflect: (str)  Name of relfectivity file for the bottom inner wall.
*                   Format [q(Angs-1) R(0-1)] (DEFAULT : no file)
*
* UIreflect: (str)  Name of relfectivity file for the top inner wall. 
*                   Format [q(Angs-1) R(0-1)]    (DEFAULT : no file)
*
* Qcxr:    [AA-1] Critical scattering vector for right vertical
*                 inner wall (DEFAULT = 0.0217)
*
* Qcxl:    [AA-1] Critical scattering vector for left vertical
*                 inner wall (DEFAULT = 0.0217)
*
* Qcyd:    [AA-1] Critical scattering vector for bottom inner wall
*                 (DEFAULT = 0.0217)
*
* Qcyu:    [AA-1] Critical scattering vector for top inner wall
*                 (DEFAULT = 0.0217)
*
* alphaxr: [AA]   Slope of reflectivity for right vertical
*                 inner wall (DEFAULT = 6.07)
*
* alphaxl: [AA]   Slope of reflectivity for left vertical
*                 inner wall (DEFAULT = 6.07)
*
* alphayd: [AA]   Slope of reflectivity for bottom inner wall
*                 (DEFAULT = 6.07)
*
* alphayu: [AA]   Slope of reflectivity for top inner wall
*                 (DEFAULT = 6.07)
*
* mxr:     [1]    m-value of material for right vertical inner wall.
*                  0 means the wall is absorbing.
*		  -1 means the wall is transparent.
*                 (DEFAULT = 3.6)
*
* mxl:     [1]    m-value of material for left vertical inner wall.
*                  0 means the wall is absorbing.
*		  -1 means the wall is transparent.
*                 (DEFAULT = 3.6)
*
* myd:     [1]    m-value of material for bottom inner wall
*                  0 means the wall is absorbing.
*		  -1 means the wall is transparent.
*                 (DEFAULT = 3.6)
*
* myu:     [1]    m-value of material for top inner wall
*                  0 means the wall is absorbing.
*		  -1 means the wall is transparent.
*                 (DEFAULT = 3.6)
*
* Wxr:       [AA-1] Width of supermirror cut-off for right inner wall
*                   (DEFAULT = 0.003)
*
* Wxl:       [AA-1] Width of supermirror cut-off for left inner wall
*                   (DEFAULT = 0.003)
*
* Wyu:       [AA-1] Width of supermirror cut-off for top inner wall
*                   (DEFAULT = 0.003)
*
* Wyd:       [AA-1] Width of supermirror cut-off for bottom inner wall
*                   (DEFAULT = 0.003)
*
* rwallthick: [m]   thickness of the right wall (DEFAULT = 0.001 m)
*
* lwallthick: [m]   thickness of the left wall  (DEFAULT = 0.001 m)
*
* uwallthick: [m]   thickness of the top wall   (DEFAULT = 0.001 m)
*
* dwallthick: [m]   thickness of the bottom wall(DEFAULT = 0.001 m)
*
*  IF NEUTRONS HIT THE WALLS FROM OUTSIDE, THE OUTER WALLS CAN BE DEFINED BY
*
* ROreflect: (str)  Name of relfectivity file for the right outer wall.
*                   Format [q(Angs-1) R(0-1)]  (DEFAULT : no file)
*
* LOreflect: (str)  Name of relfectivity file for the left outer wall.
*                   Format [q(Angs-1) R(0-1)]   (DEFAULT : no file)
* 
* DOreflect: (str)  Name of relfectivity file for the bottom outer wall.
*                   Format [q(Angs-1) R(0-1)] (DEFAULT : no file)
*
* UOreflect: (str)  Name of relfectivity file for the top outer wall.
*                   Format [q(Angs-1) R(0-1)]    (DEFAULT : no file)
*
* QcxrOW:    [AA-1] Critical scattering vector for right vertical
*                   outer wall (DEFAULT = 0.0217)
*
* QcxlOW:    [AA-1] Critical scattering vector for left vertical
*                   outer wall (DEFAULT = 0.0217)
*
* QcydOW:    [AA-1] Critical scattering vector for bottom outer wall
*                   (DEFAULT = 0.0217)
*
* QcyuOW:    [AA-1] Critical scattering vector for top outer wall
*                   (DEFAULT = 0.0217)
*
* alphaxrOW: [AA]   Slope of reflectivity for right vertical
*                   outer wall (DEFAULT = 6.07)
*
* alphaxlOW: [AA]   Slope of reflectivity for left vertical
*                   outer wall (DEFAULT = 6.07)
*
* alphaydOW: [AA]   Slope of reflectivity for bottom outer wall
*                   (DEFAULT = 6.07)
*
* alphayuOW: [AA]   Slope of reflectivity for top outer wall
*                   (DEFAULT = 6.07)
*
* mxrOW:     [1]   m-value of material for right vertical outer wall
*                  0 means the wall is absorbing. (DEFAULT)
*		  -1 means the wall is transparent.
*                  (DEFAULT = 0)
*
* mxlOW:     [1]   m-value of material for left vertical outer wall
*                  0 means the wall is absorbing.(DEFAULT)
*		  -1 means the wall is transparent.
*                  (DEFAULT = 0)
*
* mydOW:     [1]   m-value of material for bottom outer wall
*                  0 means the wall is absorbing. (DEFAULT)
*		  -1 means the wall is transparent.
*                  (DEFAULT = 0)
*
* myuOW:     [1]   m-value of material for top outer wall
*                  0 means the wall is absorbing. (DEFAULT)
*		  -1 means the wall is transparent.
*                  (DEFAULT = 0)
*
* WxrOW:       [AA-1] Width of supermirror cut-off for right outer wall
*                     (DEFAULT = 0.003)
*
* WxlOW:       [AA-1] Width of supermirror cut-off for left outer wall
*                     (DEFAULT = 0.003)
*
* WyuOW:       [AA-1] Width of supermirror cut-off for top outer wall
*                     (DEFAULT = 0.003)
*
* WydOW:       [AA-1] Width of supermirror cut-off for bottom outer wall
*                     (DEFAULT = 0.003)
*
*
* %End
*
******************************************************************************/

DEFINE COMPONENT Guide_four_side
DEFINITION PARAMETERS(string RIreflect=0, LIreflect=0, UIreflect=0, DIreflect=0, ROreflect=0, LOreflect=0, UOreflect=0, DOreflect=0)
SETTING PARAMETERS (w1l=0,w2l=0,linwl=0,loutwl=0,
                    w1r=0,w2r=0,linwr=0,loutwr=0,
                    h1u=0,h2u=0,linhu=0,louthu=0,
                    h1d=0,h2d=0,linhd=0,louthd=0,
                    l=0, R0=0.99,
                    Qcxl=0.0217,Qcxr=0.0217,Qcyu=0.0217, Qcyd=0.0217,
                    alphaxl=6.07, alphaxr=6.07, alphayu=6.07, alphayd=6.07,
                    Wxr=0.003,Wxl=0.003,Wyu=0.003,Wyd=0.003,
                    mxr=3.6, mxl=3.6, myu=3.6, myd=3.6,
                    QcxrOW=0.0217,QcxlOW=0.0217,QcyuOW=0.0217, QcydOW=0.0217,
                    alphaxlOW=6.07, alphaxrOW=6.07, alphayuOW=6.07, alphaydOW=6.07,
                    WxrOW=0.003,WxlOW=0.003,WyuOW=0.003,WydOW=0.003,
                    mxrOW=0, mxlOW=0, myuOW=0, mydOW=0,
		    rwallthick=0.001,lwallthick=0.001,uwallthick=0.001,dwallthick=0.001)
OUTPUT PARAMETERS(w1rwt,w1lwt,h1uwt,h1dwt,w2rwt,w2lwt,h2uwt,h2dwt,
                  pawr,pawl,pbwr,pbwl,pahu,pahd,pbhu,pbhd,
                  awl,bwl,awr,bwr,ahu,bhu,ahd,bhd,
                  awlwt,bwlwt,awrwt,bwrwt,ahuwt,bhuwt,ahdwt,bhdwt,
                  pawrwt,pawlwt,pbwrwt,pbwlwt,pahuwt,pahdwt,pbhuwt,pbhdwt,
                  t1,t2w1r,t2w2r,t2w1l,t2w2l,t2h1u,t2h2u,t2h1d,t2h2d,
                  t2w1rwt,t2w2rwt,t2w1lwt,t2w2lwt,t2h1uwt,t2h2uwt,t2h1dwt,t2h2dwt,
                  a2wlwt,b2wlwt,a2wrwt,b2wrwt,a2huwt,b2huwt,a2hdwt,b2hdwt,
                  a2wl,b2wl,a2wr,b2wr,a2hu,b2hu,a2hd,b2hd,
                  mru1,mru2,nru1,nru2,mrd1,mrd2,nrd1,nrd2,mlu1,mlu2,nlu1,nlu2,mld1,mld2,nld1,nld2,
                  z0wr,z0wl,z0hu,z0hd,
                  lbwl,lbwr,lbhu,lbhd,
                  p2parawr,p2parawl,p2parahu,p2parahd,
                  u1wl,u2wl,u1wr,u2wr,u1hu,u2hu,u1hd,u2h,
                  DIV1wl,DIV1wr,DIV1hu,DIV1hd,
                  p2parawrwt,p2parawlwt,p2parahuwt,p2parahdwt,
                  DIV1wlwt,DIV1wrwt,DIV1huwt,DIV1hdwt,
                  u1wrwt,u2wrwt,u1wlwt,u2wlwt,u1huwt,u2huwt,u1hdwt,u2hdwt,
                  dxwr,dzwr,dxwl,dzwl,dyhu,dzhu,dyhd,dzhd,
                  zw1r,zw2r,zw1l,zw2l,zh1u,zh2u,zh1d,zh2d,m,n,
                  pwr,qwr,pwl,qwl,phu,qhu,phd,qhd,
                  nz,nx,ny,n2,
                  pf,
                  vxin,vyin,vzin,
                  q,
                  xintersecwr,xintersecwl,yintersechd,yintersechu,
                  xlimitr,xlimitrwt,xlimitl,xlimitlwt,ylimitd,ylimitdwt,ylimitu,ylimituwt,
                  xtest,ytest,
                  pwrwt,qwrwt,pwlwt,qwlwt,phuwt,qhuwt,phdwt,qhdwt,
                  zw1rwt,zw2rwt,zw1lwt,zw2lwt,zh1uwt,zh2uwt,zh1dwt,zh2dwt,
                  xintersecwrwt,xintersecwlwt,yintersechdwt,yintersechuwt,
                  riTable, liTable, uiTable, diTable,
                  roTable, loTable, uoTable, doTable)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,sx,sy,p)

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

DECLARE
%{
  double w1rwt,w1lwt,h1uwt,h1dwt,w2rwt,w2lwt,h2uwt,h2dwt;					/* entrance and exit width for th OUTER walls of the guide*/
  double pawr,pawl,pbwr,pbwl,pahu,pahd,pbhu,pbhd;						/* parameter a and b of the parabolic curves (INNER walls)*/
  double awl,bwl,awr,bwr,ahu,bhu,ahd,bhd;							/* long and short axis a and b auf the ellipses (INNER walls)*/
  double awlwt,bwlwt,awrwt,bwrwt,ahuwt,bhuwt,ahdwt,bhdwt;					/* long and short axis a and b auf the ellipses (OUTER walls)*/
  double pawrwt,pawlwt,pbwrwt,pbwlwt,pahuwt,pahdwt,pbhuwt,pbhdwt;				/* parameter a and b of the parabolic curves for the OUTER wall*/
  double t1,t2w1r,t2w2r,t2w1l,t2w2l,t2h1u,t2h2u,t2h1d,t2h2d;					/* time variables (INNER walls)*/
  double t2w1rwt,t2w2rwt,t2w1lwt,t2w2lwt,t2h1uwt,t2h2uwt,t2h1dwt,t2h2dwt;			/* time variables (OUTER walls)*/

  double a2wlwt,b2wlwt,a2wrwt,b2wrwt,a2huwt,b2huwt,a2hdwt,b2hdwt;				/* square of long and short axis a and b auf the ellipses (OUTER walls) */
  double a2wl,b2wl,a2wr,b2wr,a2hu,b2hu,a2hd,b2hd;						/* square of long and short axis a and b auf the ellipses (INNER walls)*/
  double mru1,mru2,nru1,nru2,mrd1,mrd2,nrd1,nrd2,mlu1,mlu2,nlu1,nlu2,mld1,mld2,nld1,nld2;	/* variables the calculated the guide geometrie in the entrance and exit plane (absorbing mask given by the walls)*/
  double z0wr,z0wl,z0hu,z0hd;

  double lbwl,lbwr,lbhu,lbhd;									/* length between focal points , needed for elliptic case */
  double p2parawr,p2parawl,p2parahu,p2parahd;							/* help variables to calculate the parabolic curve parameters a and b (INNER walls)*/
  double u1wl,u2wl,u1wr,u2wr,u1hu,u2hu,u1hd,u2hd;						/* help variables to calculate the long axis a of the ellipse (INNER walls)*/
  double DIV1wl,DIV1wr,DIV1hu,DIV1hd;								/* help variables (INNER walls)*/
  double p2parawrwt,p2parawlwt,p2parahuwt,p2parahdwt;						/* help variables to calculate the parabolic curve parameters a and b for (OUTER wall)*/
  double DIV1wlwt,DIV1wrwt,DIV1huwt,DIV1hdwt;							/* help variables (OUTER wall) */
  double u1wrwt,u2wrwt,u1wlwt,u2wlwt,u1huwt,u2huwt,u1hdwt,u2hdwt;				/* help variables to calculate the long axis a of the ellipse (OUTER walls)*/
  double dxwr,dzwr,dxwl,dzwl,dyhu,dzhu,dyhd,dzhd;						/* help variables to calculate the entrance and exit witdh of the OUTER walls*/

  double zw1r,zw2r,zw1l,zw2l,zh1u,zh2u,zh1d,zh2d,m,n;						/* zcomponent of the intersection point of the neutron trajectory and the ellipse (INNER walls)*/
  double pwr,qwr,pwl,qwl,phu,qhu,phd,qhd;							/* parameters for the quadratic equations for the intersection points or the intersection time (INNER walls)*/
  double nz,nx,ny,n2;										/* component and length of the surfaces normal vector at the intersection point */
  double pf;											/* prefactor to calculate the velocity vector after the interaction */
  double vxin,vyin,vzin;									/* velocity vector components before the interaction*/
  double q;											/* q-vector for the interaction */
  double xintersecwr,xintersecwl,yintersechd,yintersechu;					/* intersection point crosscheck variables (INNER walls) */
  double xlimitr,xlimitrwt,xlimitl,xlimitlwt,ylimitd,ylimitdwt,ylimitu,ylimituwt;		/* limit variables to determine the interaction position given by the time relative to the guide walls*/
  double xtest,ytest;										/* interaction position of the neutron given by the interaction time; crosscheck with limit variables*/
  double pwrwt,qwrwt,pwlwt,qwlwt,phuwt,qhuwt,phdwt,qhdwt;					/* parameters for the quadratic equations for the intersection points or the intersection time (OUTER walls)*/
  double zw1rwt,zw2rwt,zw1lwt,zw2lwt,zh1uwt,zh2uwt,zh1dwt,zh2dwt;				/* zcomponent of the intersection point of the neutron trajectory and the ellipse (OUTER walls)*/
  double xintersecwrwt,xintersecwlwt,yintersechdwt,yintersechuwt;				/* intersection point crosscheck variables (OUTER walls) */

  t_Table riTable,liTable,uiTable,diTable;
  t_Table roTable,loTable,uoTable,doTable;
%}


INITIALIZE
%{

 int i;

if (RIreflect && strlen(RIreflect)) {
    if (Table_Read(&riTable, RIreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"right inner Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, RIreflect));
  }

if (LIreflect && strlen(LIreflect)) {
    if (Table_Read(&liTable, LIreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"left inner Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, LIreflect));
  }

if (UIreflect && strlen(UIreflect)) {
    if (Table_Read(&uiTable, UIreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"top inner Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, UIreflect));
  }

if (DIreflect && strlen(DIreflect)) {
    if (Table_Read(&diTable, DIreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"botton inner Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, DIreflect));
  }

if (ROreflect && strlen(ROreflect)) {
    if (Table_Read(&roTable, ROreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"right outer Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, ROreflect));
  }

if (LOreflect && strlen(LOreflect)) {
    if (Table_Read(&loTable, LOreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"left outer Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, LOreflect));
  }

if (UOreflect && strlen(UOreflect)) {
    if (Table_Read(&uoTable, UOreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"top outer Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, UOreflect));
  }

if (DOreflect && strlen(DOreflect)) {
    if (Table_Read(&doTable, DOreflect, 1) <= 0) 							/* read 1st block data from file into pTable */
      exit(fprintf(stderr,"botton outer Wall: %s: can not read file %s\n", NAME_CURRENT_COMP, DOreflect));
  }

if (w1r < 0)											/* Ceckup for the input variables */
  {
    fprintf(stderr,"Component: %s (Guide_four_side) w1r should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (w1l < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) w1l should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (h1u < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) h1u should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (h1d < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) h1d should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (w2r < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) w2r should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (w2l < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) w2l should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (h2u < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) h2u should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (h2d < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) h2d should  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    fprintf(stderr,"           (for negative values use the global guide position !) \n");
    exit(-1);
  }

if (mxrOW !=-1 && mxrOW < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) mxrOW must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (mxlOW !=-1 && mxlOW<0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) mxlOW must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (myuOW !=-1 && myuOW<0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) myuOW must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (mydOW !=-1 && mydOW<0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) mydOW must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (Qcxl < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Qcxl can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (Qcxr < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Qcxr can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (Qcyu < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Qcyu can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (Qcyd < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Qcyd can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (alphaxl < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) alphaxl can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (alphaxr < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) alphaxr can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (alphayu < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) alphayu can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (alphayd < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) alphayd can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (rwallthick < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) rwallthick can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (lwallthick < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) lwallthick can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }
if (uwallthick < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) uwallthick can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (dwallthick < 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) dwallthick can  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           NOT be negative \n");
    exit(-1);
  }

if (mxr < 0 && mxr!=-1)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) mxr must  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (mxl < 0 && mxl!=-1)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) mxl must  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (myu < 0 && myu!=-1)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) myu must  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (myd < 0 && myd!=-1)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) myd must  \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be -1 (transperent) or \n");
    fprintf(stderr,"           be  0 (absorbing) or \n");
    fprintf(stderr,"           be  > 0 (reflecting) \n");
    exit(-1);
  }

if (Wxr <=0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Wxr must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be positive\n");
    exit(-1);
  }

if (Wxl <=0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Wxl must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be positive\n");
    exit(-1);
  }

if (Wyu <=0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Wyu must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be positive\n");
    exit(-1);
  }

if (Wyd <=0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) Wyd must \n", NAME_CURRENT_COMP);
    fprintf(stderr,"           be positive\n");
    exit(-1);
  }

if (l <= 0)
  {
    fprintf(stderr,"Component: %s (Guide_four_side) real guide length \n",
    NAME_CURRENT_COMP);
    fprintf(stderr,"           is <= ZERO ! \n");
    exit(-1);
  }

if (mcgravitation) fprintf(stderr,"WARNING: Guide_four_side: %s: "
    "This component produces wrong results with gravitation !\n"
    "Use Guide_gravity.\n",
    NAME_CURRENT_COMP);




													/* Calculation of curve-parameters for the right side wall - negative x-axis */

													/* elliptic right side wall */

lbwr=linwr+l+loutwr;											/* lenght between the two focal points for the right wall */
z0wr=(linwr-l-loutwr)/2.0;										/* calculation of the z-component of the center of the ellipse for the right wall (with opposite sign)*/

if(loutwr!=0 && linwr!=0)										/* ? elliptic shape ? */
  { 
    u1wr=sqrt((linwr*linwr)+(w1r*w1r));									/* length between entrance focal point and starting point of the elliptic wall (INNER wall)*/
    u2wr=sqrt((w1r*w1r)+((l+loutwr)*(l+loutwr)));							/* length between exit focal point and end point of the elliptic wall (INNER wall) */
    awr=(u1wr+u2wr)/2.0;										/* long half axis a of the ellipse (INNER  wall)*/
    a2wr=awr*awr;											/* square of the long axis a (INNER wall)*/
    bwr=sqrt(a2wr-(lbwr*lbwr/4.0));									/* short half axis b of the ellipse  (INNER wall)*/
    b2wr=bwr*bwr;											/* square of short half axis b of the ellipse  (INNER wall)*/
    DIV1wr=sqrt(1.0-((lbwr/2.0-loutwr)*(lbwr/2.0-loutwr)/a2wr));					/* help variable to calculated the exit width (INNER wall)*/
    w2r=bwr*DIV1wr;											/* exit width (INNER wall)*/
    if(l<lbwr/2-loutwr){										/* if the maximum opening of  the guide is smaller than the small half axis b, the OUTER wall is defined by: */
        dxwr=rwallthick*sin(atan(a2wr*w1r/(b2wr*z0wr)));						/* help variable dx; needed because the wall is not paralell to the z-axis or the wallthickness is not perpendicular to z*/
        dzwr=rwallthick*cos(atan(a2wr*w1r/(b2wr*z0wr)));						/* help variable dz; needed because the wall is not paralell to the z-axis or the wallthickness is not perpendicular to z*/
    	u1wrwt=sqrt(((linwr+dzwr)*(linwr+dzwr))+((w1r+dxwr)*(w1r+dxwr)));				/* length between entrance focal point and starting point of the elliptic wall (OUTER wall)*/
    	u2wrwt=sqrt(((w1r+dxwr)*(w1r+dxwr))+((l+loutwr-dzwr)*(l+loutwr-dzwr)));				/* length between exit focal point and end point of the elliptic wall (OUTER wall) */
    	awrwt=(u1wrwt+u2wrwt)/2.0;									/* long half axis a of the ellipse  (OUTER wall)*/
    	a2wrwt=awrwt*awrwt;										/* square of the long axis a (OUTER wall)*/
    	bwrwt=sqrt(a2wrwt-(lbwr*lbwr/4.0));								/* short half axis b of the ellipse  (OUTER wall)*/
    	b2wrwt=bwrwt*bwrwt;										/* square of short half axis b of the ellipse  (OUTER wall)*/
    	w2rwt=bwrwt*sqrt(1.0-((lbwr/2.0-loutwr)*(lbwr/2.0-loutwr)/a2wrwt));				/* exit width for OUTER right elliptic wall */
        w1rwt=bwrwt*sqrt(1.0-((lbwr/2.0-loutwr-l)*(lbwr/2.0-loutwr-l)/a2wrwt));				/* entrance width for OUTER right elliptic wall */
    }else{												/* if the maximum opening of  the guide is the small half axis bthe OUTER wall is defined by:*/
	bwrwt=bwr+rwallthick;										/* short half axis b of the ellipse  (OUTER wall)*/
        b2wrwt=bwrwt*bwrwt;										/* square of the long axis a (OUTER wall)*/
        awrwt=sqrt(b2wrwt+(lbwr*lbwr/4.0));								/* long half axis a of the ellipse  (OUTER  wall)*/
        a2wrwt=b2wrwt+(lbwr*lbwr/4.0);									/* square of short half axis b of the ellipse  (OUTER wall)*/
        w2rwt=bwrwt*sqrt(1.0-((lbwr/2.0-loutwr)*(lbwr/2.0-loutwr)/a2wrwt));				/* exit width for OUTER right elliptic wall */
        w1rwt=bwrwt*sqrt(1.0-((lbwr/2.0-linwr)*(lbwr/2.0-linwr)/a2wrwt));				/* entrance width for OUTER right elliptic wall */
    }
   }

													/* parabolic focusing right side wall */

if(linwr==0 && loutwr!=0){										/* ? parabolic focusing shape ? */
    DIV1wr=(l+loutwr)*(l+loutwr);									/* help variable to calculate the curve parameters (INNER wall) */
    p2parawr=2.0*(sqrt(DIV1wr+(w1r*w1r))-sqrt(DIV1wr));							/* help variable to calculate the curve parameters (INNER wall) */
    w2r=sqrt(p2parawr*(loutwr+p2parawr/4.0));								/* exit width for right focusing parabolic wall (INNER wall) */
    pbwr=l+loutwr+p2parawr/4.0;										/* parameter b for parabolic equation to define the wall (INNER wall)*/
    pawr=1.0/p2parawr;											/* parameter a for parabolic equation to define the wall (INNER wall)*/
    dxwr=rwallthick*sin(atan(w1r*2*pawr));								/* help variable dx; needed because the wall is not paralell to the z-axis or the wallthickness is not perpendicular to z*/
    dzwr=rwallthick*cos(atan(w1r*2*pawr));								/* help variable dz; needed because the wall is not paralell to the z-axis or the wallthickness is not perpendicular to z*/
    DIV1wrwt=(l+loutwr-dzwr)*(l+loutwl-dzwr);								/* help variable to calculate the curve parameters (OUTER wall) */
    p2parawrwt=2.0*(sqrt(DIV1wrwt+((w1r+dxwr)*(w1r+dxwr)))-sqrt(DIV1wrwt));				/* help variable to calculate the curve parameters (OUTER wall) */
    pbwrwt=l+loutwr+p2parawrwt/4.0;									/* parameter b for parabolic equation to define the wall (OUTER wall)*/
    pawrwt=1.0/p2parawrwt;										/* parameter a for parabolic equation to define the wall (OUTER wall)*/
    w2rwt=sqrt(p2parawrwt*(loutwr+p2parawrwt/4.0));							/* exit width for right focusing parabolic wall (OUTER wall) */
    w1rwt=sqrt(p2parawrwt*(loutwr+l+p2parawrwt/4.0));							/* entrance width for right focusing parabolic wall (OUTER wall) */
   }

													/* parabolic defocusing right side wall */

if (linwr!=0 && loutwr==0)										/* ? parabolic defocusing shape ? */
    {
      DIV1wr=linwr*linwr;										/* help variable to calculate the curve parameters (INNER wall) */
      p2parawr=2.0*(sqrt(DIV1wr+(w1r*w1r))-sqrt(DIV1wr));						/* help variable to calculate the curve parameters (INNER wall) */
      w2r=sqrt(p2parawr*(l+linwr+p2parawr/4.0));							/* exit width for right focusing parabolic wall (INNER wall) */
      pbwr=-(linwr+p2parawr/4.0);									/* parameter b for parabolic equation to define the wall (INNER wall)*/
      pawr=-1.0/p2parawr;										/* parameter a for parabolic equation to define the wall (INNER wall)*/
      dxwr=rwallthick*sin(atan(-w2r*2*pawr));								/* help variable dx; needed because the wall is not paralell to the z-axis or the wallthickness is not perpendicular to z*/
      dzwr=rwallthick*cos(atan(-w2r*2*pawr));								/* help variable dz; needed because the wall is not paralell to the z-axis or the wallthickness is not perpendicular to z*/
      DIV1wrwt=(linwr+l-dzwr)*(linwr+l-dzwr);								/* help variable to calculate the curve parameters (OUTER wall) */
      p2parawrwt=2.0*(sqrt(DIV1wrwt+((w2r+dxwr)*(w2r+dxwr)))-sqrt(DIV1wrwt));				/* help variable to calculate the curve parameters (OUTER wall) */
      w1rwt=sqrt(p2parawrwt*(linwr+p2parawrwt/4.0));							/* entrance width for right focusing parabolic wall (OUTER wall) */
      w2rwt=sqrt(p2parawrwt*(linwr+l+p2parawrwt/4.0));							/* exit width for right focusing parabolic wall (OUTER wall) */
      pbwrwt=-(linwr+p2parawrwt/4.0);									/* parameter b for parabolic equation to define the wall (OUTER wall)*/
      pawrwt=-1.0/p2parawrwt;										/* parameter a for parabolic equation to define the wall (OUTER wall)*/
    }

													/* straight right side wall */

if(linwr==0 && loutwr==0){										/* ? linear shape ? */
    w1rwt=w1r+rwallthick/(cos(atan((w1r-w2r)/l)));							/* entrance width for right focusing parabolic wall (OUTER wall) */
    w2rwt=w2r+rwallthick/(cos(atan((w1r-w2r)/l)));							/* exit width for right focusing parabolic wall (OUTER wall) */
    }


													/* Calculation of curve-parameters for the left side wall - positive x-axis - analog to right side*/

lbwl=linwl+l+loutwl;
z0wl=(linwl-l-loutwl)/2.0;

if((linwl!=0) && (loutwl!=0) ) 										/* elleptic left side wall */
  { 
    u1wl=sqrt((linwl*linwl)+(w1l*w1l));
    u2wl=sqrt((w1l*w1l)+((l+loutwl)*(l+loutwl)));
    awl=(u1wl+u2wl)/2.0;
    a2wl=awl*awl;
    bwl=sqrt(a2wl-(lbwl*lbwl/4.0));
    b2wl=bwl*bwl;
    DIV1wl=sqrt(1.0-((lbwl/2.0-loutwl)*(lbwl/2.0-loutwl)/a2wl));
    w2l=bwl*DIV1wl;
    if(l<lbwl/2-loutwl){
        dxwl=lwallthick*sin(atan(a2wl*w1l/(b2wl*z0wl)));
        dzwl=lwallthick*cos(atan(a2wl*w1l/(b2wl*z0wl)));
    	u1wlwt=sqrt(((linwl+dzwl)*(linwl+dzwl))+((w1l+dxwl)*(w1l+dxwl)));								
    	u2wlwt=sqrt(((w1l+dxwl)*(w1l+dxwl))+((l+loutwl-dzwl)*(l+loutwl-dzwl)));	
    	awlwt=(u1wlwt+u2wlwt)/2.0;
    	a2wlwt=awlwt*awlwt;
    	bwlwt=sqrt(a2wlwt-(lbwl*lbwl/4.0));
    	b2wlwt=bwlwt*bwlwt;
    	w2lwt=bwlwt*sqrt(1.0-((lbwl/2.0-loutwl)*(lbwl/2.0-loutwl)/a2wlwt));
        w1lwt=bwlwt*sqrt(1.0-((lbwl/2.0-loutwl-l)*(lbwl/2.0-loutwl-l)/a2wlwt));
    }else{
	bwlwt=bwl+lwallthick;
        b2wlwt=bwlwt*bwlwt;
        awlwt=sqrt(b2wlwt+(lbwl*lbwl/4.0));
        a2wlwt=b2wlwt+(lbwl*lbwl/4.0);
        w2lwt=bwlwt*sqrt(1.0-((lbwl/2.0-loutwl)*(lbwl/2.0-loutwl)/a2wlwt));
        w1lwt=bwlwt*sqrt(1.0-((lbwl/2.0-linwl)*(lbwl/2.0-linwl)/a2wlwt));
    }
   }


if(linwl==0 && loutwl!=0)										/* parabolic focusing left side wall */
   {
    DIV1wl=(l+loutwl)*(l+loutwl);
    p2parawl=2.0*(sqrt(DIV1wl+(w1l*w1l))-sqrt(DIV1wl));
    w2l=sqrt(p2parawl*(loutwl+p2parawl/4.0));
    pbwl=l+loutwl+p2parawl/4.0;
    pawl=1.0/p2parawl;
    dxwl=lwallthick*sin(atan(w1l*2*pawl));
    dzwl=lwallthick*cos(atan(w1l*2*pawl));
    DIV1wlwt=(l+loutwl-dzwl)*(l+loutwl-dzwl);
    p2parawlwt=2.0*(sqrt(DIV1wlwt+((w1l+dxwl)*(w1l+dxwl)))-sqrt(DIV1wlwt));
    pbwlwt=l+loutwl+p2parawlwt/4.0;
    pawlwt=1.0/p2parawlwt;
    w2lwt=sqrt(p2parawlwt*(loutwl+p2parawlwt/4.0));
    w1lwt=sqrt(p2parawlwt*(loutwl+l+p2parawlwt/4.0));
   }

if (linwl!=0 && loutwl==0)										/* parabolic defocusing left side wall */
    {
      DIV1wl=linwl*linwl;
      p2parawl=2.0*(sqrt(DIV1wl+(w1l*w1l))-sqrt(DIV1wl));
      w2l=sqrt(p2parawl*(l+linwl+p2parawl/4.0));
      pbwl=-(linwl+p2parawl/4.0);
      pawl=-1.0/p2parawl;
      dxwl=lwallthick*sin(atan(-w2l*2*pawl));
      dzwl=lwallthick*cos(atan(-w2l*2*pawl));
      DIV1wlwt=(linwl+l-dzwl)*(linwl+l-dzwl);
      p2parawlwt=2.0*(sqrt(DIV1wlwt+((w2l+dxwl)*(w2l+dxwl)))-sqrt(DIV1wlwt));
      w1lwt=sqrt(p2parawlwt*(linwl+p2parawlwt/4.0));
      w2lwt=sqrt(p2parawlwt*(linwl+l+p2parawlwt/4.0));
      pbwlwt=-(linwl+p2parawlwt/4.0);
      pawlwt=-1.0/p2parawlwt;
    }

if(linwl==0 && loutwl==0){
    w1lwt=w1l+lwallthick/(cos(atan((w2l-w1l)/l)));
    w2lwt=w2l+lwallthick/(cos(atan((w2l-w1l)/l)));
    }


													/* Calculation of curve-parameters for the top wall - positive y-axis - analog  right wall*/

lbhu=linhu+l+louthu;
z0hu=lbhu/2.0-l-louthu;	

if (linhu != 0 && louthu !=0)										/* elliptic top wall */
   {
    u1hu=sqrt((linhu*linhu)+(h1u*h1u));
    u2hu=sqrt((h1u*h1u)+((l+louthu)*(l+louthu)));
    ahu=(u1hu+u2hu)/2.0;
    bhu=sqrt(ahu*ahu-(lbhu*lbhu/4.0));
    a2hu=ahu*ahu;
    b2hu=bhu*bhu;
    DIV1hu=sqrt(1.0-((lbhu/2.0-louthu)*(lbhu/2.0-louthu)/a2hu));
    h2u=bhu*DIV1hu;
    if(l<lbhu/2-louthu){
        dyhu=uwallthick*sin(atan(a2hu*h1u/(b2hu*z0hu)));
        dzhu=uwallthick*cos(atan(a2hu*h1u/(b2hu*z0hu)));
    	u1huwt=sqrt(((linhu+dzhu)*(linhu+dzhu))+((h1u+dyhu)*(h1u+dyhu)));								
    	u2huwt=sqrt(((h1u+dyhu)*(h1u+dyhu))+((l+louthu-dzhu)*(l+louthu-dzhu)));	
    	ahuwt=(u1huwt+u2huwt)/2.0;
    	a2huwt=ahuwt*ahuwt;
    	bhuwt=sqrt(a2huwt-(lbhu*lbhu/4.0));
    	b2huwt=bhuwt*bhuwt;
    	h2uwt=bhuwt*sqrt(1.0-((lbhu/2.0-louthu)*(lbhu/2.0-louthu)/a2huwt));
        h1uwt=bhuwt*sqrt(1.0-((lbhu/2.0-louthu-l)*(lbhu/2.0-louthu-l)/a2huwt));
    }else{
	bhuwt=bhu+uwallthick;
        b2huwt=bhuwt*bhuwt;
        ahuwt=sqrt(b2huwt+(lbhu*lbhu/4.0));
        a2huwt=b2huwt+(lbhu*lbhu/4.0);
        h2uwt=bhuwt*sqrt(1.0-((lbhu/2.0-louthu)*(lbhu/2.0-louthu)/a2huwt));
        h1uwt=bhuwt*sqrt(1.0-((lbhu/2.0-linhu)*(lbhu/2.0-linhu)/a2huwt));
    }
  }

if(linhu==0 && louthu!=0)										/* parabolic focusing top wall */
   { 	
    DIV1hu=(l+louthu)*(l+louthu);
    p2parahu=2.0*(sqrt(DIV1hu+(h1u*h1u))-sqrt(DIV1hu));
    h2u=sqrt(p2parahu*(louthu+p2parahu/4.0));
    pbhu=l+louthu+p2parahu/4.0;
    pahu=1.0/p2parahu;
    dyhu=uwallthick*sin(atan(h1u*2*pahu));
    dzhu=uwallthick*cos(atan(h1u*2*pahu));
    DIV1huwt=(l+louthu-dzhu)*(l+louthu-dzhu);
    p2parahuwt=2.0*(sqrt(DIV1huwt+((h1u+dyhu)*(h1u+dyhu)))-sqrt(DIV1huwt));
    pbhuwt=l+louthu+p2parahuwt/4.0;
    pahuwt=1.0/p2parahuwt;
    h2uwt=sqrt(p2parahuwt*(louthu+p2parahuwt/4.0));
    h1uwt=sqrt(p2parahuwt*(louthu+l+p2parahuwt/4.0));
   } 

if (linhu!=0 && louthu==0)										/* parabolic defocusing top wall */
    {
      DIV1hu=linhu*linhu;
      p2parahu=2.0*(sqrt(DIV1hu+(h1u*h1u))-sqrt(DIV1hu));
      h2u=sqrt(p2parahu*(l+linhu+p2parahu/4.0));
      pbhu=-(linhu+p2parahu/4.0);
      pahu=-1.0/p2parahu;
      dyhu=uwallthick*sin(atan(-h2u*2*pahu));
      dzhu=uwallthick*cos(atan(-h2u*2*pahu));
      DIV1huwt=(linhu+l-dzhu)*(linhu+l-dzhu);
      p2parahuwt=2.0*(sqrt(DIV1huwt+((h2u+dyhu)*(h2u+dyhu)))-sqrt(DIV1huwt));
      h1uwt=sqrt(p2parahuwt*(linhu+p2parahuwt/4.0));
      h2uwt=sqrt(p2parahuwt*(linhu+l+p2parahuwt/4.0));
      pbhuwt=-(linhu+p2parahuwt/4.0);
      pahuwt=-1.0/p2parahuwt;
    }

if(linhu==0 && louthu==0){
    h1uwt=h1u+uwallthick/(cos(atan((h2u-h1u)/l)));
    h2uwt=h2u+uwallthick/(cos(atan((h2u-h1u)/l)));
    }


													/* Calculation of curve-parameters for the bottom wall - negative y-axis - analog right wall */
lbhd=linhd+l+louthd;
z0hd=lbhd/2.0-l-louthd;

if (linhd != 0 && louthd !=0)										/* elliptic bottom wall */
   {
    u1hd=sqrt((linhd*linhd)+(h1d*h1d));
    u2hd=sqrt((h1d*h1d)+((l+louthd)*(l+louthd)));
    ahd=(u1hd+u2hd)/2.0;
    bhd=sqrt(ahd*ahd-(lbhd*lbhd/4.0));
    a2hd=ahd*ahd;
    b2hd=bhd*bhd;
    DIV1hd=sqrt(1.0-((lbhd/2.0-louthd)*(lbhd/2.0-louthd)/a2hd));
    h2d=bhd*DIV1hd;
    if(l<lbhd/2-louthd){
        dyhd=dwallthick*sin(atan(a2hd*h1d/(b2hd*z0hd)));
        dzhd=dwallthick*cos(atan(a2hd*h1d/(b2hd*z0hd)));
    	u1hdwt=sqrt(((linhd+dzhd)*(linhd+dzhd))+((h1d+dyhd)*(h1d+dyhd)));	
    	u2hdwt=sqrt(((h1d+dyhd)*(h1d+dyhd))+((l+louthd-dzhd)*(l+louthd-dzhd)));
    	ahdwt=(u1hdwt+u2hdwt)/2.0;
    	a2hdwt=ahdwt*ahdwt;
    	bhdwt=sqrt(a2hdwt-(lbhd*lbhd/4.0));
    	b2hdwt=bhdwt*bhdwt;
    	h2dwt=bhdwt*sqrt(1.0-((lbhd/2.0-louthd)*(lbhd/2.0-louthd)/a2hdwt));
        h1dwt=bhdwt*sqrt(1.0-((lbhd/2.0-louthd-l)*(lbhd/2.0-louthd-l)/a2hdwt));
    }else{
	bhdwt=bhd+dwallthick;
        b2hdwt=bhdwt*bhdwt;
        ahdwt=sqrt(b2hdwt+(lbhd*lbhd/4.0));
        a2hdwt=b2hdwt+(lbhd*lbhd/4.0);
        h2dwt=bhdwt*sqrt(1.0-((lbhd/2.0-louthd)*(lbhd/2.0-louthd)/a2hdwt));
        h1dwt=bhdwt*sqrt(1.0-((lbhd/2.0-linhd)*(lbhd/2.0-linhd)/a2hdwt));
    }
  }

if(linhd==0 && louthd!=0)										/* parabolic focusing bottom wall */
   {	
    DIV1hd=(l+louthd)*(l+louthd);
    p2parahd=2.0*(sqrt(DIV1hd+(h1d*h1d))-sqrt(DIV1hd));
    h2d=sqrt(p2parahd*(louthd+p2parahd/4.0));
    pbhd=l+louthd+p2parahd/4.0;
    pahd=1.0/p2parahd;
    dyhd=dwallthick*sin(atan(h1d*2*pahd));
    dzhd=dwallthick*cos(atan(h1d*2*pahd));
    DIV1hdwt=(l+louthd-dzhd)*(l+louthd-dzhd);
    p2parahdwt=2.0*(sqrt(DIV1hdwt+((h1d+dyhd)*(h1d+dyhd)))-sqrt(DIV1hdwt));
    pbhdwt=l+louthd+p2parahdwt/4.0;
    pahdwt=1.0/p2parahdwt;
    h2dwt=sqrt(p2parahdwt*(louthd+p2parahdwt/4.0));
    h1dwt=sqrt(p2parahdwt*(louthd+l+p2parahdwt/4.0));
   }

if (linhd!=0 && louthd==0)										/* parabolic defocusing bottom wall */
    {
      DIV1hd=linhd*linhd;
      p2parahd=2.0*(sqrt(DIV1hd+(h1d*h1d))-sqrt(DIV1hd));
      h2d=sqrt(p2parahd*(l+linhd+p2parahd/4.0));
      pbhd=-(linhd+p2parahd/4.0);
      pahd=-1.0/p2parahd;
      dyhd=dwallthick*sin(atan(-h2d*2*pahd));
      dzhd=dwallthick*cos(atan(-h2d*2*pahd));
      DIV1hdwt=(linhd+l-dzhd)*(linhd+l-dzhd);
      p2parahdwt=2.0*(sqrt(DIV1hdwt+((h2d+dyhd)*(h2d+dyhd)))-sqrt(DIV1hdwt));
      h1dwt=sqrt(p2parahdwt*(linhd+p2parahdwt/4.0));
      h2dwt=sqrt(p2parahdwt*(linhd+l+p2parahdwt/4.0));
      pbhdwt=-(linhd+p2parahdwt/4.0);
      pahdwt=-1.0/p2parahdwt;
    }

if(linhd==0 && louthd==0){
    h1dwt=h1d+dwallthick/(cos(atan((h2d-h1d)/l)));
    h2dwt=h2d+dwallthick/(cos(atan((h2d-h1d)/l)));
    }



mru1=(h1uwt-h1u)/(w1r-w1rwt);										/* calculation for entrance and exit absorbing mask for the right upper corner*/
nru1=h1u-mru1*(-w1r);
mru2=(h2uwt-h2u)/(w2r-w2rwt);
nru2=h2u-mru2*(-w2r);

mrd1=(-h1dwt+h1d)/(w1r-w1rwt);										/* calculation for entrance and exit absorbing mask for the right lower corner*/
nrd1=-h1d-mrd1*(-w1r);
mrd2=(-h2dwt+h2d)/(w2r-w2rwt);
nrd2=-h2d-mrd2*(-w2r);

mlu1=(h1uwt-h1u)/(-w1l+w1lwt);										/* calculation for entrance and exit absorbing mask for the left upper corner*/
nlu1=h1u-mlu1*w1l;
mlu2=(h2uwt-h2u)/(-w2l+w2lwt);
nlu2=h2u-mlu2*w2l;

mld1=(-h1dwt+h1d)/(-w1l+w1lwt);										/* calculation for entrance and exit absorbing mask for the left lower corner*/
nld1=-h1d-mld1*w1l;
mld2=(-h2dwt+h2d)/(-w2l+w2lwt);
nld2=-h2d-mld2*w2l;

%}


TRACE
%{
  
  int i;

  PROP_Z0;												/* Propagate neutron to guide entrance. */ 


  if(x <= -w1r && x >= -w1rwt && y <= mru1*x+nru1 && y>= mrd1*x+nrd1 && mxr!=-1 && mxrOW!=-1){		/* absorbing the neutron if it hit the RIGHT entrance wall and the wall is not transparent*/
      ABSORB;
    }
  if(x >= w1l && x <=  w1lwt && y <= mlu1*x+nlu1 && y>= mld1*x+nld1 && mxl!=-1 && mxlOW!=-1 ){ 		/* absorbing the neutron if it hit the LEFT entrance wall and the wall is not transparent*/
      ABSORB;
     }
  if(y<=-h1d && y >=-h1dwt && x <= (y-nld1)/mld1 && x>= (y-nrd1)/mrd1 && myd!=-1 && mydOW!=-1){		/* absorbing the neutron if it hit the BOTTOM entrance wall and the wall is not transparent*/
      ABSORB;
     }
  if(y>=h1u && y <= h1uwt && x <= (y-nlu1)/mlu1 && x>= (y-nru1)/mru1 && myu!=-1 && myuOW!=-1){		/* absorbing the neutron if it hit the TOP entrance wall and the wall is not transparent*/
    ABSORB; 
    }


do{										/* start the propagation loop inside the guide */
  t1=(l-z)/vz;									/* needed time to pass the guide (or rest of the guide without any interaction)*/

  if(loutwr==0 && linwr==0)							/* for a straight right wall */
  {
  m=(-w2r+w1r)/l;
  t2w1r=(m*z+(-w1r)-x)/(vx-m*vz);						/* time untill next interaction with this wall (INNER wall)*/
  t2w1rwt=(m*z+(-w1rwt)-x)/(vx-m*vz);						/* time untill next interaction with this wall (OUTER wall)*/
  if(t2w1r<1e-15)								/* solving the precision problem for the intersection times given by double variable type, to small times (<1e-15) gives scattering events */
    t2w1r=t1+2.0;								/* at the same postion again (time is set to zero). this results in a sign change in the velocity components, which results in a wall tunneling.*/
  if(t2w1rwt<1e-15)								/* see comments above*/
    t2w1rwt=t1+2.0; 
  }


  if(loutwr!=0 && linwr!=0)							/* for a elliptic right wall */
  { 										/* solving the elliptic equation in respect to z and the straight neutron trajectoty, only two z values possible! */
    m=vx/vz;									/* m parameter of the neutron trajectory*/
    n=-m*z+x;									/* n parameter of the neutron trajectory */
    pwr=2.0*(a2wr*m*n+b2wr*z0wr)/(a2wr*m*m+b2wr);				/* p parameter of quadratic equation for calulation the z component of the intersection point with respect to the neutron trajectory (INNER wall)*/
    qwr=(a2wr*n*n+b2wr*z0wr*z0wr-a2wr*b2wr)/(a2wr*m*m+b2wr);			/* q parameter of quadratic equation for calulation the z component of the intersection point with respect to the neutron trajectory (INNER wall)*/ 
    if ((pwr*pwr/4.0)-qwr<0){							
    t2w1r=t1+2.0;								/* if the neutron never touch the ellipse the time is set to be bigger than the time (t1) needed to pass the component */
    }else{
        zw1r=-pwr/2.0+sqrt((pwr*pwr/4.0)-qwr);					/* first solution for z (INNER wall)*/
        zw2r=-pwr/2.0-sqrt((pwr*pwr/4.0)-qwr);					/* second solution for z (INNER wall)*/
        t2w1r=(zw1r-z)/vz;							/* interaction time for first z value (INNER wall)*/
        t2w2r=(zw2r-z)/vz;							/* interactime time for second z value (INNER wall)*/
        if(t2w1r<1e-15)								/* solving the precision problem for the intersection times given by double variable type, to small times (<1e-15) gives scattering events */
        t2w1r=t1+2.0;								/* at the same postion again (time is set to zero). this results in a sign change in the velocity components, which results in a wall tunneling.*/
        if(t2w2r<1e-15)								/* see comments above*/
        t2w2r=t1+2.0;
        if(t2w2r<t2w1r)								/* choosing the smaller positive time solution (INNER wall)*/
        t2w1r=t2w2r;
        xintersecwr=m*(vz*t2w1r+z)+n;						/* crosscheck of the x-coordinate of the intersection point */
        if (xintersecwr>0){							/* for the right wall x-coordinate of the intersection point have to be negative */
        t2w1r=t1+2.0;}								/* if this is not the case the time is set to t1+2.0 (time point behind the component) */
        }
    pwrwt=2.0*(a2wrwt*m*n+b2wrwt*z0wr)/(a2wrwt*m*m+b2wrwt);			/* p parameter of quadratic equation for calulation the z component of the intersection point with respect to the neutron trajectory (OUTER wall)*/
    qwrwt=(a2wrwt*n*n+b2wrwt*z0wr*z0wr-a2wrwt*b2wrwt)/(a2wrwt*m*m+b2wrwt);	/* q parameter of quadratic equation for calulation the z component of the intersection point with respect to the neutron trajectory (OUTER wall)*/
    if ((pwrwt*pwrwt/4.0)-qwrwt<0){							
    t2w1rwt=t1+2.0;								/* if the neutron never touch the ellipse the time is set bigger than need to pass the component */
    }else{
        zw1rwt=-pwrwt/2.0+sqrt((pwrwt*pwrwt/4.0)-qwrwt);			/* first solution for z (OUTER wall) */
        zw2rwt=-pwrwt/2.0-sqrt((pwrwt*pwrwt/4.0)-qwrwt);			/* second solution for z (OUTER wall)*/
        t2w1rwt=(zw1rwt-z)/vz;							/* interaction time for first z value (OUTER wall)*/
        t2w2rwt=(zw2rwt-z)/vz;							/* interactime time for second z value (OUTER wall)*/
        if(t2w1rwt<1e-15)							/* solving the precision problem for the intersection times given by double variable type, to small times (<1e-15) gives scattering events */
        t2w1rwt=t1+2.0;								/* at the same postion again (time is set to zero). this results in a sign change in the velocity components, which results in a wall tunneling.*/
        if(t2w2rwt<1e-15)							/* see comments above*/
        t2w2rwt=t1+2.0;
        if(t2w2rwt<t2w1rwt)							/* choosing the smaller positive time solution (OUTER wall)*/
        t2w1rwt=t2w2rwt;
        xintersecwrwt=m*(vz*t2w1rwt+z)+n;					/* crosscheck of the x-coordinate of the intersection point */
        if (xintersecwrwt>0){							/* for the right wall x-coordinate of the intersection point have to be negative */
        t2w1rwt=t1+2.0;}							/* if this is not the case the time is set to t1+2.0 (time point behind the component) */
        }
    }

  if((loutwr!=0 && linwr==0)|| (loutwr==0 && linwr!=0))				/* for a parabolic focusing or defocusing right wall*/
  {										/* solving the elliptic equation in respect to the straight neutron trajectoty to get the interaction time*/
    m=vx/vz;									/* m parameter of the neutron trajectory*/
    n=-m*z+x;									/* n parameter of the neutron trajectory */
    pwr=(2.0*m*n*pawr+1.0)/(pawr*m*m);						/* p parameter of quadratic equation (INNER wall)*/
    qwr=n*n/(m*m)-pbwr/(pawr*m*m);						/* q parameter of quadratic equation (INNER wall)*/
    if(qwr>0 && qwr>(pwr*pwr/4)){						/* in the very special case of no intersection the quadratic equation has no solution (negative square root) the time is set to t1+2.0 */
    t2w1r=t1+2.0;
    }else{
    if(vx==0)									/* in the special case of vx = 0 is x a constant */
    {
      if(x<0){									/* only neutron with a negativ x-component can hit the RIGHT wall (INNER wall)*/
      t2w1r=(pbwr-pawr*x*x-z)/vz;
      }else{
      t2w1r=t1+2.0; 								/* the time solution for neutron with a positive x component is set to a time long behind the exit of the guide */
										/* (means will not  scatter with the right wall)*/
      }
    }else{									/* if vx is not zero and x is a real variable*/
      zw1r=-pwr/2.0+sqrt(pwr*pwr/4.0-qwr);					/* first z-solution for intersection (INNER wall)*/
      zw2r=-pwr/2.0-sqrt(pwr*pwr/4.0-qwr);					/* second z-solution for intersection (INNER wall)*/
      t2w1r=(zw1r-z)/vz;							/* first time solution (INNER wall)*/
      t2w2r=(zw2r-z)/vz;							/* second time solution (INNER wall)*/
      if(t2w1r<1e-15)								/* solving the precision problem for the intersection times given by double variable type, to small times (<1e-15) gives scattering events */
      t2w1r=t1+2.0;								/* at the same postion again (time is set to zero). this results in a sign change in the velocity components, which results in a wall tunneling.*/
      if(t2w2r<1e-15)								/* see comments above*/
      t2w2r=t1+2.0;
      if(t2w2r<t2w1r)								/* choosing the smaller positive time solution (INNER wall)*/	
      t2w1r=t2w2r;
    }
    xintersecwr=m*(vz*t2w1r+z)+n;						/* crosscheck of the x-coordinate of the intersection point */
    if (xintersecwr>0){								/* for the right wall the x-coordinate of the intersection point have to be negative */
    t2w1r=t1+2.0;}								/* if this is not the case the time is set to t1+2.0 (time point behind the component) */
    }
    pwrwt=(2.0*m*n*pawrwt+1.0)/(pawrwt*m*m);					/* p parameter of quadratic equation (OUTER wall)*/
    qwrwt=n*n/(m*m)-pbwrwt/(pawrwt*m*m);					/* q parameter of quadratic equation (OUTER wall)*/
    if(qwrwt>0 && qwrwt>(pwrwt*pwrwt/4)){					/* in the very special case of no intersection the quadratic equation has no solution (negative square root) and the time is set to t1+2.0 */
    t2w1rwt=t1+2.0;
    }else{
    if(vx==0)									/* in the special case of vx = 0 is x a constant */
    {
      if(x<0){
      t2w1rwt=(pbwrwt-pawrwt*x*x-z)/vz;						/* only neutron with a negativ x-component can hit the RIGHT wall (OUTER wall)*/
      }else{
      t2w1rwt=t1+2.0; 
      }
    }else{									/* if vx is not zero */
      zw1rwt=-pwrwt/2.0+sqrt(pwrwt*pwrwt/4.0-qwrwt);				/* first z-solution for intersection (OUTER wall)*/
      zw2rwt=-pwrwt/2.0-sqrt(pwrwt*pwrwt/4.0-qwrwt);				/* second z-solution for intersection (OUTER wall)*/
      t2w1rwt=(zw1rwt-z)/vz;							/* first time solution (OUTER wall)*/
      t2w2rwt=(zw2rwt-z)/vz;							/* second time solution (OUTER wall)*/
      if(t2w1rwt<1e-15)								/* solving the precision problem for the intersection times given by double variable type, to small times (<1e-15) gives scattering events */
      t2w1rwt=t1+2.0;								/* at the same postion again (time is set to zero). this results in a sign change in the velocity components, which results in a wall tunneling.*/
      if(t2w2rwt<1e-15)								/* see comments above*/
      t2w2rwt=t1+2.0;
      if(t2w2rwt<t2w1rwt)							/* choosing the smaller positive time solution (OUTER wall)*/
      t2w1rwt=t2w2rwt;
    }
    xintersecwrwt=m*(vz*t2w1rwt+z)+n;						/* crosscheck of the x-coordinate of the intersection point */
    if (xintersecwrwt>0){							/* for the right wall the x-coordinate of the intersection point have to be negative */
    t2w1rwt=t1+2.0;}								/* if this is not the case the time is set to t1+2.0 (time point behind the component) */
  }
  }

  if(loutwl==0 && linwl==0)							/* comments for the left wall are same like for the right wall*/
  {
  m=(w2l-w1l)/l;
  t2w1l=(m*z+(w1l)-x)/(vx-m*vz);
  t2w1lwt=(m*z+(w1lwt)-x)/(vx-m*vz);
  if(t2w1l<1e-15)
  t2w1l=t1+2.0;
  if(t2w1lwt<1e-15)
  t2w1lwt=t1+2.0;
  };

  if(loutwl!=0 && linwl!=0)
  { 
    m=vx/vz;
    n=-m*z+x;
    pwl=2.0*(a2wl*m*n+b2wl*z0wl)/(a2wl*m*m+b2wl);
    qwl=(a2wl*n*n+b2wl*z0wl*z0wl-a2wl*b2wl)/(a2wl*m*m+b2wl);
    if((pwl*pwl/4.0)-qwl<0){
    t2w1l=t1+2.0;
    }else{
        zw1l=-pwl/2.0+sqrt((pwl*pwl/4.0)-qwl);
        zw2l=-pwl/2.0-sqrt((pwl*pwl/4.0)-qwl);
        t2w1l=(zw1l-z)/vz;
        t2w2l=(zw2l-z)/vz;
        if(t2w1l<1e-15)
        t2w1l=t1+2.0;
        if(t2w2l<1e-15)
        t2w2l=t1+2.0;
        if(t2w2l<t2w1l)
        t2w1l=t2w2l;
        xintersecwl=m*(vz*t2w1l+z)+n;
        if (xintersecwl<0){
        t2w1l=t1+2.0;}
        }
    pwlwt=2.0*(a2wlwt*m*n+b2wlwt*z0wl)/(a2wlwt*m*m+b2wlwt);
    qwlwt=(a2wlwt*n*n+b2wlwt*z0wl*z0wl-a2wlwt*b2wlwt)/(a2wlwt*m*m+b2wlwt);
    if((pwlwt*pwlwt/4.0)-qwlwt<0){
    t2w1lwt=t1+2.0;
    }else{
        zw1lwt=-pwlwt/2.0+sqrt((pwlwt*pwlwt/4.0)-qwlwt);
        zw2lwt=-pwlwt/2.0-sqrt((pwlwt*pwlwt/4.0)-qwlwt);
        t2w1lwt=(zw1lwt-z)/vz;
        t2w2lwt=(zw2lwt-z)/vz;
        if(t2w1lwt<1e-15)
        t2w1lwt=t1+2.0;
        if(t2w2lwt<1e-15)
        t2w2lwt=t1+2.0;
        if(t2w2lwt<t2w1lwt)
        t2w1lwt=t2w2lwt;
        xintersecwlwt=m*(vz*t2w1lwt+z)+n;
        if (xintersecwlwt<0){
        t2w1lwt=t1+2.0;}
        }
  }

  if((loutwl!=0 && linwl==0) || (loutwl==0 && linwl!=0))
  {
    m=vx/vz;
    n=-m*z+x;
    pwl=(2.0*m*n*pawl+1.0)/(pawl*m*m);
    qwl=n*n/(m*m)-pbwl/(pawl*m*m);
    if(qwl>0 && qwl>(pwl*pwl/4)){
    t2w1l=t1+2.0;
    }else{
    if(vx==0)
    {
      if(x>0){
      t2w1l=(pbwl-pawl*x*x-z)/vz;
      }else{
      t2w1l=t1+2.0;
      }
    }else{
      zw1l=-pwl/2.0+sqrt(pwl*pwl/4.0-qwl);
      zw2l=-pwl/2.0-sqrt(pwl*pwl/4.0-qwl);
      t2w1l=(zw1l-z)/vz;
      t2w2l=(zw2l-z)/vz;
      if(t2w1l<1e-15)
      t2w1l=t1+2.0;
      if(t2w2l<1e-15)
      t2w2l=t1+2.0;
      if(t2w2l<t2w1l)
      t2w1l=t2w2l;
    }
    }
    xintersecwl=m*(vz*t2w1l+z)+n;
    if (xintersecwl<0){
    t2w1l=t1+2.0;}
    pwlwt=(2.0*m*n*pawlwt+1.0)/(pawlwt*m*m);
    qwlwt=n*n/(m*m)-pbwlwt/(pawlwt*m*m);
    if(qwlwt>0 && qwlwt>(pwlwt*pwlwt/4)){
    t2w1lwt=t1+2.0;
    }else{
    if(vx==0)
    {
      if(x>0){
      t2w1lwt=(pbwlwt-pawlwt*x*x-z)/vz;
      }else{
      t2w1lwt=t1+2.0;
      }
    }else{
      zw1lwt=-pwlwt/2.0+sqrt(pwlwt*pwlwt/4.0-qwlwt);
      zw2lwt=-pwlwt/2.0-sqrt(pwlwt*pwlwt/4.0-qwlwt);
      t2w1lwt=(zw1lwt-z)/vz;
      t2w2lwt=(zw2lwt-z)/vz;
      if(t2w1lwt<1e-15)
      t2w1lwt=t1+2.0;
      if(t2w2lwt<1e-15)
      t2w2lwt=t1+2.0;
      if(t2w2lwt<t2w1lwt)
      t2w1lwt=t2w2lwt;
    }
    xintersecwlwt=m*(vz*t2w1lwt+z)+n;
    if (xintersecwlwt<0){
    t2w1lwt=t1+2.0;}
  }
  }

  if(louthu==0 && linhu==0)							/* comments for the top wall are same like for the right wall (only the ralation between left and right become relation between up and down)*/
  {
  m=(h2u-h1u)/l;
  t2h1u=(m*z+(h1u)-y)/(vy-m*vz);
  t2h1uwt=(m*z+(h1uwt)-y)/(vy-m*vz);
  if(t2h1u<1e-15)
  t2h1u=t1+2.0;
  if(t2h1uwt<1e-15)
  t2h1uwt=t1+2.0;
  }


  if(louthu!=0 && linhu!=0)
  { 
    m=vy/vz;
    n=-m*z+y;
    phu=2.0*(a2hu*m*n+b2hu*z0hu)/(a2hu*m*m+b2hu);
    qhu=(a2hu*n*n+b2hu*z0hu*z0hu-a2hu*b2hu)/(a2hu*m*m+b2hu);
    if((phu*phu/4.0)-qhu<0){
    t2h1u=t1+2.0;
    }else{
        zh1u=-phu/2.0+sqrt((phu*phu/4.0)-qhu);
        zh2u=-phu/2.0-sqrt((phu*phu/4.0)-qhu);
        t2h1u=(zh1u-z)/vz;
        t2h2u=(zh2u-z)/vz;
        if(t2h1u<1e-15)
        t2h1u=t1+2.0;
        if(t2h2u<1e-15)
        t2h2u=t1+2.0;
        if(t2h2u<t2h1u)
        t2h1u=t2h2u;
        yintersechu=m*(vz*t2h1u+z)+n;
        if (yintersechu<0){
        t2h1u=t1+2.0;}
        }
    phuwt=2.0*(a2huwt*m*n+b2huwt*z0hu)/(a2huwt*m*m+b2huwt);
    qhuwt=(a2huwt*n*n+b2huwt*z0hu*z0hu-a2huwt*b2huwt)/(a2huwt*m*m+b2huwt);
    if((phuwt*phuwt/4.0)-qhuwt<0){
    t2h1uwt=t1+2.0;
    }else{
        zh1uwt=-phuwt/2.0+sqrt((phuwt*phuwt/4.0)-qhuwt);
        zh2uwt=-phuwt/2.0-sqrt((phuwt*phuwt/4.0)-qhuwt);
        t2h1uwt=(zh1uwt-z)/vz;
        t2h2uwt=(zh2uwt-z)/vz;
        if(t2h1uwt<1e-15)
        t2h1uwt=t1+2.0;
        if(t2h2uwt<1e-15)
        t2h2uwt=t1+2.0;
        if(t2h2uwt<t2h1uwt)
        t2h1uwt=t2h2uwt;
        yintersechuwt=m*(vz*t2h1uwt+z)+n;
        if (yintersechuwt<0){
        t2h1uwt=t1+2.0;}
        }
  }


  if((louthu!=0 && linhu==0)|| (louthu==0 && linhu!=0))
  {
    m=vy/vz;
    n=-m*z+y;
    phu=(2.0*m*n*pahu+1.0)/(pahu*m*m);	
    qhu=n*n/(m*m)-pbhu/(pahu*m*m);
    if(qhu>0 && qhu>(phu*phu/4)){
    t2h1u=t1+2.0;
    }else{
    if(vy==0)
    {
      if(y>0){
      t2h1u=(pbhu-pahu*y*y-z)/vz;
      }else{
      t2h1u=t1+2.0;
      }
    }else{
      zh1u=-phu/2.0+sqrt(phu*phu/4.0-qhu);
      zh2u=-phu/2.0-sqrt(phu*phu/4.0-qhu);
      t2h1u=(zh1u-z)/vz;
      t2h2u=(zh2u-z)/vz;
      if(t2h1u<1e-15)
      t2h1u=t1+2.0;
      if(t2h2u<1e-15)
      t2h2u=t1+2.0;
      if(t2h2u<t2h1u)
      t2h1u=t2h2u;
    }
    yintersechu=m*(vz*t2h1u+z)+n;
    if (yintersechu<0){
    t2h1u=t1+2.0;}
  }
    phuwt=(2.0*m*n*pahuwt+1.0)/(pahuwt*m*m);						
    qhuwt=n*n/(m*m)-pbhuwt/(pahuwt*m*m);
    if(qhuwt>0 && qhuwt>(phuwt*phuwt/4)){
    t2h1uwt=t1+2.0;
    }else{
    if(vy==0)
    {
      if(y>0){
      t2h1uwt=(pbhuwt-pahuwt*y*y-z)/vz;
      }else{
      t2h1uwt=t1+2.0;
      }
    }else{
      zh1uwt=-phuwt/2.0+sqrt(phuwt*phuwt/4.0-qhuwt);
      zh2uwt=-phuwt/2.0-sqrt(phuwt*phuwt/4.0-qhuwt);
      t2h1uwt=(zh1uwt-z)/vz;
      t2h2uwt=(zh2uwt-z)/vz;
      if(t2h1uwt<1e-15)
      t2h1uwt=t1+2.0;
      if(t2h2uwt<1e-15)
      t2h2uwt=t1+2.0;
      if(t2h2uwt<t2h1uwt)
      t2h1uwt=t2h2uwt;
    }
    yintersechuwt=m*(vz*t2h1uwt+z)+n;
    if (yintersechuwt<0){
    t2h1uwt=t1+2.0;}
  }
  }

  if(louthd==0 && linhd==0)							/* comments for the down wall are same like for the right wall (only the relation between left and right become relation between up and down)*/	
  {
  m=(-h2d+h1d)/l;
  t2h1d=(m*z+(-h1d)-y)/(vy-m*vz);
  t2h1dwt=(m*z+(-h1dwt)-y)/(vy-m*vz);
  if(t2h1d<1e-15)
  t2h1d=t1+2.0;
  if(t2h1dwt<1e-15)
  t2h1dwt=t1+2.0;
  }

  if(louthd!=0 && linhd!=0)
  { 
    m=vy/vz;
    n=-m*z+y;
    phd=2.0*(a2hd*m*n+b2hd*z0hd)/(a2hd*m*m+b2hd);
    qhd=(a2hd*n*n+b2hd*z0hd*z0hd-a2hd*b2hd)/(a2hd*m*m+b2hd);
    if((phd*phd/4.0)-qhd<0){
    t2h1d=t1+2.0;
    }else{
        zh1d=-phd/2.0+sqrt((phd*phd/4.0)-qhd);
        zh2d=-phd/2.0-sqrt((phd*phd/4.0)-qhd);
        t2h1d=(zh1d-z)/vz;
        t2h2d=(zh2d-z)/vz;
        if(t2h1d<1e-15)
        t2h1d=t1+2.0;
        if(t2h2d<1e-15)
        t2h2d=t1+2.0;
        if(t2h2d<t2h1d)
        t2h1d=t2h2d;
        yintersechd=m*(vz*t2h1d+z)+n;
        if (yintersechd>0){
        t2h1d=t1+2.0;}
        }
    phdwt=2.0*(a2hdwt*m*n+b2hdwt*z0hd)/(a2hdwt*m*m+b2hdwt);
    qhdwt=(a2hdwt*n*n+b2hdwt*z0hd*z0hd-a2hdwt*b2hdwt)/(a2hdwt*m*m+b2hdwt);
    if((phdwt*phdwt/4.0)-qhdwt<0){
    t2h1dwt=t1+2.0;
    }else{
        zh1dwt=-phdwt/2.0+sqrt((phdwt*phdwt/4.0)-qhdwt);
        zh2dwt=-phdwt/2.0-sqrt((phdwt*phdwt/4.0)-qhdwt);
        t2h1dwt=(zh1dwt-z)/vz;
        t2h2dwt=(zh2dwt-z)/vz;
        if(t2h1dwt<1e-15)
        t2h1dwt=t1+2.0;
        if(t2h2dwt<1e-15)
        t2h2dwt=t1+2.0;
        if(t2h2dwt<t2h1dwt)
        t2h1dwt=t2h2dwt;
        yintersechdwt=m*(vz*t2h1dwt+z)+n;
        if (yintersechdwt>0){
        t2h1dwt=t1+2.0;}
        }
  }

  if((louthd!=0 && linhd==0)|| (louthd==0 && linhd!=0))
  {
    m=vy/vz;
    n=-m*z+y;
    phd=(2.0*m*n*pahd+1.0)/(pahd*m*m);
    qhd=n*n/(m*m)-pbhd/(pahd*m*m);
    if(qhd>0 && qhd>(phd*phd/4)){
    t2h1d=t1+2.0;
    }else{
    if(vy==0)
    {
      if(y<0){
      t2h1d=(pbhd-pahd*y*y-z)/vz;
      }else{
      t2h1d=t1+2.0;
      }
    }else{
      zh1d=-phd/2.0+sqrt(phd*phd/4.0-qhd);
      zh2d=-phd/2.0-sqrt(phd*phd/4.0-qhd);
      t2h1d=(zh1d-z)/vz;
      t2h2d=(zh2d-z)/vz;
      if(t2h1d<1e-15)
      t2h1d=t1+2.0;
      if(t2h2d<1e-15)
      t2h2d=t1+2.0;
      if(t2h2d<t2h1d)
      t2h1d=t2h2d;
    }
    yintersechd=m*(vz*t2h1d+z)+n;
    if (yintersechd>0){
    t2h1d=t1+2.0;}
  }
    phdwt=(2.0*m*n*pahdwt+1.0)/(pahdwt*m*m);
    qhdwt=n*n/(m*m)-pbhdwt/(pahdwt*m*m);
    if(qhdwt>0 && qhdwt>(phdwt*phdwt/4)){
    t2h1dwt=t1+2.0;
    }else{
    if(vy==0)
    {
      if(y<0){
      t2h1dwt=(pbhdwt-pahdwt*y*y-z)/vz;
      }else{
      t2h1dwt=t1+2.0;
      }
    }else{
      zh1dwt=-phdwt/2.0+sqrt(phdwt*phdwt/4.0-qhdwt);
      zh2dwt=-phdwt/2.0-sqrt(phdwt*phdwt/4.0-qhdwt);
      t2h1dwt=(zh1dwt-z)/vz;
      t2h2dwt=(zh2dwt-z)/vz;
      if(t2h1dwt<1e-15)
      t2h1dwt=t1+2.0;
      if(t2h2dwt<1e-15)
      t2h2dwt=t1+2.0;
      if(t2h2dwt<t2h1dwt)
      t2h1dwt=t2h2dwt;
    }
    yintersechdwt=m*(vz*t2h1dwt+z)+n;
    if (yintersechdwt>0){
    t2h1dwt=t1+2.0;}
  }
  }

													/* TEST OF THE INNER INTERSECTION - TIMES */
													/* possible interactions outside the guide have to be eliminated*/

if(t2w1r<t1+2.0){											/* test of RIGHT INNER wall interaction time*/
if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1r+z)-h1d;								/* calculation of the lower y-limit given by a linear bottom wall and the interaction time*/
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1r+z))*(z0hd+(vz*t2w1r+z)))/a2hd);				/* calculation of the lower y-limit given by a elliptic bottom wall and the interaction time*/
      }else{
      ylimitd=-sqrt(((vz*t2w1r+z)-pbhd)/-pahd);							/* calculation of the lower y-limit given by a parabolic bottom wall and the interaction time*/
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1r+z)+h1u;								/* calculation of the upper y-limit given by a linear top wall and the interaction time*/
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1r+z))*(z0hu+(vz*t2w1r+z)))/a2hu);				/* calculation of the upper y-limit given by a elliptic top wall and the interaction time*/
     }else{
     ylimitu=sqrt(((vz*t2w1r+z)-pbhu)/-pahu);							/* calculation of the upper y-limit given by a parabolic top wall and the interaction time*/
    }
    }
  ytest=vy*t2w1r+y;											/* calculation of the y coordinate of the neutron at the interaction time */
  if (ytest<ylimitd || ytest>ylimitu){									/* if the neutron is below or above the outer bottom or top wall of the guide NO interaction is possible!*/
  t2w1r=t1+2.0;
  }
  }


if(t2w1l<t1+2.0){											/* test of LEFT INNER wall interaction time - analog to right wall*/
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1l+z)-h1d;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1l+z))*(z0hd+(vz*t2w1l+z)))/a2hd);
      }else{
      ylimitd=-sqrt(((vz*t2w1l+z)-pbhd)/-pahd);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1l+z)+h1u;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1l+z))*(z0hu+(vz*t2w1l+z)))/a2hu);
     }else{
     ylimitu=sqrt(((vz*t2w1l+z)-pbhu)/-pahu);
    }
    }
  ytest=vy*t2w1l+y;
  if (ytest<ylimitd || ytest>ylimitu){
  t2w1l=t1+2.0;
  }
  }

if(t2h1u<t1+2.0){											/* test of TOP INNER wall interaction time - analog to right wall*/
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1u+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1u+z))*(z0wr+(vz*t2h1u+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1u+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1u+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1u+z))*(z0wl+(vz*t2h1u+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1u+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1u+x;
  if (xtest<xlimitr || xtest>xlimitl){
  t2h1u=t1+2.0;
  }
  }


if(t2h1d<t1+2.0){											/* test of BOTTOM INNER wall interaction time - analog to right wall*/
if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1d+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1d+z))*(z0wr+(vz*t2h1d+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1d+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1d+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1d+z))*(z0wl+(vz*t2h1d+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1d+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1d+x;
  if (xtest<xlimitr || xtest>xlimitl)
  {
  t2h1d=t1+2.0;
  }
  }


													/* TEST OF THE OUTER INTERSECTION - TIMES */

if(t2w1rwt<t1+2.0){											/* test of RIGHT OUTER wall interaction time - analog inner wall*/
if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1rwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1rwt+z))*(z0hd+(vz*t2w1rwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1rwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1rwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1rwt+z))*(z0hu+(vz*t2w1rwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1rwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1rwt+y;
  if (ytest<ylimitd || ytest>ylimitu){
  t2w1rwt=t1+2.0;
  }
  }


if(t2w1lwt<t1+2.0){											/* test of LEFT OUTER wall interaction time - analog inner wall*/
if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1lwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1lwt+z))*(z0hd+(vz*t2w1lwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1lwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1lwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1lwt+z))*(z0hu+(vz*t2w1lwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1lwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1lwt+y;
  if (ytest<ylimitd || ytest>ylimitu){
  t2w1lwt=t1+2.0;
  }
  }

if(t2h1uwt<t1+2.0){											/* test of TOP OUTER wall interaction time - analog inner wall*/
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1uwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1uwt+z))*(z0wr+(vz*t2h1uwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1uwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1uwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1uwt+z))*(z0wl+(vz*t2h1uwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1uwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1uwt+x;
  if (xtest<xlimitr || xtest>xlimitl){
  t2h1uwt=t1+2.0;
  }
  }


if(t2h1dwt<t1+2.0){											/* test of BOTTOM OUTER wall interaction time - analog inner wall*/
if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1dwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1dwt+z))*(z0wr+(vz*t2h1dwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1dwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1dwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1dwt+z))*(z0wl+(vz*t2h1dwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1dwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1dwt+x;
  if (xtest<xlimitr || xtest>xlimitl)
  {
  t2h1dwt=t1+2.0;
  }
  }


/* printf("t1 = %f \n",t1);
printf("t2w1r = %f \n",t2w1r);
printf("t2w1l = %f \n",t2w1l);
printf("t2h1u = %f \n",t2h1u);
printf("t2h1d = %f \n",t2h1d);
*/
/*
 printf("t1 = %f \n",t1);
printf("t2w1rwt = %f \n",t2w1rwt);
printf("t2w1lwt = %f \n",t2w1lwt);
printf("t2h1uwt = %f \n",t2h1uwt);
printf("t2h1dwt = %f \n",t2h1dwt);
*/
/*
printf("vx = %f \n",vx);
printf("x = %f \n",x); */

/*printf("mxr = %f \n", mxr);
printf("mxl = %f \n", mxl);
printf("myu = %f \n", myu);
printf("myd = %f \n", myd);*/

																				/* which wall is hit first? which geometry? */

  if (t1 < t2w1r && t1 < t2w1l && t1 < t2h1u && t1 < t2h1d && t1 < t2w1rwt && t1 < t2w1lwt && t1 < t2h1uwt && t1 < t2h1dwt){					/* neutron pass the guide without interaction*/
  i=1;
  }

													/* neutron interacts with the INNER elliptic right wall and this wall is NOT transparent*/

  if (mxr!=-1 && t2w1r > 0 && t2w1r < t1 && t2w1r<t2w1l && t2w1r < t2h1u && t2w1r < t2h1d && t2w1r < t2w1rwt && t2w1r<t2w1lwt && t2w1r < t2h1uwt && t2w1r < t2h1dwt && ((linwr!=0) && (loutwr!=0) )) {	
    if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1r+z)-h1d;								/* calculation of the lower y-limit given by a linear INNER bottom wall and the interaction time*/
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1r+z))*(z0hd+(vz*t2w1r+z)))/a2hd);				/* calculation of the lower y-limit given by a elliptic INNER bottom wall and the interaction time*/
      }else{
      ylimitd=-sqrt(((vz*t2w1r+z)-pbhd)/-pahd);								/* calculation of the lower y-limit given by a parabolic INNER bottom wall and the interaction time*/
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1r+z)+h1u;								/* calculation of the upper y-limit given by a linear INNER top wall and the interaction time*/
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1r+z))*(z0hu+(vz*t2w1r+z)))/a2hu);				/* calculation of the upper y-limit given by a elliptic INNER top wall and the interaction time*/
     }else{
     ylimitu=sqrt(((vz*t2w1r+z)-pbhu)/-pahu);								/* calculation of the upper y-limit given by a parabolic INNER top wall and the interaction time*/
    }
    }
  ytest=vy*t2w1r+y;											/* calculation of y-coordinate of the neutron at the interaction time*/
  xlimitr=bwr*sqrt(1-((vz*t2w1r+z)+z0wr)*((vz*t2w1r+z)+z0wr)/a2wr);					/* calculation of right INNER limit given by a elliptic shape of the wall*/
  xtest=vx*t2w1r+x;											/* calculation of x-coordinate of the neutron at the interaction time*/
  if (ytest>=ylimitd && ytest<=ylimitu && x>=-xlimitr && mxr > 0){					/* only if the neutron is realy inside the guide (testet by the limits) and the wall is reflecting*/
  i=2;													/* the neutron will be reflected*/
  }else{
  i=18;													/* otherwise it will be absorbed*/
  }
  }

													/* neutron interacts with the parabolic right INNER wall - comments are analog to inner elliptic right wall */

  if (mxr!=-1 && t2w1r > 0 && t2w1r < t1 && t2w1r<t2w1l && t2w1r < t2h1u && t2w1r < t2h1d && t2w1r < t2w1rwt && t2w1r<t2w1lwt && t2w1r < t2h1uwt && t2w1r < t2h1dwt && ((loutwr!=0 && linwr==0) || (loutwr==0 && linwr!=0))){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1r+z)-h1d;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1r+z))*(z0hd+(vz*t2w1r+z)))/a2hd);
      }else{
      ylimitd=-sqrt(((vz*t2w1r+z)-pbhd)/-pahd);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1r+z)+h1u;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1r+z))*(z0hu+(vz*t2w1r+z)))/a2hu);
     }else{
     ylimitu=sqrt(((vz*t2w1r+z)-pbhu)/-pahu);
    }
    }
  ytest=vy*t2w1r+y;
  xtest=vx*t2w1r+x;
  xlimitr=-sqrt(pbwr-(vz*t2w1r+z)/pawr);
  if (ytest>ylimitd && ytest<ylimitu && x>=xlimitr && mxr > 0){
  i=3;
  }else{
  i=18;}
  }

													/* neutron interacts with the straight right INNER wall - comments are analog to inner elliptic right wall*/

  if (mxr!=-1 && t2w1r > 0 && t2w1r < t1 && t2w1r<t2w1l && t2w1r < t2h1u && t2w1r < t2h1d && t2w1r < t2w1rwt && t2w1r<t2w1lwt && t2w1r < t2h1uwt && t2w1r < t2h1dwt && (loutwr==0 && linwr==0)){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1r+z)-h1d;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1r+z))*(z0hd+(vz*t2w1r+z)))/a2hd);
      }else{
      ylimitd=-sqrt(((vz*t2w1r+z)-pbhd)/-pahd);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1r+z)+h1u;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1r+z))*(z0hu+(vz*t2w1r+z)))/a2hu);
     }else{
     ylimitu=sqrt(((vz*t2w1r+z)-pbhu)/-pahu);
    }
    }
  ytest=vy*t2w1r+y;
  xtest=vx*t2w1r+x;
  xlimitr=-(w2r+w1r)*(vz*t2w1r+z)/l-w1r;
  if (ytest>=ylimitd && ytest<=ylimitu && x>=xlimitr && mxr>0){
  i=4;
  }else{
  i=18;
  }
  }

													/* neutron interacts with the elliptic left INNER wall - comments are analog to inner elliptic right wall*/

  if (mxl!=-1 && t2w1l > 0 && t2w1l < t1 && t2w1l<t2w1r && t2w1l < t2h1u && t2w1l < t2h1d && t2w1l < t2w1lwt && t2w1l<t2w1rwt && t2w1l < t2h1uwt && t2w1l < t2h1dwt && ((linwl!=0) && (loutwl!=0) )){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1l+z)-h1d;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1l+z))*(z0hd+(vz*t2w1l+z)))/a2hd);
      }else{
      ylimitd=-sqrt(((vz*t2w1l+z)-pbhd)/-pahd);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1l+z)+h1u;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1l+z))*(z0hu+(vz*t2w1l+z)))/a2hu);
     }else{
     ylimitu=sqrt(((vz*t2w1l+z)-pbhu)/-pahu);
    }
    }
  ytest=vy*t2w1l+y;
  xlimitl=bwl*sqrt(1-((vz*t2w1l+z)+z0wl)*((vz*t2w1l+z)+z0wl)/a2wl);
  if (ytest>=ylimitd && ytest<=ylimitu && x<=xlimitl && mxl>0){
  i=5;
  }else{
  i=19;
  }
  }

													/* neutron interacts with the parabolic left INNER wall - comments are analog to inner elliptic right wall*/

  if (mxl!=-1 && t2w1l > 0 && t2w1l < t1 && t2w1l<t2w1r && t2w1l < t2h1u && t2w1l < t2h1d && t2w1l < t2w1lwt && t2w1l<t2w1rwt && t2w1l < t2h1uwt && t2w1l < t2h1dwt && ((loutwl!=0 && linwl==0) || (loutwl==0 && linwl!=0))){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1l+z)-h1d;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1l+z))*(z0hd+(vz*t2w1l+z)))/a2hd);
      }else{
      ylimitd=-sqrt(((vz*t2w1l+z)-pbhd)/-pahd);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1l+z)+h1u;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1l+z))*(z0hu+(vz*t2w1l+z)))/a2hu);
     }else{
     ylimitu=sqrt(((vz*t2w1l+z)-pbhu)/-pahu);
    }
    }
  ytest=vy*t2w1l+y;
  xlimitl=sqrt(pbwl-(vz*t2w1l+z)/pawl);
  if (ytest>=ylimitd && ytest<=ylimitu && x<=xlimitl && mxl>0){
  i=6;
  }else{
  i=19;
  }
  }

														/* neutron interacts with the straight left INNER wall - comments are analog to inner elliptic right wall*/

  if (mxl!=-1 && t2w1l > 0 && t2w1l < t1 && t2w1l<t2w1r && t2w1l < t2h1u && t2w1l < t2h1d && t2w1l < t2w1lwt && t2w1l < t2w1rwt && t2w1l < t2h1uwt && t2w1l < t2h1dwt && (loutwl==0 && linwl==0)){
 if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1l+z)-h1d;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhd*sqrt(1-((z0hd+(vz*t2w1l+z))*(z0hd+(vz*t2w1l+z)))/a2hd);
      }else{
      ylimitd=-sqrt(((vz*t2w1l+z)-pbhd)/-pahd);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1l+z)+h1u;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhu*sqrt(1-((z0hu+(vz*t2w1l+z))*(z0hu+(vz*t2w1l+z)))/a2hu);
     }else{
     ylimitu=sqrt(((vz*t2w1l+z)-pbhu)/-pahu);
    }
    }
  ytest=vy*t2w1l+y;
  xlimitl=(w2l-w1l)*(vz*t2w1l+z)/l+w1l;
  if (ytest>=ylimitd && ytest<=ylimitu && x<=xlimitl && mxl>0){
  i=7;
  }else{
  i=19;
  }
  }

														/* neutron interacts with the elliptic top INNER wall - comments are analog to inner elliptic right wall*/

  if (myu!=-1 && t2h1u > 0 && t2h1u < t1 && t2h1u < t2w1r && t2h1u < t2w1l && t2h1u < t2h1d && t2h1u < t2h1uwt && t2h1u < t2w1rwt && t2h1u < t2w1lwt && t2h1u < t2h1dwt && (louthu !=0 && linhu!=0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1u+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1u+z))*(z0wr+(vz*t2h1u+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1u+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1u+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1u+z))*(z0wl+(vz*t2h1u+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1u+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1u+x;
  ylimitu=bhu*sqrt(1-((vz*t2h1u+z)+z0hu)*((vz*t2h1u+z)+z0hu)/a2hu);
  if (xtest>=xlimitr && xtest<=xlimitl && y<=ylimitu && myu>0 ){
  i=8;
  }else{
  i=20;
  }
  }
	
														/* neutron interacts with the parabolic top INNER wall - comments are analog to inner elliptic right wall*/

  if (myu!=-1 && t2h1u > 0 && t2h1u < t1 && t2h1u < t2w1r && t2h1u < t2w1l && t2h1u < t2h1d && t2h1u < t2h1uwt && t2h1u < t2w1rwt && t2h1u < t2w1lwt && t2h1u < t2h1dwt && ((louthu!=0 && linhu==0) || (louthu==0 && linhu!=0))){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1u+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1u+z))*(z0wr+(vz*t2h1u+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1u+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1u+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1u+z))*(z0wl+(vz*t2h1u+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1u+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1u+x;
  ylimitu=sqrt(pbhu-(vz*t2h1u+z)/pahu);
  if (xtest>=xlimitr && xtest<=xlimitl && y<=ylimitu && myu>0){
  i=9;
  }else{
  i=20;
  }
  }


															/* neutron interacts with the straight top INNER wall - comments are analog to inner elliptic right wall*/

  if (myu!=-1 && t2h1u > 0 && t2h1u < t1 && t2h1u < t2w1r && t2h1u < t2w1l && t2h1u < t2h1d && t2h1u < t2h1uwt && t2h1u < t2w1rwt && t2h1u < t2w1lwt && t2h1u < t2h1dwt && (louthu == 0 && linhu == 0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1u+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1u+z))*(z0wr+(vz*t2h1u+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1u+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1u+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1u+z))*(z0wl+(vz*t2h1u+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1u+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1u+x;
  ylimitu=(h2u-h1u)*(vz*t2h1u+z)/l+h1u;
  if (xtest>=xlimitr && xtest<=xlimitl && y<=ylimitu && myu>0 ){
  i=10;
  }else{
  i=20;
  }
  }

															/* neutron interacts with the elliptic down INNER wall - comments are analog to inner elliptic right wall*/

  if (myd!=-1 && t2h1d > 0 && t2h1d < t1 && t2h1d < t2w1r && t2h1d < t2w1l && t2h1d < t2h1u  && t2h1d < t2h1dwt && t2h1d < t2w1rwt && t2h1d < t2w1lwt && t2h1d < t2h1uwt  && (louthd !=0 && linhd!=0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1d+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1d+z))*(z0wr+(vz*t2h1d+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1d+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1d+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1d+z))*(z0wl+(vz*t2h1d+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1d+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1d+x;
  ylimitd=-bhd*sqrt(1-((vz*t2h1d+z)+z0hd)*((vz*t2h1d+z)+z0hd)/a2hd);
  if (xtest>=xlimitr && xtest<=xlimitl && y>=ylimitd && myd>0){
  i=11;
  }else{
  i=21;
  }
  }

															/* neutron interacts with the parabolic down INNER wall - comments are analog to inner elliptic right wall*/

  if (myd!=-1 && t2h1d > 0 && t2h1d < t1 && t2h1d < t2w1r && t2h1d < t2w1l && t2h1d < t2h1u  && t2h1d < t2h1dwt && t2h1d < t2w1rwt && t2h1d < t2w1lwt && t2h1d < t2h1uwt  && ((louthd !=0 && linhd==0) || (louthd ==0 && linhd!=0))){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1d+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1d+z))*(z0wr+(vz*t2h1d+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1d+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1d+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1d+z))*(z0wl+(vz*t2h1d+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1d+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1d+x;
  ylimitd=-sqrt(pbhd-(vz*t2h1d+z)/pahd);
  if (xtest>=xlimitr && xtest<=xlimitl && y>=ylimitd && myd>0 ){
  i=12;
  }else{
  i=21;
  }
  }

															/* neutron interacts with the straight down INNER wall - comments are analog to inner elliptic right wall*/

  if (myd!=-1 && t2h1d > 0 && t2h1d < t1 && t2h1d < t2w1r && t2h1d < t2w1l && t2h1d < t2h1u  && t2h1d < t2h1dwt && t2h1d < t2w1rwt && t2h1d < t2w1lwt && t2h1d < t2h1uwt  && (louthd == 0 && linhd == 0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1d+z)-w1r;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwr*sqrt(1-((z0wr+(vz*t2h1d+z))*(z0wr+(vz*t2h1d+z)))/a2wr);
    }else{
    xlimitr=-sqrt(((vz*t2h1d+z)-pbwr)/-pawr);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1d+z)+w1l;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwl*sqrt(1-((z0wl+(vz*t2h1d+z))*(z0wl+(vz*t2h1d+z)))/a2wl);
      }else{
      xlimitl=sqrt(((vz*t2h1d+z)-pbwl)/-pawl);
      }
   }
  xtest=vx*t2h1d+x;
  ylimitd=(-h2d+h1d)*(vz*t2h1d+z)/l-h1d;
  if (xtest>=xlimitr && xtest<=xlimitl && y>=ylimitd && myd>0 ){
  i=13;
  }else{
  i=21;
  }
  }

   /* SAME INTERACTIONS BUT TRANSPERENT INNER WALLS - in these cases no limit check up is needed , because will neutrons will pass unchanged anyway*/

  if (mxr==-1 && t2w1r > 0 && t2w1r < t1 && t2w1r<t2w1l  && t2w1r < t2h1u && t2w1r < t2h1d && t2w1r < t2w1rwt && t2w1r<t2w1lwt  && t2w1r < t2h1uwt && t2w1r < t2h1dwt && ((linwr!=0) && (loutwr!=0) )) {
  i=14;
  }

  if (mxr==-1 && t2w1r > 0 && t2w1r < t1 && t2w1r<t2w1l  && t2w1r < t2h1u && t2w1r < t2h1d && t2w1r < t2w1rwt && t2w1r<t2w1lwt  && t2w1r < t2h1uwt && t2w1r < t2h1dwt && ((loutwr!=0 && linwr==0) || (loutwr==0 && linwr!=0))){
  i=14;
  }

  if (mxr==-1 && t2w1r > 0 && t2w1r < t1 && t2w1r<t2w1l  && t2w1r < t2h1u && t2w1r < t2h1d && t2w1r < t2w1rwt && t2w1r<t2w1lwt  && t2w1r < t2h1uwt && t2w1r < t2h1dwt && (loutwr==0 && linwr==0)){
  i=14;
  }

  if (mxl==-1 && t2w1l > 0 && t2w1l < t1 && ( t2w1l<t2w1r ) && t2w1l < t2h1u && t2w1l < t2h1d && t2w1l < t2w1lwt && t2w1l<t2w1rwt && t2w1l < t2h1uwt && t2w1l < t2h1dwt && ((linwl!=0) && (loutwl!=0) )){
  i=15;
  }

  if (mxl==-1 && t2w1l > 0 && t2w1l < t1 && ( t2w1l<t2w1r ) && t2w1l < t2h1u && t2w1l < t2h1d && t2w1l < t2w1lwt && t2w1l<t2w1rwt && t2w1l < t2h1uwt && t2w1l < t2h1dwt &&  ((loutwl!=0 && linwl==0) || (loutwl==0 && linwl!=0))){
  i=15;
  }

  if (mxl==-1 && t2w1l > 0 && t2w1l < t1 && ( t2w1l<t2w1r ) && t2w1l < t2h1u && t2w1l < t2h1d && t2w1l < t2w1lwt && t2w1l<t2w1rwt && t2w1l < t2h1uwt && t2w1l < t2h1dwt && (loutwl==0 && linwl==0)){
  i=15;
  }

  if (myu==-1 && t2h1u > 0 && t2h1u < t1 && t2h1u < t2w1r && t2h1u < t2w1l && t2h1u < t2h1d  && t2h1u < t2h1uwt && t2h1u < t2w1rwt && t2h1u < t2w1lwt && t2h1u < t2h1dwt  && (louthu !=0 && linhu!=0)){
  i=16;
  }

  if (myu==-1 && t2h1u > 0 && t2h1u < t1 && t2h1u < t2w1r && t2h1u < t2w1l && t2h1u < t2h1d  && t2h1u < t2h1uwt && t2h1u < t2w1rwt && t2h1u < t2w1lwt && t2h1u < t2h1dwt  && ((louthu!=0 && linhu==0) || (louthu==0 && linhu!=0))){
  i=16;
  }

  if (myu==-1 && t2h1u > 0 && t2h1u < t1 && t2h1u < t2w1r && t2h1u < t2w1l && t2h1u < t2h1d  && t2h1u < t2h1uwt && t2h1u < t2w1rwt && t2h1u < t2w1lwt && t2h1u < t2h1dwt  && (louthu == 0 && linhu == 0)){
  i=16;
  }

  if (myd==-1 && t2h1d > 0 && t2h1d < t1 && t2h1d < t2w1r && t2h1d < t2w1l && t2h1d < t2h1u  && t2h1d < t2h1dwt && t2h1d < t2w1rwt && t2h1d < t2w1lwt && t2h1d < t2h1uwt  && (louthd !=0 && linhd!=0)){
  i=17;
  }

  if (myd==-1 && t2h1d > 0 && t2h1d < t1 && t2h1d < t2w1r && t2h1d < t2w1l && t2h1d < t2h1u  && t2h1d < t2h1dwt && t2h1d < t2w1rwt && t2h1d < t2w1lwt && t2h1d < t2h1uwt  && ((louthd !=0 && linhd==0) || (louthd ==0 && linhd!=0))){
  i=17;
  }

  if (myd==-1 && t2h1d > 0 && t2h1d < t1 && t2h1d < t2w1r && t2h1d < t2w1l && t2h1d < t2h1u  && t2h1d < t2h1dwt && t2h1d < t2w1rwt && t2h1d < t2w1lwt && t2h1d < t2h1uwt  && (louthd == 0 && linhd == 0)){
  i=17;
  }

   
												/* EVERTHING AGAIN FOR THE OUTER WALLS */

												/* neutron interacts with the elliptic right OUTER wall - comments are analog to inner elliptic right wall*/

  if (mxrOW!=-1 && t2w1rwt > 0 && t2w1rwt < t1 && t2w1rwt < t2w1r && t2w1rwt <t2w1l && t2w1rwt < t2h1u && t2w1rwt < t2h1d &&  t2w1rwt <t2w1lwt && t2w1rwt < t2h1uwt && t2w1rwt < t2h1dwt && ((linwr!=0) && (loutwr!=0) )){	
   if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1rwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1rwt+z))*(z0hd+(vz*t2w1rwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1rwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1rwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1rwt+z))*(z0hu+(vz*t2w1rwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1rwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1rwt+y;
  xtest=vx*t2w1rwt+x;
  xlimitrwt=-bwrwt*sqrt(1-((vz*t2w1rwt+z)+z0wr)*((vz*t2w1rwt+z)+z0wr)/a2wrwt);
  if (ytest>=ylimitd && ytest<=ylimitu && xtest<=xlimitr && mxrOW > 0){
  i=22;
  }
  else{
  i=34;
  }
  }

													/* neutron interacts with the parabolic right OUTER wall - comments are analog to inner elliptic right wall*/

  if (mxrOW!=-1 && t2w1rwt > 0 && t2w1rwt < t1 && t2w1rwt < t2w1r && t2w1rwt <t2w1l && t2w1rwt < t2h1u && t2w1rwt < t2h1d &&  t2w1rwt <t2w1lwt && t2w1rwt < t2h1uwt && t2w1rwt < t2h1dwt && ((loutwr!=0 && linwr==0) || (loutwr==0 && linwr!=0))){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1rwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1rwt+z))*(z0hd+(vz*t2w1rwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1rwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1rwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1rwt+z))*(z0hu+(vz*t2w1rwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1rwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1rwt+y;
  xtest=vx*t2w1rwt+x;
  xlimitr=-sqrt((pbwrwt-(vz*t2w1rwt+z))/pawrwt);
  if (ytest>ylimitd && ytest<ylimitu && xtest<=xlimitr && mxrOW > 0){
  i=23; 
  }else{
  i=34;}
  }

													/* neutron interacts with the straight right OUTER wall - comments are analog to inner elliptic right wall */

  if (mxrOW!=-1 && t2w1rwt > 0 && t2w1rwt < t1 && t2w1rwt < t2w1r && t2w1rwt <t2w1l && t2w1rwt < t2h1u && t2w1rwt < t2h1d &&  t2w1rwt <t2w1lwt && t2w1rwt < t2h1uwt && t2w1rwt < t2h1dwt && (loutwr==0 && linwr==0)){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1rwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1rwt+z))*(z0hd+(vz*t2w1rwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1rwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1rwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1rwt+z))*(z0hu+(vz*t2w1rwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1rwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1rwt+y; 
  xtest=vx*t2w1rwt+x;
  xlimitr=-(-w2r+w1r)*(vz*t2w1rwt+z)/l-w1rwt;
  if(ytest>=ylimitd && ytest<=ylimitu && xtest<=xlimitr && mxrOW>0){
  i=24;
  }else{
  i=34;
  }
  }

													/* neutron interacts with the elliptic left OUTER wall - comments are analog to inner elliptic right wall*/

  if (mxlOW!=-1 && t2w1lwt > 0 && t2w1lwt < t1 && t2w1lwt < t2w1r && t2w1lwt <t2w1l && t2w1lwt < t2h1u && t2w1lwt < t2h1d &&  t2w1lwt <t2w1rwt && t2w1lwt < t2h1uwt && t2w1lwt < t2h1dwt && ((linwl!=0) && (loutwl!=0) )){
  if(linhd==0 && louthd==0)
   if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1lwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1lwt+z))*(z0hd+(vz*t2w1lwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1lwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1lwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1lwt+z))*(z0hu+(vz*t2w1lwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1lwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1lwt+y;
  xtest=vx*t2w1lwt+x;
  xlimitlwt=-bwlwt*sqrt(1-((vz*t2w1lwt+z)+z0wl)*((vz*t2w1lwt+z)+z0wl)/a2wlwt);
  if (ytest>=ylimitd && ytest<=ylimitu && xtest>=xlimitl && mxlOW>0){
  i=25;
  }else{
  i=35;
  }
  }

													/* neutron interacts with the parabolic left OUTER wall - comments are analog to inner elliptic right wall */

  if (mxlOW!=-1 && t2w1lwt > 0 && t2w1lwt < t1 && t2w1lwt < t2w1r && t2w1lwt <t2w1l && t2w1lwt < t2h1u && t2w1lwt < t2h1d &&  t2w1lwt <t2w1rwt && t2w1lwt < t2h1uwt && t2w1lwt < t2h1dwt && ((loutwl!=0 && linwl==0) || (loutwl==0 && linwl!=0))){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1lwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1lwt+z))*(z0hd+(vz*t2w1lwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1lwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1lwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1lwt+z))*(z0hu+(vz*t2w1lwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1lwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1lwt+y;
  xtest=vx*t2w1lwt+x;
  xlimitl=-sqrt((pbwlwt-(vz*t2w1lwt+z))/pawlwt);
  if (ytest>=ylimitd && ytest<=ylimitu && xtest>=xlimitl && mxlOW>0){
  i=26;
  }else{
  i=35;
  }
  }

														/* neutron interacts with the straight left OUTER wall - comments are analog to inner elliptic right wall*/

  if (mxlOW!=-1 && t2w1lwt > 0 && t2w1lwt < t1 && t2w1lwt < t2w1r && t2w1lwt <t2w1l && t2w1lwt < t2h1u && t2w1lwt < t2h1d &&  t2w1lwt <t2w1rwt && t2w1lwt < t2h1uwt && t2w1lwt < t2h1dwt && (loutwl==0 && linwl==0)){
  if(linhd==0 && louthd==0)
   {
      ylimitd=(-h2d+h1d)/l*(vz*t2w1lwt+z)-h1dwt;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ylimitd=-bhdwt*sqrt(1-((z0hd+(vz*t2w1lwt+z))*(z0hd+(vz*t2w1lwt+z)))/a2hdwt);
      }else{
      ylimitd=-sqrt(((vz*t2w1lwt+z)-pbhdwt)/-pahdwt);
      }
   }
   if(linhu==0 && louthu==0)
    {
     ylimitu=(h2u-h1u)/l*(vz*t2w1lwt+z)+h1uwt;
    }else{
     if(linhu!=0 && louthu!=0)
     {
     ylimitu=bhuwt*sqrt(1-((z0hu+(vz*t2w1lwt+z))*(z0hu+(vz*t2w1lwt+z)))/a2huwt);
     }else{
     ylimitu=sqrt(((vz*t2w1lwt+z)-pbhuwt)/-pahuwt);
    }
    }
  ytest=vy*t2w1lwt+y; 
  xtest=vx*t2w1lwt+x;
  xlimitl=(w2l-w1l)*(vz*t2w1lwt+z)/l+w1lwt;
  if (ytest>=ylimitd && ytest<=ylimitu && xtest>=xlimitl && mxlOW>0){
  i=27;
  }else{
  i=35;
  }
  }

															/* neutron interacts with the elliptic top OUTER wall - comments are analog to inner elliptic right wall*/

  if (myuOW!=-1 && t2h1uwt > 0 && t2h1uwt < t1 && t2h1uwt < t2w1r && t2h1uwt < t2w1l && t2h1uwt < t2h1d && t2h1uwt < t2h1u && t2h1uwt < t2w1rwt && t2h1uwt < t2w1lwt && t2h1uwt < t2h1dwt && (louthu !=0 && linhu!=0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1uwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1uwt+z))*(z0wr+(vz*t2h1uwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1uwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1uwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1uwt+z))*(z0wl+(vz*t2h1uwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1uwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1uwt+x;
  ytest=vy*t2w1lwt+y;
  ylimitu=bhuwt*sqrt(1-((vz*t2h1uwt+z)+z0hu)*((vz*t2h1uwt+z)+z0hu)/a2huwt);
  if (xtest>=xlimitr && xtest<=xlimitl && ytest>=ylimitu && myuOW>0 ){
  i=28;
  }else{
  i=36;
  }
  }
	
															/* neutron interacts with the parabolic top OUTER wall - comments are analog to inner elliptic right wall */

  if (myuOW!=-1 && t2h1uwt > 0 && t2h1uwt < t1 && t2h1uwt < t2w1r && t2h1uwt < t2w1l && t2h1uwt < t2h1d && t2h1uwt < t2h1u && t2h1uwt < t2w1rwt && t2h1uwt < t2w1lwt && t2h1uwt < t2h1dwt && ((louthu!=0 && linhu==0) || (louthu==0 && linhu!=0))){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1uwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1uwt+z))*(z0wr+(vz*t2h1uwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1uwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1uwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1uwt+z))*(z0wl+(vz*t2h1uwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1uwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1uwt+x;
  ytest=vy*t2w1lwt+y;
  ylimitu=sqrt((pbhuwt-(vz*t2h1uwt+z))/pahuwt);
  if (xtest>=xlimitr && xtest<=xlimitl && ytest>=ylimitu && myuOW>0){
  i=29;
  }else{
  i=36;
  }
  }


															/* neutron interacts with the straight top OUTER wall - comments are analog to inner elliptic right wall*/

  if (myuOW!=-1 && t2h1uwt > 0 && t2h1uwt < t1 && t2h1uwt < t2w1r && t2h1uwt < t2w1l && t2h1uwt < t2h1d && t2h1uwt < t2h1u && t2h1uwt < t2w1rwt && t2h1uwt < t2w1lwt && t2h1uwt < t2h1dwt && (louthu == 0 && linhu == 0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1uwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1uwt+z))*(z0wr+(vz*t2h1uwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1uwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1uwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1uwt+z))*(z0wl+(vz*t2h1uwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1uwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1uwt+x;
  ytest=vy*t2w1lwt+y;
  ylimitu=(h2u-h1u)*(vz*t2h1uwt+z)/l+h1uwt;
  if (xtest>=xlimitr && xtest<=xlimitl && ytest>=ylimitu && myuOW>0 ){
  i=30;
  }else{
  i=36;
  }
  }

															/* neutron interacts with the elliptic down OUTER wall - comments are analog to inner elliptic right wall*/

  if (mydOW!=-1 && t2h1dwt > 0 && t2h1dwt < t1 && t2h1dwt < t2w1r && t2h1dwt < t2w1l && t2h1dwt < t2h1u  && t2h1dwt < t2h1d && t2h1dwt < t2w1rwt && t2h1dwt < t2w1lwt && t2h1dwt < t2h1uwt  && (louthd !=0 && linhd!=0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1dwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1dwt+z))*(z0wr+(vz*t2h1dwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1dwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1dwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1dwt+z))*(z0wl+(vz*t2h1dwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1dwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1dwt+x;
  ytest=vy*t2w1lwt+y;
  ylimitd=-bhdwt*sqrt(1-((vz*t2h1dwt+z)+z0hd)*((vz*t2h1dwt+z)+z0hd)/a2hdwt);
  if (xtest>=xlimitr && xtest<=xlimitl && ytest<=ylimitd && mydOW>0){
  i=31;
  }else{
  i=37;
  }
  }

															/* neutron interacts with the parabolic down OUTER wall - comments are analog to inner elliptic right wall */

  if (mydOW!=-1 && t2h1dwt > 0 && t2h1dwt < t1 && t2h1dwt < t2w1r && t2h1dwt < t2w1l && t2h1dwt < t2h1u  && t2h1dwt < t2h1d && t2h1dwt < t2w1rwt && t2h1dwt < t2w1lwt && t2h1dwt < t2h1uwt  && ((louthd !=0 && linhd==0) || (louthd ==0 && linhd!=0))){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1dwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1dwt+z))*(z0wr+(vz*t2h1dwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1dwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
     xlimitl=(w2l-w1l)/l*(vz*t2h1dwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1dwt+z))*(z0wl+(vz*t2h1dwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1dwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1dwt+x;
  ytest=vy*t2w1lwt+y;
  ylimitd=-sqrt((pbhdwt-(vz*t2h1dwt+z))/pahdwt);
  if (xtest>=xlimitr && xtest<=xlimitl && ytest<=ylimitd && mydOW>0 ){
  i=32;
  }else{
  i=37;
  }
  }

															/* neutron interacts with the straight down OUTER wall - comments are analog to inner elliptic right wall*/

  if (mydOW!=-1 && t2h1dwt > 0 && t2h1dwt < t1 && t2h1dwt < t2w1r && t2h1dwt < t2w1l && t2h1dwt < t2h1u  && t2h1dwt < t2h1d && t2h1dwt < t2w1rwt && t2h1dwt < t2w1lwt && t2h1dwt < t2h1uwt  && (louthd == 0 && linhd == 0)){
  if(linwr==0 && loutwr==0)
  {
    xlimitr=(-w2r+w1r)/l*(vz*t2h1dwt+z)-w1rwt;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xlimitr=-bwrwt*sqrt(1-((z0wr+(vz*t2h1dwt+z))*(z0wr+(vz*t2h1dwt+z)))/a2wrwt);
    }else{
    xlimitr=-sqrt(((vz*t2h1dwt+z)-pbwrwt)/-pawrwt);
    }
  }
  if(linwl==0 && loutwl == 0)
   {
      xlimitl=(w2l-w1l)/l*(vz*t2h1dwt+z)+w1lwt;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xlimitl=bwlwt*sqrt(1-((z0wl+(vz*t2h1dwt+z))*(z0wl+(vz*t2h1dwt+z)))/a2wlwt);
      }else{
      xlimitl=sqrt(((vz*t2h1dwt+z)-pbwlwt)/-pawlwt);
      }
   }
  xtest=vx*t2h1dwt+x;
  ytest=vy*t2w1lwt+y;
  ylimitd=(-h2d+h1d)*(vz*t2h1dwt+z)/l-h1dwt;
  if (xtest>=xlimitr && xtest<=xlimitl && ytest<=ylimitd && mydOW>0 ){
  i=33;
  }else{
  i=37;
  }
  }

   /* SAME INTERACTIONS BUT TRANSPERENT OUTER WALLS */

  if (mxrOW==-1 && t2w1rwt > 0 && t2w1rwt < t1 && t2w1rwt < t2w1r && t2w1rwt <t2w1l && t2w1rwt < t2h1u && t2w1rwt < t2h1d &&  t2w1rwt <t2w1lwt && t2w1rwt < t2h1uwt && t2w1rwt < t2h1dwt && ((linwr!=0) && (loutwr!=0) )) {
  i=38;
  }

  if (mxrOW==-1 && t2w1rwt > 0 && t2w1rwt < t1 && t2w1rwt < t2w1r && t2w1rwt <t2w1l && t2w1rwt < t2h1u && t2w1rwt < t2h1d &&  t2w1rwt <t2w1lwt && t2w1rwt < t2h1uwt && t2w1rwt < t2h1dwt && ((loutwr!=0 && linwr==0) || (loutwr==0 && linwr!=0))){
  i=38;
  }

  if (mxrOW==-1 && t2w1rwt > 0 && t2w1rwt < t1 && t2w1rwt < t2w1r && t2w1rwt <t2w1l && t2w1rwt < t2h1u && t2w1rwt < t2h1d &&  t2w1rwt <t2w1lwt && t2w1rwt < t2h1uwt && t2w1rwt < t2h1dwt && (loutwr==0 && linwr==0)){
  i=38;
  }

  if (mxlOW==-1 && t2w1lwt > 0 && t2w1lwt < t1 && t2w1lwt <t2w1r && t2w1lwt < t2w1l && t2w1lwt < t2h1u && t2w1lwt < t2h1d &&  t2w1lwt < t2w1rwt && t2w1lwt < t2h1uwt && t2w1lwt < t2h1dwt && ((linwl!=0) && (loutwl!=0) )){
  i=39;
  }

  if (mxlOW==-1 && t2w1lwt > 0 && t2w1lwt < t1 && t2w1lwt <t2w1r && t2w1lwt < t2w1l && t2w1lwt < t2h1u && t2w1lwt < t2h1d &&  t2w1lwt < t2w1rwt && t2w1lwt < t2h1uwt && t2w1lwt < t2h1dwt &&  ((loutwl!=0 && linwl==0) || (loutwl==0 && linwl!=0))){
  i=39;
  }

  if (mxlOW==-1 && t2w1lwt > 0 && t2w1lwt < t1 && t2w1lwt <t2w1r && t2w1lwt < t2w1l && t2w1lwt < t2h1u && t2w1lwt < t2h1d &&  t2w1lwt < t2w1rwt && t2w1lwt < t2h1uwt && t2w1lwt < t2h1dwt && (loutwl==0 && linwl==0)){
  i=39;
  }

  if (myuOW==-1 && t2h1uwt > 0 && t2h1uwt < t1 && t2h1uwt < t2w1r && t2h1uwt < t2w1l && t2h1uwt < t2h1d && t2h1uwt < t2h1u && t2h1uwt < t2w1rwt && t2h1uwt < t2w1lwt && t2h1uwt < t2h1dwt && (louthu !=0 && linhu!=0)){
  i=40;
  }

  if (myuOW==-1 && t2h1uwt > 0 && t2h1uwt < t1 && t2h1uwt < t2w1r && t2h1uwt < t2w1l && t2h1uwt < t2h1d && t2h1uwt < t2h1u && t2h1uwt < t2w1rwt && t2h1uwt < t2w1lwt && t2h1uwt < t2h1dwt && ((louthu!=0 && linhu==0) || (louthu==0 && linhu!=0))){
  i=40;
  }

  if (myuOW==-1 && t2h1uwt > 0 && t2h1uwt < t1 && t2h1uwt < t2w1r && t2h1uwt < t2w1l && t2h1uwt < t2h1d && t2h1uwt < t2h1u && t2h1uwt < t2w1rwt && t2h1uwt < t2w1lwt && t2h1uwt < t2h1dwt && (louthu == 0 && linhu == 0)){
  i=40;
  }

  if (mydOW==-1 && t2h1dwt > 0 && t2h1dwt < t1 && t2h1dwt < t2w1r && t2h1dwt < t2w1l && t2h1dwt < t2h1u  && t2h1dwt < t2h1d && t2h1dwt < t2w1rwt && t2h1dwt < t2w1lwt && t2h1dwt < t2h1uwt  && (louthd !=0 && linhd!=0)){
  i=41;
  }

  if (mydOW==-1 && t2h1dwt > 0 && t2h1dwt < t1 && t2h1dwt < t2w1r && t2h1dwt < t2w1l && t2h1dwt < t2h1u  && t2h1dwt < t2h1d && t2h1dwt < t2w1rwt && t2h1dwt < t2w1lwt && t2h1dwt < t2h1uwt  && ((louthd !=0 && linhd==0) || (louthd ==0 && linhd!=0))){
  i=41;
  }

  if (mydOW==-1 && t2h1dwt > 0 && t2h1dwt < t1 && t2h1dwt < t2w1r && t2h1dwt < t2w1l && t2h1dwt < t2h1u  && t2h1dwt < t2h1d && t2h1dwt < t2w1rwt && t2h1dwt < t2w1lwt && t2h1dwt < t2h1uwt  && (louthd == 0 && linhd == 0)){
  i=41;
  }



switch(i){ /* the principal for the calculation is in every case the same: 1.) one needs the surface normal vector at the intersection point. 2.) calculation of the velocity vector after the interaction by  */
           /* vector subrtation (the basic idea and explanations can be found in the 'Mcstas component manual' in the section 'straight guide') */

  case 1: 										/* no interaction, propagation to the end of the guide */
        PROP_DT(t1);
        break;

  case 2: 
        PROP_DT(t2w1r);									/* propagation to interaction point */
        vxin=vx;									/* saving the velocity vector before the interaction*/
        vyin=vy;
        vzin=vz;
        nx=-x;										/* surface normal vector components at the intersection point */
        nz=-x*x/((a2wr/(z+z0wr))-(z0wr+z));
        n2=sqrt(nx*nx+nz*nz);								/* lenght of the surface normal */
        pf=2.0*(vx*nx+vz*nz)/n2;							/* prefactor for the calculation of the velocity vector after the interaction */
        vx-=pf*nx/n2;									/* velocity vector after the interaction*/
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));	/* calculation the q-vector to calculated the reflectivity*/
        break;

 case 3: 
        PROP_DT(t2w1r);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=-x;
        nz=-0.5/pawr;
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 4: 
        PROP_DT(t2w1r);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=l;
        nz=w2r-w1r;
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 5: 
        PROP_DT(t2w1l);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=-x;
        nz=-x*x/((a2wl/(z+z0wl))-(z0wl+z));
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        SCATTER;
        break;

 case 6: 
        PROP_DT(t2w1l);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=-x;
        nz=-0.5/pawl;
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 7: 
        PROP_DT(t2w1l);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=-l;
        nz=w2l-w1l;
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 8: 
        PROP_DT(t2h1u);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=-y;
        nz=-y*y/((a2hu/(z+z0hu))-(z0hu+z));
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

 case 9: 
        PROP_DT(t2h1u);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=-y;
        nz=-0.5/pahu;
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 10: 
        PROP_DT(t2h1u);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=-l;
        nz=h2u-h1u;
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 11: 
        PROP_DT(t2h1d);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=-y;
        nz=-y*y/((a2hd/(z+z0hd))-(z0hd+z));
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

 case 12: 
        PROP_DT(t2h1d);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=-y;
        nz=-0.5/pahd;
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 13: 
        PROP_DT(t2h1d);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=l;
        nz=h2d-h1d;
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 14:										/* transperent walls - no interaction */
	 PROP_DT(t2w1r);
         break;

  case 15:
	 PROP_DT(t2w1l);
         break;

  case 16:
	 PROP_DT(t2h1u);
         break;

  case 17:
	 PROP_DT(t2h1d);
         break;
  
   case 18:										/* absorbing walls - neutrons are absorbed at interaction point*/
	 PROP_DT(t2w1r);
         ABSORB;
         break;

  case 19:
	 PROP_DT(t2w1l);
         ABSORB;
         break;

  case 20:
	 PROP_DT(t2h1u);
         ABSORB;
         break;

  case 21:
	 PROP_DT(t2h1d);
         ABSORB;
         break;

											/* OUTER WALLS - analog to inner walls, but sign of surface normal vector is changed */

  case 22: 
        PROP_DT(t2w1rwt);								/* propagation to interaction point */
        vxin=vx;									/* saving the velocity vector before the interaction*/
        vyin=vy;
        vzin=vz;
        nx=x;										/* surface normal vector components at the intersection point */
        nz=x*x/((a2wrwt/(z+z0wr))-(z0wr+z));
        n2=sqrt(nx*nx+nz*nz);								/* lenght of the surface normal */
        pf=2.0*(vx*nx+vz*nz)/n2;							/* prefactor for the calculation of the velocity vector after the interaction */
        vx-=pf*nx/n2;									/* velocity vector after the interaction*/
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));	/* calculation the q-vector to calculated the reflectivity*/
        break;

 case 23: 
        PROP_DT(t2w1rwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=x;
        nz=0.5/pawrwt;
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 24: 
        PROP_DT(t2w1rwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=-l;
        nz=-(w2r-w1r);
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 25: 
        PROP_DT(t2w1lwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=x;
        nz=x*x/((a2wlwt/(z+z0wl))-(z0wl+z));
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

 case 26: 
        PROP_DT(t2w1lwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=x;
        nz=0.5/pawlwt;
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 27: 
        PROP_DT(t2w1lwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        nx=l;
        nz=-(w2l-w1l);
        n2=sqrt(nx*nx+nz*nz);
        pf=2.0*(vx*nx+vz*nz)/n2;
        vx-=pf*nx/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 28: 
        PROP_DT(t2h1uwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=y;
        nz=y*y/((a2huwt/(z+z0hu))-(z0hu+z));
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

 case 29: 
        PROP_DT(t2h1uwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=y;
        nz=0.5/pahuwt;
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 30: 
        PROP_DT(t2h1uwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=l;
        nz=-(h2u-h1u);
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 31: 
        PROP_DT(t2h1dwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=y;
        nz=y*y/((a2hdwt/(z+z0hd))-(z0hd+z));
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

 case 32: 
        PROP_DT(t2h1dwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=y;
        nz=0.5/pahdwt;
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 33: 
        PROP_DT(t2h1dwt);
        vxin=vx;
        vyin=vy;
        vzin=vz;
        ny=-l;
        nz=-(h2d-h1d);
        n2=sqrt(ny*ny+nz*nz);
        pf=2.0*(vy*ny+vz*nz)/n2;
        vy-=pf*ny/n2;
        vz-=pf*nz/n2;
        q=V2Q*sqrt((vxin-vx)*(vxin-vx)+(vyin-vy)*(vyin-vy)+(vzin-vz)*(vzin-vz));
        break;

  case 34:
	 PROP_DT(t2w1rwt);
         ABSORB;
         break;

  case 35:
	 PROP_DT(t2w1lwt);
         ABSORB;
         break;

  case 36:
	 PROP_DT(t2h1uwt);
         ABSORB;
         break;

  case 37:
	 PROP_DT(t2h1dwt);
         ABSORB;
         break;

  case 38:
	 PROP_DT(t2w1rwt);
         break;

  case 39:
	 PROP_DT(t2w1lwt);
         break;

  case 40:
	 PROP_DT(t2h1uwt);
         break;

  case 41:
	 PROP_DT(t2h1dwt);
         break;

 }




  if (((i==2) ||(i==3) || (i == 4 ))){								/* calculating the the probability that the neutron is reflected at the RIGHT INNER wall*/
          if (RIreflect && strlen(RIreflect))
              {
               p=Table_Value(riTable, q, 1);
              }else{
                 if(mxr > 0 && q > Qcxr){
                 double arg = (q - mxr*Qcxr)/Wxr;
                 if(arg<10){
                   p *= 0.5*(1.0-tanh(arg))*(1.0-alphaxr*(q-Qcxr));
                 }else
                 ABSORB;
                }
              }
           }

  if (((i==22) ||(i==23) || (i==24 ))){								/* calculating the the probability that the neutron is reflected at the RIGHT OUTER wall*/
          if (ROreflect && strlen(ROreflect))
              {
               p=Table_Value(roTable, q, 1);
              }else{
                 if(mxrOW > 0 && q > QcxrOW){
                   double arg = (q - mxrOW*QcxrOW)/WxrOW;
                   if(arg<10){
                      p *= 0.5*(1.0-tanh(arg))*(1.0-alphaxrOW*(q-QcxrOW));
                   }
                   else
                   ABSORB;
                  }
               }
             }

  if (((i==5) ||(i==6) || (i == 7 ) ) ){							/* calculating the the probability that the neutron is reflected at the LEFT INNER wall*/
          if (LIreflect && strlen(LIreflect))
              {
               p=Table_Value(liTable, q, 1);
              }else{
                 if(mxl > 0 && q > Qcxl){
                   double arg = (q - mxl*Qcxl)/Wxl;
                   if(arg<10){
                     p *= 0.5*(1.0-tanh(arg))*(1.0-alphaxl*(q-Qcxl));
                   }else
                   ABSORB;
                }
               }
          }

  if (((i==25) || (i==26) || (i==27 )) && mxlOW >0 && q > QcxlOW ){				/* calculating the the probability that the neutron is reflected at the LEFT OUTER wall*/
         if (LOreflect && strlen(LOreflect))
              {
               p=Table_Value(loTable, q, 1);
              }else{
                 if(mxlOW > 0 && q > QcxlOW){
                   double arg = (q - mxlOW*QcxlOW)/WxlOW;
                   if(arg<10){
                      p *= 0.5*(1.0-tanh(arg))*(1.0-alphaxlOW*(q-QcxlOW));
                   }
                   else
                   ABSORB;
                 }
               }
          }

  if (((i==8) ||(i==9) || (i == 10 ))){								/* calculating the the probability that the neutron is reflected at the TOP INNER wall*/
	  if (UIreflect && strlen(UIreflect))
              {
               p=Table_Value(uiTable, q, 1);
              }else{
                 if(myu > 0 && q > Qcyu){
                 double arg = (q - myu*Qcyu)/Wyu;
                   if(arg<10){
                     p *= 0.5*(1.0-tanh(arg))*(1.0-alphayu*(q-Qcyu));
                   }else
                   ABSORB;
                 }
              }
            }

  if (((i==28) || (i==29) || (i==30 )) ){							/* calculating the the probability that the neutron is reflected at the TOP OUTER wall*/
	  if (UOreflect && strlen(UOreflect))
              {
               p=Table_Value(uoTable, q, 1);
              }else{
                 if(myuOW > 0 && q > QcyuOW){
                   double arg = (q - myuOW*QcyuOW)/WyuOW;
                   if(arg<10){
                      p *= 0.5*(1.0-tanh(arg))*(1.0-alphayuOW*(q-QcyuOW));
                   }else
                   ABSORB;
                   }
               }
           }

  if (((i==11) ||(i==12) || (i == 13 ))){							/* calculating the the probability that the neutron is reflected at the BOTTOM INNER wall*/
          if (DIreflect && strlen(DIreflect))
              {
               p=Table_Value(diTable, q, 1);
              }else{
                 if(myd > 0 && q > Qcyd){
                   double arg = (q - myd*Qcyd)/Wyd;
                   if(arg<10){
                      p *= 0.5*(1.0-tanh(arg))*(1.0-alphayd*(q-Qcyd));
                   }else
                   ABSORB;
                   }
               }
            }

   if (((i==31) || (i==32) || (i==33 )) ){							/* calculating the the probability that the neutron is reflected at the BOTTOM OUTER wall*/
          if (DOreflect && strlen(DOreflect))
              {
               p=Table_Value(doTable, q, 1);
              }else{
                 if(mydOW > 0 && q > QcydOW){
                   double arg = (q - mydOW*QcydOW)/WydOW;
                   if(arg<10){
                      p *= 0.5*(1.0-tanh(arg))*(1.0-alphaydOW*(q-QcydOW));
                   }else
                   ABSORB;
                   }
                 }
           }


  p *= R0;
  SCATTER;

}while(z<l);												/* repeat the interaction loop untill the neutron pass the end of guide */


  if(x <= -w2r && x >= -w2rwt && y <= mru2*x+nru2 && y >= mrd2*x+nrd2 && mxr!=-1 && mxrOW!=-1)		/* absorbing the neutron if it hit the RIGHT exit wall and the wall is not transparent*/
      ABSORB;
  if(x >=  w2l && x <=  w2lwt && y <= mlu2*x+nlu2 && y >= mld2*x+nld2 && mxl!=-1 && mxlOW!=-1) 		/* absorbing the neutron if it hit the LEFT exit wall and the wall is not transparent*/
      ABSORB;
  if(y <= -h2d && y >= -h2dwt && x <= (y-nld2)/mld2 && x>= (y-nrd2)/mrd2 && myd!=-1 && mydOW!=-1) 	/* absorbing the neutron if it hit the BOTTOM exit wall and the wall is not transparent*/
      ABSORB; 
  if(y >=  h2u && y <=  h2uwt && x <= (y-nlu2)/mlu2 && x>= (y-nru2)/mru2 && myu!=-1 && myuOW!=-1)	/* absorbing the neutron if it hit the TOP exit wall and the wall is not transparent*/
      ABSORB;


%}



FINALLY
%{
printf("w2r = %f \n", w2r);
printf("w2l = %f \n", w2l);
printf("h2u = %f \n", h2u);
printf("h2d = %f \n", h2d);
%}

MCDISPLAY
%{
  int i,imax;
  double x,y,z,x1,y1,z1,z0wr,z0wl,z0hu,z0hd,xwt,ywt,x1wt,y1wt;
  double mr,ml,mu,md,nr1,nl1,nu1,nd1,nr2,nl2,nu2,nd2;
  double lbwl,lbwr,lbhu,lbhd;								/* length between focal points , needed for elliptic case */

  magnify("xy");


  imax=100;										/* maximum points for every line in z direction*/

  lbwr=linwr+l+loutwr;
  lbwl=linwl+l+loutwl;
  lbhu=linhu+l+louthu;
  lbhd=linhd+l+louthd;

  if (linwr==0 && loutwr==0){
   mr=(-w2r+w1r)/l;
   nr1=-w1r;
   nr2=-(w1rwt);
   }


if (linwl==0 && loutwl==0){
   ml=(w2l-w1l)/l;
   nl1=w1l;
   nl2=(w1lwt);
   }


if (linhu == 0 && louthu==0) 
  {
   mu=(h2u-h1u)/l;
   nu1=h1u;
   nu2=(h1uwt);
  }


if (linhd == 0 && louthd==0) 
  {
   md=(-h2d+h1d)/l;
   nd1=-h1d;
   nd2=-(h1dwt);
  }

  z0wr=(linwr-l-loutwr)/2.0;
  z0wl=(linwl-l-loutwl)/2.0;
  z0hu=lbhu/2.0-l-louthu;
  z0hd=lbhd/2.0-l-louthd;

  if(myd!=-1) line(w1l, -h1d, 0.0, -w1r, -h1d, 0.0);							/* entrance window given by the INNER walls*/
  if(myu!=-1)line(w1l, h1u, 0.0, -w1r, h1u, 0.0);
  if(mxl!=-1)line(w1l, -h1d, 0.0, w1l, h1u, 0.0);
  if(mxr!=-1)line( -w1r, h1u, 0.0, -w1r, -h1d, 0.0);

  if(myd!=-1)line(w2l, -h2d, l, -w2r, -h2d, l);								/* exit window given by the INNER walls*/
  if(myu!=-1)line(w2l, h2u,  l, -w2r, h2u, l);
  if(mxl!=-1)line(w2l, -h2d, l, w2l, h2u, l);
  if(mxr!=-1)line( -w2r, -h2d,  l, -w2r, h2u, l);

  if(mydOW!=-1) line((w1lwt), -(h1dwt), 0.0, -(w1rwt), -(h1dwt), 0.0);					/* entrance window given by the OUTER walls */
  if(myuOW!=-1)line((w1lwt), (h1uwt), 0.0, -(w1rwt), (h1uwt), 0.0);
  if(mxlOW!=-1)line((w1lwt), -(h1dwt), 0.0, (w1lwt), (h1uwt), 0.0);
  if(mxrOW!=-1)line( -(w1rwt), (h1uwt), 0.0, -(w1rwt), -(h1dwt), 0.0);

  if(mydOW!=-1)line((w2lwt), -(h2dwt), l, -(w2rwt), -(h2dwt), l);					/* exit windows given by the OUTER walls*/
  if(myuOW!=-1)line((w2lwt), (h2uwt),  l, -(w2rwt), (h2uwt), l);
  if(mxlOW!=-1)line((w2lwt), -(h2dwt), l, (w2lwt), (h2uwt), l);
  if(mxrOW!=-1)line( -(w2rwt), -(h2dwt),  l, -(w2rwt), (h2uwt), l);

  if((myd!=-1 && mydOW!=-1) || (mxl!=-1 && mxlOW!=-1)) line(w1l, -h1d, 0.0, (w1lwt), -(h1dwt), 0.0);	/* corner connection lines for the entrance windows*/	
  if((myu!=-1 && myuOW!=-1) || (mxl!=-1 && mxlOW!=-1)) line(w1l, h1u, 0.0, (w1lwt), (h1uwt), 0.0);
  if((myd!=-1 && mydOW!=-1) || (mxr!=-1 && mxrOW!=-1)) line(-w1r, -h1d, 0.0,-(w1rwt),  -(h1dwt), 0.0);
  if((myu!=-1 && myuOW!=-1) || (mxr!=-1 && mxrOW!=-1)) line( -w1r, h1u, 0.0, -(w1rwt), (h1uwt), 0.0);

  if((myd!=-1 && mydOW!=-1) || (mxl!=-1 && mxlOW!=-1)) line(w2l, -h2d, l, (w2lwt), -(h2dwt), l);	/* corner connection lines for the exit windows*/
  if((myu!=-1 && myuOW!=-1) || (mxl!=-1 && mxlOW!=-1)) line(w2l, h2u, l, (w2lwt), (h2uwt), l);
  if((myd!=-1 && mydOW!=-1) || (mxr!=-1 && mxrOW!=-1)) line(-w2r, -h2d, l,-(w2rwt),  -(h2dwt), l);
  if((myu!=-1 && myuOW!=-1) || (mxr!=-1 && mxrOW!=-1)) line( -w2r, h2u, l, -(w2rwt), (h2uwt), l);

for(i=0;i<imax;i++){											/* calculation of the points for the INNER LEFT BOTTOM line */
   z=i*l/imax;
   z1=(i+1)*l/imax;
  if(linwl==0 && loutwl == 0)
   {
     x=ml*z+nl1;
     x1=ml*z1+nl1;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      x=bwl*sqrt(1-((z0wl+z)*(z0wl+z))/(awl*awl));
      x1=bwl*sqrt(1-((z0wl+z1)*(z0wl+z1))/(awl*awl));
      }else{
      x=sqrt((z-pbwl)/-pawl);
      x1=sqrt((z1-pbwl)/-pawl);
      }
   }
  if(linhd==0 && louthd==0)
   {
     y=md*z+nd1;
     y1=md*z1+nd1;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      y=-bhd*sqrt(1-((z0hd+z)*(z0hd+z))/(ahd*ahd));
      y1=-bhd*sqrt(1-((z0hd+z1)*(z0hd+z1))/(ahd*ahd));
      }else{
      y=-sqrt((z-pbhd)/-pahd);
      y1=-sqrt((z1-pbhd)/-pahd);
      }
   }
  if(mxl!=-1 || myd!=-1)line((double)x, (double)y, (double)z, (double)x1, (double)y1, (double)z1); 
 }

for(i=0;i<imax;i++){											/* calculation of the points for the INNER RIGHT BOTTOM line */
  z=i*l/imax;
  z1=(i+1)*l/imax;
  if(linwr==0 && loutwr==0)
  {
    x=mr*z+nr1;
    x1=mr*z1+nr1;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    x=-bwr*sqrt(1-((z0wr+z)*(z0wr+z))/(awr*awr));
    x1=-bwr*sqrt(1-((z0wr+z1)*(z0wr+z1))/(awr*awr));
    }else{
    x=-sqrt((z-pbwr)/-pawr);
    x1=-sqrt((z1-pbwr)/-pawr);
    }
  }
  if(linhd==0 && louthd==0)
  {
    y=md*z+nd1;
    y1=md*z1+nd1;
  }else{
    if(linhd!=0 && louthd!=0)
    {
    y=-bhd*sqrt(1-((z0hd+z)*(z0hd+z))/(ahd*ahd));
    y1=-bhd*sqrt(1-((z0hd+z1)*(z0hd+z1))/(ahd*ahd));
    }else{
    y=-sqrt((z-pbhd)/-pahd);
    y1=-sqrt((z1-pbhd)/-pahd);
    }
  }
  if(mxr!=-1 || myd!=-1) line((double)x, (double)y, (double)z, (double)x1, (double)y1, (double)z1); 
}


for(i=0;i<imax;i++){											/* calculation of the points for the INNER RIGHT TOP line */
  z=i*l/imax;
  z1=(i+1)*l/imax;
  if(linwr==0 && loutwr==0)
  {
    x=mr*z+nr1;
    x1=mr*z1+nr1;
   }else{
     if(linwr!=0 && loutwr!=0){
     x=-bwr*sqrt(1-((z0wr+z)*(z0wr+z))/(awr*awr));
     x1=-bwr*sqrt(1-((z0wr+z1)*(z0wr+z1))/(awr*awr));
     }else{
     x=-sqrt((z-pbwr)/-pawr);
     x1=-sqrt((z1-pbwr)/-pawr);
     }
   }
   if(linhu==0 && louthu==0)
  {
    y=mu*z+nu1;
    y1=mu*z1+nu1;
  }else{
    if(linhu!=0 && louthu!=0)
    {
    y=bhu*sqrt(1-((z0hu+z)*(z0hu+z))/(ahu*ahu));
    y1=bhu*sqrt(1-((z0hu+z1)*(z0hu+z1))/(ahu*ahu));
    }else{
    y=sqrt((z-pbhu)/-pahu);
    y1=sqrt((z1-pbhu)/-pahu);
    }
  }
  if(mxr!=-1 || myu!=-1)line((double)x, (double)y, (double)z, (double)x1, (double)y1, (double)z1); 
}

for(i=0;i<imax;i++){											/* calculation of the points for the INNER LEFT TOP line */
  z=i*l/imax;
  z1=(i+1)*l/imax;
  if(linwl==0 && loutwl==0)
  {
    x=ml*z+nl1;
    x1=ml*z1+nl1;
   }else{
     if(linwl!=0 && loutwl!=0)
     {
     x=bwl*sqrt(1-((z0wl+z)*(z0wl+z))/(awl*awl));
     x1=bwl*sqrt(1-((z0wl+z1)*(z0wl+z1))/(awl*awl));
     }else{
     x=sqrt((z-pbwl)/-pawl);
     x1=sqrt((z1-pbwl)/-pawl);
     }
   }
   if(linhu==0 && louthu==0)
  {
    y=mu*z+nu1;
    y1=mu*z1+nu1;
  }else{
    if(linhu!=0 && louthu!=0)
    {
    y=bhu*sqrt(1-((z0hu+z)*(z0hu+z))/(ahu*ahu));
    y1=bhu*sqrt(1-((z0hu+z1)*(z0hu+z1))/(ahu*ahu));
    }else{
    y=sqrt((z-pbhu)/-pahu);
    y1=sqrt((z1-pbhu)/-pahu);
    }
  }
  if(mxl!=-1 || myu!=-1)line((double)x, (double)y, (double)z, (double)x1, (double)y1, (double)z1); 
}																	/* END INNER LINES*/


for(i=0;i<imax;i++){											/* calculation of the points for the OUTER LEFT BOTTOM line */
   z=i*l/imax;
   z1=(i+1)*l/imax;
  if(linwl==0 && loutwl == 0)
   {
     xwt=ml*z+nl2;
     x1wt=ml*z1+nl2;
   }else{
      if(linwl!=0 && loutwl != 0)
      {
      xwt=bwlwt*sqrt(1-((z0wl+z)*(z0wl+z))/(awlwt*awlwt));
      x1wt=bwlwt*sqrt(1-((z0wl+z1)*(z0wl+z1))/(awlwt*awlwt));
      }else{
      xwt=sqrt((z-pbwlwt)/-pawlwt);
      x1wt=sqrt((z1-pbwlwt)/-pawlwt);
      }
   }
  if(linhd==0 && louthd==0)
   {
     ywt=md*z+nd2;
     y1wt=md*z1+nd2;
    }else{
      if(linhd!=0 && louthd!=0)
      {
      ywt=-bhdwt*sqrt(1-((z0hd+z)*(z0hd+z))/(ahdwt*ahdwt));
      y1wt=-bhdwt*sqrt(1-((z0hd+z1)*(z0hd+z1))/(ahdwt*ahdwt));
      }else{
      ywt=-sqrt((z-pbhdwt)/-pahdwt);
      y1wt=-sqrt((z1-pbhdwt)/-pahdwt);
      }
   }
  if(mxlOW!=-1 || mydOW!=-1)line((double)xwt, (double)ywt, (double)z, (double)x1wt, (double)y1wt, (double)z1); 
 }

for(i=0;i<imax;i++){											/* calculation of the points for the OUTER RIGHT BOTTOM line */
  z=i*l/imax;
  z1=(i+1)*l/imax;
  if(linwr==0 && loutwr==0)
  {
    xwt=mr*z+nr2;
    x1wt=mr*z1+nr2;
  }else{
    if(linwr!=0 && loutwr!=0)
    {
    xwt=-bwrwt*sqrt(1-((z0wr+z)*(z0wr+z))/(awrwt*awrwt));
    x1wt=-bwrwt*sqrt(1-((z0wr+z1)*(z0wr+z1))/(awrwt*awrwt));
    }else{
    xwt=-sqrt((z-pbwrwt)/-pawrwt);
    x1wt=-sqrt((z1-pbwrwt)/-pawrwt);
    }
  }
  if(linhd==0 && louthd==0)
  {
    ywt=md*z+nd2;
    y1wt=md*z1+nd2;
  }else{
    if(linhd!=0 && louthd!=0)
    {
    ywt=-bhdwt*sqrt(1-((z0hd+z)*(z0hd+z))/(ahdwt*ahdwt));
    y1wt=-bhdwt*sqrt(1-((z0hd+z1)*(z0hd+z1))/(ahdwt*ahdwt));
    }else{
    ywt=-sqrt((z-pbhdwt)/-pahdwt);
    y1wt=-sqrt((z1-pbhdwt)/-pahdwt);
    }
  }
  if(mxrOW!=-1 || mydOW!=-1) line((double)xwt, (double)ywt, (double)z, (double)x1wt, (double)y1wt, (double)z1); 
}


for(i=0;i<imax;i++){											/* calculation of the points for the OUTER RIGHT TOP line */
  z=i*l/imax;
  z1=(i+1)*l/imax;
  if(linwr==0 && loutwr==0)
  {
    xwt=mr*z+nr2;
    x1wt=mr*z1+nr2;
   }else{
     if(linwr!=0 && loutwr!=0){
     xwt=-bwrwt*sqrt(1-((z0wr+z)*(z0wr+z))/(awrwt*awrwt));
     x1wt=-bwrwt*sqrt(1-((z0wr+z1)*(z0wr+z1))/(awrwt*awrwt));
     }else{
     xwt=-sqrt((z-pbwrwt)/-pawrwt);
     x1wt=-sqrt((z1-pbwrwt)/-pawrwt);
     }
   }
   if(linhu==0 && louthu==0)
  {
    ywt=mu*z+nu2;
    y1wt=mu*z1+nu2;
  }else{
    if(linhu!=0 && louthu!=0)
    {
    ywt=bhuwt*sqrt(1-((z0hu+z)*(z0hu+z))/(ahuwt*ahuwt));
    y1wt=bhuwt*sqrt(1-((z0hu+z1)*(z0hu+z1))/(ahuwt*ahuwt));
    }else{
    ywt=sqrt((z-pbhuwt)/-pahuwt);
    y1wt=sqrt((z1-pbhuwt)/-pahuwt);
    }
  }
   if(mxrOW!=-1 || myuOW!=-1)line((double)xwt, (double)ywt, (double)z, (double)x1wt, (double)y1wt, (double)z1); 
}

for(i=0;i<imax;i++){											/* calculation of the points for the OUTER LEFT TOP line */
  z=i*l/imax;
  z1=(i+1)*l/imax;
  if(linwl==0 && loutwl==0)
  {
    xwt=ml*z+nl2;
    x1wt=ml*z1+nl2;
   }else{
     if(linwl!=0 && loutwl!=0)
     {
     xwt=bwlwt*sqrt(1-((z0wl+z)*(z0wl+z))/(awlwt*awlwt));
     x1wt=bwlwt*sqrt(1-((z0wl+z1)*(z0wl+z1))/(awlwt*awlwt));
     }else{
     xwt=sqrt((z-pbwlwt)/-pawlwt);
     x1wt=sqrt((z1-pbwlwt)/-pawlwt);
     }
   }
   if(linhu==0 && louthu==0)
  {
    ywt=mu*z+nu2;
    y1wt=mu*z1+nu2;
  }else{
    if(linhu!=0 && louthu!=0)
    {
    ywt=bhuwt*sqrt(1-((z0hu+z)*(z0hu+z))/(ahuwt*ahuwt));
    y1wt=bhuwt*sqrt(1-((z0hu+z1)*(z0hu+z1))/(ahuwt*ahuwt));
    }else{
    ywt=sqrt((z-pbhuwt)/-pahuwt);
    y1wt=sqrt((z1-pbhuwt)/-pahuwt);
    }
  }
  if(mxlOW!=-1 || myuOW!=-1)line((double)xwt, (double)ywt, (double)z, (double)x1wt, (double)y1wt, (double)z1); 
}



%}

END



