插值.cpp

用来计算航天器再入质点弹道

#include"插值.h"
//**************************************************
Chazhi::Chazhi()
{
	fp		=	NULL;
	n		=	0;
	xnum	=	1;
	ynum	=	1;
	znum	=	1;
	xdata	=	NULL;
	ydata	=	NULL;
	zdata	=	NULL;
	value	=	NULL;
	data	=	0.0;
}
//***************************************************
bool Chazhi::readdata(char* filename)
{
	fp = fopen( filename , "r" ) ;
	
	if( fp == NULL )
	{
		cout<<"Can't open the file!"<<endl ;
		
		return false ;
	}
	
	char  str[20] ;
	
	fscanf( fp , "%s" , str ) ;
	
	if( !strcmp( str, "Dim" ) ) 
	{
        fscanf( fp, "%d" , &n ) ;
        
		switch(n)
		{
		  case 1:
			fscanf( fp , "%d" , &xnum ) ;
			
			xdata = new double[xnum];
			
			fscanf( fp , "%s" , str ) ;
			
			if( !strcmp( str , "Line" ))
			for( int i = 0 ; i < xnum ; i++ )
				fscanf( fp , "%lf" , xdata+i ) ;
				
			break;

		  case 2:
			fscanf( fp, "%d %d" , &xnum , &ynum ) ;
			
			xdata = new double[xnum];
			ydata = new double[ynum];
			
			fscanf( fp , "%s" , str ) ;
			
			if( !strcmp( str , "Line" ) ) 
			for( int i = 0 ; i < xnum ; i++ )
				fscanf( fp , "%lf" , xdata+i ) ;
				
			fscanf( fp , "%s" , str ) ;
			
			if( !strcmp( str , "Row" ) )
			for( int i = 0 ; i < ynum ; i++ )
				fscanf( fp , "%lf" , ydata+i ) ;
					
			break;

		  case 3:
			fscanf( fp, "%d %d %d" , &xnum , &ynum , &znum ) ;
			
			xdata = new double[xnum];
			ydata = new double[ynum];
			zdata = new double[znum];
			
			fscanf( fp , "%s" , str ) ;
			
			if( !strcmp( str , "Line" ) )
			for( int i = 0; i < xnum ; i++ )
				fscanf( fp , "%lf" , xdata+i ) ;
				
			fscanf( fp , "%s" , str ) ;
			
			if( !strcmp( str , "Row" ) )
			for( int i = 0 ; i < ynum; i++ )
				fscanf( fp , "%lf" , ydata+i ) ;
					
			fscanf( fp , "%s" , str ) ;
			
			if( !strcmp( str , "Pag" ) )
			for( int i = 0; i < znum; i++ )
				fscanf( fp , "%lf" , zdata+i ) ;
						
			break;
		}
		
	}
	
	fscanf( fp, "%s", str ) ;
	
	if( strcmp( str , "Data" ) )
	{
		cout<<"Data wrong!\n" ;
		return false ;
	}
	
	value = new double [xnum*ynum*znum] ;
	
	for( int i = 0; i < xnum*ynum*znum ; i++)
		fscanf( fp , "%lf" , value+i ) ;	
	
	fclose(fp) ;
	
	return true ;
	
}

//******************************************************************
double Chazhi::insert(char* filename,int dim,double x,double y,double z)
{  
	data	=	0.0;
	if(readdata(filename)==true)
	{
       if(dim != n)
		 cout<<"维数出现错误!"<<endl ;

       switch(n)
	   {
          case 1:
		    element1(x);
		    break;
	      case 2:
		    element2(x,y);
		    break;
	      case 3:
		    element3(x,y,z);
		    break;

	   }
	}

	return data;
}
//*******************************************************************
void Chazhi::location(double* ldata ,double l ,int* loc , int lnum  )
{
	int     k  = 0 ;
	double  t  = fabs( l - ldata[0] ) ;
	
	for(int i = 0 ; i < lnum ; i++)
	{
	   if( t >= fabs( l - ldata[i] ) )
	   {
	         t = fabs( l - ldata[i] ) ;
		     k = i ;
	   }
	}
		
	if( k == 0)
	{		
	   loc[0] = 0 ;
	   loc[1] = 1 ;
	   loc[2] = 2 ;
	}
	else 
	if( k == lnum -	1 )
	{
	   loc[0] = k - 2 ;
	   loc[1] = k - 1 ;
	   loc[2] = k ;
	} 
	else
	{
	   loc[0] = k - 1 ;
	   loc[1] = k ;
	   loc[2] = k + 1 ;			
	}
}
//*********************************************************************
void Chazhi::element1(double x)
{
   int* xloc = new int[3];
   
   location(xdata , x , xloc , xnum );
   
   double yy = 0.0 ;

   for(int i = 0 ; i < 3 ; i ++)
	{
		yy = *( value + xloc[i] ) ;
       
		for(int j = 0 ; j < 3 ; j ++)
		{					
			if( j == i ) continue ;					
			
			yy *= (( x - *(xdata+xloc[j]) )/( *(xdata+xloc[i]) - *(xdata+xloc[j]) ) ) ;
		}
		data += yy ;
	}
}
//*******************************************************************
void Chazhi::element2(double x,double y)
{	
	int* xloc = new int[3];
	int* yloc = new int[3];
	location(xdata , x , xloc , xnum );
	location(ydata , y , yloc , ynum );
	
	double zz	=	0.0;
    
	for(int i = 0 ; i < 3 ; i ++)
	{
		for(int j = 0 ; j < 3 ; j ++)
		{					
			zz = *(value + xnum*yloc[j] + xloc[i]) ;
			
			for(int k = 0 ; k < 3 ; k ++)
			{					
				if( k == j ) continue ;					
				zz *= ( (y - *(ydata+yloc[k]) )/( *(ydata+yloc[j]) - *(ydata+yloc[k]) ) ) ;
			}
			for( k = 0 ; k < 3 ; k ++)
			{
				if( k == i ) continue ;				
				zz *= ( (x - *(xdata+xloc[k]) )/( *(xdata+xloc[i]) - *(xdata+xloc[k]) ) ) ;   
			}
			data+=zz;
		}	
	} 

}
//*******************************************************
void Chazhi::element3(double x,double y,double z)
{
    int* xloc = new int[3];
	int* yloc = new int[3];
    int* zloc = new int[3];
	location(xdata , x , xloc , xnum );
	location(ydata , y , yloc , ynum );
    location(zdata , z , zloc , znum );
    
	double dd;
	for(int i=0;i<3;i++)
	{ 
    	for(int j=0;j<3;j++)
		{  
            
			for(int k=0;k<3;k++)
			{
			   dd = (*(value+zloc[k]*xnum*ynum+yloc[j]*xnum+xloc[i]));
			   
			   for(int m=0;m<3;m++)
               {
				   if( m == k ) continue ;
			       dd *= ( (z - *(zdata+zloc[m]) )/( *(zdata+zloc[k]) - *(zdata+zloc[m]) ) ) ;
               }

               for( m=0;m<3;m++)
               {
				   if( m == j )  continue;
			       dd *= ( (y - *(ydata+yloc[m]) )/( *(ydata+yloc[j]) - *(ydata+yloc[m]) ) ) ;
			   }

               for( m=0;m<3;m++)
			   {
                   if( m == i )  continue;
                   dd *= ( (x - *(xdata+xloc[m]) )/( *(xdata+xloc[i]) - *(xdata+xloc[m]) ) ) ;
			   }
               
               data+=dd;

			}	
		}  
	}			 
}