/************************************************************************************************/
/*  VITESS module 'chopper_fermi'                                                               */
/*                                                                                              */
/* The free non-commercial use of these routines is granted providing due credit is given to    */
/* the authors.                                                                                 */
/*                                                                                              */
/* 0.1  Jul 2002  G. Zsigmond	  initial version                                               */
/* 1.0  Sep 2002  G. Zsigmond	  active on both sides                                          */
/* 1.1  Apr 2003  S. Manoshin	  visualisation included                                        */
/*              + G. Zsigmond	   and position filter                                          */
/* 1.2  May 2003  G. Zsigmond	  collision file name as input parameter, coll. output changed, */
/*                              quantisation direction incl.                                    */
/* 1.3  Jun 2003  G. Zsigmond	  included if when closing coll file; softabort included        */
/* 1.4  Jan 2004  K. Lieutenant changes for 'instrument.dat'                                    */
/* 1.5  Feb 2004  S. Manoshin   Visualise only first 10000 neutrons                             */
/*                              Choose the output device : screen, file or both                 */
/*                              New external variable gselec                                    */
/* 1.6  Apr 2004  G. Zsigmond	  Visualise only first 1000 neutrons, write sm_ensemble.ps      */
/* 1.7  Feb 2008  K. Lieutenant extension to MAX_MIRR (=13) mirrors                             */
/* 1.8  Apr 2008  M. Fromme changeable constant MAX_MIRR                                        */
/************************************************************************************************/

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

#ifdef VT_GRAPH
 #include "cpgplot.h"

/* MF: visualisation for Windows and generation of file for the picture */
 #ifdef DO_WIN32 
  #define GDEV "sm_ensemble.ps"
 #else 
  #define GDEV "sm_ensemble.png"
 #endif 

 extern int gselec;
 // choose the output 1 - display only, 2 - file only,
 //                   3 - both, defined in cpgplot.c
#else
 int gselec;
#endif

/* START HEADER STORY */

#define MAX_MIRR 25

#define	STRING_BUFFER 1000

FILE	   *Par_Field, *COLLFILE;
char       Option[STRING_BUFFER], *ParameterFileName, *ReflUpFileName, *ReflDownFileName,
           *COLLFILEName="collision.dat";
int	   p=0, datanumber, vistype=0, quant_dir=2;
int        max_mirr; // highest number of used mirror
long	   User, NumWrong, Repetition, repet, nocol, nocolM = 10000, Wallonoff, NoCh, ok;
long	   number_vis_tr=0; /* current number of visualised trajectories, add SM */
double	   TOF, WL, Prob, phi, the;
double	   rupdata[1001], rdowndata[1001], OutputAngleHoriz, OutputAngleVert, RotMatrixOut[3][3];
double	   IntegralIntensity, Path, Path0;
double     prob[MAX_MIRR+1];
	
double	   Windw=-10.0, WindW=200.0, Windh=-10.0, WindH=10.0, wei_min1=0.0;
	
VectorType Pos, Dir, SpinVector, n, TranslOutput;
Neutron	   Neutrons;

VectorType WallOffset[MAX_MIRR+1], WallNormal[MAX_MIRR+1], r1[MAX_MIRR+1], r2[MAX_MIRR+1],
           r3[MAX_MIRR+1], r4[MAX_MIRR+1];
double	   RotMatrixWall[MAX_MIRR+1][3][3], WallVert[MAX_MIRR+1], WallHoriz[MAX_MIRR+1], PathA[MAX_MIRR+1],
           thetaC[MAX_MIRR+1][2], thetaCSM[MAX_MIRR+1][2], RthetaCSM[MAX_MIRR+1][2], mued[MAX_MIRR+1][4],
           mrangh[MAX_MIRR+1], mrangv[MAX_MIRR+1]; 

int	hittriangle(VectorType r1, VectorType r2, VectorType rt);
int	hitwall(VectorType r1, VectorType r2, VectorType r3, VectorType r4, VectorType rt);
double	CollideWall(double *prob, VectorType pos,VectorType dir,VectorType spin,
		    VectorType WallOffset,VectorType WallNormal,
		    double RotMatrixWall[3][3],VectorType r1,VectorType r2,VectorType r3,VectorType r4,
		    double thetaC[2],
		    double thetaCSM[2], double RthetaCSM[2], double mued[4], double mrangh, double mrangv);
