units.h

Go to the documentation of this file.

//-----------------------------------------------------------------------bl-
//--------------------------------------------------------------------------
//
// Antioch - A Gas Dynamics Thermochemistry Library
//
// Copyright (C) 2014-2016 Paul T. Bauman, Benjamin S. Kirk,
//                         Sylvain Plessis, Roy H. Stonger
//
// Copyright (C) 2013 The PECOS Development Team
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the Version 2.1 GNU Lesser General
// Public License as published by the Free Software Foundation.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc. 51 Franklin Street, Fifth Floor,
// Boston, MA  02110-1301  USA
//
//-----------------------------------------------------------------------el-
#ifndef ANTIOCH_UNITS_H
#define ANTIOCH_UNITS_H

//Antioch
#include "antioch/metaprogramming.h"
#include "antioch/antioch_asserts.h"
#include "antioch/siprefix.h"
#include "antioch/insi.h"
#include "antioch/converter.h"
#include "antioch/unit_defs.h"

//C++
#include <vector>
#include <cmath>
#include <sstream>
#include <iomanip>

namespace Antioch{

template <typename T = double>
class Units{
      public:

        Units(const std::string &sym, const std::string &na,
                const T &conva, const T &convb,
                int mi,int kgi=0, int si=0, int Ai=0, int Ki=0, int moli=0, int cdi=0, int radi=0);//constructors for known Units
        Units(const std::string &sym, const std::string &na="");  //constructors when automatic conversion
        Units(const std::string & sym, const Converter<T> &conv, const std::string &na="");//constructor for given conversion
        Units(){}

        ~Units(){}


        void print(std::ostream &os = std::cout) const;


       friend std::ostream& operator<<( std::ostream& os,
                                        const Units<T>& unit )
       {
         unit.print(os);
         return os;
       } 


        bool operator== (const Units<T> &rhs) const;
        bool operator!= (const Units<T> &rhs) const {return !(*this == rhs);}
        Units<T> &operator= (const Units<T> & rhs);
        Units<T> &operator+=(const Units<T> & rhs);
        Units<T> &operator-=(const Units<T> & rhs);
        Units<T>  operator+ (const Units<T> & rhs) const;
        Units<T>  operator- (const Units<T> & rhs) const;
        Units<T> &operator*=(int r);
        Units<T>  operator* (int r)             const;
        Units<T> &operator/=(int r);
        Units<T>  operator/ (int r)             const;

        void equalize(const Units<T> &rhs)      {*this  = rhs;}
        void equalize(Units<T> *rhs)            {*this  = *rhs;}
        void equalize(std::string unit)       {*this  = Units<T>(unit);}
        void add(const Units<T> &rhs)            {*this += rhs;}
        void add(Units<T> *rhs)                  {*this += *rhs;}
        void add(std::string unit)            {*this += Units<T>(unit);}
        void substract(const Units<T> &rhs)      {*this -= rhs;}
        void substract(Units<T> *rhs)            {*this -= *rhs;}
        void substract(const std::string &unit)      {*this -= Units<T>(unit);}
        void root(int r)                      {*this /= r;}

        void set_unit(const std::string &sym, std::string na);
        void set_unit(const std::string &sym)       {set_unit(sym,"");}
        void set_name(const std::string &na)        {name = na;}
        void set_symbol(const std::string &symb)    {symbol = symb;}
        void set_SI_converter(const Converter<T> &conv) {Antioch::init_clone(toSI,conv);}
        void set_power_to_SI(const InSI &mat)         {power = mat;}

        bool is_homogeneous(const Units<T> &rhs) const;
        bool is_homogeneous(std::string target) const;
        bool is_united() const {return !power.empty();}


        T factor_to_some_unit(const Units<T> & target)     const;
        T factor_to_some_unit(const std::string &target)    const;
        T translator_to_some_unit(const Units<T> & target)  const;
        T translator_to_some_unit(const std::string &target)  const {return translator_to_some_unit(Units<T>(target));}

