//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
//             LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION 
// 
//  ID    :LEMuSRElectricField.cc , v 1.3
//  AUTHOR: Taofiq PARAISO
//  DATE  : 2004-09-17 10:20
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
//
//  &           &&&&&&&&&&                            &&&&&&&         &&&&&&&&
//  &           &                                   &&      &&       &       &&
//  &           &                   &      &       &                &       &&
//  &           &&&&&&&            &      &         &&&&&&         &&&&&&&&  
//  &           &                 &      &&                &      &      &&
//  &           &                &&     &  &     &&      &&      &        &
//  &&&&&&&&&&  &&&&&&&&&&      &  &&&&&    &&    &&&&&&&       &        &&   
//                             &
//                            &
//                           &
//                         &  
//                               Electric Field   
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
#include"LEMuSRElectricField.hh"
#include"G4UnitsTable.hh"


LEMuSRElectricField::LEMuSRElectricField(G4double fieldval,G4String Xfile, 
					 G4String Yfile,
					 G4String Zfile,
					 G4String map_length_unit,
					 G4double Offset, G4double nbx, G4double nby, G4double nbz)
{
  FieldVal=fieldval;
  map_unit = map_length_unit;
  G4cout<<G4endl
	<< "                   ++++++ELECTRIC FIELD BEINg BUILT++++++"<<G4endl;
  G4cout<< "                            +++++IMPORTANT++++++"<<G4endl;
  G4cout<< "+ MAKE SURE UNITS IN FIELD MAP AND IN LEMuSRElectricField.cc ARE COHERENT +"<<G4endl<<G4endl;

  // open files for reading
  std::ifstream fx(Xfile);
  std::ifstream fy(Yfile);
  std::ifstream fz(Zfile);

  // Ignore first blank line
  char Xbuffer[256];
  fx.getline(Xbuffer,256);
  char Ybuffer[256];
  fy.getline(Ybuffer,256);
  char Zbuffer[256];
  fz.getline(Zbuffer,256);

  zOffset = Offset;
  nx=nbx;
  ny=nby;
  nz=nbz;


  int ix, iy, iz;  

  // Read in the data
  double xval,yval,zval,bx,by,bz;
  xval=0;
  yval=0;
  zval=0;
  bx=0;
  by=0;
  bz=0;
  ix=0;
  iy=0;
  iz=0;
  G4int firstline;
  firstline=0;

  // ATTENTION :: FIELD MAP IS ASSUMED TO BEGIN WITH XMIN YMIN ZMIN AND END WITH XMAX YMAX ZMAX
  // THIS IS A CASE FOR MAPS GENERATED BY FEMLAB

  for (ix=0; ix<nx; ix++) {
    for (iy=0; iy<ny; iy++) {
      for (iz=0; iz<nz; iz++) {

	fx >> xval >> yval >> zval >> bx;
	fy >> xval >> yval >> zval >> by;
	fz >> xval >> yval >> zval >> bz;
	      
	xField[ix][iy][iz] = bx*volt/meter;
	yField[ix][iy][iz] = by*volt/meter;
	zField[ix][iy][iz] = bz*volt/meter;
	if (firstline==0)
	  {
	    minx =xval;
	    miny =yval;
	    minz =zval;
		  
	    firstline=1;
	  }

      }
    }
  }


  fx.close(); fy.close(); fz.close();

  maxx = xval;
  maxy = yval;
  maxz = zval;
 
  G4cout<<G4endl
	<< "                (-: ++++++ELECTRIC FIELD INITIALIZED++++++ :-)";
  G4cout<<"\n from files "<< Xfile << " etc. with following statistics : \n"
	<<"\n last line: " <<  xval <<" "<< yval <<" "<< zval<<" "<<bx<<" "<<by<<" "<<bz<<" with field "<<fieldval<<"\n"
	<<"\n boundaries: " <<  minx <<" "<< maxx <<" "<< miny<<" "<<maxy<<" "<<minz<<" "<<maxz<<" (in "<<map_unit<<")\n"
	<<G4endl
	<<G4endl;
 

  //#ifdef DEBUG_INTERPOLATING_FIELD2
  G4String output_name;
  G4cout<<"Output Field Map file name?"<<G4endl;
  G4cin>>output_name;
  std::ofstream fout (output_name);     // open destination file
  fout.close ();	
 

  fout.open(output_name);
  if (!fout.is_open()) exit(8);
   
  for (ix=0; ix<nx; ix++) {
    for (iy=0; iy<ny; iy++) {
      for (iz=0; iz<nz; iz++) {
	fout<< xField[ix][iy][iz]/volt*meter <<" " << yField[ix][iy][iz]/volt*meter <<" " << zField[ix][iy][iz]/volt*meter <<G4endl;
      } 
    }
  }

  fout<< maxx<<" " << maxy<<" " <<maxz <<G4endl;
  fout<< minx<<" " << miny<<" " <<minz <<G4endl;

  fout.close();

  G4cout<<"Field Map "<<output_name<<" genarated ___ press enter"<<G4endl;

  //#endif


}
 



