// discSt.cpp assign disccrete state to each amino acid in a *.ang file

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <cmath>
#include <map>
#include <cstring>
using namespace std;

#define pi 3.1415926

typedef map <string, int > SIMAP;
typedef map <string, string > SSMAP;
SSMAP AAMap;
SIMAP AIMap;
SIMAP::iterator aiItr;
SSMAP::iterator aaItr;
void editMap();
void Usage();
void at_co(double alpha, double tau, double& x, double& y, double& z);

int main(int argc, char* argv[])
{
  int i, j, k, NC=4, PI=0,aaP,posn;
  char fileName[200],OutFile[200],temp[200],caa[5];
  string saa,fname;
  double x,y,state[20][8][3];
  strcpy(fileName, "at.dat");
  strcpy(OutFile,"out.tmp");
  editMap();
  ifstream infile;
  ofstream outfile;

  i=1;
  if(argc==1){
    Usage();
  }
  while(i<argc) {
    if(strcmp(argv[i],"-h")==0) {
      Usage();
      i++;
    }
    else if(strcmp(argv[i],"-nc")==0) {
      if(argc==(i+1))
        Usage();
      NC=atoi(argv[i+1]);
      i=i+2;
    }
    else if(strcmp(argv[i],"-if")==0) {
      if(argc==i+1)
	Usage();
      strcpy(fileName,argv[i+1]);
      fname=fileName;
	i+=2;
    }
    else if(strcmp(argv[i],"-pi")==0) {
      PI=1;
      i++;
    }
  }

  strcat(fileName,".ang");
  fname=fname+".sts";

  outfile.open(fname.c_str());
  outfile<<fname.substr(0,fname.size()-4).c_str()<<" ";
  if(NC==4)
    infile.open("/mnt/autofs/home/jinfeng/SharedData/cls_nc4.dat");
  else if(NC==5)
    infile.open("/mnt/autofs/home/jinfeng/SharedData/cls_nc5.dat");
  else if(NC==6)
    infile.open("/mnt/autofs/home/jinfeng/SharedData/cls_nc6.dat");
  if(!infile.is_open()) {
    cout<<"Erro in opening file cls_nc4.dat."<<endl;
      exit(0);
  }
  i=0;
  while(!infile.eof()&&i<20){
    infile>>saa;
    if((aiItr=AIMap.find(saa))!=AIMap.end())
      posn=(*aiItr).second;
    else {
      cout<<" can not find the amino acid type "<<saa<<endl;
    }
    for(j=0;j<NC;++j) {
      infile>>x>>y;
      at_co(x,y,state[posn][j][0],state[posn][j][1],state[posn][j][2]);
    }
    ++i;
  }
  infile.close();
  infile.open(fileName);
  if(!infile.is_open()) {
    fname=fname.substr(0,fname.size()-4);
    strcpy(fileName,"/mnt/autofs/home/jinfeng/SharedData/pc30_res2.0_ca_ang/");
    strcat(fileName,fname.c_str());
    strcat(fileName,".ca.ang");
    infile.open(fileName);
    if(!infile.is_open()) {
      cout<<"Erro in opening file "<<fileName<<endl;
      exit(0);
    }
  }
  float mindis;
  float mintemp;
  int minid;
  double tdx,tdy,tdz;
  while(!infile.eof()) {
    infile.getline(temp,200);
    if(strlen(temp)<10) break;
    sscanf(temp,"%s %d %lf %lf",caa,&aaP,&x,&y);
    saa=caa;
    if((aiItr=AIMap.find(saa))!=AIMap.end())
      posn=(*aiItr).second;
    else {
      cout<<"can not find the amino acid type "<<temp<<endl;
      exit(0);
    }
    mindis=100000;
    minid=0;
    at_co(x,y,tdx,tdy,tdz);
    for(i=0;i<NC;++i) {
	if((mintemp=(pow((tdx-state[posn][i][0]),2)+pow(tdy-state[posn][i][1],2)+pow(tdz-state[posn][i][2],2)))<mindis) {
	  mindis=mintemp;
	  minid=i;
	}
    }
    outfile<<minid;
  }
  outfile<<endl;
  return 0;
}

void at_co(double alpha, double tau, double& x, double& y, double& z)
{
  x=cos(pi-pi*alpha/180);
  y=-1*sin(pi-pi*alpha/180)*cos(pi*tau/180);
  z=-1*sin(pi-pi*alpha/180)*sin(pi*tau/180);
}

