#include <Map2D.h>
#include <ablator.h>

double xdelta = 0.010;
double ydelta = 0.010;
double xview[2] = { -1, 8 };
double yview[2] = { -1, 8 };
double xdims[2] = {  0, 7.2 };
double ydims[2] = {  0, 7.2 };
double xhole[2] = { 0.4, 6.8 };
double yhole[2] = { 0.4, 6.8 };
double thick[2] = { 300.,20. };

void run_ablator(char *file)
{
   int nxbins = (xview[1]-xview[0])/xdelta;
   int nybins = (yview[1]-yview[0])/ydelta;
   Map2D *virgin = new Map2D(TH2D("virgin","virgin surface",
                                  nxbins,xview[0],xview[1],
                                  nybins,yview[0],yview[1]));
   Map2D *target = new Map2D(TH2D("target","target surface",
                                  nxbins,xview[0],xview[1],
                                  nybins,yview[0],yview[1]));

   for (int ix = 1; ix <= nxbins; ++ix) {
      for (int iy = 1; iy <= nybins; ++iy) {
         double x = virgin->GetXaxis()->GetBinCenter(ix);
         double y = virgin->GetYaxis()->GetBinCenter(iy);
         if (x > xhole[0] && x < xhole[1] &&
             y > yhole[0] && y < yhole[1])
         {
            target->SetBinContent(ix,iy,thick[1]);
            virgin->SetBinContent(ix,iy,thick[0]);
         }
         else if (x > xdims[0] && x < xdims[1] &&
                  y > ydims[0] && y < ydims[1])
         {
            target->SetBinContent(ix,iy,thick[0]);
            virgin->SetBinContent(ix,iy,thick[0]);
         }
         else {
            virgin->SetBinContent(ix,iy,0);
            target->SetBinContent(ix,iy,thick[0]);
         }
      }
   }

   ablator abl;
   Map2D *targets = abl.smooth(target);
   abl.setRaster(xdelta,ydelta);
   abl.setInterlace(0);
   Map2D *part = abl.mill(virgin,targets,file);
}

void run_ablator_map(Map2D *part, Map2D *target, char *file)
{
   ablator abl;
   //abl.setRaster(xdelta,ydelta);
   Map2D *part = abl.mill_sparse(part,target,file);
}