/********************************************************************************************************************************************************
  VITESS module 'chopper_fermi'                                                            
                                                                                           
 The free non-commercial use of these routines is granted providing due credit is given to 
 the authors.                                                                              
                                                                                           
 1.00  Jul 2002  G. Zsigmond	initial version                                             
 1.01  Aug 2002  G. Zsigmond	forward all coordinates                                     
 1.02  Sep 2002  G. Zsigmond	included more channel windows                               
 1.03  Apr 2003  G. Zsigmond	sign correction                                             
 1.04  May 2003  G. Zsigmond	info changed                                                
 1.05  Jun 2003  G. Zsigmond	modulo function included for safety                         
 1.06  Jul 2003  G. Zsigmond	generalised for optional number of pulses; warnings included
 1.07  Oct 2003  G. Zsigmond	superfluous modulo function cancelled                       
 1.08  Nov 2003  G. Zsigmond	put 2 more windows representing channels, now 6 windows    
                               in the big IF loop ">=" changed to ">"                      
 1.09  Jan 2004  K. Lieutenant changes for 'instrument.dat'                                
 1.10  Jan 2004  G. Zsigmond  back to 4 windows representing channels
 1.11  Apr 2004  G. Zsigmond  negative time of flight defined
 1.12  Apr 2004  G. Zsigmond  negative time of flight - corrections, set zero time
 1.13  May 2004  G. Zsigmond  circular geom option and channel length included in curved fc
 1.14  JUL 2004  G. Zsigmond  small change in Init to adapt to new GUI
 1.15  OCT 2004  G. Zsigmond  changed to use both even or odd number of channels
 1.16  MAY 2005  G. Zsigmond  output changed to give trajectory coordinates at a plane crossing the center of the chopper
			      (to be compatible with zero time option	
                               zero time option fixed to get one peak 
                              shadowing cylinder opening activated 
 1.17  MAY 2005  G. Zsigmond  new option choice of 4, 6(better,slower) or 8(much better, very slow) gates, 4 gates option adjusted
 1.18  SEP 2005  G. Zsigmond  optimisations to speed up the algorithm
       Nov 2011  M. Fromem    local declarations of single letter variables
*******************************************************************************************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "general.h"
#include "init.h"
#include "intersection.h"
#include "softabort.h"

/* START HEADER STORY */

#define	STRING_BUFFER 100

int			Option, CurvGeomOption, GatesNumber, zerotime=0;
long		NumOut, Nchannels;
double		TOF, TOF_zero, WL, omega, height, width, depth, optimal_wl, radius_of_curv, main_depth, 
			diameter, Phase, shift_y=0., angle_channel, phase0, wallwidth, expon, pos[3], n[3],
			IntegralIntensity, y_ch[10][2000], x_ch[10][2000], phase[10][2000], coef_pi;
const char	*GeomFileName;
char            XFILEName[STRING_BUFFER];
FILE		*GeomFilePtr, *GatesFilePtr;
VectorType	Pos, Dir, pos_ch, Path;
Neutron		Neutrons;

void		OutputTransformations(double *tof, double *wl, double *prob, VectorType Pos, VectorType Dir, VectorType SpinVector);
void		ReadParameterFile();
void		OwnInit(int argc, char *argv[]);
void		OwnCleanup();

/* FINISH HEADER STORY */


