#include "lem4Axial2DElField.hh"
#include "G4UnitsTable.hh"
//#include "lem4Parameters.hh"

lem4Axial2DElField::lem4Axial2DElField(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume),
 ffieldValue(fieldValue)

//lem4Axial2DElField::lem4Axial2DElField(const G4String filename, double fieldValue, double lenUnit, double fieldNormalisation, double offset) : ffieldValue(fieldValue/1000.0), zOffset(offset), invertR(false), invertZ(false)
{
  //  double lenUnit   = decimeter; = 100   // Base unit mm
  //  double fieldUnit = kV;        = 0.001 // Base unit MV
  //  double fieldNorm = 0.00001 = 1E-5     // = 0.001/100
  G4cout << "\n-----------------------------------------------------------"
	 << "\n      Electric field"
	 << "\n-----------------------------------------------------------"
         << G4endl;
  // G4cout << " kilovolt = "<< kilovolt << G4endl; // = 1e-3 MV - default G4 unit
  ///G4cout << "\n Potential set to "<< fieldValue/1000.0/kilovolt << " kV"<< G4endl;
  G4cout << "\n Potential set to "<< fieldValue/kilovolt << " kV"<< G4endl;
  G4cout << "\n ---> " "Reading 2D E-field map from " << filename << " ... " << G4endl;
  std::ifstream file(filename); // Open the file for reading

  // Start reading file content
  char buffer[256];
  
  // Read table dimensions first
  file >> nr >> nz; // r - radial, z - longitudinal coordinate 

  G4cout << "      2D field table dimensions (r-z): ["
         << nr << ", " << nz << "] "
         << G4endl;

  // Set up storage space for table
  rField.resize(nr);
  zField.resize(nr);
  int ir, iz;

  for (ir = 0; ir < nr; ir++) {
    rField[ir].resize(nz);
    zField[ir].resize(nz);
  }

  // Ignore header information. All lines whose first character
  // is '%' are considered to be part of the header.
  do {
    file.getline(buffer, 256);
  } while (buffer[0] != '%');

  // Read in the data: [r, z, Er, Ez]
  double rval, zval, er, ez;
  for (ir = 0; ir < nr; ir++) {
    for (iz = 0; iz < nz; iz++) {
      file >> rval >> zval >> er >> ez;
      if (ir == 0 && iz == 0) {
	minr = rval * lenUnit;
	minz = zval * lenUnit;
      }
      rField[ir][iz] = er*fieldNormalisation;
      zField[ir][iz] = ez*fieldNormalisation;
    }
  }
  G4cout << "      Normalisation coeff.: " << fieldNormalisation << G4endl;
  file.close();

  maxr = rval * lenUnit;
  maxz = zval * lenUnit;

  G4cout << "\n ---> ... done reading (assumed order: r, z, Er, Ez)." << G4endl;
  G4cout << "\n ---> Min values r, z: "
	 << minr/cm << " " << minz/cm << " cm "
	 << "\n ---> Max values r, z: "
	 << maxr/cm << " " << maxz/cm << " cm " << G4endl;

  // Should really check that the limits are not the wrong way around.
  if (maxr < minr) {Invert("x");} ///{std::swap(maxr, minr); invertR = true;}
  if (maxz < minz) {Invert("z");} ///{std::swap(maxz, minz); invertZ = true;}
  if (maxr < minr || maxz < minz){ ///(invertR == true || invertZ == true){
  G4cout << "\n      After reordering:"
	 << "\n ---> Min values r, z: "
	 << minr/cm << " " << minz/cm << " cm "
	 << "\n ---> Max values r, z: "
	 << maxr/cm << " " << maxz/cm << " cm " << G4endl;
  }

  dr = maxr - minr;
  dz = maxz - minz;
  
  G4cout << " ---> Range of values: "
	 << dr/cm << " cm (in r) and " << dz/cm << " cm (in z)."
	 << "\n-----------------------------------------------------------\n" << G4endl;
}


///lem4Axial2DElField::~lem4Axial2DElField() /// Made virtual!
///{
///  delete [] &rField ;
///  delete [] &zField ;
///}


