Integer ParseNonbondedFile::parse(StringParser *sp)
{
ParseComment parseComment("[#;");
ParseSigmaEpsilonCharge parseSEC;

Units *DU = Units::distanceUnits();
Units *EU = Units::energyUnits();

int NParsed = 0;
do{
sp->skipChars(" \t\n");
if(sp->isOver()) return 1;

if(expectStatement(sp,&parseComment)) continue;

if(!expectStatement(sp,&parseSEC)) return NParsed>0;


SigmaEpsilonCharge sec;

sec.sigma = parseSEC.m_sigma * DU->unit2unit("nm","Angstr");
sec.epsilon = parseSEC.m_epsilon * EU->unit2unit("kJ/mol","kcal/mol");
sec.charge = parseSEC.m_charge;
	
(*m_ff)[parseSEC.m_atomType] = sec;
NParsed++;

}while(1);

}

bool
VerticalDatum::transform(const VerticalDatum* from,
                         const VerticalDatum* to,
                         double               lat_deg,
                         double               lon_deg,
                         double&              in_out_z)
{
    if ( from == to )
        return true;

    if ( from )
    {
        in_out_z = from->msl2hae( lat_deg, lon_deg, in_out_z );
    }

    Units fromUnits = from ? from->getUnits() : Units::METERS;
    Units toUnits = to ? to->getUnits() : Units::METERS;

    in_out_z = fromUnits.convertTo(toUnits, in_out_z);

    if ( to )
    {
        in_out_z = to->hae2msl( lat_deg, lon_deg, in_out_z );
    }

    return true;
}

int PSPower::GetNumberOfUnits(PSData* pData)
{
	if (!pData)
		pData = TG.GetData();

	Units vUnits = pData->GetUnits(this);
	return vUnits.size();
}

bool
is_2D_char_2D_type(MContext *mc, DNode *t)
{
    Units *units = (dale::Units*) mc->units;
    Node *n = units->top()->dnc->toNode(t);
    n = units->top()->mp->parsePotentialMacroCall(n);
    return (n && n->token && !n->token->str_value.compare("char"));
}

void Units::MergeFrom(const Units& from) {
    GOOGLE_CHECK_NE(&from, this);
    if (from._has_bits_[0 / 32] & (0xffu << (0 % 32))) {
        if (from._has_bit(0)) {
            set_ret(from.ret());
        }
    }
    mutable_unknown_fields()->MergeFrom(from.unknown_fields());
}

void MDAtomsTyped<T>::setUnits(const Units& units,const Units& MDUnits){
  double lscale=units.getLength()/MDUnits.getLength();
  double escale=units.getEnergy()/MDUnits.getEnergy();
// scalep and scaleb are used to convert MD to plumed
  scalep=1.0/lscale;
  scaleb=1.0/lscale;
// scalef and scalev are used to convert plumed to MD
  scalef=escale/lscale;
  scalev=escale;
}

BD_convert::BD_convert (const Units& units, const symbol has, const symbol want,
                        const symbol bulk_unit)
  : in (units.get_convertion (has, bulk_unit)),
    out (units.get_convertion (Units::dry_soil_fraction (), want)),
    bulk (-42.42e42)
{ 
#if 0
  std::ostringstream tmp;
  tmp << "has = " << has << ", bulk_unit = " << bulk_unit
      << ", dsf = " << Units::dry_soil_fraction () << ", want = " << want;
  Assertion::message (tmp.str ());
#endif  
}

bool
is_2D_pointer_2D_to_2D_type(MContext *mc, DNode *t, DNode *pointee)
{
    Units *units = (dale::Units*) mc->units;
    Node *n = units->top()->dnc->toNode(t);
    n = units->top()->mp->parsePotentialMacroCall(n);
    if (!n) {
        return false;
    }

    Node *n2 = units->top()->dnc->toNode(pointee);
    n2 = units->top()->mp->parsePotentialMacroCall(n2);
    if (!n2) {
        return false;
    }

    int error_count_begin =
        units->top()->ctx->er->getErrorTypeCount(ErrorType::Error);

    Type *type  = FormTypeParse(units, n,  false, false);
    Type *type2 = FormTypeParse(units, n2, false, false);

    units->top()->ctx->er->popErrors(error_count_begin);
    if (!type || !type2) {
        return false;
    }
    return (type->points_to->isEqualTo(type2));
}

