From dfae9d45af4551e38581136b3d397204f1a4ea35 Mon Sep 17 00:00:00 2001
From: Tilo Arens <tilo.arens@kit.edu>
Date: Wed, 18 Jan 2023 09:15:12 +0100
Subject: [PATCH] Program files to produce the plots added to repository.
---
src/indicator_func.m | 70 +++++++++++++++++++
src/operator_N.m | 148 +++++++++++++++++++++++++++++++++++++++++
src/operator_T_B.m | 118 ++++++++++++++++++++++++++++++++
src/plot_of_figure_1.m | 52 +++++++++++++++
src/plot_of_figure_2.m | 49 ++++++++++++++
src/plot_of_figure_3.m | 52 +++++++++++++++
6 files changed, 489 insertions(+)
create mode 100644 src/indicator_func.m
create mode 100644 src/operator_N.m
create mode 100644 src/operator_T_B.m
create mode 100644 src/plot_of_figure_1.m
create mode 100644 src/plot_of_figure_2.m
create mode 100644 src/plot_of_figure_3.m
diff --git a/src/indicator_func.m b/src/indicator_func.m
new file mode 100644
index 0000000..e34c7cd
--- /dev/null
+++ b/src/indicator_func.m
@@ -0,0 +1,70 @@
+function [indicat, Xi1, Xi2] = indicator_func(op_N, alpha, k_out, H, R, x_max, mesh_step, a, q_positive)
+%
+% Compute the indicator function for finding an obstacle for a transmission
+% problem in a Neumann waveguide of height H. The wave number in the
+% waveguide is k_out.
+%
+% The near field measurements, which are taken on vertikal lines at +-R,
+% are contained in the matrix operator_N.
+%
+% The function creates a grid of points on [-x_max, xmax] x [0, H] with
+% mesh size mesh_step. For each grid point, the operator T_B is computed for
+% B equal to a square with lateral length 2a centered at the grid point.
+% The number of negative eigenvalues of the matrix
+%
+% Re(N) - alpha T_B
+%
+% is computed. This the value of the indicator for this point.
+%
+% The parameter q_positive should be true for the case q > 0, it should be
+% false for the case -1 < q <= 0.
+
+% threshold for negative eigenvalues. For now we define this as a constant.
+if q_positive
+ delta = -2e-6; % for figure 1
+ %delta = -2e-6; % for figure 2
+else
+ delta = 2e-6; % for figure 3
+end
+
+% obtain the number of source/receiver points from the dimension of
+% operator_N
+M = size(op_N, 1) / 2 - 1;
+
+% grid in x1 direction: x_max_true is the largest value such that the
+% square centred at that x_1 coordinate is completely inside the rectangle
+% [-x_max, x_max] x [0, H]
+N_x1 = floor(2*(x_max-a) / mesh_step);
+x_max_true = mesh_step * N_x1 / 2;
+x1 = linspace(-x_max_true, x_max_true, N_x1+1);
+
+% grid in x2 direction. The centers of the first and last pixels are
+% x2_offset away from the waveguide boundary.
+bound_gap = 0.05;
+N_x2 = floor((H-2*a-2*bound_gap) / mesh_step);
+x2_offset = (H - N_x2*mesh_step) / 2;
+x2 = linspace(x2_offset, H-x2_offset, N_x2+1);
+
+[Xi1, Xi2] = meshgrid(x1, x2);
+indicat = zeros(size(Xi1));
+
+real_N = (op_N + op_N') / 2;
+
+for j=1:N_x1+1
+ fprintf('.');
+ for l=1:N_x2+1
+ T_B = operator_T_B(H, k_out, R, Xi1(l, j), Xi2(l, j), a, M);
+ the_op = real_N - alpha * T_B;
+ eigenvals = eig(the_op);
+ if q_positive
+ indicat(l, j) = sum(eigenvals < delta);
+ else
+ indicat(l, j) = sum(eigenvals > delta);
+ end
+ end
+end
+fprintf('\n\n');
+
+
+end
+
diff --git a/src/operator_N.m b/src/operator_N.m
new file mode 100644
index 0000000..b41dc40
--- /dev/null
+++ b/src/operator_N.m
@@ -0,0 +1,148 @@
+function N_mat_ct = operator_N(H, curve, k_out, k_in, R, M, N, N_quad_ct, debug_mode)
+%
+% The function generates a discretization for the near field operator N for
+% a penetrable obstacle in a waveguide with Neumann boundary conditions.
+%
+% The waveguide is \R \times (0, 'H'). The obstacle is given by the boundary
+% curve 'curve' and should be a BIEPack.geom.ParamCurve instance. It must
+% be contained in the waveguide.
+%
+% The numbers 'k_out' and 'k_in' are the wavenumbers outside and inside the
+% obstacle, respectively. The positive number 'R' specifes the
+% x_1-coordinate of the vertical lines containing the source and receiver
+% points.
+%
+% The space L^2(C_R) is discretized using the basis induced by the Neumann
+% eigenfunctions of the Laplacian on (0, H). If one of these are used as
+% the density in the conjugated single layer potential, the series for the
+% Green's function reduces to just one term. We use the first M+1
+% eigenfunctions, hence the returned matrix will be of size 2M+2 x 2M+2.
+%
+% To compute the coefficients of the scattered field with respect to the
+% basis, the composite trapezoidal rule with M+1 points is used on each
+% vertical line. As the field can be extended as an even function for
+% negative x_2 to a 2H-periodic function, the rule is highly accurate.
+%
+% The parameter 'N' specifies the number of unknowns to use in the boundary
+% integral equation method.
+%
+% The parameter debug_mode is optional and defaults to false. If it is
+% true, some tests whether the boundary data is correct, are carried out.
+
+% set default value of parameters
+if (nargin < 8)
+ N_quad_ct = 80;
+end
+
+if (nargin < 9)
+ debug_mode = false;
+end
+
+% shorten names of BIEPack classes
+import BIEPack.*
+
+% allocate memory for the output and define source/receiver positions
+%N_mat_gl = zeros(2*M+2, 2*M+2);
+N_mat_ct = zeros(2*M+2, 2*M+2);
+
+% quadrature points on C_R and corresponding weights for evaluating scalar
+% products with basis functions
+y_quad_ct = linspace(0, H, N_quad_ct+1).';
+w_ct = H/N_quad_ct * ones(1, N_quad_ct+1);
+w_ct([1, N_quad_ct+1]) = 0.5 * w_ct([1, N_quad_ct+1]);
+
+%[y_quad_gl, w_gl] = BIEPack.utils.quadGaussLegendre(40, 0, H);
+%y_quad_gl = y_quad_gl.';
+
+%y1_gl = [-R * ones(size(y_quad_gl)); R * ones(size(y_quad_gl))];
+%y2_gl = [y_quad_gl; y_quad_gl];
+y1_ct = [-R * ones(size(y_quad_ct)); R * ones(size(y_quad_ct))];
+y2_ct = [y_quad_ct; y_quad_ct];
+
+
+% define the mesh
+mesh = meshes.ParamCurveMesh(curve, N);
+space = spaces.TrigPol(mesh);
+
+% define all the integral operators
+slop_ext = ops.HelmholtzNwgSLOp(space, space, k_out, H);
+dlop_ext = ops.HelmholtzNwgDLOp(space, space, k_out, H);
+adlop_ext = ops.HelmholtzNwgADLOp(space, space, k_out, H);
+slop_int = ops.HelmholtzSLOp(space, space, k_in);
+dlop_int = ops.HelmholtzDLOp(space, space, k_in);
+adlop_int = ops.HelmholtzADLOp(space, space, k_in);
+hsing_dif_op = ops.HelmholtzNwgHsingDifOp(space, space, k_out, k_in, H);
+slpot = ops.HelmholtzNwgSLPot(space, space, k_out, H);
+dlpot = ops.HelmholtzNwgDLPot(space, space, k_out, H);
+
+% Define the block operator A in the IE system (I - A) phi = psi in
+% the Kress-Roach integral equation formulation for a transmission problem.
+A = {(-1)*dlop_ext + dlop_int, slop_ext + (-1)*slop_int;
+ (-1)*hsing_dif_op, adlop_ext + (-1)*adlop_int};
+
+% Set up the Nyström method for solving the IE system
+nystroemMethod = methods.NystroemPeriodicLog();
+nystroemMethod.setOperators(A);
+
+if (debug_mode)
+ % Operators for checking that the boundary data is sensible. It
+ % corresponds to Cauchy data of a solution of the Helmholtz equation in
+ % the obstacle, so it should be in the kernel of the Calderon projector.
+ slop_ff_out = ops.HelmholtzSLOp(space, space, k_out);
+ dlop_ff_out = ops.HelmholtzDLOp(space, space, k_out);
+ calderon_dir = 2.0*dlop_ff_out;
+ calderon_neu = (-2.0)*slop_ff_out;
+ calderon_dir_mat = calderon_dir.getImplementation('NystroemPeriodicLog');
+ calderon_neu_mat = calderon_neu.getImplementation('NystroemPeriodicLog');
+ fprintf('Exterior wave number is %6.4f\n', k_out);
+end
+
+% for each source point, assign right hand side and solve
+import BIEPack.utils.PeriodicGreensFunc
+for j=1:2*M+2
+ mesh.evalDerivs(1);
+
+ k_n = mod(j-1, M+1) * pi / H;
+ beta_n_conj = sqrt(k_out^2 - k_n^2)';
+
+ if j == 1
+ ui_dirichlet = -0.5i / (sqrt(H)*beta_n_conj) * exp(-1i * beta_n_conj * (R + mesh.nodesX.'));
+ ui_neumann = -0.5 / sqrt(H) * exp(-1i * beta_n_conj * (R + mesh.nodesX.')) .* mesh.normalAtNodeX.';
+ elseif 2 <= j && j <= M+1
+ ui_dirichlet = -0.5i / (sqrt(H/2)*beta_n_conj) * exp(-1i * beta_n_conj * (R + mesh.nodesX.')) .* cos(k_n * mesh.nodesY.');
+ ui_neumann = -0.5 / sqrt(H/2) * exp(-1i * beta_n_conj * (R + mesh.nodesX.')) ...
+ .* (cos(k_n * mesh.nodesY.') .* mesh.normalAtNodeX.' - 1i * k_n / beta_n_conj * sin(k_n * mesh.nodesY.') .* mesh.normalAtNodeY.');
+ elseif j == M+2
+ ui_dirichlet = -0.5i / (sqrt(H)*beta_n_conj) * exp(-1i * beta_n_conj * (R - mesh.nodesX.'));
+ ui_neumann = 0.5 / sqrt(H) * exp(-1i * beta_n_conj * (R - mesh.nodesX.')) .* mesh.normalAtNodeX.';
+ else
+ ui_dirichlet = -0.5i / (sqrt(H/2)*beta_n_conj) * exp(-1i * beta_n_conj * (R - mesh.nodesX.')) .* cos(k_n * mesh.nodesY.');
+ ui_neumann = 0.5 / sqrt(H/2) * exp(-1i * beta_n_conj * (R - mesh.nodesX.')) ...
+ .* (cos(k_n * mesh.nodesY.') .* mesh.normalAtNodeX.' + 1i * k_n / beta_n_conj * sin(k_n * mesh.nodesY.') .* mesh.normalAtNodeY.');
+ end
+
+ if (debug_mode)
+ % check that the boundary data is sensible.
+ residual = ui_dirichlet + calderon_dir_mat * ui_dirichlet + calderon_neu_mat * ui_neumann;
+ fprintf('mode = %2d, k_n = %6.4f, boundary data Calderon residual: %e\n', j, k_n, norm(residual));
+ end
+
+ nystroemMethod.setRhs({-ui_dirichlet; -ui_neumann});
+ phi = nystroemMethod.apply;
+
+ u_scat_ct = dlpot.evalPotential(y1_ct, y2_ct, phi(1:N)) - slpot.evalPotential(y1_ct, y2_ct, phi(N+1:2*N));
+
+ theta_l = 1/sqrt(H) * ones(size(y_quad_ct));
+ N_mat_ct(1, j) = w_ct * (theta_l .* u_scat_ct(1:length(y_quad_ct)));
+ N_mat_ct(M+2, j) = w_ct * (theta_l .* u_scat_ct(length(y_quad_ct)+1:end));
+ for l=2:M+1
+ k_l = (l-1)*pi/H;
+ theta_l = sqrt(2/H) * cos(k_l * y_quad_ct);
+ N_mat_ct(l, j) = w_ct * (theta_l .* u_scat_ct(1:length(y_quad_ct)));
+ N_mat_ct(M+1+l, j) = w_ct * (theta_l .* u_scat_ct(length(y_quad_ct)+1:end));
+ end
+end
+
+% fprintf('Difference %e\n', norm(N_mat_gl - N_mat_ct));
+
+end
\ No newline at end of file
diff --git a/src/operator_T_B.m b/src/operator_T_B.m
new file mode 100644
index 0000000..3db3d2f
--- /dev/null
+++ b/src/operator_T_B.m
@@ -0,0 +1,118 @@
+function T_B = operator_T_B(H, k_out, R, xi1, xi2, a, M)
+%
+% Set up the operator T_B where B is a square centered at (xi_1, xi_2) with
+% lateral length 2a. The operator is discretized by using the basis in
+% L^2(C_R) induced by the Neumann eigenfunctions of the Laplacian,
+%
+% \theta_n(t) = sqrt(2/H) cos( n pi / H * t ) for n \geq 1
+% \theta_0(t) = sqrt(1/H)
+%
+% The corresponding basis is
+%
+% g^1_n = ( theta_n, 0 ) , g^2_n = ( 0, theta_n )
+%
+% where the first component function is defined on C_{-, R}, the second one
+% on C_{+, R}. Fixing the number M of the highest eigen function used, the
+% resulting matrix is of size 2M+2 x 2M+2 and contains the scalar products
+%
+% A^(pq)_{jn} = ( T_B g^q_n , g^p_j )_{L^2(0, H)} p,q = 1,2, j,n = 0,...,M.
+%
+% As the Green's function is also expanded in the theta_n, only one term in
+% the double series in the expression of T_B remains. This still envolves
+% an integral over two exponential functions. This is computed by a
+% Gauss-Legendre quadrature rule.
+
+T_B = zeros(2*M+2, 2*M+2);
+
+k_n = (0:M) * pi / H;
+beta_n = sqrt(k_out^2 - k_n.^2).';
+max_real_index = sum(imag(beta_n) == 0);
+
+% % The number of Gauss-Legendre points for quadratur with respect to y_1
+% N_gauss = 15;
+% [y1, w] = BIEPack.utils.quadGaussLegendre(N_gauss, xi1-a, xi1+a);
+% W = diag(w);
+%
+% pre_factors = 0.25 * k_out^2 ./ beta_n ./ beta_n' .* gamma_n_m_eval(M, k_n, xi2, a, H);
+% exp_plus = exp(1i * beta_n * (R + y1));
+% exp_minus = exp(1i * beta_n * (R - y1));
+%
+% T_B(1:M+1, 1:M+1) = exp_plus * W * exp_plus' .* pre_factors;
+% T_B(1:M+1, M+2:end) = exp_plus * W * exp_minus' .* pre_factors;
+% T_B(M+2:end, 1:M+1) = exp_minus * W * exp_plus' .* pre_factors;
+% T_B(M+2:end, M+2:end) = exp_minus * W * exp_minus' .* pre_factors;
+
+one_vec = ones(size(beta_n));
+gamma_n_m = gamma_n_m_eval(M, k_n, xi2, a, H);
+pre_factors = 0.25 * k_out^2 ./ beta_n ./ beta_n' .* gamma_n_m ...
+ .* (exp(1i * R * (beta_n * one_vec.' - one_vec * beta_n')));
+
+% block nu = mu = 1
+denom = 1i * (beta_n * one_vec.' - one_vec * beta_n');
+exp_plus = exp(denom * (xi1 + a));
+exp_minus = exp(denom * (xi1 - a));
+T_B_block = (exp_plus - exp_minus) ./ denom .* pre_factors;
+diag_indeces = (1:M+2:(M+2)*(max_real_index-1)+1).';
+T_B_block(diag_indeces) = 0.5 * a * k_out^2 ./ beta_n(1:max_real_index).^2 .* gamma_n_m(diag_indeces);
+T_B(1:M+1, 1:M+1) = T_B_block;
+
+% block nu = 2, mu = 2
+denom = -1i * (beta_n * one_vec.' - one_vec * beta_n');
+exp_plus = exp(denom * (xi1 + a));
+exp_minus = exp(denom * (xi1 - a));
+T_B_block = (exp_plus - exp_minus) ./ denom .* pre_factors;
+T_B_block(diag_indeces) = 0.5 * a * k_out^2 ./ beta_n(1:max_real_index).^2 .* gamma_n_m(diag_indeces);
+T_B(M+2:end, M+2:end) = T_B_block;
+
+% block nu = 1, mu = 2
+denom = 1i * (beta_n * one_vec.' + one_vec * beta_n');
+exp_plus = exp(denom * (xi1 + a));
+exp_minus = exp(denom * (xi1 - a));
+T_B_block = (exp_plus - exp_minus) ./ denom .* pre_factors;
+diag_indeces = ( (M+2)*max_real_index + (1:M+2:(M+2)*(M + 1 - max_real_index)) ).';
+T_B_block(diag_indeces) = 0.5 * a * k_out^2 ./ abs(beta_n(max_real_index+1:end)).^2 .* gamma_n_m(diag_indeces) ...
+ .* exp(-2 * R * abs(beta_n(max_real_index+1:end)));
+T_B(1:M+1, M+2:end) = T_B_block;
+
+% block nu = 2, mu = 1
+denom = -1i * (beta_n * one_vec.' + one_vec * beta_n');
+exp_plus = exp(denom * (xi1 + a));
+exp_minus = exp(denom * (xi1 - a));
+T_B_block = (exp_plus - exp_minus) ./ denom .* pre_factors;
+T_B_block(diag_indeces) = 0.5 * a * k_out^2 ./ abs(beta_n(max_real_index+1:end)).^2 .* gamma_n_m(diag_indeces) ...
+ .* exp(-2 * R * abs(beta_n(max_real_index+1:end)));
+T_B(M+2:end, 1:M+1) = T_B_block;
+
+end
+
+
+% The matrix (gamma_n^m)_{nm} contains the scalar products of two Neumann
+% functions over the interval (\xi_2-a, \xi_2 + a) \subseteq (0, H). It is
+% symmetric.
+function g_n_m = gamma_n_m_eval(M, k_n, xi2, a, H)
+
+g_n_m = zeros(M+1, M+1);
+k_n_nat = k_n(2:end);
+
+% entry for n = m = 0
+g_n_m(1, 1) = 2*a / H;
+
+% other entries for first row and first column
+v = sqrt(2) ./ (H * k_n_nat) .* (sin(k_n_nat*(xi2+a)) - sin(k_n_nat*(xi2-a)));
+g_n_m(1, 2:end) = v;
+g_n_m(2:end, 1) = v.';
+
+% assign the diagonal of the matrix except for first entry
+g_n_m(M+3:M+2:end) = 2*a / H + (sin(2*k_n_nat*(xi2+a)) - sin(2*k_n_nat*(xi2-a))) ./ (2 * k_n_nat * H);
+
+% assign all the sub-diagonals
+for l=1:M
+ k_diff = k_n(2+l:end) - k_n(2:end-l);
+ k_sum = k_n(2+l:end) + k_n(2:end-l);
+ v = sin(k_diff*(xi2+a)) / (k_diff * H) + sin(k_sum*(xi2+a)) / (k_sum * H) ...
+ - sin(k_diff*(xi2-a)) / (k_diff * H) - sin(k_sum*(xi2-a)) / (k_sum * H);
+ g_n_m((M+1)*(l+1)+2:M+2:(M+1)^2-l) = v;
+ g_n_m((M+2)+l+1:M+2:(M+1-l)*(M+1)) = v;
+end
+
+end
diff --git a/src/plot_of_figure_1.m b/src/plot_of_figure_1.m
new file mode 100644
index 0000000..378ff23
--- /dev/null
+++ b/src/plot_of_figure_1.m
@@ -0,0 +1,52 @@
+function [] = plot_of_figure_1()
+
+H = 5;
+R = 6;
+x_max = 5;
+
+curve = BIEPack.geom.Ellipse([0; 0], 0.6, 0.2);
+curve.rotate(pi/6);
+curve.translate([-0.4; 2.6]);
+
+k_out = 6;
+k_in = sqrt(2) * k_out;
+M = 10;
+
+N_q = operator_N(H, curve, k_out, k_in, R, M, 512, 80, true);
+
+% Define the plots
+a = 0.01;
+%alpha = [0.1 0.2 0.5 0.75];
+alpha = [0.05 0.06 0.1 0.2];
+fig_nr = 1:length(alpha);
+indicat = cell(1, length(fig_nr));
+%step = 0.1;
+step = 0.05;
+s = linspace(-pi, pi, 201);
+obs_height = 25;
+
+% lambda = eig(N_q);
+% figure(length(alpha)+1);
+% semilogy(abs(lambda));
+
+for j = 1:length(fig_nr)
+ [indicat{j}, Xi1, Xi2] = indicator_func(N_q, alpha(j), k_out, H, R, x_max, step, a, true);
+
+ figure(fig_nr(j))
+ colormap(parula);
+ surf(Xi1, Xi2, indicat{j});
+ hold on
+ [x, y] = curve.eval(s, 0);
+ plot3(x, y, obs_height*ones(size(x)), 'w');
+ hold off
+ view(2);
+ axis equal
+ axis([-5 5 0 H]);
+ shading interp
+ colorbar
+ set(gcf, 'Color', 'w');
+ export_fig(sprintf('../figures/ellipse_central_k_6.0_alpha_%04.2f.png', alpha(j)));
+end
+
+end
+
diff --git a/src/plot_of_figure_2.m b/src/plot_of_figure_2.m
new file mode 100644
index 0000000..1caf0be
--- /dev/null
+++ b/src/plot_of_figure_2.m
@@ -0,0 +1,49 @@
+function [] = plot_of_figure_2()
+
+H = 5;
+R = 6;
+
+curve = BIEPack.geom.Circle([1.5; 3.5], 0.35);
+
+k_out = 11.0;
+k_in = sqrt(2)*k_out;
+M = 17;
+
+N_q = operator_N(H, curve, k_out, k_in, R, M, 512, 80, true);
+%lambda = eig(N_q);
+
+% Define the plots
+a = 0.01;
+alpha = [0.01 0.025 0.0375 0.05];
+fig_nr = 1:length(alpha);
+indicat = cell(1, length(fig_nr));
+x_max = 5;
+%step = 0.1;
+step = 0.05;
+s = linspace(-pi, pi, 201);
+obs_height = 25;
+
+% figure(length(alpha)+1);
+% semilogy(abs(lambda));
+
+for j = 1:length(fig_nr)
+ [indicat{j}, Xi1, Xi2] = indicator_func(N_q, alpha(j), k_out, H, R, x_max, step, a, true);
+
+ figure(fig_nr(j))
+ colormap(parula);
+ surf(Xi1, Xi2, indicat{j});
+ hold on
+ [x, y] = curve.eval(s, 0);
+ plot3(x, y, obs_height*ones(size(x)), 'w');
+ hold off
+ view(2);
+ axis equal
+ axis([-5 5 0 H]);
+ shading interp
+ colorbar
+ set(gcf, 'Color', 'w');
+ export_fig(sprintf('../figures/circle_off_center_k_%4.2f_alpha_%4.2f.png', k_out, alpha(j)));
+end
+
+end
+
diff --git a/src/plot_of_figure_3.m b/src/plot_of_figure_3.m
new file mode 100644
index 0000000..ec9a9e6
--- /dev/null
+++ b/src/plot_of_figure_3.m
@@ -0,0 +1,52 @@
+function [] = plot_of_figure_3()
+
+H = 5;
+R = 6;
+
+curve = BIEPack.geom.Ellipse([0; 0], 0.5, 0.2);
+curve.rotate(pi/6);
+curve.translate([-1; 2.5]);
+
+k_out = 6.0;
+M = 9;
+
+%alpha = [-0.01 -0.008 -0.007 -0.006 ];
+alpha = [-0.01 -0.02 -0.03 -0.04];
+q = -0.5 * ones(size(alpha));
+fig_nr = 1:length(q);
+a = 0.01;
+
+indicat = cell(1, length(fig_nr));
+x_max = 5;
+step = 0.05;
+s = linspace(-pi, pi, 201);
+obs_height = 25;
+
+for j = 1:length(fig_nr)
+
+ N_q = operator_N(H, curve, k_out, sqrt((1+ q(j)))*k_out, R, M, 512, true);
+ [indicat{j}, Xi1, Xi2] = indicator_func(N_q, alpha(j), k_out, H, R, x_max, step, a, false);
+
+ figure(fig_nr(j))
+ colormap(parula);
+ surf(Xi1, Xi2, indicat{j});
+ hold on
+ [x, y] = curve.eval(s, 0);
+ plot3(x, y, obs_height*ones(size(x)), 'w');
+ hold off
+ view(2);
+ axis equal
+ axis([-5 5 0 5]);
+ shading interp
+ colorbar
+ set(gcf, 'Color', 'w');
+ export_fig(sprintf('../figures/ellipse_neg_q_alpha_%05.3f.png', -alpha(j)))
+end
+
+% lambda = eig(N_q);
+% figure(length(alpha)+1);
+% semilogy(abs(lambda));
+
+
+end
+
--
GitLab