/*******************************************************************************
*
* 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;
        while (!feof(num))
        {
          char *ret = fgets(segment[ii].st,128,num);
          if (ii >  799 || !ret) {
             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 < segno; 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