#include "iostream.h"
#include "math.h"
#include "stdio.h"
#define p 206264.806247   //一弧度對應的秒值 
#define pi 3.1415926  //圓周率 
#define a_  298.3     //北京54橢球的扁率倒數 
#define a  6378245.00000   //北京54橢球的長半軸 


main()
{     double B,L,x,y,X1;
      int i;
      double m0,m2,m4,m6,m8,q0,q2,q4,q6,q8,b,e2,t,g2,N,l,x1,x2,x3,y1,y2,y3;
      double X,Y,Z;
	  double p1,r,A1,A2,A3,A4;
	 
     //輸入空間直角坐標系中的坐標 
  	cout<<"請輸入已知點在空中直角坐標系中的坐標："<<endl;
	cin>>X>>Y>>Z;
	cout<<endl;
	
	b=a*(1-1/a_);
	e2=1-b*b/(a*a);
	
	//坐標轉換的直接解法求解經緯度 
	p1=atan(Z/sqrt(X*X+Y*Y));
	r=a*sqrt(1-e2)/((1-e2*sin(p1)*sin(p1))*(1-e2*sin(p1)*sin(p1)));
    r=sqrt(X*X+Y*Y+Z*Z);
	A1=a*tan(p1)/r;
    A2=sin(p1)*sin(p1)+2*(a/r)*cos(p1)*cos(p1);
    A3=3*sin(p1)*sin(p1)*sin(p1)*sin(p1)+16*((a/r)*sin(p1)*sin(p1)*cos(p1)*cos(p1)+4*(a/r)*(a/r)*cos(p1)*cos(p1)*(2-5*sin(p1)*sin(p1)));
	A4=5*sin(p1)*sin(p1)*sin(p1)*sin(p1)*sin(p1)*sin(p1)+48*((a/r)*sin(p1)*sin(p1)*sin(p1)*sin(p1)*cos(p1)*cos(p1)+20*(a/r)*(a/r)*sin(p1)*sin(p1)*cos(p1)*cos(p1)*(4-7*sin(p1)*sin(p1))+16*(a/r)*(a/r)*(a/r)*cos(p1)*cos(p1)*(1-7*sin(p1)*sin(p1)+8*sin(p1)*sin(p1)*sin(p1)*sin(p1)));
	
	B=atan(tan(p1)+A1*e2*(1+e2*A2/2*(1+e2*A3/4*(1+A4*e2/2))));

	L=atan(Y/X);
	if(L<0)
    	L=pi+atan(Y/X);
	  
     
     L=L*180/pi;
       

     //計算高斯投影帶的帶號 
	  for(i=0;fabs(3*i-L)>1.5;i++);
      printf("高斯投影帶的帶號為：");
      cout<<i<<endl;
      
      //計算經差，并用弧度的形式表示 
      l=(L-3*i)*pi/180;
      
      t=tan(B);
      g2=e2*cos(B)*cos(B)/(1-e2);
      N=a/sqrt(1-e2*sin(B)*sin(B));
	
     //計算子午弧長 X1
      m0=a*(1-e2);
      m2=3*e2*m0/2;
      m4=5*e2*m2/4;
      m6=7*e2*m4/6;
      m8=9*e2*m6/8;
      
      q0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
      q2=m2/2+m4/2+15*m6/32+7*m8/16;
      q4=m4/8+3*m6/16+7*m8/32;
      q6=m6/32+m8/16;
      q8=m8/128;
      
      X1=q0*B-q2*sin(2*B)/2+q4*sin(4*B)/4-q6*sin(6*B)/6+q8*sin(8*B)/8;

     
     //使用實用電算公式計算進行高斯投影正算
      
      x1=N*t*cos(B)*cos(B)*l*l/2;
      x2=N*t*cos(B)*cos(B)*cos(B)*cos(B)*l*l*l*l*(5-t*t+9*g2+4*g2*g2)/24;
      x3=N*t*cos(B)*cos(B)*cos(B)*cos(B)*cos(B)*cos(B)*l*l*l*l*l*l*(61-58*t*t+t*t*t*t)/720;

        
      y1=N*cos(B)*l;
      y2=N*cos(B)*cos(B)*cos(B)*l*l*l*(1-t*t+g2)/6;
      y3=N*cos(B)*cos(B)*cos(B)*cos(B)*cos(B)*l*l*l*l*l*(5-18*t*t+t*t*t*t+14*g2-58*g2*t*t)/120;

      x=X1+x1+x2+x3;
      y=y1+y2+y3;
      
      printf("對應的高斯投影坐標系中的坐標為：\n");
      printf("x=%10.6f,y=%10.6f",x,y);
      cin>>x;
      
}