        const std::string get_name()            const {return name;}
        const std::string get_symbol()          const {return symbol;}
        const Converter<T> get_SI_coef() const {return toSI;}
        T get_SI_factor()                const {return toSI.geta();}
        T get_SI_translator()            const {return toSI.getb();}
        const InSI get_power()                  const {return power;}
        const std::string get_SI_symb()     const;
        int get_SI_power(const std::string &SIask) const;
        const std::string get_SI_convenient_symb() const;
        const std::string harmonized_symbol(const std::string &input = "")    const {return manipulate_symbol(input,false);}
        const std::string contracted_symbol(const std::string &input = "")    const {return manipulate_symbol(input,true);}

        void clear(){symbol.clear();name.clear();toSI.clear();power.clear();}

        void showAll(std::ostream &out = std::cout);

      private:

        const std::string manipulate_symbol(std::string input, bool contract) const;


        void fill_in_power(bool doConv);
        bool parse_single_unit(int signe, std::string unit, bool doConv);
        void parse_prefix_unit(int &iUnit,int &iPre, const std::string &unit)  const;
        int indexUnit(std::string unit)                               const;
        int parse_power(std::string unit, int &nc)                    const;
        bool is_number(char c)                                   const {return ((int)c <= 57 && (int)c >= 48);} //ascii table

        Converter<T> raise(const Converter<T> &tbm,int power)       const;
        Converter<T> raise(const Converter<T> &tbm,double power)       const;
        void develop_symbol(std::string &subsymbol);
        void reverse_power_symbol(std::string &subsymbol);
        bool is_in_symb(char c)               const;
        int n_dimension_of_units()             const;
        const std::string add_SI(int pow,std::string SIsymb) const;
        void symbol_to_the_power(int r, const int &key);
        int get_integer_power(int unit,int r, const int &key);
        void check_if_name(std::string &nameOut, const Units<T> &un) const;

