step/stepcore

particle.cc

00001 /* This file is part of StepCore library.
00002    Copyright (C) 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
00003 
00004    StepCore library is free software; you can redistribute it and/or modify
00005    it under the terms of the GNU General Public License as published by
00006    the Free Software Foundation; either version 2 of the License, or
00007    (at your option) any later version.
00008 
00009    StepCore library is distributed in the hope that it will be useful,
00010    but WITHOUT ANY WARRANTY; without even the implied warranty of
00011    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012    GNU General Public License for more details.
00013 
00014    You should have received a copy of the GNU General Public License
00015    along with StepCore; if not, write to the Free Software
00016    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017 */
00018 
00019 #include "particle.h"
00020 #include "types.h"
00021 #include <cstring>
00022 #include <cmath>
00023 
00024 namespace StepCore
00025 {
00026 
00027 STEPCORE_META_OBJECT(Particle, "Simple zero-size particle", 0,
00028         STEPCORE_SUPER_CLASS(Item) STEPCORE_SUPER_CLASS(Body),
00029         STEPCORE_PROPERTY_RW_D(StepCore::Vector2d, position, "m", "position", position, setPosition)
00030         STEPCORE_PROPERTY_RW_D(StepCore::Vector2d, velocity, "m/s", "velocity", velocity, setVelocity)
00031         STEPCORE_PROPERTY_R_D(StepCore::Vector2d, acceleration, STEPCORE_FROM_UTF8("m/s²"),
00032                                                             "acceleration", acceleration)
00033         STEPCORE_PROPERTY_R_D(StepCore::Vector2d, force, "N", "force", force)
00034         STEPCORE_PROPERTY_RW(double, mass, "kg", "mass", mass, setMass)
00035         STEPCORE_PROPERTY_RWF(StepCore::Vector2d, momentum, "kg m/s", "momentum",
00036                         StepCore::MetaProperty::DYNAMIC, momentum, setMomentum)
00037         STEPCORE_PROPERTY_RWF(double, kineticEnergy, "J", "kinetic energy",
00038                         StepCore::MetaProperty::DYNAMIC, kineticEnergy, setKineticEnergy))
00039 
00040 STEPCORE_META_OBJECT(ParticleErrors, "Errors class for Particle", 0, STEPCORE_SUPER_CLASS(ObjectErrors),
00041         STEPCORE_PROPERTY_RW_D(StepCore::Vector2d, positionVariance, "m",
00042                     "position variance", positionVariance, setPositionVariance)
00043         STEPCORE_PROPERTY_RW_D(StepCore::Vector2d, velocityVariance, "m/s",
00044                     "velocity variance", velocityVariance, setVelocityVariance)
00045         STEPCORE_PROPERTY_R_D(StepCore::Vector2d, accelerationVariance, STEPCORE_FROM_UTF8("m/s²"),
00046                     "acceleration variance", accelerationVariance)
00047         STEPCORE_PROPERTY_R_D(StepCore::Vector2d, forceVariance, "N",
00048                     "force variance", forceVariance)
00049         STEPCORE_PROPERTY_RW(double, massVariance, "kg",
00050                     "mass variance", massVariance, setMassVariance )
00051         STEPCORE_PROPERTY_RWF(StepCore::Vector2d, momentumVariance, "kg m/s",
00052                     "momentum variance", StepCore::MetaProperty::DYNAMIC, momentumVariance, setMomentumVariance)
00053         STEPCORE_PROPERTY_RWF(double, kineticEnergyVariance, "J",
00054                     "kinetic energy variance", StepCore::MetaProperty::DYNAMIC, kineticEnergyVariance, setKineticEnergyVariance))
00055 
00056 STEPCORE_META_OBJECT(ChargedParticle, "Charged zero-size particle", 0, STEPCORE_SUPER_CLASS(Particle),
00057         STEPCORE_PROPERTY_RW(double, charge, "C", "charge", charge, setCharge))
00058 
00059 STEPCORE_META_OBJECT(ChargedParticleErrors, "Errors class for ChargedParticle", 0,
00060         STEPCORE_SUPER_CLASS(ParticleErrors),
00061         STEPCORE_PROPERTY_RW(double, chargeVariance, "kg",
00062                     "charge variance", chargeVariance, setChargeVariance ))
00063 
00064 Particle* ParticleErrors::particle() const
00065 {
00066     return static_cast<Particle*>(owner());
00067 }
00068 
00069 Vector2d ParticleErrors::accelerationVariance() const
00070 {
00071     return _forceVariance/square(particle()->mass()) +
00072         _massVariance*(particle()->force()/square(particle()->mass())).cSquare();
00073 }
00074 
00075 Vector2d ParticleErrors::momentumVariance() const
00076 {
00077     return _velocityVariance * square(particle()->mass()) +
00078            particle()->velocity().cSquare() * _massVariance;
00079 }
00080 
00081 void ParticleErrors::setMomentumVariance(const Vector2d& momentumVariance)
00082 {
00083     _velocityVariance = (momentumVariance - particle()->velocity().cSquare() * _massVariance) /
00084                         square(particle()->mass());
00085 }
00086 
00087 double ParticleErrors::kineticEnergyVariance() const
00088 {
00089     return _velocityVariance.innerProduct(particle()->velocity().cSquare()) * square(particle()->mass()) +
00090            square(particle()->velocity().norm2()/2) * _massVariance;
00091 }
00092 
00093 void ParticleErrors::setKineticEnergyVariance(double kineticEnergyVariance)
00094 {
00095     _velocityVariance = (kineticEnergyVariance - square(particle()->velocity().norm2()/2) * _massVariance) /
00096                         square(particle()->mass()) / 2 *
00097                         Vector2d(1,1).cDivide(particle()->velocity().cSquare());
00098     if(!std::isfinite(_velocityVariance[0]) || _velocityVariance[0] < 0 ||
00099        !std::isfinite(_velocityVariance[1]) || _velocityVariance[1]) {
00100         _velocityVariance.setZero();
00101     }
00102 }
00103 
00104 ChargedParticle* ChargedParticleErrors::chargedParticle() const
00105 {
00106     return static_cast<ChargedParticle*>(owner());
00107 }
00108 
00109 Particle::Particle(Vector2d position, Vector2d velocity, double mass)
00110     : _position(position), _velocity(velocity), _force(0), _mass(mass)
00111 {
00112 }
00113 
00114 void Particle::getVariables(double* position, double* velocity,
00115                      double* positionVariance, double* velocityVariance)
00116 {
00117     std::memcpy(position, _position.array(), 2*sizeof(*position));
00118     std::memcpy(velocity, _velocity.array(), 2*sizeof(*velocity));
00119     if(positionVariance) {
00120         ParticleErrors* pe = particleErrors();
00121         std::memcpy(positionVariance, pe->_positionVariance.array(), 2*sizeof(*positionVariance));
00122         std::memcpy(velocityVariance, pe->_velocityVariance.array(), 2*sizeof(*velocityVariance));
00123     }
00124 }
00125 
00126 void Particle::setVariables(const double* position, const double* velocity,
00127                const double* positionVariance, const double* velocityVariance)
00128 {
00129     std::memcpy(_position.array(), position, 2*sizeof(*position));
00130     std::memcpy(_velocity.array(), velocity, 2*sizeof(*velocity));
00131     _force.setZero();
00132     if(positionVariance) {
00133         ParticleErrors* pe = particleErrors();
00134         std::memcpy(pe->_positionVariance.array(), positionVariance, 2*sizeof(*positionVariance));
00135         std::memcpy(pe->_velocityVariance.array(), velocityVariance, 2*sizeof(*positionVariance));
00136         pe->_forceVariance.setZero();
00137     }
00138 }
00139 
00140 void Particle::getAccelerations(double* acceleration, double* accelerationVariance)
00141 {
00142     acceleration[0] = _force[0] / _mass;
00143     acceleration[1] = _force[1] / _mass;
00144     if(accelerationVariance) {
00145         ParticleErrors* pe = particleErrors();
00146         accelerationVariance[0] = pe->_forceVariance[0]/square(_mass) +
00147                                         square(_force[0]/square(_mass))*pe->_massVariance;
00148         accelerationVariance[1] = pe->_forceVariance[1]/square(_mass) +
00149                                         square(_force[1]/square(_mass))*pe->_massVariance;
00150     }
00151 }
00152 
00153 void Particle::addForce(const double* force, const double* forceVariance)
00154 {
00155     _force[0] += force[0];
00156     _force[1] += force[1];
00157     if(forceVariance) {
00158         ParticleErrors* pe = particleErrors();
00159         pe->_forceVariance[0] += forceVariance[0];
00160         pe->_forceVariance[1] += forceVariance[1];
00161     }
00162 }
00163 
00164 void Particle::resetForce(bool resetVariance)
00165 {
00166     _force.setZero();
00167     if(resetVariance) particleErrors()->_forceVariance.setZero();
00168 }
00169 
00170 void Particle::getInverseMass(GmmSparseRowMatrix* inverseMass,
00171                         GmmSparseRowMatrix* variance, int offset)
00172 {
00173     inverseMass->row(offset).w(offset, 1/_mass);
00174     inverseMass->row(offset+1).w(offset+1, 1/_mass);
00175     if(variance) {
00176         double v = particleErrors()->_massVariance / square(square(_mass));
00177         variance->row(offset).w(offset, v);
00178         variance->row(offset+1).w(offset+1, v);
00179     }
00180 }
00181 
00182 void Particle::setKineticEnergy(double kineticEnergy)
00183 {
00184     double v = _velocity.norm();
00185     _velocity = sqrt(kineticEnergy*2/_mass) * (v>0 ? _velocity/v : Vector2d(1,0));
00186 }
00187 
00188 } // namespace StepCore
00189
step/stepcore

particle.cc

step/stepcore

kdeedu
