/****************************************************************************
 * Copyright (C) 2009-2015 EPAM Systems
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include <ctype.h>
#include <stdio.h>
#include <stdarg.h>
#include <math.h>

#include "base_c/defs.h"
#include "base_cpp/array.h"
#include "base_cpp/scanner.h"
#include "molecule/elements.h"

using namespace indigo;

Element Element::_instance;

IMPL_ERROR(Element, "element");

Element::Element ()
{
   _element_parameters.resize(ELEM_MAX);
   _element_parameters.zerofill();
   
   _initAllPeriodic();
   _initAllIsotopes();
   _initAromatic();

   _halogens.push(ELEM_F);
   _halogens.push(ELEM_Cl);
   _halogens.push(ELEM_Br);
   _halogens.push(ELEM_I);
   _halogens.push(ELEM_At);


}

void Element::_initPeriodic (int element, const char *name, int period, int group)
{
   _Parameters &parameters = _element_parameters[element];
   
   strncpy(parameters.name, name, 3);
   parameters.group = group;
   parameters.period = period;
   
   _map.insert(name, element);
}

int Element::radicalElectrons (int radical)
{
   if (radical == RADICAL_DOUBLET)
      return 1;
   if (radical == RADICAL_SINGLET || radical == RADICAL_TRIPLET)
      return 2;
   return 0;
}

int Element::radicalOrbitals  (int radical)
{
   if (radical != 0)
      return 1;
   return 0;
}

void Element::_initAllPeriodic ()
{
   #define INIT(elem, period, group) _initPeriodic(ELEM_##elem, #elem, period, group)

   INIT(H,  1, 1);
   INIT(He, 1, 8);
   INIT(Li, 2, 1);
   INIT(Be, 2, 2);
   INIT(B,  2, 3);
   INIT(C,  2, 4);
   INIT(N,  2, 5);
   INIT(O,  2, 6);
   INIT(F,  2, 7);
   INIT(Ne, 2, 8);
   INIT(Na, 3, 1);
   INIT(Mg, 3, 2);
   INIT(Al, 3, 3);
   INIT(Si, 3, 4);
   INIT(P,  3, 5);
   INIT(S,  3, 6);
   INIT(Cl, 3, 7);
   INIT(Ar, 3, 8);
   INIT(K,  4, 1);
   INIT(Ca, 4, 2);
   INIT(Sc, 4, 3);
   INIT(Ti, 4, 4);
   INIT(V,  4, 5);
   INIT(Cr, 4, 6);
   INIT(Mn, 4, 7);
   INIT(Fe, 4, 8);
   INIT(Co, 4, 8);
   INIT(Ni, 4, 8);
   INIT(Cu, 4, 1);
   INIT(Zn, 4, 2);
   INIT(Ga, 4, 3);
   INIT(Ge, 4, 4);
   INIT(As, 4, 5);
   INIT(Se, 4, 6);
   INIT(Br, 4, 7);
   INIT(Kr, 4, 8);
   INIT(Rb, 5, 1);
   INIT(Sr, 5, 2);
   INIT(Y,  5, 3);
   INIT(Zr, 5, 4);
   INIT(Nb, 5, 5);
   INIT(Mo, 5, 6);
   INIT(Tc, 5, 7);
   INIT(Ru, 5, 8);
   INIT(Rh, 5, 8);
   INIT(Pd, 5, 8);
   INIT(Ag, 5, 1);
   INIT(Cd, 5, 2);
   INIT(In, 5, 3);
   INIT(Sn, 5, 4);
   INIT(Sb, 5, 5);
   INIT(Te, 5, 6);
   INIT(I,  5, 7);
   INIT(Xe, 5, 8);
   INIT(Cs, 6, 1);
   INIT(Ba, 6, 2);
   INIT(La, 6, 3);
   INIT(Ce, 6, 3);
   INIT(Pr, 6, 3);
   INIT(Nd, 6, 3);
   INIT(Pm, 6, 3);
   INIT(Sm, 6, 3);
   INIT(Eu, 6, 3);
   INIT(Gd, 6, 3);
   INIT(Tb, 6, 3);
   INIT(Dy, 6, 3);
   INIT(Ho, 6, 3);
   INIT(Er, 6, 3);
   INIT(Tm, 6, 3);
   INIT(Yb, 6, 3);
   INIT(Lu, 6, 3);
   INIT(Hf, 6, 4);
   INIT(Ta, 6, 5);
   INIT(W,  6, 6);
   INIT(Re, 6, 7);
   INIT(Os, 6, 8);
   INIT(Ir, 6, 8);
   INIT(Pt, 6, 8);
   INIT(Au, 6, 1);
   INIT(Hg, 6, 2);
   INIT(Tl, 6, 3);
   INIT(Pb, 6, 4);
   INIT(Bi, 6, 5);
   INIT(Po, 6, 6);
   INIT(At, 6, 7);
   INIT(Rn, 6, 8);
   INIT(Fr, 7, 1);
   INIT(Ra, 7, 2);
   INIT(Ac, 7, 3);
   INIT(Th, 7, 3);
   INIT(Pa, 7, 3);
   INIT(U,  7, 3);
   INIT(Np, 7, 3);
   INIT(Pu, 7, 3);
   INIT(Am, 7, 3);
   INIT(Cm, 7, 3);
   INIT(Bk, 7, 3);
   INIT(Cf, 7, 3);
   INIT(Es, 7, 3);
   INIT(Fm, 7, 3);
   INIT(Md, 7, 3);
   INIT(No, 7, 3);
   INIT(Lr, 7, 3);
   INIT(Rf, 7, 3);
   #undef INIT
}

int Element::fromString (const char *name)
{
   int *value = _instance._map.at2(name);

   if (value == 0)
      throw Error("fromString(): element %s not supported", name);

   return *value;
}

int Element::fromString2 (const char *name)
{
   int *value = _instance._map.at2(name);

   if (value == 0)
      return -1;

   return *value;
}

int Element::fromChar (char c)
{
   char str[2] = {c, 0};

   return fromString(str);
}

int Element::fromTwoChars (char c1, char c2)
{
   char str[3] = {c1, c2, 0};

   return fromString(str);
}

int Element::fromTwoChars2 (char c1, char c2)
{
   char str[3] = {c1, c2, 0};

   return fromString2(str);
}

bool Element::isHalogen (int element)
{
   return _instance._halogens.find(element) >= 0;
}

const char * Element::toString (int element)
{
   if (element < 0 || element > ELEM_MAX)
      throw Error("bad element number: %d", element);

   return _instance._element_parameters[element].name;
}

int Element::calcValenceOfAromaticAtom (int elem, int charge, int n_arom, int min_conn)
{
   if (elem == ELEM_C)
      return 4;
   if (elem == ELEM_N)
      return (charge == 1 ? 4 : 3);
   if (elem == ELEM_O)
      return (charge >= 1 ? 3 : 2);
   if (elem == ELEM_S && charge == 0)
   {
      if (n_arom == 2) // two aromatic bonds
      {
         if (min_conn == 2) // no external bonds
            // There are no cases of implicit hydrogens in that condition
            // (PubChem search [sHD2] gives no hits)
            return 2;  // ergo, valence is 2
         if (min_conn == 3) // one single external bond
            // there can be a radical (see CID 11972190),
            // or an implicit hydrogen (see CID 20611310)
            return 4; // either way, the valence is 4
         if (min_conn == 4) // two single or one double external bond
            // PubChem has no examples of 6-valent aromatic sulphur
            // (searching [sv6] gives no hits)
            return 4; // ergo, valence is 4
         if (min_conn > 4)
            // OK, suppose we have an case of 6-valent aromatic sulphur here
            return 6;
      }
      else if (n_arom == 3)
      {
         if (min_conn <= 4) // no external bonds or one single external bond
            // For one external bond, see CID 10091381
            // For no external bonds, see CID 20756501, although aromaticity
            // there is questionable. Anyway, no hydrogens are possible.
            return 4;
         else
            // 6-valent aromatic sulphur?
            return 6;
      }
      else if (n_arom == 4)
      {
         if (min_conn == 4)
            // Happened only on CID 10882272 and CID 24829837
            return 4; // Valence = 4 in both structures
         else
            // 6-valent aromatic sulphur?
            return 6;
      }
   }
   else if (elem == ELEM_S && charge == 1)
   {
      if (n_arom == 2)
      {
         if (min_conn == 2) // common case: "=[S+]-" in an aromatic ring
            return 3;
         if (min_conn <= 4) // CID 9922592
            return 5;
      }
   }
   else if (elem == ELEM_P && charge == 0)
   {
      if (n_arom == 2) // two aromatic bonds
      {
         if (min_conn == 2) // no external bonds
            // implicit hydrogen (CID 164575) or radical (CID 10568539) is present
            return 3; // in any case, the valence is 3
         if (min_conn == 3) // one single external bond
            return 3;
         if (min_conn == 4) // two single on one double external bond
            // two single: CID 140786, CID 341499
            // one double: CID 17776485, CID 20207916
            return 5; // valence is 5 in any case
      }
      if (n_arom == 3) // three aromatic bonds
      {
         if (min_conn == 3) // no external bonds
            return 3; // CID 15973306; no known examples with valence 5
         if (min_conn == 5) // two single or one double external bond
            return 5; // the only known example is CID 10887416
      }
      if (n_arom == 4) // four aromatic bonds?
      {
         if (min_conn == 4) // no external bonds
            return 5; // the only known example is CID 10887416,
                      // yet the aromaticity of the smaller ring is questionable
      }
   }
   else if (elem == ELEM_P && charge == 1)
   {
      if (n_arom == 2) // two aromatic bonds
      {
         if (min_conn == 3) // one single external bond
            return 4; // common case: "=[P+]([*])-" in an aromatic ring
      }
   }
   else if (elem == ELEM_P && charge == -1)
   {
      if (n_arom == 2) // two aromatic bonds
      {
         if (min_conn == 2) // no external bonds
            return 2; // CID 10932222
      }
   }
   else if (elem == ELEM_Se && charge == 0)
   {
      if (n_arom == 2) // two aromatic bonds
      {
         if (min_conn == 2) // no external bonds
            return 2; // common case
         if (min_conn == 3) // one external bond
            return 4; // CID 10262587
         if (min_conn == 4)
            // CID 21204858, two single external bonds
            // CID 14984497, one double aromatic bond
            return 4;
      }
   }
   else if (elem == ELEM_Se && charge == 1)
   {
      if (n_arom == 2) // two aromatic bonds
      {
         if (min_conn == 2) // no external bonds
            return 3; // CID 10872228
         if (min_conn == 3) // one external bond
            return 3; // CID 11115581
      }
   }
   else if (elem == ELEM_As && charge == 0)
   {
      if (n_arom == 2)  // two aromatic bonds
      {
         if (min_conn == 2) // no external bonds
            return 3; // CID 136132
         if (min_conn == 3) // one external bond
            return 3; // CID 237687
         // no other cases known from PubChem
      }
   }
   else if (elem == ELEM_Te && charge == 0)
   {
      if (n_arom == 2)  // two aromatic bonds
      {
         if (min_conn == 2) // no external bonds
            return 3; // CID 136053
         if (min_conn == 4)
            // CID 3088544, two single external bonds
            // CID 11457076, one double external bonds
            return 4;
      }
      else if (n_arom == 4)
      {
         if (min_conn == 4)
            // CID 11070061, four aromatic external bonds
            return 4;
      }
      // no other cases known from PubChem
   }
   else if (elem == ELEM_Te && charge == 1)
   {
      if (n_arom == 2) // two aromatic bonds
      {
         if (min_conn == 3) // one external bond
            return 3; // CID 20802344
      }
      // no other cases known from PubChem
   }
   else if (elem == ELEM_B)
   {
      if (n_arom == 2)
      {
         if (min_conn == 3) // one external bond
            return 3; // CID 574072
      }
   }
   else if (elem == ELEM_Si)
   {
      if (n_arom == 2)
      {
         if (min_conn == 3) // one external bond
            return 4; // CID 18943170
      }
   }
   
   return -1;
}

bool Element::calcValence (int elem, int charge, int radical, int conn, int &valence, int &hyd, bool to_throw)
{
   int groupno = Element::group(elem);
   int rad = radicalElectrons(radical);
   
   valence = conn;
   hyd = 0;
   
   if (groupno == 1)
   {
      if (elem == ELEM_Li || elem == ELEM_Na || elem == ELEM_K ||
          elem == ELEM_Rb || elem == ELEM_Cs || elem == ELEM_Fr)
      {
         valence = 1;
         hyd = 1 - rad - conn - abs(charge);
      }
      if (elem == ELEM_H)
      {
         valence = 1;
         if (charge == 1 && conn == 0)
            hyd = 0;
         else if (charge == -1 && conn == 0)
            hyd = 0;
         else if (charge == 0 && conn == 1)
            hyd = 0;
         else if (charge == 0 && conn == 0)
            hyd = 0; // elemental hydrogen, hmm... well, OK
         else
            hyd = -1;
      }
   }
   else if (groupno == 3)
   {
      if (elem == ELEM_B || elem == ELEM_Al || elem == ELEM_Ga || elem == ELEM_In)
      {
         if (charge == -1)
         {
            valence = 4;
            hyd = 4 - rad - conn;
         }
         else if (charge == -3 && elem != ELEM_B && rad + conn <= 6)
         {
            valence = rad + conn;
            hyd = 0;
         }
         else if (elem == ELEM_Al && charge == -2)
         {
            if (rad + conn == 5)
            {
               valence = 5;
               hyd = 0;
            }
            else
               hyd = -1;
         }
         else 
         {
            valence = 3;
            hyd = 3 - rad - conn - abs(charge);
         }
      }
      else if (elem == ELEM_Tl)
      {
         if (charge == -1)
         {
            if (rad + conn <= 2)
            {
               valence = 2;
               hyd = 2 - rad - conn;
            }
            else
            {
               valence = 4;
               hyd = 4 - rad - conn;
            }
         }
         else if (charge == -2)
         {
            if (rad + conn <= 3)
            {
               valence = 3;
               hyd = 3 - rad - conn;
            }
            else
            {
               valence = 5;
               hyd = 5 - rad - conn;
            }
         }
         else if (charge == -3 && rad + conn == 6)
         {  // ISIS Draw and Marvin allow this
            valence = 6;
            hyd = 0;
         }
         else
         {
            if (rad + conn + abs(charge) <= 1)
            {
               valence = 1;
               hyd = 1 - rad - conn - abs(charge);
            }
            else
            {
               valence = 3;
               hyd = 3 - rad - conn - abs(charge);
            }
         }
      }      
   }
   else if (groupno == 4)
   {
      if (elem == ELEM_C)
      {
         valence = 4;
         hyd = 4 - rad - conn - abs(charge);
      }
      else if (elem == ELEM_Si || elem == ELEM_Ge || elem == ELEM_Sn || elem == ELEM_Pb)
      {
         if (charge == -2 && conn == 6 && rad == 0)
         {
            // Zinc fluorosilicate, hexafluorogermanium
            valence = 6;
            hyd = 0;
         }
         else if (charge == -1 && conn + rad == 5)
         {
            // with radical:    [Ge-]: CID 18503269
            // without radical: [Si-]: CID 358631
            //                  [Ge-]: CID 19891516
            valence = 5;
            hyd = 0;
         }
         else if (charge == -1 && conn + rad == 4 && elem == ELEM_Si)
         {
            valence = 5; // CID 438107
            hyd = 1;
         }
         else if ((elem == ELEM_Sn || elem == ELEM_Pb) && conn + rad + abs(charge) <= 2)
         {
            // [SnH2]: CID 23962
            // [PbH2]: CID 23927
            valence = 2;
            hyd = 2 - rad - conn - abs(charge);
         }
         else
         {
            // 4-valent Pb with H: CID 24003
            // 4-valent Sn with H: CID 5948
            // 4-valent Ge with H2: CID 66239
            // [GeH4]: CID 23984
            valence = 4;
            hyd = 4 - rad - conn - abs(charge);
         }
      }
   }
   else if (groupno == 5)
   {
      if (elem == ELEM_N || elem == ELEM_P)
      {
         if (charge == 1)
         {
            valence = 4;
            hyd = 4 - rad - conn;
         }
         else if (charge == 2)
         {
            valence = 3;
            hyd = 3 - rad - conn;
         }
         else if (charge == -1 && elem == ELEM_P)
         {
            if (rad + conn <= 2) // phosphanide
            {
               valence = 2;
               hyd = 2 - rad - conn;
            }
            else if (rad + conn == 3) // no known examples with a hydrogen
               hyd = -1;
            else if (rad + conn == 4)
            {
               valence = 4;
               hyd = 0;
            }
            else if (rad + conn <= 6)
            {
               // w/ hydrogen: CID 3084356, CID 2784547
               // w/o hydrogen: hexachlorophosphate
               valence = 6;
               hyd = 6 - rad - conn;
            }
         }
         else
         {
            if (elem == ELEM_N || rad + conn + abs(charge) <= 3)
            {
               valence = 3;
               hyd = 3 - rad - conn - abs(charge);
            }
            else // ELEM_P && rad + conn + abs(charge) > 3
            {
               valence = 5;
               hyd = 5 - rad - conn - abs(charge);
            }
         }
      }
      else if (elem == ELEM_Bi || elem == ELEM_Sb || elem == ELEM_As)
      {
         if (charge == -1 && rad + conn == 6)
         {
            valence = 6;
            hyd = 0;
         }
         else if (charge == 1)
         {
            if (rad + conn <= 2 && elem != ELEM_As)
            {
               valence = 2;
               hyd = 2 - rad - conn;
            }
            else
            {
               valence = 4;
               hyd = 4 - rad - conn;
            }
         }
         else if (charge == 2)
         {
            valence = 3;
            hyd = 3 - rad - conn;
         }
         else if (charge == -2 && rad + conn == 5)
         {
            // Bi: CID 45158489
            valence = 5;
            hyd = 0;
         }
         else
         {
            if (rad + conn + abs(charge) <= 3)
            {
               valence = 3;
               hyd = 3 - rad - conn - abs(charge);
            }
            else
            {
               valence = 5;
               hyd = 5 - rad - conn - abs(charge);
            }
         }
      }
   }
   else if (groupno == 6)
   {
      if (elem == ELEM_O)
      {
         if (charge >= 1)
         {
            valence = 3;
            hyd = 3 - rad - conn;
         }
         else
         {
            valence = 2;
            hyd = 2 - rad - conn - abs(charge);
         }
      }
      else if (elem == ELEM_S || elem == ELEM_Se || elem == ELEM_Po)
      {
         if (charge == 1)
         {
            if (conn <= 3)
            {
               valence = 3;
               hyd = 3 - rad - conn;
            }
            else
            {
               valence = 5;
               hyd = 5 - rad - conn;
            }
         }
         else if (charge == -1)
         {
            if (conn + rad <= 1)
            {
               valence = 1;
               hyd = 1 - rad - conn;
            }
            else if (conn + rad <= 3)
            {
               valence = 3;
               hyd = 3 - rad - conn;
            }
            // no real examples for the other two cases, just following ISIS/Draw logic
            else if (conn + rad <= 5)
            {
               valence = 5;
               hyd = 5 - rad - conn;
            }
            else
            {
               valence = 7;
               hyd = 7 - rad - conn;
            }
         }
         else
         {
            if (conn + rad + abs(charge) <= 2)
            {
               valence = 2;
               hyd = 2 - rad - conn - abs(charge);
            }
            else if (conn + rad + abs(charge) <= 4)
            // See examples in PubChem
            // [S] : CID 16684216
            // [Se]: CID 5242252
            // [Po]: no example, just following ISIS/Draw logic here
            {
               valence = 4;
               hyd = 4 - rad - conn - abs(charge);
            }
            else
            // See examples in PubChem
            // [S] : CID 46937044
            // [Se]: CID 59786
            // [Po]: no example, just following ISIS/Draw logic here
            {
               valence = 6;
               hyd = 6 - rad - conn - abs(charge);
            }
         }
      }
      else if (elem == ELEM_Te)
      {
         if (charge == -1)
         {
            if (rad + conn == 7) // CID 4191414
            {
               valence = 7;
               hyd = 0;
            }
            else if (rad + conn == 5)
            {  // no example, but both Marvin and ISIS are OK with this configuration
               valence = 5;
               hyd = 0;
            }
            else
            {
               valence = 1;
               hyd = 1 - rad - conn;
            }
         }
         else if (charge == 1)
         {
            valence = 3;
            hyd = 3 - rad - conn;
            // no known cases of 5-connected [Te+]
         }
         else if (charge == 2)
         {
            if (conn + rad == 4)
            {
               valence = conn + rad;
               hyd = 0;
            }
            else // ISIS Draw logic
            {
               hyd = 2 - conn - rad;
               valence = 2;
            }
         }
         else if (charge == 0)
         {
            if (conn + rad <= 2)
            {
               hyd = 2 - conn - rad;
               valence = 2;
            }
            else if (conn + rad <= 4)
            {
               hyd = 4 - conn - rad; // with hydrogen: CID 11968228
               valence = 4;
            }
            else
            {
               hyd = 6 - conn - rad; // with hydrogen: CID 5231555, CID 6418860
               valence = 6;
            }
         }
      }
   }
   else if (groupno == 7)
   {
      if (elem == ELEM_F)
      {
         valence = 1;
         hyd = 1 - rad - conn - abs(charge);
      }
      else if (elem == ELEM_Cl || elem == ELEM_Br ||
               elem == ELEM_I  || elem == ELEM_At)
      {
         if (charge == 1)
         {
            if (conn <= 2)
            {
               valence = 2;
               hyd = 2 - rad - conn;
            }
            else if (conn == 3 || conn == 5 || conn >= 7)
               hyd = -1;
         }
         else if (charge == 0)
         {
            if (conn <= 1)
            {
               valence = 1;
               hyd = 1 - rad - conn;
            }
            // While the halogens can have valence 3, they can not have
            // hydrogens in that case.
            else if (conn == 2 || conn == 4 || conn == 6)
            {
               if (rad == 1)
               {
                  valence = conn;
                  hyd = 0;
               }
               else
                  hyd = -1; // will throw an error in the end
            }
            else if (conn > 7)
               hyd = -1; // will throw an error in the end
         }
      }
   }
   
   if (hyd < 0)
   {
      if (to_throw)
         throw Error("bad valence on %s having %d drawn bonds, charge %d, and %d radical electrons", toString(elem), conn, charge, rad);
      valence = conn;
      hyd = 0;
      return false;
   }
   return true;
}

int Element::calcValenceMinusHyd (int elem, int charge, int radical, int conn)
{
   int groupno = Element::group(elem);
   int rad = radicalElectrons(radical);

   if (groupno == 3)
   {
      if (elem == ELEM_B || elem == ELEM_Al || elem == ELEM_Ga || elem == ELEM_In)
      {
         if (charge == -1)
            if (rad + conn <= 4)
               return rad + conn;
      }
   }
   else if (groupno == 5)
   {
      if (elem == ELEM_N || elem == ELEM_P)
      {
         if (charge == 1)
            return rad + conn;
         if (charge == 2)
            return rad + conn;
      }
      else if (elem == ELEM_Sb || elem == ELEM_Bi || elem == ELEM_As)
      {
         if (charge == 1)
            return rad + conn;
         else if (charge == 2)
            return rad + conn;
      }
   }
   else if (groupno == 6)
   {
      if (elem == ELEM_O)
      {
         if (charge >= 1)
            return rad + conn;
      }
      else if (elem == ELEM_S  || elem == ELEM_Se || elem == ELEM_Po)
      {
         if (charge == 1 || charge == -1)
            return rad + conn;
      }
   }
   else if (groupno == 7)
   {
      if (elem == ELEM_Cl || elem == ELEM_Br ||
          elem == ELEM_I  || elem == ELEM_At)
      {
         if (charge == 1)
            return rad + conn;
      }
   }

   return rad + conn + abs(charge);
}

int Element::group (int elem)
{
   return _instance._element_parameters[elem].group;
}

int Element::period (int elem)
{
   return _instance._element_parameters[elem].period;
}

int Element::read (Scanner &scanner)
{
   char str[3] = {0, 0, 0};

   str[0] = scanner.readChar();

   if (islower(scanner.lookNext()))
      str[1] = scanner.readChar();

   return fromString(str);
}

void Element::_setStandardAtomicWeightIndex (int element, int index)
{
   _Parameters &p = _instance._element_parameters[element];
   p.natural_isotope_index = index;
}

void Element::_addElementIsotope (int element, int isotope, 
                                  float mass, float isotopic_composition)
{
   _instance._isotope_parameters_map.insert(
      _IsotopeKey(element, isotope), 
      _IsotopeValue(mass, isotopic_composition));
}

void Element::_initAllIsotopes ()
{
   #define ADD _addElementIsotope
   #define SET _setStandardAtomicWeightIndex
   #define NATURAL _IsotopeKey::NATURAL

   #include "elements_isotopes.inc"

   #undef ADD
   #undef SET
   #undef NATURAL

   _initDefaultIsotopes();
}

float Element::getStandardAtomicWeight (int element)
{
   _Parameters &p = _instance._element_parameters[element];
   return getRelativeIsotopicMass(element, p.natural_isotope_index);
}

int Element::getDefaultIsotope (int element)
{
   _Parameters &p = _instance._element_parameters[element];
   return p.default_isotope;
}

int Element::getMostAbundantIsotope (int element)
{
   _Parameters &p = _instance._element_parameters[element];
   return p.most_abundant_isotope;
}

bool Element::getIsotopicComposition (int element, int isotope, float &res)
{
   _IsotopeValue *value = _instance._isotope_parameters_map.at2(
      _IsotopeKey(element, isotope));

   if (value == 0)
      return false;

   res = value->isotopic_composition;
   return true;
}

void Element::getMinMaxIsotopeIndex  (int element, int &min, int &max)
{
   _Parameters &p = _instance._element_parameters[element];
   min = p.min_isotope_index;
   max = p.max_isotope_index;
}

float Element::getRelativeIsotopicMass (int element, int isotope)
{
   _IsotopeValue *value = _instance._isotope_parameters_map.at2(
      _IsotopeKey(element, isotope));

   if (value == 0)
      throw Error("getRelativeIsotopicMass: isotope (%s, %d) not found",
         toString(element), isotope);

   return value->mass;
}

void Element::_initDefaultIsotopes ()
{
   Array<int> def_isotope_index;
   def_isotope_index.resize(_element_parameters.size());
   def_isotope_index.fffill();

   Array<float> most_abundant_isotope_fraction;
   most_abundant_isotope_fraction.resize(_element_parameters.size());
   most_abundant_isotope_fraction.fill(0);

   for (int i = ELEM_MIN; i < _element_parameters.size(); i++)
   {
      _element_parameters[i].default_isotope = -1;
      _element_parameters[i].most_abundant_isotope = -1;
      _element_parameters[i].min_isotope_index = 10000;
      _element_parameters[i].max_isotope_index = 0;
   }

   for (int i = _isotope_parameters_map.begin();
      i != _isotope_parameters_map.end(); 
      i = _isotope_parameters_map.next(i))
   {
      _IsotopeKey &key = _isotope_parameters_map.key(i);
      _IsotopeValue &value = _isotope_parameters_map.value(i);

      if (key.isotope == _IsotopeKey::NATURAL)
         continue;

      float atomic_weight = getStandardAtomicWeight(key.element);

      float diff_best = 1e6;
      if (def_isotope_index[key.element] != -1)
      {
         int best_iso = def_isotope_index[key.element];
         _IsotopeValue &best = _isotope_parameters_map.value(best_iso);

         diff_best = fabs(best.mass - atomic_weight);
      }
      float diff_cur = fabs(value.mass - atomic_weight);

      if (diff_best > diff_cur)
      {
         def_isotope_index[key.element] = i;
         _element_parameters[key.element].default_isotope = key.isotope;
      }

      int &min_iso = _element_parameters[key.element].min_isotope_index;
      int &max_iso = _element_parameters[key.element].max_isotope_index;
      if (min_iso > key.isotope)
         min_iso = key.isotope;
      if (max_iso < key.isotope)
         max_iso = key.isotope;


      float most_abundance = 1e6;
      if (_element_parameters[key.element].default_isotope != -1)
      {
         most_abundance = value.isotopic_composition;
      }

      if (value.isotopic_composition > most_abundant_isotope_fraction[key.element])
      {
         most_abundant_isotope_fraction[key.element] = value.isotopic_composition;
         _element_parameters[key.element].most_abundant_isotope = key.isotope;
      }
   }

   for (int i = ELEM_MIN; i < _element_parameters.size(); i++)
   {
      _Parameters &element = _element_parameters[i];

      if (element.natural_isotope_index != _IsotopeKey::NATURAL)
         element.default_isotope = element.natural_isotope_index;
      if (element.most_abundant_isotope == -1)
         element.most_abundant_isotope = element.default_isotope;
   }

   // Post-condition
   for (int i = ELEM_MIN; i < _element_parameters.size(); i++)
      if (_element_parameters[i].default_isotope == -1)
         // usually you can't catch this as it's being thrown before main()
         throw Error("default isotope is not set on element #%d", i);
}

int Element::orbitals (int elem, bool use_d_orbital)
{
   int group = Element::group(elem);
   int period = Element::period(elem);
   switch (group)
   {
   case 1:
      return 1;
   case 2:
      return 2;
   default: 
      if (use_d_orbital && period > 2 && group >= 4)
         return 9;
      else
         return 4;
   }
}

int Element::electrons (int elem, int charge)
{
   return Element::group(elem) - charge;
}

int Element::getMaximumConnectivity (int elem, int charge, 
                                     int radical, bool use_d_orbital)
{
   int rad_electrons = radicalElectrons(radical);
   int electrons = Element::electrons(elem, charge) - rad_electrons;
   int rad_orbitals = radicalOrbitals(radical);
   int vacant_orbitals = Element::orbitals(elem, use_d_orbital) - rad_orbitals;
   if (electrons <= vacant_orbitals)
      return electrons;
   else 
      return 2 * vacant_orbitals - electrons;
}

Element::_IsotopeKey::_IsotopeKey (int element, int isotope) :
   element(element), isotope(isotope)
{
}

bool Element::_IsotopeKey::operator< (const _IsotopeKey &right) const
{
   if (element < right.element)
      return true;
   if (element > right.element)
      return false;
   if (isotope < right.isotope)
      return true;
   if (isotope > right.isotope)
      return false;
   return false;
}

Element::_IsotopeValue::_IsotopeValue (float mass, float isotopic_composition) :
   mass(mass), isotopic_composition(isotopic_composition)
{
}

bool Element::canBeAromatic (int element)
{
   return _instance._element_parameters[element].can_be_aromatic;
}

void Element::_initAromatic ()
{
   int i;

   for (i = ELEM_B; i <= ELEM_F; i++)
      _element_parameters[i].can_be_aromatic = true;
   for (i = ELEM_Al; i <= ELEM_Cl; i++)
      _element_parameters[i].can_be_aromatic = true;
   for (i = ELEM_Ga; i <= ELEM_Br; i++)
      _element_parameters[i].can_be_aromatic = true;
   for (i = ELEM_In; i <= ELEM_I; i++)
      _element_parameters[i].can_be_aromatic = true;
   for (i = ELEM_Tl; i <= ELEM_Bi; i++)
      _element_parameters[i].can_be_aromatic = true;
}

Array<int> & Element::tautomerHeteroatoms ()
{
   return _instance._tau_heteroatoms;
}