|
"Jane " <j.l.terry@hotmail.co.uk> wrote in message <idlsbt$9ip$1@fred.mathworks.com>...
> Hi, I have written a function (probably not very efficiently) to measure the distance along a line, defined by the angle from the negative y-axis and passing through a known point.
>
> The problem I am having is coping with instance where the user selects the angle phi = [0, 90, 180, 270,...] the problem becoems ill-conditioned. I'm sure there must be an easy work around for this, but I've become so bogged down in the details I can't see it!
>
> I would appreciate any help.
>
> gridSize = [20 20];
> phi = 90;
> es = [10 10];
>
> function [virtualx, virtualy] = arrayDistances(gridSize, phi, es)
> % arrayDistances - For an (n x m) grid, find the orthogonal distance from
> % the centre of each cell to the centre plume line of the emission source.
> % Also find the distance along the centre plume line between the emission
> % source and the intersection with the orthogonal line.
> %
> % Uses the principle that for determining the gradient of the orthogonal
> % to a straight-line; gradient = -1/m, where the line has a gradient of m.
> %
> % Input Parameters:
> % gridSize - [n m] the size of the grid, it is assumed that each
> % cell in the grid is 1m^2
> % phi - The angle, in degrees, specifying the direction of
> % wind. Taken as the angle between the centre
> % plume line and the positive y-axis
> % es - [x y] the (x,y) position of the emission source
> %
> % Output Parameters:
> % virtualx - An array of size (n,m) containing the distance
> % between the emission source and the
> % intersection of the centre plume line and the
> % orthogonal line that passes through the centre
> % point of the (i,j)th cell of the grid
> % virtualy - An array of size (n,m) containing the orthogonal
> % distance between the centre of each cell in the
> % grid and the centre plume line
> %
> % Created on: 6 December 2010
> % -------------------------------------------------------------------------
>
> % % x = [0 gridSize(1)]; % x-positions of line ends for centre plume
> % % % line - only used in plotting centre plume
> % % % line
>
> phi = pi/180 * phi; % Convert angle to radians
>
> c = es(2) - es(1)/tan(phi); % Calculate constant in equation of the
> % centre plume line, knowing it passes
> % through point es and is at an angle phi
> % from the positive y-axis
>
> % Initialise array to store virtual distance information
> virtualx = nan(gridSize(2), gridSize(1));
> virtualy = nan(gridSize(2), gridSize(1));
>
> for i = 1:gridSize(1)
> xpoint = i - 0.5; % x-position to measure distance for
> for j = 1:gridSize(2)
> ypoint = j - 0.5; % y-position to measure ditance for
>
> cprimed = ypoint + xpoint*tan(phi);
> % Intercept of orthogonal to centre plume
> % line
>
> cross = [tan(phi) 1; -1/tan(phi) 1]\[cprimed; c];
> % (x,y) position of intercept between
> % centre plume line and its orthogonal
>
> virtualy(j,i) = sign(ypoint - cross(2))*...
> sqrt( (cross(1) - xpoint)^2 + (cross(2) - ypoint)^2 );
> % orthogonal distance from point of
> % interest to centre plume line
>
> virtualx(j,i) = -1*sign(cross(1) - es(1))*sign(sin(phi))*...
> sqrt( (cross(1) - es(1))^2 + (es(2) - cross(2))^2 );
> % distance along centre plume line between
> % emitter and the point the orthogonal line
> % intersects it plume line
> end
> end
if any(phi==[0 90 180 270])
error('Ill conditioned');
end
If they can have numbers outside that range, use MOD or REM to convert. Also, be careful of floating point issues. E.g. What do you want to do if they enter 90.00002?
|