ellipsoid.cpp

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

#if !defined(__ellipsoid_cpp)
#define __ellipsoid_cpp

//Difinition of ellipsoid class    19990318

#include <math.h>
#pragma hdrstop

#include "ellipsoid.h"

////Construcor and destructor
ellipsoid::ellipsoid(){

	cx=134.6465624; //[mm]
	cy=134.6465624; //[mm]
	cz=1009.024131; //[mm]

   dx=0.; //[mm]
   dy=-53.7931; //[mm]
   dz=-925.; //[mm]

   rx_design=13.; //rotation around x-axis(form)  [deg]
   rx_setup =45.;//ritation araound x-asis(setup)  [deg]

   h_plane=11.96;//=11.6-1.64+2.;

   z_top=34.;//=38.86-5;
};

ellipsoid::~ellipsoid(){

};
//equations for calcuration
inline bool ellipsoid::InMirrorArea(double x,double y){

   return  InMirrorBody(x,y,PlaneZ(x,y));

};

inline bool ellipsoid::InMirrorBody(double x,double y,double z){

    double t=(rx_setup-rx_design)*M_PI/180.;

    double xo=x;
    double yo=y*cos(t)+z*sin(t)+dy;
    double zo=-y*sin(t)+z*cos(t)+dz;

    return ((xo/cx)*(xo/cx)+(yo/cy)*(yo/cy)+(zo/cz)*(zo/cz)<=1);
};

inline int ellipsoid::AreaJudgement(double x,double y, double &z){

   //儚乕僋偺僄儕傾傪敾掕偟丄僄儕傾ID偲偦偺揰偱偺Z嵗昗傪曉偡
   //AreaID=   0 :Mirror
   //          1 :Plane
   //          2 :ztop

   int judge;

   if(InMirrorArea(x,y)) {

       z=MirrorZ(x,y);
       judge=0;

   }else{

       z=PlaneZ(x,y);
       judge=1;
   };


   if (z>=z_top) {

           z=z_top;
           judge=2;
   };

   return judge;
};

inline double ellipsoid::PlaneZ(double x,double y){

    double b=rx_setup*M_PI/180.;
    if (b==0) {
       return 1000.;
    }else{
       return ((h_plane-y*cos(b))/sin(b));
    };
};

