//  Geant4 simulation for MuSR
//  AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
//  DATE  : 2008-05
//

#include "lem4TabulatedElementField3D.hh"
#include "lem4Parameters.hh"

#include "G4UnitsTable.hh"
#include <string>

lem4TabulatedElementField3D::lem4TabulatedElementField3D(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume),
fldType(fieldType), ffieldValue(fieldValue)
{
  G4cout << "\n-----------------------------------------------------------";

  // The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla)
  // NOTE: Should be the same as EMfieldUnit defined in DetectorConstruction.cc!
  if      (fldType == 'E') {G4cout << ("\n      Electric field 3D"); fUnit = "kV/mm"; fieUnit= kilovolt/mm;}
  else if (fldType == 'B') {G4cout << ("\n      Magnetic field 3D"); fUnit = "T";     fieUnit = tesla;}

  G4cout << "\n-----------------------------------------------------------" << G4endl;
  G4cout << "\n ---> Reading the "<< fldType <<"-field map from " << filename << " ... " << G4endl;

  // Open the file for reading
  std::ifstream file( filename );
  
  // Read the first line of the file header (contains the metadata)
  char buffer[256];
  file.getline(buffer,256);
  
  //G4cout << "The first line of file "<< filename << " is:\n"<< buffer << G4endl;
  // Read the file header
  ///file >> nx >> ny >> nz >> lUnit >> fieldNormalisation; // OLD VERSION
  
  // Read the number of arguments and decide filetype - 3 or 6 columns
  int n_arg = sscanf (buffer,"%d %d %d %s %lf %lf %lf %lf %lf %lf %lf",
                      &nx, &ny, &nz, lUnit, &fieldNormalisation,
                      &minimumx, &maximumx, &minimumy, &maximumy, &minimumz, &maximumz);
  // Add manually the length unit and norm. factor in the field-map file if missing!

  //printf ("The following data were read in: \n");
  //printf ("nx = %d, ny = %d, nz = %d, Unit = %s, Norm = %f,\n minx = %f, maxx = %f, miny = %f, maxy = %f, minz = %f, maxz = %f\n",nx,ny,nz,lUnit,fieldNormalisation,minimumx, maximumx, minimumy, maximumy, minimumz, maximumz);
  ///printf ("nx = %d, ny = %d, nz = %d, Unit = %s, Norm = %f,\n minx = %f\n",nx,ny,nz,lUnit,fieldNormalisation,minimumx);
  //printf ("Number of args read: %d\n", n_arg);

  G4cout << "      The grid consists of [" << nx << " x " << ny << " x " << nz << "] x, y, z values" << G4endl;
 
  // Get the field map LENGTH unit and its value in G4 native units.
  lenUnit = G4UnitDefinition::GetValueOf(lUnit);
  G4cout << "      Field length unit = " << lUnit << ", in G4 native units = " << lenUnit << G4endl;
  G4cout << "      Field map normalisation factor = " << fieldNormalisation << G4endl;

  // Set up storage space for the 3D table: notice the order!
  xField.resize(nx);
  yField.resize(nx);
  zField.resize(nx);
  int ix, iy, iz;
  for (ix=0; ix<nx; ix++) {
    xField[ix].resize(ny);
    yField[ix].resize(ny);
    zField[ix].resize(ny);
    for (iy=0; iy<ny; iy++) {
      xField[ix][iy].resize(nz);
      yField[ix][iy].resize(nz);
      zField[ix][iy].resize(nz);
    }
  }

  
  // Ignore header information. All lines whose first character
  // is '%' are considered to be part of the header.
  do {
     file.ignore(256, '\n');
   } while (file.peek() == '%');
  
  /// NOTE: The old version works only ONCE, i.e. at a STOP signal, e.g. !
  //do {
  //     file.getline(buffer, 256, '!');
  //} while (buffer[0] != '!');

  
  // Read in the data: [x, y, z, EMx, EMy, EMz] or [EMx, EMy, EMz]
  double xval, yval, zval; // Coordinates
  double fx,   fy,   fz;   // Field values
  //std::ofstream fileout("three_columns.txt"); // Useful for rewriting files
  
  for (ix=0; ix<nx; ix++) { // notice the order!
    for (iy=0; iy<ny; iy++) {
      for (iz=0; iz<nz; iz++) {
        
	if (n_arg == 5) { // Read BOTH coordinates and field values
	  file >> xval >> yval >> zval >> fx >> fy >> fz;
	  //fileout << fx <<"  \t"<< fy << "  \t"<< fz << std::endl; // for rewriting files
	  if ( ix==0 && iy==0 && iz==0 ) {
            minimumx = xval * lenUnit;
            minimumy = yval * lenUnit;
            minimumz = zval * lenUnit;
          }
	}
	else if (n_arg == 11) {
	  file >> fx >> fy >> fz; // Read ONLY field values
	}

        xField[ix][iy][iz] = fx*fieldNormalisation;
        yField[ix][iy][iz] = fy*fieldNormalisation;
        zField[ix][iy][iz] = fz*fieldNormalisation;
      }
    }
  }
  //fileout.close(); // for rewriting files
  file.close();

  if (n_arg == 5) { 
    maximumx = xval * lenUnit;
    maximumy = yval * lenUnit;
    maximumz = zval * lenUnit;
  }
  else if (n_arg == 11) { 
    minimumx = minimumx * lenUnit;
    minimumy = minimumy * lenUnit;
    minimumz = minimumz * lenUnit;
    maximumx = maximumx * lenUnit;
    maximumy = maximumy * lenUnit;
    maximumz = maximumz * lenUnit;
  }

    
  G4String volumeName = logVolume->GetName().substr(4);
  
  G4cout << "      ... done reading " << G4endl;
  G4cout << "\n ---> " << fldType << "-field in volume "<< volumeName 
         << " set to: " << fieldValue/fieUnit << " " << fUnit << G4endl;
         if (n_arg == 5) { G4cout << "\n ---> Assumed order (6 col.):  x, y, z, ";}
   else if (n_arg == 11) { G4cout << "\n ---> Assumed order (3 col.):  ";}
  G4cout << fldType <<"x, "<< fldType <<"y, "<< fldType <<"z "
	 << "\n      Min values x, y, z: "
	 << minimumx/cm << " " << minimumy/cm << " " << minimumz/cm << " cm "
	 << "\n      Max values x, y, z: "
	 << maximumx/cm << " " << maximumy/cm << " " << maximumz/cm << " cm" << G4endl;

  // Should really check that the limits are not the wrong way around.
  if (maximumx < minimumx) {Invert("x");}
  if (maximumy < minimumy) {Invert("y");}
  if (maximumz < minimumz) {Invert("z");}

  if ((maximumx < minimumx) || (maximumy < minimumy)|| (maximumz < minimumz)) {
    G4cout << "\n ---> After reordering:"
	 << "\n      Min values x, y, z: "
	 << minimumx/cm << " " << minimumy/cm << " " << minimumz/cm << " cm "
	 << "\n      Max values x, y, z: "
	 << maximumx/cm << " " << maximumy/cm << " " << maximumz/cm << " cm " << G4endl;
  }

  dx = maximumx - minimumx;
  dy = maximumy - minimumy;
  dz = maximumz - minimumz;

  G4cout << "      Range of values: "
  	 << dx/cm << " cm (in x), " << dy/cm << " cm (in y) and " << dz/cm << " cm (in z)."
  	 << "\n-----------------------------------------------------------\n" << G4endl;

}

