simulations contra the terroristic nimbus of the second law

On Mar 2, 11:29 pm, galath...@veawb.coop (galathaea) wrote:


still being very naive about this whole crackpot / crank thing
  i accidentally let the engineer inside think too hard about this

since my degree is in physics simulations
  and my career is in programming
  and i happened to have a simulation generator i have been building
i generated a simulation for a gas in a column with gravity

the simulation generator is my attempt at building
  a commercially viable product at home
and although it still needs much work
and doesn't generate the prettiest code yet
  i couldn't help seeing what it would give me here

i have run some simulations
  playing around with the specifications
and all my tests have validated loschmidt so far

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
THE TEMPERATURE AT THE BOTTOM IS HOTTER THAN AT TOP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

of course this seemed obvious to me
  because gravity accelerates particles moving down the column
  and slows them moving up
but there is nothing better than seeing it validated in a simulation

now
  currently particle collisions are very rare
    ( due to size of particles and iteration delta )
  and i haven't fully tested that handling yet
so the values are not mixing the z components with the x,y

collisions with the wall only affect x,y
collisions with the floor or ceiling only affect z

but i will play with that some more

i have hand entered some comments
and added output of the simulation
  ( the generator creates blanket output for the specification
    but litle descriptive output currently )

anyways
  i cannot stress that i am well aware
  of many of the refactorings needed in this code
as the generator placed code in places i normally wouldn't

also since it was working with a larger model
it has some stuff that can be completely refactored out

despite that
here is my generated loschmidt simulation
  ( in c++, from a generator written in OCaml )
in one large file format (the only working generation right now)

this may get mercilessly chopped by google (sorry)

-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

#include <iostream>
#include <vector>
#include <cmath>
#include <string>
#include <sstream>
#include <exception>

#include <boost/shared_ptr.hpp>
#include <boost/random.hpp>

namespace
{
  // math
  double pi(3.14159265358);

  // particle descriptors
  unsigned numberOfParticles(1000);
  double particleMass(1.67372e-24); // (kg)
  double particleRadius(2.4e-11); // (m)

  // environment
  double columnRadius(20.); // (m)
  double columnHeight(100.); // (m)

  // energetics
  double temperature(273.15); // (K) = 0(C)
  double boltzmann_k(1.3806505e-23); // (J/K)
  double sigma(std::sqrt(boltzmann_k * temperature /
particleMass)); // ~ 47.468

  // gravity
  double g(-9.807); // (m/s)
}

double distance(double x, double y)
{
  return std::sqrt(x * x + y * y);
}

double distance(double x, double y, double z)
{
  return std::sqrt(x * x + y * y + z * z);
}

struct Point3D
{
  Point3D() :
    x_(),
    y_(),
    z_()
  {}

  Point3D(double x, double y, double z) :
    x_(x),
    y_(y),
    z_(z)
  {}

  double x_;
  double y_;
  double z_;
};

struct Velocity3D : Point3D
{
  Velocity3D()
  {}

  Velocity3D(double x, double y, double z) :
    Point3D(x, y, z)
  {}
};

enum CollisionType
{
  column,
  particle
};

class Existent
{
public:
  Existent(Point3D const& initialPosition,
           Velocity3D const& initialVelocity,
           double mass,
           std::string const& name) :
    position_(initialPosition),
    velocity_(initialVelocity),
    mass_(mass),
    name_(name)
  {}

  double x() const
  {
    return position_.x_;
  }

  void x(double x)
  {
    position_.x_ = x;
  }

  double y() const
  {
    return position_.y_;
  }

  void y(double y)
  {
    position_.y_ = y;
  }

  double z() const
  {
    return position_.z_;
  }

  void z(double z)
  {
    position_.z_ = z;
  }

  double vx() const
  {
    return velocity_.x_;
  }

  void vx(double vx)
  {
    velocity_.x_ = vx;
  }

  double vy() const
  {
    return velocity_.y_;
  }

  void vy(double vy)
  {
    velocity_.y_ = vy;
  }