void lem4Axial2DElField::addFieldValue(const G4double  point[4], 
                                             G4double *field) const
{
  G4double B[3]; // Field value obtained from the field table
  ///TS
  ///G4cout<<"TTfileTT Global coord. ("<<point[0]/mm<<", "<<point[1]/mm<<", "<<point[2]/mm<<")."<<G4endl;
  G4ThreeVector global(point[0],point[1],point[2]);
  G4ThreeVector local;

  local = global2local.TransformPoint(global);
  ///TS
  ///G4cout<<"TTfileTT Local coord. ("<<local[0]/mm<<", "<<local[1]/mm<<", "<<local[2]/mm<<")."<<G4endl;
  ///G4cout<<"The LOGICAL volume is named "<<lvolume->GetName()<<G4endl;
    
  double r = sqrt(local.x()*local.x()+local.y()*local.y());
  double z = fabs(local.z());
  double z_sign = (local.z()>0) ? 1.:-1.;
  
  // Check that the point is within the defined region 
  if ( r < maxr && z < maxz ) {
    // Relative point position within region, normalized to [0,1].
    double rfraction = (r - minr) / dr;
    double zfraction = (z - minz) / dz;
    
    // Need addresses of these to pass to modf below.
    // modf uses its second argument as an OUTPUT argument.
    double rdindex, zdindex;
    
    // Position of the point within the cuboid defined by the
    // nearest surrounding tabulated points
    double rlocal = ( modf(rfraction*(nr-1), &rdindex));
    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 rindex = static_cast<int>(rdindex);
    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 ((rindex<0)||(rindex>(nr-2))) {
      std::cout<<"SERIOUS PROBLEM:  rindex out of range!  rindex="<<rindex<<"   r="<<r<<"  rfraction="<<rfraction<<std::endl;
      if (rindex<0) rindex=0;
      else rindex=nr-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;
    }

    //    G4cout<<"xField["<<xindex<<"]["<<zindex<<"]="<<xField[xindex  ][zindex  ]<<G4endl;
    //    G4cout<<"zField["<<xindex<<"]["<<zindex<<"]="<<zField[xindex  ][zindex  ]<<G4endl;
   
    // Interpolate between the neighbouring points
    double Bfield_R =
      rField[rindex  ][zindex  ] * (1-rlocal) * (1-zlocal) +
      rField[rindex  ][zindex+1] * (1-rlocal) *    zlocal  +
      rField[rindex+1][zindex  ] *    rlocal  * (1-zlocal) +
      rField[rindex+1][zindex+1] *    rlocal  *    zlocal  ;
    B[0] = (r>0) ? Bfield_R * (local.x() /r) : 0.;
    B[1] = (r>0) ? Bfield_R * (local.y() /r) : 0.;
    B[2] =
      zField[rindex  ][zindex  ] * (1-rlocal) * (1-zlocal) +
      zField[rindex  ][zindex+1] * (1-rlocal) *    zlocal  +
      zField[rindex+1][zindex  ] *    rlocal  * (1-zlocal) +
      zField[rindex+1][zindex+1] *    rlocal  *    zlocal  ;
    
    B[0] *= ffieldValue * z_sign;
    B[1] *= ffieldValue * z_sign;
    B[2] *= ffieldValue;

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

    field[0] += finalField.x();
    field[1] += finalField.y();
    field[2] += finalField.z();
  }
  //  G4cout<<"Kamil: Field: ("<<field[0]/tesla<<","<<field[1]/tesla<<","<<field[2]/tesla<<")"<<G4endl;

}



