高斯正算结果不对，求查错

边走边创

void CmapgisDoc::Calculateellipse2plane(double B,double L,double&x,double&y,int index,int tuoqiu)//如使用6°带index请输入1，使用3°带请输入2
{
int ProjNo=0;
int ZoneWide; ////带宽
double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
double a,f, e2,ee, NN, T,C,A, M, iPI;
iPI =  0.0174532925199433; ////3.1415926535898/180.0;  
if (tuoqiu==0)
{
  a=6378245.0;
  f=1.0/298.3; //54年北京坐标系参数
}
else
{
  a=6378140.0;
  f=1/298.257; //80年西安坐标系参数
}
if(index==1)   //已知a点在6°带的带号和中央子午线经度
{
  ProjNo =int(L/6);
  longitude0=6*ProjNo+3;
  ZoneWide=6; ////6度带宽
}
if(index==0)  //已知a点在3°带的带号和中央子午线经度
{
  ProjNo=int((L+1.5)/3);
  longitude0=3*ProjNo;
  ZoneWide=3; ////3度带宽
}
longitude0 =longitude0 * iPI ;
latitude0=0;
longitude1 = L * iPI ; //经度转换为弧度
latitude1 =B * iPI ; //纬度转换为弧度
e2=2*f-f*f;
ee=e2*(1.0-e2);
NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));
T=tan(latitude1)*tan(latitude1);
C=ee*cos(latitude1)*cos(latitude1);
A=(longitude1-longitude0)*cos(latitude1);
M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2*e2
  *e2/1024)*sin(2*latitude1)+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-
  (35*e2*e2*e2/3072)*sin(6*latitude1));
xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);
yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
  +(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);
X0 = 1000000L*(ProjNo+1)+500000L;
Y0 = 0;
xval = xval+X0;
yval = yval+Y0;
y=xval;
x=yval;
}

NN卯酉圈曲率半径，测量学里面用N表示

M为子午线弧长，测量学里用大X表示

fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示

R为底点所对的曲率半径，测量学里用Nf表示