definefunction.h

计算子午圈半径M 计算卯酉圈半径N 计算玮圈半径r 计算经线弧长Sm 计算纬线弧长Sn 计算球面梯形面积 计算最大角度变形 计算球面极坐标

#include <process.h>
#include <fstream.h>
#include <conio.h>
#include <stdlib.h>
# include <math.h>
double rad()
{
	int a,b=0,c=0;
	cout<<"请输入角度： 度/分/秒\n";
	cin>>a;
	cout<<"度\n";
	cin>>b;
	cout<<"分\n";
	cin>>c;
	cout<<"秒\n";
	double rad;
	rad=a*PI/180+b*PI/10800+c*PI/648000;
	cout<<rad<<"rad"<<endl;
	return rad;
}
double rad(int a)
{
	double rad;
	rad=a*PI/180;
	return rad;
}
double rad1( double a)
{
	double rad1;
	rad1=a*180/PI;
	return rad1;
}	

double choice()
{
	char input;
	double a;
	cout<<"请选择角度/弧度，1为角度，2为弧度：\n";
	cin>>input;
	switch (input)
	{
	case '1': 
		a=rad();break;
	case '2':
		cout<<"请输入弧度：";
        cin>>a;
		break;
	}
	return a;
}
double M(double q)
{
	double M;
	double c;
	double d;
	d=1-e2*sin(q)*sin(q);
	c=pow(d,1.5);
	M=a*(1-e2)/c;
	return M;
}
double N(double q)
{
	double N;
	double c;
	double d;
	d=1-e2*sin(q)*sin(q);
	c=sqrt(d);
	N=a/c;
	return N;
}
double r(double q)
{
	double d,c,r;
	d=1-e2*sin(q)*sin(q);
	c=sqrt(d);
	r=a*cos(q)/c;
	return r;
}
double Sn(double r, double l)
{
	double Sn;
	Sn=r*l;
	return Sn;
}
double Sm(double rad) //用以计算赤道到rad的经线弧长
{
	double A=1.0050517739;
	double B=0.00506237764;
	double C=0.0000106245;
	double D=0.00000002081;
	double s;
	s=a*(1-e2)*(A*rad-B*sin(2*rad)/2+C*sin(4*rad)/4-D*sin(6*rad)/6);
	return s;
}

double Sm(double rad1, double rad2)
{
    double d=0;
	double rad;
	double f;
	double b,Drad;
	Drad=(rad2-rad1)/10000;
	for(int i=0;i<10000;i++)
	{
	rad=rad1+Drad*i;
	b=1/(1-e2*sin(rad)*sin(rad));
    f=pow(b,3/2);
	d+=f*Drad;
	}
	d=a*(1-e2)*d;
	return d;
}
double s(double rad1, double rad2,double M1,double M2)
{
    double d=0;
	double rad;
	double f;
	double b,Drad1,Drad2;
	Drad1=fabs(rad2-rad1)/1000;
    Drad2=fabs(M2-M1)/1000;
	for( int j=0;j<1000;j++)
	{
	for(int i=0;i<1000;i++)
	{
	 rad=rad1+Drad1*i;
	 f=cos(rad)/(1-e2*sin(rad)*sin(rad));
	 d+=f*Drad1;
	}
	d+=d*Drad2;
    }
	d=a*(1-e2)*d;
	return d;
}
double w(double a,double b)
{
	double w;
	w=acos(2*sqrt(a*b)/(a+b));
    w=2*w;
	return w;
}
double w(double Z)
{
	double w;
	w=asin((Z-sin(Z))/(Z+sin(Z)));
	w=2*w;
	return w;
}
position getp(double M0,double N0,double M,double N)
{
	position gp;
	if(M==M0&&N==N0)
	{
		gp.a=0;
		gp.Z=0;
	}
	else 
	{
	N=rad(N);
    M=rad(M);
	M0=rad(M0);
	N0=rad(N0);
    gp.Z=acos(sin(N)*sin(N0)+cos(N)*cos(N0)*cos(M-M0));
	if (sin(gp.Z)==0)
		gp.a=atan(sin(M-M0)/tan(N));
	else
	gp.a=asin(cos(N)*sin(M-M0)/sin(gp.Z));
	}
	return gp;
}
position0 getp(double a,double Z)
{
	position0 gp0;
	double R=63.710;
	gp0.a=a;
	gp0.p=R*Z;
	return gp0;
}
position1 getp(double a,double p,int)
{
	position1 gp1;
	gp1.x=p*cos(a);
	gp1.y=p*sin(a);
	return gp1;
}
double max(double a,double b)
{
	double max;
	max=a>b?a:b;
	return max;
}
double min(double a,double b)
{
	double min;
	min=a<b?a:b;
	return min;
}