void editMap()
{
    AAMap.insert(SSMAP::value_type("Ala", "A"));
    AAMap.insert(SSMAP::value_type("Cys", "C"));
    AAMap.insert(SSMAP::value_type("Asp", "D"));
    AAMap.insert(SSMAP::value_type("Glu", "E"));
    AAMap.insert(SSMAP::value_type("Phe", "F"));
    AAMap.insert(SSMAP::value_type("Gly", "G"));
    AAMap.insert(SSMAP::value_type("His", "H"));
    AAMap.insert(SSMAP::value_type("Ile", "I"));
    AAMap.insert(SSMAP::value_type("Lys", "K"));
    AAMap.insert(SSMAP::value_type("Leu", "L"));
    AAMap.insert(SSMAP::value_type("Met", "M"));
    AAMap.insert(SSMAP::value_type("Asn", "N"));
    AAMap.insert(SSMAP::value_type("Pro", "P"));
    AAMap.insert(SSMAP::value_type("Gln", "Q"));
    AAMap.insert(SSMAP::value_type("Arg", "R"));
    AAMap.insert(SSMAP::value_type("Ser", "S"));
    AAMap.insert(SSMAP::value_type("Thr", "T"));
    AAMap.insert(SSMAP::value_type("Val", "V"));
    AAMap.insert(SSMAP::value_type("Trp", "W"));
    AAMap.insert(SSMAP::value_type("Tyr", "Y"));
    AAMap.insert(SSMAP::value_type("A", "Ala"));
    AAMap.insert(SSMAP::value_type("C", "Cys"));
    AAMap.insert(SSMAP::value_type("D", "Asp"));
    AAMap.insert(SSMAP::value_type("E", "Glu"));
    AAMap.insert(SSMAP::value_type("F", "Phe"));
    AAMap.insert(SSMAP::value_type("G", "Gly"));
    AAMap.insert(SSMAP::value_type("H", "His"));
    AAMap.insert(SSMAP::value_type("I", "Ile"));
    AAMap.insert(SSMAP::value_type("K", "Lys"));
    AAMap.insert(SSMAP::value_type("L", "Leu"));
    AAMap.insert(SSMAP::value_type("M", "Met"));
    AAMap.insert(SSMAP::value_type("N", "Asn"));
    AAMap.insert(SSMAP::value_type("P", "Pro"));
    AAMap.insert(SSMAP::value_type("Q", "Gln"));
    AAMap.insert(SSMAP::value_type("R", "Arg"));
    AAMap.insert(SSMAP::value_type("S", "Ser"));
    AAMap.insert(SSMAP::value_type("T", "Thr"));
    AAMap.insert(SSMAP::value_type("V", "Val"));
    AAMap.insert(SSMAP::value_type("W", "Trp"));
    AAMap.insert(SSMAP::value_type("Y", "Tyr"));

  // initialize AIMap
    AIMap.insert(SIMAP::value_type("A",0));
    AIMap.insert(SIMAP::value_type("C",1));
    AIMap.insert(SIMAP::value_type("D",2));
    AIMap.insert(SIMAP::value_type("E",3));
    AIMap.insert(SIMAP::value_type("F",4));
    AIMap.insert(SIMAP::value_type("G",5));
    AIMap.insert(SIMAP::value_type("H",6));
    AIMap.insert(SIMAP::value_type("I",7));
    AIMap.insert(SIMAP::value_type("K",8));
    AIMap.insert(SIMAP::value_type("L",9));
    AIMap.insert(SIMAP::value_type("M",10));
    AIMap.insert(SIMAP::value_type("N",11));
    AIMap.insert(SIMAP::value_type("P",12));
    AIMap.insert(SIMAP::value_type("Q",13));
    AIMap.insert(SIMAP::value_type("R",14));
    AIMap.insert(SIMAP::value_type("S",15));
    AIMap.insert(SIMAP::value_type("T",16));
    AIMap.insert(SIMAP::value_type("V",17));
    AIMap.insert(SIMAP::value_type("W",18));
    AIMap.insert(SIMAP::value_type("Y",19));

    AIMap.insert(SIMAP::value_type("ALA",0));
    AIMap.insert(SIMAP::value_type("CYS",1));
    AIMap.insert(SIMAP::value_type("ASP",2));
    AIMap.insert(SIMAP::value_type("GLU",3));
    AIMap.insert(SIMAP::value_type("PHE",4));
    AIMap.insert(SIMAP::value_type("GLY",5));
    AIMap.insert(SIMAP::value_type("HIS",6));
    AIMap.insert(SIMAP::value_type("ILE",7));
    AIMap.insert(SIMAP::value_type("LYS",8));
    AIMap.insert(SIMAP::value_type("LEU",9));
    AIMap.insert(SIMAP::value_type("MET",10));
    AIMap.insert(SIMAP::value_type("ASN",11));
    AIMap.insert(SIMAP::value_type("PRO",12));
    AIMap.insert(SIMAP::value_type("GLN",13));
    AIMap.insert(SIMAP::value_type("ARG",14));
    AIMap.insert(SIMAP::value_type("SER",15));
    AIMap.insert(SIMAP::value_type("THR",16));
    AIMap.insert(SIMAP::value_type("VAL",17));
    AIMap.insert(SIMAP::value_type("TRP",18));
    AIMap.insert(SIMAP::value_type("TYR",19));

}

void Usage()
{
  cout<<"assign_disAngle, program to assign discrete state to a protein structure."<<endl;
  cout<<"options: "<<endl;
  cout<<"\t-h output this help file."<<endl;
  cout<<"\t-nc number of discrete state."<<endl;
  cout<<"\t-if input filename."<<endl;
  exit(0);
}