/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: multi_pipe
*
* %I
* Written by: Uwe Filges
* Date: March, 2005
*
* multi pipe circular slit.
*
* %D
* No transmission around the slit is allowed.
*
* example for an input file
* # user defined geometry
* # x(i) y(i) r(i) [m]
*   0.02 0.03 0.01
*  -0.04 0.015 0.005
*
*  warning: at least two values must be in the file
*
* Example: multi_pipe(xmin=-0.01, xmax=0.01, ymin=-0.01, ymax=0.01)
*
* %P
* INPUT PARAMETERS
*
* radius: []       radius of a single hole
* xmin: [m]        Lower x bound 
* xmax: [m]        Upper x bound 
* ymin: [m]        Lower y bound 
* ymax: [m]        Upper y bound 
* xwidth: [m]      Width of slit plate. Overrides xmin,xmax. 
* yheight: [m]     Height of slit plate. Overrides ymin,ymax. 
* gap: [m]         distance between holes 
* thickness: [m]   thickness of the pipe 
* filename: [str]  define user table of holes 
*
* %E
*******************************************************************************/

DEFINE COMPONENT multi_pipe

SETTING PARAMETERS (string filename=0, xmin,xmax,ymin,ymax,radius=0.0,gap=0.0,thickness=0.0, xwidth=0, yheight=0)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "read_table-lib"
%}

DECLARE
%{
t_Table pTable;
int user_file;
int ynum;
int xnum;
double diam;
double w_pipe;
double h_pipe;
int h_number;
%}


INITIALIZE
%{

char file_name[1024];
char *pos;
char *fu;
int  check;

if (xwidth  > 0) { xmax = xwidth/2;  xmin = -xmax; }
  if (yheight > 0) { ymax = yheight/2; ymin = -ymax; }

  if (xmin == 0 && xmax == 0 && ymin == 0 && ymax == 0 && radius == 0 && gap==0)
    { fprintf(stderr,"multi_pipe: %s: Error: give geometry\n", NAME_CURRENT_COMP); exit(-1); }
  user_file=0;
  if (filename != NULL) {
     fu=(char*)malloc(sizeof(char)*(strlen(filename)+1));
     strcpy(fu,filename);

     user_file=1;

     Table_Read(&pTable, fu, 1);  /* read 1st block data from file into pTable */
     if (pTable.rows < 2) Table_Free(&pTable);
     Table_Info(pTable);
     free(fu);

  } else {
    /* no user file is set */
    /* calculate the number of holes */
    diam=2*radius;
    w_pipe=xmax-xmin;
    h_pipe=ymax-ymin;
    xnum = (int) (w_pipe/(diam+gap));
    ynum = (int) (h_pipe/(diam+gap));
    h_number= xnum*ynum;
    fprintf(stderr,"multi_pipe-hole number: %d\n", h_number);
    user_file=2;
  }
%}

TRACE
%{
    int check1,check2;
    double test_pipes;
    double test_pipe2;
    double xl,yl,dt0;
    int i,ii;
    int valuex,valuey;
    int even_x,even_y;

    PROP_Z0;
    if (x<xmin || x>xmax || y<ymin || y>ymax) {
       /* neutrons outside of component */
       ABSORB;
    } else {
      if (user_file==1) {
         /*  use file data */
         if (pTable.data != NULL)
         {
             double uradius;

             check1=0;
             check2=0;

             for (i=0;i<pTable.rows;i++) {
                xl = Table_Index(pTable, i,   0); /* 1st column */
                yl = Table_Index(pTable, i,   1); /* 2nd column */
                uradius = Table_Index(pTable, i,   2); /* 3rd column */
                test_pipes=((x-xl)*(x-xl))+((y-yl)*(y-yl));
                if (test_pipes < (uradius*uradius)) {
                   check1=1;
                   if (thickness==0) {
                      /* nothing to do */
                   } else {
                      /* check if neutron absorb inside the component */
                      dt0=thickness/vz;
                      PROP_DT(dt0);
                      check2=0;
                      test_pipe2=((x-xl)*(x-xl))+((y-yl)*(y-yl));
                      if (test_pipe2 > (uradius*uradius)) {
                         check2=1;
                         ABSORB;
                      }

                   }
                   break;
                }
             }
         }
      } else if (user_file==2) {
         /* automatic generated holes */
         /* inside component - proof whether the neutrons come into a hole */
         valuex=0;
         even_x = (xnum % 2);
         /* check if xnum is odd or even */
         if (even_x==0) {
            valuex=1;
         } else {
            valuex=2;
         }
         valuey=0;
         even_y = (ynum % 2);
         /* check if ynum is odd or even */
         if (even_y==0) {
            /* even */
            valuey=1;
         } else {
            valuey=2;
         }
         test_pipes=0;
         check1=0;
         for(i=-(xnum/2);i<(xnum/2+1);i++)
         {
            /* calculate local coordiantes */
            if (valuex==1) {
               if (i==xnum/2) {
                  /* outside */
               } else {
                  xl=i*(diam+gap)+radius+gap/2.0;
               }
             } else {
               xl=i*(diam+gap);
             }
             /* fprintf(stderr,"xl: %f\n", xl); */
             for (ii=-(ynum/2);ii<(ynum/2+1);ii++)
             {
                if (valuey==1) {
                   if (ii==ynum/2) {
                      /* ouside */
                   } else {
                      yl=ii*(diam+gap)+radius+gap/2;
                      test_pipes=((x-xl)*(x-xl))+((y-yl)*(y-yl));
                      if (test_pipes < (radius*radius)) {
                         check1=1;
                         if (thickness==0) {
                            /* nothing to do */
                         } else {
                            /* check if neutron absorb inside the component */
                            dt0=thickness/vz;
                            PROP_DT(dt0);
                            check2=0;
                            test_pipe2=((x-xl)*(x-xl))+((y-yl)*(y-yl));
                            if (test_pipe2 > (radius*radius)) {
                               check2=1;
                               ABSORB;
                            }
                         }
                         break;
                      }
                   }
                } else {
                   yl=ii*(diam+gap);
                   test_pipes=((x-xl)*(x-xl))+((y-yl)*(y-yl));
                   if (test_pipes < (radius*radius)) {
                      check1=1;
                      if (thickness==0) {
                         /* nothing to do */
                      } else {
                         /* check if neutron absorb inside the component */
                         dt0=thickness/vz;
                         PROP_DT(dt0);
                         check2=0;
                         test_pipe2=((x-xl)*(x-xl))+((y-yl)*(y-yl));
                         if (test_pipe2 > (radius*radius)) {
                            check2=1;
                            ABSORB;
                         }
                      }
                      break;
                   }
                }
             }
         }
      }
      if (check1==1) {
         SCATTER;
      } else {
         ABSORB;
      }
   }
%}

