heightanomaly.cpp

计算高程异常比较好的方法

//---------------------------------------------------------------------------

#include <vcl.h>


#pragma hdrstop

#include "HA.h"


//---------------------------------------------------------------------------

#pragma package(smart_init)

#include <iostream.h>   //输入输出流体系
#include <fstream.h>   //定义了流类系中支持以外部文件为存储介质的标准流的操作类
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <iomanip.h>  //定义了若干个只接受一个入口参数的操作符(manipulator)



double THA::DuiShu(double H,double Hp,double Hr,double x,double y,double x0,double y0)
        {
		double m,n,r0;
		r0=sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));

		if(r0==0)    H=H+0.1;

		m=H-Hp;
		m+=sqrt(m*m+r0*r0);

		n=Hr-Hp;
		n+=sqrt(n*n+r0*r0);

		return log(m/n);

}


double THA:: JiFenTi( double * DEM, DEM_INFO head,double Hr,double Hp,double x0,double y0)
{
		//DEM_INFO head;
	    int D=99999999;
		int HangHP,LieHP;
		int HangMinY,HangMaxY,LieMinX,LieMaxX;

		LieHP= (int)((x0-head.MINX)/head.DistanceX); //左右
		HangHP=(int)((y0-head.MINY)/head.DistanceY);   //上下       //目标点的近似点的点号
	    HangMinY=HangHP-D;	//Y轴方向

        if(HangMinY>=head.NumberHang)
		    ShowMessage("HangMinY 不可用");
        else if(HangMinY<0)
			 {
            HangMinY=0;
			 }

		    HangMaxY=HangHP+D;
			if(HangMaxY<0)
			 {
             ShowMessage("HangMaxY 不可用");
             }
			 else if(HangMaxY>=head.NumberHang)
			 {
				 HangMaxY=head.NumberHang-1;
			 }

		    LieMinX=LieHP-D;
		   if(LieMinX>=head.NumberLie)
           {
		     ShowMessage("LieMinX 不可用");
           }
           else if(LieMinX<0)
			 {
				 LieMinX=0;
			 }
		LieMaxX=LieHP+D;

		if(LieMaxX<0)
		{
            ShowMessage("LieMaxX 不可用");  //cout<<""<<endl;
		}
        else if(LieMaxX>=head.NumberLie)
                LieMaxX=head.NumberLie-1;

		double sum=0;
		for(int j=HangMinY;j<=HangMaxY;j++)
		{
			for(int i=LieMinX;i<=LieMaxX;i++)
			{
			double H=DEM[j*head.NumberLie+i];
			if(H<head.MINZ||H>head.MAXZ)
				continue;

		double x=head.MINX+i*head.DistanceX;
		double y=head.MINY+j*head.DistanceY;

		sum+=DuiShu(H,Hp,Hr,x,y,x0,y0)*head.s;
			}
		}
		return (sum) ;

}
struct aim
{           double x0;
			double y0;
			double Hp;
}point[30]

void THA::calMain ()

{
        FILE *fp;
        double g,g0,h,Tc;
		double *DEM;

		double Hr;
		double Hp,x0,y0;

		double test_score;
		char ch[15];
        int i,j,k;

	   const double G=6.673E-8,L=2.67,k=0.72E-7,t=5.8E-6,pi(3.1415926);
            double sum=0;

 		ifstream test_file;
		test_file.open(InFileName.c_str(),ios::in);
			if(test_file.fail())
			{
				ShowMessage(InFileName+"不能打开");
				return;
            }    //end if
		test_file>>ch;
		cout<<"测试的DEM图幅号为"<<ch<<endl;           //READ "DSAA"

		test_file>>test_score;
		head.NumberLie=(int)test_score;   //行号
		cout<<"列号"<<head.NumberLie<<endl;

		test_file>>test_score;
		head.NumberHang=(int)test_score;   //列号
		cout<<"行号"<<head.NumberHang<<endl;

		test_file>>test_score;
		head.MINX=test_score;
		cout<<"x的最小值"<<head.MINX<<endl;
		test_file>>test_score;
		head.MAXX=test_score;
		cout<<"x的最大值"<<head.MAXX<<endl;     //X的最大值和最小值

		test_file>>test_score;
		head.MINY=test_score;
		cout<<"y的最小值"<<head.MINY<<endl;
		test_file>>test_score;
		head.MAXY=test_score;
		cout<<"y的最大值"<<head.MAXY<<endl;    //Y的最大值和最小值

		test_file>>test_score;
		head.MINZ=test_score;
		cout<<"z的最小值"<<head.MINZ<<endl;
		test_file>>test_score;	
		head.MAXZ=test_score;
		cout<<"z的最大值"<<head.MAXZ<<endl;    //Z的最大值和最小值
		///read file HEAD

		DEM=new double[head.NumberHang*head.NumberLie];
			if(DEM==NULL) 
			{
				cout<<"出现错误"<<endl;
				   return;
			}

		for(int j=0;j<head.NumberHang;j++)
		{
			for(int i=0;i<head.NumberLie;i++)
			{
						 test_file>>test_score;
						 //cout<<test_score<<endl;
						 DEM[j*head.NumberLie+i]=test_score;
			}
		}
		test_file.close();
		
		head.DistanceX=(head.MAXX-head.MINX)/(head.NumberLie-1);
		head.DistanceY=(head.MAXY-head.MINY)/(head.NumberHang-1);
		head.s=head.DistanceX*head.DistanceY;
		cout<<head.DistanceX<<"++++"<<head.DistanceY<<"+++"<<head.s<<"+++++++"<<endl;

        if((fp=fopen("point","w"))==NULL)
            {
            cout<<"cannot open files\n">>;
            exit(0);
            }
         for(k=0;k<30;k++)
            {
            fseek(fp,k*sizeof(struct aim),0);
            fscanf("%f,%f,%f,\n",point[k].x0,point[k].y0,point[k].HP);

            }
            fclose(fp);
            
		const double a;
            cout<<"请输入纬度a "<<endl;
	        cin>>a;
		    cout<<"请输入参考面高程Hr"<<endl;
            cin>>Hr;
        h=a*pi/180;
		g0=978.0327*(1+0.0053024*sin(h)*sin(h))-t*sin(2*h)*sin(2*h);
		/*////////////////////////////////////////////////////////////计算点值
		int m=5,n=5;
        double Hc=DEM[n*head.NumberLie+m];
        double xc=head.MINX+m*head.DistanceX;
	    double yc=head.MINY+n*head.DistanceY;
		cout<<":::Hc:::"<<Hc<<":::xc:::"<<xc<<":::yc:::"<<yc<<":::"<<endl;
		////////////////////////////////////////////////////////////////////*/
		double CANKAO[4]={0,1500,2500,3500};            //参考面高程
		for(int j=0;j<4;j++)
		{
			Hr=CANKAO[j];
			//cout<<"..................................................................."<<endl;
			//cout<<"参考面高程Hr为"<<Hr<<"时的各点高程异常为："<<endl;
		  for(int i=0;i<12;i++)
		  {
			  x0=point[i].x0;
			  y0=point[i].y0;
			  Hp=point[i].Hp;

		g=g0-0.3086*(Hp/1000)+k*(Hp/1000)*(Hp/1000);   //Hp---Km

	   sum=JiFenTi(DEM,head,Hr,Hp,x0,y0);


	   Tc=100*(G*L*sum)/g;
	   cout<</*"目标点"<<i<<"的高程异常是"<<*/Tc<<endl;
		  }
		}
}