/*-------------------------------------------------------------------------
| The software accompanies the paper                                      |
|                                                                         |
| Classes and Objects of Chemical Thermodynamics in Object-Oriented       |
| Programming. 2. A Class of Chemical Species                             |
|                                                                         |
| E.B. Rudnyi                                                             |
| E-mail: rudnyi@comp.chem.msu.su                                         |
| Homepages: http://www.chem.msu.su/~rudnyi/welcome.html                  |
|                                                                         |
| presented at Second Electronic Computational Chemistry Conference,      |
| November 1995, http://hackberry.chem.niu.edu/ECCC2/                     |
|                                                                         |
| E.B. Rudnyi, Copyright (c) 1995                                         |
|                                                                         |
| Permission to use, copy, modify, distribute and sell this software and  |
| its documentation for any purpose is hereby granted without fee,        |
| provided that the above copyright notice with the name of the paper and |
| the conference appear in all copies and that both that copyright notice |
| and this permission notice appear in supporting documentation.          |
| E.B. Rudnyi makes no representations about the suitability of this      |
| software for any purpose. It is provided "as is" without expressed or   |
| implied warranty.                                                       |
--------------------------------------------------------------------------*/

#include <iostream.h>
#include <strstrea.h>
#include <ctype.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "species.h"

char elem::el_name[NELEMENTS][2] =
{ 0 , 0 ,'e', 0 ,'O', 0 ,'H', 0 ,'D', 0 ,'T', 0 ,'F', 0 ,'C','l','B','r','I', 0 ,
 'A','t','H','e','N','e','A','r','K','r','X','e','R','n','S', 0 ,'S','e','T','e',
 'P','o','N', 0 ,'P', 0 ,'A','s','S','b','B','i','C', 0 ,'S','i','G','e','S','n',
 'P','b','B', 0 ,'A','l','G','a','I','n','T','l','Z','n','C','d','H','g','C','u',
 'A','g','A','u','F','e','C','o','N','i','R','u','R','h','P','d','O','s','I','r',
 'P','t','M','n','T','c','R','e','C','r','M','o','W', 0 ,'V', 0 ,'N','b','T','a',
 'T','i','Z','r','H','f','S','c','Y', 0 ,'L','a','C','e','P','r','N','d','P','m',
 'S','m','E','u','G','d','T','b','D','y','H','o','E','r','T','m','Y','b','L','u',
 'A','c','T','h','P','a','U', 0 ,'N','p','P','u','A','m','C','m','B','k','C','f',
 'E','s','F','m','M','d','N','o','B','e','M','g','C','a','S','r','B','a','R','a',
 'L','i','N','a','K', 0 ,'R','b','C','s','F','r'};

float elem::el_mass[NELEMENTS] = {
// NEL
0,
// e,      O ,  H ,   D ,    T ,    F ,    Cl,    Br,     I ,   At,
5.485e-4,15.9994,1.0079,2.014,3.016,18.9984,35.453,79.904,126.904,210.,
// He,   Ne,     Ar,   Kr,   Xe,   Rn,   S ,  Se,  Te,   Po,
4.0026,20.179,39.948,83.80,131.30,222.,32.06,78.9,127.6,209.,
// N ,     P ,      As,     Sb,    Bi,     C ,    Si,    Ge,   Sn,    Pb,
14.0067,30.97376,74.9216,121.75,208.9804,12.011,28.0855,72.59,118.69,207.2,
// B ,   Al,     Ga,  In,     Tl,   Zn,   Cd,     Hg,    Cu,     Ag,
10.81,26.98154,69.72,114.82,204.37,65.38,112.41,200.59,63.546,107.868,
// Au,     Fe,     Co,   Ni,    Ru,     Rh,    Pd,   Os,    Ir,    Pt,
196.9665,55.847,58.9332,58.70,101.07,102.9055,106.4,190.2,192.22,195.09,
// Mn,  Tc,   Re,     Cr,   Mo,   W ,     V ,    Nb,     Ta,     Ti,
54.938,97.,186.207,51.996,95.94,183.85,50.9414,92.9064,180.9479,47.90,
// Zr,  Hf,    Sc,     Y ,      La,     Ce,     Pr,    Nd,   Pm,   Sm,
91.22,178.49,44.9559,88.9059,138.9055,140.12,140.9077,144.24,145.,150.4,
// Eu,    Gd,     Tb,     Dy,     Ho,     Er,     Tm,     Yb,   Lu,   Ac,
151.96, 157.25,158.9254,162.50,164.9304,167.26,168.9342,173.04,174.97,227.,
// Th,      Pa,       U ,   Np,     Pu, Am,  Cm,   Bk, Cf,  Es,
232.0381,231.0359,238.029,237.0482,244.,243.,247.,247.,251.,254.,
// Fm, Md, No,    Be,     Mg,   Ca,   Sr,   Ba,     Ra,     Li,
257.,258.,259.,9.01218,24.305,40.08,87.62,137.33,226.0254,6.941,
// Na,      K ,     Rb,     Cs,   Fr
22.98977,39.0983,85.4678,132.9054,223.};

