/////////////////////////////////////////////////////////////////////////////
/// movemol.cpp - modify 3d coords
///
/// Author: Jeremy Yang
///  
/////////////////////////////////////////////////////////////////////////////
#include "openeye.h"
#include <iostream>
#include <vector>
#include <math.h>
#include "oesystem.h"
#include "oechem.h"

using namespace OESystem;
using namespace OEChem;
using namespace OEMath;
using namespace std;

static string ifile="";
static string ofile="";
static bool center=false;
static bool mcmols=false;
static float rotx=0.0;
static float roty=0.0;
static float rotz=0.0;
static double trans[3]={0.0,0.0,0.0};
static int verbose=0;


int ParseCommandLine(int argc,char *argv[]);
double *RotMatAxial(char axis,double a);

void Help(string msg)
{
cerr << msg << endl
     <<"\t************************************************************"<<endl
     <<"\t| movemol [options]                                        |"<<endl
     <<"\t|   options:                                               |"<<endl
     <<"\t|     -i INFILE                                            |"<<endl
     <<"\t|     -o OUTFILE                                           |"<<endl
     <<"\t|     -center   ... center mol at origin                   |"<<endl
     <<"\t|     -rotx AX ... rotate about x axis AX degrees          |"<<endl
     <<"\t|     -roty AY ... rotate about y axis AY degrees          |"<<endl
     <<"\t|     -rotz AZ ... rotate about z axis AZ degrees          |"<<endl
     <<"\t|     -trans X,Y,Z ... translate (X,Y,Z)                   |"<<endl
     <<"\t| Note that multiple operations will be done in the order  |"<<endl
     <<"\t| listed here.                                             |"<<endl
     <<"\t|     -mcmols ... handle as multiconf mols                 |"<<endl
     <<"\t|     -v .......  verbose                                  |"<<endl
     <<"\t|     -vv ......  very verbose                             |"<<endl
     <<"\t|     -help                                                |"<<endl
     <<"\t************************************************************"<<endl;
  exit(1);
}


int main(int argc,char *argv[])
{
  oemolistream ims;
  oemolostream oms;
  // OEMol mol;
  OEGraphMol mol;
  string smi;
  double *rot;
  double xyz[3];
  int molcount;

  ParseCommandLine(argc,argv);

  if (!ims.open(ifile))
    OEThrow.Fatal("Unable to open \"%s\".\n",ifile.c_str());
  if (!ofile.empty() && !oms.open(ofile))
    OEThrow.Fatal("Unable to open \"%s\".\n",ofile.c_str());

  if (mcmols) ims.SetConfTest(OEAbsoluteConfTest(false));

  for (molcount=0; OEReadMolecule(ims,mol); ++molcount)
  {
    if (verbose)
    {
      OECreateIsoSmiString(smi,mol);
      cerr << smi <<endl;
      cerr << "formal dimension of mol: " << mol.GetDimension() <<endl;
    }

    OESetDimensionFromCoords(mol);
    if (mol.GetDimension()<1) 
    {
      cerr << "No coordinates -- skipping mol." <<endl;
      continue;
    }

    if (center)
    {
      OECenter(mol,xyz);  // OECenter() centers and fills xyz.
      if (verbose)
        fprintf(stderr,"centered from: %.2f,%.2f,%.2f\n",xyz[0],xyz[1],xyz[2]);
    }

    if (rotx!=0.0)
    {
      if (verbose) fprintf(stderr,"x-rotating: %f deg\n",rotx);
      rot = RotMatAxial('x',rotx*OEMath::Pi/180.0);
      OERotate(mol,rot);
    }
    if (roty!=0.0)
    {
      if (verbose) fprintf(stderr,"y-rotating: %f deg\n",roty);
      rot = RotMatAxial('y',roty*OEMath::Pi/180.0);
      OERotate(mol,rot);
    }
    if (rotz!=0.0)
    {
      if (verbose) fprintf(stderr,"z-rotating: %f deg\n",rotz);
      rot = RotMatAxial('z',rotz*OEMath::Pi/180.0);
      OERotate(mol,rot);
    }

    if (trans[0]!=0.0 || trans[1]!=0.0 || trans[2]!=0.0)
    {
      if (verbose) fprintf(stderr,"trans: %f,%f,%f\n",trans[0],trans[1],trans[2]);
      OETranslate(mol,trans);
    }


    OEWriteMolecule(oms,mol);
  }
  return 0;
}


/////////////////////////////////////////////////////////////////////////////
int ParseCommandLine(int argc,char *argv[])
{
  string a;
  if (argc==1) Help("");
  for (++argv;--argc;++argv)
  {
    a=*argv;
    if (a=="-i")
    {
      if (--argc==0) Help("ERROR: "+a+" requires value.");
      ifile=*++argv;
    }
    else if (a=="-o")
    {
      if (--argc==0) Help("ERROR: "+a+" requires value.");
      ofile=*++argv;
    }
    else if (a=="-center") center=true;
    else if (a=="-mcmols") mcmols=true;
    else if (a=="-rotx")
    {
      if (--argc==0) Help("ERROR: "+a+" requires value.");
      rotx=atof(*++argv);
    }
    else if (a=="-roty")
    {
      if (--argc==0) Help("ERROR: "+a+" requires value.");
      roty=atof(*++argv);
    }
    else if (a=="-rotz")
    {
      if (--argc==0) Help("ERROR: "+a+" requires value.");
      rotz=atof(*++argv);
    }
    else if (a=="-trans")
    {
      if (--argc==0) Help("ERROR: "+a+" requires value.");
      sscanf(*++argv,"%lf,%lf,%lf",trans,trans+1,trans+2);
    }
    else if (a=="-v") verbose=1;
    else if (a=="-vv") verbose=2;
    else if (a=="-h") Help("");
    else Help("ERROR: Bad option: "+a);
  }
  return 1;
}
  
/////////////////////////////////////////////////////////////////////////////
///  RotMatAxial - returns axial rotation matrices (clockwise looking
///  into axis).
/// 
/// Refs: http://mathworld.wolfram.com/RotationMatrix.html,
/// Goldstein 1980, pp. 146-147 and 608; Arfken 1985, pp. 199-200
/////////////////////////////////////////////////////////////////////////////
double *RotMatAxial(char axis,double a)
{
  static double mat[9];

  switch (axis) 
  {
    case 'x':
      mat[0]=1.0;      mat[1]=1.0;       mat[2]=0.0;
      mat[3]=0.0;      mat[4]=cos(a);    mat[5]=sin(a);
      mat[6]=0.0;      mat[7]=-sin(a);   mat[8]=cos(a);
      break;
    case 'y':
      mat[0]=cos(a);   mat[1]=0.0;       mat[2]=-sin(a);
      mat[3]=0.0;      mat[4]=1.0;       mat[5]=0.0;    
      mat[6]=sin(a);   mat[7]=0.0;       mat[8]=cos(a);
      break;
    default:
      mat[0]=cos(a);   mat[1]=sin(a);    mat[2]=0.0;    
      mat[3]=-sin(a);  mat[4]=cos(a);    mat[5]=0.0;    
      mat[6]=0.0;      mat[7]=0.0;       mat[8]=1.0;   
  }
  return mat;
}