bool
types_2D_equal(MContext *mc, DNode *t1, DNode *t2)
{
    Units *units = (dale::Units*) mc->units;

    Node *n = units->top()->dnc->toNode(t1);
    n = units->top()->mp->parsePotentialMacroCall(n);
    if (!n) {
        return false;
    }

    Node *n2 = units->top()->dnc->toNode(t2);
    n2 = units->top()->mp->parsePotentialMacroCall(n2);
    if (!n2) {
        return false;
    }

    int error_count_begin =
        units->top()->ctx->er->getErrorTypeCount(ErrorType::Error);

    Type *type1 = FormTypeParse(units, n,  false, false);
    Type *type2 = FormTypeParse(units, n2, false, false);

    units->top()->ctx->er->popErrors(error_count_begin);
    if (!type1 || !type2) {
        return 0;
    }
    return (type1->isEqualTo(type2));
}

void SED::write(const QString& filename) const
{
    WavelengthGrid* lambdagrid = find<WavelengthGrid>();
    Units* units = find<Units>();

    ofstream file(filename.toLocal8Bit().constData());
    file << setprecision(8) << scientific;
    for (int ell=0; ell<lambdagrid->Nlambda(); ell++)
    {
        double lambda = lambdagrid->lambda(ell);
        double dlambda = lambdagrid->dlambda(ell);
        file << units->owavelength(lambda)
             << '\t'
             << _Lv[ell]/dlambda*lambda
             << endl;
    }
    file.close();
}

// class methods
double Units::convertValue(double value, const Units& fromUnits, const Units& toUnits) {
  if (fromUnits.isSameDimensionAs(toUnits)) {
    return (value * fromUnits._conversion / toUnits._conversion);
  }
  else {
    cerr << "Units are not dimensionally consistent" << endl;
    return 0.;
  }
}

ContourFindingBase::ContourFindingBase(const ActionOptions&ao):
Action(ao),
ActionWithInputGrid(ao),
mymin(this),
mydata(NULL)
{
  if( ingrid->noDerivatives() ) error("cannot find contours if input grid has no derivatives");
  parse("CONTOUR",contour); 
  log.printf("  calculating dividing surface along which function equals %f \n", contour);

  if( keywords.exists("FILE") ){
      std::string file; parse("FILE",file);
      if(file.length()==0 && keywords.style("FILE","compulsory") ) error("name out output file was not specified");
      else if( file.length()>0 ){
         std::string type=Tools::extension(file);
         log<<"  file name "<<file<<"\n";
         if(type!="xyz") error("can only print xyz file type with contour finding");

         fmt_xyz="%f";
         std::string precision; parse("PRECISION",precision);
         if(precision.length()>0){
           int p; Tools::convert(precision,p);
           log<<"  with precision "<<p<<"\n";
           std::string a,b;
           Tools::convert(p+5,a);
           Tools::convert(p,b);
           fmt_xyz="%"+a+"."+b+"f";
         }

         std::string unitname; parse("UNITS",unitname);
         if(unitname!="PLUMED"){
           Units myunit; myunit.setLength(unitname);
           lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
         }
         else lenunit=1.0; 

         of.link(*this); of.open(file); 

         // Now create a store data vessel to hold the points
         mydata=buildDataStashes( NULL );
      }
  }
}

	QString FrameCd::sizeDescription( const Units& units ) const
	{
		if ( units.toEnum() == Units::IN )
		{
			QString dStr = StrUtil::formatFraction( 2 * mR1.in() );

			return QString().sprintf( "%s %s %s",
						  qPrintable(dStr),
						  qPrintable(units.toTrName()),
						  qPrintable(tr("diameter")) );
		}
		else
		{
			return QString().sprintf( "%.5g %s %s",
						  2 * mR1.inUnits(units),
						  qPrintable(units.toTrName()),
						  qPrintable(tr("diameter")) );
		}
	}

bool
VCheck::Compatible::valid (const Units& units, symbol value, Treelog& msg) const
{
  if (units.can_convert (dimension, value))
    return true;
  
  std::ostringstream tmp;
  tmp << "Cannot convert [" << dimension << "] to [" << value << "]";
  msg.error (tmp.str ());
  return false;
}