const char* elem::anal(const char *in) {
  el = NEL;
  if (*in == '\0') return in;
  while (isspace(*in)) in++;
  if (!isupper(*in) && *in != 'e') return in;
  int find = 0;
  for (int i = 1; i < NELEMENTS; i++) {
    if (in[0] == el_name[i][0]) {
      if (el_name[i][1]) {
        if (el_name[i][1] == in[1]) {
          el = (elements)i;
          return in += 2;
        }
      }
      else {
        if (!islower(in[1])) {
          el = (elements) i;
          return ++in;
        }
        else
          find = i;
      }
    }
  }
  if (find) {
    el = (elements) find;
    in++;
  }
  return in;
}

istream& operator>>(istream &in, elem &to) {
  char buf[2];
  in >> ws;
  in.get(buf[0]);
  if (buf[0] == EOF) return in;
  in.get(buf[1]);
  to.anal(buf);
  if (!to.el) {
    in.putback(buf[0]);
    if (buf[1] != EOF)
      in.putback(buf[1]);
  }
  else if (to.len() == 1)
    in.putback(buf[1]);
  return in;
}

void formula::sort(int &ne, elem *ele, int *c) {
  if (ne == 1) return;
  elem tele;
  int tc;
  int i, j;
  int notSorted;
  do {
    notSorted = 0;
    for (i = 0; i + 1 < ne; i++) {
      if (ele[i] < ele[i + 1]) {
        notSorted = 1;
        tele = ele[i];
        ele[i] = ele[i + 1];
        ele[i + 1] = tele;
        tc = c[i];
        c[i] = c[i + 1];
        c[i + 1] = tc;
      }
      else if (ele[i] == ele[i + 1]) {
        c[i] += c[i + 1];
        for (j = i + 1; j < ne; j++) {
          ele[j] = ele[j + 1];
          c[j] = c[j + 1];
        }
        i--;
        ne--;
      }
    }
  }
  while (notSorted);
}

void formula::copy(const formula &old) {
  clear();
  if ((ne = old.ne) != 0) {
    if ((elems = new elem[ne]) == NULL || (coef = new int[ne]) == NULL
                    || (frml = new char[strlen(old.frml) + 1]) == NULL)
    {
      clear();
      return;
    }
    memcpy(elems, old.elems, sizeof(elem)*ne);
    memcpy(coef, old.coef, sizeof(int)*ne);
    strcpy(frml, old.frml);
  }
}

const char* formula::anal(const char *in) {
  clear();
  if (*in == '\0') return in;
  elem ele[BUF_EL];
  int c[BUF_EL];
  const char *start = in;
  int i = 0;
  while (isspace(*in)) in++;
  if (!read(i, in, ele, c)) {
    clear();
    return in;
  }
  ne = i + 1;
  if (!ne) return in;
  sort(ne, ele, c);
  if ((elems = new elem[ne]) == NULL || (coef = new int[ne]) == NULL
              || (frml = new char[in - start + 1]) == NULL) {
    clear();
    return in;
  }
  memcpy(elems, ele, sizeof(elem)*ne);
  memcpy(coef, c, sizeof(int)*ne);
  memcpy(frml, start, in - start);
  frml[in - start] = '\0';
  return in;
}

