UnitIO.h

// -*- C++ -*-
//
// UnitIO.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_UnitIO_H
#define ThePEG_UnitIO_H
// This is the declaration of the IUnit and OUnit classes and
// associated templated functions.

#include <complex>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <cmath>

namespace ThePEG {

template <typename T, typename UT>
struct OUnit {

  OUnit(const T & t, const UT & u): theX(t), theUnit(u) {}

  OUnit(const OUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}

  const T & theX;

  const UT & theUnit;
};

template <typename T, typename UT>
struct IUnit {

  IUnit(T & t, const UT & u): theX(t), theUnit(u) {}

  IUnit(const IUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}

  T & theX;

  const UT & theUnit;

};

template <typename T, typename UT>
inline OUnit<T,UT> ounit(const T & t, const UT & ut) {
  return OUnit<T,UT>(t, ut);
}

template <typename T, typename UT>
inline IUnit<T,UT> iunit(T & t, const UT & ut) {
  return IUnit<T,UT>(t, ut);
}

template <typename OStream, typename T, typename UT>
void ounitstream(OStream & os, const T & t, UT & u) {
  os << t/u;
}

template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, T & t, UT & u) {
  double d;
  is >> d;
  t = d*u;;
}

template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, std::complex<T> & t, UT & u) {
  std::complex<double> d;
  is >> d;
  t = d*u;;
}

template <typename OStream, typename T, typename UT>
OStream & operator<<(OStream & os, const OUnit<T,UT> & u) {
  ounitstream(os, u.theX, u.theUnit);
  return os;
}

template <typename IStream, typename T, typename UT>
IStream & operator>>(IStream & is, const IUnit<T,UT> & u) {
  iunitstream(is, u.theX, u.theUnit);
  return is;
}

template <typename T, typename UT>
struct OUnitErr {

  OUnitErr(const T & t, const T & dt, const UT & u): x(t/u), dx(dt/u) {}

  double x;

  double dx;
  
};

template <typename T, typename UT>
inline OUnitErr<T,UT> ouniterr(const T & t, const T & dt, const UT & ut) {
  return OUnitErr<T,UT>(t, dt, ut);
}

inline OUnitErr<double,double> ouniterr(double t, double dt) {
  return OUnitErr<double,double>(t, dt, 1.0);
}

template <typename OStream, typename T, typename UT>
OStream & operator<<(OStream & os, const OUnitErr<T,UT> & u) {
  if ( ! isfinite(u.x) ) return os << u.x;
  if ( ! isfinite(u.dx) ) {
    ostringstream out;
    out << u.x << '(' << u.dx << ')';
    return os << out.str();
  }
  double dx = min(u.dx, abs(u.x));
  if ( dx <= 0.0 ) return os << u.x;
  double x = abs(u.x);
  ostringstream osse;
  osse << std::scientific << setprecision(0) << dx;
  string sse = osse.str();
  string::size_type ee = sse.find('e');
  long m = static_cast<long>(round(abs(x)/std::pow(10.0,std::atoi(sse.substr(ee + 1).c_str()))));
  int powx = m <= 0? os.precision(): int(log10(double(m)));
  if ( m <= 0 || powx > os.precision() ) sse[0]='0';  
  ostringstream oss;
  oss << std::scientific << setprecision(powx) << x;
  string ss = oss.str();
  string::size_type e = ss.find('e');
  ostringstream out;
  int pp = std::atoi(ss.substr(e + 1).c_str());
  if ( pp%3 == 0 )
    out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
  else if ( (pp - 1)%3 == 0 ) {
    ostringstream oss;
    oss << std::scientific << setprecision(powx) << x/10.0;
    string ss = oss.str();
    string::size_type e = ss.find('e');
    if ( powx == 0 )
      out << ss.substr(0, e) << "0(" << sse[0] << "0)" << ss.substr(e);
    else if ( powx == 1 )
      out << ss.substr(0, ss.find('.'))
          << ss.substr(ss.find('.') + 1, e - ss.find('.') - 1)
          << "(" << sse[0] << ")" << ss.substr(e);
    else {
      swap(ss[ss.find('.')], ss[ss.find('.') + 1]);
      out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
    }
  }
  else {
    ostringstream oss;
    oss << std::scientific << setprecision(powx) << x*10.0;
    string ss = oss.str();
    string::size_type e = ss.find('e');
    if ( powx == 0 )
      out << "0." << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
    else {
      swap(ss[ss.find('.')], ss[ss.find('.') - 1]);
      out << ss.substr(0, ss.find('.')) << "0" << ss.substr(ss.find('.'), e)
          << "(" << sse[0] << ")" << ss.substr(e);
    }
  }
  string res = out.str();
  if ( u.x < 0.0 )
    res = "-" + res;
  return os << res;
}

template <typename T, typename UT>
struct IUnitErr {

  IUnitErr(T & t, T & dt, const UT & u): x(t), dx(dt), ut(u) {}

  T & x;

  T & dx;
  
  UT ut;

};

template <typename T, typename UT>
inline IUnitErr<T,UT> iuniterr(T & t, T & dt, const UT & ut) {
  return IUnitErr<T,UT>(t, dt, ut);
}

inline IUnitErr<double,double> iuniterr(double & t, double & dt) {
  return IUnitErr<double,double>(t, dt, 1.0);
}

template <typename IStream, typename T, typename UT>
IStream & operator>>(IStream & is, const IUnitErr<T,UT> & u) {
  string s;
  double x = 0.0;
  double dx = 0.0;
  double ex = 1.0;
  is >> s;
  string::size_type open = s.find('(');
  string::size_type close = s.find(')');
  string se = "0";
  string sp = "1";
  double pe = 1.0;
  if ( open != string::npos && close != string::npos ) {
    se = s.substr(open + 1);
    sp += s.substr(close + 1);
    string::size_type dot = s.find('.');
    if ( dot != string::npos && dot < open ) pe = std::pow(10.0, 1.0 - (open - dot));
  }

  istringstream(s) >> x;
  istringstream(se) >> dx;
  istringstream(sp) >> ex;

  u.x = x*ex*u.ut;
  u.dx = dx*ex*pe*u.ut;

  return is;
}

}

#endif /* ThePEG_UnitIO_H */