inline double ellipsoid::MirrorZ(double x,double y){

   double a=cx;
   double b=cy;
   double c=cz;
   double t=(rx_setup-rx_design)*M_PI/180.;

   double z=-(2*pow(a,2)*(pow(c,2)*dy*sin(t) + cos(t)*(pow(b,2)*dz - (pow(b,2) - pow(c,2))*y*sin(t))) +
      sqrt(-4*pow(a,2)*(pow(b,2)*pow(cos(t),2) + pow(c,2)*pow(sin(t),2))*
         (-(pow(a,2)*pow(b,2)*pow(c,2)) + pow(a,2)*pow(c,2)*pow(dy,2) + pow(a,2)*pow(b,2)*pow(dz,2) + pow(b,2)*pow(c,2)*pow(x,2) +
           2*pow(a,2)*pow(c,2)*dy*y*cos(t) + pow(a,2)*pow(c,2)*pow(y,2)*pow(cos(t),2) - 2*pow(a,2)*pow(b,2)*dz*y*sin(t) +
           pow(a,2)*pow(b,2)*pow(y,2)*pow(sin(t),2)) + 4*pow(a,4)*
         pow(pow(c,2)*dy*sin(t) + cos(t)*(pow(b,2)*dz - (pow(b,2) - pow(c,2))*y*sin(t)),2)))/
   (2.*pow(a,2)*(pow(b,2)*pow(cos(t),2) + pow(c,2)*pow(sin(t),2)));

   return z;
};
inline double ellipsoid::MirrorgradX(double x,double y){
   double a=cx;
   double b=cy;
   double c=cz;
   double t=(rx_setup-rx_design)*M_PI/180.;

   double z=(2*pow(b,2)*pow(c,2)*x)/sqrt(-4*pow(a,2)*(pow(b,2)*pow(cos(t),2) + pow(c,2)*pow(sin(t),2))*
      (-(pow(a,2)*pow(b,2)*pow(c,2)) + pow(a,2)*pow(c,2)*pow(dy,2) + pow(a,2)*pow(b,2)*pow(dz,2) + pow(b,2)*pow(c,2)*pow(x,2) +
        2*pow(a,2)*pow(c,2)*dy*y*cos(t) + pow(a,2)*pow(c,2)*pow(y,2)*pow(cos(t),2) - 2*pow(a,2)*pow(b,2)*dz*y*sin(t) +
        pow(a,2)*pow(b,2)*pow(y,2)*pow(sin(t),2)) + 4*pow(a,4)*pow(pow(c,2)*dy*sin(t) + cos(t)*(pow(b,2)*dz - (pow(b,2) - pow(c,2))*y*sin(t)),2));

   return z;

};
inline double ellipsoid::MirrorgradY(double x,double y){
   double a=cx;
   double b=cy;
   double c=cz;
   double t=(rx_setup-rx_design)*M_PI/180.;

   double z=-(-(pow(a,2)*(pow(b,2) - pow(c,2))*sin(2*t)) - (2*sqrt(2)*pow(a,4)*pow(b,2)*pow(c,2)*(y + dy*cos(t) - dz*sin(t)))/
       sqrt(pow(a,2)*pow(b,2)*pow(c,2)*(pow(a,2)*pow(b,2) + pow(a,2)*pow(c,2) - pow(a,2)*pow(dy,2) - pow(a,2)*pow(dz,2) -
           pow(b,2)*pow(x,2) - pow(c,2)*pow(x,2) - 2*pow(a,2)*pow(y,2) - 4*pow(a,2)*dy*y*cos(t) +
           (pow(a,2)*pow(b,2) - pow(a,2)*pow(c,2) - pow(a,2)*pow(dy,2) + pow(a,2)*pow(dz,2) - pow(b,2)*pow(x,2) + pow(c,2)*pow(x,2))*
            cos(2*t) + 4*pow(a,2)*dz*y*sin(t) + 2*pow(a,2)*dy*dz*sin(2*t))))/(2.*pow(a,2)*(pow(b,2)*pow(cos(t),2) + pow(c,2)*pow(sin(t),2)));

   return z;

};
inline double ellipsoid::MirrorgradXX(double x,double y){
   double a=cx;
   double b=cy;
   double c=cz;
   double t=(rx_setup-rx_design)*M_PI/180.;

   double z=(sqrt(2)*pow(a,4)*pow(b,4)*pow(c,4)*(pow(b,2) + pow(c,2) - pow(dy,2) - pow(dz,2) - 2*pow(y,2) - 4*dy*y*cos(t) +
       (pow(b,2) - pow(c,2) - pow(dy,2) + pow(dz,2))*cos(2*t) + 4*dz*y*sin(t) + 2*dy*dz*sin(2*t)))/
   pow(pow(a,2)*pow(b,2)*pow(c,2)*(pow(a,2)*pow(b,2) + pow(a,2)*pow(c,2) - pow(a,2)*pow(dy,2) - pow(a,2)*pow(dz,2) - pow(b,2)*pow(x,2) -
       pow(c,2)*pow(x,2) - 2*pow(a,2)*pow(y,2) - 4*pow(a,2)*dy*y*cos(t) +
       (pow(a,2)*pow(b,2) - pow(a,2)*pow(c,2) - pow(a,2)*pow(dy,2) + pow(a,2)*pow(dz,2) - pow(b,2)*pow(x,2) + pow(c,2)*pow(x,2))*cos(2*t) +
       4*pow(a,2)*dz*y*sin(t) + 2*pow(a,2)*dy*dz*sin(2*t)),1.5);

   return z;

};
inline double ellipsoid::MirrorgradYY(double x,double y){
   double a=cx;
   double b=cy;
   double c=cz;
   double t=(rx_setup-rx_design)*M_PI/180.;

   double z=(2*sqrt(2)*pow(a,4)*pow(b,4)*pow(c,4)*(pow(a,2) - pow(x,2)))/
   pow(pow(a,2)*pow(b,2)*pow(c,2)*(pow(a,2)*pow(b,2) + pow(a,2)*pow(c,2) - pow(a,2)*pow(dy,2) - pow(a,2)*pow(dz,2) - pow(b,2)*pow(x,2) -
       pow(c,2)*pow(x,2) - 2*pow(a,2)*pow(y,2) - 4*pow(a,2)*dy*y*cos(t) +
       (pow(a,2)*pow(b,2) - pow(a,2)*pow(c,2) - pow(a,2)*pow(dy,2) + pow(a,2)*pow(dz,2) - pow(b,2)*pow(x,2) + pow(c,2)*pow(x,2))*cos(2*t) +
       4*pow(a,2)*dz*y*sin(t) + 2*pow(a,2)*dy*dz*sin(2*t)),1.5);

   return z;

};

