[地图投影之C实现]BASIC_FU.C

double Eccentricity()                             /*THE FIRST ECCENTRICITY*/
{
   double a,b,c;
   a=pow(6378245.00,2);
   b=pow(6356863.00,2);
   c=(a-b)/a;
   c=sqrt(c);
   return c;
}

#define eccentricity Eccentricity()

double Eccentricity1()                           /*THE SECOND ECCENTRICITY*/
{
   double x,y;
   x=eccentricity*eccentricity;
   y=x/(1-x);
   return sqrt(y);
}

#define eccentricity1 Eccentricity1()

double Meridian_radius(double latitude)
{
   double x,y;
   x=sin(latitude);
   x=pow((1-eccentricity*eccentricity*x*x),1.5);
   y=major_radius*(1-eccentricity*eccentricity);
   return y/x;
}

double Prim_vertical_radius(double latitude)
{
   double x;
   x=sin(latitude);
   x=pow(1-eccentricity*eccentricity*x*x,0.5);
   return major_radius/x;
}

double Parallel_radius(double latitude)
{
   double x,y;
   x=cos(latitude);
   y=Prim_vertical_radius(latitude);
   return x*y;
}

double Meridian_length(double lowest_latitude,double highest_latitude)
{
   double x,h,meridian_length=0;
   int n=10000;
   h=(highest_latitude-lowest_latitude)/n;
   for(x=lowest_latitude;x      meridian_length+=Meridian_radius(x)*h;
   return meridian_length;
}

double Parallel_length(double latitude,double lowest_longitude,double highest_longitude)
{
   return Prim_vertical_radius(latitude)*cos(latitude)*(highest_longitude-lowest_longitude);
}

double Area(double lowest_latitude,double highest_latitude,double lowest_longitude,double highest_longitude)
{
   double x=highest_longitude-lowest_longitude;
   double y=(highest_latitude+lowest_latitude)/2;
   double z=(highest_latitude-lowest_latitude)/2;
   double k=2*pow(major_radius,2)*(1-pow(eccentricity,2))*x;
   double a=1+pow(eccentricity,2)/2+3*pow(eccentricity,4)/8+5*pow(eccentricity,6)/16;
   double b=pow(eccentricity,2)/6+3*pow(eccentricity,4)/16+3*pow(eccentricity,6)/16;
   double c=3*pow(eccentricity,4)/80+pow(eccentricity,6)/16;
   double d=pow(eccentricity,6)/112;
   double area=0;
   area=k*(a*sin(z)*cos(y)-b*sin(3*z)*cos(3*y)+c*sin(5*z)*cos(5*y)-d*sin(7*z)*cos(7*y));
   return area;
}

double Radiam(double dms)                         /*CHANGE DGREE TO ARC*/
{
   double p=3.14159265358979/180;
   return dms*p;
}

double U(double latitude)                  /*IMPORTANT CONST BEGIN FROM EX2*/
{
   double a=Radiam(45),e=Eccentricity(),u;
   u=asin(e*sin(latitude));
   u=tan(a+latitude/2)/pow(tan(a+u/2),e);
   return u;
}

double U1(double latitude)            /*FOR GAMA IN MAIN OF EX4*/
{
   double x;
   x=eccentricity1*eccentricity1*cos(latitude)*cos(latitude);
   return sqrt(x);
}

double Radiam1(double x,double y,double z)          /*CHANGE D.M.S. TO ARC*/
{
   double m;
   m=x+y/60+z/3600;
   return Radiam(m);
}