int formula::read(int &i, const char* &in, elem *ele, int *c) {
  for (i; ;i++) {
    if (i == BUF_EL - 1) {
      return 0;
    }
    in = ele[i].anal(in);
    if (!ele[i]) {
      if (*in == '(') {
        int starti = i;
        in++;
        if (isspace(*in)) return 0;
        if (!read(i, in, ele, c)) {
          return 0;
        }
        if (*in != ')')
          return 0;
        in++;
        if (isdigit(*in) || (in[0] == '-' && isdigit(in[1]))) {
          int num = atoi(in);
          in++;
          while (isdigit(*in)) in++;
          for (int j = starti; j <= i; j++)
            c[j] *= num;
          if (!num) i = starti - 1;
        }
        continue;
      }
      else
      {
        i--;
        return 1;
      }
    }
    if (isdigit(*in) || (in[0] == '-' && isdigit(in[1]))) {
      c[i] = atoi(in);
      in++;
      while (isdigit(*in)) in++;
      if (!c[i]) i--;
    }
    else
      c[i] = 1;
    if (*in == '-') {
      in++;
      i++;
      ele[i] = elem::e;
      c[i] = 1;
      return 1;
    }
    if (*in == '+') {
      in++;
      i++;
      ele[i] = elem::e;
      c[i] = -1;
      return 1;
    }
    if (isspace(*in)) return 1;
  }
}

char* formula::to_str(char *buf) const {
  strcpy(buf, frml);
/* molecular formula
  char *str = buf;
  for (int i = 0; i < ne; i++) {
    if (i == ne - 1 && elems[i] == elem::e) {
      if (coef[i] == 1) {
        *str++ = '-';
        break;
      }
      if (coef[i] == -1) {
        *str++ = '+';
        break;
      }
    }
    elems[i].to_str(str);
    str += strlen(str);
    if (coef[i] != 1)
      str += strlen(itoa(coef[i], str, 10));
  }
  *str = '\0';
*/
  return buf;
}

int formula::compare(const formula &second) const {
  int diff;
  int i = 0;
  while (i < ne && i < second.ne) {
    if ((diff = elems[i] - second.elems[i]) != 0) return diff;
    if ((diff = coef[i] - second.coef[i]) != 0) return diff;
    i++;
  }
  if (ne - second.ne)
    return (ne - second.ne);
  else
    return strcmp(frml, second.frml);
}

ostream& operator<<(ostream &out, const formula &old) {
  return out << old.frml;
}

istream& operator>>(istream &in, formula &to) {
  char buf[BUF_EL*5];
  char *ptr = buf;
  char ch;
  in >> ws;
  while (ch = in.peek(),
         (!isspace(ch) &&
          ch != EOF &&
          ch != ';' &&
          ch != '"' &&
          ch != '[' &&
          ptr < buf + BUF_EL*5 - 1))
    in.get(*ptr++);
  *ptr = '\0';
  const char *buf1 = to.anal(buf);
  for (; buf1 < ptr; buf1++)
    in.putback(*buf1);
  return in;
}

int formula::len() const {
/* len of molecular formula
  char buffer[17];
  int len = 0;
  for (int i = 0; i < ne; i++) {
    if (i == ne - 1 && elems[i] == elem::e) {
      if (coef[i] == 1 || coef[i] == -1) {
        len++;
        break;
      }
    }
    len += elems[i].len();
    if (coef[i] != 1) len += strlen(itoa(coef[i], buffer, 10));
  }
  return len;
*/
  return frml ? strlen(frml) : 0;
}

int formula::noa(elem ele) {
  for (int i = 0; i < ne; i++)
    if (elems[i] == ele) return coef[i];
  return 0;
}

elem formula::el(unsigned int i) {
  return (i < ne) ? elems[i] : elem();
}

double formula::mass() {
  double mass = 0.;
  for (int i = 0; i < ne; i++)
    mass += coef[i]*elems[i].mass();
  return mass;
}