inline double ellipsoid::GetZ( double x, double y){

   double z;

   AreaJudgement(x,y,z);

   return z;

   //z=0 ,when x=y=0.

};

inline vector ellipsoid::GetPosition( double x, double y){
	return ( vector( x, y, GetZ(x,y) ) );
};

inline double ellipsoid::gradX(double x, double y){

   double z;
   double value;

   switch (AreaJudgement(x,y,z)){

       case 0:           //Mirror
           value = MirrorgradX(x,y);
           break;

       case 1:           //Plane
           value = 0;
           break ;

       case 2:           //ztop
           value = 0;
           break ;

       default:
           value = 0;

   };

   return value;

};

inline double ellipsoid::gradY(double x, double y){

   double z;
   double value;

   double b;

   switch (AreaJudgement(x,y,z)){

       case 0:
           value = MirrorgradY(x,y);
           break;

       case 1:
           b=rx_setup*M_PI/180.;

           if (tan(b)!=0){
               value= -1./tan(b);     //PlanegradY
           }else{
               value= 0;
           };

           break ;

       case 2:
           value = 0;   //ztop
           break ;

       default:
           value = 0;

   };

   return value;

};

inline double ellipsoid::gradXX(double x, double y){

   double z;
   double value;

   switch (AreaJudgement(x,y,z)){

       case 0:           //Mirror
           value = MirrorgradXX(x,y);
           break;

       case 1:           //Plane
           value = 0;
           break ;

       case 2:           //ztop
           value = 0;
           break ;

       default:
           value = 0;

   };

   return value;

};

inline double ellipsoid::gradYY(double x, double y){


   double z;
   double value;

   switch (AreaJudgement(x,y,z)){

       case 0:           //Mirror
           value = MirrorgradYY(x,y);
           break;

       case 1:           //Plane
           value = 0;
           break ;

       case 2:           //ztop
           value = 0;
           break ;

       default:
           value = 0;

   };

   return value;
};

// This is the normal vector of surface on the point(x,y).
// This is toward the inside of ellipsoid.
// This.length is 1.
vector ellipsoid::NormalVector ( double x, double y){
	return ( vector( -gradX(x,y), -gradY(x,y), 1).normalize() );  //OK 980520
};
vector ellipsoid::NormalVector ( vector position){
	return ( NormalVector( position.getX(), position.getY() ) );
};

////GradientVector////////////////////////////////////////////
// This is one of tangential vactors on the tangential plane.
// This is toward direction which has the largest gradient.
// This length is 1.
vector ellipsoid::GradientVector(double x, double y){
	double gx = gradX( x,y);
	double gy = gradY( x,y);

	return ( vector( gx, gy, gx*gx+gy*gy ).normalize() );
};
vector ellipsoid::GradientVector( vector position){
	return( GradientVector( position.getX(), position.getY() ) );
};

