tp@0: function [rs,thetas,zs,rr,thetar,zr] = EDB1coordtrans2(xsou,xrec,xwedge,nvec1)
tp@0: % EDB1coordtrans2 - Transforms two sets 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: %	xrec		Matrix, [n2,3] of cartesian coordinates of n2 other 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: %	rr, thetar, zr		cyl. coord. of the points in xrec
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) 20050922
tp@0: %
tp@0: % [rs,thetas,zs,rr,thetar,zr] = EDB1coordtrans2(xsou,xrec,xwedge,nvec1)
tp@0: 
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: 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: [npoints,slask] = size(xrec);
tp@0: if npoints >0
tp@0: 	xrec =    (A*( xrec.' - xneworigo(ones(npoints,1),:).' )).';
tp@0: 	rr = sqrt( sum(xrec(:,1:2).'.^2) ).';
tp@0: 	zr = xrec(:,3);
tp@0: 	thetar = zeros(npoints,1);
tp@0: 	iv = find(rr>0);
tp@0: 	if ~isempty(iv),	
tp@0: 		thetar(iv) = real( acos( xrec(iv,1)./rr(iv) ).*( xrec(iv,2) ~= 0) );
tp@0: 		thetar(iv) = thetar(iv) + pi*( (xrec(iv,2)==0) & xrec(iv,1) < 0 );
tp@0: 		thetar(iv) = thetar(iv).*( xrec(iv,2) >=0 ) + (2*pi - thetar(iv)).*( xrec(iv,2) < 0 );
tp@0: 	end
tp@0: else
tp@0: 	rr = []; thetar =[]; zr = [];
tp@0: end