角度逼近圆弧插补法MATLAB与C语言实现插补与仿真

插补算法

插补（Interpolation），即机床数控系统依照一定方法确定刀具运动轨迹的过程。也可以说，已知曲线上的某些数据，按照某种算法计算已知点之间的中间点的方法，也称为“数据点的密化”；数控装置根据输入的零件程序的信息，将程序段所描述的曲线的起点、终点之间的空间进行数据密化，从而形成要求的轮廓轨迹，这种“数据密化”机能就称为“插补”。

在这里插入图片描述
在这里插入图片描述

圆弧插补的基本概念

如图2所示，圆弧P 0​ P e 是圆心在坐标原点的逆圆弧段，P 0 (x 0 ,y 0​ ) 是圆弧起点，P e (x e​ ,y e​ )是圆弧终点，R是圆弧半径。F是刀具沿圆弧段的移动速度。圆弧插补的要求是根据编程指定的刀具移动速度F在圆弧P 0 P e​ 上找出若干个插补点P 1​ ,P 2​ ,P 3​ ,…,P i​ ,…用若干个微小直线段来代替圆弧，使刀具在每个插补周期T内沿小直线段从一个插补点移动到另一个插补点。
起始角：

终止角：

每个插补周期半径矢量移动角:

每个插补周期半径矢量与x轴正向夹角:

则每个插补点的计算为:

角度逼近圆弧插补算法

    为了求⁡cos⁡θ和sin⁡θ，可以先求

    根据角度逼近定理，将角度θ近似分解成代数和的形式

    算法插补点的核心求法:

    KaTeX parse error: Can't use function '$' in math mode at position 3: $̲\left.x_{n}^{2}…

    \begin{array}{c}
    r=\sqrt{\left(x_{0}{2}+y_{0}{2}\right) K_{n}} \
    x=R \cos \theta \approx R \cos \theta_{n}=R k x_{n}, y=R \sin \theta \approx R \sin \theta_{n}=R k y_{n}
    \end{array}
    $$

    在这里插入图片描述

    插补算法的比较

    除了角度逼近圆弧插补算法外，通常使用的还有EDDA算法（数值积分法）和直接函数法。

    在这里插入图片描述

    表中提到角度逼近算法的速度误差易控制是指该误差有明确的误差公式，且可以通过限制与输入参考量无关的迭代次数n来控制。表中没有直接比较三种算法的运算速度，这是因为三种算法所用的计算方法不同，所以在不同的硬件支持条件下，它们的运算速度是不同的。例如在有快速算术协处理器的系统中，EDDA或直接函数法运算速度要快一些，而在没有硬件乘除部件的系统中，则角度逼近法运算速度要快一些。

    在象限、整圆、终点判别处理中，角度逼近法较其它两种算法要方便一些。这是因为角度逼近法是以角度为插补参量，而其它两种算法是以坐标(x，y)为插补参量的。

    角度逼近圆弧算法

    程序的局限性及展望

    目前该程序存在一定的局限性，它的插补圆心目前是固定在原点，并不能随着插补起点终点和半径自动算出，还有改程序需要输入角度，也与实际数控系统的输入有所不同，算法内的迭代次数也不能根据实际的弧长自动选择合适的迭代次数。目前的改进思路是根据现有的起点坐标，输入的终点坐标和圆弧半径，代入圆的函数解析式，求出圆心坐标(a,b)，再进行坐标变换，计算插补圆弧的角度；同时可以建立一个圆弧角度与迭代次数的关系，让迭代次数伴随着圆弧角度自动适配。

    参考文献

    [1] 除书法,周建来,唐学飞. 基于角度逼近法的圆弧插补算法[J]. 现代制造工程,2002,(08):13-14.
    [2] 张力. 一种新的圆弧插补算法——角度逼近法[J]. 北京航空航天大学学报,1991,(01):77-84.

    /***********************************************************************
角度逼近圆弧数控插补算法
   Editor：Ezekiel Mok/JyLi/Tiger chen
   Function:          void AXY_interpretion(double xs, double ys, double xe, 
                                            double ye, double x0, double y0,double ANG)
   Description:       在X-Y轴所构成的平面上，利用角度逼近原理计算圆弧的插补值系列
   Parameters:        xs, ys     起点 
                     xe ye   终点
                    x0,y0    圆心坐标
						ANG	插补圆弧角度
   Return Value:      大于或等余0      给予此运动命令的编码
                      小于0            失败，传回值的意义可参考错误信息代码
***************************************************************************/
#include "stdio.h"
#include "math.h"
#include "conio.h"
#define n 1000//算法迭代次数
#define F 0.25//进给速度
## 参考文献：

