Rotation of curves , intersection of curve and for -loop application [Floating point comparison problem]]

Question

The problem scenario is defined as follows problem statement (Click here)

The information what happens in one single iteration is provided here in this link

I am supposed to get one small arc touching the circle at both ends (curve 1 rotates circle about point A and then rotates touching circle at point B so both ends of arc touch the circle). I am getting the result when I do this manually where I find intersection by using polyxpoly() , place the arc on intersection and keep rotating it till it touches the circle. But when I try to automate the process for many arcs by using a for loop I do not get desired result.

Since rotation of arc takes place , I have given small rotation step size but I get different results . Sometimes the for-loop works normally and gives desired result and sometimes gives abnormal results for certain iteration. I am unable to figure out the the problem with my approach.

Results are as follows click here to see the result

These are additional details which shows difference between desired result and undesired result Desired Result Vs Undesired Result for theta=500 & i=2 ,3 Desired Result Vs Undesired Result for theta=1800 & i=4

It can be seen for theta = 500 & i = 2 the arc does not rotate and since it dose not rotate it gives same result for theta = 500 & i = 3 but when theta = 1800 it gives different results but fails at iteration 4 (for theta = 1800 i = 1,2,3 gives desired result but fails at i = 4). Check the result image.

• The parameter theta can be found in rotation()

The code I used is as follows but For-loop is giving abnormal results for certain iteration depending on theta

clc,clear
r = 10;                     % arc length
R = 55;                     % radius of a circle
aa = 60*pi/180;              % arc angle
ap = 0*pi/180;             % arc position angle
k=0;

t = linspace(0,pi/2);
[x,y] = pol2cart(t,R);      % circle data
t1 = linspace(0,aa)-aa/2+ap;
[x1,y1] = pol2cart(t1,r); % arc data
% Intersection of x axis (y=0) and quarter circle
x_axis_range=1:1000;
y_val_on_x_axis=zeros(1,length(x_axis_range));
[x1rc,y1rc]=polyxpoly(x_axis_range,y_val_on_x_axis,x,y);
xc_s(1)=x1rc;
yc_s(1)=y1rc;
% x1rc=55;
% y1rc=0;
for z=1:3
[x1_p,y1_p]= placement(x,y,x1,y1,x1rc,y1rc);
[x1_rotated,y1_rotated,xc,yc]=rotate(x,y,x1_p,y1_p,x1rc,y1rc);
%[x1_rotated,y1_rotated,xc,yc]=rotate(x,y,x1,y1,x1rc,y1rc);
x1rc=xc;
y1rc=yc;
end
% havent written the plot command

Placement()- The purpose of placement () is to shift the new arc to point of intersection of previous arc and the circle

function [x1_p,y1_p] = placement(x,y,x1,y1,x1rc,y1rc)

xcen=x1rc;
ycen=y1rc;

% shifting arc-lower-end to (t(1),0) or (
delx=xcen-x1(1); % Finding the x difference between arc-lower-end x-coordinate & x(1)[circle]
dely=ycen-y1(1); % Finding the y difference between arc-lower-end y-coordinate & y(1) [circle]

x1_p=x1+delx;
y1_p=y1+dely;
end

The rotation function():

 function [x1_rotated,y1_rotated,xc,yc] = rotate(x,y,x1_p,y1_p,x1rc,y1rc)
    % [x1_p,y1_p]= placement(x,y,x1,y1,x1rc,y1rc);
    theta =linspace(0,pi,500);
    i=1;
    xc=[];
    yc=[];
    xcen=x1rc;
    ycen=y1rc;
    while isempty(xc)||xc==xcen && isempty(yc)||yc==ycen
    % create a matrix of these points, which will be useful in future calculations
    v = [x1_p;y1_p];
    % choose a point which will be the center of rotation
    x_center = xcen;
    y_center = ycen;
    % create a matrix which will be used later in calculations
    center = repmat([x_center; y_center], 1, length(x1_p));
    % define a theta degree counter-clockwise rotation matrix
    R = [cos(theta(i)) -sin(theta(i)); sin(theta(i)) cos(theta(i))];
    % do the rotation...
    s = v - center;     % shift points in the plane so that the center of rotation is at the origin
    so = R*s;           % apply the rotation about the origin
    vo = so + center;   % shift again so the origin goes back to the desired center of rotation
    % this can be done in one line as:
    % vo = R*(v - center) + center
    % pick out the vectors of rotated x- and y-data
    x1_rotated = vo(1,:);
    y1_rotated = vo(2,:);
    [xc,yc] = polyxpoly(x1_rotated,y1_rotated,x,y);
    [xc1,yc1] = polyxpoly(x1_p,y1_p,x,y);
    if length(xc)==2
    xc=xc(2);
    yc=yc(2);
    end
    i=i+1;
    end
    end

I tried modifying the "if statement" as following but it seems after a certain iteration it dosen't work. The new modification in if statement are as follows

   if length(xc)==2
   ubx=xc(1); lbx=xc(2);
   uby=yc(2);lby=yc(1);
   xc=xc(2);
   yc=yc(2)
    
  ub=[ubx uby];
  lb=[lbx-4 lby];
  end

Then I was trying to return the ub and lb and got the following error message Output argument "ub" (and maybe others) not assigned during call to "rotate". after fifth iteration of for loop . The if statement worked fine for 1,2,3,4 iteation

Answer 1

It's a floating point comparison problem. Long story short, both 3*(4/3 -1) == 1 and sin(pi)==0 return false (Avoiding Common Problems with Floating-Point Arithmetic).

Change your while condition to something like this:

while (isempty(xc)||abs(xc-xcen)<=10*eps(xcen)) && (isempty(yc)||abs(yc-ycen)<=10*eps(ycen))

BTW, your assumption that when polyxpoly returns 2 point, the second point is your result, is not true. Try to find the furthest point from start point of the arc:

if length(xc)>1
    [~, i] = max((xc-xcen).^2+(yc-ycen).^2);
    xc=xc(i);
    yc=yc(i);
end