recipe_sp2.cpp

曲面加工程序源码,用C＋＋编写的加工曲线曲面的程序

//---------------------------------------------------------------------------
#include <vcl\vcl.h>
#include <fstream.h>
#include <math.h>
#pragma hdrstop

#include "recipe_sp2.h"


//僌儘乕僶儖曄悢
int kx; //Jisu+1
int ky;

double *data_x;
double *data_y;
double **data_z;
int nx;
int ny;
//data_x[0..nx-1]
//data_y[0..ny-1]
//data_z[0..nx-1,0..ny-1]

double **data_z2a;

void initial(void){

        nx=9;
        ny=7;

	data_x=new double[nx];
	data_y=new double[ny];

	data_z=new double*[nx];
	for (int i=0;i<nx;i++) data_z[i]=new double[ny];

        // for recipe
	data_z2a=new double*[nx];
	for (int i=0;i<nx;i++) data_z2a[i]=new double[ny];



        double xxx[]={-65,-50,-35,-20,0,20,35,50,65};
	double yyy[]={-8,-4,-2,0,2,4,8};
	double zzz[9][7]={ { 5.60, 1.40, 0.35, 0.00, 0.35, 1.40, 5.60},
	          	 { 4.86, 1.22, 0.30, 0.00, 0.30, 1.22, 4.86},
		         { 2.92, 0.73, 0.18, 0.00, 0.18, 0.73, 2.92},
		         { 1.50, 0.37, 0.09, 0.00, 0.09, 0.37, 1.50},
		         { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00},
		         {-1.71,-0.43,-0.11, 0.00,-0.11,-0.43,-1.71},
		         {-3.28,-0.82,-0.21, 0.00,-0.21,-0.82,-3.28},
		         {-4.73,-1.18,-0.30, 0.00,-0.30,-1.18,-4.73},
		         {-5.32,-1.33,-0.33, 0.00,-0.33,-1.33,-5.32}};



        for(int i=0;i<nx;i++) data_x[i]=xxx[i];
        for(int j=0;j<ny;j++) data_y[j]=yyy[j];
        for(int i=0;i<nx;i++)
             for(int j=0;j<ny;j++)
                  data_z[i][j]=zzz[i][j]*1.e-3;


};
void final(void){

	delete data_x;
	delete data_y;
	for (int i=0;i<nx;i++) delete data_z[i];
	delete data_z;

        //for recipe
	for (int i=0;i<nx;i++) delete data_z2a[i];
	delete data_z2a;

};
void open_file(AnsiString filename){

    //open file
    ifstream fin;
    fin.open(filename.c_str());

  if(fin!=NULL){

    int nox,noy;
    double x,y,z;

    fin >>nox;
    fin >>noy;

    final();

    nx=nox;
    ny=noy;

    data_x=new double[nx];
    data_y=new double[ny];

    data_z=new double*[nx];
    for (int i=0;i<nx;i++) data_z[i]=new double[ny];

    // for recipe
    data_z2a=new double*[nx];
    for (int i=0;i<nx;i++) data_z2a[i]=new double[ny];

    //char temp[255];

    for (int i=0;i<nox;i++){

        fin>>x;
        data_x[i]=x;
    };

    for (int i=0;i<noy;i++){

        fin>>y;
        data_y[i]=y;
    };

    for (int i=0;i<nox;i++){
        for( int j=0;j<noy;j++){

             fin>>z;
             data_z[i][j]=z;
        };
    };

    //close file
    fin.close();
  };

};

//---------------------------------------------------------------------------
void splie2(double *x1a,double *x2a,double **ya,int m,int n,double **y2a)
{
 //娭悢抣偺昞ya[0..m-1][0..n-1]偲撈棫曄悢偺昞x1a[0..m-1],x2a[0..n-1]偑梌偊傜傟偨偲偒
 //偙偺儖乕僠儞偼ya偺峴偺1師尦帺慠3師僗僾儔僀儞傪峔惉偟丄2奒摫娭悢傪攝楍
 //y2a[0..m-1][0..n-1]偵偄傟偰栠傞丅

 int j;

 for(j=0;j<m;j++)
         spline(x2a,ya[j],n,1.0e300,1.0e300,y2a[j]);

};