        std::string symbol;
        std::string name;
        Converter<T> toSI;
        InSI power;

};


template<typename T>
inline
Units<T>::Units(const std::string &sym, const std::string &na,
             const T &conva, const T &convb,
             int mi,int kgi, int si, int Ai, int Ki, int moli, int cdi, int radi): //fully descriptive constructor
  symbol(sym),
  name(na),
  toSI(conva,convb),
  power(mi,kgi,si,Ai,Ki,moli,cdi,radi)
{
  develop_symbol(symbol);
}

template<typename T>
inline
Units<T>::Units(const std::string &sym, const std::string &na):  //constructors when automatic conversion
       symbol(sym),
       name(na)
{
  fill_in_power(true);
}

template<typename T>
inline
Units<T>::Units(const std::string &sym, const Converter<T> &conv, const std::string &na)://constructor for given conversion
       symbol(sym),
       name(na),
       toSI(conv)
{
  fill_in_power(false);
}

template <typename T>
void Units<T>::fill_in_power(bool doConv)
{
  if(symbol.empty())return; // no unity
  develop_symbol(symbol);

  std::string tmp(""),symboltmp(symbol);
  int signe(1),istart(0);

  while(symboltmp != contracted_symbol(symboltmp))symboltmp = contracted_symbol(symboltmp);
  if(symboltmp.empty())return; // no unity

  if(symboltmp[0] == '/')
  {
    signe = -1;
    istart = 1;
  }
  for(unsigned int i = istart ; i < symboltmp.size() ; i++)
  {
    if(symboltmp[i] == '.')
    {
      if(!parse_single_unit(signe,tmp,doConv))
      {
        antioch_unit_error("In symbol " + symboltmp + ", unit \"" + tmp + "\" does not ring a bell");
        break;
      }
      signe = 1;
      tmp.clear();
      }else if(symboltmp[i] == '/')
      {
       if(!parse_single_unit(signe,tmp,doConv))
       {
         antioch_unit_error("In symbol " + symboltmp + ", unit \"" + tmp + "\" does not ring a bell");
         break;
       }
       signe = -1;
       tmp.clear();
       }else
       {
         tmp += symboltmp[i];
       }
   }
   if(!parse_single_unit(signe,tmp,doConv))
      antioch_unit_error("In symbol " + symboltmp + ", unit \"" + tmp + "\" does not ring a bell");
}

template <typename T>
void Units<T>::develop_symbol(std::string &subsymb)
{
  if(subsymb == "no unit" || subsymb == "No unit" || subsymb == "NO UNIT")
  {
     subsymb.clear();
     return;
  }
  if(subsymb.find("(") == std::string::npos)return;

//model is a serie ().()/()..., we find it
// we count the opening parenthesises (no) and
// the closing parenthesises (nc). We need to
// keep track of the position of opening (po)
// and closing (pc) parenthesises. We have a pair
// when no == nc.
  unsigned int no(0),nc(0);
  unsigned int po(0),pc(subsymb.size() - 1);
  for(unsigned int cc = 0; cc < subsymb.size(); cc++)
  {
    if(subsymb[cc] == '(')
    {
      no++;
      if(po == 0 && cc != 0)po = cc;
    }
    if(subsymb[cc] == ')')
    {
      nc++;
      pc = cc;
    }
    if(no == 0)continue;

    if(no == nc)//found a pair
    {
     //develop it
      if(pc == po + 1)
      {
        unsigned int off(0);
        if(po == 0)off = 1;//if ( is the first character or not
        subsymb.erase(po - 1 + off,3 - off);//if yes, suppress '()', if not suppress '.()' or '/()'
      }else
      {
        std::string insideStr = subsymb.substr(po + 1,pc - po - 1);
        develop_symbol(insideStr);
        if(po != 0)if(subsymb[po - 1] == '/')reverse_power_symbol(insideStr);
        subsymb.replace(po,pc - po + 1,insideStr);
      }
     //reset the system
      no = 0;
      nc = 0;
      po = 0;
      pc = subsymb.size() - 1;
      cc -= 2;
    }
  }

//if first character is a power determinant ('/' or '.')
  if(subsymb[0] == '/') //we change only the first atomic unit
  {
    subsymb.erase(0,1);
    po = subsymb.find(".");
    if(po > subsymb.find("/"))po = subsymb.find("/");
    std::string curUnit = subsymb.substr(0,po);
    int nc(0),pow = - parse_power(curUnit,nc); //power
    if(pow != 1)
    {
      std::ostringstream np;
      np << pow;
      if(po < subsymb.size())
      {
         subsymb.replace(po - nc,nc,np.str());
      }else
      {
        subsymb += np.str();
      }
    }else
    {
      subsymb.erase(po - nc,nc);
    }
  }else if(subsymb[0] == '.')
  {
    subsymb.erase(0,1);
  }

}

template <typename T>
T Units<T>::factor_to_some_unit(const Units<T> &target) const
{
  if(is_homogeneous(target))
  {
    return get_SI_factor()/target.get_SI_factor();
  }else
  {
    antioch_unit_error("Units are not homogeneous:\n\"" + symbol + "\" and \"" + target.get_symbol() + "\".");
    return -1.;
  }
}

template<typename T>
T Units<T>::factor_to_some_unit(const std::string &target) const 
{
  return factor_to_some_unit(Units<T>(target));
}

template <typename T>
T Units<T>::translator_to_some_unit(const Units<T> & target)  const
{
  if(is_homogeneous(target))
  {
    return ((this->get_SI_translator() - target.get_SI_translator())/target.get_SI_factor());
  }
  else
  {
    antioch_unit_error("Units are not homogeneous.");
    return -1.;
  }
}

template <typename T>
void Units<T>::parse_prefix_unit(int &iUnit,int &iPre,const std::string& unit) const
{
  iPre = -1;

  iUnit = UnitBaseStorage::known_units().stored_index(unit);
//prefix, if exists, is one or two character
   if(iUnit == -1 && unit.size() > 1)
   {
     std::string pre = unit.substr(0,1);
     std::string un  = unit.substr(1,std::string::npos);
     iPre = UnitBaseStorage::known_prefixes().stored_index(pre);
     iUnit = UnitBaseStorage::known_units().stored_index(un);
     if(iPre == -1 && iUnit == -1 && unit.size() > 2)
     {
       pre = unit.substr(0,2);
       un  = unit.substr(2,std::string::npos);
       iPre = UnitBaseStorage::known_prefixes().stored_index(pre);
       iUnit = UnitBaseStorage::known_units().stored_index(un);
     }
   }
   return;
}

template <typename T>
bool Units<T>::parse_single_unit(int signe,std::string unit,bool doConv)
{
  int iUnit = -1, nc = 0;
  int ipower = parse_power(unit,nc); // find power
  int iPre = -1;

  unit = unit.substr(0,unit.length()-nc); //unit without power

  parse_prefix_unit(iUnit,iPre,unit);//find prefix and unit

  if(iUnit == -1)return false; //not found

  T pre = 1.L;
  if(iPre != -1)pre = UnitBaseStorage::known_prefixes().stored(iPre).value<T>();
  if(UnitBaseStorage::known_units().stored(iUnit).symbol() == "kg" && 
     unit != "kg")pre *= 1e-3L; //rescale
  InSI powerTmp = UnitBaseStorage::known_units().stored(iUnit).power_array();
  power += (powerTmp * signe * ipower);
  if(doConv)
  {
    Converter<T> convTmp = UnitBaseStorage::known_units().stored(iUnit).converter() * pre;
    toSI += raise(convTmp,signe * ipower);
  }

  return true;
}

template <typename T>
int Units<T>::indexUnit(std::string unit) const
{
  int iUnit;

  if(unit.empty())return -1;

  if(unit == "g")//special case to adapt SI
  {
    unit = "kg";
  }
  for (iUnit = 0; iUnit < UnitBaseStorage::known_units().n_known_units() ; iUnit++)
  {
    if(unit == UnitBaseStorage::known_units().stored(iUnit).symbol())break;
  }

  if(iUnit >= UnitBaseStorage::known_units().n_known_units())iUnit=-1;


  return iUnit;
}

template <typename T>
int Units<T>::parse_power(std::string unit,int &nc) const
{
  int ip = 1, loc = unit.length();
  char c = unit[loc-1];
  nc = 0;
  while(is_number(c)){
    nc++;
    loc--;
    c = unit[loc-1];
  }

  if(unit[loc-1] == '-'){
        loc--;
        nc++;
  }

  std::string power = unit.substr(loc,std::string::npos);
  if(power.size() > 0)
  {
     std::stringstream p;
     p << power;
     p >> ip;
  }

  if(ip == 0)antioch_unit_error("Invalid power found: " + unit);

  return ip;
}

template <typename T>
bool Units<T>::is_homogeneous(const Units<T> &rhs) const 
{
  return (power == rhs.get_power());
}

template <typename T>
bool Units<T>::is_homogeneous(std::string target) const
{
  if(target.empty())
  {
    return (!this->is_united());
  }else
  {
    Units<T> rhs(target);
    return is_homogeneous(rhs);
  }
  
}

template <typename T>
int Units<T>::get_SI_power(const std::string &SIask) const
{
  if(SIask == "m")  return power.get_m();
  if(SIask == "kg") return power.get_kg();
  if(SIask == "s")  return power.get_s();
  if(SIask == "A")  return power.get_A();
  if(SIask == "K")  return power.get_K();
  if(SIask == "mol")return power.get_mol();
  if(SIask == "cd") return power.get_cd();
  if(SIask == "rad")return power.get_rad();

  antioch_unit_error(SIask + " is not a SI symbol");

  return 0;
}

template <typename T>
std::string const Units<T>::get_SI_symb() const
{
 std::string SISymb;
 SISymb.clear();

 if(power.get_m() != 0)SISymb += add_SI(power.get_m(),"m");

 if(power.get_kg() != 0)
 {
   if(!SISymb.empty())SISymb += ".";
   SISymb += add_SI(power.get_kg(),"kg");
 }
 if(power.get_s() != 0)
 {
   if(!SISymb.empty())SISymb += ".";
   SISymb += add_SI(power.get_s(),"s");
 }
 if(power.get_A() != 0)
 {
   if(!SISymb.empty())SISymb += ".";
   SISymb += add_SI(power.get_A(),"A");
 }
 if(power.get_K() != 0)
 {
   if(!SISymb.empty())SISymb += ".";
   SISymb += add_SI(power.get_K(),"K");
 }
 if(power.get_mol() != 0)
 {
   if(!SISymb.empty())SISymb += ".";
   SISymb += add_SI(power.get_mol(),"mol");
 }
 if(power.get_cd() != 0)
 {
   if(!SISymb.empty())SISymb += ".";
   SISymb += add_SI(power.get_cd(),"cd");
 }
 if(power.get_rad() != 0)
 {
   if(!SISymb.empty())SISymb += ".";
   SISymb += add_SI(power.get_rad(),"rad");
 }

 return SISymb;

}

template <typename T>
const std::string Units<T>::add_SI(int pow,std::string symb) const
{
 std::string out("");

 if(pow == 0)return out;

 out = symb;
 if(pow != 1)
 {
   std::stringstream po;
   po << pow;
   out += po.str();
 }
 
 return out;

}

template <typename T>
std::string const Units<T>::manipulate_symbol(std::string input, bool contract) const
{
  std::string harmSymb("");
  if(input.empty())input = symbol;

  std::vector<std::string> unitvec;
  std::vector<int> powervec;
  std::string curUnit(""),interUnit(".");//,strPower;

  for(unsigned int i = 0; i < input.size(); i++)
  {
    if(input[i] != '.' && input[i] != '/')curUnit += input[i];
    if(input[i] == '.' || input[i] == '/' || i == input.size() - 1 )
    {
      if(curUnit.empty())continue;
//first, parsing the unit
      int nc(0);
      int iUnit(-1),iPre(-1);
      int curPower = parse_power(curUnit,nc); //power
//      strPower = curUnit.substr(curUnit.length() - nc,std::string::npos); //resulting power
      curUnit = curUnit.substr(0,curUnit.length() - nc); //resulting unit
      parse_prefix_unit(iUnit,iPre,curUnit);
 // unit is UnitBaseStorage::Prefixes[iPre].symbol() + UnitBaseStorage::knownUnits[iUnit].get_symbol() + strPower
      if(iUnit == -1)
      {
        antioch_unit_error("The unit \"" + curUnit + "\" is unknown. No harmonized symbol will be produced.");
        harmSymb.clear();
        return harmSymb;
      }
//updating variables:
//  - curUnit is the unit (std::string)
//  - curPower is the power (int)
      if(interUnit == "/")curPower *= -1;

// checking what we've got
      unsigned int j;
      for(j = 0; j < unitvec.size(); j++)
      {
         Units<T> tmp(unitvec[j]);
         bool same(tmp.get_symbol() == curUnit);//only contraction: strong condition
         if(!contract && !same)same = tmp.is_homogeneous(curUnit); //harmonizing: weak condition
         if(same)//if in there, update the power, 
         {
            powervec[j] += curPower;
            break;
         }
      }
      if(j >= unitvec.size())//if not, add
      {
        unitvec.push_back(curUnit);
        powervec.push_back(curPower);
      }
      interUnit = input[i];
      curUnit.clear();
    }
  }

//litre fix: only if harmonizing
//  if there is l AND m in the vector, all is converted to m. 
  if(!contract)
  {
    bool ism(false),isl(false);
    unsigned int herem(0),herel(0);
    for(unsigned int i = 0; i < unitvec.size(); i++)
    {
      if(unitvec[i] == "m")
      {
        ism = true;
        herem = i;
      }
      if(unitvec[i] == "l")
      {
        isl = true;
        herel = i;
      }
    }
    if(ism && isl)
    {
      powervec[herem] += 3 * powervec[herel];
      unitvec.erase(unitvec.begin() + herel);
      powervec.erase(powervec.begin() + herel);
    }
  }
//now writing harmsymb

  std::ostringstream outsym;
  int k(0);
  for(unsigned int i = 0; i < unitvec.size(); i++)
  {
    if(powervec[i] == 0)continue; //ignore deleted unit
    if(k != 0) // need a symbol, '.' or '/'
    {
      if(powervec[i] > 0) // '.'
      {
        outsym << ".";
      }else if(powervec[i] < 0) // '/' and reverse the power
      {
        outsym << "/";
        powervec[i] *= -1;
      }
    }
    outsym << unitvec[i];
    if(powervec[i] != 1)outsym << powervec[i];
    k++;
  }

  harmSymb = outsym.str();

  return harmSymb;
}

template <typename T>
int Units<T>::n_dimension_of_units() const
{
  int ndim(0);

  if(power.get_m()   != 0)ndim++;
  if(power.get_kg()  != 0)ndim++;
  if(power.get_s()   != 0)ndim++;
  if(power.get_A()   != 0)ndim++;
  if(power.get_K()   != 0)ndim++;
  if(power.get_mol() != 0)ndim++;
  if(power.get_cd()  != 0)ndim++;
  if(power.get_rad() != 0)ndim++;

  return ndim;
}

template <typename T>
const std::string Units<T>::get_SI_convenient_symb() const
{
  if(this->n_dimension_of_units() == 0)
  {
    std::string outStr("");
    return outStr;
  }else
  {
    for(int iUnit = 1; iUnit < UnitBaseStorage::known_units().n_known_units(); iUnit++)
    {
      if(is_homogeneous(UnitBaseStorage::known_units().stored(iUnit).symbol()) && 
         UnitBaseStorage::known_units().stored(iUnit).converter().geta() == 1. && 
         UnitBaseStorage::known_units().stored(iUnit).converter().getb() == 0.)return UnitBaseStorage::known_units().stored(iUnit).symbol();
    }
    return get_SI_symb();
  }
}

template <typename T>
void Units<T>::set_unit(const std::string &sym, std::string na)
{
  symbol = sym;
  name = na;
  toSI.clear();
  power.clear();
  fill_in_power(true);
}

template <typename T>
void Units<T>::reverse_power_symbol(std::string &subsymbol)
{
  unsigned int curInter;
  curInter = subsymbol.find(".");
  if(subsymbol.find("/") < curInter)curInter = subsymbol.find("/");
  while(curInter < subsymbol.size())
  {
   (subsymbol[curInter] == '.')?
      subsymbol.replace(curInter,1,"/"):
      subsymbol.replace(curInter,1,".");
   curInter++;
   unsigned int tmp = subsymbol.find(".",curInter);
   if(subsymbol.find("/",curInter) < tmp)tmp = subsymbol.find("/",curInter);
   curInter = tmp;
  }
}

template <typename T>
void Units<T>::symbol_to_the_power(int r,const int &key)
{
  if(n_dimension_of_units() == 0)return;
  std::string curUnit("");
  std::string tmpSymbol = contracted_symbol();
  for(unsigned int i = 0; i < tmpSymbol.size(); i++)
  {
    if(tmpSymbol[i] != '.' && tmpSymbol[i] != '/')curUnit += tmpSymbol[i];
    if(tmpSymbol[i] == '.' || tmpSymbol[i] == '/' || i == tmpSymbol.size() - 1 )
    {
      int nc(0);
      std::ostringstream po;
      int resultPower = get_integer_power(parse_power(curUnit,nc),r,key); //power
      if(resultPower == 0)
      {
        symbol = "failed";
        return;
      }
      po << resultPower;
      std::string postr = po.str();
      if(resultPower == 1)postr = "";
      if(nc != 0)
      {
        (i != tmpSymbol.size() - 1)?tmpSymbol.replace(i - nc,nc,postr):
                                    tmpSymbol.replace(i - nc + 1,nc,postr); //resulting unit
      }else
      {
        (i != tmpSymbol.size() - 1)?tmpSymbol.insert(i,postr):tmpSymbol.insert(i + 1,postr);
      }
      curUnit.clear();
      i += postr.size();
    }
  }

  symbol = tmpSymbol;
}

template <typename T>
int Units<T>::get_integer_power(int unit,int r, const int &key)
{
  if(key == 1)//multiplication
  {
     return unit * r;
  }else if(key == -1)//division
  {
    if(unit%r != 0)return 0;
    return unit / r;
  }else
  {
    std::cerr << "Key is not acceptable. This is a private method, there is a big problem..." << std::endl;
    antioch_error();
    return 0;
  }


}

template <typename T>
Converter<T> Units<T>::raise(const Converter<T> &tbm,int power) const
{
  return Converter<T>(std::pow(tbm.geta(),power),(power != 1)?0.:tbm.getb());
}

template <typename T>
Converter<T> Units<T>::raise(const Converter<T> &tbm, double power) const
{
  return Converter<T>(std::pow(tbm.geta(),power),(power != 1.)?0.:tbm.getb());
}

template <typename T>
bool Units<T>::is_in_symb(char c) const
{

  return (c != '/' && 
          c != '.' && 
          c != '-' &&
          !this->is_number(c));
}

template <typename T>
void Units<T>::showAll(std::ostream &out)
{
  out << "Unit description:"  << std::endl;
  out << "name: "             << name   << std::endl;
  out << "symbol: "           << symbol << std::endl;
  out << "SI decomposition: " << power  << std::endl;
  out << "SI converter: "     << toSI   << std::endl << std::endl;
}

template <typename T>
Units<T>& Units<T>::operator=(const Units<T> & rhs)
{
  if(this == &rhs){return *this;}
  name = rhs.get_name();
  symbol = rhs.get_symbol();
  toSI.clear();
  power.clear();
  fill_in_power(true);
  return *this;
}

template <typename T>
Units<T> & Units<T>::operator+=(const Units<T> & rhs)
{

// the name
  if(!rhs.get_name().empty())name  += " " + rhs.get_name();

// the symbol
  if(rhs.get_symbol().empty())return *this;
  if(!symbol.empty())
  {
    symbol  += ".(" + rhs.get_symbol() + ")";
  }else
  {
    symbol = rhs.get_symbol();
  }
  toSI  *= rhs.get_SI_coef();
  power += rhs.get_power();
  return *this;
}

template <typename T>
Units<T> & Units<T>::operator-=(const Units<T> & rhs)
{
  if(!rhs.get_name().empty())name += " / " + rhs.get_name();
  if(!rhs.get_symbol().empty())
  {
    symbol  += "/(" + rhs.get_symbol() + ")";
  }
  this->develop_symbol(symbol);
  toSI  /= rhs.get_SI_coef();
  power -= rhs.get_power();

  return *this;
}

template <typename T>
Units<T> & Units<T>::operator*=(int r)
{
   power *= r;
   toSI = this->raise(toSI,r);
   this->symbol_to_the_power(r,1);
   if(symbol == "failed")
   {
     symbol = this->get_SI_symb();
   }

   return *this;
}

template <typename T>
Units<T> & Units<T>::operator/=(int r)
{
  power /= r;//check consistency of root
  this->symbol_to_the_power(r,-1);
  if(symbol == "failed")
  {
     symbol = this->get_SI_symb();
  }
  toSI = this->raise(toSI,1./(T)r);

  return *this;
}


template <typename T>
Units<T> Units<T>::operator+(const Units<T> & rhs) const
{
  return (Units<T>(*this) += rhs);
}

template <typename T>
Units<T> Units<T>::operator-(const Units<T> & rhs) const
{
  return (Units<T>(*this) -= rhs);
}

template <typename T>
Units<T> Units<T>::operator*(int r) const
{
  return (Units<T>(*this) *= r);
}

template <typename T>
Units<T> Units<T>::operator/(int r) const
{
  return (Units<T>(*this) /= r);
}

  template <typename T>
  void Units<T>::print(std::ostream &os) const
  {
     os << name << " (" << symbol << "), "
        << toSI << ", "
        << power << std::endl;
  }


} //Antioch namespace

#endif