12 std::vector<Matrix<T>> velocities;
13 std::vector<Matrix<T>> positions;
57 unsigned int int_order;
58 unsigned long num_particles;
62 T margin, ac_length, mu_gas, T_gas;
64 void (*local_kernel)(T*, T*, T*, T*, T*);
106 this->surface = surface;
114 roots_2d_r = std::get<1>(roots_tuple);
115 roots_2d_c = std::get<0>(roots_tuple);
116 weights_2d = weights.
outer(weights);
120 gamma_2d = roots_2d_r;
121 gamma_dot_2d = roots_2d_c;
132 p_gamma = exp(gamma * gamma * (-1.0) / 2.0) * 1.0 / std::sqrt(2.0 * M_PI);
133 p_gamma_2d = exp(gamma_2d * gamma_2d * (-1.0) / 2.0) * 1.0 / std::sqrt(2.0 * M_PI);
134 p_gamma_dot = exp(gamma_dot * gamma_dot * (-1.0) * ac_length * ac_length / 4.0) * ac_length / std::sqrt(4.0 * M_PI);
135 p_gamma_dot_2d = exp(gamma_dot_2d * gamma_dot_2d * (-1.0) * ac_length * ac_length / 4.0) * ac_length / std::sqrt(4.0 * M_PI);
140 this->local_kernel = local_kernel;
142 this->incident_velocity = incident_velocity;
143 this->num_particles = num_particles;
144 this->sim_name = sim_name;
161 this->local_kernel =
nullptr;
204 std::tuple<T, T>
random_sampler_2D(T (
PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[]);
219 std::tuple<T, T>
random_sampler_angles(T (
PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[]);
229 T
mixture_pdf(T gamma, T gamma_dot, std::vector<T> parameters);
239 T
kirchhoff_pdf(T theta_r1, T theta_r2, std::vector<T> parameters);
249 T
shadow(T theta_r1, T height,
int direction);
260 T
occlusion(T theta_r1, T height, T distance,
int direction);
355 std::vector<trajectory<T>>
import_data(std::string filename);
377 unsigned long steps = 100;
379 T x = x0, x_p = 0.0, acceptance = 0.0;
381 std::random_device rd;
382 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
383 std::mt19937 rng(fullSeed);
384 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
385 std::normal_distribution<T> normDist(0.0f, 1.0f);
387 for(
auto i = 0; i < steps; i ++) {
389 T step = normDist(rng) * size;
390 T u_sample = uniformDist(rng);
393 if(x_p >= domain_x[0] && x_p <= domain_x[1]) {
394 acceptance = std::min(1.0, (this->*pdf)(x_p, parameters) / (this->*pdf)(x, parameters));
396 if(acceptance > u_sample) {
405std::tuple<T, T>
PG_Kernel<T>::random_sampler_2D(T (
PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[]) {
407 unsigned long steps = 100;
413 std::random_device rd;
414 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
415 std::mt19937 rng(fullSeed);
416 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
417 std::normal_distribution<T> normDist(0.0f, 1.0f);
419 for(
auto i = 0; i < steps; i ++) {
421 T step_x = normDist(rng) * size_x;
422 T step_y = normDist(rng) * size_y;
423 T u_sample = uniformDist(rng);
427 if(x_p >= domain_x[0] && x_p <= domain_x[1] && y_p >= domain_y[0] && y_p <= domain_y[1]) {
428 acceptance = std::min(1.0, (this->*pdf)(x_p, y_p, parameters) / (this->*pdf)(x, y, parameters));
430 if(acceptance > u_sample) {
436 return std::make_tuple(x, y);
440std::tuple<T, T>
PG_Kernel<T>::random_sampler_angles(T (
PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T theta_r1, T theta_r2, T size_x, T size_y, T domain_x[], T domain_y[]) {
442 unsigned long steps = 100;
444 std::random_device rd;
445 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
446 std::mt19937 rng(fullSeed);
447 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
448 std::normal_distribution<T> normDist(0.0f, 1.0f);
452 for(
auto i = 0; i < steps; i ++) {
454 T step_x = normDist(rng) * size_x;
455 T step_y = normDist(rng) * size_y;
456 T u_sample = uniformDist(rng);
457 T theta_r1_p = theta_r1 + step_x;
458 T theta_r2_p = theta_r2 + step_y;
460 if(theta_r1_p < 0.0) {
461 theta_r1_p = std::abs(theta_r1_p);
464 if(theta_r2_p < 0.0) {
465 theta_r2_p += 2.0 * M_PI;
467 if(theta_r2_p > 2.0 * M_PI) {
468 theta_r2_p -= 2.0 * M_PI;
471 if(theta_r1_p < domain_x[1]) {
472 acceptance = std::min(1.0, (this->*pdf)(theta_r1_p, theta_r2_p, parameters) / (this->*pdf)(theta_r1, theta_r2, parameters));
475 if(acceptance > u_sample) {
476 theta_r1 = theta_r1_p;
477 theta_r2 = theta_r2_p;
481 return std::make_tuple(theta_r1, theta_r2);
487 bool first = (bool)parameters[0];
488 T theta_i = parameters[1];
489 T distance = parameters[2];
496 T w = std::sqrt(2.0 * (sigma_local * sigma_local / ac_length / ac_length) + sigma_gamma_local * sigma_gamma_local * gamma_dot_local * gamma_dot_local);
497 T eta = 1.0 / std::tan(theta_i);
498 T slope_term = 0.5 * (1.0 + std::erf((eta + mu_dot_local) / (w * std::sqrt(2.0)))) * (std::cos(theta_i) + mu_dot_local * std::sin(theta_i)) +
499 w * std::sin(theta_i) / std::sqrt(2.0 * M_PI) * std::exp(- (eta + mu_dot_local) * (eta + mu_dot_local) / 2.0 / w / w);
500 T p_gamma_local = 1.0 / std::sqrt(2.0 * M_PI) * std::exp( - gamma_local * gamma_local / 2.0);
501 T p_gamma_dot_local = ac_length / std::sqrt(4.0 * M_PI) * std::exp(- gamma_dot_local * gamma_dot_local * ac_length * ac_length / 4.0);
504 return p_gamma_local * p_gamma_dot_local * slope_term * shadow(theta_i, mu_local, -1);
507 return p_gamma_local * p_gamma_dot_local * slope_term * occlusion(theta_i, mu_local, distance, -1);
515 T theta_i = parameters[0];
516 T gamma = parameters[1];
517 T gamma_dot_x = parameters[2];
518 T gamma_dot_y = parameters[3];
519 T slope_x = parameters[4];
520 T slope_y = parameters[5];
521 T tau = parameters[6];
523 T C = std::exp(- tau * tau / ac_length / ac_length);
529 T mu_kirchhoff_x = mu_gamma * gamma_dot_x;
530 T mu_kirchhoff_y = mu_gamma * gamma_dot_y;
532 T vx = std::sin(theta_i) - std::sin(theta_r1) * std::cos(theta_r2);
533 T vy = - std::sin(theta_r1) * std::sin(theta_r2);
534 T vz = - std::cos(theta_i) - std::cos(theta_r1);
536 T F = (1.0 + std::cos(theta_i) * std::cos(theta_r1) - std::sin(theta_i) * std::sin(theta_r1) * std::cos(theta_r2)) /
537 (std::cos(theta_i) * std::cos(theta_i) + std::cos(theta_r1) * std::cos(theta_i));
539 T sigma_kirchhoff_x = std::sqrt((sigma * sigma / ac_length / ac_length) + 0.5 * sigma_gamma * sigma_gamma * gamma_dot_x * gamma_dot_x);
540 T sigma_kirchhoff_y = std::sqrt((sigma * sigma / ac_length / ac_length) + 0.5 * sigma_gamma * sigma_gamma * gamma_dot_y * gamma_dot_y);
542 T vx_k = (vx + mu_kirchhoff_x * vz);
543 T vy_k = (vy + mu_kirchhoff_y * vz);
545 T p_kirchhoff_x = std::sqrt(2.0 * M_PI) * F / vz / std::sqrt(1.0 - C * C) *
546 std::exp(- (vx_k * vx_k - 2.0 * C * vx_k * (slope_x - mu_kirchhoff_x) * vz + vz * vz * (slope_x - mu_kirchhoff_x) * (slope_x - mu_kirchhoff_x))
547 / (4.0 * vz * vz * sigma_kirchhoff_x * sigma_kirchhoff_x * (1.0 - C * C)));
548 T p_kirchhoff_y = std::sqrt(2.0 * M_PI) * F / vz / std::sqrt(1.0 - C * C) *
549 std::exp(- (vy_k * vy_k - 2.0 * C * vy_k * (slope_y - mu_kirchhoff_y) * vz + vz * vz * (slope_y - mu_kirchhoff_y) * (slope_y - mu_kirchhoff_y))
550 / (4.0 * vz * vz * sigma_kirchhoff_y * sigma_kirchhoff_y * (1.0 - C * C)));
552 return p_kirchhoff_x * p_kirchhoff_y * std::sin(theta_r1);
559 if(theta_r1 > M_PI / 2.0) {
563 theta_r1 = std::max(EPSILON, theta_r1);
564 T eta = 1.0 / std::tan(theta_r1);
565 Matrix<T> w_2d = (sigma_2d * sigma_2d / ac_length / ac_length * 2.0 + sigma_gamma_2d * sigma_gamma_2d * gamma_dot_2d * gamma_dot_2d) ^ 0.5;
566 Matrix<T> Delta = (w_2d / std::sqrt(2.0 * M_PI) * exp(((mu_gamma_2d * gamma_dot_2d * (-1.0) * direction + eta) ^ 2.0) / ((w_2d * w_2d) * 2.0) * (-1.0)) -
567 (mu_gamma_2d * gamma_dot_2d * (-1.0) * direction + eta) * 0.5 * (erfc((mu_gamma_2d * gamma_dot_2d * (-1.0) * direction + eta) / (w_2d * std::sqrt(2.0)))));
568 T exponent = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * Delta) / eta;
569 T base = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
570 base = std::min(1.0, base);
571 return std::pow(base, exponent);
577 theta_r1 = std::max(1e-12, theta_r1);
578 T eta = 1.0 / std::tan(theta_r1);
579 T height_new = height + distance * eta;
580 Matrix<T> w_2d = (sigma_2d * sigma_2d / ac_length / ac_length * 2.0 + sigma_gamma_2d * sigma_gamma_2d * gamma_dot_2d * gamma_dot_2d) ^ 0.5;
581 Matrix<T> Delta = (w_2d / std::sqrt(2.0 * M_PI) * exp(((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) ^ 2.0) / ((w_2d * w_2d) * 2.0) * (-1.0)) -
582 (mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) * 0.5 * (erfc((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) / (w_2d * std::sqrt(2.0)))));
583 T exponent = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * Delta) / eta;
584 T base_old = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
585 T base_new = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height_new) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
587 return std::pow(base_old / base_new, exponent);
593 T theta_r1 = std::max(parameters[0], 1e-12);
594 T height_0 = parameters[1];
595 T gamma_0 = parameters[2];
596 T gamma_dot_0 = parameters[3];
598 T eta = 1.0 / std::tan(theta_r1);
599 T height_new = height_0 + eta * distance;
602 Matrix<T> w_2d = (sigma_2d * sigma_2d / ac_length / ac_length * 2.0 + sigma_gamma_2d * sigma_gamma_2d * gamma_dot_2d * gamma_dot_2d) ^ 0.5;
603 Matrix<T> Delta = (w_2d / std::sqrt(2.0 * M_PI) * exp(((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) ^ 2.0) / ((w_2d * w_2d) * 2.0) * (-1.0)) -
604 (mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) * 0.5 * (erfc((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) / (w_2d * std::sqrt(2.0)))));
605 T exponent = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * Delta) / eta;
606 T base_old = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height_0) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
607 T base_new = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height_new) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
608 T p_h = 1.0 / std::sqrt(2.0 * M_PI) / sigma_0 * std::exp(- (height_new - mu_0) * (height_new - mu_0) / 2.0 / sigma_0 / sigma_0);
609 T S = std::pow(base_old, exponent) / std::pow(base_new, exponent + 1.0);
617 T theta_i = parameters[0];
618 T gamma_0 = parameters[1];
619 T gamma_dot_0 = parameters[2];
624 T p_h = 1.0 / std::sqrt(2.0 * M_PI) / sigma_0 * std::exp(- std::pow(height_0 - mu_0, 2.0) / 2.0 / sigma_0 / sigma_0);
626 return p_h * shadow(theta_i, height_0, -1);
632 Matrix<T> unit_incident = global_vector / global_vector.
norm();
634 z_local(0, 0) = std::sin(theta_r1) * std::cos(theta_r2) - std::sin(theta_i);
635 z_local(1, 0) = std::sin(theta_r1) * std::sin(theta_r2);
636 z_local(2, 0) = std::cos(theta_r1) + std::cos(theta_i);
637 z_local /= z_local.
norm();
639 x_local = unit_incident - z_local * unit_incident.
tr().
dot(z_local)(0, 0);
640 x_local /= x_local.
norm();
642 y_local = x_local.
cross(z_local);
643 y_local /= y_local.
norm();
645 T res_x = global_vector.
tr().
dot(x_local)(0, 0);
646 T res_y = global_vector.
tr().
dot(y_local)(0, 0);
647 T res_z = global_vector.
tr().
dot(z_local)(0, 0);
649 Matrix<T> res(3, 1, {res_x, res_y, res_z});
657 Matrix<T> unit_incident(3, 1, {std::sin(theta_i), 0.0, -std::cos(theta_i)});
659 z_local(0, 0) = std::sin(theta_r1) * std::cos(theta_r2) - std::sin(theta_i);
660 z_local(1, 0) = std::sin(theta_r1) * std::sin(theta_r2);
661 z_local(2, 0) = std::cos(theta_r1) + std::cos(theta_i);
662 z_local /= z_local.
norm();
664 x_local = unit_incident - z_local * unit_incident.
tr().
dot(z_local)(0, 0);
665 x_local /= x_local.
norm();
667 y_local = x_local.
cross(z_local);
668 y_local /= y_local.
norm();
670 return x_local * local_vector(0, 0) + y_local + local_vector(1, 0) + z_local * local_vector(2, 0);
676 Matrix<T> unit_incident(3, 1, {std::sin(theta_i), 0.0, -std::cos(theta_i)});
677 Matrix<T> unit_reflected(3, 1, {std::sin(theta_r1) * std::cos(theta_r2), std::sin(theta_r1) * std::sin(theta_r2), std::cos(theta_r1)});
679 normal = unit_reflected - unit_incident;
680 normal /= std::abs(normal(2, 0));
688 Matrix<T> incident_velocity_copy = this->incident_velocity;
690 T incident_norm = incident_velocity_copy.
norm();
691 T theta_i = std::acos(- incident_velocity_copy(2, 0) / incident_norm);
692 T theta_r1 = 0.0, theta_r2 = 0.0, theta_r2_total = 0.0, distance = 1e8;
693 T gamma_0 = 0.0, gamma_dot_0_x = 0.0, gamma_dot_0_y = 0.0, slope_x = 0.0, slope_y = 0.0, slope_x_new =0.0, slope_y_new = 0.0;
696 bool collided =
true;
698 std::vector<Matrix<T>> p_positions;
699 std::vector<Matrix<T>> p_velocities;
701 std::random_device rd;
702 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
703 std::mt19937 rng(fullSeed);
704 std::normal_distribution<T> normDist(0.0f, 1.0f);
706 T gamma, gamma_dot_x, gamma_dot_y;
707 T domain[] = {- margin, margin};
708 std::tie(gamma, gamma_dot_x) = random_sampler_2D(&
PG_Kernel::mixture_pdf, {1.0, theta_i, 0.0}, normDist(rng), normDist(rng), 0.5, 0.5, domain, domain);
709 gamma_dot_y = normDist(rng) * std::sqrt(2.0) / ac_length;
711 T height_domain[] = {-margin, margin};
712 T height_0 = random_sampler_1D(&
PG_Kernel::height_0_pdf, {theta_i, gamma, gamma_dot_x}, normDist(rng), 0.5, height_domain);
713 p_positions.push_back(
Matrix<T>(3, 1, {height_0, 0.0, 0.0}));
717 T angle_domain_x[] = {0.0, M_PI - theta_i};
718 T angle_domain_y[] = {-M_PI, M_PI};
719 std::tie(theta_r1, theta_r2) = random_sampler_angles(&
PG_Kernel::kirchhoff_pdf, {theta_i, gamma, gamma_dot_x, gamma_dot_y, slope_x, slope_y, distance}, 1e-5, 1e-5,
720 (M_PI - theta_i) * 0.25, (M_PI - theta_i) * 0.25, angle_domain_x, angle_domain_y);
721 Matrix<T> slopes_matrix = angles_to_slopes(theta_i, theta_r1, theta_r2);
722 slope_x_new = - slopes_matrix(0, 0);
723 slope_y_new = - slopes_matrix(1, 0);
725 Matrix<T> incident_velocity_local = global_to_local(incident_velocity_copy, theta_i, theta_r1, theta_r2);
726 T vel_in_array[] = {incident_velocity_local(0, 0), incident_velocity_local(1, 0), incident_velocity_local(2, 0)};
727 T vel_refl_array[] = {0.0, 0.0, 0.0};
732 local_kernel(vel_in_array, vel_refl_array, gsi_parameters, gas_parameters, surface_parameters);
734 Matrix<T> reflected_velocity_local(3, 1, {vel_refl_array[0], vel_refl_array[1], vel_refl_array[2]});
735 Matrix<T> reflected_velocity = local_to_global(reflected_velocity_local, theta_i, theta_r1, theta_r2);
737 theta_r1 = std::acos(reflected_velocity(2, 0) / reflected_velocity.
norm());
738 theta_r2 = std::atan2(reflected_velocity(1, 0), reflected_velocity(0, 0));
740 std::bernoulli_distribution shadow_pdf(1.0 - shadow(theta_r1, height_0, 1));
741 collided = shadow_pdf(rng);
743 T distance_domain[] = {0.0, margin * 5.0};
744 distance = random_sampler_1D(&
PG_Kernel::distance_pdf, {theta_r1, height_0, gamma, gamma_dot_x}, 0.0, 0.5, distance_domain);
745 height_0 += distance / std::tan(theta_r1);
746 theta_i = M_PI - theta_r1;
747 theta_r2_total += theta_r2;
749 incident_norm = reflected_velocity.
norm();
750 incident_velocity_copy(0, 0) = incident_norm * std::sin(theta_i);
751 incident_velocity_copy(1, 0) = 0.0;
752 incident_velocity_copy(2, 0) = - incident_norm * std::cos(theta_i);
754 Matrix<T> global_reflected_velocity(3, 1, {std::sin(theta_r1) * std::cos(theta_r2_total),
755 std::sin(theta_r1) * std::sin(theta_r2_total),
756 std::cos(theta_r1)});
757 global_reflected_velocity *= incident_norm;
758 p_velocities.push_back(global_reflected_velocity);
760 Matrix<T> position_change(3, 1, {std::sin(theta_r1) * std::cos(theta_r2_total),
761 std::sin(theta_r1) * std::sin(theta_r2_total), std::cos(theta_r1)});
762 position_change *= distance / std::sin(theta_r1);
763 if (collided) p_positions.push_back(p_positions.back() + position_change);
765 slope_x = - std::cos(theta_r2) * slope_x_new - std::sin(theta_r2) * slope_y_new;
766 slope_y = std::sin(theta_r2) * slope_x_new - std::cos(theta_r2) * slope_y_new;
768 std::tie(gamma, gamma_dot_x) = random_sampler_2D(&
PG_Kernel::mixture_pdf, {0.0, theta_i, distance}, normDist(rng), normDist(rng), 0.5, 0.5, domain, domain);
769 gamma_dot_y = normDist(rng) * std::sqrt(2.0) / ac_length;
780 std::vector<trajectory<T>> trajectories;
782 for(
auto i = 0; i < this->num_particles; i ++) {
784 trajectories.push_back(p_trajectory);
792 std::ofstream myfile;
793 myfile.open(filename);
795 unsigned long num_particles = trajectories.size();
797 for(
auto i = 0; i < num_particles; i ++) {
800 unsigned long num_bounces = p_trajectory.positions.size();
802 myfile << num_bounces <<
"\n";
804 for(
auto j = 0; j < 3; j ++) {
805 for(
auto k = 0; k < num_bounces; k ++) {
806 myfile << p_trajectory.positions[k](j, 0) <<
" ";
811 for(
auto j = 0; j < 3; j ++) {
812 for(
auto k = 0; k < num_bounces; k ++) {
813 myfile << p_trajectory.velocities[k](j, 0) <<
" ";
826 std::ifstream linecounter(filename);
828 std::vector<trajectory<T>> trajectories;
832 while( std::getline(linecounter, line)) {
836 long num_particles = num_lines / 7;
839 std::ifstream myfile(filename);
841 for(
long i = 0; i < num_particles; i ++) {
843 myfile >> num_bounces;
846 std::vector<Matrix<T>> p_positions(num_bounces,
Matrix<T>(3, 1));
847 std::vector<Matrix<T>> p_velocities(num_bounces,
Matrix<T>(3, 1));
849 for(
long k = 0; k < 3; k ++) {
850 for(
long j = 0; j < num_bounces; j ++) {
851 myfile >> p_positions[j](k, 0);
854 for(
long k = 0; k < 3; k ++) {
855 for(
long j = 0; j < num_bounces; j ++) {
856 myfile >> p_velocities[j](k, 0);
869 std::vector<trajectory<T>> trajectories_total;
871 for(
long i = 0; i < num_procs; i ++) {
872 std::vector<trajectory<T>> trajectories = import_data(filename +
"_" + std::to_string(i) +
".dat");
873 trajectories_total.insert(trajectories_total.end(), trajectories.begin(), trajectories.end());
874 std::remove((filename +
"_" + std::to_string(i) +
".dat").c_str());
876 save(trajectories_total, filename +
".dat");
Class representing a Gas with properties such as density, temperature, molar mass,...
Definition Gas.h:10
T get_molar_mass() const
Returns the molar mass of the gas.
Definition Gas.h:141
T get_density() const
Returns the density of the gas.
Definition Gas.h:134
T get_temperature() const
Returns the temperature of the gas.
Definition Gas.h:127
T get_speed() const
Returns the speed of the gas.
Definition Gas.h:148
Class for computing Hermite polynomials and related operations.
Definition Hermite_Tools.h:11
Matrix< T > evaluate(Matrix< T > coefficients, Matrix< T > values)
Evaluates the Hermite polynomial at specified values.
Definition Hermite_Tools.h:136
Matrix< T > get_roots(unsigned int order)
Retrieves the roots of the Hermite polynomial of specified order.
Definition Hermite_Tools.h:276
Matrix< T > get_weights(unsigned int order)
Retrieves the weights of the Hermite polynomial of specified order.
Definition Hermite_Tools.h:299
Matrix< T > evaluate_deriv(Matrix< T > coefficients, Matrix< T > values)
Evaluates the derivative of the Hermite polynomial at specified values.
Definition Hermite_Tools.h:215
This class implements a matrix object used for linear algebra and vectorized operations.
Definition Matrix.h:25
Matrix cross(const Matrix &other)
Calculates the cross product of two matrices.
Definition Matrix.h:624
std::tuple< Matrix< T >, Matrix< T > > meshgrid(const Matrix &other)
Calculates the meshgrid of 2 1D matrices.
Definition Matrix.h:656
Matrix dot(const Matrix &other)
Calculates the dot product of two matrices.
Definition Matrix.h:604
T norm()
Calculates the Frobenius norm of the matrix.
Definition Matrix.h:674
Matrix outer(const Matrix &other)
Calculates the outer product of two matrices.
Definition Matrix.h:639
Matrix tr()
Transposes the matrix.
Definition Matrix.h:705
Represents a class for the proposed rough surface kernel applied to Poly-Gaussian surfaces.
Definition PG_Kernel.h:26
trajectory< T > sample_particle()
Sample a single particle trajectory.
Definition PG_Kernel.h:686
void set_num_particles(unsigned long num_particles)
Set the number of particles.
Definition PG_Kernel.h:347
T distance_pdf(T distance, std::vector< T > parameters)
Calculate the PDF for the distance parameter.
Definition PG_Kernel.h:591
unsigned long get_num_particles()
Get the number of particles.
Definition PG_Kernel.h:369
void save(std::vector< trajectory< T > > trajectories, std::string filename)
Save the trajectories to a file.
Definition PG_Kernel.h:790
PG_Kernel(Surface< T > surface, Gas< T > gas, void(*local_kernel)(T *, T *, T *, T *, T *), Matrix< T > incident_velocity, unsigned long num_particles, std::string sim_name)
Constructor for PG_Kernel.
Definition PG_Kernel.h:101
std::vector< trajectory< T > > import_data(std::string filename)
Import trajectory data from a file.
Definition PG_Kernel.h:824
T mixture_pdf(T gamma, T gamma_dot, std::vector< T > parameters)
Calculate the PDF of the control process gamma and its derivative gamma_dot, for Poly-Gaussian surfac...
Definition PG_Kernel.h:485
Matrix< T > local_to_global(Matrix< T > local_vector, T theta_i, T theta_r1, T theta_r2)
Convert local coordinates to global coordinates.
Definition PG_Kernel.h:655
Matrix< T > angles_to_slopes(T theta_i, T theta_r1, T theta_r2)
Convert angles to slopes.
Definition PG_Kernel.h:674
std::string & get_sim_name()
Get the simulation name.
Definition PG_Kernel.h:362
std::tuple< T, T > random_sampler_angles(T(PG_Kernel::*pdf)(T, T, std::vector< T >), std::vector< T > parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[])
Perform random sampling of angles in 2D using the Metropolis-Hastings algorithm.
Definition PG_Kernel.h:440
T shadow(T theta_r1, T height, int direction)
Calculate the probability of a reflected particle to escape the surface.
Definition PG_Kernel.h:556
T random_sampler_1D(T(PG_Kernel::*pdf)(T, std::vector< T >), std::vector< T > parameters, T x0, T size, T domain_x[])
Perform random sampling in 1D using the Metropolis-Hastings algorithm.
Definition PG_Kernel.h:375
T occlusion(T theta_r1, T height, T distance, int direction)
Calculate the probability of a reflected particle to escape the surface, given a certain travelled ho...
Definition PG_Kernel.h:575
~PG_Kernel()
Destructor for PG_Kernel.
Definition PG_Kernel.h:153
T height_0_pdf(T height_0, std::vector< T > parameters)
Calculate the PDF for the height_0 parameter.
Definition PG_Kernel.h:615
Matrix< T > global_to_local(Matrix< T > global_vector, T theta_i, T theta_r1, T theta_r2)
Convert global coordinates to local coordinates.
Definition PG_Kernel.h:630
std::vector< trajectory< T > > sample_batch()
Sample a batch of particle trajectories.
Definition PG_Kernel.h:778
std::tuple< T, T > random_sampler_2D(T(PG_Kernel::*pdf)(T, T, std::vector< T >), std::vector< T > parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[])
Perform random sampling in 2D using the Metropolis-Hastings algorithm.
Definition PG_Kernel.h:405
Surface< T > & get_surface()
Get the Surface object.
Definition PG_Kernel.h:177
T kirchhoff_pdf(T theta_r1, T theta_r2, std::vector< T > parameters)
Calculate the reflected particle angular PDF according to the Kirchhoff model.
Definition PG_Kernel.h:513
Gas< T > & get_gas()
Get the Gas object.
Definition PG_Kernel.h:170
PG_Kernel()
Default constructor for PG_Kernel.
Definition PG_Kernel.h:76
void combine_files(std::string filename, int num_procs)
Combine multiple trajectory files into one.
Definition PG_Kernel.h:867
A class representing a surface with geometry and statistical properties.
Definition Surface.h:11
T autocorrelation_length()
Retrieves the autocorrelation length of the surface.
Definition Surface.h:100
T * get_local_parameters()
Retrieves the local parameters of the surface.
Definition Surface.h:107
Matrix< T > get_sigma_coefficients()
Retrieves the sigma coefficients for the surface.
Definition Surface.h:86
Matrix< T > get_mu_coefficients()
Retrieves the mu coefficients for the surface.
Definition Surface.h:79
Represents a trajectory of a certain type.
Definition PG_Kernel.h:10