step/stepcore/coulombforce.cc

0001 /*
0002     SPDX-FileCopyrightText: 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
0003
0004     SPDX-License-Identifier: GPL-2.0-or-later
0005 */
0006
0007 #include "coulombforce.h"
0008 #include "particle.h"
0009
0010 #include <cmath>
0011
0012 namespace StepCore {
0013
0014 STEPCORE_META_OBJECT(CoulombForce, QT_TRANSLATE_NOOP("ObjectClass", "CoulombForce"), QT_TRANSLATE_NOOP("ObjectDescription", "Coulomb force"), 0,
0015     STEPCORE_SUPER_CLASS(Item) STEPCORE_SUPER_CLASS(Force),
0016     STEPCORE_PROPERTY_RW(double, coulombConst, QT_TRANSLATE_NOOP("PropertyName", "coulombConst"), STEPCORE_FROM_UTF8(QT_TRANSLATE_NOOP("Units", "N m²/C²")),
0017                 QT_TRANSLATE_NOOP("PropertyDescription", "Coulomb constant"), coulombConst, setCoulombConst))
0018
0019 STEPCORE_META_OBJECT(CoulombForceErrors, QT_TRANSLATE_NOOP("ObjectClass", "CoulombForceErrors"), QT_TRANSLATE_NOOP("ObjectDescription", "Errors class for CoulombForce"), 0,
0020     STEPCORE_SUPER_CLASS(ObjectErrors),
0021     STEPCORE_PROPERTY_RW(double, coulombConstVariance, QT_TRANSLATE_NOOP("PropertyName", "coulombConstVariance"), STEPCORE_FROM_UTF8(QT_TRANSLATE_NOOP("Units", "N m²/C²")),
0022                 QT_TRANSLATE_NOOP("PropertyDescription", "Coulomb constant variance"), coulombConstVariance, setCoulombConstVariance))
0023
0024 CoulombForce* CoulombForceErrors::coulombForce() const
0025 {
0026     return static_cast<CoulombForce*>(owner());
0027 }
0028
0029 CoulombForce::CoulombForce(double coulombConst)
0030     : Force()
0031     , _coulombConst(coulombConst)
0032 {
0033     coulombForceErrors()->setCoulombConstVariance(square(Constants::CoulombError));
0034 }
0035
0036 void CoulombForce::calcForce(bool calcVariances)
0037 {
0038     const BodyList::const_iterator end = world()->bodies().end();
0039     for(BodyList::const_iterator b1 = world()->bodies().begin(); b1 != end; ++b1) {
0040         if(!(*b1)->metaObject()->inherits<ChargedParticle>()) continue;
0041         for(BodyList::const_iterator b2 = b1+1; b2 != end; ++b2) {
0042             if(!(*b2)->metaObject()->inherits<ChargedParticle>()) continue;
0043             ChargedParticle* p1 = static_cast<ChargedParticle*>(*b1);
0044             ChargedParticle* p2 = static_cast<ChargedParticle*>(*b2);
0045
0046             Vector2d r = p2->position() - p1->position();
0047             double rnorm2 = r.squaredNorm();
0048             Vector2d force = _coulombConst* p1->charge() * p2->charge() * r / (rnorm2*sqrt(rnorm2));
0049             p2->applyForce(force);
0050             force = -force;
0051             p1->applyForce(force);
0052
0053             if(calcVariances) {
0054                 // XXX: CHECKME
0055                 ChargedParticleErrors* pe1 = p1->chargedParticleErrors();
0056                 ChargedParticleErrors* pe2 = p2->chargedParticleErrors();
0057                 Vector2d rV = pe2->positionVariance() + pe1->positionVariance();
0058                 double chargeVariance = coulombForceErrors()->_coulombConstVariance / square(_coulombConst) +
0059                                         pe1->chargeVariance() / square(p1->charge()) +
0060                                         pe2->chargeVariance() / square(p2->charge());
0061                 Vector2d forceV = force.array().square()* (
0062                         Vector2d(chargeVariance, chargeVariance) +
0063                         Vector2d(rV[0] * square(1/r[0] - 3*r[0]/rnorm2) + rV[1] * square(3*r[1]/rnorm2),
0064                                  rV[1] * square(1/r[1] - 3*r[1]/rnorm2) + rV[0] * square(3*r[0]/rnorm2))).array();
0065                 pe1->applyForceVariance(forceV);
0066                 pe2->applyForceVariance(forceV);
0067             }
0068         }
0069     }
0070 }
0071
0072 } // namespace StepCore
0073