void splin2(double *x1a,double *x2a,double **ya,double **y2a,int m,int n,double x1,double x2,double *y)
{
 //偙偙偱偺x1a,x2a,ya,m,n偺堄枴偼splie2偺拲庍傪尒傛丅y2a偼splie2偑栠偡昞偱偁傞丅
 //x1,x2偼曗娫抣傪媮傔偨偄揰偺嵗昗偱偁傞丅偙傟傜傪偡傋偰梌偊傞偲丄偙偺儖乕僠儞偼
 //憃3師僗僾儔僀儞曗娫偱曗娫娭悢抣y傪寁嶼偟偰栠傞丅
 int j;
 double *ytemp,*yytemp;

 ytemp = new double [m];     //980817
 yytemp = new double [m];    //980817 debugged

 for(j=0;j<m;j++)
         splint(x2a,ya[j],y2a[j],n,x2,&yytemp[j]);
 spline(x1a,yytemp,m,1.0e300,1.0e300,ytemp);
 splint(x1a,yytemp,ytemp,m,x1,y);

 delete[] ytemp;
 delete[] yytemp;

};
void dy_splin2(double *x1a,double *x2a,double **ya,double **y2a,int m,int n,double x1,double x2,double *y)
{
 //偙偙偱偺x1a,x2a,ya,m,n偺堄枴偼splie2偺拲庍傪尒傛丅y2a偼splie2偑栠偡昞偱偁傞丅
 //x1,x2偼曗娫抣傪媮傔偨偄揰偺嵗昗偱偁傞丅偙傟傜傪偡傋偰梌偊傞偲丄偙偺儖乕僠儞偼
 //憃3師僗僾儔僀儞曗娫偱曗娫堦師旝暘抣y傪寁嶼偟偰栠傞丅
 int i,j;
 double *ytemp,*yytemp;

 ytemp = new double [n];
 yytemp = new double [n];

 double **y2a_temp,**ya_temp;

 y2a_temp=new double*[n];
 for (i=0;i<n;i++) y2a_temp[i]=new double[m];

 ya_temp=new double*[n];
 for (i=0;i<n;i++) ya_temp[i]=new double[m];

 for (i=0;i<m;i++) {
     for (j=0;j<n;j++){

         y2a_temp[j][i]=y2a[i][j];
         ya_temp[j][i]=ya[i][j];
     };
 };

 for (j=0;j<n;j++)
     splint(x1a,ya_temp[j] ,y2a_temp[j],m,x1,&yytemp[j]);
 spline(x2a,yytemp,n,1.0e300,1.0e300,ytemp);
 d_splint(x2a,yytemp,ytemp,n,x2,y);

 delete[] ytemp;
 delete[] yytemp;

 for (i=0;i<n;i++) delete[] y2a_temp[i];
 delete[] y2a_temp;
 for (i=0;i<n;i++) delete[] ya_temp[i];
 delete[] ya_temp;

};

void dx_splin2(double *x1a,double *x2a,double **ya,double **y2a,int m,int n,double x1,double x2,double *y)
{
 //偙偙偱偺x1a,x2a,ya,m,n偺堄枴偼splie2偺拲庍傪尒傛丅y2a偼splie2偑栠偡昞偱偁傞丅
 //x1,x2偼曗娫抣傪媮傔偨偄揰偺嵗昗偱偁傞丅偙傟傜傪偡傋偰梌偊傞偲丄偙偺儖乕僠儞偼
 //憃3師僗僾儔僀儞曗娫偱曗娫堦師旝暘抣y傪寁嶼偟偰栠傞丅
 int j;
 double *ytemp,*yytemp;

 ytemp = new double [m];
 yytemp = new double [m];

 for(j=0;j<m;j++)
         splint(x2a,ya[j],y2a[j],n,x2,&yytemp[j]);
 spline(x1a,yytemp,m,1.0e300,1.0e300,ytemp);
 d_splint(x1a,yytemp,ytemp,m,x1,y);

 delete[] ytemp;
 delete[] yytemp;

};

