#include "ideal_defden.h"

double IdealDefdenGenerator::defden(V3 pos,CEnvironment &belo, QList<Matrix> &mi){
  double rho=0.0;
  double argu;
  for (int i=0; i<belo.size(); i++){
    //V3 frac = V3(fmod(9999.+belo.at(i).frac.x,1.0), fmod(9999.+belo.at(i).frac.y,1.0), fmod(9999.+belo.at(i).frac.z,1.0));
    V3 pos0 = pos - belo.at(i).frac;
    argu = pos0 * (mi.at(i) * pos0);
    if (argu>-37.0)// this speeds it up!
        rho += belo.at(i).sof_org * exp(argu);
  }
  return rho;
}

void IdealDefdenGenerator::beloWriteCube(){
  if (!chgl->hideBeLo->isVisible()) {qDebug()<<"BeLo: 'rrrrrrr!'"; return; }
  CEnvironment belo,bu;
  double stepWitdth=0.11;
  int stepx = chgl->mol->cell.a / stepWitdth;
  int stepy = chgl->mol->cell.b / stepWitdth;
  int stepz = chgl->mol->cell.c / stepWitdth;  
  size_t ntot = stepx * stepy * stepz;
  float *dat = NULL;
  do {
	dat = (float*) malloc(ntot * sizeof (float));
	if (dat == NULL){	  //
		if (stepWitdth > 0.3) return; // stop if this gets to coarse
		stepWitdth *= 1.1;//try a larger stepWitdth
		stepx = chgl->mol->cell.a / stepWitdth;
		stepy = chgl->mol->cell.b / stepWitdth;
		stepz = chgl->mol->cell.c / stepWitdth;  
		ntot = stepx * stepy * stepz;
		dat = (float*) malloc(ntot * sizeof (float));    
	}  
  } while (dat == NULL);  
  V3 _dx = V3(1.0 / stepx, 0, 0), _dy = V3(0, 1.0 / stepy, 0), _dz = V3(0, 0, 1.0 / stepz);
  V3 org = 0.5 * _dx + 0.5 * _dy + 0.5 * _dz;
  for (int bl = 0; bl < chgl->mol->asymm.size(); bl++){
    if (chgl->mol->asymm.at(bl).an==-42){  //bede lone objects
      bu.append(chgl->mol->asymm[bl]);
    }
  }
  belo = chgl->mol->packInLimits(bu,-0.1,1.1,-0.1,1.1,-0.1,1.1);
  emit message(QString("Computing on %1 BEDE LONE objects.").arg(belo.size()));
  bu.clear();
  double midd=9.9e37,madd=-9.9e37,rho,summ=0.0,sum1=0.0, sum2=0.0;
  QList<Matrix> mi;
  const double pisqrt = sqrt(M_PI);
  const double pi2=M_PI*M_PI;
  const double pisqrt3=M_PI*pisqrt;
  Matrix M,Mi;
  double d;
  for (int i=0; i<belo.size(); i++){
    M = (chgl->mol->cell.Gi* (belo.at(i).peakHeight*4)) + belo.at(i).uf;
    //M = (chgl->mol->cell.Gi* (belo.at(i).uf*belo.at(i).peakHeight*.25));
    //M = (chgl->mol->cell.Gi* ((belo.at(i).uf)*(belo.at(i).peakHeight*0.25)));//uf is Bij here
    //M = (chgl->mol->cell.Gi* (belo.at(i).peakHeight*0.25)) * belo.at(i).uf;//uf is Bij here
    Mi = inverse(M);
    //sof==a peakHeight=b
    d = determinant(Mi);
    belo[i].sof_org = sqrt(d) * belo.at(i).sof  / chgl->mol->cell.V;// * pisqrt3;
    Mi = Mi *(-pi2);
    mi.append(Mi);
  } 
  uint pro=0,aproc,bproc=10;
  for (int iz=0; iz<stepz; iz++){
    for (int iy=0; iy<stepy; iy++){
      for (int ix=0; ix<stepx; ix++){
        rho=defden((ix*_dx)+(iy*_dy)+(iz*_dz)+org,belo,mi);//+minp
        midd=qMin(rho,midd);
        madd=qMax(rho,madd);
        if (rho>0.0) sum1+=rho;
        else sum2+=rho;
        summ+=rho;
        if (fabs(rho)<=1.0e-98) rho=0.0;
        dat[ix + iy*stepx + iz*stepx*stepy]=(float)rho;
        pro++;
      }
    }
    aproc = (pro * 100) / ntot;
    if (aproc!=bproc){
        emit message(QString("Calculating deformation density %1% (steps:%2)").arg(aproc).arg(ntot));
        bproc = aproc;
    }    
  }
  mi.clear();
  belo.clear();
  printf("%f %f %f \n",midd,madd,summ);
  bool gut = fxle->loadMap(stepx, stepy, stepz, dat);//dat lives as datfo_fc its life in fxle 
  if (!gut) {
	  emit message("Oh no! Something went wrong while creating BODD maps. Perhaps we are out of memory. Sorry!");
	  free(dat);
	  dat=NULL;
	  return;
  }
  printf("%f %f %f %d\n",midd, madd, summ, gut);
  chgl->dimc=QColor("darkorange");
  chgl->dipc=QColor("teal");
  emit message(QString("finished! minRho %1 e/A^3 maxRho %2 e/A^3 sum %3 e/A^3 OK %4  %5 %6 ").arg(midd).arg(madd).arg(summ).arg(gut).arg(sum1).arg(sum2));
}