airmodel.cpp

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

#include "AirModel.h"
#include<math.h>
#include<stdio.h>


/* 1976 COESA atmosphere model */

//分段数
#define NUM1976PTS 8

//海拔向量
static double altitude76[NUM1976PTS] = {  /* in meters (m) */
    0.0, 11000.0, 20000.0, 32000.0, 47000.0, 51000.0, 71000.0, 84852.0 };
//温度梯度向量
static double tempGradient76[NUM1976PTS] = {  /* in K/m  */
    (-0.0065), 0.0, 0.0010, 0.0028, 0.0, -0.0028, -0.0020, -0.0020 };
//温度向量
static double temperature76[NUM1976PTS];
//压力向量
static double pressureRatio76[NUM1976PTS];

/*==================================================================*/
//大气常数
#define PRESSURE0        101325.0     /*  N/m^2                  */
#define TEMPERATURE0     288.15       /*  K                      */
#define GRAV_CONST       9.80665      /*  m/s^2                  */
#define MOL_WT           28.9644      /*  kg/kgmol (air)         */
#define R_HAT            8314.32      /*  J/kgmol.K (气体常数)   */
#define GAMMA            1.4          /*  (specific heat ratio)  */
#define GMR              ( GRAV_CONST * MOL_WT / R_HAT ) 
/*===================================================================*/

/*====================================================================
  函数: InitCalcAtmosCOESA 
  功能: 初始化温度和压力表
======================================================================*/
void InitCalcAtmosCOESA(void)
{
    int k;

    temperature76[0]   = TEMPERATURE0;
    pressureRatio76[0] = 1.0;

    for (k=0; k<(NUM1976PTS-1); k++) 
	{
        if (tempGradient76[k] != 0.0) 
		{
            temperature76[k+1]   = temperature76[k] + 
                tempGradient76[k]*(altitude76[k+1] - altitude76[k]);
            pressureRatio76[k+1] = pressureRatio76[k] *
                exp( log(temperature76[k]/temperature76[k+1]) * 
                     GMR/tempGradient76[k] );
        } 
		else 
		{
            temperature76[k+1]   = temperature76[k];
            pressureRatio76[k+1] = pressureRatio76[k] *
              exp( (-GMR)*(altitude76[k+1] - altitude76[k])/temperature76[k] );
        }
    }
}

/*====================================================================
  函数: CalcAtmosCOESA 
  功能: 给定海拔高度，计算大气密度和当地声速
  参数:
     输入参数:  
	     altitude: 海拔高度
     输出参数:
	     density: 大气密度
         speedofsound: 当地声速
         pressure: 大气压强
======================================================================*/
void CalcAtmosCOESA(double altitude, 
		            double &density,
                    double &speedofsound,
					double &pressure)
{
	double temp=0; 

	int bottom = 0;
	int top    = NUM1976PTS-1;
	int idx;
	
	if (altitude<= altitude76[bottom]) 
	{
		idx = bottom;
	} 
	else if (altitude>= altitude76[top]) 
	{
		idx = NUM1976PTS-2;
	} 
	else 
	{
		for (;;)
		{
			idx = (bottom + top)/2;
			if (altitude< altitude76[idx]) 
			{
				top = idx - 1;
			} 
			else if (altitude>= altitude76[idx+1])
			{
				bottom = idx + 1;
			} 
			else 
			{
				break;
			}
		} 
	}

	if ( tempGradient76[idx] != 0.0 ) 
	{
		temp = temperature76[idx] + 
			tempGradient76[idx] * (altitude - altitude76[idx]);
		pressure = PRESSURE0 * pressureRatio76[idx] * 
			pow(temperature76[idx]/temp, GMR/tempGradient76[idx] );
	} 
	else 
	{
		temp = temperature76[idx];
		pressure = PRESSURE0 * pressureRatio76[idx] *
			exp( (-GMR)*(altitude - altitude76[idx]) 
			/ temperature76[idx] );
	}
	density = pressure / ((R_HAT/MOL_WT)*temp);
	speedofsound = sqrt(GAMMA*temp*(R_HAT/MOL_WT));
}