#ifndef PDB_H
#define PDB_H

#include <iostream>
#include <fstream>
#include <math.h>
#include <string>
#include <vector>
#include <map>

#define PI              3.14159265358979323846
#define RAD             57.29577951
#define MAXNUM          9999

#define max(a,b) (a>b)?a:b
#define min(a,b) (a>b)?b:a

using namespace std;
using std::string;

typedef float Vec3[3];
typedef float Mat3[3][3];
typedef map<int, map<int, double> > DoubleD2;
typedef map<int, double>  DoubleD1;

typedef struct CORdata {
  int atomNum, resNum;
  string atomName, resName, chainName;
  float X, Y, Z, B_Factor;
  Vec3 Coord;
} PDB_Entry;


typedef struct {
	int resID, atomNo;
	string resName;
	PDB_Entry atoms[6];

	Vec3 coordA[6];
	Vec3 center, norm;
	Mat3 transM;
	float rad;

	float ringFact;
} RingData;

typedef map<int, PDB_Entry> Mol; //atom list indexed by atom number
typedef map<int, Mol> Mols;


class PDB
{
	long sgn(float x);
	float arccos_(float x);
	void Vec3Zero(Vec3 v);
	void Vec3Copy(Vec3 v1, Vec3 v2);
	float Vec3Abs(Vec3 v);
	float Vec3DiffAbs(Vec3 v1, Vec3 v2);
	void Vec3Norm(Vec3 v);
	void Vec3Scale(Vec3 v, float s);
	void Vec3Add(Vec3 v1, Vec3 v2);
	void Vec3Sub(Vec3 v1, Vec3 v2);
	float Vec3Scalar(Vec3 v1, Vec3 v2);
	void Vec3Cross(Vec3 v1, Vec3 v2);
	void Mat3VecMult(Vec3 v, Mat3 m);

public:

	string PDBfileName;

	string Format;
	int VarsNumber;

	RingData Rings[200]; //Aromatic ring list
	int RingNo;

	PDB_Entry EMPTY;

	//molecules
	Mols Conformers;
	//atom list indexed by conformer ID + atom number
	map<int, map<int, map<string, PDB_Entry> > > ATOMS; 
	//atom list indexed by conformer ID + residue number + atom name

	map<int, string> residList, residListOne; //resdiue lists with three-letter-aa and one-letter-aa name
	int r1, rN;

	vector<int> acceptorList; // Acceptor List for H-bond indexed by atom index
	map<int, int> donorList; // Donor List for H-bond indexed by atom indeice of itself and its connected heavy atom
	map<int, map<string, float> > HBDistList; // HBond distance indexed by residue number and atom name

	PDB();
	PDB(const string& fileName);

	//amino acid name convertion
	string getThreeAAName(char a);
	string getOneAAName(const string& a);

	//load pdb coordinates
	void loadPDB(const string &fileName);
	void loadPDB_Entry(const string &str, PDB_Entry &entry);
	string getField(const string &str, int index);

	//get atom
	PDB_Entry getEntry(int conformerID, int rNum, const string &aName); //get atom by residue number and atom name
	PDB_Entry getEntry(int conformerID, int aNum); //get atom by atom number

	//get angles and distances
	float getBondAngle(Vec3 A, Vec3 B, Vec3 C);
	float getBondAngle(PDB_Entry a, PDB_Entry b, PDB_Entry c);
	float getDihedralAngle(PDB_Entry a, PDB_Entry b, PDB_Entry c, PDB_Entry d);
	float getPhi(int conformerID, int resNum);
	float getPsi(int conformerID, int resNum);
	float getOmega(int conformerID, int resNum);
	float getChi1(int conformerID, int resNum);
	float getChi2(int conformerID, int resNum);
	float getDist(Vec3 A, Vec3 B);
	float getDist(PDB_Entry A, PDB_Entry B);

	bool isSSBonded(int conformerID, int resNum);

	//get ring current shifts
	void initOrbitalShift();
	float getOrbitalShift(int conformerID, int resNum, const string &aName);
	void calcPlane(RingData *ringP);

	//get H-bond information
	void initHBond();
	float getHBondDist(PDB_Entry D);
	float getHBondDist(int resNum, string atomName);
	bool isAcceptor(PDB_Entry A); //check if an atom is an Acceptor for a H-Hond
	PDB_Entry isDonor(PDB_Entry D); //check if an aotm is a Donor for a H-hond, and return the connected heavy atom

	//string functions
	static string simplifyWhiteSpace(const string &str);
	static bool isSpace( const char &c );

	int getVectors(DoubleD2 &vec, int rr1, int rrN, vector<string> &atomNames, vector<int> &segList);
	//get coordinates for a given frament, return the number of selected atoms
	int getVectors(DoubleD2 &vec, vector<string> &atomNames, vector<int> &segList);

	double r_sign(double x, double y);
	int svd_(DoubleD2 &a, DoubleD1 &w, DoubleD2 &v);
	double rmsd(DoubleD2 &a, DoubleD2 &b, int n);
	double rmsd(vector<Vec3> &a, vector<Vec3> &b);

};

#endif