void InstrumentFrame::calibrateAndWriteDataFrames(int ell, QList<Array*> farrays, QStringList fnames)
{
    Units* units = find<Units>();
    WavelengthGrid* lambdagrid = find<WavelengthGrid>();

    // conversion from bolometric luminosities (units W) to monochromatic luminosities (units W/m)
    // --> divide by delta-lambda
    double dlambda = lambdagrid->dlambda(ell);

    // correction for the area of the pixels of the images; the units are now W/m/sr
    // --> divide by area
    double xpresang = 2.0*atan(_xpres/(2.0*_distance));
    double ypresang = 2.0*atan(_ypres/(2.0*_distance));
    double area = xpresang*ypresang;

    // calibration step 3: conversion of the flux per pixel from monochromatic luminosity units (W/m/sr)
    // to flux density units (W/m3/sr) by taking into account the distance
    // --> divide by fourpid2
    double fourpid2 = 4.0*M_PI*_distance*_distance;

    // conversion from program SI units (at this moment W/m3/sr) to the correct output units
    // --> multiply by unit conversion factor
    double unitfactor = units->osurfacebrightness(lambdagrid->lambda(ell), 1.);

    // Perform the conversion, in place
    foreach (Array* farr, farrays)
    {
        (*farr) *= (unitfactor / (dlambda * area * fourpid2));
    }

    // Write a FITS file for each array
    for (int q = 0; q < farrays.size(); q++)
    {
        QString filename = _instrument->instrumentName() + "_" + fnames[q] + "_" + QString::number(ell);
        QString description = fnames[q] + " flux " + QString::number(ell);

        // Create the image and save it
        Image image(this, _Nxp, _Nyp, 1, _xpres, _ypres, "surfacebrightness");
        image.saveto(this, *(farrays[q]), filename, description);
    }
}

void DustSystem::writequality() const
{
    Log* log = find<Log>();
    Units* units = find<Units>();
    Parallel* parallel = find<ParallelFactory>()->parallel();

    // Density metric

    log->info("Calculating quality metric for the grid density...");
    DustSystemDensityCalculator calc1(this, _Nrandom, _Ncells/5);
    parallel->call(&calc1, _Nrandom);

    log->info("  Mean value of density delta: "
              + QString::number(units->omassvolumedensity(calc1.meanDelta()*1e9))
              + " nano" + units->umassvolumedensity());
    log->info("  Standard deviation of density delta: "
              + QString::number(units->omassvolumedensity(calc1.stddevDelta()*1e9))
              + " nano" + units->umassvolumedensity());

    // Optical depth metric

    log->info("Calculating quality metric for the optical depth in the grid...");
    DustSystemDepthCalculator calc2(this, _Nrandom, _Ncells/50, _Nrandom*10);
    parallel->call(&calc2, _Nrandom);

    log->info("  Mean value of optical depth delta: " + QString::number(calc2.meanDelta()));
    log->info("  Standard deviation of optical depth delta: " + QString::number(calc2.stddevDelta()));

    // Output to file

    QString filename = find<FilePaths>()->output("ds_quality.dat");
    log->info("Writing quality metrics for the grid to " + filename + "...");
    ofstream file(filename.toLocal8Bit().constData());
    file << "Mean value of density delta: "
         << units->omassvolumedensity(calc1.meanDelta()) << ' '
         << units->umassvolumedensity().toStdString() << '\n'
         << "Standard deviation of density delta: "
         << units->omassvolumedensity(calc1.stddevDelta()) << ' '
         << units->umassvolumedensity().toStdString() << '\n';
    file << "Mean value of optical depth delta: "
         << calc2.meanDelta() << '\n'
         << "Standard deviation of optical depth delta: "
         << calc2.stddevDelta() << '\n';
    file.close();
    log->info("File " + filename + " created.");
}

s16 Battle::AIAttackPosition(Arena & arena, const Unit & b, const Indexes & positions)
{
    s16 res = -1;

    if(b.isMultiCellAttack())
    {
        res = AIMaxQualityPosition(positions);
    }
    else
    if(b.isDoubleCellAttack())
    {
        Indexes results;
        results.reserve(12);

	const Units enemies(arena.GetForce(b.GetColor(), true), true);

	if(1 < enemies.size())
	{
	    for(Units::const_iterator
		it1 = enemies.begin(); it1 != enemies.end(); ++it1)
	    {
		const Indexes around = Board::GetAroundIndexes(**it1);

		for(Indexes::const_iterator
		    it2 = around.begin(); it2 != around.end(); ++it2)
		{
		    const Unit* unit = Board::GetCell(*it2)->GetUnit();
		    if(unit && enemies.end() != std::find(enemies.begin(), enemies.end(), unit))
			results.push_back(*it2);
		}
	    }

    	    if(results.size())
    	    {
        	// find passable results
        	Indexes passable = Arena::GetBoard()->GetPassableQualityPositions(b);
        	Indexes::iterator it2 = results.begin();

        	for(Indexes::const_iterator
		    it = results.begin(); it != results.end(); ++it)
            	    if(passable.end() != std::find(passable.begin(), passable.end(), *it))
                	*it2++ = *it;

        	if(it2 != results.end())
            	    results.resize(std::distance(results.begin(), it2));

        	// get max quality
        	if(results.size())
            	    res = AIMaxQualityPosition(results);
    	    }
	}
    }

    return 0 > res ? AIShortDistance(b.GetHeadIndex(), positions) : res;
}

