计算椭球上大地线长度

CSharpBoy · 发表于 2016-10-18 20:47:00

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？立即注册

×

#include "stdio.h"
#include "9simp.c"
#include "math.h"
#define rou 206264.8062
double a=6378245;
double b=6356863.0;
#define GS 21
main()
{ static double /*要输入弧度*/
B[GS]={ },
L[GS]={ };
double B0=0.629026214,a0,b0,eps,sum=0.0,t[GS-1],simpf(double);
double L_L[GS-1],B_B[GS-1],Bm[GS-1];
double c,E,e,e1,ec[GS-1],es[GS-1],Vm[GS-1],Nm[GS-1],tm[GS-1];
double u0,v0,u1,v1,aa,bb,suml=0.0;
int i,j=0;
/*求大地多边形面积*/
a0=B0; eps=0.000001;
for(i=1;i<GS;i++)
{b0=(B+B[i-1])/2;
t[i-1]=b*b*(L-L[i-1])*simp(a0,b0,eps,simpf);
}
printf("\n与x轴面积为 ");
for(i=0;i<GS-1;i++)
sum=sum+t;
printf("%.4f ",sum);
/* 计算两点间大地线长度的程序*/
for(i=1;i<GS;i++)
{L_L[i-1]=(L-L[i-1])*rou;
B_B[i-1]=(B-B[i-1])*rou;
Bm[i-1]=(B+B[i-1])/2;
}
c=a*a/b;E=sqrt(a*a-b*b);e=E/a;e1=E/b;
for(i=0;i<GS-1;i++)
{
ec=e1*cos(Bm);
es=e*sin(Bm);
Vm=sqrt(1+ec*ec);
Nm=a/sqrt(1-es*es);
tm=tan(Bm);
u1=(L_L/rou)*Nm*cos(Bm);
v1=(B_B/rou)*Nm/(Vm*Vm);
do
{j=j+1;
u0=u1;v0=v1;
bb=(u0*u0*(2+3*tm*tm+2*ec*ec)+3*v0*v0*ec*ec*(-1+tm*tm-ec*ec-4*tm*tm*ec*ec))/(24*Nm*Nm);
aa=(u0*u0*tm*tm-v0*v0*(1+ec*ec-9*tm*tm*ec*ec))/(24*Nm*Nm);
u1=L_L*Nm*cos(Bm)/(rou*(1+bb));
v1=B_B*Nm/(rou*Vm*Vm*(1+aa));
}while((fabs(u1-u0)>=1e-10)&&(fabs(v1-v0)>=1e-10));
suml=suml+sqrt(u1*u1+v1*v1);
}
/*printf("运行次数j=%d \n",j);*/
printf("大地线长度 %.4f \n\n",suml);
}
double simpf(x)
double x;
{ double y,e;
e=sqrt(a*a-b*b)/a;
y=cos(x)/(1+e*e*sin(x)*sin(x))*(1+e*e*sin(x)*sin(x));
return(y);
}

#include "stdio.h" 
  #include "9simp.c" 
  #include "math.h" 
  #define rou 206264.8062 
  double  a=6378245; 
  double  b=6356863.0; 
  #define  GS  21 
  main() 
  { static double /*要输入弧度*/ 
   
  B[GS]={ }, 
  L[GS]={ }; 
   
  double B0=0.629026214,a0,b0,eps,sum=0.0,t[GS-1],simpf(double); 
  double L_L[GS-1],B_B[GS-1],Bm[GS-1]; 
  double c,E,e,e1,ec[GS-1],es[GS-1],Vm[GS-1],Nm[GS-1],tm[GS-1]; 
  double u0,v0,u1,v1,aa,bb,suml=0.0; 
  int i,j=0; 
 /*求大地多边形面积*/ 
  a0=B0;  eps=0.000001; 
  for(i=1;i<GS;i++)          
  {b0=(B+B[i-1])/2; 
    t[i-1]=b*b*(L-L[i-1])*simp(a0,b0,eps,simpf); 
    
  } 
  printf("\n与x轴面积为  "); 
  for(i=0;i<GS-1;i++) 
   sum=sum+t; 
  printf("%.4f  ",sum); 
 
  /* 计算两点间大地线长度的程序*/ 
  for(i=1;i<GS;i++) 
  {L_L[i-1]=(L-L[i-1])*rou; 
    B_B[i-1]=(B-B[i-1])*rou; 
    Bm[i-1]=(B+B[i-1])/2; 
  } 
 
  c=a*a/b;E=sqrt(a*a-b*b);e=E/a;e1=E/b; 
  for(i=0;i<GS-1;i++) 
  { 
  ec=e1*cos(Bm);    
  es=e*sin(Bm); 
  Vm=sqrt(1+ec*ec); 
  Nm=a/sqrt(1-es*es); 
  tm=tan(Bm); 
  u1=(L_L/rou)*Nm*cos(Bm); 
  v1=(B_B/rou)*Nm/(Vm*Vm); 
  do 
  {j=j+1; 
  u0=u1;v0=v1; 
  bb=(u0*u0*(2+3*tm*tm+2*ec*ec)+3*v0*v0*ec*ec*(-1+tm*tm-ec*ec-4*tm*tm*ec*ec))/(24*Nm*Nm); 
  aa=(u0*u0*tm*tm-v0*v0*(1+ec*ec-9*tm*tm*ec*ec))/(24*Nm*Nm); 
  u1=L_L*Nm*cos(Bm)/(rou*(1+bb)); 
  v1=B_B*Nm/(rou*Vm*Vm*(1+aa)); 
  }while((fabs(u1-u0)>=1e-10)&&(fabs(v1-v0)>=1e-10)); 
  suml=suml+sqrt(u1*u1+v1*v1); 
  } 
  /*printf("运行次数j=%d \n",j);*/ 
  printf("大地线长度 %.4f \n\n",suml); 
 
  } 
 
   
  double simpf(x) 
  double x; 
  { double y,e; 
   e=sqrt(a*a-b*b)/a; 
    y=cos(x)/(1+e*e*sin(x)*sin(x))*(1+e*e*sin(x)*sin(x)); 
    return(y); 
  }

账号		自动登录	找回密码
密码			立即注册

[代码] 计算椭球上大地线长度

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。