//
// ablator: ROOT c++ class implementing the capability to analyze a surface
//          profile image of a diamond sample (eg. from a Zygo microscope)
//          and create a laser ablation sequence file that can be used to
//          mill that sample down to match an arbitrary target surface.
//
// author: richard.t.jones at uconn.edu
// version: february 10, 2015
//

#include "Map2D.h"

class ablator {
 public:
   ablator():
      // general spot model guidance parameters
      fSpotHsize(0.350),
      fSpotVsize(0.350),
      fSpotHresol(0.010), // mm
      fSpotVresol(0.010), // mm

      // parameters for a single 2d gaussian model for the laser ablation spot
      // from a single laser pulse, averaged over a large number of pulses
      fSpotHrms(0.067), // mm
      fSpotVrms(0.047), // mm
      fSpotDepth(0.183), // microns

      // parameters for an empirical 2d model for the laser ablation spot
      // from a single laser pulse, averaged over a large number of pulses
      fSpotPitDepth(0.130), // microns
      fSpotHbase(0.030), // microns
      fSpotHsigma(0.077), // microns
      fSpotVsigma0(0.013), // mm
      fSpotVsigma1(0.048), // mm
      fSpotVampli0(0.072), // microns
      fSpotVampli1(0.058), // microns
      fSpotVsigma2(0.036), // mm
      fSpotVoffset2(0.079), // mm
      fSpotVampli2(0.023), // mm
      fSpotVbase2(0.047), // mm

      // parameters that control aspects of the ablator algorithms
      fSpotDepthRMS(0.15), // relative units
      fXedgeParam(0.070), // mm
      fYedgeParam(0.050), // mm
      fRaster_dx(0.010), // mm
      fRaster_dy(0.010), // mm
      fTaperFactor(1.0), // multiplies steady-state pulse energy
      fTaperLength(10.), // pulses
      fReps(10) {}

   int getReps() const {
      return fReps;
   }

   double getRaster_dx() const {
      return fRaster_dx;
   }
   double getRaster_dy() const {
      return fRaster_dy;
   }
   int getSpotHsize() const {
      return fSpotHsize;
   }
   int getSpotVsize() const {
      return fSpotVsize;
   }
   double getSpotHresol() const {
      return fSpotHresol;
   }
   double getSpotVresol() const {
      return fSpotVresol;
   }

   double getSpotHrms() const {
      return fSpotHrms;
   }
   double getSpotVrms() const {
      return fSpotVrms;
   }
   double getSpotDepth() const {
      return fSpotDepth;
   }

   double getSpotPitDepth() const {
      return fSpotPitDepth;
   }
   double getSpotHbase() const {
      return fSpotHbase;
   }
   double getSpotHsigma() const {
      return fSpotHsigma;
   }
   double getSpotVsigma0() const {
      return fSpotVsigma0;
   }
   double getSpotVsigma1() const {
      return fSpotVsigma1;
   }
   double getSpotVampli0() const {
      return fSpotVampli0;
   }
   double getSpotVampli1() const {
      return fSpotVampli1;
   }
   double getSpotVsigma2() const {
      return fSpotVsigma2;
   }
   double getSpotVoffset2() const {
      return fSpotVoffset2;
   }
   double getSpotVampli2() const {
      return fSpotVampli2;
   }
   double getSpotVbase2() const {
      return fSpotVbase2;
   }

   double getSpotDepthRMS() const {
      return fSpotDepthRMS;
   }
   double getXedgeParam() const {
      return fXedgeParam;
   }
   double getYedgeParam() const {
      return fYedgeParam;
   }
   double getTaperFactor() const {
      return fTaperFactor;
   }
   double getTaperLength() const {
      return fTaperLength;
   }
   void setReps(int reps) {
      fReps = reps;
   }
   void setRaster(double dx, double dy) {
      fRaster_dx = dx;
      fRaster_dy = dy;
   }

   void setSpotSize(double xsize, double ysize) {
      fSpotHsize = xsize;
      fSpotVsize = ysize;
   }
   void setSpotResol(double hresol, double vresol) {
      fSpotHresol = hresol;
      fSpotVresol = vresol;
   }