double	reflectivity(double angle, double CritUp, double WL, double smearf, double ReflUp);
void	OutputTransformations(double *tof, double *wl, double *prob, VectorType Pos, VectorType Dir,
			      VectorType SpinVector);
void	ReadParameterFile();
void	OwnInit(int argc, char *argv[]);
void	OwnCleanup();
void	CartesianToSpherical2(VectorType Vector, double *Theta, double *Phi);

int     useQuantDir=1;

/* FINISH HEADER STORY */

#ifdef VT_GRAPH
static void moveIt(VectorType P) {
  switch (p) {
  case 2:
    cpgmove((float) P[0], (float) P[1]); 
    cpgpt1 ((float) P[0], (float) P[1], -2);
    break;
  case 3:
    cpgmove((float) P[0], (float) P[2]); 
    cpgpt1 ((float) P[0], (float) P[2], -2);
    break;
  case 4:
    cpgmove((float) P[1], (float) P[2]);
    cpgpt1 ((float) P[1], (float) P[2], -2);
  }
}

static void drawIt(VectorType P) {
  switch (p) {
  case 2:
    if (vistype) cpgpt1 ((float) P[0], (float) P[1], -2);
    else         cpgdraw((float) P[0], (float) P[1]); 
    break;
  case 3:
    if (vistype) cpgpt1 ((float) P[0], (float) P[2], -2);
    else         cpgdraw((float) P[0], (float) P[2]); 
    break;
  case 4:
    if (vistype) cpgpt1 ((float) P[1], (float) P[2], -2);
    else         cpgdraw((float) P[1], (float) P[2]);
  }
}
#endif

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

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

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

  if (p==1) {
    COLLFILE = fopen(COLLFILEName, "w");
    fprintf(COLLFILE, "     ID      debug color          no  sp wall         x/cm          y/cm          z/cm           dir y/°       dir z/° \n\n");
  }

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

  while (ReadNeutrons()) {
    int i;

    for (i=0; i < NumNeutGot; i++) { 
      int j,m;
      CHECK;
      
      if (number_vis_tr == 1000)
      	p = 100;

      InputNeutrons[i].Vector[0] = (double) sqrt(1 - sq(InputNeutrons[i].Vector[1]) - sq(InputNeutrons[i].Vector[2]));

      TOF  = InputNeutrons[i].Time;
      WL   = InputNeutrons[i].Wavelength;
      Prob = InputNeutrons[i].Probability;

      CopyVector(InputNeutrons[i].Position, Pos);
      CopyVector(InputNeutrons[i].Vector, Dir);
      CopyVector(InputNeutrons[i].Spin, SpinVector); 

      /*  compute hit positions on the walls */

      Path0 = 0.; m = nocol = 0;
      {
	VectorType pos[MAX_MIRR+1], dir[MAX_MIRR+1], spin[MAX_MIRR+1]; 
 
        if (p==1)
	  fprintf(COLLFILE,
		  "     %c%c%07ld %c %5d    %10d  %2d  0     %12.5f  %12.5f  %12.5f     %12.5f  %12.5f \n",
		  InputNeutrons[i].ID.IDGrp[0], InputNeutrons[i].ID.IDGrp[1], InputNeutrons[i].ID.IDNo,
		  InputNeutrons[i].Debug, InputNeutrons[i].Color, i, ((int) SpinVector[quant_dir]),
		  Pos[0], Pos[1], Pos[2], 180./M_PI * atan2(Dir[1],Dir[0]), 180./M_PI * atan2(Dir[2],Dir[0]));

#ifdef VT_GRAPH
        else if (p>=2 && p<=4) {
	  cpgsci((int)(InputNeutrons[i].Color));
	  if(Prob > wei_min1)
	    moveIt(Pos);
	}
#endif

        for (j=0; j<1000; j++) { // loop over at most 1000 collisions
	  int im, l;
          for (l=1; l<=max_mirr; l++) {
            CopyVector(Pos, pos[l]);
	    CopyVector(Dir, dir[l]);
	    CopyVector(SpinVector, spin[l]);
	    prob[l]= Prob; 
            if (m != l)
	      PathA[l] = CollideWall(&prob[l], pos[l], dir[l], spin[l], WallOffset[l], WallNormal[l],
				     RotMatrixWall[l], r1[l], r2[l], r3[l], r4[l], thetaC[l], thetaCSM[l],
				     RthetaCSM[l], mued[l], mrangh[l], mrangv[l]);
	    else
	      PathA[l] = 99999;
          }

	  for (im=1; im <= max_mirr; im++) 
	    if (PathA[im] != 99999.0)
	      break;
	  if (im > max_mirr) {
	    // all PathA are 99999
	    Path = 99999.0;
	    goto conti;
          }

          for (l=1; l<=max_mirr; l++) { // loop over mirrors
	    if (m == l) continue;
	    for (im=1; im<=max_mirr; im++) 
	      if (im != l && PathA[l] > PathA[im])
		break;
	    if (im <= max_mirr) continue; // because PathA[l] > PathA[im]

	    CopyVector(pos[l], Pos); CopyVector(dir[l], Dir); Path = PathA[l]; Prob = prob[l]; 
	    m=l; nocol +=1; 
	    if (p==1)
	      fprintf(COLLFILE,
		      "     %c%c%07ld %c %5d    %10d  %2d  %d     %12.5f  %12.5f  %12.5f     %12.5f  %12.5f \n", 
		      InputNeutrons[i].ID.IDGrp[0], InputNeutrons[i].ID.IDGrp[1],
		      InputNeutrons[i].ID.IDNo, InputNeutrons[i].Debug,  InputNeutrons[i].Color, i,
		      ((int) SpinVector[quant_dir]), 
		      m, Pos[0], Pos[1], Pos[2], 180./M_PI * atan2(Dir[1],Dir[0]), 180./M_PI * atan2(Dir[2],Dir[0]));
	    
#ifdef VT_GRAPH
	    else if(Prob > wei_min1)
	      drawIt(Pos);
#endif
	    
          }  // end loop over mirrors
	  
          if(nocol == nocolM) goto conti; 
	
          if(Prob < wei_min) goto getlost; 

          Path0 += Path;

        }  // end loop over collisions
conti:;

      }    // end of range for local definitions
      if(Prob < wei_min) goto getlost; 


      /* transform into output frame  */
	
      SubVector(Pos, TranslOutput);
      RotVector(RotMatrixOut, Pos);
      RotVector(RotMatrixOut, Dir);

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

      {
        VectorType Path ;
      
        if (Pos[0] > 0.0) goto getlost; /*	Filter if collision after output YZ plane */
      
        TOF += - Pos[0] / fabs(Dir[0]) / V_FROM_LAMBDA(WL) ;
      	
        CopyVector(Dir, Path) ;
      
        MultiplyByScalar(Path, - Pos[0]/ Dir[0] ) ;
      
        AddVector(Pos, Path) ;  
	
      }			/* Path = displacement vector */

      /* transform into input frame to give the exit positions here */

      if (p>0) {
        VectorType posex, direx ;
        CopyVector(Pos, posex) ;
        CopyVector(Dir, direx) ;
        RotBackVector(RotMatrixOut, direx);
        RotBackVector(RotMatrixOut, posex);
        AddVector(posex, TranslOutput);

	if (p == 1) 
	  fprintf(COLLFILE,
		  "     %c%c%07ld %c %5d    %10d  %2d  0     %12.5f  %12.5f  %12.5f     %12.5f  %12.5f \n",
		  InputNeutrons[i].ID.IDGrp[0], InputNeutrons[i].ID.IDGrp[1], InputNeutrons[i].ID.IDNo,          
		  InputNeutrons[i].Debug,       InputNeutrons[i].Color, i, ((int) SpinVector[quant_dir]),
		  posex[0], posex[1], posex[2],
		  180./M_PI * atan2(direx[1],direx[0]), 180./M_PI * atan2(direx[2], direx[0]));
#ifdef VT_GRAPH
        else if (Prob > wei_min1)
	  drawIt(posex);
#endif
      }
	
      /* transmit coordinates which were not changed, overwrite the rest below */
      Neutrons = InputNeutrons[i]; 
      
      Neutrons.Time = TOF + Path0 / V_FROM_LAMBDA(WL);
      Neutrons.Probability = Prob;
      
      CopyVector(Pos, Neutrons.Position);
      CopyVector(Dir, Neutrons.Vector);
      
      Neutrons.Vector[0]	= (double) sqrt(1 - sq(Neutrons.Vector[1]) - sq(Neutrons.Vector[2]));
      
      CopyVector(SpinVector, Neutrons.Spin);

      number_vis_tr = number_vis_tr + 1;

      /* write output binary file ps(nocol);*/
      WriteNeutron(&Neutrons);  

getlost:;
      if(p==1) fprintf(COLLFILE, " \n");

    }  // end loop over trajectories
  }
   
