Memahami istilah dalam rumus Panjang Gelar?

Kalkulator online seperti http://www.csgnetwork.com/degreelenllavcalc.html (lihat sumber halaman) menggunakan rumus di bawah ini untuk mendapatkan meter per derajat. Saya mengerti secara umum bagaimana jarak per derajat bervariasi tergantung pada lokasi lintang, tetapi saya tidak mengerti bagaimana itu diterjemahkan ke bawah. Lebih khusus lagi, dari mana konstanta, 3 "cos" istilah dalam setiap rumus, dan koefisien (2, 4, 6; 3, dan 5) untuk "lat" berasal?

    // Set up "Constants"
    m1 = 111132.92;     // latitude calculation term 1
    m2 = -559.82;       // latitude calculation term 2
    m3 = 1.175;         // latitude calculation term 3
    m4 = -0.0023;       // latitude calculation term 4
    p1 = 111412.84;     // longitude calculation term 1
    p2 = -93.5;         // longitude calculation term 2
    p3 = 0.118;         // longitude calculation term 3

    // Calculate the length of a degree of latitude and longitude in meters
    latlen = m1 + (m2 * Math.cos(2 * lat)) + (m3 * Math.cos(4 * lat)) +
            (m4 * Math.cos(6 * lat));
    longlen = (p1 * Math.cos(lat)) + (p2 * Math.cos(3 * lat)) +
                (p3 * Math.cos(5 * lat));

Jari-jari utama dari spheroid WGS84 adalah a = 6378137 meter dan perataan terbaliknya adalah f = 298.257223563, di mana eksentrisitas kuadrat adalah

e2 = (2 - 1/f)/f = 0.0066943799901413165.

Jari-jari lengkung melidional pada garis lintang phi adalah

M = a(1 - e2) / (1 - e2 sin(phi)^2)^(3/2)

dan jari-jari kelengkungan di sepanjang paralel adalah

N = a / (1 - e2 sin(phi)^2)^(1/2)

Selanjutnya, jari-jari paralelnya adalah

r = N cos(phi)

Ini adalah koreksi berganda untuk nilai-nilai bola M dan N , yang keduanya sama dengan jari - jari bola a , yang mereduksi menjadi ketika e2 = 0.

Pada titik kuning pada 45 derajat lintang utara, piringan biru jari-jari M adalah lingkaran osculating ("berciuman") ke arah meridian dan piringan merah jari-jari N adalah lingkaran terukur dalam arah paralel: keduanya cakram berisi arah "turun" pada titik ini. Angka ini melebih-lebihkan perataan bumi dengan dua urutan besarnya.

Jari-jari kelengkungan menentukan panjang derajat: ketika sebuah lingkaran memiliki jari-jari R , perimeter panjangnya 2 pi R mencakup 360 derajat, di mana panjang satu derajat adalah pi * R / 180. Mengganti M dan r untuk R - - yaitu, mengalikan M dan r dengan pi / 180 - memberikan rumus tepat sederhana untuk panjang derajat.

Formula-formula ini - yang semata-mata didasarkan pada nilai a dan f (yang dapat ditemukan di banyak tempat ) dan deskripsi spheroid sebagai ellipsoid rotasi - setuju dengan perhitungan dalam pertanyaan hingga dalam 0,6 bagian per juta (beberapa sentimeter), yang kira-kira sama besarnya besarnya koefisien terkecil dalam pertanyaan, menunjukkan mereka setuju. (Perkiraan selalu sedikit rendah.) Dalam plot kesalahan relatif dalam panjang derajat lintang adalah hitam dan garis bujur putus-putus merah:

Dengan demikian, kita dapat memahami perhitungan dalam pertanyaan sebagai perkiraan (melalui deret trigonometri terpotong) dengan rumus yang diberikan di atas.

Koefisien dapat dihitung dari deret Fourier cosinus untuk M dan r sebagai fungsi garis lintang. Mereka diberikan dalam hal fungsi elips e2, yang akan terlalu berantakan untuk mereproduksi di sini. Untuk spheroid WGS84, perhitungan saya memberi

  m1 = 111132.95255
  m2 = -559.84957
  m3 = 1.17514
  m4 = -0.00230
  p1 = 111412.87733
  p2 = -93.50412
  p3 = 0.11774
  p4 = -0.000165

(Anda mungkin menebak bagaimana p4memasukkan rumus. :) Kedekatan nilai-nilai ini dengan parameter dalam kode membuktikan kebenaran interpretasi ini. Peningkatan perkiraan ini akurat hingga jauh lebih baik daripada satu bagian per miliar di mana-mana.

Untuk menguji jawaban ini saya mengeksekusi Rkode untuk melakukan kedua perhitungan:

#
# Radii of meridians and parallels on a spheroid.  Defaults to WGS84 meters.
# Input is latitude (in degrees).
#
radii <- function(phi, a=6378137, e2=0.0066943799901413165) {
  u <- 1 - e2 * sin(phi)^2
  return(cbind(M=(1-e2)/u, r=cos(phi)) * (a / sqrt(u))) 
}
#
# Approximate calculation.  Same interface (but no options).
#
m.per.deg <- function(lat) {
  m1 = 111132.92;     # latitude calculation term 1
  m2 = -559.82;       # latitude calculation term 2
  m3 = 1.175;         # latitude calculation term 3
  m4 = -0.0023;       # latitude calculation term 4
  p1 = 111412.84;     # longitude calculation term 1
  p2 = -93.5;         # longitude calculation term 2
  p3 = 0.118;         # longitude calculation term 3

  latlen = m1 + m2 * cos(2 * lat) + m3 * cos(4 * lat) + m4 * cos(6 * lat);
  longlen = p1 * cos(lat) + p2 * cos(3 * lat) + p3 * cos(5 * lat);
  return(cbind(M.approx=latlen, r.approx=longlen))
}
#
# Compute the error of the approximation `m.per.deg` compared to the 
# correct formula and plot it as a function of latitude.
#
phi <- pi / 180 * seq(0, 90, 10)
names(phi) <- phi * 180 / pi
matplot(phi * 180 / pi, 10^6 * ((m.per.deg(phi) - radii(phi) * pi / 180) / 
       (radii(phi) * pi / 180)),
        xlab="Latitude (degrees)", ylab="Relative error * 10^6",lwd=2, type="l")

Perhitungan yang tepat dengan radiidapat digunakan untuk mencetak tabel dengan panjang derajat, seperti pada

zapsmall(radii(phi) * pi / 180)

Outputnya dalam meter dan terlihat seperti ini (dengan beberapa garis dihapus):

          M         r
0  110574.3 111319.49
10 110607.8 109639.36
20 110704.3 104647.09
...
80 111659.9  19393.49
90 111694.0      0.00

Referensi

LM Bugayevskiy dan JP Snyder, Proyeksi Peta - Manual Referensi. Taylor & Francis, 1995. (Lampiran 2 dan Lampiran 4)

JP Snyder, Proyeksi Peta - Manual Kerja. Kertas Profesional USGS 1395, 1987. (Bab 3)

Itulah rumus Haversine , meskipun diungkapkan dengan cara yang aneh.