void lem4TabulatedElementField3D::addFieldValue(const G4double point[4],
				      G4double *field ) const
{
  G4double EMf[3];  // EM field value obtained from the field map

  G4ThreeVector global(point[0],point[1],point[2]);
  G4ThreeVector local;

  local = global2local.TransformPoint(global);

  double x = local.x();
  double y = local.y();
  double z = local.z();

  //  G4cout<<"Global points= "<<point[0]<<", "<<point[1]<<", "<<point[2]<<",  Local point= "<<x<<", "<<y<<", "<<z<<G4endl;


  // Check that the point is within the defined region
  if ( x>minimumx && x<maximumx &&
       y>minimumy && y<maximumy &&
       z>minimumz && z<maximumz ) {

    // Position of given point within region, normalized to the range
    // [0,1]
    double xfraction = (x - minimumx) / dx;
    double yfraction = (y - minimumy) / dy;
    double zfraction = (z - minimumz) / dz;

    // Need addresses of these to pass to modf below.
    // modf uses its second argument as an OUTPUT argument.
    double xdindex, ydindex, zdindex;

    // Position of the point within the cuboid defined by the
    // nearest surrounding tabulated points
    double xlocal = ( modf(xfraction*(nx-1), &xdindex));
    double ylocal = ( modf(yfraction*(ny-1), &ydindex));
    double 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);

    //cks  The following check is necessary - even though xindex and zindex should never be out of range,
    //     it may happen (due to some rounding error ?).  It is better to leave the check here.
    if ((xindex<0)||(xindex>(nx-2))) {
      std::cout<<"SERIOUS PROBLEM:  xindex out of range!  xindex="<<xindex<<"   x="<<x<<"  xfraction="<<xfraction<<std::endl;
      if (xindex<0) xindex=0;
      else xindex=nx-2;
    }
    if ((yindex<0)||(yindex>(ny-2))) {
      std::cout<<"SERIOUS PROBLEM:  yindex out of range!  yindex="<<yindex<<"   y="<<y<<"  yfraction="<<yfraction<<std::endl;
      if (yindex<0) yindex=0;
      else yindex=ny-2;
    }
    if ((zindex<0)||(zindex>(nz-2))) {
      std::cout<<"SERIOUS PROBLEM:  zindex out of range!  zindex="<<zindex<<"   z="<<z<<"  zfraction="<<zfraction<<std::endl;
      if (zindex<0) zindex=0;
      else zindex=nz-2;
    }

        // Full 3-dimensional version
    EMf[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 ;
    EMf[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 ;
    EMf[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 ;

    EMf[0] *= ffieldValue;
    EMf[1] *= ffieldValue;
    EMf[2] *= ffieldValue;

    G4ThreeVector finalField(EMf[0],EMf[1],EMf[2]);
    finalField = global2local.Inverse().TransformAxis(finalField);

  if (fldType == 'B') {
    field[0] += finalField.x();
    field[1] += finalField.y();
    field[2] += finalField.z();
    }
  else if (fldType == 'E') {
    field[3] += finalField.x();
    field[4] += finalField.y();
    field[5] += finalField.z();
    }

  }
  //  G4cout<<"Kamil: Field: ("<<field[0]/tesla<<","<<field[1]/tesla<<","<<field[2]/tesla<<")"<<G4endl;

}



G4double lem4TabulatedElementField3D::GetNominalFieldValue() {
  return ffieldValue;
}

void lem4TabulatedElementField3D::SetNominalFieldValue(G4double newFieldValue) {
  //  // Rescale the magnetic field for a new value of the magnetic field
  ffieldValue=newFieldValue;
  G4cout<<"lem4TabulatedElementField3D.cc:  ffieldValue changed to="<< ffieldValue/fieUnit << fUnit << G4endl;
}

void lem4TabulatedElementField3D::Invert(const char* indexToInvert) {
  // This function inverts the indices of the field table for a given axis (x or z).
  // It should be called in the case when the x or z coordinate in the initial
  // field table is ordered in the decreasing order.
  std::vector< std::vector< std::vector< double > > > xFieldTemp(xField);
  std::vector< std::vector< std::vector< double > > > yFieldTemp(yField);
  std::vector< std::vector< std::vector< double > > > zFieldTemp(zField);
  G4bool invertX=false;
  G4bool invertY=false;
  G4bool invertZ=false;

  G4cout<<"Check that the lem4TabulatedElementField3D::Invert() function works properly!"<<G4endl;
  G4cout<<"It has not been tested yet!"<<G4endl;

  if (strcmp(indexToInvert,"x")==0) {invertX=true; std::swap(maximumx,minimumx);}
  if (strcmp(indexToInvert,"y")==0) {invertY=true; std::swap(maximumx,minimumx);}
  if (strcmp(indexToInvert,"z")==0) {invertZ=true; std::swap(maximumz,minimumz);}

  for (int ix=0; ix<nx; ix++) {
    for (int iy=0; iy<ny; iy++) {
      for (int iz=0; iz<nz; iz++) {
	if (invertX) {
	  xField[ix][iy][iz] = xFieldTemp[nx-1-ix][iy][iz];
	  yField[ix][iy][iz] = yFieldTemp[nx-1-ix][iy][iz];
	  zField[ix][iy][iz] = zFieldTemp[nx-1-ix][iy][iz];
	}
	else if(invertY) {
	  xField[ix][iy][iz] = xFieldTemp[ix][ny-1-iy][iz];
	  yField[ix][iy][iz] = yFieldTemp[ix][ny-1-iy][iz];
	  zField[ix][iy][iz] = zFieldTemp[ix][ny-1-iy][iz];
	}
	else if(invertZ) {
	  xField[ix][iy][iz] = xFieldTemp[ix][iy][nz-1-iz];
	  yField[ix][iy][iz] = yFieldTemp[ix][iy][nz-1-iz];
	  zField[ix][iy][iz] = zFieldTemp[ix][iy][nz-1-iz];
	}
      }
    }
  }


}