my_exit:;

  /* Do the cleanup */

  if(p==1) fclose(COLLFILE);	
	
  OwnCleanup(); 

  Cleanup(TranslOutput[0], TranslOutput[1], TranslOutput[2], OutputAngleHoriz,OutputAngleVert);	

  return 0;
}


/* own initialization of the ensemble module */

void OwnInit(int argc, char *argv[])
{
  int j;
  fprintf(LogFilePtr," \n");
  print_module_name("supermirror_ensemble 1.8");

  for(j=0;j<3;j++) TranslOutput[j]=0.;
  OutputAngleHoriz=OutputAngleVert=0.;

  gselec = 1 ; /* Activate visualisation device -screen */

  while(argc>1) {
    switch(argv[1][1]) {
				
    case 'P':
      if((Par_Field = fopen(&argv[1][2],"r"))==NULL) {
	fprintf(LogFilePtr,"\nERROR: Parameter file '%s' not found.\n",&argv[1][2]);
	exit(0);
      }
      ParameterFileName=&argv[1][2];
      break;
		
    case 'C':
      COLLFILEName=&argv[1][2];
      break;
		
    case 'M':
      sscanf(&argv[1][2], "%ld", &nocolM); fprintf(LogFilePtr,"Stops at %ld collisions.\n", nocolM);
      break;

    case 'T':
      sscanf(&argv[1][2], "%d", &p); 
      break;

    case 'Q':
      sscanf(&argv[1][2], "%d", &quant_dir); 
      if ((quant_dir !=0)&&(quant_dir !=1)&&(quant_dir !=2)&&(quant_dir !=-1)) {
	fprintf(LogFilePtr,"\nERROR:  wrong quantization direction definition. \n\n");
	exit(0);
      }
      else if (quant_dir == -1) {
	useQuantDir = 0;
      }
      break;

    case 'r':
      sscanf(&argv[1][2], "%lf", &TranslOutput[0]); 
      break;

    case 's':
      sscanf(&argv[1][2], "%lf", &TranslOutput[1]); 
      break;

    case 't':
      sscanf(&argv[1][2], "%lf", &TranslOutput[2]); 
      break;

    case 'h':
      sscanf(&argv[1][2], "%lf", &OutputAngleHoriz); 
      break;

    case 'v':
      sscanf(&argv[1][2], "%lf", &OutputAngleVert); 
      break;

    /* Visual data */					
	
    case 'w':
      sscanf(&argv[1][2], "%lf", &Windw); 
      break;
					
    case 'W':
      sscanf(&argv[1][2], "%lf", &WindW); 
      break;
										
    case 'a':
      sscanf(&argv[1][2], "%lf", &Windh); 
      break;
					
    case 'A':
      sscanf(&argv[1][2], "%lf", &WindH); 
      break;
					
    case 'b':
      sscanf(&argv[1][2], "%lf", &wei_min1); 
      break;
			
    case 'c':
      sscanf(&argv[1][2], "%d", &vistype); 
      break;

    case 'o':
      gselec = atol(&argv[1][2]);
      break;
    }
    argc--;
    argv++;
  }
	

  if (p==1)  
    fprintf(LogFilePtr,"Prints the coordinates of collisions to '%s'. \n", COLLFILEName);
	
#ifdef VT_GRAPH
  else if (p >= 2 && p <= 4) {
	
    /* initialize pgplot with device GraphDev */
    
    if (gselec == 1)
      fprintf(LogFilePtr,"Open visual output device - display\n");
    else if (gselec == 2)
      fprintf(LogFilePtr,"Open visual output device - file \n");
    else if (gselec == 3)
      fprintf(LogFilePtr,"Open visual output device - display+file \n");
    else {
      fprintf(LogFilePtr,"ERROR: Incorrect output device, correct option -o, value 1,2 or 3 \n");
      exit(-1);
    }

    fprintf(LogFilePtr,"Visual activation for the first  %ld  trajectories. \n", BufferSize);
    /*	    initialize pgplot with device GraphDev  */
    if (cpgopen(GDEV) < 1) {
      fprintf(LogFilePtr,"ERROR: I cannot open graphical output device. \n");    
      exit(-1);
    }
    fprintf(LogFilePtr,"Minimum neutron weigth for visualisation: %f. \n", wei_min1);

    if (p == 2) {
      cpgenv((float) Windw, (float) WindW, (float) Windh, (float) WindH, 0, 0);
      cpgsfs((int) 2);
      cpgsch((float) 1.2);
      cpglab("X, cm","Y, cm","NEUTRON VISUALISATION - PLANE X0Y");
    } else if (p == 3) {
      cpgenv((float) Windw, (float) WindW, (float) Windh, (float) WindH, 0, 0);	
      cpgsfs((int) 2);
      cpgsch((float) 1.2);
      cpglab("X, cm","Z, cm","NEUTRON VISUALISATION - PLANE X0Z");
    } else if (p == 4) {
      cpgenv((float) (float) Windw, (float) (float) WindW, (float) (float) Windh, (float) (float) WindH, 0, 0);	
      cpgsfs((int) 2);
      cpgsch((float) 1.2);
      cpglab("Y, cm","Z, cm","NEUTRON VISUALISATION - PLANE Y0Z");
    }
  }
#endif	

  /* init for ASCII output */

  NumWrong=0;

  IntegralIntensity = 0.; Path = 0.; 

  OutputAngleHoriz	*= M_PI/180.;
  OutputAngleVert	*= M_PI/180.;
  FillRotMatrixZY(RotMatrixOut, OutputAngleVert, OutputAngleHoriz);
	
  ReadParameterFile(); 

  if (Par_Field != NULL) fclose(Par_Field);

} /* End OwnInit */


