#include<stdio.h>
#include<math.h>
#define a 6378245
#define e2 0.006693421622966
#define PI 3.1415926535898
#define z (PI/180)
#define lambda0 116
double fun1(double phi)/*计算U值*/
{
double psi,e,U;
e=sqrt(e2);
psi=asin(e*sin(phi*z));
U=(tan((45+phi/2)*z))/pow((tan((45+psi/2)*z)),e);
return U;
}
double fun2(double phi)/*计算r值*/
{
double r;
r=(a*cos(phi*z))/sqrt(1-e2*pow(sin(phi*z),2));
return r;
}
void fun(double phi,double lambda,double alpha,double K,double rhos,double *t)/*计算坐标*/
{
double delta,rho;
delta=alpha*(lambda-lambda0);
rho=K/pow(fun1(phi),alpha);
t[0]=rhos-rho*cos(delta*z);
t[1]=rho*sin(delta*z);
}
void main()
{
double phi1,phi2,phi,phin=0,phis=0,alpha,U1,U2,r1,r2,lambda1=0,lambda2=0,lambda=0,K,rhos;
long int M0;
double t[2]={0,0};
printf("请输入标准纬线φ1,φ2\t");
/*phi1=37.5;phi2=41;lambda1=112;lambda2=120;phin=36;phis=43;M0=5000000;*/
scanf("%lf%lf",&phi1,&phi2);
printf("请输入经线区间λw、λz(角度值,以度为单位)\t");
scanf("%lf%lf",&lambda1,&lambda2);
printf("请输入纬线区间φn、φs");
scanf("%lf%lf",&phin,&phis);
printf("请输入比例尺分母,如1比1000输入:1000\t");
scanf("%ld",&M0);
U1=fun1(phi1);
U2=fun1(phi2);
r1=fun2(phi1);
r2=fun2(phi2);
alpha=(log10(r1)-log10(r2))/(log10(U2)-log10(U1));
K=(r1*pow(U1,alpha)*100)/(alpha*M0);
printf("\t经纬度\t\t\t地图上坐标\n");
rhos=K/pow(fun1(phin),alpha);
for(phi=phin;phi<=phis;phi++)
for(lambda=lambda1;lambda<=lambda2;lambda++)
{fun(phi,lambda,alpha,K,rhos,t);
printf("%lf,%lf\t\t%lf,%lf\n",phi,lambda,t[0],t[1]);}
}
#include<math.h>
#define a 6378245
#define e2 0.006693421622966
#define PI 3.1415926535898
#define z (PI/180)
#define lambda0 116
double fun1(double phi)/*计算U值*/
{
double psi,e,U;
e=sqrt(e2);
psi=asin(e*sin(phi*z));
U=(tan((45+phi/2)*z))/pow((tan((45+psi/2)*z)),e);
return U;
}
double fun2(double phi)/*计算r值*/
{
double r;
r=(a*cos(phi*z))/sqrt(1-e2*pow(sin(phi*z),2));
return r;
}
void fun(double phi,double lambda,double alpha,double K,double rhos,double *t)/*计算坐标*/
{
double delta,rho;
delta=alpha*(lambda-lambda0);
rho=K/pow(fun1(phi),alpha);
t[0]=rhos-rho*cos(delta*z);
t[1]=rho*sin(delta*z);
}
void main()
{
double phi1,phi2,phi,phin=0,phis=0,alpha,U1,U2,r1,r2,lambda1=0,lambda2=0,lambda=0,K,rhos;
long int M0;
double t[2]={0,0};
printf("请输入标准纬线φ1,φ2\t");
/*phi1=37.5;phi2=41;lambda1=112;lambda2=120;phin=36;phis=43;M0=5000000;*/
scanf("%lf%lf",&phi1,&phi2);
printf("请输入经线区间λw、λz(角度值,以度为单位)\t");
scanf("%lf%lf",&lambda1,&lambda2);
printf("请输入纬线区间φn、φs");
scanf("%lf%lf",&phin,&phis);
printf("请输入比例尺分母,如1比1000输入:1000\t");
scanf("%ld",&M0);
U1=fun1(phi1);
U2=fun1(phi2);
r1=fun2(phi1);
r2=fun2(phi2);
alpha=(log10(r1)-log10(r2))/(log10(U2)-log10(U1));
K=(r1*pow(U1,alpha)*100)/(alpha*M0);
printf("\t经纬度\t\t\t地图上坐标\n");
rhos=K/pow(fun1(phin),alpha);
for(phi=phin;phi<=phis;phi++)
for(lambda=lambda1;lambda<=lambda2;lambda++)
{fun(phi,lambda,alpha,K,rhos,t);
printf("%lf,%lf\t\t%lf,%lf\n",phi,lambda,t[0],t[1]);}
}