LEMuSRElectricField::LEMuSRElectricField(G4double fieldval,G4String file,
					 G4String map_length_unit,
					 G4double Offset, 
					 G4double nbx, 
					 G4double nby, 
					 G4double nbz)
{

  FieldVal= fieldval;
  map_unit = map_length_unit;
  G4cout<<G4endl
	<< "++++++LEMUSR ELECTRIC FIELD INITIALIZATION+++++\n";
  G4cout<< "                 +++++IMPORTANT++++++"<<G4endl;
  G4cout<< "+ MAKE SURE UNITS IN FIELD MAP AND IN LEMuSRElectricField.cc ARE COHERENT +"<<G4endl;
  G4cout<< "++++++ELECTRIC FIELD BEINg BUILT++++++"<<G4endl<<G4endl;
 
  // open files for reading
  std::ifstream fmap(file);

  // Ignore first blank line
  //  char buffer[256];
  //  fmap.getline(buffer,256);

  zOffset = Offset;
  nx=nbx;
  ny=nby;
  nz=nbz;

  int ix, iy, iz;

 
  // Read in the data
  double bx,by,bz;
  for (ix=0; ix<nx; ix++) {
    for (iy=0; iy<ny; iy++) {
      for (iz=0; iz<nz; iz++) {
        fmap >>bx >>by >>bz;

        xField[ix][iy][iz] = bx*volt/meter;
	yField[ix][iy][iz] = by*volt/meter;
	zField[ix][iy][iz] = bz*volt/meter;

      }
    }
  }
  
  fmap>> maxx >> maxy >> maxz ;
  fmap>> minx >> miny >> minz ;
  G4cout<< maxx<<" " << maxy<<" " <<maxz <<G4endl;
  G4cout<< minx<<" " << miny<<" " <<minz <<G4endl;
  
  fmap.close();
  
  G4cout<< "++++++ELECTRIC FIELD INITIALIZED++++++";
  G4cout<<"\n from file "<< file <<G4endl<<G4endl;

}


LEMuSRElectricField::~LEMuSRElectricField()
{
  delete  []  xField ;
  delete  []  yField ;
  delete  []  zField ;
}

 

