/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Guide_tapering.comp
*
* Component: Guide_tapering
*
* %I
* Written by: Uwe Filges
* Date: 22/10/2003
* Origin: PSI
* Modified by: Rob Dalgliesh, ISIS, 2007
*
* Models a rectangular tapered guide (many shapes)
*
* %D
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane.
* The guide may be tapered.
*
* The component includes a feature to read in self-defined functions for
* guide tapering. Under the parameter 'option' the KEYWORD 'file=' offers
* the possibility to read in parameters from an ASC-file. The file structure
* is shown below in the example. It is important to know that the first
* 3 lines will be interpreted as comments.
* Afterwards the dimension of each guide segment must be defined. The
* length of each segment is constant l(i)=l/segments . The number of
* segments is defined through the number of lines minus the first 3
* lines taken from the Input-File.
* The guide can be made curved both horizontally and vertically (not shown in 3d
* view), and m-coating, when set negative, is varied in 1/width.
*
* Example Input-File:
*
* c Guide_tapering.comp
* c i = 0 - 199 segments
* c h1(i)     h2(i)   w1(i)    w2(i)
* 0.120000 0.119850 0.020000 0.020000
* 0.119850 0.119700 0.020000 0.020000
* 0.119700 0.119550 0.020000 0.020000
* 0.119550 0.119400 0.020000 0.020000
* 0.119400 0.119250 0.020000 0.020000
* 0.119250 0.119100 0.020000 0.020000
* ...
*
* Example1: Guide_tapering(w1=0.1, h1=0.18, linw=0.1, loutw=0.1, linh=0.1, louth=0.1, l=1.5, option="elliptical", R0=0.99, Qcx=0.021, Qcy=0.021, alphax=6.07, alphay=6.07, W=0.003, mx=1, my=1, segno=800)
*
* Example2: Guide_tapering(w1=0, h1=0, linw=0, loutw=0, linh=0, louth=0, l=1.5, option="file=ownfunction.txt", R0=0.99, Qcx=0.021, Qcy=0.021, alphax=6.07, alphay=6.07, W=0.003, mx=1, my=1)
*
* %BUGS
* This component does not work with gravitation on. Use component Guide_gravity then.
*
* %P
* INPUT PARAMETERS:
*
* w1: [m]      Width at the guide entry
* h1: [m]      Height at the guide entry
* linw: [m]    distance from 1st focal point to real guide entry - left and right horizontal mirrors
* loutw: [m]   distance from real guide exit to 2nd focal point - left and right horizontal mirrors
* l: [m]       length of guide
* linh: [m]         distance from 1st focal point to real guide entry - top and bottom vertical mirrors
* louth: [m]        distance from real guide exit to 2nd focal point - top and bottom vertical mirrors
* option: [str]   define the input function for the curve of the guide walls options are: "elliptical" - define elliptical function of guide walls "parabolical" - define parabolical function of guide walls "straight"    - define a straight elements guide"file=[filename]" - read in ASC-file with arbitrary definition for the curve of the guide walls
* R0: [1]           Low-angle reflectivity
* Qcx: [AA-1]       Critical scattering vector for left and right vertical mirrors in each channel
* Qcy: [AA-1]       Critical scattering vector for top and bottom mirrors
* alphax: [AA]      Slope of reflectivity for left and right vertical mirrors in each channel
* alphay: [AA]      Slope of reflectivity for top and bottom mirrors
* mx: [1]           m-value of material for left and right vertical mirrors in each channel. Zero means completely absorbing. Negative  value will adapt coating as e.g. m=mx*w1/w
* my: [1]           m-value of material for top and bottom mirrors. Zero means completely absorbing. Negative value will adapt coating as e.g. m=my*h1/h
* W: [AA-1]         Width of supermirror cut-off for all mirrors
* segno: [1]        number of segments (z-axis) for cutting the tube
* curvature: [m]    guide horizontal radius of curvature. Zero means straight.
* curvature_v: [m]  guide vertical radius of curvature. Zero means straight.
*
* %End
*
******************************************************************************/