double phase_k_j(double x_ch_k_j, double y_ch_k_j)
{
  double sq_x_ch_k_j, Denom_k, Arg_k, arg_k, pha_k_j, y_ch_new_k_j ;

  double sq_D_0_1, sq_term, omega_fact, dirpos, vz_pos;
 
  sq_D_0_1 = sq(Dir[0]) + sq(Dir[1]);
  dirpos = Dir[0]*Pos[1] - Dir[1]*Pos[0];
  sq_term  = sq(dirpos) / sq_D_0_1;
  omega_fact = omega / (V_FROM_LAMBDA(WL) * Dir[0]);
  vz_pos = Pos[1] > 0.0 ? 1.0 : -1.0;

  sq_x_ch_k_j = sq(x_ch_k_j);
  Denom_k     = sqrt( sq_D_0_1 * (sq_x_ch_k_j + sq(y_ch_k_j)) );

  Arg_k = dirpos / Denom_k;

  if (fabs(Arg_k) > 1.) {
    Arg_k = vz_pos;
    y_ch_new_k_j = Arg_k * sqrt(sq_term - sq_x_ch_k_j);
  } else
    y_ch_new_k_j = y_ch_k_j; /* no intersection with trajectory */

  Denom_k = sqrt( sq_D_0_1 * (sq_x_ch_k_j + sq(y_ch_new_k_j)) );

  arg_k = (Dir[0]*y_ch_new_k_j - Dir[1]*x_ch_k_j) / Denom_k;

  if (fabs(arg_k) > 1.) return  777;
			  
  else {
    pha_k_j = asin(Arg_k) - asin(arg_k); 

    if(x_ch_k_j < 0.) pha_k_j = - pha_k_j; 
								
    return  pha_k_j - omega_fact * (x_ch_k_j * cos(pha_k_j) - y_ch_new_k_j * sin(pha_k_j) - Pos[0]);
  }
}

int searchgates(int gates, double phase0)
{ 
  int j, scs=0;
  int m = -1; 

  phase[0][0] = phase_k_j(x_ch[0][0], y_ch[0][0]);

  for(j=1;j<2*Nchannels+2; j++) {
    phase[0][j] = phase_k_j(x_ch[0][j], y_ch[0][j]);

    if(gates==4){

      if((m == 1)&&(phase[0][j-1] < phase0 )&&(phase0 < phase[0][j]))
        {					  
          phase[1][j-1] = phase_k_j(x_ch[1][j-1], y_ch[1][j-1]); 
          phase[1][j] = phase_k_j(x_ch[1][j], y_ch[1][j]);

          if((phase[1][j-1] < phase0 )&&(phase0 < phase[1][j]))
            {
              phase[2][j-1] = phase_k_j(x_ch[2][j-1], y_ch[2][j-1]); 
              phase[2][j] = phase_k_j(x_ch[2][j], y_ch[2][j]);

              if((phase[2][j-1] > phase0 )&&(phase0 > phase[2][j]))
                {
                  phase[3][j-1] = phase_k_j(x_ch[3][j-1], y_ch[3][j-1]); 
                  phase[3][j] = phase_k_j(x_ch[3][j], y_ch[3][j]);

                  if((phase[3][j-1] > phase0 )&&(phase0 > phase[3][j])) scs = 1;
                }
            }
        } 			
    }
    if(gates==6){

      if((m == 1)&&(phase[0][j-1] < phase0 )&&(phase0 < phase[0][j]))
        {					  
          phase[1][j-1] = phase_k_j(x_ch[1][j-1], y_ch[1][j-1]); 
          phase[1][j] = phase_k_j(x_ch[1][j], y_ch[1][j]);

          if((phase[1][j-1] < phase0 )&&(phase0 < phase[1][j]))
            {
              phase[2][j-1] = phase_k_j(x_ch[2][j-1], y_ch[2][j-1]); 
              phase[2][j] = phase_k_j(x_ch[2][j], y_ch[2][j]);

              if((phase[2][j-1] < phase0 )&&(phase0 < phase[2][j]))
                {
                  phase[3][j-1] = phase_k_j(x_ch[3][j-1], y_ch[3][j-1]); 
                  phase[3][j] = phase_k_j(x_ch[3][j], y_ch[3][j]);

                  if((phase[3][j-1] > phase0 )&&(phase0 > phase[3][j]))
                    {
                      phase[4][j-1] = phase_k_j(x_ch[4][j-1], y_ch[4][j-1]); 
                      phase[4][j] = phase_k_j(x_ch[4][j], y_ch[4][j]);

                      if((phase[4][j-1] > phase0 )&&(phase0 > phase[4][j]))
                        {
                          phase[5][j-1] = phase_k_j(x_ch[5][j-1], y_ch[5][j-1]); 
                          phase[5][j] = phase_k_j(x_ch[5][j], y_ch[5][j]);

                          if((phase[5][j-1] > phase0 )&&(phase0 > phase[5][j])) scs = 1;
                        }
                    }
                }
            }
        } 			
    }
    if(gates==8){

      if((m == 1)&&(phase[0][j-1] < phase0 )&&(phase0 < phase[0][j]))
        {					  
          phase[1][j-1] = phase_k_j(x_ch[1][j-1], y_ch[1][j-1]); 
          phase[1][j] = phase_k_j(x_ch[1][j], y_ch[1][j]);

          if((phase[1][j-1] < phase0 )&&(phase0 < phase[1][j]))
            {
              phase[2][j-1] = phase_k_j(x_ch[2][j-1], y_ch[2][j-1]); 
              phase[2][j] = phase_k_j(x_ch[2][j], y_ch[2][j]);

              if((phase[2][j-1] < phase0 )&&(phase0 < phase[2][j]))
                {
                  phase[3][j-1] = phase_k_j(x_ch[3][j-1], y_ch[3][j-1]); 
                  phase[3][j] = phase_k_j(x_ch[3][j], y_ch[3][j]);

                  if((phase[3][j-1] < phase0 )&&(phase0 < phase[3][j]))
                    {
                      phase[4][j-1] = phase_k_j(x_ch[4][j-1], y_ch[4][j-1]); 
                      phase[4][j] = phase_k_j(x_ch[4][j], y_ch[4][j]);

                      if((phase[4][j-1] > phase0 )&&(phase0 > phase[4][j]))
                        {
                          phase[5][j-1] = phase_k_j(x_ch[5][j-1], y_ch[5][j-1]); 
                          phase[5][j] = phase_k_j(x_ch[5][j], y_ch[5][j]);

                          if((phase[5][j-1] > phase0 )&&(phase0 > phase[5][j]))
                            {
                              phase[6][j-1] = phase_k_j(x_ch[6][j-1], y_ch[6][j-1]); 
                              phase[6][j] = phase_k_j(x_ch[6][j], y_ch[6][j]);

                              if((phase[6][j-1] > phase0 )&&(phase0 > phase[6][j]))
                                {
                                  phase[7][j-1] = phase_k_j(x_ch[7][j-1], y_ch[7][j-1]); 
                                  phase[7][j] = phase_k_j(x_ch[7][j], y_ch[7][j]);

                                  if((phase[7][j-1] > phase0 )&&(phase0 > phase[7][j])) scs = 1;
                                }
                            }
                        }
                    }
                }
            }
        } 			
    }
    m = m * (-1); 
  }
  return scs;
}


