1

在此处输入图像描述

在此处输入图像描述

这是我的代码：

public class Test {

  public static void main(String[] args) {
    // TODO Auto-generated method stub
    System.out.println("Test ");
    gaussPositiveFormula(36.082846593,119.5520122887);
  }

  /* gauss-kruger projection forward calculation
   * params:B latitude，L longitude
   */
  public static void gaussPositiveFormula (double B,double L){

    double L0 = (L%15)>7.5?Math.floor(L/15+1)*15:Math.floor(L/15)*15;    // longitude of central meridian
    double l =  L - L0;

    int [] array_B = DmsRad.dec2Dms(B);
    B = DmsRad.dms2rad(array_B[0], array_B[1], array_B[2]);
    int [] array_L = DmsRad.dec2Dms(L);
    L = DmsRad.dms2rad(array_L[0], array_L[1], array_L[2]);
    int [] array_L0 = DmsRad.dec2Dms(L0);
    L0 = DmsRad.dms2rad(array_L0[0], array_L0[1], array_L0[2]);



    System.out.println("B = "+ B);
    System.out.println("L = " + L);

    //WGS-84 spheroid params
    double a = 6378137.0;   //major semi axis
    double efang = 0.0066943799013;     //square of e
    double e2fang = 0.00673949674227;   //suqre of e2
    double c = 6399593.6258;
    double sinB = Math.sin(B);
    double cosB = Math.cos(B);
    double t = Math.tan(B);
    double p = 180*3600/Math.PI;
    double nfang = e2fang*cosB*cosB;



    double m = l*cosB;
    double V = Math.pow((1+e2fang*cosB*cosB), 1/2.0);
    double N = a*Math.pow((1-efang*sinB*sinB),-1/2);   
    N = c/V;
    System.out.println("L0 = "+ L0);
    System.out.println("sinB = "+ sinB);
    System.out.println("cosB = "+ cosB);
    System.out.println("t = "+ t);
    System.out.println("p = "+ p);
    System.out.println("nfang = "+ nfang);
    System.out.println("l = "+ l);

    System.out.println("N = "+ N);

    //principal curvatures
    double m0,m2,m4,m6,m8;
    m0 = a*(1-efang);
    m2 = 3.0/2.0*efang*m0;
    m4 = efang*m2*5.0/4.0;
    m6 = efang*m4*7.0/6.0;
    m8 = efang*m6*9.0/8.0;

    //meridional curvature
    double a0,a2,a4,a6,a8;
    a0 = m0 + m2/2.0 + m4*3.0/8.0 + m6*5.0/16.0 + m8*35.0/128.0; 
    a2 = m2/2.0 + m4/2.0 + m6*15.0/32.0 + m8*7.0/16.0;
    a4 = m4/8.0 + m6*3.0/16.0 + m8*7.0/32.0;
    a6 = m6/32.0 + m8/16.0;
    a8 = m8/128.0;


    //meridional arc length calculate

    double X0 = a0*B -a2*Math.sin(2*B)/2.0 + a4*Math.sin(4*B)/4.0 - a6*Math.sin(6*B)/6.0 + a8*Math.sin(8*B)/8.0;
    System.out.println("X0 = "+ X0);

    double x = X0 + N*sinB*cosB*l*l/(p*p*2.0) + N*sinB*Math.pow(cosB, 3)*(5-t*t+9*nfang+4*nfang*nfang)*Math.pow(l, 4)/(24*Math.pow(p, 4)) + N*sinB*Math.pow(cosB, 5)*(61-58*t*t + Math.pow(t, 4))*Math.pow(l, 6)/(720*Math.pow(p, 6)); 
    System.out.println("x = "+ x);
    double y = N*cosB*l/p  +  N*Math.pow(cosB, 3)*(1-t*t + nfang)*Math.pow(l, 3)/(6*Math.pow(p, 3)) + N*Math.pow(cosB, 5)*(5-18*t*t +Math.pow(t, 4) + 14*nfang -58*nfang*t*t)*Math.pow(l, 5)/(120*Math.pow(p, 5));



    System.out.println("y = "+ y);
    System.out.println(N*cosB*l/p );

  }

  //transform coordinates to degrees
  public static double formToDegress (double x){

    double x1 = Math.abs(x);

    double degree,min,sec;

    degree = Math.floor(x1);
    x1 = x1-degree;
    min =  Math.floor(x1*60);
    x1 = (x1*60-min);
    sec = x1*60;

    double s = degree*3600+min*60+sec;

    return s;
  }

}

这是我的代码，计算结果相当正确。

我发现参数l是主要因素。例如坐标是(36.0,119.0),L = 119.0，我已经计算了L0 = 120.0，对吗？并计算l = L -L0，我应该转换l度数吗？

m0~m8是常数，是由WGS84参数计算出来的，我觉得是对的。a0~a8使用方法相同。

我的问题是：如何使用高斯公式？公式的参数，需要我转换哪种格式？</p>

4

0 回答 0