static bool
get_type(MContext *mc, DNode *dnode, Type **type)
{
    Units *units = (dale::Units*) mc->units;
    Node *n = units->top()->dnc->toNode(dnode);
    n = units->top()->mp->parsePotentialMacroCall(n);
    if (!n) {
        return false;
    }

    int error_count_begin =
        units->top()->ctx->er->getErrorTypeCount(ErrorType::Error);

    *type = FormTypeParse(units, n, false, false);

    units->top()->ctx->er->popErrors(error_count_begin);
    if (!*type) {
        return false;
    }

    return true;
}

CString PSPower::GetAllInfo()
{
	CString strInfo, strSCs, strArmies, strFleets;
	int nArmies = 0;
	int nFleets = 0;
	SCs vSCs = TG.GetSCs(this);
	strSCs.Format(L"%d%s", vSCs.size(), (vSCs.size() != 1) ? L" SC-s:" : L" SC:");
	for (int i = 0; i < vSCs.size(); i++)
		strSCs += " " + vSCs[i]->m_pProvince->m_strName;

	Units vUnits = TG.GetUnits(this);
	for (int i = 0; i < vUnits.size(); i++)
	{
		if (vUnits[i]->GetType() == ARMY)
		{
			nArmies++;
			strArmies += " " + vUnits[i]->GetLocation()->m_strName;
		}
		else
		{
			nFleets++;
			strFleets += " " + vUnits[i]->GetLocation()->m_strName;
			CString strName = vUnits[i]->GetCoast()->m_strName;
			if (strName != "Central")
			{
				strFleets += "(" + strName + ")";
			}
		}
	}

	if (!nArmies && ! nFleets && vSCs.empty())
		strInfo.Format(L"%s%s%s", m_strName, L":\r\n", L"Eliminated\r\n\r\n");// TODO: Conditions may vary.
	else
		strInfo.Format(L"%s%s%s%s%d%s%s%s%d%s%s%s", m_strName, L":\r\n", strSCs, L"\r\n", 
			nArmies, (nArmies != 1) ? L" Armies: " : L" Army: ", strArmies, L"\r\n", 
			nFleets, (nFleets != 1) ? L" Fleets: " : L" Fleet: ", strFleets, "\r\n\r\n");

	return strInfo;
}

void PerspectiveInstrument::write()
{
    Units* units = find<Units>();
    WavelengthGrid* lambdagrid = find<WavelengthGrid>();
    int Nlambda = find<WavelengthGrid>()->Nlambda();

    // Put the data cube in a list of f-array pointers, as the sumResults function requires
    QList< Array* > farrays;
    farrays << &_ftotv;

    // Sum the flux arrays element-wise across the different processes
    sumResults(farrays);

    // Multiply each sample by lambda/dlamdba and by the constant factor 1/(4 pi s^2)
    // to obtain the surface brightness and convert to output units (such as W/m2/arcsec2)

    double front = 1. / (4.*M_PI*_s*_s);
    for (int ell=0; ell<Nlambda; ell++)
    {
        double lambda = lambdagrid->lambda(ell);
        double dlambda = lambdagrid->dlambda(ell);
        for (int i=0; i<_Nx; i++)
        {
            for (int j=0; j<_Ny; j++)
            {
                int m = i + _Nx*j + _Nx*_Ny*ell;
                _ftotv[m] = units->osurfacebrightness(lambda, _ftotv[m]*front/dlambda);
            }
        }
    }

    // Write a FITS file containing the data cube

    QString filename = _instrumentname + "_total";
    Image image(this, _Nx, _Ny, Nlambda, _s, _s, "surfacebrightness");
    image.saveto(this, _ftotv, filename, "total flux");
}