
#include "jform.h"

#define MAXCLE          703   /* 26*26+26+1      */

static char* TM[MAXELM] = { "@@",
  "H@", "HE", "LI", "BE",   "B@", "C@", "N@", "O@",
  "F@", "NE", "NA", "MG",   "AL", "SI", "P@", "S@",
  "CL", "AR", "K@", "CA",   "SC", "TI", "V@", "CR",
  "MN", "FE", "CO", "NI",   "CU", "ZN", "GA", "GE",
  "AS", "SE", "BR", "KR",   "RB", "SR", "Y@", "ZR",
  "NB", "MO", "TC", "RU",   "RH", "PD", "AG", "CD",
  "IN", "SN", "SB", "TE",   "I@", "XE", "CS", "BA",
  "LA", "CE", "PR", "ND",   "PM", "SM", "EU", "GD",
  "TB", "DY", "HO", "ER",   "TM", "YB", "LU", "HF",
  "TA", "W@", "RE", "OS",   "IR", "PT", "AU", "HG",
  "TL", "PB", "BI", "PO",   "AT", "RN", "FR", "RA",
  "AC", "TH", "PA", "U@",   "OX"
};

// tableau de valeurs positives ou nulles sauf l'��ent 0 symbolique = -1
static double MASSE[MAXELM] = { -1,
    1.008,  4.003,  6.939,  9.012,   10.811, 12.011, 14.007, 15.999,
   18.998, 20.183, 22.991, 24.312,   26.982, 28.086, 30.974, 32.064,
   35.453, 39.948, 39.102, 40.08 ,   44.956, 47.90 , 50.942, 51.996,
   54.938, 55.847, 58.933, 58.71 ,   63.54 , 65.37 , 69.72 , 72.59 ,
   74.922, 78.96 , 79.909, 83.80 ,   85.47 , 87.62 , 88.905, 91.22 ,
   92.906, 95.94 , 98.0  ,101.07 ,  102.905, 106.4 ,107.870,112.40 ,
   114.82,118.69 ,121.75 ,127.60 ,  126.904, 131.30,132.905,137.34 ,
   138.91,140.12 ,140.907,144.24 ,  147.0  ,150.35 ,151.96 ,157.25 ,
  158.924,162.50 ,164.930,167.26,   168.934, 173.04, 174.97, 178.49,
  180.948,183.85 ,186.2  ,190.2  ,  192.2, 195.09 , 196.967, 200.59,
  204.37 ,207.19 , 208.980, 210.0,  210.0,  222.0, 223.0, 226.0,
  227.0  ,232.038, 231.0 ,238.03 ,  18.0
};

// tableau de valeurs positives ou nulles sauf l'��ent 0 symbolique = -1
// et ceux pour lesquels la valeur est inconnue, on fait aussi -1
// le calcul de la valeur formulaire retournera une valeur <0 si un
// ��ent est inconnue ou si la propri��de cet ��ent est inconnue.
//
static double ENTROPIE[MAXELM] = {     -1,
   65.340,     -1, 29.120,  9.540,     -1,  5.740, 95.806, 102.573,
  101.398,     -1, 51.300, 32.680, 28.350, 18.810,     -1,  32.054,
  111.538,     -1, 64.680, 41.600,     -1, 30.720,     -1,      -1,
   32.010, 27.276, 30.040, 29.870, 33.140, 41.630,     -1,      -1,
   35.690,     -1,     -1,     -1,     -1,     -1,     -1,      -1,
     0   , 28.570,     -1,     -1,     -1,     -1, 42.550,      -1,
       -1, 51.180,     -1,     -1,     -1,     -1,     -1,  62.500,
       -1,     -1,     -1,     -1,     -1,     -1,     -1,      -1,
       -1,     -1,     -1,     -1,     -1,     -1,     -1,      -1,
       -1,     -1,     -1,     -1,     -1,     -1, 47.490,  75.900,
       -1, 64.800,     -1,     -1,     -1,     -1,     -1,      -1,
       -1,     -1,     -1, 50.210,     -1
};

