this may be
in secret
constructivisms greatest defense
against the defeatism of ultrafinitism
-+-+-+-
i performed the checks on particle-particlelast night
and noticed that when collisions occurred
there was a boost in the energy of the system
it turns out that the spherical collision equations
it uses as reference
are in a canonical basis in relation to the particles
but that i only coded it to boost one of the velocities to the frame
i've changed it to now properly perform the coordinate transformations
completely in both directions
i generated the program
asking it to also make a few code motions useful to debugging
to illustrate the usefulness and speed of generation
i have verified that the energetics of collisions conserve energy
and i have tracked several collisions
and verified realistic final conditions
( angles work out to maintain momentum )
now
velocities are mixing properly
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
AND THE TEMPERATURE STRATIFICATION REMAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
i have tested radii of particles up to 1m
and a variety of deltas as low as .0001s
i have extended the times of the simulation
everything currently still supports loschmidt
here is the dump of the current generation:
-+-+-+-+-+-+-+-+-+-+-+-+-++++++...
#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
double particleMass(1.67372e-24); // (kg)
//double particleRadius(2.4e-11); // (m)
double particleRadius(.5); // (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)
// simulation controls
unsigned numberOfParticles(2000);
double simulationTime(5.);
double simulationDelta(.01);
// debug
bool collisionData(true);
bool stateData(true);
bool temperatureStratumData(true);
bool particleCollisionCalculations(true);
}
template <typename Metric>
Metric distance(Metric const& x, Metric const& y)
{
return std::sqrt(x * x + y * y);
}
template <typename Metric>
Metric distance(Metric const& x, Metric const& y, Metric const& 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)->radius()))
{
if (collisionData && (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 (collisionData && (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 (collisionData && (otherEntity->name() == "0"))
std::cout << "collision with ceiling" << std::endl;
return true;
}
else if (otherEntity->z() <
boost::static_pointer_cast<Particle>(otherEntity)->radius())
{
if (collisionData && (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;
squareAccumulatorX_ += tx * tx;
squareAccumulatorY_ += ty * ty;
squareAccumulatorZ_ += tz * tz;
++count_;
}
double tx() const
{
return temperatureAccumulatorX_ / count_;
}
double ty() const
{
return temperatureAccumulatorY_ / count_;
}
double tz() const
{
return temperatureAccumulatorZ_ / count_;
}
double mean() const
{
return (tx() + ty() + tz()) / 3.;
}
double standardDeviation() const
{
return std::sqrt(((squareAccumulatorX_ + squareAccumulatorY_ +
squareAccumulatorZ_) / (3. * count_)) - mean() * mean());
}
double temperatureAccumulatorX_;
double temperatureAccumulatorY_;
double temperatureAccumulatorZ_;
double squareAccumulatorX_;
double squareAccumulatorY_;
double squareAccumulatorZ_;
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);
if (stateData)
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(static_cast<unsigned>(((*existent)->z() * 5.0) / columnHeight));
if (particularStrata >= strata)
particularStrata = strata - 1;
if (particularStrata < 0)
particularStrata = 0;
double tx((*existent)->m() * (*existent)->vx() * (*existent)->vx() / boltzmann_k);
double ty((*existent)->m() * (*existent)->vy() * (*existent)->vy() / boltzmann_k);
double tz((*existent)->m() * (*existent)->vz() * (*existent)-
vz() / boltzmann_k);
temperatures[particularStrata].add(tx, ty, tz);
++existent;
}
for (unsigned stratum(0); stratum < strata; ++stratum)
{
if (temperatureStratumData)
std::cout << "in stratum " << stratum << " there were " <<
temperatures[stratum].count_ << " particles" << std::endl;
if (temperatures[stratum].count_)
{
if (temperatureStratumData)
std::cout << " tx: " << temperatures[stratum].tx()
<< " ty: " << temperatures[stratum].ty()
<< " tz: " << temperatures[stratum].tz()
<< " Tavg: " << temperatures[stratum].mean()
<< " std dev: +/-" <<
temperatures[stratum].standardDeviation() << std::endl;
if (0 == stratum)
bottom.add(temperatures[stratum].tx(),
temperatures[stratum].ty(), temperatures[stratum].tz());
else if ((strata - 1) == stratum)
top.add(temperatures[stratum].tx(),
temperatures[stratum].ty(), temperatures[stratum].tz());
}
else if (temperatureStratumData)
std::cout << "no temperature information possible" <<
std::endl;
}
}
// dump full run statistics
std::cout << "over the simulation we have"
<< "\ntop strata of the column - tx: " << top.tx() << " ty: " <<
top.ty() << " tz: " << top.tz()
<< " Tavg: " << top.mean() << " std dev: +/-" <<
top.standardDeviation()
<< "\nbottom strata of the column - tx: " << bottom.tx() << "
ty: " << bottom.ty() << " tz: " << bottom.tz()
<< " Tavg: " << bottom.mean() << " std dev: +/-" <<
bottom.standardDeviation()
<< 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)
{
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "particles colliding:"
<< "\n 1 - x: " << (*object1)->x() << " y: " << (*object1)->y() << " z: " << (*object1)->z()
<< "\n 1 - vx: " << (*object1)->vx() << " vy: " <<
(*object1)->vy() << " vz: " << (*object1)->vz()
<< "\n 2 - x: " << (*object2)->x() << " y: " << (*object2)->y() << " z: " << (*object2)->z()
<< "\n 2 - vx: " << (*object2)->vx() << " vy: " <<
(*object2)->vy() << " vz: " << (*object2)->vz() << std::endl;
// calculate relative distances and speed
double deltaX((*object2)->x() - (*object1)->x());
double deltaY((*object2)->y() - (*object1)->y());
double deltaZ((*object2)->z() - (*object1)->z());
double deltaVX((*object2)->vx() - (*object1)->vx());
double deltaVY((*object2)->vy() - (*object1)->vy());
double deltaVZ((*object2)->vz() - (*object1)->vz());
double relativeDistance(distance(deltaX, deltaY, deltaZ));
double relativeVelocity(distance(deltaVX, deltaVY, deltaVZ));
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "relativeDistance: " << relativeDistance << "
relativeVelocity: " << relativeVelocity << std::endl;
// find relative angle
double theta(std::acos(deltaZ / relativeDistance));
double phi(((0 == deltaX) && (0 == deltaY)) ? 0. :
std::atan2(deltaY, deltaX));
if ((*object1)->name() == "0")
std::cout << "theta: " << theta << " phi: " << phi <<
std::endl;
// rotate to canonical coordinates
double rotatedVelocityX(- std::cos(theta) * std::cos(phi) *
deltaVX - std::cos(theta) * std::sin(phi) * deltaVY + std::sin(theta)
* deltaVZ);
double rotatedVelocityY(std::sin(phi) * deltaVX - std::cos(phi)
* deltaVY);
double rotatedVelocityZ(- std::sin(theta) * std::cos(phi) *
deltaVX - std::sin(theta) * std::sin(phi) * deltaVY - std::cos(theta)
* deltaVZ);
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "rotatedVelocityX: " << rotatedVelocityX << "
rotatedVelocityY: "
<< rotatedVelocityY << " rotatedVelocityZ: " <<
rotatedVelocityZ << std::endl;
// find relative velocity angles
double velocityTheta(std::acos(rotatedVelocityZ /
relativeVelocity));
double velocityPhi(((0 == rotatedVelocityX) && (0 ==
rotatedVelocityY)) ? 0. : std::atan2(rotatedVelocityY,
rotatedVelocityX));
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "velocityTheta: " << velocityTheta << "
velocityPhi: " << velocityPhi << std::endl;
// classic impact parameter
double impact(relativeDistance * std::sin(velocityTheta) /
(boost::static_pointer_cast<Particle>(*object1)->radius() +
boost::static_pointer_cast<Particle>(*object2)->radius()));
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "impact: " << impact << std::endl;
// reverse motion to collision time
double collisionTime((relativeDistance*std::cos(velocityTheta) -
(boost::static_pointer_cast<Particle>(*object1)->radius() +
boost::static_pointer_cast<Particle>(*object2)->radius()) *
std::sqrt(1. - impact * impact)) / relativeVelocity);
(*object1)->x((- deltaVX + (*object2)->vx()) * collisionTime +
(*object1)->x());
(*object1)->y((- deltaVY + (*object2)->vy()) * collisionTime +
(*object1)->y());
(*object1)->z((- deltaVZ + (*object2)->vz()) * collisionTime +
(*object1)->z());
(*object2)->x(deltaX + deltaVX * collisionTime + (*object1)->x());
(*object2)->y(deltaY + deltaVY * collisionTime + (*object1)->y());
(*object2)->z(deltaZ + deltaVZ * collisionTime + (*object1)-
z());
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "collisionTime: " << collisionTime << std::endl;
// calculate collision angles
double impactAngle(std::asin(-impact));
double velocityAngleTotal(std::tan(velocityTheta +
impactAngle));
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "impactAngle: " << impactAngle << "
velocityAngleTotal: " << velocityAngleTotal << std::endl;
// mass impact
double mu((*object2)->m() / (*object1)->m());
double diffZ(2. * (rotatedVelocityZ + velocityAngleTotal *
(std::cos(velocityPhi) * rotatedVelocityX +
std::sin(velocityPhi) * rotatedVelocityY)) /
((1 + velocityAngleTotal * velocityAngleTotal) * (1
+ mu)));
// boost coordinate frame for impact
double boostedVelocityX(velocityAngleTotal *
std::cos(velocityPhi) * diffZ);
double boostedVelocityY(velocityAngleTotal *
std::sin(velocityPhi) * diffZ);
double boostedVelocityZ(diffZ);
rotatedVelocityX -= mu * boostedVelocityX;
rotatedVelocityY -= mu * boostedVelocityY;
rotatedVelocityZ -= mu * boostedVelocityZ;
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "boostedVelocityX: " << boostedVelocityX << "
boostedVelocityY: "
<< boostedVelocityY << " boostedVelocityZ: " <<
boostedVelocityZ
<< "\n and updated rotatedVelocityX: " << rotatedVelocityX
<< " rotatedVelocityY: "
<< rotatedVelocityY << " rotatedVelocityZ: " <<
rotatedVelocityZ << std::endl;
// and finally transform velocities
(*object1)->vx(std::cos(theta) * std::cos(phi) *
rotatedVelocityX - std::sin(phi) * rotatedVelocityY
+ std::sin(theta) * std::cos(phi) *
rotatedVelocityZ + (*object2)->vx());
(*object1)->vy(std::cos(theta) * std::sin(phi) *
rotatedVelocityX + std::cos(phi) * rotatedVelocityY
+ std::sin(theta) * std::sin(phi) *
rotatedVelocityZ + (*object2)->vy());
(*object1)->vz(- std::sin(theta) * rotatedVelocityX +
std::cos(theta) * rotatedVelocityZ + (*object2)->vz());
(*object2)->vx(std::cos(theta) * std::cos(phi) *
boostedVelocityX - std::sin(phi) * boostedVelocityY
+ std::sin(theta) * std::cos(phi) *
boostedVelocityZ + (*object2)->vx());
(*object2)->vy(std::cos(theta) * std::sin(phi) *
boostedVelocityX + std::cos(phi) * boostedVelocityY
+ std::sin(theta) * std::sin(phi) *
boostedVelocityZ + (*object2)->vy());
(*object2)->vz(- std::sin(theta) * boostedVelocityX +
std::cos(theta) * boostedVelocityZ + (*object2)->vz());
if (particleCollisionCalculations && ((*object1)->name() ==
"0"))
std::cout << "particles leaving:"
<< "\n 1 - x: " << (*object1)->x() << " y: " << (*object1)->y() << " z: " << (*object1)->z()
<< "\n 1 - vx: " << (*object1)->vx() << " vy: " <<
(*object1)->vy() << " vz: " << (*object1)->vz()
<< "\n 2 - x: " << (*object2)->x() << " y: " << (*object2)->y() << " z: " << (*object2)->z()
<< "\n 2 - vx: " << (*object2)->vx() << " vy: " <<
(*object2)->vy() << " vz: " << (*object2)->vz() << std::endl;
}
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)->vy() == 0) ? 0. : std::atan2((*object2)->vy(), (*object2)->vx()));
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, simulationTime,
simulationDelta);
engine.run();
}
catch (std::exception & blowUp)
{
std::cout << "exception at lowest level: " << blowUp.what() <<
std::endl;
}
catch (...)
{
std::cout << "blowUps happen" << std::endl;
}
}
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
galathaea: prankster, fablist, magician, liar