tp@0: function [rs,thetas,zs] = EDB1coordtrans1(xsou,xwedge,nvec1)
tp@0: % EDB1coordtrans1 - Transforms one set of cartesian coordinates to edge-related cylindrical coordinates.
tp@0: % The cyl. coord. system is defined so that:
tp@0: %     A z-axis is placed along the edge, from the given endpoint 1 to the given
tp@0: %     endpoint 2.
tp@0: %     The origo of the cyl. syst. will be edge endpoint 1.
tp@0: %     The theta-angles of the cyl. coord. syst. will refer to the
tp@0: %     reference plane of the edge.
tp@0: % The ref. plane of the edge is described by its normal vector.
tp@0: % NB! The order of the edge points is important!!
tp@0: % The vector going from xwedge(1,:) to xwedge(2,:) must be
tp@0: % oriented so that if the RH thumb is along this vector, the tips
tp@0: % of the fingers "come out of" the open face of plane1, i.e. where nvec1
tp@0: % is the normal vector.
tp@0: %
tp@0: % Input parameters:
tp@0: %	xsou		Matrix, [n1,3] of cartesian coordinates of n1 points. 
tp@0: %	xwedge 		Matrix, [2,3], with the cartesian coordinates of the two 
tp@0: %				wedge end points: [xw1 yw1 zw1;xw2 yw2 zw2].
tp@0: %   nvec1       List, [1,3], with the normal vector of the reference plane
tp@0: %               of the edge.
tp@0: %
tp@0: % Output parameters:
tp@0: %	rs, thetas, zs		cyl. coord. of the points in xsou
tp@0: %
tp@0: % Uses the function EDB1cross
tp@0: %
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: %   This file is part of the Edge Diffraction Toolbox by Peter Svensson.                       
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify       
tp@0: %   it under the terms of the GNU General Public License as published by the Free Software     
tp@0: %   Foundation, either version 3 of the License, or (at your option) any later version.        
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is distributed in the hope that it will be useful,       
tp@0: %   but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS  
tp@0: %   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.             
tp@0: %                                                                                              
tp@0: %   You should have received a copy of the GNU General Public License along with the           
tp@0: %   Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.                 
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20061118
tp@0: %
tp@0: % [rs,thetas,zs] = EDB1coordtrans1(xsou,xwedge,nvec1)
tp@0: 
tp@0: npoints = size(xsou,1);
tp@0: if npoints == 1
tp@0:     xneworigo = xwedge(1,:);
tp@0: 
tp@0:     xknown1 = xwedge(2,:) - xneworigo;
tp@0:     xknown1 = xknown1 / norm(xknown1);
tp@0: 
tp@0:     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:     A = inv([xknown1.' A nvec1.']);
tp@0: 
tp@0:     xsou =    (A([2 3 1],:)*( xsou.' - xneworigo.' )).';
tp@0: 
tp@0:     rs = norm(xsou(1:2));
tp@0:     zs = xsou(:,3);
tp@0:     thetas = 0;
tp@0:     if rs > 0
tp@0:         thetas = real( acos( xsou(1)./rs ).*( xsou(2) ~= 0) );
tp@0:         thetas = thetas + pi*( (xsou(2)==0) & xsou(1) < 0 );
tp@0:         thetas = thetas.*( xsou(2) >=0 ) + (2*pi - thetas).*( xsou(2) < 0 );
tp@0:     end
tp@0:     
tp@0: else
tp@0:     xneworigo = xwedge(1,:);
tp@0: 
tp@0:     xknown1 = xwedge(2,:) - xneworigo;
tp@0:     xknown1 = xknown1 / sqrt( sum( xknown1.^2 ));
tp@0: 
tp@0: 
tp@0:     A = [0 1 0;0 0 1;1 0 0]*inv([xknown1.' EDB1cross(nvec1.',xknown1.') nvec1.']);
tp@0: 
tp@0:     [npoints,slask] = size(xsou);
tp@0:     xsou =    (A*( xsou.' - xneworigo(ones(npoints,1),:).' )).';
tp@0: 
tp@0:     rs = sqrt( sum(xsou(:,1:2).'.^2) ).';
tp@0:     zs = xsou(:,3);
tp@0:     thetas = zeros(npoints,1);
tp@0:     iv = find(rs>0);
tp@0:     if ~isempty(iv)
tp@0:         thetas(iv) = real( acos( xsou(iv,1)./rs(iv) ).*( xsou(iv,2) ~= 0) );
tp@0:         thetas(iv) = thetas(iv) + pi*( (xsou(iv,2)==0) & xsou(iv,1) < 0 );
tp@0:         thetas(iv) = thetas(iv).*( xsou(iv,2) >=0 ) + (2*pi - thetas(iv)).*( xsou(iv,2) < 0 );
tp@0:     end
tp@0: 
tp@0: end