   void setSpotRMS(double xrms, double yrms) {
      fSpotHrms = xrms;
      fSpotVrms = yrms;
   }
   void setSpotDepth(double depth) {
      fSpotDepth = depth;
   }

   void setSpotPitDepth(double depth) {
      fSpotPitDepth = depth;
   }
   void setSpotHbase(double base) {
      fSpotHbase = base;
   }
   void setSpotHsigma(double sigma) {
      fSpotHsigma = sigma;
   }
   void setSpotVsigma0(double sigma) {
      fSpotVsigma0 = sigma;
   }
   void setSpotVsigma1(double sigma) {
      fSpotVsigma1 = sigma;
   }
   void setSpotVampli0(double amp) {
      fSpotVampli0 = amp;
   }
   void setSpotVampli1(double amp) {
      fSpotVampli1 = amp;
   }
   void setSpotVsigma2(double sigma) {
      fSpotVsigma2 = sigma;
   }
   void setSpotVoffset2(double offset) {
      fSpotVoffset2 = offset;
   }
   void setSpotVampli2(double amp) {
      fSpotVampli2 = amp;
   }
   void setSpotVbase2(double base) {
      fSpotVbase2 = base;
   }

   void setSpotDepthRMS(double fractional_rms) {
      fSpotDepthRMS = fractional_rms;
   }
   void setSlopeParams(double xwidth, double ywidth) {
      fXedgeParam = xwidth;
      fYedgeParam = ywidth;
   }
   void setTaperFactor(double taper_factor) {
      fTaperFactor = taper_factor;
   }
   void setTaperLength(double taper_pulses) {
      fTaperLength = taper_pulses;
   }

   double get_cut_depth() const;
   double get_max_pit_depth() const;
   Map2D *make_beam_spot(const Map2D *templ) const;
   Map2D *make_gaussian_spot(const Map2D *templ, 
                             double sigx, double sigy) const;
   Map2D *smooth(Map2D *target_map, double smear_factor = 0);
   Map2D *unsmooth(const Map2D *target_map, double noise_cutoff = 1e+12);
   Map2D *unsmooth2(const Map2D *target_map, double eigenvalue_cutoff = 1e-12);
   Map2D *fillvoids(const Map2D *part_map, double max_void_sigmas = 3);
   Map2D *mill(const Map2D *part_map, Map2D *target_map);
   Map2D *discretize(const Map2D *pulse_map);
   Map2D *detrench(const Map2D *program_map);
   Map2D *simulate_mill(const Map2D *program_map, Map2D *part_map);
   int rasterize(const Map2D *program_map, const char *seqfile);

 private:
   // general spot model guidance
   double fSpotHsize;
   double fSpotVsize;
   double fSpotHresol;
   double fSpotVresol;

   // naive gaussian beam spot model parameters
   double fSpotHrms;
   double fSpotVrms;
   double fSpotDepth;

   // empirical non-gaussian beam spot parameterization
   double fSpotPitDepth;
   double fSpotHbase;
   double fSpotHsigma;
   double fSpotVsigma0;
   double fSpotVsigma1;
   double fSpotVampli0;
   double fSpotVampli1;
   double fSpotVsigma2;
   double fSpotVoffset2;
   double fSpotVampli2;
   double fSpotVbase2;

   // parameters for ablator algorithms
   double fSpotDepthRMS;
   double fXedgeParam;
   double fYedgeParam;
   double fRaster_dx;
   double fRaster_dy;
   double fTaperFactor;
   double fTaperLength;
   int fReps;

   void check_orthogonality(const TMatrixD& matr);
};

class Pulse {
 public:
   Pulse() : ix(0), iy(0) {}
   int ix;
   int iy;
};

inline bool operator<(Pulse &p1, Pulse &p2) {
   if (p1.iy != p2.iy) {
      return (p1.iy < p2.iy);
   }
   else {
      return (p1.ix < p2.ix);
   }
}