tp@0
|
1 function [rs,thetas,zs] = EDB1coordtrans1(xsou,xwedge,nvec1)
|
tp@0
|
2 % EDB1coordtrans1 - Transforms one set of cartesian coordinates to edge-related cylindrical coordinates.
|
tp@0
|
3 % The cyl. coord. system is defined so that:
|
tp@0
|
4 % A z-axis is placed along the edge, from the given endpoint 1 to the given
|
tp@0
|
5 % endpoint 2.
|
tp@0
|
6 % The origo of the cyl. syst. will be edge endpoint 1.
|
tp@0
|
7 % The theta-angles of the cyl. coord. syst. will refer to the
|
tp@0
|
8 % reference plane of the edge.
|
tp@0
|
9 % The ref. plane of the edge is described by its normal vector.
|
tp@0
|
10 % NB! The order of the edge points is important!!
|
tp@0
|
11 % The vector going from xwedge(1,:) to xwedge(2,:) must be
|
tp@0
|
12 % oriented so that if the RH thumb is along this vector, the tips
|
tp@0
|
13 % of the fingers "come out of" the open face of plane1, i.e. where nvec1
|
tp@0
|
14 % is the normal vector.
|
tp@0
|
15 %
|
tp@0
|
16 % Input parameters:
|
tp@0
|
17 % xsou Matrix, [n1,3] of cartesian coordinates of n1 points.
|
tp@0
|
18 % xwedge Matrix, [2,3], with the cartesian coordinates of the two
|
tp@0
|
19 % wedge end points: [xw1 yw1 zw1;xw2 yw2 zw2].
|
tp@0
|
20 % nvec1 List, [1,3], with the normal vector of the reference plane
|
tp@0
|
21 % of the edge.
|
tp@0
|
22 %
|
tp@0
|
23 % Output parameters:
|
tp@0
|
24 % rs, thetas, zs cyl. coord. of the points in xsou
|
tp@0
|
25 %
|
tp@0
|
26 % Uses the function EDB1cross
|
tp@0
|
27 %
|
tp@0
|
28 % ----------------------------------------------------------------------------------------------
|
tp@0
|
29 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
tp@0
|
30 %
|
tp@0
|
31 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
tp@0
|
32 % it under the terms of the GNU General Public License as published by the Free Software
|
tp@0
|
33 % Foundation, either version 3 of the License, or (at your option) any later version.
|
tp@0
|
34 %
|
tp@0
|
35 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
tp@0
|
36 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
tp@0
|
37 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
tp@0
|
38 %
|
tp@0
|
39 % You should have received a copy of the GNU General Public License along with the
|
tp@0
|
40 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
tp@0
|
41 % ----------------------------------------------------------------------------------------------
|
tp@0
|
42 % Peter Svensson (svensson@iet.ntnu.no) 20061118
|
tp@0
|
43 %
|
tp@0
|
44 % [rs,thetas,zs] = EDB1coordtrans1(xsou,xwedge,nvec1)
|
tp@0
|
45
|
tp@0
|
46 npoints = size(xsou,1);
|
tp@0
|
47 if npoints == 1
|
tp@0
|
48 xneworigo = xwedge(1,:);
|
tp@0
|
49
|
tp@0
|
50 xknown1 = xwedge(2,:) - xneworigo;
|
tp@0
|
51 xknown1 = xknown1 / norm(xknown1);
|
tp@0
|
52
|
tp@0
|
53 A = [nvec1(2)*xknown1(3)-nvec1(3)*xknown1(2) ; nvec1(3)*xknown1(1)-nvec1(1)*xknown1(3) ; nvec1(1)*xknown1(2)-nvec1(2)*xknown1(1)];
|
tp@0
|
54 A = inv([xknown1.' A nvec1.']);
|
tp@0
|
55
|
tp@0
|
56 xsou = (A([2 3 1],:)*( xsou.' - xneworigo.' )).';
|
tp@0
|
57
|
tp@0
|
58 rs = norm(xsou(1:2));
|
tp@0
|
59 zs = xsou(:,3);
|
tp@0
|
60 thetas = 0;
|
tp@0
|
61 if rs > 0
|
tp@0
|
62 thetas = real( acos( xsou(1)./rs ).*( xsou(2) ~= 0) );
|
tp@0
|
63 thetas = thetas + pi*( (xsou(2)==0) & xsou(1) < 0 );
|
tp@0
|
64 thetas = thetas.*( xsou(2) >=0 ) + (2*pi - thetas).*( xsou(2) < 0 );
|
tp@0
|
65 end
|
tp@0
|
66
|
tp@0
|
67 else
|
tp@0
|
68 xneworigo = xwedge(1,:);
|
tp@0
|
69
|
tp@0
|
70 xknown1 = xwedge(2,:) - xneworigo;
|
tp@0
|
71 xknown1 = xknown1 / sqrt( sum( xknown1.^2 ));
|
tp@0
|
72
|
tp@0
|
73
|
tp@0
|
74 A = [0 1 0;0 0 1;1 0 0]*inv([xknown1.' EDB1cross(nvec1.',xknown1.') nvec1.']);
|
tp@0
|
75
|
tp@0
|
76 [npoints,slask] = size(xsou);
|
tp@0
|
77 xsou = (A*( xsou.' - xneworigo(ones(npoints,1),:).' )).';
|
tp@0
|
78
|
tp@0
|
79 rs = sqrt( sum(xsou(:,1:2).'.^2) ).';
|
tp@0
|
80 zs = xsou(:,3);
|
tp@0
|
81 thetas = zeros(npoints,1);
|
tp@0
|
82 iv = find(rs>0);
|
tp@0
|
83 if ~isempty(iv)
|
tp@0
|
84 thetas(iv) = real( acos( xsou(iv,1)./rs(iv) ).*( xsou(iv,2) ~= 0) );
|
tp@0
|
85 thetas(iv) = thetas(iv) + pi*( (xsou(iv,2)==0) & xsou(iv,1) < 0 );
|
tp@0
|
86 thetas(iv) = thetas(iv).*( xsou(iv,2) >=0 ) + (2*pi - thetas(iv)).*( xsou(iv,2) < 0 );
|
tp@0
|
87 end
|
tp@0
|
88
|
tp@0
|
89 end |