void spline(double *x,double *y,int n,double yp1,double ypn,double *y2)
{
//x[0],x[1],x[2],...傪枮偨偡揰偲娭悢抣y[i]=f(x[i])傪擖傟偨攝楍x[0..n-1],y[0..n-1]
//偑梌偊傜傟丄揰1,n偱偺曗娫娭悢偺1奒摫娭悢偺抣yp1,ypn偑梌偊傜傟偨偲偒丄偙偺儖乕僠儞偼
//奺揰x[i]偱偺曗娫娭悢偺2奒摫娭悢偺抣傪擖傟偨攝楍y2[0..n-1]傪曉偡丅傕偟yp1,ypn偺
//曅曽傑偨偼椉曽偑1e300埲忋側傜丄懳墳偡傞抂揰偱帺慠僗僾儔僀儞偺嫬奅忦審傪梫媮偟偨
//偙偲偲側傝丄偦偺抂揰偱偺2師摫娭悢偑0偵側傞丅
 int i,k;
 double p,qn,sig,un,*u;

 u = new double [n-1];
 if (yp1 > 0.99e300)
    y2[0] = u[0]= 0.0;
 else {
    y2[0] = -0.5;
    u[0] = (3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
 };

 for(i=1;i<=n-2;i++){
    sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
    p = sig*y2[i-1]+2.0;
    y2[i] = (sig-1.0)/p;
    u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
    u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
 };

 if (ypn > 0.99e300)
    qn = un = 0.0;
 else {
    qn =0.5;
    un =(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
 };
 y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);

for (k=n-2; k>=0; k--)
    y2[k]=y2[k]*y2[k+1]+u[k];

delete[] u;

};

void splint(double *xa,double *ya,double *y2a,int n,double x,double *y)
{
//彫偝偄弴偵暲傫偩昗杮揰xai偲娭悢抣傪擖傟傜傟偨攝楍xa[0..n-1],ya[0..n-1]
//偑梌偊傜傟丄忋偺spline偑弌椡偡傞攝楍y2a[0..n-1]偑梌偊傜傟丄x偺抣偑梌偊傜傟偨偲偒丄
//偙偺娭悢偼3師僗僾儔僀儞偵傛傞曗娫抣y傪曉偡丅

 int klo,khi,k;
 double h,b,a;

 klo = 0;
 khi = n-1;
 while (khi-klo>1){
    k=(khi+klo)>>1;
    if (xa[k] > x) khi=k;
    else klo = k;
 };

 h= xa[khi]-xa[klo];
 if (h==0.0) throw "Bad xa input routine splint";
 a=(xa[khi]-x)/h;
 b=(x-xa[klo])/h;
 *y=a*ya[klo]+b*ya[khi]+((a*a-1.)*a*y2a[klo]-(b*b-1.)*b*y2a[khi])*(h*h)/6.0;
};
void d_splint(double *xa,double *ya,double *y2a,int n,double x,double *y)
{
//彫偝偄弴偵暲傫偩昗杮揰xai偲娭悢抣傪擖傟傜傟偨攝楍xa[0..n-1],ya[0..n-1]
//偑梌偊傜傟丄忋偺spline偑弌椡偡傞攝楍y2a[0..n-1]偑梌偊傜傟丄x偺抣偑梌偊傜傟偨偲偒丄
//偙偺娭悢偼3師僗僾儔僀儞偵傛傞堦師旝暘抣y(dy/dx)傪曉偡丅

 int klo,khi,k;
 double h,b,a;

 klo = 0;
 khi = n-1;
 while (khi-klo>1){
    k=(khi+klo)>>1;
    if (xa[k] > x) khi=k;
    else klo = k;
 };

 h= xa[khi]-xa[klo];
 if (h==0.0) throw "Bad xa input routine d_splint";
 a=(xa[khi]-x)/h;
 b=(x-xa[klo])/h;
 *y=(ya[khi]-ya[klo])/h -(3.*a*a-1.)*y2a[klo]*h/6.0
                        +(3.*b*b-1.)*y2a[khi]*h/6.0;
};