#define T 0.5// 插补周期
//角度均以弧度制表示
/*************************************************************************
**************************************************************************/
void AXY_interpretion(double xs, double ys, double xe, double ye, double x0, double y0, double ANG);

int main()
{ double xs;
/**********************
	xs=atan(1/(2^0));
	xs=atan(1/pow(2,0));
	xs=1;
	ys=1;
	xs=atan(ys/xs);
	xs=atan(1/(2^0));
  	double k=2,Kn=1;
	xs=(Kn*(1+1/pow(2,2*k)));
  **********************/
	xs=pow(2,8);
	printf("int is:%f\n",xs);
	
    AXY_interpretion(1.5,0,-1.5,0,0,0,180);//xs, ys(起点)xe ye(终点) x0,y0(圆心坐标) ANG(角度)
    return 0;
}
void AXY_interpretion(double xs, double ys, double xe, double ye, double x0, double y0, double ANG)
{   
     FILE *f1;
	double i=0,j=0,flag=1;
	double R; 
	double theta_i,theta_S,theta_E,deta_theta,R0,S_0,S_i,deta_x,deta_y,A_Judge,A_Judge1,x_calulator,y_calulator;//角度逼近圆弧算法参数
	double x_temp,y_temp,S_j;//存放i-1的数据 
	double x_inter,y_inter,x_i,y_i,deta_xi,deta_yi;//插补操作参数
	double Kn=1,k=1;//计算逼近常数
/*******************************************************************
初始化参数
********************************************************************/
    R=(sqrt(pow(xs,2)+pow(ys,2)));
	S_0=atan(1/pow(2,0));
	theta_S=atan(ys/xs);
	printf("theta_S:%f\n",theta_S);
    ANG=3.14159*ANG/180;
	theta_E=ANG;
//	theta_E=atan(ye/xe);
	theta_i=theta_S;
	deta_theta=(F*T)/R;
	printf("deta_theta:%f\n",deta_theta);
   // f1=fopen("D:\My Documents\桌面\INTERP\c.txt","w+");
     f1=fopen("C:\\c.txt","w+");
	for(k=1;k<=n;k++)
		{
			Kn=(Kn*(1+1/pow(2,2*k)));
		}
			printf("Kn:%f\n",Kn);
 	R0=R/(sqrt(2*Kn));
	S_0=atan((1/pow(2,0)));
	printf("i:%f\n",theta_i);
	printf("E:%f\n",theta_E);	
	   i=1;
	while(theta_i<=theta_E&&flag==1)//避免死循环 
		{
			theta_i=theta_i+deta_theta;
			if(theta_i>theta_E)
	    	  {
	    	   theta_i=theta_E;
	    	   flag=0;
			  }
			printf("thi:%f\n",theta_i);
			if(theta_i<=1.5708&&theta_i>0)
			{
				    S_0=atan((1/pow(2,0)));
					A_Judge=theta_i-S_0; 
					x_inter=R0;
					y_inter=R0;
					j=1;
					while(j<=n)
					{ 
					 S_j=atan((1/pow(2,j)));
					 x_calulator=(1/pow(2,j))*y_inter;
					 y_calulator=(1/pow(2,j))*x_inter;
					 if(A_Judge>=0)
					 {
					 	y_inter=y_inter+y_calulator;
					 	x_inter=x_inter-x_calulator;
					 	A_Judge=A_Judge-S_j;
					 }
					else
					 {
					 	y_inter=y_inter-y_calulator;
					 	x_inter=x_inter+x_calulator;
					 	A_Judge=A_Judge+S_j;
					 }
					 j++;	
					}
					x_i=x_inter;
					y_i=y_inter;
					deta_xi=x_i-x_temp;
					x_temp=x_i;
		            deta_yi=y_i-y_temp;
					y_temp=y_i;
				  fprintf(f1,"%f,",x_i);
                  fprintf(f1,"%f\n",y_i);
				}
			if(theta_i>1.5708&&theta_i<=3.14159)
			  {
			  		S_0=atan((1/pow(2,0)));
					A_Judge=theta_i-S_0-1.5708; 
					x_inter=-R0;
					y_inter=R0;
					j=1;
					while(j<=n)
					{ 
					 S_j=atan((1/pow(2,j)));
					 x_calulator=(1/pow(2,j))*y_inter;
					 y_calulator=(1/pow(2,j))*x_inter;
					 if(A_Judge>=0)
					 {
					 	y_inter=y_inter+y_calulator;
					 	x_inter=x_inter-x_calulator;
					 	A_Judge=A_Judge-S_j;
					 }
					else
					 {
					 	y_inter=y_inter-y_calulator;
					 	x_inter=x_inter+x_calulator;
					 	A_Judge=A_Judge+S_j;
					 }
					 j++;	
					}
					x_i=-x_inter;
					y_i=-y_inter;
					deta_xi=x_i-x_temp;
					x_temp=x_i;
		            deta_yi=y_i-y_temp;
					y_temp=y_i;
				fprintf(f1,"%f,",-x_i);
                fprintf(f1,"%f\n",-y_i);
			  	
			  }     
			if(theta_i>3.14159&&theta_i<=4.7124)
			  {
			  		S_0=atan((1/pow(2,0)));
					A_Judge=theta_i-S_0-3.14159; 
					x_inter=-R0;
					y_inter=-R0;
					j=1;
					while(j<=n)
					{ 
					 S_j=atan((1/pow(2,j)));
					 x_calulator=(1/pow(2,j))*y_inter;
					 y_calulator=(1/pow(2,j))*x_inter;
					 if(A_Judge>=0)
					 {
					 	y_inter=y_inter+y_calulator;
					 	x_inter=x_inter-x_calulator;
					 	A_Judge=A_Judge-S_j;
					 }
					else
					 {
					 	y_inter=y_inter-y_calulator;
					 	x_inter=x_inter+x_calulator;
					 	A_Judge=A_Judge+S_j;
					 }
					 j++;	
					}
					x_i=x_inter;
					y_i=y_inter;
					deta_xi=x_i-x_temp;
					x_temp=x_i;
		            deta_yi=y_i-y_temp;
					y_temp=y_i;
				fprintf(f1,"%f,",x_i);
                fprintf(f1,"%f\n",y_i);
			  	
			  }                    
           	if(theta_i>4.7124&&theta_i<=6.283)
			  {
			  		S_0=atan((1/pow(2,0)));
					A_Judge=theta_i-S_0-4.71; 
					x_inter=R0;
					y_inter=-R0;
					j=1;
					while(j<=n)
					{ 
					 S_j=atan((1/pow(2,j)));
					 x_calulator=(1/pow(2,j))*y_inter;
					 y_calulator=(1/pow(2,j))*x_inter;
					 if(A_Judge>=0)
					 {
					 	y_inter=y_inter+y_calulator;
					 	x_inter=x_inter-x_calulator;
					 	A_Judge=A_Judge-S_j;
					 }
					else
					 {
					 	y_inter=y_inter-y_calulator;
					 	x_inter=x_inter+x_calulator;
					 	A_Judge=A_Judge+S_j;
					 }
					 j++;	
					}
					x_i=-x_inter;
					y_i=-y_inter;
					deta_xi=x_i-x_temp;
					x_temp=x_i;
		            deta_yi=y_i-y_temp;
					y_temp=y_i;
				fprintf(f1,"%f,",-x_i);
                fprintf(f1,"%f\n",-y_i);
			  	
			  }               
		   i+=1;
	
		}
	
}  


    clc,clear all;
plot_data=load('c:\\c.txt');
%  plot_data=load('D:\My Documents\Desktop\INTERP\My Documents桌面INTERPc.txt');
x=0:0.01:1;
y=sqrt(1-x.^2);
z=x.*0;
plot(x,y,'r');
hold on;
plot(-x,y,'r');
plot(-x,-y,'r');
plot(x,-y,'r');
grid on;
for i=1:1:size(plot_data)-1
   A1=plot_data(i,1);
   B1=plot_data(i,2);
   A2=plot_data(i+1,1);
   B2=plot_data(i+1,2);
   hold on
   if(rem(i,2)==0)
   line([A1,A2],[B1,B1]);
   pause(0.3);
   line([A2,A2],[B1,B2]);
   pause(0.3);
   else
   line([A1,A1],[B1,B2]);
   pause(0.3);
   line([A1,A2],[B2,B2]);
   pause(0.3);
   end
end
hold on;
plot_data=load('c:\\c.txt');
plot(plot_data(:,1),plot_data(:,2),'*')