//
// File: mixture.cpp
// Purpose: Represent a material mixture
// Created: 27-OCT-1998 by Stephen Kahn
//

#include "material/Mixture.hpp"
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

using namespace D0Material;

#ifndef _WIN32
using std::abs;
#endif

Mixture::Mixture(
  const std::string& name,
  const std::string& symbol,
  const std::vector<Material> elements, 
  const std::vector<double> weights,
  const bool weights_are_fractions /* true */
)
: Material(name, symbol, 1.0, 1.0, 1.0, 1.0)
, _weights(weights)
, _elements(elements)
, _nlmat(elements.size())
, _wtot(0.0)
, _wtype(!weights_are_fractions)
{
  int i;

  for (i = 0; i < elements.size(); i++)
  {
    _wtot += weights[i];
  }

  // Initialize the attributes in Material
  calc_eff_material();
}


Mixture::Mixture(
  const std::string &name,
  const std::string &symbol
)
: Material(name, symbol, 1.0, 1.0, 1.0, 1.0)
, _nlmat(0)
, _wtot(0.0)
{
  _weights.clear();
  _elements.clear(); 
}  

// empty constructor.  Is this smart?  
// Use dummy values to fill
// elements required in Material.  The material here is undefined until
// calc_eff_material is called.

void
Mixture::add_material_by_weight(const double wt, const Material& mat)
{
  /**  A material element is added with wt being the fraction by weight */
  _wtot += wt;
  _weights.push_back(wt);
  _elements.push_back(mat);
  _nlmat++;
  _wtype = true;
}

void  Mixture::add_material_by_fraction(const double wt, const Material& mat)
{
  /** A material element is added with wt being the fraction by number */
  _wtot += wt;
  _weights.push_back(wt);
  _elements.push_back(mat);
  _nlmat++;
  _wtype = false;
}

void
Mixture::calc_eff_material()
{
  int i;
  double aeff = 0.0;
  double zeff = 0.0;    
  double xinv = 0.0;
  double dedx = 0.0;
  double dens = 0.0;
  const double eps = 0.001;
  int nlmat = (int) _elements.size();

  if (nlmat != _nlmat) 
  {
    std::cerr
      << " Mixture : inconsistent number of elements for "
      << d_name << nlmat << _nlmat
      << std::endl;
  }

  if (_wtype) // fraction by weight
  {
    if (abs(_wtot - 1.0) > eps) // renormalize
    {
      std::vector<double>::iterator wti = _weights.begin();
      while (wti != _weights.end())
      {
        double wtel = (double) *wti;
        wtel /= _wtot;
        wti = _weights.erase(wti);
        _weights.insert(wti, wtel);
        wti++;
      }
    _wtot = 1.0;
    }
  }
  else // fraction by number
  {
    double atot = 0.0;
    double wtot = 0.0;
    for (i = 0; i < nlmat; i++) //  calc atot
    {
      double wt = _weights[i];
      atot += wt*_elements[i].getA();
      wtot += wt;
    }
    std::vector<double> new_weights;
    for (i = 0; i < nlmat; i++)
    {
      double wt = _weights[i];
      wt *= _elements[i].getA() / atot;
      wt *= wtot;
      new_weights.push_back(wt);
    }
    _weights = new_weights;
  }
         
  for (i = 0; i < nlmat; i++)
  {
    Material m = _elements[i]; 
    aeff += _weights[i]*m.getA();
    zeff += _weights[i]*m.getZ();
    dens += _weights[i]*m.getDensity();
    double x0 = m.getX0_gm();
    if (x0 > 0.0)
    {
      xinv += _weights[i] / x0;
    }
    else
    {
      std::cerr << " Mixture : error calculating rad length for "
                << m.getName()
                << std::endl;
      break;
    }
  }

  d_zeff = zeff;
  d_aeff = aeff;
  d_rho = dens;
  d_x0 = 1.0 / xinv;
}

double
Mixture::atomic_fraction(int i) const
{
  double wt;
  if (!_wtype) //   weight stored by atomic fraction
  {
    wt = _weights[i];
  }
  else //   weight stored by weight fraction
  {
    wt = _weights[i] / _elements[i].getA();
    double wtot = 0.0;
    int nelements = _weights.size();
    int j;
    for (j = 0; j < nelements; j++)
    {
      wtot += _weights[j] / _elements[j].getA();
    }
    wt *= 1.0 / wtot;
  }
  return wt;
}

double
Mixture::weight_fraction(int i) const
{
  double  wt;
  if (!_wtype) //  weight stored by atomic fraction
  {
    wt = _weights[i] * _elements[i].getA();
    double atot = 0.0;
    int nelements = _weights.size();
    int j;
    for (j = 0; j < nelements; j++)
    {
      atot += _weights[i] * _elements[j].getA();
    }
    wt *= 1.0 / atot;
   }
   else
   {
     wt = _weights[i];
    }
   return wt;
}

std::ostream&
D0Material::operator<<(std::ostream& os, const Mixture& me) 
{
  os
    << "Mixture "
    << me.getName()
    << " consists of "
    << me.numb_elements()
    << " elements"
    << std::endl;

  os.setf(std::ios::fixed, std::ios::floatfield);
  int oldprec = os.precision();

  os
    << "Mixture: " << me.getName() 
    << " (" << me.getSymbol() << ")"
    << ", Z = " << std::setprecision(2) << me.getZ() 
    << ", A = " << std::setprecision(2) << me.getA()
    << ", density = " << std::setprecision(2) << me.getDensity() << " g/cm**3"
    << ", X0 = " << std::setprecision(2) << me.getX0_gm() << " g/cm**2"
    << std::endl;
      
  os.unsetf(std::ios::fixed);
  os.precision(oldprec);
    
  int i;
  for (i = 0; i < me.numb_elements(); i++)
  {
    double wt = me.atomic_fraction(i);
    double wt1 = me.weight_fraction(i);
    Material mat = me.get_material(i);
    os
      << "Fraction " << mat.getName()
      << " atom frac= " << wt 
      << ", wt_frac= " << wt1
      << std::endl;
  }

  return os;
}