/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Guide_tapering.comp
*
* Component: Guide_tapering
*
* %I
* Written by: Uwe Filges
* Date: 22/10/2003
* Version: $Revision: 1.6 $
* Origin: McStas 1.9.1
*
* 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.
*
* 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 of 1. focal point and real guide entry -
*                 left and right horizontal mirrors
* loutw    [m]    distance of real guide exit and 2nd focal point -
*                 left and right horizontal mirrors
* l:       [m]    length of guide
* linh:    [m]    distance of 1st focal point and real guide entry -
*                 top and bottom vertical mirrors
* louth:   [m]    distance of real guide exit and 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
*                 "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.
* my:      [1]    m-value of material for top and bottom mirrors. Zero
*                 means completely absorbing.
* W:       [AA-1] Width of supermirror cut-off for all mirrors
* segno:   [1]    number of segments (z-axis) for cutting the tube
*
* %End
*
******************************************************************************/

DEFINE COMPONENT Guide_tapering
DEFINITION PARAMETERS (option=0,segno=800)
SETTING PARAMETERS (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)
OUTPUT PARAMETERS(w1c,w2c,ww,hh,whalf,hhalf,lwhalf,lhhalf,h1_in,h2_out,w1_in,w2_out,l_seg)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,sx,sy,p)

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

%}

INITIALIZE
%{

  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;
  if (option != NULL)
  {
     fu = (char*)malloc(sizeof(char)*(strlen(option)+1));
     strcpy(fu,option);
  } else {
     exit(-1);
  }
  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)*((u1 + u2)/2);
        b_ell_q = a_ell_q - ((lbh/2)*(lbh/2));
        /* calculate heigth of guide exit (h2) */
        div1 = ((lbh/2-louth)*(lbh/2-louth))/a_ell_q;
        h2 = sqrt(b_ell_q*(1-div1));
        h2 = h2*2;
     }
  } 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", file_name);
        free(fu);exit(-1);
     }
     if (louth == 0 && linh == 0)
     {
        /* plane mirrors (vertical) */
        h2 = h1;
     } else {
        /* parabolical mirrors */
        if (linh == 0)
        {
           Div1=((2*louth+2*l)*(2*louth+2*l))/4;
           p2_para=((sqrt(Div1+(h12*h12)))-(louth+l))*2;
           /* calculate heigth of guide exit (h2) */
           h2 = sqrt(p2_para*(louth+p2_para/4));
           h2 = h2*2;
         } else {
            /* anti-trompete */
           Div1=((2*linh)*(2*linh))/4;
           p2_para=((sqrt(Div1+(h12*h12)))-linh)*2;
           /* calculate heigth of guide exit (h2) */
           h2 = sqrt(p2_para*(l+linh+p2_para/4));
           h2 = h2*2;
         }
     }
  } 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))
        {
          fgets(segment[ii].st,128,num);
          if (ii >  799) {
             fprintf(stderr,"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 {
     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)-linh) - (l_seg * i));
         h1_in[i] = (sqrt((1-((mxi*mxi)/a_ell_q))*b_ell_q))*2;
       }
     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+louth+(l_seg*ii))*p2_para))*2;
             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+linh+(l_seg*i))*p2_para))*2;
             h2_out[i-1] = h1_in[i];
           }
        }
     }
  }
  /* compute each value for horizontal mirrors */
  w12 = w1/2;
  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 = ((lbw/2)-linw)*((lbw/2)-linw);
       u2 = (lbw/2)*(lbw/2);
       div1 = u1/u2;
       b_ell_q = (w12*w12)/(1-div1);
       u1 = sqrt((linw*linw)+(w12*w12));
       u2 = sqrt((w12*w12) + ((l+loutw)*(l+loutw)));
       a_ell_q = ((u1 + u2)/2)*((u1 + u2)/2);
       b_ell_q = a_ell_q - ((lbw/2)*(lbw/2));
       /* calculate weigth of guide exit (w2) */
       div1 = ((lbw/2-loutw)*(lbw/2-loutw))/a_ell_q;
       w2 = sqrt(b_ell_q*(1-div1));
       w2 = w2*2;
     }
  } 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", file_name);
         free(fu);exit(-1);
     }
     if (loutw == 0 && linw == 0)
     {
        /* plane mirrors (horizontal) */
        w2 = w1;
     } else {
       if (linw == 0)
       {
          /* parabolical mirrors */
          Div1=((2*loutw+2*l)*(2*loutw+2*l))/4;
          p2_para=((sqrt(Div1+(w12*w12)))-(loutw+l))*2;
          /* calculate weigth of guide exit (w2) */
          w2 = sqrt(p2_para*(loutw+p2_para/4));
          w2 = w2*2;
       } else {
          /* anti-trompete */
          Div1=((2*linw)*(2*linw))/4;
          p2_para=((sqrt(Div1+(w12*w12)))-linw)*2;
          /* calculate heigth of guide exit (w2) */
          w2 = sqrt(p2_para*(l+linw+p2_para/4));
          w2 = w2*2;
       }
     }
  }
  fprintf(stderr,"Component: %s (Guide_tapering)\n",NAME_CURRENT_COMP);
  fprintf(stderr,"           Weigth at the guide exit (w2): %lf \n", w2);
  if (w2 <= 0)
  {
   fprintf(stderr,"Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP);
   fprintf(stderr,"           Weigth 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;
     a_ell_q = (lbw/2)*(lbw/2);
     for (i=1;i<seg;i++)
     {
       if (b_ell_q == 0)
       {
         w1_in[i]=w1;
       } else {
         mxi = (((lbw/2)-linw) - (l_seg * i));
         w1_in[i] = (sqrt((1-((mxi*mxi)/a_ell_q))*b_ell_q))*2;
       }
       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];
  }
%}

TRACE
%{
  double t1,t2;                                 /* 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 */
  double sz=0;
  int ii;
  Coords mctc1,mctc2;
  Rotation mctr1;

  for (ii=0;ii<seg;ii++)
  {
    /* Propagate neutron to guide entrance. */
    PROP_Z0;
    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. */
      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*ww[ii]; cv2 = x*l_seg;
      ch1 = -hhalf[ii]*l_seg - z*hh[ii]; ch2 = y*l_seg;
      /* Compute intersection times. */
      t1 = (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;
      }
      if((i <= 2 && q > Qcx) || (i > 2 && q > Qcy))
      {
        if (i <= 2)
        {
          double arg = (q - mx*Qcx)/W;
          if(arg < 10)
            p *= .5*(1-tanh(arg))*(1-alphax*(q-Qcx));
          else
          {
            x += hadj; /* Re-adjust origin */
            ABSORB;                               /* Cutoff ~ 1E-10 */
          }
        } else {
          double arg = (q - my*Qcy)/W;
          if(arg < 10)
            p *= .5*(1-tanh(arg))*(1-alphay*(q-Qcy));
          else
          {
            x += hadj; /* Re-adjust origin */
            ABSORB;                               /* Cutoff ~ 1E-10 */
          }
        }
      }
      p *= R0;
    x += hadj; SCATTER; x -= hadj;
    }
    x += hadj; /* Re-adjust origin */
    /* store current neutron state  */
    if (ii<(seg-1))
    {
        z=z-l_seg;
        STORE_NEUTRON(INDEX_CURRENT_COMP,x,y,z,vx,vy,vz,t,sx,sy,sz,p);
    }
  }
  if (seg > 1)
  {
  z=z+l_seg*(seg-1);
  }

%}

MCDISPLAY
%{
  double x;
  int i,ii;

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

END