static int cle[MAXCLE];

static void
ini_cle()
{
int i, j;
int was_init = 0;
//
  if (was_init++) {
    return;
  }
  else {
    for ( i=0; i<MAXELM; i++) {
      j = (TM[i][0] - '@') * 26 + TM[i][1] - '@';
      cle[j] = i;
    }
  }
}

void
symbel(int *za, char *symb){
//
  if ((*za<0)||(*za>=MAXELM)) za[0]=0;
  strncpy(symb,TM[*za],sizeof symb);
  if (symb[1]=='@')
    symb[1]=' ';
  else
    symb[1]=tolower(symb[1]);
}

/* function sur le mod�e xxx strtox(char*, char**) */
int zat(char *s, char **t){
  int r, j, k0, k1;
//
  ini_cle();
  *t = s;
  if (!isupper(*s)) return 0;
  k0 = s[0]-'@';
  *t += 1;
  if (islower(**t)){
    k1 = toupper(s[1])-'@';
    *t += 1;
  }
  else
    k1 = 0;
  j = 26*k0+k1;
  r = ((j>=0)&&(j<MAXCLE))?cle[j]:0;
  return r;
}

/* fonctionne comme strtod...  */
static double
wat(double *pelm, char *s, char **t, set sys){
  card z;
// pelm est le vecteur des masses atomiques ou des entropies
// rtourne le z-atom. dans *z
// type card h�it�de set.h
  z = zat(s,t);
  if (sys!=NULL) set_add(sys,z);
  return pelm[z];
}

static char
fw1(double *pelm, char *f, char **g, double *w, set sys){
  double coeff;
  double aw;
//
  if isupper(f[0]){ // double emploi avec test dans zat ? NON
    aw = wat(pelm, f, &f, sys);
    if ((*w=aw)<0.0) return EOF; // ��ent ou valeur inconnue
                                 // la valeur 0 est possible
                                 // utilis� pour calcul de stoechio.
// symbole valide
    if (isdigit(*f))
      coeff = strtod(f,&f);
    else
      coeff = 1.0;
// masse atomique * coeff. steochio
    *w = aw*coeff;
  }
  else {  // symbole non valide prem. car. n'est pas une majuscule
    *w = 0.0;
  }
  *g = f;
  return f[0];
}

//----------------------------------------------

double
fchtod(double *pel, char *f, char **g, set sys){
  double coeff, w, w0;
  char c;
  pil p;
// syst=NULL si on ne veut pas sys, l'ensemble des ��ents
// qui composent la formule f.
// ACHTUNG les propr. elem. DOIVENT �re >= 0
// la valeur -1.0 pour l'��ent fictif de z=0
// est utilis� comme sortie pour symbole erron�// ou valeur de la propri��inconnue
  if (f==NULL) return -1.0;
  while (isspace(f[0])) ++f;
  p=pil_new(sizeof(double));
  w = 0.0;
  while(1){
    c=fw1(pel,f,&f,&w0,sys);
    if (c=='('){
      w += w0;
      pil_emp(p,&w);
      w = 0.0;
      ++f;
    }
    else if (c==')'){
      w += w0;
      ++f;
      if
        isdigit(f[0]) coeff = strtod(f,&f);
      else
        coeff = 1.0;
      w *= coeff;
      pil_dep(p,&w0);
      w += w0;
    }
    else {  // majuscule ou ponctuation ou fin ou erreur
      w += w0;
      if (!isupper(c)){
        if (isspace(c)||(c=='\0')){ // sortie normale
          break;
        }
        else if ((c==',')||(c==';')) { // , ; admis comme s�arateurs
          ++f;      // peut retourner 0 mais on peut continuer
          break;
        }
        else{  // ni une majuscule ni un blanc ni la fin de chaine
          w = -1.0; // on ne peut pas continuer
          break;
        }
      }
    }
  }
  if(pil_pro(p)!=0) return -2.0;
  pil_free(p); // il faudrait la repasser !
  *g=f; // aspirer les blancs s'il y a lieu ?
  return (w);
}

