#if !defined (SMP_H)
 #define SMP_H
 #include <fstream>
 #include "dataset.h"
 #include "vect.h"
 #include "box.h"
 #include <string>
 #include "cpgplot.h"
 using namespace std;
 
 
 struct pixel {
         double x,y,z;
 };
 
 
 class SMPMap
 {
  private:
  protected:
   int verbosity;
   int nx,ny;
   double maxx,maxy,minx,miny;
   double* mapdata;
   double maxvalue,minvalue;
   int min_max_flag;
         string id;
         bool IsCont(const int& x,const int& y,const double& val) const;
  public:
   static const int ok=0;
   static const int out_of_bounds=1;
   static const int could_not_open=2;
   static const int write_error=4;
   static const int could_not_close=8;
   static const int read_error=16;
   static const int could_not_create=32;
   SMPMap();
   SMPMap(const SMPMap &map);
   SMPMap(const int& xin,const int& yin,const double& sxin,const double& syin,const vect& origin);
   SMPMap(const vect& midin,const int& xin,const int& yin,const double& sxin,const double& syin);
   SMPMap(const SMPMap* mapin);
   SMPMap(const double& fsize,const double& fres,const SMPMap* mapin);
   SMPMap (const int& xin,const int& yin,const double& xmaxin,const double& ymaxin);
   SMPMap (const int& xin,const int& yin,const double& xminin,const double& xmaxin,const double& yminin,const double& ymaxin);
   SMPMap (const char* filename);
   SMPMap (const string& filename);
         SMPMap (const int& nxi,const int& nyi,const vector<pixel>& p);
   SMPMap& operator=(const SMPMap &map);
   double& operator[](const int &i) {return mapdata[i];}
   void SetID(const string& ins) {id=ins;}
   string& GetID() {return id;}
   double* pmapdata() {return mapdata;}
   friend SMPMap operator+(const SMPMap &map1,const SMPMap &map2);
   SMPMap operator+(const SMPMap &map2);
   friend SMPMap operator-(const SMPMap &map1,const SMPMap &map2);
   SMPMap operator-(const SMPMap &map2);
   friend SMPMap operator+(const SMPMap &map1,const double &fact);
   SMPMap operator+(const double &fact);
   friend SMPMap operator-(const SMPMap &map1,const double &fact);
   SMPMap operator-(const double &fact);
   friend SMPMap operator+(const double &fact,const SMPMap &map1);
   friend SMPMap operator-(const double &fact,const SMPMap &map1);
   friend SMPMap operator*(const SMPMap &map1,const SMPMap &map2);
   SMPMap operator*(const SMPMap &map2);
   friend SMPMap operator/(const SMPMap &map1,const SMPMap &map2);
   SMPMap operator/(const SMPMap &map2);
   friend SMPMap operator*(const SMPMap &map1,const double &fact);
   SMPMap operator*(const double &fact);
   friend SMPMap operator/(const SMPMap &map1,const double &fact);
   SMPMap operator/(const double &fact);
   friend SMPMap operator*(const double &fact,const SMPMap &map1);
   friend SMPMap operator/(const double &fact,const SMPMap &map1);
   friend ostream& operator<<(ostream& out,const SMPMap& m);
   friend istream& operator>>(istream& in,SMPMap& m);
   friend SMPMap Square(const SMPMap &map1);
   friend SMPMap Sqrt(const SMPMap &map1);
   friend SMPMap Log(const SMPMap &map1);
   friend SMPMap Abs(const SMPMap &map1);
   friend double Difference(const SMPMap &map1,const SMPMap &map2);
   friend double DifferenceSq(const SMPMap &map1,const SMPMap &map2);
   double GetSizeX() const {return (maxx-minx);}
   double GetSizeY() const {return (maxy-miny);}
   double size() const {return nx*ny;}
   int GetXPixelSize() const {return nx;}
   int GetYPixelSize() const {return ny;}
   int GetYYPixelSize() const {return ny;}
   void SetOrigin(const vect& origin);
   vect GetOrigin() const {return (vect(minx,miny));}
   void SetRange(const double& x1,const double& x2,const double& y1,const double& y2) {minx=x1;maxx=x2;miny=y1;maxy=y2;}
   vect GetMiddle() const {return (vect(minx+(maxx-minx)/2.0,miny+(maxy-miny)/2.0));}
   vect GetPixelSize() const;
   vect FindCentre() const;
   vect Pixel2vect(const int& x,const int& y) const;
   vect PixelMidpoint(const int& x,const int& y) const;
   vect PixelMidpoint(const int& i) const {return PixelMidpoint(i % nx,int(i/nx));}
   int Vect2pixel_x(const vect& position) const;
   int Vect2pixel_y(const vect& position) const;
   bool Inside(vect p) const {
     int x=Vect2pixel_x(p),y=Vect2pixel_y(p);
     if (x>=0 && y>=0 && x<nx && y<ny) return true;
     return false;
   }
         vector<pixel> getPixellist(const int& n);
         
   double Area(const double& val) const;
   bool SameRes(const SMPMap& map2) const;
   void SetMiddle(const vect& middle){SetOrigin(middle-vect((maxx-minx)/2.0,(maxy-miny)/2.0));}
   void SetCentre(const vect& centre);
   void FindMinMax();
   double GetMax();
   double GetMax(vect& vmax);
   double GetMin();
   double GetMin(vect& vmax);
   double GetXMax() const {return maxx;}
   double GetXMin() const {return minx;}
   double GetYMax() const {return maxy;}
   double GetYMin() const {return miny;}
   box getBounds() {return box(GetMiddle(),maxx-minx,maxy-miny);}
         double GetPixel(const int& x,const int& y) const {if (x>=0 && x<nx && y>=0 && y<ny) return mapdata[y*nx+x]; else return 0;}
   double GetPixel(const int& i) const {if (i>=0 && i<nx*ny) return mapdata[i]; else return 0;}
   double SetPixel(const int& x,const int& y,const double& val) {if (x>=0 && x<nx && y>=0 && y<ny) mapdata[y*nx+x]=val;}
   double SetPixel(const int& i,const double& val) {if (i>=0 && i<nx*ny) {mapdata[i]=val;return ok;} else return out_of_bounds;}
   int SetPixelsBox(vect bl,vect tr,double val);
   int SetPixelsCircle(const double& x,const double& y,const double& r,const double& val);
         double GetPoint(const vect& position) const;
   int addPixels(const vector<pixel>& pl);
   int setPixels(const vector<pixel>& pl);
   void TextDump() const;
   Dataset CutX(const double& val) const;
   Dataset CutY(const double& val) const;
   Dataset LineCut(const vect& start,const vect& end,const int& np) const;
   Dataset SumX() const;
   Dataset SumY() const;
   double FractVal(const double& minin,const double& maxin) const;
   SMPMap Threshold(const double& val) const;
   Dataset FollowContour(const double& val) const;
   void AddGaussianError(const double& stdev);
   void AddPoissonError(const double& stdev);
   void AddGaussianNoise(const double& stdev);
   int SaveMap (const char* name) const;
   int SaveBinary (const char* name) const;
   int LoadMap (const char* name);
   int LoadBinary (const char* name);
   int SaveFITS (const char* name) const;
   int SaveFITS (const char* name,const int& n) const;
   int SaveFITS (const char* name,const int& n,const char* extension) const;
   int SaveMatLab (const char* name) const;
   int LoadFITS (const char* name,const int& n);
   int LoadFITS (const char* name,const int& n,const char* extension);
   int LoadFITS (const char* name);
         int ReadSubimage(const char* name,int x,int y,int nx,int ny);
         int ReadThumb(const char* name,int sx,int sy);
   void Fill(const double& fillwith);
   void Fill(const double& fillwhat,const double& fillwith);
   void Saturate(const double& thresh);
   void CutLow(const double& thresh,const double& val);
         void PositiveDef();
   void FixRange(const double& low,const double& thresh);
   int FillGaps(const double& fillwhat,int limit=10);
   int SetPoint(const vect& pos,const double& value);
   int AddToPoint(const vect& pos,const double& value);
   void ShiftMap(const int& sx,const int& sy);
   int TakeLog10(const int& flag);
   Dataset GetHistogram(const int& np) const;
   Dataset GetHistogram(const int& np,const double& minimum,const double& maximum) const;
   Dataset GetWeighedHistogram(const int& np,const double& minimum,const double& maximum,const SMPMap& weigh) const;
   Dataset GetValuegram(const int& np,const double& minimum,const double& maximum) const;
   Dataset GetRadiagram(const vect& pos,const int& np,const double& maximum) const;
   Dataset GetRadiagram(const vect& pos,const int& np,const double& minimum,const double& maximum) const;
   double Average() const;
   double StdDev() const;
   double Sum() const;
   double Sum(const double& thresh) const;
   double SumInCircle(const vect& position,const double& radius) const;
   double FindSumLevel(const double& val,const double& precin);
   double AverageInCircle(const vect& position,const double& radius) const;
   double GetMaxAverage(const double& radius,vect& where,int stepsize=-1) const;
   int CopyTo(SMPMap& target,const vect& pos) const;
   int BltTo(SMPMap& target,const vect& pos) const;
   int BltTo(SMPMap& target,const int& sx,const int& sy) const;
   int RemoveSpikes(const double& thresh);
   friend SMPMap Subsection(const SMPMap& mapin,const double& lx,const double& ly,const vect& mid);
   friend SMPMap Subsection(const SMPMap& mapin,const int& sx,const int& ey,const int& sy,const int& ey);
   SMPMap subsection(const int& sx,const int& ey,const int& sy,const int& ey);
   friend SMPMap Regrid(const SMPMap& mapin,const int& newnx,const int& newny,int flag=0);
   SMPMap Regrid(const int& newnx,const int& newny,int flag=0);
   friend SMPMap Overlap(const SMPMap &map1,const SMPMap &map2);
   void AttachMap(const SMPMap& map);
         void AddMap(const SMPMap& map);
   int SetVerbosity(const int& level) {verbosity=level;}
   int GetVerbosity() {return verbosity;}
   //friend SMPMap TangentialComponent(SMPMap& map1,SMPMap& map2,vect v=vect(0.0,0.0));
         int Rotate(const double& phi,int flag=0);
         int Stretch(vect amount,int flag);
         int Stretch(const double& amount,const double& phi,int flag=0);
         int Stretch2(const double& amount,const double& phi,int flag=0);
 /*
 
 
 ********** These are declaration that should be in an interface class or so.
   friend SMPMap ConvolveMaps(SMPMap& map1,SMPMap& map2);
   friend SMPMap PowerMap(SMPMap& map1);
   friend SMPMap PowerMap(SMPMap& map1,SMPMap& map2);
   friend void GetRedShears(const SMPMap& mapin,int res,SMPMap** pe1,SMPMap** pe2);
   friend SMPMap GetKappa(SMPMap& mapin,int res);
   friend SMPMap GetKappa(SMPMap& e1,SMPMap& e2);
   friend SMPMap GetTPC(const SMPMap& mapin,int res,double min,double max);
   friend SMPMap GetKappaWithNoise(SMPMap& mapin,int res,double noise);
   friend SMPMap TwoPointCorrelation(SMPMap& map1);
   friend SMPMap TwoPointCorrelation(SMPMap& map1,SMPMap& map2);
   friend int CalculateF(double& result,SMPMap& data,SMPMap& predict,SMPMap& errors,SMPMap& model,SMPMap& defaultmap,double alpha);
   friend int CalculateXi(double& result,SMPMap& data,SMPMap& predict,SMPMap& errors); 
   friend int MateMaps(SMPMap& father,SMPMap& mother,SMPMap& child,double mutprop,double mutsig); 
   SMPMap SecondOrderMoments();
   void GetPowerMap(int inverse);
   void PowerMake(double  (*Pfunction)(double k));
 
   */
   ~SMPMap();
 };
 
 
 inline vect SMPMap::Pixel2vect(const int& x,const int& y) const
 {
   vect result;
   result.xval=(x*(maxx-minx))/nx+minx;
   result.yval=(y*(maxy-miny))/ny+miny;
   return result;
 }
 
 inline vect SMPMap::PixelMidpoint(const int& x,const int& y) const
 {
   vect result;
   result.xval=((x+0.5)*(maxx-minx))/nx+minx;
   result.yval=((y+0.5)*(maxy-miny))/ny+miny;
   return result;
 }
 
 inline int SMPMap::Vect2pixel_x(const vect& position) const
 {
   return ((int)floor(nx*(position.xval-minx)/(maxx-minx)));
 }
 inline int SMPMap::Vect2pixel_y(const vect& position) const
 {
   return ((int)floor(ny*(position.yval-miny)/(maxy-miny)));
 }
 
 inline double SMPMap::GetPoint(const vect& position) const {
   double result;
   int x=Vect2pixel_x(position);
   int y=Vect2pixel_y(position);
   result=GetPixel(x,y);
   return result;
 }
 
 inline int SMPMap::SetPoint(const vect& pos,const double& value)
 {
   int x,y;
   x=Vect2pixel_x(pos);
   y=Vect2pixel_y(pos);
   if (x>=0 && x<nx && y>=0 && y<ny) {
     mapdata[y*nx+x]=value;
     return ok;
   }
   return out_of_bounds; // point was outside map
 }
 inline int SMPMap::AddToPoint(const vect& pos,const double& value)
 {
   int x,y;
   x=Vect2pixel_x(pos);
   y=Vect2pixel_y(pos);
   if (x>=0 && x<nx && y>=0 && y<ny) {
     mapdata[y*nx+x]+=value;
     return 0;
   }
   return 1; // point was outside map
 }
 
 #endif
 
 
 
 
 
 