////NoGradientVector///////////////////////////////////////////
// This is one of tangential vactors on the tangential plane.
// This is toward direction which has no gradient.
//// This length will be 1,
//// and this Z-coordinate will be Zero.
vector ellipsoid::NoGradientVector(double x,double y){
	return( NormalVector(x,y)%GradientVector(x,y) );
};
vector ellipsoid::NoGradientVector(vector position){
	return( NormalVector(position)%GradientVector(position) );
};


int ellipsoid::WriteFaceName(ofstream& fout){

	fout << "(Face name is : ellipsoid);" << endl;
	return(1);
};
int ellipsoid::WriteParameters(ofstream& fout){

	fout << "(cx[mm] is :" << cx <<");" << endl;
	fout << "(cy[mm] is :" << cy <<");" << endl;
	fout << "(cz[mm] is :" << cz <<");" << endl;

	fout << "(dx[mm] is :" << dx <<");" << endl;
	fout << "(dy[mm] is :" << dy <<");" << endl;
	fout << "(dz[mm] is :" << dz <<");" << endl;

	fout << "(rx_design is[deg] :" << rx_design <<");" << endl;
	fout << "(rx_setup[deg] is :" << rx_setup <<");" << endl;

	fout << "(h_plane[mm] is :" << h_plane <<");" << endl;
	fout << "(z_top[mm] is :" << z_top <<");" << endl;

	return(1);
};


void ellipsoid::SetParameter(TStringList *str)
{
   for (int i=0; i<str->Count; i++){

         AnsiString temp;

         temp = "cx";
         if (str->Names[i]==temp){
              cx = str->Values[str->Names[i]].ToDouble();
         };
         temp = "cy";
         if (str->Names[i]==temp){
              cy = str->Values[str->Names[i]].ToDouble();
         };
         temp = "cz";
         if (str->Names[i]==temp){
              cz = str->Values[str->Names[i]].ToDouble();
         };

         temp = "dx";
         if (str->Names[i]==temp){
              dx = str->Values[str->Names[i]].ToDouble();
         };
         temp = "dy";
         if (str->Names[i]==temp){
              dy = str->Values[str->Names[i]].ToDouble();
         };
         temp = "dz";
         if (str->Names[i]==temp){
              dz = str->Values[str->Names[i]].ToDouble();
         };

         temp = "rx_design";
         if (str->Names[i]==temp){
              rx_design = str->Values[str->Names[i]].ToDouble();
         };
         temp = "rx_setup";
         if (str->Names[i]==temp){
              rx_setup = str->Values[str->Names[i]].ToDouble();
         };

         temp = "h_plane";
         if (str->Names[i]==temp){
              h_plane = str->Values[str->Names[i]].ToDouble();
         };
         temp = "z_top";
         if (str->Names[i]==temp){
              z_top = str->Values[str->Names[i]].ToDouble();
         };

         //isConcave
   };
   face::SetFaceParameter(str);
};

void ellipsoid::GetParameter(TStringList *str)
{

   str->Add(AnsiString("---ellipsoid parameter(s)---"));
   str->Add(AnsiString("FaceType=ellipsoid"));

   str->Add(AnsiString("---(x/cx)^2+(y/cy)^2+(z/cz)^2=1---"));

   str->Add(AnsiString("cx=")+ cx);
   str->Add(AnsiString("cy=")+ cy);
   str->Add(AnsiString("cz=")+ cz);

   str->Add(AnsiString("dx=")+ dx);
   str->Add(AnsiString("dy=")+ dy);
   str->Add(AnsiString("dz=")+ dz);

   str->Add(AnsiString("rx_design=")+ rx_design);
   str->Add(AnsiString("rx_setup=")+ rx_setup);

   str->Add(AnsiString("h_plane=")+ h_plane);
   str->Add(AnsiString("z_top=")+ z_top);

   face::GetFaceParameter(str);
};

#endif