void OwnCleanup()
{
#ifdef VT_GRAPH
  if ((p == 2)||(p == 3)||(p == 4)||(p == 100))  cpgclos();    
#endif
  if(NumWrong !=0) fprintf(LogFilePtr,"\nERROR:  %ld collided trajectories have wrong input spin. \n\n", NumWrong);
}/* End OwnCleanup */


/* ReadParameterFile() reads the parameters from a file */

void ReadParameterFile()
{
  // You may give more than MAX_MIRR description lines in the parameter file,
  // but only the first MAX_MIRR mirrors which are marked valid are taken into account.
  // If less than MAX_MIRR description lines are given, max_mirr will be less 
  // than MAX_MIRR.
  int l = 1;
  while (1) {
    int use_this_mirror;
    if (feof(Par_Field)) break;
    r1[l][0] = r2[l][0] = r3[l][0] = r4[l][0] = 0;
    use_this_mirror = ReadParI(Par_Field);
    r1[l][1]=ReadParF(Par_Field); r1[l][2]=ReadParF(Par_Field);
    r2[l][1]=ReadParF(Par_Field); r2[l][2]=ReadParF(Par_Field);
    r3[l][1]=ReadParF(Par_Field); r3[l][2]=ReadParF(Par_Field);
    r4[l][1]=ReadParF(Par_Field); r4[l][2]=ReadParF(Par_Field); 
    WallOffset[l][0]=ReadParF(Par_Field); WallOffset[l][1]=ReadParF(Par_Field); WallOffset[l][2]=ReadParF(Par_Field); 
    WallHoriz[l]=ReadParF(Par_Field); WallVert[l]=ReadParF(Par_Field); 
    mrangh[l]=ReadParF(Par_Field); mrangv[l]=ReadParF(Par_Field); 
    thetaC[l][0]=ReadParF(Par_Field); thetaCSM[l][0]=ReadParF(Par_Field); RthetaCSM[l][0]=ReadParF(Par_Field);
    mued[l][0]=ReadParF(Par_Field); mued[l][1]=ReadParF(Par_Field); 
    thetaC[l][1]=ReadParF(Par_Field); thetaCSM[l][1]=ReadParF(Par_Field); RthetaCSM[l][1]=ReadParF(Par_Field);
    mued[l][2]=ReadParF(Par_Field); mued[l][3]=ReadParF(Par_Field);  
    ReadParComment(Par_Field); 

    if(use_this_mirror) {
      fprintf(LogFilePtr,
	      "\n%c:  %10.5f%10.5f  %10.5f%10.5f  %10.5f%10.5f  %10.5f%10.5f  %10.5f%10.5f%10.5f  "
	      "%10.5f%10.5f  %10.5f%10.5f  %10.5f%10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f%10.5f \n", 
	      'A'+l-1,r1[l][1],r1[l][2],r2[l][1],r2[l][2],r3[l][1],r3[l][2],r4[l][1],r4[l][2],
	      WallOffset[l][0],WallOffset[l][1],WallOffset[l][2], WallHoriz[l],WallVert[l],
	      mrangh[l], mrangv[l], thetaC[l][0], thetaCSM[l][0], RthetaCSM[l][0],
	      mued[l][0], mued[l][1], thetaC[l][1], thetaCSM[l][1], RthetaCSM[l][1], mued[l][2], mued[l][3]);
      WallHoriz[l] *= M_PI/180.;
      WallVert[l]  *= M_PI/180.;
      mrangh[l]	 *= M_PI/180.;
      mrangv[l]	 *= M_PI/180.;
      FillRotMatrixZY(RotMatrixWall[l],  WallVert[l],  WallHoriz[l]);
      EulerToCartesianZY(WallNormal[l], &WallVert[l], &WallHoriz[l]);
      l++;
      // terminate if we have MAX_MIRR valid mirrors
      if (l > MAX_MIRR)
	break;
    }
    // else this mirror is to be skipped
  }
  max_mirr = l - 1;
  
  fprintf(LogFilePtr,"\n");   
  for (l=1; l<=max_mirr; l++)
    fprintf(LogFilePtr,"%c", 'A'+l-1);   

  fprintf(LogFilePtr,"\n");   
  
}/* End ReadParameterFile  */





