simulations contra the terroristic nimbus of the second law
On Mar 7, 7:46 pm, Publius <m.publ...@nospam.comcast.net> wrote:
"galathaea" <galath...@gmail.com> wrote in news:1173314396.031294.195110
@q40g2000cwq.googlegroups.com:
at this point
i have enough of the "temperature" to show
that a rarefied heat bath at equilibrium
will have a delta capable of driving a heat engine
which violates the second law
Interesting stuff, g. Probably worth pointing out, though, that the
"violation" is only apparent. Gravity produces a temperature stratification
(a delta T) from which energy may (in principle) be extracted. But unless
that energy is replaced, e.g., by solar heating, the delta T will quickly
become too small to extract any further energy. Basically we are "stealing"
from the system heat which would otherwise drive convection.
absolutely
that is the point
there is no violation
of the conservation of energy
extracting heat from a heat bath cools it
the friction and other dissipative forces
near turbulence and free energy
generate heat
so we have a completely reversible system
if you can turn a heat bath
into something that generates work
you have a perpetuum mobile of the second kind
all energy becomes potentially infinitely reusable
this is actually a problem in information theory as well
distinguishing force and work
from heat
is actually a distinction of _use_ or _understanding_
that is why thermodynamics was so important to the industrial
revolution
because it allowed us to automate uses of energy
understanding engines
meant understanding how to make the universe
do what we _intend_
gravity stratification
( or general field stratification )
means we can recover energy for use
in a physics that associates information with energy
in the manner of maxwell's demon solutions
field stratification means we can extract
a potentially infinite amount of energy from the universe
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