const char* label::anal(const char *in) {
  clear();
  while (isspace(*in)) in++;
  if (*in == '\0') return in;
  if (*in != del1) return in;
  in++;
  const char *s = in;
  while (*in != del2) {
    if (in == '\0') return in;
    in++;
  }
  str = new char[in - s + 1];
  if (str) {
    char *s1 = str;
    while (s < in) *s1++ = *s++;
    *s1 = '\0';
  }
  return ++in;
}

istream& operator>>(istream &in, label &to) {
  ostrstream out;
  char ch;
  to.clear();
  in >> ws;
  if (in.peek() != to.del1) return in;
  in.get(ch);
  while (in.get(ch), ch != to.del2) {
    if (ch == EOF) {
      delete out.str();
      return in;
    }
    out << ch;
  }
  in.get(ch);
  out << ends;
  if (!strlen(out.str()))
  {
    delete out.str();
    return in;
  }
  to.str = new char[strlen(out.str()) + 1];
  if (to.str) {
    strcpy(to.str, out.str());
  }
  delete out.str();
  return in;
}

species::species_err species::errno;

void species::copy(const species &old) {
  errno = OKAY;
  formula::copy(old);
  function::copy(old);
  if (function::errno != function::OKAY) errno = WRFNC;
  isomer.copy(old.isomer);
  phase.copy(old.phase);
}

const char* species::anal(const char *in) {
  const char *tmp;
  errno = OKAY;
  in = formula::anal(in);
  tmp = isomer.anal(in);
  if (tmp > in + 2 && !isomer)
  {
    errno = WRLBL;
    clear();
    return in;
  }
  in = tmp;
  tmp = phase.anal(in);
  if (tmp > in + 2 && !phase)
  {
    errno = WRLBL;
    clear();
    return in;
  }
  in = tmp;
  while (isspace(*in)) in++;
  in = function::anal(in);
  if (function::errno != function::OKAY) {
    errno = WRFNC;
    clear();
  }
  return in;
}

char* species::to_str(char *buf) const {
  errno = OKAY;
  char *str = buf;
  formula::to_str(str);
  str += strlen(str);
  if (isomer) {
    isomer.to_str(str);
    str += strlen(str);
  }
  if (phase) {
    phase.to_str(str);
    str += strlen(str);
  }
  function::to_str(str);
  if (function::errno != function::OKAY) errno = WRFNC;
  return buf;
}

int species::compare(const species &second) const {
  int is;
  if ((is = formula::compare(second)) != 0) return is;
  if ((is = isomer.compare(second.isomer)) != 0) return is;
  if ((is = phase.compare(second.phase)) != 0) return is;
  return 0;
}

int species::len() const {
  int len = 0;
  len += formula::len();
  if (isomer) len += isomer.len();
  if (phase) len += phase.len();
  len += function::len();
  return len;
}

ostream& operator<<(ostream &out, const species &old) {
  species::errno = species::OKAY;
  out << *(formula*)&old;
  out << old.isomer;
  out << old.phase;
  if (old.function::operator int()) {
    out << endl;
    out << *(function*)&old;
    if (function::errno != function::OKAY) species::errno = species::WRFNC;
  }
  return out;
}

istream& operator>>(istream &in, species &to) {
  species::errno = species::OKAY;
  in >> *(formula*)&to;
  in >> to.isomer;
  in >> to.phase;
  in >> *(function*)&to;
  if (function::errno != function::OKAY) {
    species::errno = species::WRFNC;
    to.clear();
  }
  return in;
}

double species::H(const double &T, const double &p) {
  double G = est(T, p);
  if (G == HUGE_VAL) return HUGE_VAL;
  double S = df(global::dT, T, p);
  if (S == HUGE_VAL) return HUGE_VAL;
  return G - T*S;
}

double species::S(const double &T, const double &p) {
  double S = df(global::dT, T, p);
  if (S == HUGE_VAL) return HUGE_VAL;
  return -S;
}

double species::Cp(const double &T, const double &p) {
  double C = df(global::dT2, T, p);
  if (C == HUGE_VAL) return HUGE_VAL;
  return -C*T;
}