1+ % Evaluate partial derivatives of wavefunction expansion series.
2+ %
3+ % [dr,dtheta,er,etheta,rad] = dsumcof(x,x0,k,c,'H') returns the [partial
4+ % derivatives dr and dtheta of the radiating wavefunction expansion with
5+ % coefficients c, centre x0 and wavenumber k at points x. The unit
6+ % vectors er and etheta associated with dr and dtheta are also computed.
7+ %
8+ % [...] = dsumcof(x,x0,k,c,'J') returns the values z of the regular
9+ % wavefunction expansion with coefficients c, centre x0 and wavenumber k
10+ % at points x.
11+ %
12+ % Note: in the above vectors in the plane are represented by
13+ % complex numbers.
14+ %
15+ % Stuart C. Hawkins - 7 May 2024
16+
17+ % Copyright 2014, 2015, 2016, 2017, 2018, 2022, 2023, 2024 Stuart C. Hawkins and M. Ganesh.
18+ %
19+ % This file is part of TMATROM.
20+ %
21+ % TMATROM is free software: you can redistribute it and/or modify
22+ % it under the terms of the GNU General Public License as published by
23+ % the Free Software Foundation, either version 3 of the License, or
24+ % (at your option) any later version.
25+ %
26+ % TMATROM is distributed in the hope that it will be useful,
27+ % but WITHOUT ANY WARRANTY; without even the implied warranty of
28+ % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29+ % GNU General Public License for more details.
30+ %
31+ % You should have received a copy of the GNU General Public License
32+ % along with TMATROM. If not, see <http://www.gnu.org/licenses/>.
33+
34+
35+ function [dr ,dtheta ,er ,etheta ,rad ] = dsumcof(points ,centre ,kwave ,cof ,type )
36+
37+ if strcmp(type ,' F' )
38+ error(' Type F not supported' )
39+ end
40+
41+ % -------------------------------------------------
42+ % setup
43+ % -------------------------------------------------
44+
45+ % make sure the coefficient vector is a column vector
46+ cof= cof(: );
47+
48+ % determine the maximum order from the length of the coefficient vector
49+ nmax= 0.5 *(length(cof )-1 );
50+
51+ % create a vector of indexes.... helps to vectorize the computation
52+ n= -nmax : nmax ;
53+
54+ % -------------------------------------------------
55+ % turn points into a vector
56+ % -------------------------------------------------
57+
58+ % get the shape of points so we can restore it later
59+ [np ,mp ]=size(points );
60+
61+ % reshape into a vector
62+ p= reshape(points - centre ,np * mp ,1 );
63+
64+ % -------------------------------------------------
65+ % compute the field
66+ % -------------------------------------------------
67+
68+ % convert to polar coordinates
69+ theta= angle(p );
70+ rad= abs(p );
71+
72+ % make a matrix from n and rad
73+ [nd ,rd ]=meshgrid(n ,kwave * rad );
74+
75+ % get Bessel/Hankel/far field values as appropriate
76+ if strcmp(type ,' J' )
77+
78+ bess= besselj(abs(nd ),rd );
79+ bessd= kwave * besseljd(abs(nd ),rd );
80+
81+ elseif strcmp(type ,' H' )
82+
83+ bess= besselh(abs(nd ),rd );
84+ bessd= kwave * besselhd(abs(nd ),rd );
85+
86+ end
87+
88+ % compute the angular part
89+ Y= exp(1i * theta * n );
90+ Yd= Y * diag(1i * n );
91+
92+ % put it together
93+ M = bess .* Y ;
94+ Mdr = bessd .* Y ;
95+ Mdtheta = bess .* Yd ;
96+
97+ % -------------------------------------------------
98+ % make the return value the same shape as the original
99+ % array of points
100+ % -------------------------------------------------
101+
102+ % compute the sum of the wavefunctions and reshape
103+ dr= reshape(Mdr * cof ,np ,mp );
104+ dtheta= reshape(Mdtheta * cof ,np ,mp );
105+
106+ if nargout > 2
107+ er = reshape(p ./ abs(p ),np ,mp );
108+ etheta = reshape(exp(1i * pi / 2 )*er ,np ,mp );
109+ rad = reshape(rad ,np ,mp );
110+ end
0 commit comments