  double vz() const
  {
    return velocity_.z_;
  }

  void vz(double vz)
  {
    velocity_.z_ = vz;
  }

  double m() const
  {
    return mass_;
  }

  std::string name() const
  {
    return name_;
  }

  virtual bool isCollided(boost::shared_ptr<Existent> otherEntity)
const = 0;
  virtual CollisionType collisionType() const = 0;

private:
  Point3D position_;
  Velocity3D velocity_;

  double mass_;

  std::string name_;
};

class Particle : public Existent
{
public:
  Particle(
      Point3D const& initialPosition,
      Velocity3D const& initialVelocity,
      double radius,
      double mass,
      std::string const& name) :
    Existent(initialPosition, initialVelocity, mass, name),
    radius_(radius)
  {}

  double radius() const
  {
    return radius_;
  }

  virtual bool isCollided(boost::shared_ptr<Existent> otherEntity)
const
  {
    // a particle can only be a first entity if the second is also a
particle
    if (distance(x() - otherEntity->x(), y() - otherEntity->y(), z() -
otherEntity->z()) <
        (radius() + boost::static_pointer_cast<Particle>(otherEntity)-


    {
      if (name() == "0")
        std::cout << "particle on particle action" << std::endl;
      return true;
    }
    return false;
  }

  virtual CollisionType collisionType() const
  {
    return particle;
  }

private:
  double radius_;
};

class EnclosingColumn : public Existent
{
public:
  EnclosingColumn(double radius, double height) :
    Existent(Point3D(), Velocity3D(), 0., "column"),
    radius_(radius),
    height_(height)
  {}

  virtual bool isCollided(boost::shared_ptr<Existent> otherEntity)
const
  {
    // otherEntity must always be a particle!
    if (distance(otherEntity->x(), otherEntity->y()) > (radius_ -
boost::static_pointer_cast<Particle>(otherEntity)->radius()))
    {
      if (otherEntity->name() == "0")
        std::cout << "collision with wall of column" << std::endl;
      return true;
    }
    else if (otherEntity->z() > (height_ -
boost::static_pointer_cast<Particle>(otherEntity)->radius()))
    {
      if (otherEntity->name() == "0")
        std::cout << "collision with ceiling" << std::endl;
      return true;
    }
    else if (otherEntity->z() <
boost::static_pointer_cast<Particle>(otherEntity)->radius())
    {
      if (otherEntity->name() == "0")
        std::cout << "collision with floor" << std::endl;
      return true;
    }
    return false;
  }

  virtual CollisionType collisionType() const
  {
    return column;
  }

private:
  double radius_;
  double height_;
};

class PhysicalWorld
{
public:
  typedef std::vector<boost::shared_ptr<Existent> > Existents;

  Existents::iterator begin()
  {
    return existents_.begin();
  }

  Existents::iterator end()
  {
    return existents_.end();
  }

  void pushBack(boost::shared_ptr<Existent> existent)
  {
    existents_.push_back(existent);
  }

private:
  Existents existents_;
};

struct StrataProfile
{
  StrataProfile() :
    temperatureAccumulatorX_(),
    temperatureAccumulatorY_(),
    temperatureAccumulatorZ_(),
    count_()
  {}

  void add(double tx, double ty, double tz)
  {
    temperatureAccumulatorX_ += tx;
    temperatureAccumulatorY_ += ty;
    temperatureAccumulatorZ_ += tz;

    ++count_;
  }

  double temperatureAccumulatorX_;
  double temperatureAccumulatorY_;
  double temperatureAccumulatorZ_;

  unsigned count_;
};

class PhysicsEngine
{
public:
  PhysicsEngine(PhysicalWorld & world, double final, double delta) :
    world_(world),
    time_(0.),
    final_(final),
    delta_(delta)
  {}

  double time() const
  { return time_; }

  double delta() const
  { return delta_; }

  void cycle()
  {
    time_ += delta_;
  }

  void run()
  {
    std::cout << "entering the run loop" << std::endl;

    // full run statistics
    StrataProfile top;
    StrataProfile bottom;

    // main event loop
    while (time_ < final_)
    {
      std::cout << "^~^~^~^~^~^~^~^ loop time " << time_ << "
^~^~^~^~^~^~^~^" << std::endl;

      // check for collisions first as initial adjustment of velocity
      for (PhysicalWorld::Existents::iterator
initialIntersector(world_.begin());
           initialIntersector != world_.end();
           ++initialIntersector)
      {
        if ((initialIntersector + 1) != world_.end())
        {
          //std::cout << "-" << std::flush;

          for (PhysicalWorld::Existents::iterator
secondIntersector(initialIntersector + 1);
               secondIntersector != world_.end();
               ++secondIntersector)
          {
            //std::cout << "+" << std::flush;

            if (collision(initialIntersector, secondIntersector))
            {
              //std::cout << "collision" << std::endl;
              adjustFromCollision(initialIntersector,
secondIntersector);
            }
          }
        }
      }

      // then cycle through each and adjust their final positions and
force changes to velocity
      for (PhysicalWorld::Existents::iterator
existent(world_.begin());
            existent != world_.end();
            ++existent)
      {
        adjustFreeMoving(existent);
      }

      // and cycle the time
      cycle();

      // now we can snapshot the new state
      // follow a particle
      PhysicalWorld::Existents::iterator existent(world_.begin() + 1);

      std::cout << "particle 1's iterated state:"
        << "\n x: " << (*existent)->x()
        << "\n y: " << (*existent)->y()
        << "\n z: " << (*existent)->z()
        << "\n vx: " << (*existent)->vx()
        << "\n vy: " << (*existent)->vy()
        << "\n vz: " << (*existent)->vz() << std::endl;

      // and build temperature profile
      unsigned strata(5);
      std::vector<StrataProfile> temperatures(strata);

      while (existent != world_.end())
      {
        unsigned particularStrata(((*existent)->z() * 5.0) /
columnHeight);
        if (particularStrata >= strata)
          particularStrata = strata - 1;
        if (particularStrata < 0)
          particularStrata = 0;

        double tx((*existent)->m() * (*existent)->vx() * (*existent)-


        double ty((*existent)->m() * (*existent)->vy() * (*existent)-


        double tz((*existent)->m() * (*existent)->vz() * (*existent)-


        temperatures[particularStrata].add(tx, ty, tz);

        ++existent;
      }

      for (unsigned stratum(0); stratum < strata; ++stratum)
      {
        std::cout << "in stratum " << stratum << " there were " <<
temperatures[stratum].count_ << " particles" << std::endl;

        if (temperatures[stratum].count_)
        {
          double tx(temperatures[stratum].temperatureAccumulatorX_ /
temperatures[stratum].count_);
          double ty(temperatures[stratum].temperatureAccumulatorY_ /
temperatures[stratum].count_);
          double tz(temperatures[stratum].temperatureAccumulatorZ_ /
temperatures[stratum].count_);

          std::cout << " tx: " << tx
            << " ty: " << ty
            << " tz: " << tz
            << " Tavg: " << ((tx + ty + tz) / 3.) << std::endl;

          if (0 == stratum)
            bottom.add(tx, ty, tz);
          else if ((strata - 1) == stratum)
            top.add(tx, ty, tz);
        }
        else std::cout << "no temperature information possible" <<
std::endl;
      }
    }

    // dump full run statistics
    double txTop(top.temperatureAccumulatorX_ / top.count_);
    double tyTop(top.temperatureAccumulatorY_ / top.count_);
    double tzTop(top.temperatureAccumulatorZ_ / top.count_);
    double txBottom(bottom.temperatureAccumulatorX_ / bottom.count_);
    double tyBottom(bottom.temperatureAccumulatorY_ / bottom.count_);
    double tzBottom(bottom.temperatureAccumulatorZ_ / bottom.count_);

    std::cout << "over the simulation we have"
      << "\ntop strata of the column - tx: " << txTop << " ty: " <<
tyTop << " tz: " << tzTop << " Tavg: " << ((txTop + tyTop + tzTop) /
3.)
      << "\nbottom strata of the column - tx: " << txBottom << " ty: "
<< tyBottom << " tz: " << tzBottom << " Tavg: " << ((txBottom +
tyBottom + tzBottom) / 3.)
      << std::endl;
  }

private:
  bool collision(PhysicalWorld::Existents::iterator object1,
PhysicalWorld::Existents::iterator object2)
  {
    return (*object1)->isCollided(*object2);
  }

  virtual void adjustFromCollision(PhysicalWorld::Existents::iterator
object1, PhysicalWorld::Existents::iterator object2) = 0;
  virtual void adjustFreeMoving(PhysicalWorld::Existents::iterator
object) = 0;

  PhysicalWorld & world_;

  double time_;
  double final_;
  double const delta_;
};

class LinearGravityEngine : public PhysicsEngine
{
public:
  LinearGravityEngine(PhysicalWorld & world, double final, double
delta) :
    PhysicsEngine(world, final, delta)
  {}

private:
  virtual void adjustFromCollision(PhysicalWorld::Existents::iterator
object1, PhysicalWorld::Existents::iterator object2)
  {
    if ((*object1)->collisionType() == particle)
    {
      // calculate relative distances and speed
      double deltaX((*object1)->x() - (*object2)->x());
      double deltaY((*object1)->y() - (*object2)->y());
      double deltaZ((*object1)->z() - (*object2)->z());
      double deltaVX((*object1)->vx() - (*object2)->vx());
      double deltaVY((*object1)->vy() - (*object2)->vy());
      double deltaVZ((*object1)->vz() - (*object2)->vz());

      double relativeDistance(std::sqrt(deltaX * deltaX + deltaY *
deltaY + deltaZ * deltaZ));
      double relativeVelocity(std::sqrt(deltaVX * deltaVX + deltaVY *
deltaVY + deltaVZ * deltaVZ));

      // find relative angle
      double theta(std::acos(deltaZ / relativeDistance));
      double phi((deltaX == 0 && deltaY == 0) ? 0. :
std::atan2(deltaY, deltaX));

      // ch-ch-ch-changes
      double mu((*object2)->m() / (*object1)->m());
      double diffZ(2. * (deltaVZ + std::tan(theta) *
                         (std::cos(phi) * deltaVX + std::sin(phi) *
deltaVY)) /
                   (1 + std::tan(theta) * std::tan(theta)) * (1 +
mu));

      (*object2)->vz((*object2)->vz() + diffZ);
      (*object2)->vx((*object2)->vx() + std::tan(theta) *
std::cos(phi) * diffZ);
      (*object2)->vy((*object2)->vy() + std::tan(theta) *
std::sin(phi) * diffZ);
      (*object1)->vz((*object1)->vz() - mu * diffZ);
      (*object1)->vx((*object1)->vx() - std::tan(theta) *
std::cos(phi) * mu * diffZ);
      (*object1)->vy((*object1)->vy() - std::tan(theta) *
std::sin(phi) * mu * diffZ);
    }
    else if ((*object1)->collisionType() == column)
    {
      if ((*object2)->z() < 0.)
      {
        // reflect off bottom
        (*object2)->vz(-(*object2)->vz());
      }
      else if ((*object2)->z() > columnHeight)
      {
        // reflect off top
        (*object2)->vz(-(*object2)->vz());
      }

      if (distance((*object2)->x(), (*object2)->y()))
      {
        // reflect off sides
        // calculate the collision theta
        double theta(((*object2)->x() == 0 && (*object2)->y() == 0) ?
0. : std::atan2((*object2)->y(), (*object2)->x()));
        double velocityTheta(((*object2)->vx() == 0 && (*object2)-


        double velocity(distance((*object2)->vx(), (*object2)->vy()));

        (*object2)->vx(velocity * std::cos(pi + 2 * theta -
velocityTheta));
        (*object2)->vy(velocity * std::sin(pi + 2 * theta -
velocityTheta));
      }
    }
  }

  virtual void adjustFreeMoving(PhysicalWorld::Existents::iterator
object)
  {
    // simple linear free with vertical gravity
    (*object)->x((*object)->x() + (*object)->vx() * delta());
    (*object)->y((*object)->y() + (*object)->vy() * delta());

    (*object)->z((*object)->z() + (*object)->vz() * delta() + (1./2.)
* g * delta() * delta());
    (*object)->vz((*object)->vz() + g * delta());
  }
};

int main()
{
  try
  {
    std::cout << "Functional Simulations Engine 0.39 - starting" <<
std::endl;
    std::cout << "generated from source loschmidt.ml" << std::endl;

    // build up the world
    PhysicalWorld world;

    std::cout << "created the world" << std::endl;

    // add the column (radius 1
    boost::shared_ptr<Existent> column(new EnclosingColumn(columnRadius,
columnHeight));
    world.pushBack(column);

    std::cout << "added the column" << std::endl;

    // rng for simulations
    boost::mt19937 rng;
    rng.seed(static_cast<unsigned int>(std::time(0)));

    // create particles
    for (unsigned particleCounter(0); particleCounter <
numberOfParticles; ++particleCounter)
    {
    if (! particleCounter)
      std::cout << "creating particle 1" << std::endl;

    // find a random position in the column
    // height in [0, columnHeight]
    boost::uniform_real<> possibleHeight(0., columnHeight);
    boost::variate_generator<boost::mt19937 &, boost::uniform_real<> >
heightGenerator(rng, possibleHeight);

    double particleHeight(heightGenerator());

    double particleX(0.);
    double particleY(0.);

    // find an x, y position in the column
    // just loop until valid in a circle (to maintain distribution)
    do
    {
      // build generator
      boost::uniform_real<> possibleXY(0., columnRadius);
      boost::variate_generator<boost::mt19937 &, boost::uniform_real<> >
tranverseGenerator(rng, possibleXY);

      // pull an x, y
      particleX = tranverseGenerator();
      particleY = tranverseGenerator();
    }
    while (distance(particleX, particleY) >= columnRadius);

    // now get a velocity in the maxwellian-boltzmann distribution for
the temperature
    // which is just a normal distribution
    // with mean = 0
    // variance sigma^2 = k T / m
    boost::normal_distribution<> possibleVelocity(0., sigma);
    boost::variate_generator<boost::mt19937 &,
boost::normal_distribution<> > velocityGenerator(rng,
possibleVelocity);

    double velocityX(velocityGenerator());
    double velocityY(velocityGenerator());
    double velocityZ(velocityGenerator());

    if (! particleCounter)
      std::cout << "initial state position: "
        << "\n x: " << particleX
        << " y: " << particleY
        << " z: " << particleHeight
        << "\n vx: " << velocityX
        << " vy: " << velocityY
        << " vz: " << velocityZ << std::endl;

    std::stringstream converter;
    converter << particleCounter;

    // now build the particle
    boost::shared_ptr<Existent> particle(new Particle(
      Point3D(particleX, particleY, particleHeight),
      Velocity3D(velocityX, velocityY, velocityZ),
      particleRadius,
      particleMass,
      converter.str()));

    world.pushBack(particle);
    }

    std::cout << "all particles created!" << std::endl;

    // now that we have the world built
    // build the physics engine
    // attach the world
    // and configure to run x deltas
    LinearGravityEngine engine(world, 5., 0.1);
    engine.run();
  }
  catch (std::exception & blowUp)
  {
    std::cout << "exception at lowest level: " << blowUp.what() <<
std::endl;
  }
  catch (...)
  {
    std::cout << "blowUps happen" << std::endl;
  }
}

-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

i apologise if my hand editing inserted any tabs or other errata

you will need to have a download of boost

on a UN*X/Linux/MacOSX machine simply type:
g++ -I<boost location> main.cpp

where <boost location> is the path to the outermost boost folder
  from the distribution

enjoy the crankdom!

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
galathaea: prankster, fablist, magician, liar