int main(int argc, char **argv)
{

  /* Initialize the program according to the parameters given  */

  Init(argc, argv, VT_CHOP_FERMI);
  OwnInit(argc, argv);


  /* Get the neutrons from the file */
  DECLARE_ABORT;

  while((ReadNeutrons())!= 0) {
    int i;
    /* here is what happens to the neutron */

    for(i=0;i<NumNeutGot;i++)	{

      CHECK;

      TOF = InputNeutrons[i].Time;

      WL = InputNeutrons[i].Wavelength;

      CopyVector(InputNeutrons[i].Position, Pos); 

      CopyVector(InputNeutrons[i].Vector, Dir); 

      Dir[0]	= (double)sqrt(1 - sq(Dir[1]) - sq(Dir[2])); 


      /* shift to center of Fermi-Chopper */

      SubVector(Pos, pos_ch);


      /*trajectories which do not intersect the entrance and exit window */


      n[1]=n[2]=0.; n[0]=1.;
      if((PlaneLineIntersect(Pos, Dir, n, - diameter/2., pos))==1)
        {
          if((pos[2]>=height/2.)||(pos[2]<= - height/2.)||(pos[1]>=diameter/2.)||(pos[1]<= - diameter/2.)) goto getlost;
        }
      else goto getlost;

      if((PlaneLineIntersect(Pos, Dir, n, diameter/2., pos))==1)
        {
          if((pos[2]>=height/2.)||(pos[2]<= - height/2.)||(pos[1]>=diameter/2.)||(pos[1]<= - diameter/2.)) goto getlost;
        }
      else goto getlost;
	

      /* translates neutron variables for X'= - diameter/2.  */

      TOF = TOF + (- diameter/2. - Pos[0]) / fabs(Dir[0]) / V_FROM_LAMBDA(WL); 
			
      if (TOF<0 && Nchannels==1) {
        fprintf(LogFilePtr,"\nERROR: Single-slit Fermi chopper needs positive flight time at the chopper position! \n");
        goto my_exit;
      }
				
      CopyVector(Dir, Path);

      MultiplyByScalar(Path, (- diameter/2. - Pos[0])/ Dir[0] );

      AddVector(Pos, Path);  /*	 Path = displacement vector */
							
	
      /* calculate time entering-edge and exiting-edge of gates along the channels */

      phase0 = fmod(Phase + omega*TOF, coef_pi*M_PI); 

      if(searchgates(GatesNumber, phase0)==1) goto happyend;

      /* also tries one turn earlier  */
		  
      if((phase0 > 0)&&(omega > 0)){ 
			  
        phase0 +=  - coef_pi*M_PI;

        if(searchgates(GatesNumber, phase0)==1) goto happyend;
      }

      if((phase0 < 0)&&(omega < 0)){
			  
        phase0 +=  coef_pi*M_PI;

        if(searchgates(GatesNumber, phase0)==1) goto happyend;
      }

      goto getlost;

    happyend:;

      /* Output matters */

      /* transmit coordinates which were not changed, the rest overwrite below */

      Neutrons = InputNeutrons[i]; 

      /* translates neutron variables for output - X'= 0. . */

      Neutrons.Time = TOF + (- Pos[0]) / Dir[0] / V_FROM_LAMBDA(WL); 
			
      if(zerotime==1)
        { 
          Neutrons.Time = fabs(fmod(Neutrons.Time + Phase/omega + coef_pi*M_PI/omega/2., coef_pi*M_PI/omega)) - coef_pi*M_PI/2./omega ;
        }

      CopyVector(Dir, Path);

      MultiplyByScalar(Path, (- Pos[0])/ Dir[0] );

      AddVector(Pos, Path);  /*Path = displacement vector */
		
      CopyVector(Pos, Neutrons.Position);


      /* writes output binary file */

      WriteNeutron(&Neutrons);

    getlost:;
    }
  }

  /* Do the general cleanup */

 my_exit:;

  OwnCleanup();

  fprintf(LogFilePtr," \n");

  Cleanup(pos_ch[0],pos_ch[1],pos_ch[2], 0.0,0.0);	

  return 0;
}


