proj.cpp

一个VC写的直接读取MIF文件的GIS系统

#include "stdafx.h"
#include "proj.h"
#include <math.h>

void CMapProj::rhumb_line1(double bb1, double ll1, double bb2, double ll2, double* angle, double* length )
{
	rhumb_line(bb1 * 0.0174532925199433, ll1 * 0.0174532925199433, bb2 * 0.0174532925199433, ll2 * 0.0174532925199433, angle, length );
	*angle *= 57.29577951308232;
	for(; *angle<0; *angle += 360);
	for(; *angle > 360; *angle -= 360);
}

// This subroutine is used for calculating angle and length of long line on the earth
void CMapProj::rhumb_line(double bb1, double ll1, double bb2, double ll2, double* angle, double* length )
{
	double q1,q2,angl,len,bb0,dbb,dll;

	bb0=(bb1+bb2)/2.0; 
	dbb=(bb2-bb1); 
	dll=(ll2-ll1);
	q1=q(bb1);
	q2=q(bb2);
	angl=arctg( dll, (q2-q1) );

	if( fabs(fabs(angl)-PAI/2.0) > PAI/180.0 )  
		{   len=m(bb0)*dbb/cos(angl); }
	else
		{   len=r(bb0)*dll/sin(angl); } 
	*angle=angl; *length=len;
}

// This subroutine is used for next point on the earth accoding to start point(b1,l1),angle and length
// (b1, l1) and angle in rad; length in meter;
void CMapProj::inverse_rhumb1( double b1,double l1,double angl,double len, double* b2,double *l2 )
{
	inverse_rhumb(b1 * 0.0174532925199433, l1 * 0.0174532925199433, angl * 0.0174532925199433, len,  b2, l2);
	*b2 *= 57.29577951308232;
	*l2 *= 57.29577951308232;
}

void CMapProj::inverse_rhumb( double b1,double l1,double angl,double len, double* b2,double *l2 )
{
	double sig0, sigi, dsig, sigii, err, qi, bi,li, dqi,q1, dli;
	short i=0;
	if( len < 1852 )
		err=0.001;
	else
		err=1.0;
	q1=q(b1);
	double tmp = invq(q1);
	sig0=len;
	sigi=len;
	sigii=len;
	do {
		sigii=sigi;

		dli=sigi/3437.746771/1852*sin(angl);
		li=l1+dli;

		dqi=sigi/3437.746771/1852*cos(angl);
		qi=q1+dqi;
		bi=invq(qi);

		rhumb_line(b1,l1,bi,li,&angl,&sigi);
		dsig=sig0-sigi;

		sigi=sigii+dsig;
		i=i+1;
	} while( dsig>err ) ;

	*b2=bi; *l2=li;

	return;
}
double CMapProj::arctg( double dy, double dx)
{
	if( dx==0.0 && dy==0.0 ) return(0.0);
	return ( atan2( dy, dx ) );
}

double CMapProj::atanhhh(double x)
{
	double y;
    y=log((1.0+x)/(1.0-x));
	y=y*0.5 ;
	return(y);
}

double CMapProj::q ( double bb )
{
	double tempq,tempe;
	tempe=sqrt(EARTH_E12);
	tempq=atanhhh(sin(bb))-tempe*atanhhh(tempe*sin(bb));
	return(tempq);
}

double CMapProj::invq(double q0 )
{
	double q1,q10,b0,dq;
	q10=q0;
	b0=asin(tanh(q0));
	q1=q(b0);
	do 
	{  dq=q0-q1;
	   q10=q10+dq;
	   b0=asin(tanh(q10));
	   q1=q(b0);
	}  while( (fabs(dq)) > 0.0000000001);
	return(b0);
}

double CMapProj::r( double b )
{
	 return(EARTH_A*cos(b)/sqrt(1.0-EARTH_E12*sin(b)*sin(b)));
}

double CMapProj::m( double bb )
{
	 double mm;
	 mm=sin(bb);
	 mm=sqrt(1.0-EARTH_E12*mm*mm);
	 return( EARTH_A*(1.0-EARTH_E12)/(mm*mm*mm) );
}