source: src/analysis_correlation.cpp@ 6b919f8

/*
 * analysis.cpp
 *
 *  Created on: Oct 13, 2009
 *      Author: heber
 */

#include <iostream>

#include "analysis_correlation.hpp"
#include "element.hpp"
#include "molecule.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "vector.hpp"


/** Calculates the pair correlation between given elements.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * \param *out output stream for debugging
 * \param *mol molecule with atoms
 * \param *type1 first element or NULL (if any element)
 * \param *type2 second element or NULL (if any element)
 * \return Map of doubles with values the pair of the two atoms.
 */
PairCorrelationMap *PairCorrelation( ofstream *out, molecule *mol, element *type1, element *type2 )
{
  PairCorrelationMap *outmap = NULL;
  double distance = 0.;

  if ((mol == NULL)) {
    cout << "No molecule given." << endl;
    return outmap;
  }
  outmap = new PairCorrelationMap;
  atom *Walker = mol->start;
  while (Walker->next != mol->end) {
    Walker = Walker->next;
    if ((type1 == NULL) || (Walker->type == type1)) {
      atom *OtherWalker = mol->start;
      while (OtherWalker->next != mol->end) { // only go up to Walker
        OtherWalker = OtherWalker->next;
        if (Walker->nr < OtherWalker->nr)
          if ((type2 == NULL) || (OtherWalker->type == type2)) {
            distance = Walker->node->Distance(OtherWalker->node);
            //cout << "Inserting " << *Walker << " and " << *OtherWalker << endl;
            outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
          }
      }
    }
  }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a point.
 * \param *out output stream for debugging
 * \param *mol molecule with atoms
 * \param *type element or NULL (if any element)
 * \param *point vector to the correlation point
 * \return Map of dobules with values as pairs of atom and the vector
 */
CorrelationToPointMap *CorrelationToPoint( ofstream *out, molecule *mol, element *type, Vector *point )
{
  CorrelationToPointMap *outmap = NULL;
  double distance = 0.;

  if ((mol == NULL)) {
    cout << "No molecule given." << endl;
    return outmap;
  }
  outmap = new CorrelationToPointMap;
  atom *Walker = mol->start;
  while (Walker->next != mol->end) {
    Walker = Walker->next;
    if ((type == NULL) || (Walker->type == type)) {
      distance = Walker->node->Distance(point);
      outmap->insert ( pair<double, pair<atom *, Vector*> >(distance, pair<atom *, Vector*> (Walker, point) ) );
    }
  }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a surface.
 * \param *out output stream for debugging
 * \param *mol molecule with atoms
 * \param *type element or NULL (if any element)
 * \param *Surface pointer to Tesselation class surface
 * \param *LC LinkedCell structure to quickly find neighbouring atoms
 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
 */
CorrelationToSurfaceMap *CorrelationToSurface( ofstream *out, molecule *mol, element *type, Tesselation *Surface, LinkedCell *LC )
{

  CorrelationToSurfaceMap *outmap = NULL;
  double distance = 0.;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;

  if ((Surface == NULL) || (LC == NULL) || (mol == NULL)) {
    cout << "No Tesselation, no LinkedCell or no molecule given." << endl;
    return outmap;
  }
  outmap = new CorrelationToSurfaceMap;
  atom *Walker = mol->start;
  while (Walker->next != mol->end) {
    Walker = Walker->next;
    if ((type == NULL) || (Walker->type == type)) {
      triangle = Surface->FindClosestTriangleToPoint(out, Walker->node, LC );
      if (triangle != NULL) {
        distance = DistanceToTrianglePlane(out, Walker->node, triangle);
        outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
      }
    }
  }

  return outmap;
};

/** Returns the start of the bin for a given value.
 * \param value value whose bin to look for
 * \param BinWidth width of bin
 * \param BinStart first bin
 */
double GetBin ( double value, double BinWidth, double BinStart )
{
  double bin =(double) (floor((value - BinStart)/BinWidth));
  return (bin*BinWidth+BinStart);
};


/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelation( ofstream *file, BinPairMap *map )
{
  *file << "# BinStart\tCount" << endl;
  for (BinPairMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << runner->first << "\t" << runner->second << endl;
  }
};

/** Prints correlation (double, (atom*,atom*) ) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputPairCorrelation( ofstream *file, PairCorrelationMap *map )
{
  *file << "# BinStart\tAtom1\tAtom2" << endl;
  for (PairCorrelationMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
  }
};

/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelationToPoint( ofstream *file, CorrelationToPointMap *map )
{
  *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
  for (CorrelationToPointMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << runner->first;
    for (int i=0;i<NDIM;i++)
      *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
    *file << endl;
  }
};

/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelationToSurface( ofstream *file, CorrelationToSurfaceMap *map )
{
  *file << "# BinStart\tTriangle" << endl;
  for (CorrelationToSurfaceMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << runner->first << "\t" << *(runner->second.second) << endl;
  }
};