/* Theta and Phi of a unit vector if Theta is the angle with the axis of lowest index  */

void CartesianToSpherical2(VectorType Vector, double *Theta, double *Phi)
{
	*Theta	= (double) acos(Vector[0]);

	if(Vector[2]>=0) *Phi	= (double) atan2(Vector[2], Vector[1]);
	if(Vector[2]< 0) *Phi	= 2. * M_PI + (double) atan2(Vector[2], Vector[1]);

}



/* hittriangle */

int hittriangle(VectorType r1, VectorType r2, VectorType rt)
{
	double the1, the2, thet, phi1, phi2, phit;

	CartesianToSpherical2(r1, &the1, &phi1);
	CartesianToSpherical2(r2, &the2, &phi2);
	CartesianToSpherical2(rt, &thet, &phit);

	if(phi1 < phi2) {
		if((phi1 < phit)&&(phit < phi2)&&(Area(r1,r2) > (Area(r1,rt) + Area(r2, rt)))) return 1; 
		else return 0;
	}
	if(phi1 > phi2) {
		if(((phi1 < phit)&&(phit < 2*M_PI)||(0. < phit)&&(phit < phi2))&&(Area(r1,r2) > (Area(r1,rt) + Area(r2, rt)))) return 1; 
		else return 0;
	}

	else return 0;
}

