//whihr_radproj.h 头文件 #if !defined(AFX_RADPROJECTION_H__B0D62835_14C8_41CB_A125_5FE50444D39E__INCLUDED_)#define AFX_RADPROJECTION_H__B0D62835_14C8_41CB_A125_5FE50444D39E__INCLUDED_#if !defined (__PI_VALUE__)#define __PI_VALUE__#define PI 3.1415927#endif//-------------------- 地球坐标系定位计算 ---------------------------//根据经纬度值来计算方位和地心角void geo2polar(double phy,do
uble delta_numda,do
uble &beta,do
uble &theta);void polar2geo(double beta,do
uble theta,do
uble &phy,do
uble &delta_numda);void arcS_of_gps(double Az,do
uble S,double& lon,do
uble& lat);void gps_to_Azimuth(double lon,do
uble lat,do
uble& Az,do
uble& beta);void gps_to_arcS(double lon,do
uble lat,do
uble& Az,do
uble& S);void gps_to_arcS(double lon,double lat, int& sx, int& sy);//-------------------- 等效地球代替地球半径近似定位计算 --------------------------- void Initialized(double dNdz, int x0, int y0, int res); //根据雷达回波坐标计算目标物的地理坐标 void tilt_of_gps(double e,do
uble dl,do
uble& dh,do
uble& dS); void tilt_of_gps(double e, do
uble az,do
uble dl,/ double& lon,do
uble& lat,double& dh,double& dS ); void tilt_of_gps(double e, int i, int j,do
uble& lon,do
uble& lat,/ double& dh,do
uble& dS) ; //根据经纬度值和仰角来计算方位、距离 void gps_to_tilt(double lon,do
uble lat,double e,do
uble& Az,do
uble& L,double& H ); void gps_to_tilt(double lon,double lat,double e,do
uble& x,do
uble& y); void DEM_to_tilt(double lon,do
uble lat,do
uble h, / double& ele,do
uble& az ,double& L,double& S) ; void ES_of_H(double e,do
uble dS,do
uble& dh); void HS_of_E(double H,do
uble dS,do
uble& ele);//根据计算波束高度 void RayPath_of_SH(double e, int Bins, short *S, short *H ); //constructor void CRadEPrj(double lon,do
uble lat,do
uble h);#endif // !defined(AFX_RADPROJECTION_H__B0D62835_14C8_41CB_A125_5FE50444D39E__INCLUDED_)//RadProjection.cpp cpp文件// RadProjection.cpp: implementation of the CRadEPrj class.////////////////////////////////////////////////////////////////////////#include <stdlib>#include <stdio>#include <math.h>#include "whihr_radproj.h"double r_numda, r_phy, r_h;double EarthRadius=6371.1 , r_res= 1.000 ;int rx= 0, ry= 0;double Km;//根据经纬度值来计算方位和地心角void geo2polar(double phy,do
uble delta_numda,do
uble &beta,do
uble &theta){ double cos_beta = sin(phy)*sin(r_phy)+cos(phy)*cos(r_phy)*cos(delta_numda); double sin_beta = sqrt(1-cos_beta*cos_beta); if( sin_beta>1 ) sin_beta=1; else
if (sin_beta<-1) sin_beta=-1; beta = asin(sin_beta); double sin_theta = cos(phy)*sin(delta_numda)/sin_beta; if( sin_theta>1 ) sin_theta = 1; else
if (sin_theta<-1) sin_theta = -1; theta = asin(sin_theta); double lat90 = asin(cos_beta*sin(r_phy)); //方位判断 if( phy<lat90) theta = PI-theta; else
if(delta_numda<0) theta = 2*PI+theta; return ;}void polar2geo(double beta,do
uble theta,do
uble &phy,do
uble &delta_numda){ double sin_phy = cos(beta)*sin(r_phy)+sin(beta)*cos(r_phy)*cos(theta); double cos_phy = sqrt(1-sin_phy*sin_phy); phy = asin(sin_phy); double sin_numda=sin(theta)*sin(beta)/cos_phy; if( sin_numda>1 ) sin_numda=1; else
if (sin_numda<-1) sin_numda=-1; delta_numda = asin(sin_numda); return ;}void arcS_of_gps(double az,do
uble S,double& lon,do
uble& lat){ double beta = S/EarthRadius; double theta = az*PI/180.0; double phy, delta_numda; polar2geo(beta, theta, phy, delta_numda); lat = phy*180/PI; lon = (delta_numda+r_numda)*180/PI; return ;}void gps_to_Azimuth(double lon,do
uble lat,do
uble& az,do
uble& beta){ double phy = lat*PI/180.0; double delta_numda = lon*PI/180.0-r_numda; double theta; geo2polar(phy,delta_numda,beta, theta); az = theta*180.0/PI; return ;}void gps_to_arcS(double lon,do
uble lat,do
uble& az,do
uble& S){ gps_to_Azimuth( lon, lat, az, S); S *= EarthRadius; return ;}void gps_to_arcS(double lon,double lat, int& sx, int& sy){ double phy = lat*PI/180.0; double delta_numda = lon*PI/180.0-r_numda; double beta, theta; geo2polar(phy,delta_numda,beta, theta); sx=int(beta*EarthRadius*sin(theta)/r_res+rx+0.5); sy=int(ry-beta*EarthRadius*cos(theta)/r_res+0.5); return ;}//constructor//-------------------- 等效地球代替地球半径近似定位计算 ------------------------------//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////void CRadEPrj(double lon,do
uble lat,do
uble h){ r_numda = lon*PI/180.0;;
//unit:弧度 r_phy = lat*PI/180.0;
//unit:弧度 r_h = h*0.001;}void Initialized(double dNdz, //线性梯度 N/Km int x0, int y0, //雷达站参考坐标 int res// 格点分辨率(米) ){ r_res = res*0.001 ;
//Unit:Kilometer rx = x0;
ry = y0 ;
//radar position Km = 1/(1-EarthRadius*dNdz*0.000001) ; Km=4.0/3.0 ;}//根据雷达回波坐标计算目标物的地理坐标//根据雷达回波坐标计算目标物的地理坐标void tilt_of_gps(double e, //仰角度数 do
uble dl, //雷达探测距离(千米) do
uble& dh, //拔海高度返回值(千米) do
uble& dS //水平距离 ){ if(dl < r_res) { dh=r_h; dS=dl; return ; } double Rm=Km*EarthRadius;
//等效地球半径 double alpha = e*PI/180.0; double beta=Km*atan(dl*cos(alpha)/(Rm+r_h+dl*sin(alpha))); dS=EarthRadius*beta;
//水平距离 dh = sqrt((Rm+r_h)*(Rm+r_h)+dl*dl+2*(Rm+r_h)*dl*sin(alpha))-Rm;
//高度 return ;}void tilt_of_gps(double e, //仰角度数 do
uble az, //雷达方位度数 do
uble dl, //雷达探测距离(千米) do
uble& lon, //经度返回值 do
uble& lat, //纬度返回值 do
uble& dh, //拔海高度返回值(千米) do
uble& dS //水平距离 ){ tilt_of_gps(e, dl, dh, dS ); if(dl < r_res) { lon = r_numda*180.0/PI ; lat = r_phy*180.0/PI ; return ; } arcS_of_gps(az, dS, lon, lat); return ;}//根据回波格点坐标计算目标物的地理坐标void tilt_of_gps(do
uble e, //仰角度数 int i, //东西向格点坐标 int j, //南北向格点坐标 do
uble& lon, //经度返回值 do
uble& lat, //纬度返回值 do
uble& dh, //拔海高度返回值 do
uble& dS //经度返回值 ){ do
uble dl = sqrt((i-rx)*(i-rx)+(ry-j)*(ry-j))*r_res;
do
uble az = asin((i-rx)*r_res/dl) ; if( i >= rx && j <= ry ) az=az*180.0/PI; else
if( i < rx && j <= ry ) az = 360.0+az*180.0/PI; else
az = 180.0-az*180.0/PI; tilt_of_gps(e, az, dl, lon, lat, dh, dS);}//根据地面距离计算波束高度void ES_of_H(do
uble e, //elevation angle do
uble dS, //经度返回值 do
uble& dh //return value og height ){ if( dS <= r_res ) dh = r_h; else
{ double alpha = e*PI/180.0; double Rm=EarthRadius*Km ; double betam = dS/Rm; dh = (Rm+r_h)*cos(alpha)/cos(alpha+betam)-Rm;
//高度 }}//根据水平距离和高度计算仰角void HS_of_E(double H,do
uble dS,do
uble &ele){ double Rm = EarthRadius*Km; double beta = dS/Rm; double cos_beta = cos(beta); if(cos_beta>1.0) cos_beta=1.0; double sin_beta=sin(beta); double tan_alpha = (cos_beta*(Rm+H)-(Rm+r_h))/(sin_beta*(Rm+H)); ele = atan(tan_alpha)*180.0/PI; return ;}//根据经纬度值和仰角来计算方位、距离void gps_to_tilt(double lon, //经度值 do
uble lat, //纬度值 do
uble e, //仰角度数 do
uble& Az, //雷达方位度数返回值 do
uble& L, //雷达探测距离返回值(千米) do
uble& H //高度返回值(千米) ){ double beta0; gps_to_Azimuth(lon, lat, Az, beta0); beta0 = beta0/Km; double alpha = e*PI/180.0; L=(Km*EarthRadius+r_h)*sin(beta0)/cos(alpha+beta0); double Rm=EarthRadius*Km ; H = sqrt((Rm+r_h)*(Rm+r_h)+L*L+2*(Rm+r_h)*L*sin(alpha))-Rm;
//高度// H=r_h+L*sin(alpha)+L*L*cos(alpha)*cos(alpha)/(2*Rm) ;}void gps_to_tilt(double lon, //经度值 do
uble lat, //纬度值 do
uble e, //仰角度数 do
uble &x,do
uble &y //返回的格点坐标 ) { double Az, L, H; gps_to_tilt(lon,lat,e, Az,L,H ); double sin_az=sin(Az*PI/180); double cos_az=cos(Az*PI/180); x=L*sin_az/r_res; y=L*cos_az/r_res;}//根据经纬度和拔海高度计算仰角、方位和距离;void DEM_to_tilt(double lon, //经度值 double lat, //纬度值 double H, //高度值(千米) double &e, //仰角返回度数 double &Az, //雷达方位度数返回值 double &L, //雷达探测距离返回值(千米) double &S //水平距离返回值(千米) ){ double beta; gps_to_Azimuth(lon,lat,Az, beta);
S=beta*EarthRadius; if( H<=r_h) //没有遮挡 { e=0; L=S; return ; } beta=beta/Km; double Rm=Km*EarthRadius; double cos_beta = cos(beta); if(cos_beta>1.0) cos_beta=1.0; double sin_beta=sin(beta); e=atan(((Rm+H)*cos_beta-(Rm+r_h))/((Rm+H)*sin_beta));
L=(Rm+H)*sin_beta/cos(e); e=e*180.0/PI ;}//根据计算波束高度void RayPath_of_SH(double e, //elevation angle int Bins, //total numbers of range gates short *S, short *H )//return value og height { double dh,dS; double dl=0;
//雷达探测距离(千米) for( int i=0;i<Bins;
++i) { tilt_of_gps(e, dl, dh, dS ); S=short(dS/r_res+0.5); H=short(dh*1000+0.5); dl += r_res; } return;}你是VC编译的吧。要引出那几个函数,自己加一下就行了。