FYS3150/Project-5
3
0
Fork 0
Code Issues Pull Requests Packages Projects 1 Releases Wiki Activity

Make some changes

Browse Source
This commit is contained in:
Cory Balaton
2023-12-19 16:08:58 +01:00
parent 3e991f9918
commit e22ceb04d7
3 changed files with 107 additions and 65 deletions
Download Patch File Download Diff File Expand all files Collapse all files

83
lib/WaveSimulation.cpp
View File

@@ -18,17 +18,36 @@
#include <cstdint>
#include <cstdlib>
WaveSimulation::WaveSimulation(double h, double dt, double T)
WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c,
double y_c, double sigma_x, double sigma_y,
double p_x, double p_y, double thickness,
double pos_x, double ap_sep, double ap,
uint32_t slits)
{
this->dt = dt;
this->h = h;
this->T = T;
this->M = 1. / h;
this->N = M - 2;
this->V.set_size(this->N, this->N);
this->V.fill(0.);
this->U.set_size(this->N, this->N);
this->U.fill(0.);
this->build_V(thickness, pos_x, ap_sep, ap, slits);
this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
this->build_A();
this->build_B();
}
WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c, double y_c,
double sigma_x, double sigma_y, double p_x, double p_y)
{
this->dt = dt;
this->h = h;
this->T = T;
this->M = 1. / h;
this->N = M - 2;
this->build_V(0.,0.,0.,0., 0);
this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
this->build_A();
this->build_B();
}
void WaveSimulation::solve(std::ofstream &ofile)
@@ -43,7 +62,6 @@ void WaveSimulation::solve(std::ofstream &ofile)
void WaveSimulation::step()
{
DEBUG("Inside step");
arma::cx_vec tmp = this->B * this->U.as_col();
arma::spsolve(this->U, this->A, tmp);
}
@@ -111,6 +129,7 @@ void WaveSimulation::build_B()
void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x,
double sigma_y, double p_x, double p_y)
{
this->U.set_size(this->N, this->N);
double x, y, diff_x, diff_y;
std::complex<double> sum = 0.;
for (size_t j = 0; j < this->U.n_cols; j++) {
@@ -145,25 +164,42 @@ void WaveSimulation::build_V(double thickness, double pos_x,
double aperture_separation, double aperture,
uint32_t slits)
{
uint32_t mid_y = this->N / 2 - (this->N % 2 == 0);
this->V.set_size(this->N, this->N);
this->V.fill(0.);
if (slits == 0) {
return;
}
arma::cx_vec res;
if (slits % 2 == 0) {
res = arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10));
for (size_t i=0; i < slits; i+=2) {
res = arma::join_cols(res, arma::cx_vec(aperture/this->h,arma::fill::zeros));
res = arma::join_cols(arma::cx_vec(aperture/this->h,arma::fill::zeros), res);
res = arma::join_cols(res, arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)));
res = arma::join_cols(arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)), res);
res = arma::cx_vec(aperture_separation / this->h,
arma::fill::value(1e10));
for (size_t i = 0; i < slits; i += 2) {
res = arma::join_cols(
res, arma::cx_vec(aperture / this->h, arma::fill::zeros));
res = arma::join_cols(
arma::cx_vec(aperture / this->h, arma::fill::zeros), res);
res =
arma::join_cols(res, arma::cx_vec(aperture_separation / this->h,
arma::fill::value(1e10)));
res = arma::join_cols(arma::cx_vec(aperture_separation / this->h,
arma::fill::value(1e10)),
res);
}
}
else {
res = arma::cx_vec(aperture/this->h,arma::fill::value(0));
for (size_t i=0; i < slits-1; i+=2) {
res = arma::join_cols(res, arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)));
res = arma::join_cols(arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)), res);
res = arma::join_cols(res, arma::cx_vec(aperture/this->h,arma::fill::zeros));
res = arma::join_cols(arma::cx_vec(aperture/this->h,arma::fill::zeros), res);
res = arma::cx_vec(aperture / this->h, arma::fill::value(0));
for (size_t i = 0; i < slits - 1; i += 2) {
res =
arma::join_cols(res, arma::cx_vec(aperture_separation / this->h,
arma::fill::value(1e10)));
res = arma::join_cols(arma::cx_vec(aperture_separation / this->h,
arma::fill::value(1e10)),
res);
res = arma::join_cols(
res, arma::cx_vec(aperture / this->h, arma::fill::zeros));
res = arma::join_cols(
arma::cx_vec(aperture / this->h, arma::fill::zeros), res);
}
}
if (res.size() > this->N) {
@@ -171,10 +207,11 @@ void WaveSimulation::build_V(double thickness, double pos_x,
}
uint32_t fill = (this->N - res.size()) / 2;
res = arma::join_cols(arma::cx_vec(fill, arma::fill::value(1e10)), res);
res = arma::join_cols(res, arma::cx_vec(fill + ((this->N - res.size()) % 2), arma::fill::value(1e10)));
res = arma::join_cols(res, arma::cx_vec(fill + ((this->N - res.size()) % 2),
arma::fill::value(1e10)));
uint32_t start = pos_x/this->h - thickness/this->h/2;
for (size_t i=0; i < thickness/this->h; i++) {
this->V.col(start+i) = res;
uint32_t start = pos_x / this->h - thickness / this->h / 2;
for (size_t i = 0; i < thickness / this->h; i++) {
this->V.col(start + i) = res;
}
}