/* OLD IMPLEMENTATION - CHANGED WITH GLOBALFIELD! TS
void lem4Axial2DElField::GetFieldValue(const double point[4], double *Efield) const
{
  //  musrEventAction* myEventAction= musrEventAction::GetInstance();
  //  G4int evNr=myEventAction->GetLatestEventNr();
  //  G4int evNrKriz=6795;  //457;  //14250
  //  if (evNr==evNrKriz) {
  //    std::cout<<"evNr="<<evNr<<std::endl;
  //    printf ("for point= %f %f %f   B= %10.10f %10.10f %10.10f \n",
  //    	    point[0],point[1],point[2],Bfield[0]/tesla,Bfield[1]/tesla,Bfield[2]/tesla);
  //  }
  Efield[0]=0.; Efield[1]=0.; Efield[2]=0.;

  double r = sqrt(point[0]*point[0]+point[1]*point[1]);
  //Apply a fixed offset (in mm) to the z - coordinate.
  double rel_z = (point[2] - zOffset)/mm;
  double z = fabs(rel_z);
  double z_sign = (rel_z>0) ? 1.:-1.;
  //if (evNr==evNrKriz) std::cout<<"r ="<<r<<"   maxr="<<maxr<<std::endl;

#ifdef DEBUG_INTERPOLATING_FIELD
double dst = sqrt(point[0]*point[0] + point[1]*point[1]);
G4cout<<"POSITION:              "<< G4BestUnit(point[0], "Length")<<" "<< G4BestUnit(point[1], "Length") <<" "<< G4BestUnit(dst, "Length") <<" "<< G4BestUnit(point[2], "Length") <<G4endl;
  G4cout<<"POSITION in field map: "<< G4BestUnit(r, "Length") <<" "<< G4BestUnit(z, "Length")<<G4endl;
#endif


  // Check that the point is within the defined region
  if ( r<maxr && z<maxz ) {
    // if (evNr>evNrKriz) std::cout<<"bol som tu"<<std::endl;

    // Position of given point within region, normalized to the range
    // [0,1]
    double rfraction = (r - minr) / dr;
    double zfraction = (z - minz) / dz;

    if (invertR) { rfraction = 1 - rfraction;}
    if (invertZ) { zfraction = 1 - zfraction;}
    // Need addresses of these to pass to modf below.
    // modf uses its second argument as an OUTPUT argument.
    double rdindex, zdindex;

    // Position of the point within the cuboid defined by the
    // nearest surrounding tabulated points
    double rlocal = ( modf(rfraction*(nr-1), &rdindex));
    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 rindex = static_cast<int>(rdindex);
    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 ((rindex<0)||(rindex>(nr-2))) {
      std::cout<<"SERIOUS PROBLEM:  rindex out of range!  rindex="<<rindex<<"   r="<<r<<"  rfraction="<<rfraction<<std::endl;
      if (rindex<0) rindex=0;
      else rindex=nr-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;
    }

#ifdef DEBUG_INTERPOLATING_FIELD
      G4cout << "Local r, z: " << rlocal << " " << zlocal << G4endl;
      G4cout << "Index r, z: " << rindex << " " << zindex << G4endl;
#endif

    // Interpolate between the neighbouring points
    double Efield_R =
      rField[rindex  ][zindex  ] * (1-rlocal) * (1-zlocal) +
      rField[rindex  ][zindex+1] * (1-rlocal) *    zlocal  +
      rField[rindex+1][zindex  ] *    rlocal  * (1-zlocal) +
      rField[rindex+1][zindex+1] *    rlocal  *    zlocal  ;
    Efield[0] = (r>0) ? Efield_R * (point[0]/r) : 0.;
    Efield[1] = (r>0) ? Efield_R * (point[1]/r) : 0.;
    Efield[2] =
      zField[rindex  ][zindex  ] * (1-rlocal) * (1-zlocal) +
      zField[rindex  ][zindex+1] * (1-rlocal) *    zlocal  +
      zField[rindex+1][zindex  ] *    rlocal  * (1-zlocal) +
      zField[rindex+1][zindex+1] *    rlocal  *    zlocal  ;
// ATTENTION if dealing with MAGNETIC FIELDS reverse x, y and z signs!!
    Efield[0] = Efield[0] * ffieldValue; // * z_sign;
    Efield[1] = Efield[1] * ffieldValue; // * z_sign;
    Efield[2] = Efield[2] * ffieldValue * z_sign;
  }



  // Set some small field if field is almost zero (to avoid internal problems of Geant).
  if (sqrt(Efield[0]*Efield[0]+Efield[1]*Efield[1]+Efield[2]*Efield[2])<1.0E-12*volt) {
    //    if (evNr>evNrKriz) std::cout<<"bol som tuna"<<std::endl;
    //    Bfield[0] = 0.0;
    //    Bfield[1] = 0.0;
    //    Bfield[2] = Bfield[2] + 0.00001*tesla;
    #ifdef DEBUG_INTERPOLATING_FIELD
    G4cout<<"Zero field case"<<G4endl;
    #endif
    Efield[2] = 1.0E-12*volt;
  }
  //  Print the field (for debugging)
  //  if ((point[2]>-1*cm)&&(point[2]<1*cm)) {
  //    printf ("for point= %f %f %f   B= %10.10f %10.10f %10.10f \n",
  //  	    point[0],point[1],point[2],Bfield[0]/tesla,Bfield[1]/tesla,Bfield[2]/tesla);
  //  }
  //  if (evNr>evNrKriz) std::cout<<"this is the end"<<std::endl;

#ifdef DEBUG_INTERPOLATING_FIELD
//G4cout << " Volt/meter factor = "<< volt/meter << G4endl; // = 1e-9
G4cout<<"Field Value " << G4BestUnit(Efield[0]*meter,"Electric potential") <<"/m " <<Efield[1]/volt*meter <<" V/m "<< Efield[2]/volt*meter <<" V/m " <<G4endl;
G4cout<<"FieldVal " << ffieldValue/kilovolt << G4endl;
#endif

}

*/

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

void lem4Axial2DElField::SetNominalFieldValue(G4double newFieldValue) {
  // Rescale the EM field to account for the new applied field value
  ffieldValue=newFieldValue;
  G4cout<<"lem4Axial2DElField.cc:   ffieldValue changed to = "<< ffieldValue/kilovolt<<" kV."<<G4endl;
}

void lem4Axial2DElField::Invert(const char* indexToInvert) {
  // This function inverts the field table indices along a given axis (r or z).
  // It should be called whenever x or z coordinates in the initial field table
  // are ordered in a decreasing order.
  std::vector< std::vector< double > > rFieldTemp(rField);
  std::vector< std::vector< double > > zFieldTemp(zField);
  G4bool invertR=false;
  G4bool invertZ=false;

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

  if (strcmp(indexToInvert,"r")==0) {invertR=true; std::swap(maxr, minr);}
  if (strcmp(indexToInvert,"z")==0) {invertZ=true; std::swap(maxz, minz);}

  for (int ir=0; ir<nr; ir++) {
    for (int iz=0; iz<nz; iz++) {
      if (invertR) {
	rField[ir][iz] = rFieldTemp[nr-1-ir][iz];
      }
      else if(invertZ) {
	rField[ir][iz] = rFieldTemp[ir][nz-1-iz];
      }
    }
  }
  
}