void LEMuSRElectricField::GetFieldValue(const G4double point[4],
					G4double *Bfield ) const
{
  // Values of x y z must be scaled according to the field map unit! important
  // be carefull and use the debuging to verify it

  G4double length_ratio, field_factor;
  length_ratio=1;
  field_factor=0;
 
  if (map_unit=="m")
    {
      length_ratio=1000;
      field_factor=1;
    }
  
  else if (map_unit=="dm")
    {
      length_ratio=100;
      field_factor=10;
    }
  else if (map_unit=="cm")
    {
      length_ratio=10;
      field_factor=100;
    }
  else if (map_unit=="mm")
    {
      length_ratio=1;
      field_factor=1000;
    }
  else 
    {
      G4cout<<"ERROR!!!! UNRECOGNIZED LENGTH UNIT FOR FIELD MAP ==>> zöro";
    }
  

  //! Scaling the distances acording to the unit.    
  G4double x = (point[0]/mm)/length_ratio;
  G4double y = (point[1]/mm)/length_ratio;
  G4double z = ((point[2] - zOffset)/mm)/length_ratio;


#ifdef DEBUG_INTERPOLATING_FIELD

  G4cout<<"POSITION:: " <<point[0] <<" mm "<< G4BestUnit(point[1],"Length") <<" "<<G4BestUnit(point[2],"Length")<<G4endl;
    
  G4cout<<"POSITION in the field map:: " <<x <<map_unit<<"  "<< y <<map_unit<<"  "<< z <<map_unit<< "  "<<G4endl;
#endif
 
 
  //check that the point is within the defined region 
  if ( x>=minx && x<=maxx && y>=miny && y<=maxy && (z)>=minz && (z)<=maxz )
    { 
      
      // Position of given point within region, normalized to the range
      // [0,1]
      G4double xfraction = (x - minx) / (maxx-minx);
      G4double yfraction = (y - miny) / (maxy-miny); 
      G4double zfraction = (z - minz) / (maxz-minz);//!!!!! PROBLEM MAY BE HERE
       
      // Need addresses of these to pass to modf below.
      // modf uses its second argument as an OUTPUT argument.
      G4double xdindex, ydindex, zdindex;
      
      // Position of the point within the cuboid defined by the
      // nearest surrounding tabulated points
       
      G4double xlocal = ( modf(xfraction*(nx-1), &xdindex));
      G4double ylocal = ( modf(yfraction*(ny-1), &ydindex));
      G4double zlocal = ( modf(zfraction*(nz-1), &zdindex));

      // The indices of the nearest tabulated point whose coordinates
      // are all less than those of the given point
      /*      int xindex = static_cast<int>(xdindex);
	      int yindex = static_cast<int>(ydindex);
	      int zindex = static_cast<int>(zdindex);
      */
      G4int xindex = (G4int)(xdindex);
      G4int yindex = (G4int)(ydindex);
      G4int zindex = (G4int)(zdindex);
      
       
#ifdef DEBUG_INTERPOLATING_FIELD
      G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
      G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
#endif
 
      //********************* FIELD INTERPOLATION***************************//
      //*********** FROM GEANT4 EXTENDED EXAMPLE PURGING MAGNET ************//

      // Full 3-dimensional version

      // example: if map_unit= 'dm' multiply by 10 since field map generated by femlab is in volt/decimeter! in order to get the field in V/m.
      
      Bfield[0] =
	(xField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
	 xField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
	 xField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
	 xField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
	 xField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
	 xField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
	 xField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
	 xField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal);
      Bfield[1] =
	(yField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
	 yField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
	 yField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
	 yField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
	 yField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
	 yField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
	 yField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
	 yField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal);
      Bfield[2] =
	(zField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
	 zField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
	 zField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
	 zField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
	 zField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
	 zField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
	 zField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
	 zField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal);
      
    } else {
      //      G4cout<<"RANGE ERROR=>field 0:: RANGE is " << minx <<" "  << maxx <<" " << miny <<" " << maxy <<" " << minz <<" "<< maxz <<" ";
      Bfield[0] = 0.0;
      Bfield[1] = 0.0;
      Bfield[2] = 0.0;
    }



  Bfield[0] = Bfield[0]*FieldVal*field_factor;
  Bfield[1] = Bfield[1]*FieldVal*field_factor;
  Bfield[2] = Bfield[2]*FieldVal*field_factor;


#ifdef DEBUG_INTERPOLATING_FIELD
  G4cout<<"Field Value " << G4BestUnit(Bfield[0]*meter,"Electric potential") <<"/m " <<Bfield[1]/volt*meter <<" V/m "<< Bfield[2]/volt*meter <<" V/m " <<G4endl;
  G4cout<<"FieldVal " <<FieldVal <<" coef "  << field_factor<<" \n";
#endif

}