/* own initialization of the module */

void OwnInit(int argc, char *argv[])
{
  const char *intForm = "%d";

  fprintf(LogFilePtr," \n");
  print_module_name("Fermi-Chopper 1.18");

  /*    INPUT  */

  GeomFileName = "chopper_fermi_g.dat";

  CurvGeomOption = 1;

  GatesNumber = 4;

  expon = 0.3333;

  while(argc>1) {
    char *arg = &argv[1][2];

    switch (argv[1][1])	{
    case 'X':
      sscanf(arg, "%lf", &pos_ch[0]);
      break;
      
    case 'Y':
      sscanf(arg, "%lf", &pos_ch[1]);
      break;

    case 'V':
      sscanf(arg, "%lf", &pos_ch[2]);
      break;

    case 'a':
      sscanf(arg, "%lf", &height);
      break;

    case 'b':
      sscanf(arg, "%lf", &width);
      break;

    case 'c':
      sscanf(arg, "%lf", &depth); 
      break;

    case 'L':
      sscanf(arg, "%lf", &optimal_wl);
      break;

    case 'l':
      sscanf(arg, "%ld", &Nchannels);
      break;

    case 'm':
      sscanf(arg, "%lf", &wallwidth);
      break;

    case 'n':
      sscanf(arg, "%lf", &omega);
      if(omega == 0) omega = 1.E-6;
      omega = omega * 2. * M_PI / 1000.;
      break;

    case 'q':
      sscanf(arg, "%lf", &Phase);
      Phase = Phase * M_PI / 180.;
      break;

    case 'O':
      sscanf(arg, intForm, &Option);
      if((Option!=1)&&(Option!=2)){fprintf(LogFilePtr,"\nERROR: Wrong option! Good options: 1-straight, 2-curved \n"); exit(-1); }
      break;

    case 'g':
      sscanf(arg, intForm, &CurvGeomOption);
      break;

    case 'p':
      sscanf(arg, intForm, &GatesNumber);
      if((GatesNumber!=4)&&(GatesNumber!=6)&&(GatesNumber!=8)){fprintf(LogFilePtr,"\nERROR: Good options: 4, 6 or 8 gates  !\n"); exit(-1); }
      break;

    case 'r':
      sscanf(arg, "%lf", &diameter);
      break;

    case 'G':
      GeomFileName=arg;
      break;

    case 'z':
      sscanf(arg, intForm, &zerotime);
      break;


    }
    argc--;
    argv++;
  }

  if(pos_ch[0] < diameter/2.) {fprintf(LogFilePtr,"\nERROR: Minimum position is diameter/2 \n"); exit(-1); }

  if(Nchannels==1) wallwidth=0.;

  GatesFilePtr = fopen("gates.dat","w");


	
  /* calculate edge positions */ 
  
  if(Option==1)
    {
      fprintf(LogFilePtr,"\nStraight Fermi chopper option activated\n");

      /* diameter matter */

      main_depth = 2. * sqrt(sq(diameter/2.) - sq(width/2.));
      if(depth > main_depth) {	fprintf(LogFilePtr,"\nERROR: Diameter too small - not compatible with 'channel length' and 'width'!\nTake min %f cm\n", 2. * sqrt(sq(depth/2.) + sq(width/2.)));	exit(-1);       }

      {
	double add, w_ch;
	long j, k, m;

	if(( w_ch = ( width - (Nchannels + 1) * wallwidth ) / Nchannels ) <= 0.){fprintf(LogFilePtr,"\nERROR: Channel width =< 0 !\n"); exit(-1);}

	fprintf(LogFilePtr,"Channel width: %f cm\nChannels represented by %d gates.\n",w_ch, GatesNumber);

	for(k=0;k<GatesNumber; k++) y_ch[k][0] =  - width/2.;

	m = 1;

	for(j=1; j< 2*Nchannels+2; j++)
	  {
	    if(m == 1) add = wallwidth; else add = w_ch;

	    for(k=0;k<GatesNumber/2; k++) 
	      {
		x_ch[k][j] = - depth/2. * (1. - pow(2.*k/GatesNumber, expon)); x_ch[GatesNumber - 1 - k][j] = depth/2. * (1. - pow(2.*k/GatesNumber, expon));
	      }

	    for(k=0;k<GatesNumber; k++) y_ch[k][j] = y_ch[k][j-1] + add;

	    m = m * (-1);
	  }

	for(k=0;k<GatesNumber/2; k++) 
	  {
	    x_ch[k][0] = x_ch[k][1] = x_ch[k][2*Nchannels] = x_ch[k][2*Nchannels+1] = - depth/2. * (1. - pow(2.*k/GatesNumber, expon)); 
	    x_ch[GatesNumber - 1 - k][0] = x_ch[GatesNumber - 1 - k][1] = x_ch[GatesNumber - 1 - k][2*Nchannels] = x_ch[GatesNumber - 1 - k][2*Nchannels+1] 
	      = depth/2. * (1. - pow(2.*k/GatesNumber, expon));
	  }

	/* activating shadowing cylinder */

	x_ch[0][0] = x_ch[0][1] = x_ch[0][2*Nchannels] = x_ch[0][2*Nchannels+1] = - main_depth/2.;
	
	x_ch[GatesNumber - 1][0] = x_ch[GatesNumber - 1][1] = x_ch[GatesNumber - 1][2*Nchannels] = x_ch[GatesNumber - 1][2*Nchannels+1] = main_depth/2.;
	
		
	/* print out gate coordinates */

	for(j=0; j< 2*Nchannels+2; j++) for(k=0;k<GatesNumber; k++) fprintf(GatesFilePtr,"%f   %f  \n", x_ch[k][j] , y_ch[k][j]);

      }

    }

  if(Option==2)
    {	
      if(CurvGeomOption==1)
	{
	  fprintf(LogFilePtr,"\nCurved Fermi chopper activated \n");

	  fprintf(LogFilePtr,"\nGeometry option: ideally shaped (close to parabolic) long channels ('channel length' inactive parameter) \n");

	  fprintf(LogFilePtr,"\nRadius of curvature (parabolic approximation):\n" 
		  "	%f cm at center\n" 
		  "	%f cm at circumference\n", 
		  V_FROM_LAMBDA(optimal_wl)/2./omega, V_FROM_LAMBDA(optimal_wl)/2./omega*pow((1 + sq(omega*diameter/V_FROM_LAMBDA(optimal_wl))), 1.5));


	  GeomFilePtr = fopen(GeomFileName,"w");
	  {
	    double add, w_ch, xx[500], yy[500], tt;

	    long j, m, k;

	    if((     w_ch = ( width - (Nchannels + 1) * wallwidth ) / Nchannels    ) <= 0.)
	      {fprintf(LogFilePtr,"\nChannel width =< 0 !\n"); exit(-1);}

	    fprintf(LogFilePtr,"Channel width: %f cm\nChannels represented by %d gates.\n",w_ch, GatesNumber);


	    y_ch[0][0] =  - width/2.;
	    x_ch[0][0] = sqrt(sq(diameter/2.) - sq(y_ch[0][0]));

	    /* the circle */for(k=0; k<500; k++)
	      { tt = k /500. * 2 * M_PI; xx[k] = diameter/2. * cos(tt); yy[k] = diameter/2. * sin(tt); 
		fprintf(GeomFilePtr, " %d %lf  %lf\n", 0, xx[k], yy[k]);
	      }

	    m = 1;

	    for(j=1; j< 2*Nchannels+2; j++)
	      {
		if(m == 1) add = wallwidth; else add = w_ch;
							
		y_ch[0][j] = y_ch[0][j-1] + add;

		x_ch[0][j] = - sqrt(sq(diameter/2.) - sq(y_ch[0][j]));

		m = m * (-1);
	      }

	    for(j=0; j< 2*Nchannels+2; j++)
	      {							
		for(k=0; k<500; k++)
		  {
		    tt= k /499. * 2. * fabs(x_ch[0][j])/V_FROM_LAMBDA(optimal_wl);
		    xx[k] = y_ch[0][j] * sin(omega*tt) + (V_FROM_LAMBDA(optimal_wl)*tt - fabs(x_ch[0][j])) * cos(omega*tt);
		    yy[k] = y_ch[0][j] * cos(omega*tt) - (V_FROM_LAMBDA(optimal_wl)*tt - fabs(x_ch[0][j])) * sin(omega*tt);

		    fprintf(GeomFilePtr, "%ld   %lf  %lf\n", j, xx[k], yy[k]);		
		  }

		for(k=0;k<GatesNumber/2; k++) 
		  { int index;
										
		    index = (int)(pow(2.*k/GatesNumber, expon)*250.);
		    x_ch[k][j] = xx[index];
		    y_ch[k][j] = yy[index];
		    x_ch[GatesNumber - 1 - k][j] = xx[499 - index];
		    y_ch[GatesNumber - 1 - k][j] = yy[499 - index];
		  }
	      }
	    /* print out gate coordinates */

	    for(j=0; j< 2*Nchannels+2; j++) for(k=0;k<GatesNumber; k++) fprintf(GatesFilePtr,"%f   %f  \n", x_ch[k][j] , y_ch[k][j]);
	  }

	}
      if(CurvGeomOption==2)
	{
	  /* diameter matter */

	  main_depth = 2. * sqrt(sq(diameter/2.) - sq(width/2.));
	  if(depth > main_depth) {
	    fprintf(LogFilePtr,"\nERROR: Diameter too small - not compatible with 'channel length' and 'width'!\nTake min %f cm\n",
		    2. * sqrt(sq(depth/2.) + sq(width/2.))); exit(-1);}

		  
	  fprintf(LogFilePtr,"\nCurved Fermi chopper activated \n");

	  fprintf(LogFilePtr,"\nGeometry option: circular shaped channels with fixed length (via 'channel length') \n");

	  radius_of_curv = V_FROM_LAMBDA(optimal_wl)/2./omega;

	  angle_channel = atan(depth/2./radius_of_curv);

	  fprintf(LogFilePtr,"\nRadius of curvature (parabolic approximation):\n" 
		  "optimal_velocity/2./omega = %f cm \n" 
		  "Angle channel: %f deg \n", 
		  radius_of_curv, 180./M_PI*angle_channel);


	  GeomFilePtr = fopen(GeomFileName,"w");
	  {
	    double add, w_ch, xx[500], yy[500], tt;

	    long j, m, k;

	    if((     w_ch = ( width - (Nchannels + 1) * wallwidth ) / Nchannels    ) <= 0.)
	      {fprintf(LogFilePtr,"\nChannel width =< 0 !\n"); exit(-1);}

	    fprintf(LogFilePtr,"Channel width: %f cm\nChannels represented by %d gates.\n",w_ch, GatesNumber);

	    for(k=0; k<500; k++)
	      { tt = k /500. * 2 * M_PI; xx[k] = diameter/2. * cos(tt); yy[k] = diameter/2. * sin(tt); /* the circle */
		fprintf(GeomFilePtr, " %d %lf  %lf\n", -1, xx[k], yy[k]);
	      }

	    m = 1;

	    for(j=0; j< 2*Nchannels+2; j++)
	      {
		if(m == 1) add = wallwidth; else add = w_ch;
								
		for(k=0; k<500; k++)
		  {
		    xx[k]= (2.*k /499.  - 1.) * depth/2.;
		    yy[k] = - width/2. - sqrt(sq(radius_of_curv) - sq(depth/2.)) + sqrt(sq(radius_of_curv) - sq(xx[k])) + shift_y;

		    fprintf(GeomFilePtr, "%ld   %lf  %lf\n", j, xx[k], yy[k]);		
		  }

		for(k=0;k<GatesNumber/2; k++) 
		  { int index;
							
		    index = (int)(pow(2.*k/GatesNumber, expon)*250.);
		    x_ch[k][j] = xx[index];
		    y_ch[k][j] = yy[index];
		    x_ch[GatesNumber - 1 - k][j] = xx[499 - index];
		    y_ch[GatesNumber - 1 - k][j] = yy[499 - index];
		  }													  
		shift_y += add;
		m = m * (-1);
	      }
	    /* activating shadowing cylinder */

	    x_ch[0][0] = x_ch[0][1] = x_ch[0][2*Nchannels] = x_ch[0][2*Nchannels+1] = - main_depth/2.;
				
	    x_ch[GatesNumber - 1][0] =
	      x_ch[GatesNumber - 1][1] =
	      x_ch[GatesNumber - 1][2*Nchannels] =
	      x_ch[GatesNumber - 1][2*Nchannels+1] =
	      main_depth/2.;
	
	    /* print out gate coordinates */

	    for(j=0; j< 2*Nchannels+2; j++) for(k=0;k<GatesNumber; k++) fprintf(GatesFilePtr,"%f   %f  \n", x_ch[k][j] , y_ch[k][j]);

	  }
	}
    }
		
  if(Option==1) coef_pi=1.; else coef_pi=2.;  fprintf(LogFilePtr,"Phase set is %f°.\n", 180./M_PI*fmod(Phase , coef_pi*M_PI));   


  /*	 init for ASCII output */

  NumOut=0;

  IntegralIntensity = 0.;

}/* End OwnInit */



void OwnCleanup()
{
  if (GeomFilePtr) fclose(GeomFilePtr);
  if (GatesFilePtr) fclose(GatesFilePtr);
}/* End OwnCleanup */
