计算大地水准面高程序

计算大地水准面高程序，C语言编写，EGM96模型，和大家分享一下。。。

#include 
#include 
#include 
using namespace std;

double Arr[721];
int SignArr[181];
double T360[361];
double CosMLamta[361];
double SinMLamta[361];
double CnmSnm[65338][2];
double A360[361];
int IndCS[361];
const double a = 6378136.3;
const double GM = 3986004.415E+8;



void Init（）
{
	Arr[0] = 1;
	double Val=1;
	for（ int i=1; i<721; ++i）
		Arr[i]=i*Arr[i-1];

	SignArr[0]=1;
	for（ int j=1; j<181; ++j）
		SignArr[j]=（-1）*SignArr[j-1];
}

void InitT360（double t）
{
	T360[0]=1;
	for（ int i=1; i<361; ++i）
		T360[i]=t*T360[i-1];
}

void InitA360（double r）
{
	A360[0] = 1;
	A360[1] = a/r;
	for（ int n=2; n<=360; ++n）
		A360[n] = A360[n-1]*A360[1];
}

void InitCosSin（double lamta）
{
	for（ int m=0; m<361; ++m）
	{
		CosMLamta[m]=cos（m*lamta）;
		SinMLamta[m]=sin（m*lamta）;
	}
}

double Pnm（int nint m double t）
{
	double Ret = 1;
	Ret = Ret*pow（2n*（-1））;
	Ret = Ret*pow（1-t*t（double）m/2）;

	double Sum = 0;
	int r = （n-m）/2;
	for（int k=0; k<=r; ++k）
	{
		Sum+=SignArr[k]*Arr[2*n-2*k]/Arr[k]/Arr[n-k]/Arr[n-m-2*k]*T360[n-m-2*k];

		/*cout<<“Sum=“<		cout<<“Arr[2n-2k]=“<		cout<<“T360=“<	}
	
	//cout<<“Sum=“<	Ret = Ret* Sum;

	return Ret;
}

void CreateIndexForCS（）
{
	IndCS[0] = 0;
	IndCS[1] = 0;
	IndCS[2] = 0;
	int Ind = 3;
	for（ ; Ind<=360; ++Ind）
		IndCS[Ind]=IndCS[Ind-1]+Ind;
}

double CaculateV（double rdouble thetadouble lamta）
{
	double t = cos（theta）;
	double Ret = GM/r;
	double EmSum = 1;
	for（ int n=2;n<=5;++n）
	{
		double InSum = 0;
		double Temp = 0;
		for（ int m=0; m<=n; ++m）
		{
			int Index = IndCS[n]+m;
			Temp = CnmSnm[Index][0]*CosMLamta[m]+CnmSnm[Index][1]*SinMLamta[m];

			InSum += Temp*Pnm（nmt）;
			
            if（Temp>1e10）
    		{
    			cout<    			cout<<“Temp=“<    			cout<<“Pnm=“<    			//cout<<“A360=“<    		}
		}

		EmSum += InSum*A360[n];
	}

	return Ret*EmSum;
}

int main（）
{
	Init（）;
	CreateIndexForCS（）;

	double rthetalamta;
	cout<<“Input rthetalamta:“;
	cin>>r>>theta>>lamta;
	
	double t = cos（theta）;
	InitCosSin（lamta）;
	InitA360（r）;
	InitT360（t）;
    
   	ifstream infile（“egm96_to360.ascii“）;
	long LineCount =0;
	int Count =0;
	double temp =0;
	while（infile>>temp）
	{
		if（Count>=2）
			CnmSnm[LineCount][Count-2]=temp;
		if（Count==3）
			LineCount+=1;
	
		Count += 1;
		if（Count == 6）
			Count = 0;
	}
	infile.close（）;

	cout<	
    //system（“pause“）;

	return 0;
}

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2659  2011-05-13 20:09  EGM96\CalGravPotential.cpp

     文件     347388  2011-05-13 20:09  EGM96\Debug\CalGravPotential.obj

     文件          0  2011-05-13 20:09  EGM96\Debug\CalGravPotential.sbr

     文件     295936  2011-05-13 20:09  EGM96\Debug\EGM96.bsc

     文件     577619  2011-05-13 20:09  EGM96\Debug\EGM96.exe

     文件     834696  2011-05-13 20:09  EGM96\Debug\EGM96.ilk

    I.A....   2177648  2011-05-13 20:08  EGM96\Debug\EGM96.pch

     文件    1156096  2011-05-13 20:09  EGM96\Debug\EGM96.pdb

     文件     115712  2011-05-13 20:10  EGM96\Debug\vc60.idb

     文件     110592  2011-05-13 20:09  EGM96\Debug\vc60.pdb

     文件       4287  2009-10-26 20:04  EGM96\EGM96.dsp

     文件        535  2009-10-26 19:40  EGM96\EGM96.dsw

     文件      41984  2011-05-13 20:12  EGM96\EGM96.ncb

     文件      48640  2011-05-13 20:12  EGM96\EGM96.opt

     文件       1048  2011-05-13 20:09  EGM96\EGM96.plg

     文件    5292378  2009-03-25 09:01  EGM96\egm96_to360.ascii

     文件         67  2011-05-13 20:13  EGM96\程序说明.txt

     目录          0  2011-05-13 20:09  EGM96\Debug

     目录          0  2011-05-13 20:12  EGM96

----------- ---------  ---------- -----  ----

             11007285                    19