大地測量學高斯投影正反算
阿新 • • 發佈:2018-11-26
namespace 大地測量學2
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
private void button1_Click(object sender, EventArgs e)
{
double ρ = 206264.806247096355,
a = 6378245, b = 6356863.0187730473 , c = 6399698.9017827110, α = 1 / 298.3, e2 = 0.006693421622966, e3 = 0.006738525414683;//先BLA是角度,BmLmAm是弧度直接進行運算,可以直接進行三角函式運算,感覺在最上邊宣告MATH後就不用再一個一個的代用
double B = 51.0 + 38.0 / 60.0 + 43.9023 / 3600.0, L = 111.0 + 2.0 / 60.0 + 13.136 / 3600.0;//角度
B = B / 180 * Math.PI;
L = L / 180 * Math.PI;//弧度
double W = Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B));
double V = Math.Sqrt(1 +e3 * Math.Cos(B) * Math.Cos(B));
// double N = a / W;
double M = a * (1 - e2) / (W * W * W);
double t = Math.Tan(B);
double η = e3 * Math.Cos(B);
B = 51.0 + 38.0 / 60.0 + 43.9023 / 3600.0;
L = 111 + 2.0 / 60.0+ 13.136 / 3600;//重新化為角度 B= 51.012195083333332
int n6,n3;
n6 = (int)(L / 6) + 1;//投影帶編號19
n3 = (int)(L / 3 + 0.5);
double L06 = 6 * n6 - 3;//L06即使這樣存在,是否有需要強制轉換 此時為111
double l2 = L - L06;//此時是角度
double zc = 2 * Math.PI * b + 4 * (a - b);//橢圓的周長40026876.244215891
double X6=zc/60/2*(2*(n6-1)+1);//12341620.1752999
B = B / 180 * Math.PI;
L = L / 180 * Math.PI;//弧度
l2 = l2 / 180 * Math.PI;//沒卵用
//double x = X6 + N / 2 / ρ / ρ * Math.Sin(B) * Math.Cos(B) * l2 * l2 + N / 24 / ρ / ρ / ρ / ρ * Math.Sin(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (5 - t * t + 9 * η * η + 4 * η * η * η * η)*l2*l2
// * l2 * l2 + N / 720 / ρ / ρ / ρ / ρ / ρ / ρ * Math.Sin(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (61 - 58 * t * t + t * t * t * t) * l2 * l2 * l2 * l2 * l2 * l2;
//double y = N / ρ * Math.Cos(B) * l2 + N / 6 / ρ / ρ / ρ * Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (1 - t * t + η * η) * l2 * l2 * l2 + N / 120 / ρ / ρ / ρ / ρ / ρ * Math.Cos(B) * Math.Cos(B)
//* Math.Cos(B) * Math.Cos(B) * Math.Cos(B) * (5 - 18 * t * t + t * t * t * t + 14 * η * η - 58 * η * η * t * t) * l2 * l2 * l2 * l2 * l2;
double N=6399698.902-(21562.267-(108.973-0.62*Math.Cos(B) * Math.Cos(B)) * Math.Cos(B)* Math.Cos(B)) * Math.Cos(B)* Math.Cos(B) ;
double B2=(51.0 + 38.0 / 60.0 + 43.9023 / 3600.0)*3600;//這個為秒的格式
double a0 = (32140.404 - (135.3302 - (0.7092 - 0.004 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B))/*a0*/;
double a4 = (0.25 + 0.00252 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B) - 0.04166;
double a6 = (0.166*Math.Cos(B) * Math.Cos(B) - 0.084) * Math.Cos(B) * Math.Cos(B);
double a3 = (0.3333333 + 0.001123 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B) - 0.1666667;
double a5 = 0.0083 - (0.1667 - (0.1968 + 0.004 * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B)) * Math.Cos(B) * Math.Cos(B);
double L0 = 6 * 19 - 3;//已知帶號為19
double l = L + L0;//角度
l=l*3600/ρ;
double x=6367558.4969*B2/ρ-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*Math.Sin (B)*Math .Cos (B);
double y = (1 + (a3 + a5 * l * l) * l * l) * N * l * Math.Cos(B);
Console.WriteLine("{0} {1} ", x, y);
textBox1.Text = Convert.ToString(x);
textBox2.Text = Convert.ToString(y);
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = "";
textBox2.Text = "";
}
private void groupBox1_Enter(object sender, EventArgs e)
{
}
private void button3_Click(object sender, EventArgs e)
{
double ρ = 206264.806247096355,
a = 6378245, b = 6356863.0187730473, c = 6399698.9017827110, α = 1 / 298.3, e2 = 0.006693421622966, e3 = 0.006738525414683;//先BLA是角度,BmLmAm是弧度直接進行運算,可以直接進行三角函式運算,感覺在最上邊宣告MATH後就不用再一個一個的代用
//計運算元午面的弧長
double Round = 2 * Math.PI * b + 4 * (a - b);
double X=5724004.82211011;
double Y=502559.918404555;
//以度°為單位的bf
double Bf = 27.11115372595 + 9.02468257083 * (X / 1000000 - 3) - 0.00579740442 * (X / 1000000 - 3) * (X / 1000000 - 3) - 0.00043532572 * (X / 1000000 - 3) * (X / 1000000 - 3) * (X / 1000000 - 3) + 0.00004857285 * (X / 1000000 - 3) * (X / 1000000 - 3) * (X / 1000000 - 3) * (X / 1000000 - 3);
double W = Math.Sqrt(1 - e2 * Math.Sin(Bf) * Math.Sin(Bf));
double V = Math.Sqrt(1 + e3 * Math.Cos(Bf) * Math.Cos(Bf));
double N = a / W;
double M = a * (1 - e2) / (W * W * W);
double t = Math.Tan(Bf);
double η = e3 * Math.Cos(Bf);
double B = Bf - t / 2 / M / N * Y * Y + t / 24 / M / N / N / N * (5 + 3 * t * t + η * η - 9 * η * η * t * t) * Y * Y * Y * Y;
Bf=Bf/180*Math .PI;
double l = 1 / N / Math.Cos(Bf) * Y - 1 / 6 / N / N / N / Math.Cos(Bf) * (1 + 2 * t * t + η * η) * Y * Y * Y + 1 / 120 / N / N / N / N / N / Math.Cos(Bf) * (5 + 28 * t * t + 24 * t * t * t * t) * Y * Y * Y * Y * Y;
double L0 = 6 * 19 - 3;//已知帶號為19
double y = Y - N - 500000;
double L = L0 + l;//角度
double D=(int)(L);
double F=(int)((L-(int)(L))*60);
double MM=(L-D-F/60)*3600;
textBox3.Text = Convert.ToString(B);
textBox4.Text = Convert.ToString(L);
Console.WriteLine("{0} {1} {2}", D, F, MM);//111° 7分 35.6402392184452秒
}
private void button4_Click(object sender, EventArgs e)
{
textBox3.Text = "";
textBox4.Text = "";
}
}
}
實驗截圖如上,在這裡為了簡單隻用了一個固定資料,介面好不人性化。。。需要的讀者可以自己改換為相應的資料,需要完整程式的可以聯絡我。
感謝閱讀