// Prototypes compatibles avec Scilab

bool
sysapp(char *ph, char *sy, int *out){
  static set s1;
  static set s2;
//
  if (s1==NULL)
    s1=set_new(MAXELM);
  else
    set_mt(s1);
  if (s2==NULL)
    s2=set_new(MAXELM);
  else
    set_mt(s2);
  fchtod(MASSE,ph,&ph,s1);
  fchtod(MASSE,sy,&sy,s2);
  return *out=set_incl(s2,s1);
}

extern void         // p1:symbole -> p2:Numero atomique
za(char *f, int *z){
   *z = zat(f,&f);  // voir si f est modifi�au dessus ?
}

void
fw(char *f, double *w){
  *w=fchtod(MASSE,f,&f,NULL);
}

void
fs(char *f, double *s){
// S_formulaire_standard en kJ
  *s=fchtod(ENTROPIE,f,&f,NULL)*1.0e-3;
}

void
fs298(char *f, double *s){
// 298.15*S_formulaire_standard en kJ
  *s=fchtod(ENTROPIE,f,&f,NULL)*0.29815;
}

void
stoe(char *f, int *z, double *s){
  static double stoechio[MAXELM];
  static int was_init = 0;
//
//in   f : formule
//in   z : z-atom de l'��ent
//out  s : coeff. stoechio de l'��ent dans la formule
//
  //? je ne sais plus si c'est n�essaire ou si init 0 par def. ?
  assert((*z>0)&&(*z<MAXELM));
  if (~(was_init++)){
    memset(stoechio,0,sizeof(double)*MAXELM);
    stoechio[0]=-1.0;
  }
  stoechio[*z]=1.0;
  *s=fchtod(stoechio,f,&f,NULL);
  stoechio[*z]=0.0;
}

void
sysadd(char *f, int *z, int *n, int *zz, int *nn){
// ajouter au syst�e z(n) les ��ents de la formule f
  static int was_init=0;
  static set s1;
  static set s2;
  static set s3;
  static vec v3;
  static double pelm[MAXELM]; // je suppose init 0 �v�ifier
                              // �voir aussi dans stoe
  double r;
  int i;
  card k; // sic ACHTUNG au type casting
//
  if (~(was_init++)){
    s1=set_new(MAXELM);
    s2=set_new(MAXELM);
    s3=set_new(MAXELM);
    v3=vec_new(sizeof(card));
  }
  else {
    set_mt(s1);
    set_mt(s2);
    set_mt(s3);
    vec_mt(v3);
  }
  pelm[0]=-1.0;
  r=fchtod(pelm,f,&f,s1); // s1 ensemble des elem. de f
//  assert(r==0); // scilab explose avec le message dans la console
  if (r!=0) { // moins violent que assert !
    set_mt(s1); set_add(s1,0);
  }
  for (i=0;i<*n;i++) set_add(s2,(card)z[i]);
  set_union(s1,s2,s3);
  // conversion de s3 -> int[] merdique la fonction set_2vec
  // exige un vecteur de card (unsigned short sur 2 octets)
  set_2vec(s3,v3);
//
  *nn=(int)vec_nbr(v3);
  for (i=0;i<*nn;i++){
    vec_elm(v3,i,&k);
    zz[i]=(int)k;
  }
}

int  // f appartient au syst�e z
sysin(char *f, int *z, int *n, int *ret){
  static double junk[MAXELM];
  static set s1;
  static set s2;
  static int was_init=0;
  double r;
  int i;
//
  if (~was_init++){
    junk[0]=-1.0; // cf. fchtod()
    s1=set_new(MAXELM);
    s2=set_new(MAXELM);
  }
  else{
    set_mt(s1);
    set_mt(s2);
  }
  r=fchtod(junk,f,&f,s1);
  assert(r==0);
  for (i=0;i<*n;i++) set_add(s2,(card)z[i]);
  return *ret=set_incl(s2,s1);
}
