#ifndef LEM3AXIAL2DELFIELD_H
#define LEM3AXIAL2DELFIELD_H 1

#include "F04ElementField.hh"
#include "F04GlobalField.hh" /// ATTENTION: USE GENERAL EM GLOBAL AND NOT MAGN. FIELD ONLY

///#include "G4MagneticField.hh"
///#include "G4ElectricField.hh"
#include "globals.hh"
#include "G4ios.hh"

#include <fstream>
#include <vector>
#include <cmath>

class lem4Axial2DElField : public F04ElementField ///G4ElectricField
{

//! Contructor from a single field map
//  lem4ElectricField(G4double fieldval,G4String file,G4String map_length_unit,
		      //G4double Offset, G4double nx, G4double ny, G4double nz);
  ///lem4Axial2DElField(const G4String filename, double fieldValue, double lenUnit, double fieldNormalisation, double offset);

public: 
  // Class constructor for 2D axial field map (x, z, Er, Ez)
    lem4Axial2DElField(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter);
  //  "lenUnit" is the unit in which the grid coordinates are specified in the table
  //  "fieldNormalisation" is the normalisation that has to be applied on the field values in the table
  //                       such that when applying V_L3 = 1 kV the E values coincide with those in the table
  //  "fieldValue" is the field value (in kV) set by the user (i.e. values normalised to 1 kV will be
  //                       multiplied by this value).

  // Virtual destructor
  virtual ~lem4Axial2DElField() {};
  
  // addFieldValue() adds the field for THIS particular map into field[].
  // point[] is expressed in GLOBAL coordinates.
  void  addFieldValue(const G4double point[4], G4double* field) const;
  
  // Usual Set and Get functions
  G4double GetNominalFieldValue();
  void SetNominalFieldValue(G4double newFieldValue);
  
  //  getWidth(), getHeight(), getLength(),  return the dimensions of the field 
  // (used to define the boundary of the field)
  virtual G4double getWidth()  { return 2*dr; }   // x coordinate
  virtual G4double getHeight() { return 2*dr; }   // y coordinate
  virtual G4double getLength() { return 2*dz; }   // z coordinate
  
  ///void GetFieldValue(const double Point[4], double *Efield) const;
  ///void SetFieldValue(double newFieldValue);
  /// G4double GetFieldSetValue();
  

private:
  // Storage space for the table
  std::vector < std::vector < double > > rField;
  std::vector < std::vector < double > > zField;
  // The dimensions of the table
  int nr, nz;
  // The physical limits of the defined region and the global offset
  double minr, maxr, minz, maxz;
  // The physical extent of the defined region
  double dr, dz;
  double ffieldValue; ///, zOffset; included in log_vol offset!
  ///bool   invertR, invertZ; // substituted by the Invert function
  // Utility function for inverting field map
  void   Invert(const char* indexToInvert); 

};

#endif