/*************************************************************************
 *
 * HTCd - Copyright (C) 1998-2006 Henrik Rydberg
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include <mt/config.h>
#include <fstream>
#include <iostream>
#include <unistd.h>
#include <fcntl.h>
#include <signal.h>
#include <sys/stat.h>
using namespace std;

const mstring LABAT=htchome()+sref("/bin/labat");

const int SLOTS=8;
const int RUNLIMIT=600;
static int SLOT;

extern "C" {
static void FreeSlot(int dummy)
{
  static int ok=0;
  static char buf[256];
  if(!ok) {
    ok=1;
    sprintf(buf,"rm -rf /tmp/labat%d/",SLOT);
    system(buf);
  }
}
}

static int AllocSlot()
{
  struct stat FS;
  char buf[256];
  for(SLOT=0;SLOT<SLOTS;SLOT++) {
    sprintf(buf,"/tmp/labat%d/",SLOT);
    if(stat(buf,&FS)!=0) {
      signal(SIGTERM,FreeSlot);
      signal(SIGINT,FreeSlot);
      signal(SIGPIPE,FreeSlot);
      sprintf(buf,"mkdir /tmp/labat%d/",SLOT);
      system(buf);
      return 1;
    }
  }
  return 0;
}

///////////////////////////////////////////////////////////////////////////////

struct CTLBasis {
  int n,l,occ;
  float e;
  void SetEnergy(int Z) {
    e=-0.5*float(Z*Z)/float(0.1+n*n*n+l*(l+1)*(l+2));
  }
};

struct CTLData {
  int atomz,its,basisn,exc;
  float mix;
  CTLBasis left[24],right[24];
};

static ostream& operator<<(ostream& out,const CTLData& ctl)
{
  out<<ctl.atomz<<" "<<0<<"\n";
  out<<10<<" "<<40<<" "<<0.2<<"\n"; // Hermann-Skillman mesh
  out<<ctl.mix<<" 0.00001 0.005\n"; // mix,error,?
  out<<ctl.its<<" "<<0<<" "<<1<<" "<<0<<"\n"; // its ? step ?
  out<<ctl.basisn<<" "<<2<<" "<<ctl.exc<<"\n"; // exc (C-A, B-H, GGA)
  for(int i=0;i<ctl.basisn;i++) {
    out<<ctl.left[i].n<<ctl.left[i].l<<0<<" "<<ctl.left[i].e
       <<" "<<ctl.left[i].occ<<" "<<1<<" "<<0<<"\n"; 
  }
  for(int i=0;i<ctl.basisn;i++) {
    out<<ctl.right[i].n<<ctl.right[i].l<<0<<" "<<ctl.right[i].e
       <<" "<<ctl.right[i].occ<<" "<<1<<" "<<0<<"\n"; 
  }
  return out;
}

static void SetupLabat(mstring& resfile,int argc,char* argv[])
{
  CTLData ctl;
  ctl.atomz=atoi(argv[1]);
  ctl.its=atoi(argv[2]);
  ctl.exc=atoi(argv[3]);
  ctl.mix=atof(argv[4]);
  resfile=argv[5];
  ctl.basisn=0;
  int nl,nr,elecs=0;
  for(int i=0;i<6;i++) {
    for(int k=0;k<3;k++) {
      nl=int(argv[6+i][2*k]-'0');
      nr=int(argv[6+i][2*k+1]-'0');
      if(nl>0||nr>0) {
	ctl.left[ctl.basisn].n=i+1;
	ctl.left[ctl.basisn].l=k;
	ctl.left[ctl.basisn].occ=nl;
	ctl.left[ctl.basisn].SetEnergy(ctl.atomz);
	ctl.right[ctl.basisn].n=i+1;
	ctl.right[ctl.basisn].l=k;
	ctl.right[ctl.basisn].occ=nr;
	ctl.right[ctl.basisn].SetEnergy(ctl.atomz);
	ctl.basisn++;
      }
      elecs+=nl+nr;
    }
  }
  ctl.atomz+=elecs;
  char buf[1024];
  sprintf(buf,"/tmp/labat%d/atomctrl.dat",SLOT);
  ofstream os(buf);
  os<<ctl;
  cout<<"Atomnumber: "<<ctl.atomz<<" Electrons: "<<elecs<<"<BR>\n";
  cout<<"NOTE: if the program does >="<<ctl.its
      <<" iterations, you must increase the number\n";
  cout<<"of iterations and re-run, to assure convergence.<BR>\n";
  cout<<"Also note that there is currently a maximum time limit of "<<
    RUNLIMIT<<" seconds.<BR>\n";
}

static void RunLabat()
{
  char buf[1024];
  sprintf(buf,"cd /tmp/labat%d/\nulimit -t %d\n%s",SLOT,RUNLIMIT,LABAT.c_str());
  system(buf);
}

static void SaveLabat(const mstring& resfile)
{
  const char* docroot=getenv("DOCROOT");
  if(docroot&&resfile.size()) {
    char buf[1024];
    sprintf(buf,"cp -rfp /tmp/labat%d/atomdens.dat %s/%s",
            SLOT,docroot,resfile.c_str());
    system(buf);
    cout<<"<p>\n"
        <<"Here you may download the density profile.\n"
        <<"Example: Save the file as 'dens.dat', and load it from\n"
        <<"matlab with 'load dens.dat'.\n"
        <<"Then you can look at it with 'plot(dens(:,1),dens(:,2))'.\n"
        <<"<p>\n"<<"<a href="<<resfile.c_str()<<">Get atomic density</a>\n";
  }
}

main(int argc,char* argv[])
{
  if(argc<12) {
    cerr<<"Usage: "<<argv[0]<<" <ion> <its> <exc> <mix> <resfile> <ssppdd>x6\n";
    exit(-1);
  }
  mstring resfile;
  char buf[1024];
  if(AllocSlot()) {
    cout<<"<h4>Running</h4>\n"
        <<"<PRE>\n";
    cout.flush();
    SetupLabat(resfile,argc,argv);
    RunLabat();
    cout<<"</PRE>\n"
        <<"<h4>Done</h4>\n";
    SaveLabat(resfile);
    FreeSlot(0);
  }
  else {
    cout<<"No vacant slots at the moment. Try again later.\n";  
  }
}