FINALLY
%{
  Table_Free(&pTable);
%}

MCDISPLAY
%{
  int even_x,even_y;
  int valuex,valuey;
  double xl,yl;
  int i,ii;
  t_Table plot_table;
  double uradius;
  char *fu2;

  
  multiline(5,
            xmin, ymin, 0.0,
            xmax, ymin, 0.0,
            xmax, ymax, 0.0,
            xmin, ymax, 0.0,
            xmin, ymin, 0.0);
  if (thickness > 0) {
     multiline(5,
            xmin, ymin, (double)thickness,
            xmax, ymin, (double)thickness,
            xmax, ymax, (double)thickness,
            xmin, ymax, (double)thickness,
            xmin, ymin, (double)thickness);
     line(xmin, ymin, 0, xmin, ymin, (double)thickness);
     line(xmax, ymin, 0, xmax, ymin, (double)thickness);
     line(xmax, xmax, 0, xmax, ymax, (double)thickness);
     line(xmin, xmax, 0, xmin, ymax, (double)thickness);
  }

  if (user_file==1) {
     /* open user file */
      if (filename != NULL) {
         fu2=(char*)malloc(sizeof(char)*(strlen(filename)+1));
         strcpy(fu2,filename);
         Table_Read(&plot_table, fu2, 1);  /* read 1st block data from file into pTable */
         if (plot_table.rows < 2) Table_Free(&plot_table);
         Table_Info(plot_table);
         free(fu2);
          for (i=0;i<plot_table.rows;i++) {
             xl = Table_Index(plot_table, i,   0); /* 1st column */
             yl = Table_Index(plot_table, i,   1); /* 2nd column */
             uradius = Table_Index(plot_table, i,   2); /* 3rd column */
             circle("xy", xl, yl, 0, uradius);
             if (thickness > 0 ) {
               circle("xy", xl, yl, (double)thickness, uradius);
             }
          }
       }
     Table_Free(&plot_table);

   } else if (user_file==2) {
      valuex=0;
      even_x = (xnum % 2);
      /* check if xnum is odd or even */
      if (even_x==0) {
         valuex=1;
      } else {
         valuex=2;
      }
      valuey=0;
      even_y = (ynum % 2);
      /* check if ynum is odd or even */
      if (even_y==0) {
         /* even */
         valuey=1;
      } else {
         valuey=2;
      }
      for(i=-(xnum/2);i<(xnum/2+1);i++)
      {
         /* calculate local coordiantes */
         if (valuex==1) {
            if (i==xnum/2)
            {
               /* outside */
            } else {
              xl=i*(2*radius+gap)+radius+gap/2.0;
            }
          } else {
              xl=i*(2*radius+gap);
          }
          for (ii=-(ynum/2);ii<(ynum/2+1);ii++)
          {
             if (valuey==1) {
                if (ii==ynum/2) {
                   /* ouside */
                } else {
                   yl=ii*(2*radius+gap)+radius+gap/2;
                   circle("xy", xl, yl, 0, radius);
                   if (thickness > 0 ) {
                      circle("xy", xl, yl, (double)thickness, radius);
                   }
                }
             } else {
                yl=ii*(2*radius+gap);
                circle("xy", xl, yl, 0, radius);
                if (thickness > 0 ) {
                   circle("xy", xl, yl, (double)thickness, radius);
                }
             }
          }
      }
   }
%}

END