Home > atmlab > geodetic > circle_intersect.m

circle_intersect

PURPOSE ^

CIRCLE_INTERSECT returns intersection points of two circles in two dimensions

SYNOPSIS ^

function [y]=circle_intersect(P1,P2,R1,R2);

DESCRIPTION ^

 CIRCLE_INTERSECT returns intersection points of two circles in two dimensions

 Format [y]=circle_intersect(P1,P2,r1,r2);
      
      OUT  y coordinates of the intersection points
              y(1,:) are closest to origin (x=0,y=0) 
              if y=nan there is no intersection
     
     IN    P1 coordinates(x,y) of the circle 1 center
             P2 coordinates(x,y) of the circle 2 center
             r1 radius of circle 1
             r2 radius of circle 2

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

circle_intersect.m

SOURCE CODE ^

0001 % CIRCLE_INTERSECT returns intersection points of two circles in two dimensions
0002 %
0003 % Format [y]=circle_intersect(P1,P2,r1,r2);
0004 %
0005 %      OUT  y coordinates of the intersection points
0006 %              y(1,:) are closest to origin (x=0,y=0)
0007 %              if y=nan there is no intersection
0008 %
0009 %     IN    P1 coordinates(x,y) of the circle 1 center
0010 %             P2 coordinates(x,y) of the circle 2 center
0011 %             r1 radius of circle 1
0012 %             r2 radius of circle 2
0013 %
0014 
0015 % Histiry: created by Bengt Rydberg 2011-11-15
0016 function [y]=circle_intersect(P1,P2,R1,R2);
0017 
0018 Pv=P1-P2;
0019 d=sqrt(Pv*Pv');
0020 
0021 if d>(R1+R2)
0022  %no solutions (the circles are separate.)
0023  y=nan;
0024 elseif d<abs(R1-R2)
0025  %no solutions (one circle is contained within the other.)
0026  y=nan;
0027 elseif d==0 & R1==R2
0028  %infinite number of solutions (circles are coincident)
0029  y=inf;
0030 else
0031 
0032  a=(R1^2-R2^2+d^2)/(2*d);
0033  h=sqrt(R1^2-a^2);
0034  P3=P1+a*(P2-P1)/d;
0035  x4a=P3(1)+h*(P2(2)-P1(2))/d;
0036  x4b=P3(1)-h*(P2(2)-P1(2))/d;
0037  y4a=P3(2)-h*(P2(1)-P1(1))/d;
0038  y4b=P3(2)+h*(P2(1)-P1(1))/d;
0039 
0040  P4a=[x4a,y4a];
0041  P4b=[x4b,y4b];
0042  
0043  d4a=sqrt(P4a*P4a');
0044  d4b=sqrt(P4b*P4b');
0045   
0046  if d4a==d4b
0047   %only one intersection
0048   y=[P4a;P4b];
0049  elseif d4a>d4b
0050   y=[P4b;P4a]; 
0051  else
0052   y=[P4a;P4b]; 
0053  end
0054 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005