-
Notifications
You must be signed in to change notification settings - Fork 0
/
example3_bounded_inverse.m
99 lines (96 loc) · 2.79 KB
/
example3_bounded_inverse.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
% Copyright Mohamed Nasser 2017
% Please cite this MATLAB functions as:
%
% When citing this software please mention the URL of the master repository
% (https://github.com/mmsnasser/CircularMap), and the paper
% M.M.S. Nasser,Fast Computation of the Circular Map, Computational Methods
% and Function Theory, 15 (2015) 187-223.
%
%
% PLEASE note that this toolbox contains the files:
% zfmm2dpart.m
% fmm2d_r2012a.mexw32
% fmm2d_r2012a.mexw64
% pthreadGC2-w32.dll
% pthreadGC2-w64.dll
% From the Toolbox:
% L. G REENGARD AND Z. G IMBUTAS , FMMLIB2D: A MATLAB toolbox for
% fast multipole method in two dimensions, Version 1.2, 2012.
% http://www.cims.nyu.edu/cmcl/fmm2dlib/fmm2dlib.html
% PLEASE also cite the FMMLIB2D toolbox.
clear
clc
%%
n = 2^10
ratio = 0.2;
t = (0:2*pi/n:2*pi-2*pi/n).';
%%
Mat_coef = [0 4 2
1.0 1.8 -0.08
-1.0 0.08 -1.6
];
cent = Mat_coef(:,1);
radx = Mat_coef(:,2);
rady = Mat_coef(:,3);
m = length(cent)-1
%%
alphain = 0;
for k=1:m+1
et(1+(k-1)*n:k*n,1) = cent(k)+0.5.*(+radx(k).*cos(t)+i*rady(k).*sin(t));
etp(1+(k-1)*n:k*n,1) = 0.5.*(-radx(k).*sin(t)+i*rady(k).*cos(t));
end
%%
%%
[zet,zetp,cntd,rad]=circmap(et,etp,alphain,n);
%%
%%
[xeh,yeh]=meshgrid(-1:0.001:1,[-0.8,-0.6,-0.4,-0.37,-0.34,-0.2,0,0.2,0.34,0.37,0.4,0.6,0.8]);
weh=xeh+i.*yeh;
weho = weh-cntd(1);
weh(abs(weho)>=rad(1))=NaN+i*NaN;
weho = weh-cntd(2);
weh(abs(weho)<=rad(2))=NaN+i*NaN;
weho = weh-cntd(3);
weh(abs(weho)<=rad(3))=NaN+i*NaN;
wehv = weh(abs(weh)>=0).';
%%
[xev,yev]=meshgrid([-0.995,-0.98,-0.965,-0.8,-0.62,-0.6,-0.57,-0.4,-0.2,0,0.14,0.17,0.2,0.4,0.6,0.8,0.85,0.88],-1:0.001:1);
wev=xev+i.*yev;
wevo = wev-cntd(1);
wev(abs(wevo)>=rad(1))=NaN+i*NaN;
wevo = wev-cntd(2);
wev(abs(wevo)<=rad(2))=NaN+i*NaN;
wevo = wev-cntd(3);
wev(abs(wevo)<=rad(3))=NaN+i*NaN;
wevv = wev(abs(wev)>=0).';
%%
figure;
hold on
box on
plot(real(wehv),imag(wehv),'or','MarkerSize',1)
plot(real(wevv),imag(wevv),'ob','MarkerSize',1)
for k=1:m+1
c_cr = zet((k-1)*n+1:k*n,1); c_cr(n+1) = c_cr(1);
plot(real(c_cr),imag(c_cr),'k','LineWidth',1.5)
end
plot(0,0,'ok')
axis equal
axis([-1.05 1.05 -1.05 1.05])
%%
%%
zehv = fcau (zet,zetp,et,wehv);
zevv = fcau (zet,zetp,et,wevv);
%%
figure;
hold on
box on
plot(real(zehv),imag(zehv),'or','MarkerSize',1)
plot(real(zevv),imag(zevv),'ob','MarkerSize',1)
for k=1:m+1
c_cr = et((k-1)*n+1:k*n,1); c_cr(n+1) = c_cr(1);
plot(real(c_cr),imag(c_cr),'k','LineWidth',1.5)
end
plot(real(alphain),imag(alphain),'ok')
axis equal
axis([-2.05 2.05 -1.05 1.05])
%%