/* hitwall  calculate coordinates where will be hit the wall */

int hitwall(VectorType r1, VectorType r2, VectorType r3, VectorType r4, VectorType rt)
{
	if((hittriangle(r1, r2, rt) == 1)||(hittriangle(r2, r3, rt) == 1)||(hittriangle(r3, r4, rt) == 1)||(hittriangle(r4, r1, rt) == 1)) return 1;
	else return 0;
}

/* routine which computes the collision with a wall */

double CollideWall(double *prob, VectorType pos,VectorType dir,VectorType spin,
		   VectorType WallOffset,VectorType WallNormal, double  RotMatrixWall[3][3],
		   VectorType r1,VectorType r2,VectorType r3,VectorType r4,
		   double thetaC[2], double thetaCSM[2], double RthetaCSM[2], double mued[4], double mrangh, double mrangv)
{
	VectorType rt; double path;  double rangh=0., rangv=0., RotMatrixRang[3][3];	

	/* random angular spread */ 
	{
/*	rangh=MonteCarlo(0., mrangh);
	rangv=MonteCarlo(0., mrangv);
*/
/*      nonsense to call MonteCarlo if the mrang value is 0
	rangh=MonteCarlo(- mrangh/2., mrangh/2.);
	rangv=MonteCarlo(- mrangv/2., mrangv/2.);
*/
        rangh = mrangh == 0 ? 0 : MonteCarlo(- mrangh/2., mrangh/2.);
        rangv = mrangv == 0 ? 0 : MonteCarlo(- mrangv/2., mrangv/2.);
	FillRotMatrixZY(RotMatrixRang, rangv, rangh); /*fprintf(LogFilePtr,"%e  %e \n", rangh, rangv);*/
	}
	
	if (PlaneLineIntersect(pos, dir, WallNormal, (double) ScalarProduct(WallOffset, WallNormal), rt)== 0)
	  return 99999;

	{ VectorType replacement;

	CopyVector(rt, replacement);

	SubVector(replacement, pos);

	if((ScalarProduct(replacement, dir)< 0.)/*||(ScalarProduct(replacement, WallNormal)< 0.)*/) return 99999;	
	}


	/* transforms into frame of the wall  */
	
	SubVector(rt, WallOffset);

	RotVector(RotMatrixWall, rt);

	RotVector(RotMatrixWall, dir);

	if((mrangh!=0.)||(mrangv!=0.))
	{
	RotVector(RotMatrixRang, rt);

	RotVector(RotMatrixRang, dir);
	}



	/* here comes to collision etc. */
	{double the, Choise, Refl[2], expon[2];
		
		if(hitwall(r1, r2, r3, r4, rt)==0) 
		{		
		RotBackVector(RotMatrixWall, dir);

		return 99999;
		}
		
		Choise = MonteCarlo(0,1);
		the = M_PI_2 - (double) acos(fabs(dir[0])); 

		if(fabs(SpinVector[quant_dir])!=1. && useQuantDir) {NumWrong+=1; *prob = 0. ;}

		/* here if spin up */
		if(SpinVector[quant_dir]==1. || !useQuantDir)
		{
			if(the <= thetaC[0] * WL) 
			{	dir[0] *= -1.;
			}
			else
			{	
				expon[0] = (mued[0] * WL + mued[1]) /(double) sqrt(sq((double) sin(the)) - sq((double) sin(thetaC[0] * WL)));
 
//				expon[0] = mued[0] * WL + mued[1];
				
				if(expon[0] > 500) expon[0] = 500;
				if((the > thetaC[0] * WL)&&(the <= thetaCSM[0] * WL))
				{
					Refl[0] = RthetaCSM[0] + (1. - RthetaCSM[0])/(thetaCSM[0] * WL - thetaC[0] * WL)*(thetaCSM[0] * WL - the);
					if(Choise < Refl[0]) dir[0] *= -1.; 
					else *prob *= (double) exp(- expon[0]);
				}
				if(the >  thetaCSM[0] * WL) *prob *= (double) exp(- expon[0]);
			}
		}

		/* here if spin down */
		if(SpinVector[quant_dir]==-1.)
		{

			if(the <= thetaC[1] * WL) 
			{	dir[0] *= -1.;
			}
			else
			{
				expon[1] = (mued[2] * WL + mued[3]) /(double) sqrt(sq((double) sin(the)) - sq((double) sin(thetaC[1] * WL)));
 
//				expon[1] = mued[2] * WL + mued[3];
				
				if(expon[1] > 500) expon[1] = 500;
				if((the > thetaC[1] * WL)&&(the <= thetaCSM[1] * WL))
				{
					Refl[1] = RthetaCSM[1] + (1. - RthetaCSM[1])/(thetaCSM[1] * WL - thetaC[1] * WL)*(thetaCSM[1] * WL - the);
					if(Choise < Refl[1]) dir[0] *= -1.; 
					else *prob *= (double) exp(- expon[1]);
				}
				if(the >  thetaCSM[1] * WL) *prob *= (double) exp(- expon[1]);
			}
		}

	}
/* transforms back into original frame  */
	
	if((mrangh!=0.)||(mrangv!=0.))
	{
	RotBackVector(RotMatrixRang, rt);

	RotBackVector(RotMatrixRang, dir);
	}

	RotBackVector(RotMatrixWall, rt);

	AddVector(rt, WallOffset);

	RotBackVector(RotMatrixWall, dir);

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


/* computes pathlength until collision */

	path = DistVector(rt, pos); 
	
	CopyVector(rt, pos); 
	return path; 

}








