Select Git revision
Particle.cpp 8.21 KiB
///
/// @file Particle.cpp
/// @author Armin Co
///
/// @brief Particle class implemantation and main().
///
#include <iostream>
#include <thread>
#include <vector>
#include "Common.hpp"
#include "Particle.hpp"
/// @brief Planck constant.
///
constexpr double hbar = 0.0002;
const std::complex<double> Ci(0.0, 1.0);
/// @brief Takes the row index of an array and turns it
/// into a position.
///
const double xAt(int x)
{
return 2.0 * x / QPong::arrayWidth - 1.0;
}
/// @brief Takes the line index of an array and turns it
/// into a position.
///
const double yAt(int y)
{
return 2.0 * y / QPong::arrayHeight - 1.0;
}
void Particle::initMomentum()
{
double sumAtInit = 0.0;
for (unsigned iy = 0; iy < QPong::arrayHeight; iy++)
{
double y = yAt(iy);
for (unsigned ix = 0; ix < QPong::arrayWidth; ix++)
{
double x = xAt(ix);
unsigned index = QPong::arrayWidth * iy + ix;
constexpr double k_px0 = 0.063;
constexpr double k_py0 = 0.042;
constexpr double k_x0 = 0.0;
constexpr double k_y0 = -0.0;
constexpr double lambda = 60.0;
constexpr double A0 = 0.5;
m_psi[index] = A0 * exp(Ci / hbar * (x * k_px0 + y * k_py0) - lambda * (sqr(x - k_x0) + sqr(y - k_y0)));
sumAtInit += sqr(real(m_psi[index])) + sqr(imag(m_psi[index]));
}
}
for (unsigned iy = 0; iy < QPong::arrayHeight; iy++)
{
double y = yAt(iy);
for (unsigned ix = 0; ix < QPong::arrayWidth; ix++)
{
double x = xAt(ix);
unsigned index = QPong::arrayWidth * iy + ix;
m_psi[index] *= 1.0 / sqrt(sumAtInit);
}
}
}
Particle::Particle(ImageBuffer *momentum, ImageBuffer *position, std::map<Player::Id, std::shared_ptr<Player>> players)
: m_bufMomentumRepresentation{momentum}
, m_bufPositionRepresentation{position}
, m_ready{false}
, m_players {players}
{
fftw_init_threads();
fftw_plan_with_nthreads(QPong::numberOfThreads);
m_psi = static_cast<std::complex<double> *>(fftw_malloc(sizeof(fftw_complex) * QPong::arraySize));
m_planForward = fftw_plan_dft_2d(QPong::arrayHeight, QPong::arrayWidth, (fftw_complex *)m_psi, (fftw_complex *)m_psi, FFTW_FORWARD, FFTW_MEASURE);
m_planBackward = fftw_plan_dft_2d(QPong::arrayHeight, QPong::arrayWidth, (fftw_complex *)m_psi, (fftw_complex *)m_psi, FFTW_BACKWARD, FFTW_MEASURE);
initMomentum();
m_ready = true;
}
Particle::~Particle()
{
fftw_free((fftw_complex *)m_psi);
fftw_destroy_plan(m_planForward);
fftw_destroy_plan(m_planBackward);
}
bool Particle::notReady()
{
return !m_ready;
}
double Particle::removeParticle(int index, double cx)
{
double normPre = sqr(std::real(m_psi[index])) + sqr(std::imag(m_psi[index]));
double e = sqr(cos(cx));
m_psi[index] *= e;
double dif = normPre - (sqr(std::real(m_psi[index])) + sqr(std::imag(m_psi[index])));
return dif;
}
const double batPotential(double x, double y, Player::Position pos)
{
y -= pos.y;
double E = 0.0;
constexpr double scale = 0.000015;
constexpr double height = 0.31;
constexpr double cutoff = 0.0015;
if (y >= height)
{
double d = sqr(x-pos.x) + sqr(y-height) + cutoff;
E = 1.0 / d * scale;
}
else if (y <= (-height))
{
double d = sqr(x-pos.x) + sqr((-y)-height) + cutoff;
E = 1.0 / d * scale;
}
else
{
double d = sqr(pos.x-x) + cutoff;
E = 1.0 / d * scale;
}
return E;
}
constexpr double dt = 0.15;
void Particle::manipulateParticleInPosition(int indexFrom, int indexTo){
auto pOneX = m_players[Player::Id::One]->getPosition().x;
auto pTwoX = m_players[Player::Id::Two]->getPosition().x;
float sumScorePlayerOne = 0.0;
float sumScorePlayerTwo = 0.0;
for (int index = indexFrom; index < indexTo; index++)
{
double x = xAt(index % QPong::arrayWidth);
double y = yAt(index / QPong::arrayWidth);
if (x < pOneX) //< expects p1.x < 0
{
double cx = -pOneX + x;
cx = cx * (M_PI / 2.0) / (1.0 + pOneX);
double dif = removeParticle(index, cx);
sumScorePlayerTwo += dif;
}
if (x > pTwoX) //< expects p2.x > 0
{
double cx = x - pTwoX;
cx = cx * (M_PI / 2.0) / (1.0 - pTwoX);
double dif = removeParticle(index, cx);
sumScorePlayerOne += dif;
}
double playerOneBat = batPotential(x, y, m_players[Player::Id::One]->getPosition());
double playerTwoBat = batPotential(x, y, m_players[Player::Id::Two]->getPosition());
double horizontalBorders = sqr(y) * sqr(y) * sqr(y) * sqr(y) * 0.01;
auto potentialSum = playerOneBat + playerTwoBat + horizontalBorders;
auto tmpPsi = m_psi[index] * exp(-Ci * dt / hbar * potentialSum);
m_psi[index] = tmpPsi;
m_bufPositionRepresentation->updateAt(index, tmpPsi * 125.0, potentialSum * 150);
}
m_players[Player::Id::One]->addScore(sumScorePlayerOne);
m_players[Player::Id::Two]->addScore(sumScorePlayerTwo);
}
double flattenMomentum(double x, double y)
{
double distanceToCenter = sqr(x) + sqr(y);
constexpr double offset = 0.012;
if (distanceToCenter > offset)
{
return 1.0 / (1.0 + (distanceToCenter * 2.5));
}
else
{
return 1.0;
}
}
const double pxAt(int ix)
{
if (ix > QPong::arrayWidth / 2)
{
ix -= QPong::arrayWidth;
}
return ix * M_PI * hbar;
}
const double pyAt(int iy)
{
if (iy > QPong::arrayHeight / 2)
{
iy -= QPong::arrayHeight;
}
return iy * M_PI * hbar;
}
void Particle::moveStep(int indexFrom, int indexTo)
{
for (int i = indexFrom; i < indexTo; i++)
{
double py = pyAt(i / QPong::arrayWidth);
double px = pxAt(i % QPong::arrayWidth);
constexpr double m = 1.0;
constexpr double N = 1.0 / QPong::arraySize;
double E = (sqr(px) + sqr(py)) / (2.0 * m);
auto f = flattenMomentum(px, py);
m_psi[i] *= exp(-Ci * dt / hbar * E) * N * f;
// m_psi[i] = flattenMomentum(px, py) * 0.00001;
}
}
void Particle::propagate()
{
for (auto &player : m_players)
{
player.second->update(dt);
}
std::vector<std::thread> positionThreads;
for (int i = 0; i < QPong::numberOfThreads; i++)
{
int from = i * (QPong::arraySize/QPong::numberOfThreads);
int to = (i+1) * (QPong::arraySize/QPong::numberOfThreads);
positionThreads.emplace_back(&Particle::manipulateParticleInPosition, this, i * (QPong::arraySize/QPong::numberOfThreads), (i+1) * (QPong::arraySize/QPong::numberOfThreads));
}
for (auto &thread : positionThreads)
{
thread.join();
}
fftw_execute (m_planForward); // transform into momentum repr.
std::vector<std::thread> momentumThreads;
for (int i = 0; i < QPong::numberOfThreads; i++)
{
momentumThreads.emplace_back(&Particle::moveStep, this, i * (QPong::arraySize / QPong::numberOfThreads), (i + 1) * (QPong::arraySize / QPong::numberOfThreads));
}
for (auto &thread : momentumThreads)
{
thread.join();
}
/// @todo Pull update Momentum image into moveStep()
updateMomentumImage(); //momentum
fftw_execute(m_planBackward); // transform into position repr.
}
void Particle::updateMomentumImage()
{
for (int index = 0; index < QPong::arraySize; index++)
{
auto y = index / QPong::arrayWidth;
unsigned py = 0;
if (y >= QPong::arrayHeight / 2)
{
py = y - QPong::arrayHeight / 2;
}
else
{
py = y + QPong::arrayHeight / 2;
}
auto x = index % QPong::arrayWidth;
unsigned px = 0; if (x >= QPong::arrayWidth / 2)
{
px = x - QPong::arrayWidth / 2;
}
else
{
px = x + QPong::arrayWidth / 2;
}
constexpr double scale = 55580.0;;
auto array_index = y * QPong::arrayWidth + x;
auto momentum_index = py * QPong::arrayWidth + px;
m_bufMomentumRepresentation->updateAt(array_index, (m_psi[momentum_index] * scale), 0.0);
}
}