From 3b1f0bf1fae0188b61a70ce2d17bec33a3f27db4 Mon Sep 17 00:00:00 2001
From: Bram Veenboer <bram.veenboer@gmail.com>
Date: Fri, 24 Jan 2020 14:51:48 +0100
Subject: [PATCH] Add oskar spherical wave model
The implementation is not complete, the coeffficients are missing.
---
Station.cc | 6 +
oskar/CMakeLists.txt | 5 +-
oskar/OSKARElementResponse.cc | 24 ++++
oskar/OSKARElementResponse.h | 10 ++
oskar/oskar_evaluate_dipole_pattern.cc | 1 -
oskar/oskar_evaluate_spherical_wave_sum.cc | 131 +++++++++++++++++++++
6 files changed, 175 insertions(+), 2 deletions(-)
create mode 100644 oskar/oskar_evaluate_spherical_wave_sum.cc
diff --git a/Station.cc b/Station.cc
index 8395b8ca..f5b6d1ea 100644
--- a/Station.cc
+++ b/Station.cc
@@ -60,8 +60,14 @@ Station::Station(
}
break;
case OSKAR:
+ std::clog << "initialize oskar!" << std::endl;
+ #if 0
DualDipoleAntenna::itsElementResponse.reset(new OSKARElementResponseDipole());
TileAntenna::itsElementResponse.reset(new OSKARElementResponseDipole());
+ #else
+ DualDipoleAntenna::itsElementResponse.reset(new OSKARElementResponseSphericalWave());
+ TileAntenna::itsElementResponse.reset(new OSKARElementResponseSphericalWave());
+ #endif
break;
default:
std::stringstream message;
diff --git a/oskar/CMakeLists.txt b/oskar/CMakeLists.txt
index 33a47ca7..e1cd3d46 100644
--- a/oskar/CMakeLists.txt
+++ b/oskar/CMakeLists.txt
@@ -1,5 +1,8 @@
+# build liboskar.so
add_library(oskar SHARED
OSKARElementResponse.cc
- oskar_evaluate_dipole_pattern.cc)
+ oskar_evaluate_dipole_pattern.cc
+ oskar_evaluate_spherical_wave_sum.cc)
+# install liboskar.so
install(TARGETS oskar DESTINATION lib)
\ No newline at end of file
diff --git a/oskar/OSKARElementResponse.cc b/oskar/OSKARElementResponse.cc
index 390c05e9..eff96442 100644
--- a/oskar/OSKARElementResponse.cc
+++ b/oskar/OSKARElementResponse.cc
@@ -27,3 +27,27 @@ void OSKARElementResponseDipole::element_response(
oskar_evaluate_dipole_pattern_double(1, &theta, &phi_x, freq, dipole_length_m, response_ptr);
oskar_evaluate_dipole_pattern_double(1, &theta, &phi_y, freq, dipole_length_m, response_ptr + 2);
}
+
+void oskar_evaluate_spherical_wave_sum_double(
+ const int num_points,
+ const double* theta,
+ const double* phi_x,
+ const double* phi_y,
+ const int l_max,
+ const std::complex<double>* alpha,
+ std::complex<double>* pattern);
+
+void OSKARElementResponseSphericalWave::element_response(
+ double freq,
+ double theta,
+ double phi,
+ std::complex<double> (&response)[2][2]) const
+{
+ int l_max = 1; // TODO
+ std::complex<double>* response_ptr = (std::complex<double> *) response;
+ std::complex<double>* alpha_ptr; // TODO
+
+ double phi_x = phi;
+ double phi_y = phi + M_PI/2;
+ oskar_evaluate_spherical_wave_sum_double(1, &theta, &phi_x, &phi_y, l_max, alpha_ptr, response_ptr);
+}
\ No newline at end of file
diff --git a/oskar/OSKARElementResponse.h b/oskar/OSKARElementResponse.h
index 5fa12389..80125a93 100644
--- a/oskar/OSKARElementResponse.h
+++ b/oskar/OSKARElementResponse.h
@@ -15,3 +15,13 @@ public:
double phi,
std::complex<double> (&response)[2][2]) const final override;
};
+
+class OSKARElementResponseSphericalWave : public OSKARElementResponse
+{
+public:
+ virtual void element_response(
+ double freq,
+ double theta,
+ double phi,
+ std::complex<double> (&response)[2][2]) const final override;
+};
\ No newline at end of file
diff --git a/oskar/oskar_evaluate_dipole_pattern.cc b/oskar/oskar_evaluate_dipole_pattern.cc
index 68adae88..0debc51f 100644
--- a/oskar/oskar_evaluate_dipole_pattern.cc
+++ b/oskar/oskar_evaluate_dipole_pattern.cc
@@ -28,7 +28,6 @@
#include <cmath>
#include <complex>
-#include <iostream>
#include "oskar_vector_types.h"
#include "oskar_helper.h"
diff --git a/oskar/oskar_evaluate_spherical_wave_sum.cc b/oskar/oskar_evaluate_spherical_wave_sum.cc
new file mode 100644
index 00000000..2cb2b4c2
--- /dev/null
+++ b/oskar/oskar_evaluate_spherical_wave_sum.cc
@@ -0,0 +1,131 @@
+/*
+ * Copyright (c) 2019, The University of Oxford
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of the University of Oxford nor the names of its
+ * contributors may be used to endorse or promote products derived from this
+ * software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <cmath>
+#include <complex>
+
+#include "oskar_vector_types.h"
+#include "oskar_helper.h"
+
+template<typename FP, typename FP2, typename FP4c>
+void oskar_evaluate_spherical_wave_sum(
+ int num_points, const FP* theta,
+ const FP* phi_x, const FP* phi_y, int l_max,
+ const FP4c* alpha, int offset, FP4c* pattern)
+{
+ FP2 Xp, Xt, Yp, Yt;
+ make_zero2(Xp); make_zero2(Xt);
+ make_zero2(Yp); make_zero2(Yt);
+
+ #pragma omp parallel for
+ for (int i = 0; i < num_points; i++) {
+ const FP phi_x_ = phi_x[i], theta_ = theta[i];
+ const FP phi_y_ = phi_y[i];
+ /* Propagate NAN. */
+ if (phi_x_ != phi_x_) {
+ Xp.x = Xp.y = Xt.x = Xt.y = phi_x_;
+ Yp.x = Yp.y = Yt.x = Yt.y = phi_x_;
+ }
+ else {
+ FP sin_t, cos_t;
+ oskar_sincos(theta_, &sin_t, &cos_t);
+ for (int l = 1; l <= l_max; ++l) {
+ const int ind0 = l * l - 1 + l;
+ const FP f_ = (2 * l + 1) / (4 * ((FP)M_PI) * l * (l + 1));
+ for (int abs_m = l; abs_m >=0; --abs_m) {
+ FP p, pds, dpms, sin_p, cos_p;
+ oskar_legendre2(l, abs_m, cos_t, sin_t, p, pds, dpms);
+ if (abs_m == 0) {
+ sin_p = (FP)0; cos_p = sqrt(f_);
+ const FP4c alpha_ = alpha[ind0];
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.a, alpha_.b, Xt, Xp);
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.c, alpha_.d, Yt, Yp);
+ }
+ else {
+ FP d_fact = (FP)1, s_fact = (FP)1;
+ const int d_ = l - abs_m, s_ = l + abs_m;
+ for (int i_ = 2; i_ <= d_; ++i_) d_fact *= i_;
+ for (int i_ = 2; i_ <= s_; ++i_) s_fact *= i_;
+ const FP ff = f_ * d_fact / s_fact;
+ const FP nf = sqrt(ff);
+ const FP4c alpha_m = alpha[ind0 - abs_m];
+ const FP4c alpha_p = alpha[ind0 + abs_m];
+ p = -abs_m * phi_x_;
+ oskar_sincos(p, &sin_p, &cos_p);
+ sin_p *= nf; cos_p *= nf;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, -abs_m, alpha_m.a, alpha_m.b, Xt, Xp);
+ sin_p = -sin_p;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, abs_m, alpha_p.a, alpha_p.b, Xt, Xp);
+ p = -abs_m * phi_y_;
+ oskar_sincos(p, &sin_p, &cos_p);
+ sin_p *= nf; cos_p *= nf;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, -abs_m, alpha_m.c, alpha_m.d, Yt, Yp);
+ sin_p = -sin_p;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, abs_m, alpha_p.c, alpha_p.d, Yt, Yp);
+ }
+ }
+ }
+ }
+ pattern[i + offset].a = Xt;
+ pattern[i + offset].b = Xp;
+ pattern[i + offset].c = Yt;
+ pattern[i + offset].d = Yp;
+ }
+}
+
+void oskar_evaluate_spherical_wave_sum_double(
+ const int num_points,
+ const double* theta,
+ const double* phi_x,
+ const double* phi_y,
+ const int l_max,
+ const std::complex<double>* alpha,
+ std::complex<double>* pattern)
+{
+ double4c* alpha_ptr = (double4c *) alpha;
+ double4c* pattern_ptr = (double4c *) pattern;
+
+ oskar_evaluate_spherical_wave_sum<double, double2, double4c>(
+ num_points, theta, phi_x, phi_y, l_max, alpha_ptr, 0, pattern_ptr);
+}
+
+void oskar_evaluate_spherical_wave_sum_float(
+ const int num_points,
+ const float* theta,
+ const float* phi_x,
+ const float* phi_y,
+ const int l_max,
+ const std::complex<float>* alpha,
+ std::complex<float>* pattern)
+{
+ float4c* alpha_ptr = (float4c *) alpha;
+ float4c* pattern_ptr = (float4c *) pattern;
+
+ oskar_evaluate_spherical_wave_sum<float, float2, float4c>(
+ num_points, theta, phi_x, phi_y, l_max, alpha_ptr, 0, pattern_ptr);
+}
--
GitLab