Page MenuHomec4science
Files F76607821

wavemap.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Aug 9, 06:00

wavemap.cpp
View Options

#include "wavemap.h"
#include "debug.hpp"
wavemap::wavemap(const SVector<double, 2>& _bounds, size_t _N, double dt, const std::function<double(double)>& V)
: bounds(_bounds), N(_N),
dx((bounds[1] - bounds[0]) / N), map(N+1),
H(N+1), A(N+1), B(N+1)
{
// init H
for (size_t k = 0; k < N; ++k)
{
H[k][k] = hbar * hbar / (m * dx * dx) + V(supp(k));
H[k+1][k] = H[k][k+1] = - hbar * hbar / (2 * m * dx * dx);
}
H[N][N] = (hbar * hbar) / (m * dx * dx) + V(bounds[1]);
//npdebug("H = ")
//npdebug(H)
using namespace std::complex_literals;
// init A and B
A = H * (1i * dt / (2 * hbar)) + cmpx(1.0, 0);
B = H * (1i * dt / (-2 * hbar)) + cmpx(1.0, 0);
A[0][0] = A[N][N] = 1;
A[0][1] = A[N-1][N] = A[1][0] = A[N][N-1] = 0;
B[0][0] = B[N][N] = B[0][1] = B[N-1][N] = B[1][0] = B[N][N-1] = 0;
//npdebug()
//npdebug("A = ")
//npdebug(A)
//npdebug()
//npdebug("B = ")
//npdebug(B)
}
void wavemap::evolve()
{
map = A.solve(B * map);
/*npdebug("psi[0] = ", map[0])
npdebug("psi[1] = ", map[1])
npdebug("psi[N-1] = ", map[N-1])
npdebug("psi[N] = ", map[N])*/
}
double wavemap::probability(double xl, double xr) const
{
return spec_value([&](size_t k, const map_t& psi)
{
return psi[k]; // identity operator
}
, (size_t) std::round((xl - bounds[0]) / dx)
, (size_t) std::round((xr - bounds[0]) / dx));
}
double wavemap::energy() const
{
return spec_value(matrix_operator(H));
}
double wavemap::position(size_t m) const
{
return spec_value([&](size_t k, const map_t& psi)
{
return pow(supp(k), m) * psi[k];
});
}
// !!! k > 0, k < N
double wavemap::momentum(size_t m) const
{
using namespace std::complex_literals;
return spec_value([&](size_t k, const map_t& psi)
{
cmpx deriv(0);
switch (m)
{
case 1:
deriv = (psi[k+1] - (k != 0 ? psi[k-1] : 0)) / (2 * dx);
break;
case 2:
deriv = (psi[k+1] - cmpx(2) * psi[k] + (k != 0 ? psi[k-1] : 0)) / (dx*dx);
break;
default:
break;
}
return pow(-1i * hbar, m) * deriv;
});
}
double wavemap::incert_pos() const
{
return sqrt(position(2) - pow(position(1), 2));
}
double wavemap::incert_mom() const
{
return sqrt(momentum(2) - pow(momentum(1), 2));
}

Event Timeline

c4science · Help