DEFINE COMPONENT Guide_tapering

SETTING PARAMETERS (string option=0, w1=0,h1=0,l,linw=0,loutw=0,linh=0,louth=0, R0=0.99,
Qcx=0.021,Qcy=0.021, alphax=6.07, alphay=6.07, W=0.003,
mx=1, my=1,segno=800,curvature=0,curvature_v=0)

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

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

DECLARE
%{
  double* w1c;
  double* w2c;
  double* ww;
  double* hh;
  double* whalf;
  double* hhalf;
  double* lwhalf;
  double* lhhalf;
  double* h1_in;
  double* h2_out;
  double* w1_in;
  double* w2_out;
  double l_seg;
  double h12;
  double h2;
  double w12;
  double w2;
  double a_ell_q;
  double b_ell_q;
  double lbw;
  double lbh;
  double mxi;
  double u1;
  double u2;
  double div1;
  double p2_para;
  double test;
  double Div1;
  int seg;
  char* fu;
  char* pos;
  char file_name[1024];
  char* ep;
  FILE* num;
  double rotation_h;
  double rotation_v;
%}

INITIALIZE
%{
  int i, ii;

  rotation_h = 0;
  rotation_v = 0;

  // dynamic memory allocation is good
  w1c = (double*)malloc (sizeof (double) * segno);
  w2c = (double*)malloc (sizeof (double) * segno);
  ww = (double*)malloc (sizeof (double) * segno);
  hh = (double*)malloc (sizeof (double) * segno);
  whalf = (double*)malloc (sizeof (double) * segno);
  hhalf = (double*)malloc (sizeof (double) * segno);
  lwhalf = (double*)malloc (sizeof (double) * segno);
  lhhalf = (double*)malloc (sizeof (double) * segno);
  h1_in = (double*)malloc (sizeof (double) * (segno + 1));
  h2_out = (double*)malloc (sizeof (double) * (segno + 1));
  w1_in = (double*)malloc (sizeof (double) * (segno + 1));
  w2_out = (double*)malloc (sizeof (double) * (segno + 1));

  struct para {
    char st[128];
  } segment[800];

  if (W <= 0) {
    fprintf (stderr, "Component: %s (Guide_tapering) W must \n", NAME_CURRENT_COMP);
    fprintf (stderr, "           be positive\n");
    exit (-1);
  }
  if (l <= 0) {
    fprintf (stderr, "Component: %s (Guide_tapering) real guide length \n", NAME_CURRENT_COMP);
    fprintf (stderr, "           is <= ZERO ! \n");
    exit (-1);
  }
  if (mcgravitation)
    fprintf (stderr,
             "WARNING: Guide_tapering: %s: "
             "This component produces wrong results with gravitation !\n"
             "Use Guide_gravity.\n",
             NAME_CURRENT_COMP);
  seg = segno;
  l_seg = l / (seg);
  h12 = h1 / 2.0;
  if (option != NULL) {
    fu = (char*)malloc (sizeof (char) * (strlen (option) + 1));
    strcpy (fu, option);
  } else {
    exit (-1);
  }
  /* handle guide geometry ================================================== */
  if (!strcmp (fu, "elliptical")) {
    /* calculate parameter b of elliptical equestion - vertical mirrors */
    /* (l+linh+louth) -> distance between focal points */
    /*  printf("A1 \n"); */
    lbh = l + linh + louth;
    if (linh == 0 && louth == 0) {
      /* plane mirrors (vertical) */
      b_ell_q = 0;
      h2 = h1;
    } else {
      /* elliptical mirrors */
      u1 = sqrt ((linh * linh) + (h12 * h12));
      u2 = sqrt ((h12 * h12) + ((l + louth) * (l + louth)));
      a_ell_q = ((u1 + u2) / 2.0) * ((u1 + u2) / 2.0);
      b_ell_q = a_ell_q - ((lbh / 2.0) * (lbh / 2.0));
      /* calculate heigth of guide exit (h2) */
      div1 = ((lbh / 2.0 - louth) * (lbh / 2.0 - louth)) / a_ell_q;
      h2 = sqrt (b_ell_q * (1.0 - div1));
      h2 = h2 * 2.0;
    }
  } else if (!strcmp (fu, "parabolical")) {
    if ((linh > 0) && (louth > 0)) {
      fprintf (stderr, "Component: %s (Guide_tapering) Two focal\n", NAME_CURRENT_COMP);
      fprintf (stderr, "            points lout and linh are not allowed! \n");
      free (fu);
      exit (-1);
    }
    if (louth == 0 && linh == 0) {
      /* plane mirrors (vertical) */
      h2 = h1;
    } else {
      /* parabolical mirrors */
      if (linh == 0) {
        Div1 = ((2.0 * louth + 2.0 * l) * (2.0 * louth + 2.0 * l)) / 4.0;
        p2_para = ((sqrt (Div1 + (h12 * h12))) - (louth + l)) * 2.0;
        /* calculate heigth of guide exit (h2) */
        h2 = sqrt (p2_para * (louth + p2_para / 4.0));
        h2 = h2 * 2.0;
      } else {
        /* anti-trompete */
        Div1 = ((2.0 * linh) * (2.0 * linh)) / 4.0;
        p2_para = ((sqrt (Div1 + (h12 * h12))) - linh) * 2.0;
        /* calculate heigth of guide exit (h2) */
        h2 = sqrt (p2_para * (l + linh + p2_para / 4.0));
        h2 = h2 * 2.0;
      }
    }
  } else if (!strncmp (fu, "file", 4)) {
    pos = strtok (fu, "=");
    while ((pos = strtok (0, "="))) {
      strcpy (file_name, pos);
    }
    if ((num = fopen (file_name, "r")) == NULL) {
      fprintf (stderr, "Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
      fprintf (stderr, "           File %s not found! \n", file_name);
      free (fu);
      exit (-1);
    } else {
      ii = 0;
      /* Use ret==NULL as termination criterion on while, otherwise
         the reading may on some systems continue until "800" segments */
      char dymmy = 1;
      char* ret = &dymmy;
      while (!feof (num) && ret) {
        ret = fgets (segment[ii].st, 128, num);
        if (ii > 799) {
          fprintf (stderr, "%s: Number of segments is limited to 800 !! \n", NAME_CURRENT_COMP);
          free (fu);
          exit (-1);
        }
        ii++;
      }
      fclose (num);
      ii--;
    }
    seg = ii - 3;
    l_seg = l / seg;
    for (i = 3; i < ii; i++) {
      if (strlen (segment[i].st) < 4) {
        fprintf (stderr, "Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
        fprintf (stderr, "           Data Format Error! \n");
        free (fu);
        exit (-1);
      }
      h1_in[i - 3] = strtod (strtok (segment[i].st, " "), &ep);
      h2_out[i - 3] = strtod (strtok (0, " "), &ep);
      w1_in[i - 3] = strtod (strtok (0, " "), &ep);
      w2_out[i - 3] = strtod (strtok (0, " "), &ep);
    }
    h1 = h1_in[0];
    h2 = h2_out[seg - 1];
    w1 = w1_in[0];
    w2 = w2_out[seg - 1];
    for (i = 0; i < seg; i++) {
      fprintf (stderr, "%d: %lf %lf %lf %lf \n", i, h1_in[i], h2_out[i], w1_in[i], w2_out[i]);
    }
  } else if (!strcmp (fu, "straight")) {
    for (i = 0; i < seg; i++) {
      h1_in[i] = h2_out[i] = h2 = h1;
      w1_in[i] = w2_out[i] = w2 = w1;
    }
  } else {
    fprintf (stderr, "Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
    fprintf (stderr, "           Unknown KEYWORD: %s \n", fu);
    free (fu);
    exit (-1);
  }
  fprintf (stderr, "Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
  fprintf (stderr, "           Height at the guide exit (h2): %lf \n", h2);
  if (h2 <= 0) {
    fprintf (stderr, "Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
    fprintf (stderr, "           Height at the guide exit (h2) was calculated\n");
    fprintf (stderr, "           <=0; Please change the parameter h1 and/or\n");
    fprintf (stderr, "           linh and/or louth! \n");
    free (fu);
    exit (-1);
  }
  if (!strcmp (fu, "elliptical")) {
    h1_in[0] = h1;
    for (i = 1; i < seg; i++) {
      if (b_ell_q == 0) {
        h1_in[i] = h1;
      } else {
        mxi = (((lbh / 2.0) - linh) - (l_seg * i));
        h1_in[i] = (sqrt ((1.0 - ((mxi * mxi) / a_ell_q)) * b_ell_q)) * 2.0;
      }
      h2_out[i - 1] = h1_in[i];
    }
    h2_out[seg - 1] = h2;
  } else if (!strcmp (fu, "parabolical")) {
    h1_in[0] = h1;
    ii = seg - 1;
    if (louth == 0 && linh == 0) {
      for (i = 1; i < (seg + 1); i++) {
        h1_in[i] = h1;
        ii = ii - 1;
        h2_out[i - 1] = h1_in[i];
      }
    } else {
      if ((linh == 0) && (louth > 0)) {
        for (i = 1; i < (seg + 1); i++) {
          h1_in[i] = (sqrt ((p2_para / 4.0 + louth + (l_seg * ii)) * p2_para)) * 2.0;
          ii = ii - 1;
          h2_out[i - 1] = h1_in[i];
        }
      } else {
        for (i = 1; i < (seg + 1); i++) {
          h1_in[i] = (sqrt ((p2_para / 4.0 + linh + (l_seg * i)) * p2_para)) * 2.0;
          h2_out[i - 1] = h1_in[i];
        }
      }
    }
  }
  /* compute each value for horizontal mirrors */
  w12 = w1 / 2.0;
  if (!strcmp (fu, "elliptical")) {
    /* calculate lbw the distance between focal points of horizontal mirrors */
    lbw = l + linw + loutw;
    /* calculate parameter b of elliptical equestion - horizontal mirrors */
    if (linw == 0 && loutw == 0) {
      /* plane mirrors (horizontal) */
      b_ell_q = 0;
      w2 = w1;
    } else {
      /* elliptical mirrors */
      u1 = sqrt ((linw * linw) + (w12 * w12));
      u2 = sqrt ((w12 * w12) + ((l + loutw) * (l + loutw)));
      a_ell_q = ((u1 + u2) / 2.0) * ((u1 + u2) / 2.0);
      b_ell_q = a_ell_q - ((lbw / 2.0) * (lbw / 2.0));
      /* calculate weigth of guide exit (w2) */
      div1 = ((lbw / 2.0 - loutw) * (lbw / 2.0 - loutw)) / a_ell_q;
      w2 = sqrt (b_ell_q * (1.0 - div1));
      w2 = w2 * 2.0;
    }
  } else if (!strcmp (fu, "parabolical")) {
    if ((linw > 0) && (loutw > 0)) {
      fprintf (stderr, "Component: %s (Guide_tapering) Two focal\n", NAME_CURRENT_COMP);
      fprintf (stderr, "           points linw and loutw are not allowed! \n");
      free (fu);
      exit (-1);
    }
    if (loutw == 0 && linw == 0) {
      /* plane mirrors (horizontal) */
      w2 = w1;
    } else {
      if (linw == 0) {
        /* parabolical mirrors */
        Div1 = ((2.0 * loutw + 2.0 * l) * (2.0 * loutw + 2.0 * l)) / 4.0;
        p2_para = ((sqrt (Div1 + (w12 * w12))) - (loutw + l)) * 2.0;
        /* calculate weigth of guide exit (w2) */
        w2 = sqrt (p2_para * (loutw + p2_para / 4.0));
        w2 = w2 * 2.0;
      } else {
        /* anti-trompete */
        Div1 = ((2.0 * linw) * (2.0 * linw)) / 4.0;
        p2_para = ((sqrt (Div1 + (w12 * w12))) - linw) * 2.0;
        /* calculate heigth of guide exit (w2) */
        w2 = sqrt (p2_para * (l + linw + p2_para / 4.0));
        w2 = w2 * 2.0;
      }
    }
  }
  fprintf (stderr, "Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
  fprintf (stderr, "           Width at the guide exit (w2): %lf \n", w2);
  if (w2 <= 0) {
    fprintf (stderr, "Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
    fprintf (stderr, "           Width at the guide exit (w2) was calculated\n");
    fprintf (stderr, "           <=0; Please change the parameter w1 and/or\n");
    fprintf (stderr, "           l! \n");
    free (fu);
    exit (-1);
  }
  if (!strcmp (fu, "elliptical")) {
    w1_in[0] = w1;
    for (i = 1; i < seg; i++) {
      if (b_ell_q == 0) {
        w1_in[i] = w1;
      } else {
        mxi = (((lbw / 2.0) - linw) - (l_seg * i));
        w1_in[i] = (sqrt ((1.0 - ((mxi * mxi) / a_ell_q)) * b_ell_q)) * 2.0;
      }
      w2_out[i - 1] = w1_in[i];
    }
    w2_out[seg - 1] = w2;
  } else if (!strcmp (fu, "parabolical")) {
    w1_in[0] = w1;
    ii = seg - 1;
    if (loutw == 0 && linw == 0) {
      for (i = 1; i < (seg + 1); i++) {
        w1_in[i] = w1;
        ii = ii - 1;
        w2_out[i - 1] = w1_in[i];
      }
    } else {
      if ((linw == 0) && (loutw > 0)) {
        for (i = 1; i < (seg + 1); i++) {
          w1_in[i] = (sqrt ((p2_para / 4 + loutw + (l_seg * ii)) * p2_para)) * 2;
          ii = ii - 1;
          w2_out[i - 1] = w1_in[i];
        }
      } else {
        for (i = 1; i < (seg + 1); i++) {
          w1_in[i] = (sqrt ((p2_para / 4 + linw + (l_seg * i)) * p2_para)) * 2;
          w2_out[i - 1] = w1_in[i];
        }
      }
    }
  }
  free (fu);
  for (i = 0; i < seg; i++) {
    w1c[i] = w1_in[i];
    w2c[i] = w2_out[i];
    ww[i] = .5 * (w2c[i] - w1c[i]);
    hh[i] = .5 * (h2_out[i] - h1_in[i]);
    whalf[i] = .5 * w1c[i];
    hhalf[i] = .5 * h1_in[i];
    lwhalf[i] = l_seg * whalf[i];
    lhhalf[i] = l_seg * hhalf[i];
  }
  /* guide curvature: rotation angle [rad] between each guide segment */
  if (curvature && l && segno)
    rotation_h = l / curvature / segno;
  if (curvature_v && l && segno)
    rotation_v = l / curvature_v / segno;
%}

TRACE
%{
  double t1, t2, ts, zr;                         /* Intersection times. */
  double av, ah, bv, bh, cv1, cv2, ch1, ch2, dd; /* Intermediate values */
  double vdotn_v1, vdotn_v2, vdotn_h1, vdotn_h2; /* Dot products. */
  int i;                                         /* Which mirror hit? */
  double q;                                      /* Q [1/AA] of reflection */
  double vlen2, nlen2;                           /* Vector lengths squared */
  double edge;
  double hadj; /* Channel displacement */
  int ii;

  /* Propagate neutron to guide entrance. */
  PROP_Z0;
  for (ii = 0; ii < seg; ii++) {
    zr = ii * l_seg;
    /* Propagate neutron to segment entrance. */
    ts = (zr - z) / vz;
    PROP_DT (ts);
    if (x <= w1_in[ii] / -2.0 || x >= w1_in[ii] / 2.0 || y <= -hhalf[ii] || y >= hhalf[ii])
      ABSORB;
    /* Shift origin to center of channel hit (absorb if hit dividing walls) */
    x += w1_in[ii] / 2.0;
    edge = floor (x / w1c[ii]) * w1c[ii];
    if (x - edge > w1c[ii]) {
      x -= w1_in[ii] / 2.0; /* Re-adjust origin */
      ABSORB;
    }
    x -= (edge + (w1c[ii] / 2.0));
    hadj = edge + (w1c[ii] / 2.0) - w1_in[ii] / 2.0;
    for (;;) {
      /* Compute the dot products of v and n for the four mirrors. */
      ts = (zr - z) / vz;
      av = l_seg * vx;
      bv = ww[ii] * vz;
      ah = l_seg * vy;
      bh = hh[ii] * vz;
      vdotn_v1 = bv + av; /* Left vertical */
      vdotn_v2 = bv - av; /* Right vertical */
      vdotn_h1 = bh + ah; /* Lower horizontal */
      vdotn_h2 = bh - ah; /* Upper horizontal */
      /* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */
      cv1 = -whalf[ii] * l_seg - (z - zr) * ww[ii];
      cv2 = x * l_seg;
      ch1 = -hhalf[ii] * l_seg - (z - zr) * hh[ii];
      ch2 = y * l_seg;
      /* Compute intersection times. */
      t1 = (zr + l_seg - z) / vz;
      i = 0;
      if (vdotn_v1 < 0 && (t2 = (cv1 - cv2) / vdotn_v1) < t1) {
        t1 = t2;
        i = 1;
      }
      if (vdotn_v2 < 0 && (t2 = (cv1 + cv2) / vdotn_v2) < t1) {
        t1 = t2;
        i = 2;
      }
      if (vdotn_h1 < 0 && (t2 = (ch1 - ch2) / vdotn_h1) < t1) {
        t1 = t2;
        i = 3;
      }
      if (vdotn_h2 < 0 && (t2 = (ch1 + ch2) / vdotn_h2) < t1) {
        t1 = t2;
        i = 4;
      }
      if (i == 0) {
        break; /* Neutron left guide. */
      }
      PROP_DT (t1);
      switch (i) {
      case 1: /* Left vertical mirror */
        nlen2 = l_seg * l_seg + ww[ii] * ww[ii];
        q = V2Q * (-2) * vdotn_v1 / sqrt (nlen2);
        dd = 2 * vdotn_v1 / nlen2;
        vx = vx - dd * l_seg;
        vz = vz - dd * ww[ii];
        break;
      case 2: /* Right vertical mirror */
        nlen2 = l_seg * l_seg + ww[ii] * ww[ii];
        q = V2Q * (-2) * vdotn_v2 / sqrt (nlen2);
        dd = 2 * vdotn_v2 / nlen2;
        vx = vx + dd * l_seg;
        vz = vz - dd * ww[ii];
        break;
      case 3: /* Lower horizontal mirror */
        nlen2 = l_seg * l_seg + hh[ii] * hh[ii];
        q = V2Q * (-2) * vdotn_h1 / sqrt (nlen2);
        dd = 2 * vdotn_h1 / nlen2;
        vy = vy - dd * l_seg;
        vz = vz - dd * hh[ii];
        break;
      case 4: /* Upper horizontal mirror */
        nlen2 = l_seg * l_seg + hh[ii] * hh[ii];
        q = V2Q * (-2) * vdotn_h2 / sqrt (nlen2);
        dd = 2 * vdotn_h2 / nlen2;
        vy = vy + dd * l_seg;
        vz = vz - dd * hh[ii];
        break;
      }
      /* Now compute reflectivity. */
      if ((i <= 2 && mx == 0) || (i > 2 && my == 0)) {
        x += hadj; /* Re-adjust origin */
        ABSORB;
      } else {
        double ref = 1;
        if (i <= 2) {
          double m = (mx > 0 ? mx : fabs (mx * w1 / w1_in[ii]));
          double par[] = { R0, Qcx, alphax, m, W };
          StdReflecFunc (q, par, &ref);
          if (ref > 0)
            p *= ref;
          else {
            x += hadj; /* Re-adjust origin */
            ABSORB;    /* Cutoff ~ 1E-10 */
          }
        } else {
          double m = (my > 0 ? my : fabs (my * h1 / h1_in[ii]));
          double par[] = { R0, Qcy, alphay, m, W };
          StdReflecFunc (q, par, &ref);
          if (ref > 0)
            p *= ref;
          else {
            x += hadj; /* Re-adjust origin */
            ABSORB;    /* Cutoff ~ 1E-10 */
          }
        }
      }
      x += hadj;
      SCATTER;
      x -= hadj;
    } /* loop on reflections inside segment */
    x += hadj; /* Re-adjust origin */

    /* rotate neutron according to actual guide curvature */
    if (rotation_h) {
      double nvx, nvy, nvz;
      rotate (nvx, nvy, nvz, vx, vy, vz, -rotation_h, 0, 1, 0);
      vx = nvx;
      vy = nvy;
      vz = nvz;
    }
    if (rotation_v) {
      double nvx, nvy, nvz;
      rotate (nvx, nvy, nvz, vx, vy, vz, -rotation_v, 1, 0, 0);
      vx = nvx;
      vy = nvy;
      vz = nvz;
    }
  } /* loop on segments */
%}

FINALLY
%{
  free (w1c);
  free (w2c);
  free (ww);
  free (hh);
  free (whalf);
  free (hhalf);
  free (lwhalf);
  free (lhhalf);
  free (h1_in);
  free (h2_out);
  free (w1_in);
  free (w2_out);
%}

MCDISPLAY
%{
  int i, ii;

  for (ii = 0; ii < seg; ii++) {
    multiline (5, -w1_in[ii] / 2.0, -h1_in[ii] / 2.0, l_seg * (double)ii, -w2_out[ii] / 2.0, -h2_out[ii] / 2.0, l_seg * ((double)ii + 1.0), -w2_out[ii] / 2.0,
               h2_out[ii] / 2.0, l_seg * ((double)ii + 1.0), -w1_in[ii] / 2.0, h1_in[ii] / 2.0, l_seg * (double)ii, -w1_in[ii] / 2.0, -h1_in[ii] / 2.0,
               l_seg * (double)ii);
    multiline (5, w1_in[ii] / 2.0, -h1_in[ii] / 2.0, l_seg * (double)ii, w2_out[ii] / 2.0, -h2_out[ii] / 2.0, l_seg * ((double)ii + 1.0), w2_out[ii] / 2.0,
               h2_out[ii] / 2.0, l_seg * ((double)ii + 1.0), w1_in[ii] / 2.0, h1_in[ii] / 2.0, l_seg * (double)ii, w1_in[ii] / 2.0, -h1_in[ii] / 2.0,
               l_seg * (double)ii);
  }
  line (-w1 / 2.0, -h1 / 2.0, 0.0, w1 / 2.0, -h1 / 2.0, 0.0);
  line (-w1 / 2.0, h1 / 2.0, 0.0, w1 / 2.0, h1 / 2.0, 0.0);
  for (i = 0; i < segno; i++) {
    line (-w2_out[i] / 2.0, -h2_out[i] / 2.0, l_seg * (double)(i + 1), w2_out[i] / 2.0, -h2_out[i] / 2.0, l_seg * (double)(i + 1));
    line (-w2_out[i] / 2.0, h2_out[i] / 2.0, l_seg * (double)(i + 1), w2_out[i] / 2.0, h2_out[i] / 2.0, l_seg * (double)(i + 1));
  }
%}

END