GPS簡介講義

2007年2月28日 10:18

一份為簡介GPS所做的講義.

1. GPS簡介

2. 內容

3. 歷史

4. 系統架構

5. 太空衛星

6. 地面控制

7. 使用者

8. 定位原理

9. 平面三點定位

10. 空間多點定位

11. 定位識別

12. 延遲判斷

13. 都卜勒效應

14. 座標系統

15. 誤差來源

16. DGPS與WAAS

17. 其他定位系統

18. 應用與發展
(繼續閱讀)

::意見(1)

coLinux 介紹與應用

2005年6月29日 07:12

1. coLinux 簡介

Cooperative Linux, 簡稱 coLinux, 是一種對 Linux kernel 的移植, 讓一台機器可以協同運作不同的作業系統. coLinux 的前身 UMLWin32 最早是由 Dan Aloni 在 2000 年所開發, 當時的目的是為了將 User Mode Linux 移植到 Cygwin 上. 在 2003 年時, Dan Aloni 運用了不同以往的想法與做法, 於是, 便產生了 coLinux.

coLinux 不同於 VMware 等模擬器, coLinux 本身並不是模擬一台PC, 而是透過它本身的特殊設計, 讓在其中運作的 Linux kernel 直接使用 Windows 的硬體資源. 因此, 對大部份的使用者來說, coLinux 本身的速度與真實機器的速度相差無幾.

更重要的是, coLinux 是完全免費的. (GPL2)



2. coLinux 系統需求
  • 在 Intel 386 以上相容機器運作
  • 在 Windows 2000, Windows XP, Windows 2003 下執行
  • 最低需求: Windows 2000, PIII 450Mhz, 64MB RAM

3. coLinux 功能

A. 現有功能

  • 支援 Linux 2.4.x, Linux 2.6.x
  • 支援各種 Linux Distribution (Debian/Gentoo/Fedora/TopologiLinux/Mandrake/Slackware/Knoppix)
  • 可存取 Windows 與 Linux Partition
  • 可存取 Windows 檔案及目錄
  • 可存取 CD/DVD
  • 可存取 RAID/LVM
  • 支援網路 (透過 TAP 或 Bridge)
  • 支援聲音 (透過 ESD)

B. 未來功能

  • 支援 Frame Buffer
  • 支援 Native X Server

4. 安裝 coLinux

coLinux 目前最新的版本是 0.6.2 版, 可以到 coLinux Download 區選取下載或者直接下載 coLinux 0.6.2.

安裝過程如下:

  • 安裝 coLinux 0.6.2

  • coLinux 是採用 GPL 授權

  • 安裝時不需下載 Root Filesystem image

  • 安裝位置通常都是用 C:\coLinux

  • 如果需要使用 Bridged Ethernet, 需自行下載安裝 WinPCAP

  • coLinux 的安裝過程很簡單

  • coLinux 附的 TAP Driver, 如果要使用 TAP 的 Network, 就一定要安裝

  • 安裝完成

5. coLinux 基本設定

coLinux 的所有基本設定, 都可以在 xml 設定檔中定義. 底下以 coLinux 預附的 default.colinux.xml 來說明設定檔中, 各項參數的意義.

  • <block_device enabled="true" path="\DosDevices\c:\coLinuxroot_fs" index="0" />
  • <block_device enabled="true" path="\DosDevices\c:\coLinuxswap_device" index="1" />

coLinux 是使用它特有的 device (cobd), 這裡是在定義 /dev/cobd0 與 /dev/cobd1. cobd0 指向 c:\coLinux\root_fs, cobd1 則指向 c:\coLinux\swap_device. 在這裡, 這兩個 device 都是 image file. coLinux 下預設的 boot image 是 cobd0, 而預設的 swap 則是 cobd1.

  • root=/dev/cobd0

這是用來設定 linux kernel 開機時所需傳遞的參數. 參數設定的方法與 grub 及 lilo 相同.

  • <bootparams>root="/dev/cobd0"</bootparams>

這是指定用來做為 initrd 的 image file. coLinux 本身有它自己專門的 initrd.

  • <image path="vmlinux" />

這是指定用來開機的 linux kernel image. 預設是使用 coLinux 附的 kernel image.

  • <memory size="64" />

這是分配給 coLinux 的記憶體大小, 如果只有跑 linux console, 至少需 64MB, 如果有執行 X, 則最好能配置 128MB 以上.

  • <network index="0" type="tap" />

這是 coLinux 所要使用的網路型態, 有 TAP 與 Bridged 兩種.

以上是基本設定, 其他進階設定, 請參照 coLinux 官方文件.

6. 執行 coLinux

當所有的設定都沒問題時, 就可以執行 coLinux .coLinux 有兩種執行方法, 透過 Command-Line 或透過 System Service. 底下分別說明兩種方式的不同.

7. coLinux 下的網路

A. TAP

在 TAP 模式下, coLinux 是透過 TAP Driver 與 ICS(Internet Connection Sharing)的方式來使用 Windows 的網路. 這時, Windows 的 ip 預設為 192.168.0.1, 而 coLinux 的預設 ip 則通常為 192.168.0.2.

要順利使用 TAP 模式, 有幾個動作要先完成.

  • 選擇目前對外的網路卡, 允許 TAP 使用這張網路卡的共用連線

  • 檢視 TAP 網路卡的 TCP/IP

  • TAP 的 IP 應該會被設為 192.168.0.1

順利完成以上動作後, coLinux 執行後, TAP 就會正常運作.

B. Bridged

若要使用 Bridged 模式, 設定檔中網路的部份要改成

coLinux 所偵測到的網路卡則應該是 "Realtek RTL8139/810x Family Fast Ethernet NIC".

8. 在 coLinux 上運行 Fedora - 使用 Image File

如果想要在 Windows 上, 透過 coLinux 來運行一個全新的 Fedora, 使用 image file 是個不錯的方式.

首先, 必須先下載安裝 coLinux, 方法如上. 另外, 還需要下載 Fedora 的 image file. coLinux 官方附的是 Fedora Core 1 的 image, 這個 image 可以使用在 linux 2.4.x 與 linux 2.6.x 上.

另外, 如果想使用 linux swap, 可以下載預先製作好的 Swap Image.

底下是 coLinux 使用 Fedora 的設定檔範例:

<?xml version="1.0" encoding="UTF-8?>
<colinux />
<block_device enabled="true" path="\DosDevices\c:\coLinuxfc1_2GB_root" index="0" />
<block_device enabled="true" path="\DosDevices\c:\coLinuxswap_64MB" index="7" />
<bootparams />ro root=/dev/cobd0 3 fastboot nogui</bootparams />
<initrd path="initrd.gz" />
<image path="vmlinux" />
<memory size="64" />
<network index="0" type="tap" />
</colinux />

有幾點要注意.

  • Fedora 的 swap device 並不是在 cobd1, 而是在 cobd7.
  • 在 bootparams 所下的參數, 不能使用 root=/dev/cobd0, 而要使用 ro root=/dev/cobd0.
  • 第一次執行 coLinux Fedora 時, 可能會因為找不到 kernel module 而造成執行失敗. 這時再重新執行一次就可以了.
  • 當 Fedora Core 1 順利開機後, 要將 coLinux 所附的 kernel module 檔 vmlinux-modules.tar.gz 解壓縮到相對應的位置.
  • coLinux 官方所附的 Fedora Core 1, 一開始的網路設定是有問題的, 需要手動改 /etc/resolv.conf 與 /etc/sysconfig/network-scripts/ifcfg-eth0 內的設定.
  • 開機完成的 Fedora Core 1 若要升級, 請按照 Fedora 本身的方式升級即可.

9. 在 coLinux 上運行 Fedora - 使用 Native Partition

如果已經在電腦中將 Fedora 安裝在單獨的 Partition 中, 通常都會使用 Dual Boot 來做作業系統中的切換. 有時不斷地切換開機是一件很煩人的事, 而 VMware 等模擬器又無法存取 Native Fedora.

用 coLinux 就沒有這種煩惱! coLinux 可以將原本需要 Dual Boot 開機執行的 Fedora Partition, 直接在 Windows 中開機執行, 而且就像是一台完整的機器一般, 依然保有所有原來的設定.

這又是怎麼做到的呢?

首先要轉換原本 Fedora 的 Partition. 要在 /dev 建立 cobdN(0..7). 這可以使用底下的 script 做到:

  for i in 0 1 2 3 4 5 6 7   do      mknod /dev/cobd$i b 117 $i   done

然後要將會因 coLinux 與 Fedora 環境不同而受到影響的檔案做自動判斷. 總共會有以下幾個檔案

/etc/fstab
/etc/hostname
/etc/resolv.conf
/etc/sysconfig/network-scripts/ifcfg-eth0

可以將這幾個檔案依不同環境, 分別製作 colinux 與 fedora 專用. 並且利用bootparams傳入COLINUX這個環境變數. 自動判斷的部份則可以使用以下的 script:

if [[ $COLINUX ]]; then
SUFFIX=colinux
else
SUFFIX=fedora
fi

for conf_file in
"/etc/fstab"
"/etc/hostname"
"/etc/resolv.conf"
"/etc/sysconfig/network-scripts/ifcfg-eth0"
; do
cp -f $conf_file-$SUFFIX $conf_file
done

另外要記得將 coLinux 所會用到的 kernel module 裝到 Fedora 的 Partition.

接著便是 coLinux 的設定檔. 假設電腦上的硬碟是用以下的方法來分割:

那麼, 在 Fedora 上看到的裝置就分別是 hda1(C:), hda2(D:), hda5(/boot), hda6(/), hda7(swap). 但是, 在 coLinux 上的裝置則是有所不同, 分別為 partition1(C:), partition2(D:), partition3(/boot), partition4(/), partition5(swap). 關於在 Windows 上檢視 Partition, 請使用 dd for Windows.

所以, coLinux 的 Native Boot 的設定檔應該如以下的範例:

<?xml version="1.0" encoding="UTF-8"?>
<colinux>
<block_device index="5" path="\Device\Harddisk0\Partition3" enabled="true" alias=hda5 />
<block_device index="6" path="\Device\Harddisk0\Partition4" enabled="true" alias=hda6 />
<block_device index="7" path="\Device\Harddisk0\Partition5" enabled="true" alias=hda7 />
<bootparams>ro root=/dev/hda6 3 fastboot COLINUX=1</bootparams>
<initrd path="initrd.gz" />
<image path="vmlinux"></image>
<memory size="256"></memory>
<network index="0" type="tap" />
</colinux>
實際上 Native Fedora 在 coLinux 開機之後, 會像是以下的畫面:


10. 在 coLinux 上執行 X

coLinux 目前還不能直接支援 X. 要想在 coLinux 的環境下執行 X, 主要有兩種方法. 第一種是在 Windows 安裝 XFree86 Server, 再將 coLinux 上的 X 轉向到 Windows 上的 XFree86. 第二種則是透過 VNC 來執行 X. 這裡要介紹的是第二種方法.

步驟如下:

  • 首先要安裝 X Server 所需的所有套件, 並且安裝 vnc-server.
  • 將 /root/.vnc/xstartup 增加 startkde.
  • 執行 vncserver -geometry 800x600
  • 在 Windows 上執行 vnc-client, 連接到 192.168.0.2:5901

當連接成功後, 看到的是以下的畫面:


如果要在 Fedora 一開始時就自動執行 VNC, 請參考 VNC On Boot.

如果想嘗試其他在 coLinux 上執行 X 的方法, 請參考 XCoLinux.

11. 相關資源

A. Command-Line

範例: colinux-daemon.exe -c default.colinux.xml -t nt

參數:

  • -c 執行 coLinux 時使用的設定檔
  • -t Console 的型態, 有 nt 與 fltk 兩種

B. System Service

要先將 coLinux 安裝成 System Service. 方法如下:

colinux-daemon.exe -c default.colinux.xml --install-service "Serive Name"

參數:

  • -c 安裝 coLinux 系統服務時使用的設定檔
  • --install-service 要安裝的系統服務名稱

系統服務安裝成功後, 可以使用系統管理工具來啟動或停止 coLinux. 也可以使用 net start/stop 的方法來管理.

另外, 也可以使用 coLinuxManager 來管理 coLinux 的系統服務.



(繼續閱讀)

::意見(4)

A* 演算法簡介 (A* Algorithm Brief)

2004年12月20日 09:44

A* (A-Star) 演算法是在Game中通常用來解決最短路徑(Shortest Path)問題的一種演算法. 相對於另一個知名的 Dijkstra 演算法來說, Dijkstra演算法雖然可以保證找到一條最短的路徑, 但不如A* 演算法這樣簡捷快速. 同時, Dijkstra的搜尋深度在某些情形下也容易顯得不適用. A* 演算法便是為了這些情形而出現的, 可以算是 Dijkstra演算法的一種改良.

底下用幾張圖來表現Dijkstra演算法與A* 演算法的不同:






Dijkstra 演算法A* 演算法
無障礙
有障礙
(圖片取自 Amit's Thoughts on Path-Finding and A-Star)

從以上的比較圖可以看出, A* 演算法的搜尋深度雖然不如 Dijkstra演算法, 但其結果仍然是很令人滿意的. 為什麼A* 演算法可以達到這樣的結果呢? 這是因為A* 演算法採用了一套特殊的啟發式評價(Heuristic Estimate)公式將許多明顯為壞的路徑排除考慮, 進而快速計算出一條滿意的路徑.

公式如下:

f(n) = g(n) + h(n)

g(n): 從啟始點到目前節點的距離
h(n): 預測目前節點到結束點的距離(此為 A* 演算法的主要評價公式)
f(n): 目前結點的評價分數

其中, h(n) 主導著A* 演算法的表現方式. 有以下幾種情形:
1. h(n)=0: A* 演算法這時等同 Dijkstra演算法, 並且保證能找到最短路徑
2. h(n)<目前節點到結束點的距離: A* 演算法保證找到最短路徑, h(n)越小, 搜尋深度越深
3. h(n)=目前節點到結束點的距離: A* 演算法僅會尋找最佳路徑, 並且能快速找到結果
4. h(n)>目前節點到結束點的距離: 不保證能找到最短路徑, 但計算比較快
5. h(n)與g(n)高度相關: A* 演算法此時成為 BFS (Best-First Search)

A* 演算法的虛擬碼如下:

Add START to OPEN list

while OPEN not empty

get node n from OPEN that has the lowest f(n)

if n is GOAL then return path

move n to CLOSED

for each n' = CanMove(n, direction)

g(n') = g(n) + 1

calculate h(n')

if n' in OPEN list and new n' is not better, continue

if n' in CLOSED list and new n' is not better, continue

remove any n' from OPEN and CLOSED

add n as n's parent

add n' to OPEN

end for

end while

if we get to here, then there is No Solution


以上虛擬碼適用於一般遊戲中方格圖(Grid Map)的情形. 當要在真實地圖的路網(Road Map)上使用A* 演算法時, g(n)與h(n)的計算方式便會有所不同.

相關參考資料:

網頁:
Amit's Thoughts on Path-Finding and A-Star
Creating an A* algorithm Robot for QUT SoKoBan

維基百科:
Shortest Path Problem
A* Search Algorithm
Dijkstra Algorithm

書籍:
Core Techniques and Algorithms in Game Programming (英文)
大師談遊戲程式設計：核心技術與演算法 (中文)
(繼續閱讀)

::意見(4)

四色定理 (Four Color Theorem)

2004年8月23日 23:59

原本, 這叫做四色猜想(Four Color Problem), 不過現在已被證實, 成為一個定理了.

今天剛好遇上要討論一個四色灰階的LCD是否足夠表現一張完整的地圖, 因此, 就開始稍稍說明了一下四色定理.

四色猜想被認為是費馬最後定理相媲美的難題. 猜想的內容很簡單, 但是證明卻很難. 目前的證明極為複雜, 是由電腦所算出的, 仍不斷的有數學家在找尋更簡單的證明方法.



西元1852年，一位英國的業餘數學家弗朗西斯．格斯裡閒來沒事，拿起色筆替一份英國的分郡地圖著色的時候，突然異想天開：『如果要替所有想像得到的地圖著色，而且有共同邊相鄰的區域都不同色的話，最多需要幾種顏色呢？』

這個問題流傳到數學界，許多數學家深入地思考與嘗試之後，發現找得到的例子裡，都只需要四種顏色就可以了！但是，這不夠，必須找出一種嚴謹的數學證明，可以涵蓋任何地圖才行。

到了1879年，當時英國的數學家肯普提出一份論文，似乎證明了這個『四色猜想』，而大家也都以為這個問題已經解決了。沒想到十一年後的1890年，數學家希伍德找出了肯普的錯誤，推翻了他的證明。但是希伍德自己卻證明出『五色定理』，也就是說最多不會超過五種顏色！

不過後來的進展就非常緩慢了！一直到了1970年，數學家才證明出所有少於三十九個區域的地圖，『四色猜想』是對的。

但是，如果有一千個區域，要等到哪一年才能證明出來呢？於是，有人從不同的方向著手，並成功地將無限多的地圖簡化成1482種基本圖。

問題是：每種基本圖的顏色組合，就幾乎已經等於無限多了，想要以人工來驗證這一千多種基本圖，根本是不可能的！

還好，電腦的出現，解決了這個難題！在1975年，數學家利用三台當時最先進的大型電腦，總共花了一千兩百小時的計算，分析驗證了1482種基本圖之後，終於證明成功，而使『四色猜想』正式成為『四色定理』！

相關參考網址:

四色問題
四色定理
The Four Color Theorem
Four-Color Theorem
(繼續閱讀)

::意見(5)

Taiwan Datum (TWD67, TWD97, TM2 Source)

2004年6月3日 18:28
台灣座標系統轉換原始碼 Taiwan Datum Convert Source(TWD67, TWD97, TM2 Source)


Related Site:

洪朝貴的首頁 (Chinese)

積丹尼 ＝ Dan Jacobson (Chinese/English)

Dan Jacobson (English)

Source Below:




/*
TWD67/TWD97/TM2 datum convert functions

based on: GFC (http://www.samblackburn.com/gfc/)
proj (http://www.remotesensing.org/proj)

rewrite: minstrel(http://www.minstrel.idv.tw), 2004/6/1

license: GNU General Public License
*/

#include <math.h>

// from long-lat to tm2
void toTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y);

// from tm2 to long-lat
void fromTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y);

// used by toTM2 and fromTM2
double mercator(double y, double a, double ecc);

// from TWD97 to TWD67
void toTWD67(double* x, double* y, double* z);

// from TWD67 to TWD97
void toTWD97(double* x, double* y, double* z);

/*
Definition of math related value
*/

#define COS67_5 0.3826834323650897717284599840304e0
#define PI 3.14159265358979323e0
#define HALF_PI 1.570796326794896615e0
#define DEG_RAD 0.01745329251994329572e0
#define RAD_DEG 57.295779513082321031e0

/*
Definition of datum related value
*/

#define AD_C 1.0026000e0

#define TWD67_A 6378160.0e0
#define TWD67_B 6356774.7192e0
#define TWD67_ECC 0.00669454185458e0
#define TWD67_ECC2 0.00673966079586e0
#define TWD67_DX -752.32e0 // different from garmin and already knowned value, but those all value only
#define TWD67_DY -361.32e0 // got 5-15m accuracy. the real offical value is holded by somebody and not
#define TWD67_DZ -180.51e0 // release to public. if can got more enough twd67/twd97 control point coordinare,
#define TWD67_RX -0.00000117e0 // then we can calculate a better value than now.
#define TWD67_RY 0.00000184e0 //
#define TWD67_RZ 0.00000098e0 // and, also lack twd67/twd97 altitude convertion value...
#define TWD67_S 0.00002329e0 //


#define TWD97_A 6378137.0e0
#define TWD97_B 6356752.3141e0
#define TWD97_ECC 0.00669438002290e0
#define TWD97_ECC2 0.00673949677556e0

#define TWD67_TM2 0.9999e0 // TWD67->TM2 scale
#define TWD97_TM2 0.9999e0 // TWD97->TM2 scale

/*
datum convert function
*/

void toTWD97(double* x, double* y, double* z)
{
double r, pole, sin_lat, cos_lat;
double lat, lon, height;
double x1, y1, z1, x2, y2, z2;
double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;

lon = *x*DEG_RAD;
lat = *y*DEG_RAD;
height = *z*DEG_RAD;

if(lat<-HALF_PI&&lat>-1.001*HALF_PI)
lat = -HALF_PI;
else if(lat>HALF_PI&&lat<1.001*HALF_PI)
lat = HALF_PI;
else if((lat<-HALF_PI)||(lat>HALF_PI))
return;

if(lon>PI)
lon -= (2*PI);

sin_lat = sin(lat);
cos_lat = cos(lat);
r = TWD67_A/(sqrt(1.0-TWD67_ECC*sin_lat*sin_lat));
x1 = (r+height)*cos_lat*cos(lon);
y1 = (r+height)*cos_lat*sin(lon);
z1 = ((r*(1-TWD67_ECC))+height)*sin_lat;

x2 = x1+TWD67_DX+TWD67_S*(lon+TWD67_RZ*lat-TWD67_RY*height);
y2 = y1+TWD67_DY+TWD67_S*(-TWD67_RZ*lon+lat+TWD67_RX*height);
z2 = z1+TWD67_DZ+TWD67_S*(TWD67_RY*lon-TWD67_RX*lat+height);

pole = 0;

if(x2!=0.0)
{
lon=atan2(y2,x2);
}
else
{
if(y2>0)
{
lon = HALF_PI;
}
else if(y2<0)
{
lon=-HALF_PI;
}
else
{
pole=1;

lon=0;

if(z2>0)
{
lat=HALF_PI;
}
else if(z2<0)
{
lat=-HALF_PI;
}
else
{
lat=HALF_PI;

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = -TWD97_B;

return;
}
}
}

q2 = x2*x2+y2*y2;
q = sqrt(q2);
t = z2*AD_C;
s = sqrt(t*t+q2);
sin_b = t/s;
cos_b = q/s;
t1 = z2+TWD97_B*TWD97_ECC2*sin_b*sin_b*sin_b;
sum = q-TWD97_A*TWD97_ECC*cos_b*cos_b*cos_b;
s1 = sqrt(t1*t1+sum*sum);
sin_p = t1/s1;
cos_p = sum/s1;
r = TWD97_A/sqrt(1.0-TWD97_ECC*sin_p*sin_p);

if(cos_p>=COS67_5)
{
height=q/cos_p-r;
}
else if(cos_p<=-COS67_5)
{
height=q/-cos_p-r;
}
else
{
height=z2/sin_p+r*(TWD97_ECC-1.0);
}

if(!pole)
{
lat = atan(sin_p/cos_p);
}

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = height;
}

void toTWD67(double* x, double* y, double* z)
{
double r, pole, sin_lat, cos_lat;
double lat, lon, height;
double x1, y1, z1, x2, y2, z2;
double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;

lon = *x*DEG_RAD;
lat = *y*DEG_RAD;
height = *z*DEG_RAD;

if(lat<-HALF_PI&&lat>-1.001*HALF_PI)
lat = -HALF_PI;
else if(lat>HALF_PI&&lat<1.001*HALF_PI)
lat = HALF_PI;
else if((lat<-HALF_PI)||(lat>HALF_PI))
return;

if(lon>PI)
lon -= (2*PI);

sin_lat = sin(lat);
cos_lat = cos(lat);
r = TWD97_A/(sqrt(1.0-TWD97_ECC*sin_lat*sin_lat));
x1 = (r+height)*cos_lat*cos(lon);
y1 = (r+height)*cos_lat*sin(lon);
z1 = ((r*(1-TWD97_ECC))+height)*sin_lat;

x2 = x1-TWD67_DX-TWD67_S*(lon+TWD67_RZ*lat-TWD67_RY*height);
y2 = y1-TWD67_DY-TWD67_S*(-TWD67_RZ*lon+lat+TWD67_RX*height);
z2 = z1-TWD67_DZ-TWD67_S*(TWD67_RY*lon-TWD67_RX*lat+height);

pole = 0;

if(x2!=0.0)
{
lon=atan2(y2,x2);
}
else
{
if(y2>0)
{
lon = HALF_PI;
}
else if(y2<0)
{
lon=-HALF_PI;
}
else
{
pole=1;

lon=0;

if(z2>0)
{
lat=HALF_PI;
}
else if(z2<0)
{
lat=-HALF_PI;
}
else
{
lat=HALF_PI;

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = -TWD67_B;

return;
}
}
}

q2 = x2*x2+y2*y2;
q = sqrt(q2);
t = z2*AD_C;
s = sqrt(t*t+q2);
sin_b = t/s;
cos_b = q/s;
t1 = z2+TWD67_B*TWD67_ECC2*sin_b*sin_b*sin_b;
sum = q-TWD67_A*TWD67_ECC*cos_b*cos_b*cos_b;
s1 = sqrt(t1*t1+sum*sum);
sin_p = t1/s1;
cos_p = sum/s1;
r = TWD67_A/sqrt(1.0-TWD67_ECC*sin_p*sin_p);

if(cos_p>=COS67_5)
{
height=q/cos_p-r;
}
else if(cos_p<=-COS67_5)
{
height=q/-cos_p-r;
}
else
{
height=z2/sin_p+r*(TWD67_ECC-1.0);
}

if(!pole)
{
lat = atan(sin_p/cos_p);
}

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = height;
}

void toTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y)
{
double x0, y0, x1, y1, m0, m1;
double n, t, c, A;

x0 = *x*DEG_RAD;
y0 = *y*DEG_RAD;

x1 = lon*DEG_RAD;
y1 = lat*DEG_RAD;

m0 = mercator(y1, a, ecc);
m1 = mercator(y0, a, ecc);

n = a/sqrt(1-ecc*pow(sin(y0), 2.0));
t = pow(tan(y0), 2.0);
c = ecc2*pow(cos(y0), 2.0);
A = (x0-x1)*cos(y0);

*x = scale*n*(A + (1.0 - t + c)*A*A*A/6.0 + (5.0 - 18.0*t + t*t + 72.0*c - 58.0*ecc2)*pow(A, 5.0)/120.0);
*y = scale*(m1 - m0 + n*tan(y0)*(A*A/2.0 + (5.0 - t + 9.0*c + 4*c*c)*pow(A, 4.0)/24.0
+ (61.0 - 58.0*t + t*t + 600.0*c - 330.0*ecc2)*pow(A, 6.0)/720.0));
}

void fromTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y)
{
double x0, y0, x1, y1, phi, m, m0, mu, e1;
double c1, t1, n1, r1, d;

x0 = *x;
y0 = *y;

x1 = lon*DEG_RAD;
y1 = lat*DEG_RAD;

m0 = mercator(y1, a, ecc);
m = m0 + y0/scale;

e1 = (1.0-sqrt(1.0-ecc))/(1.0+sqrt(1.0-ecc));
mu = m / (a * (1.0 - ecc/4.0 - 3.0*ecc*ecc/64.0 - 5.0*ecc*ecc*ecc/256.0));

phi = mu + (3.0*e1/2.0 - 27.0*pow(e1, 3.0)/32.0)*sin(2.0*mu)
+ (21.0*e1*e1/16.0 - 55.0*pow(e1,4.0)/32.0)*sin(4.0*mu)
+ 151.0*pow(e1,3.0)/96.0*sin(6.0*mu) + 1097.0*pow(e1,4.0)/512.0*sin(8.0*mu);

c1 = ecc2*pow(cos(phi),2.0);
t1 = pow(tan(phi),2.0);
n1 = a/sqrt(1-ecc*pow(sin(phi),2.0));
r1 = a*(1.0-ecc)/pow(1.0-ecc*pow(sin(phi),2.0), 1.5);
d = x0/(n1*scale);

*x = (x1+(d-(1.0+2.0*t1+c1)*pow(d,3.0)/6.0
+(5.0-2.0*c1+28.0*t1-3.0*c1*c1+8.0*ecc2+24.0*t1*t1)*pow(d,5.0)/120.0)/cos(phi))*RAD_DEG;
*y = (phi-n1*tan(phi)/r1*(d*d/2.0-(5.0+3.0*t1+10.0*c1-4.0*c1*c1-9.0*ecc2)
*pow(d,4.0)/24.0+(61.0+90.0*t1+298.0*c1+45.0*t1*t1-252.0*ecc2-3.0*c1*c1)*pow(d,6.0)/72.0))*RAD_DEG;
}

double mercator(double y, double a, double ecc)
{
if(y == 0.0)
{
return 0.0;
}
else
{
return a * (
( 1.0 - ecc/4.0 - 3.0*ecc*ecc/64.0 - 5.0*ecc*ecc*ecc/256.0 ) * y -
( 3.0*ecc/8.0 + 3.0*ecc*ecc/32.0 + 45.0*ecc*ecc*ecc/1024.0 ) * sin(2.0 * y) +
( 15.0*ecc*ecc/256.0 + 45.0*ecc*ecc*ecc/1024.0 ) * sin(4.0 * y) -
( 35.0*ecc*ecc*ecc/3072.0 ) * sin(6.0 * y) );
}
}

/*************************************************************************************

Sample code below, using coordinate in Dan Jacob's website.

*/

#include <stdio.h>
#include <stdlib.h>

int main()
{
double x1, y1, z1, x2, y2, z2;
double tx1, ty1, tx2, ty2;
double dx, dy, dz, dx1, dy1;

x1 = 120.85788004; // TWD67
y1 = 24.18347242;
z1 = 777;

x2 = 120.86603958; // TWD97
y2 = 24.18170479;
z2 = 777;

tx1 = 235561; // TWD67->TM2
ty1 = 2675359;

tx2 = 236389.849; // TWD97->TM2
ty2 = 2675153.168;

/////////////////////////////////////////////
//
//
// convert TWD67->TM2
//
//
/////////////////////////////////////////////

dx = x1;
dy = y1;

toTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx, &dy);
// center longitude of taiwan is 121, for penghu is 119

dx += 250000; // TM2 in Taiwan should add 250000

printf("TWD67->TM2\nTWD67 (%f, %f)\nConvert (%.3f, %.3f)\nOrigin (%.3f, %.3f)\n", x1, y1, dx, dy, tx1, ty1);
printf("Acuuracy (%.3f, X:%.3f, Y:%.3f)\n\n", sqrt((dx-tx1)*(dx-tx1)+(dy-ty1)*(dy-ty1)), (dx-tx1), (dy-ty1));

/////////////////////////////////////////////
//
//
// convert TWD97->TM2
//
//
/////////////////////////////////////////////

dx = x2;
dy = y2;

toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx, &dy);
// center longitude of taiwan is 121, for penghu is 119

dx += 250000; // TM2 in Taiwan should add 250000

printf("TWD97->TM2\nTWD97 (%f, %f)\nConvert (%.3f, %.3f)\nOrigin (%.3f, %.3f)\n", x2, y2, dx, dy, tx2, ty2);
printf("Acuuracy (%.3f, X:%.3f, Y:%.3f)\n\n", sqrt((dx-tx2)*(dx-tx2)+(dy-ty2)*(dy-ty2)), (dx-tx2), (dy-ty2));

/////////////////////////////////////////////
//
//
// convert TM2->TWD67
//
//
/////////////////////////////////////////////

dx = tx1-250000; // should minus 250000 first in Taiwan
dy = ty1;

fromTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx, &dy);

printf("TM2->TWD67\nTM2 (%f, %f)\nConvert (%.9f, %.9f)\nOrigin (%.9f, %.9f)\n", tx1, ty1, dx, dy, x1, y1);
printf("Acuuracy (%.9f, X:%.9f, Y:%.9f)\n\n", sqrt((dx-x1)*(dx-x1)+(dy-y1)*(dy-y1)), (dx-x1), (dy-y1));

/////////////////////////////////////////////
//
//
// convert TM2->TWD97
//
//
/////////////////////////////////////////////

dx = tx2-250000; // should minus 250000 first in Taiwan
dy = ty2;

fromTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx, &dy);

printf("TM2->TWD97\nTM2 (%f, %f)\nConvert (%.9f, %.9f)\nOrigin (%.9f, %.9f)\n", tx2, ty2, dx, dy, x2, y2);
printf("Acuuracy (%.9f, X:%.9f, Y:%.9f)\n\n", sqrt((dx-x2)*(dx-x2)+(dy-y2)*(dy-y2)), (dx-x2), (dy-y2));

/////////////////////////////////////////////
//
//
// convert TWD67->TWD97
//
//
/////////////////////////////////////////////

dx = x1;
dy = y1;
dz = z1;

toTWD97(&dx, &dy, &dz);

dx1 = dx;
dy1 = dy;

toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx1, &dy1);

dx1 += 250000; // TM2 in Taiwan should add 250000

printf("TWD67->TWD97\nTWD67 (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x1, y1, z1, tx1, ty1);
printf("Convert (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", dx, dy, dz, dx1, dy1);
printf("Origin (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x2, y2, z2, tx2, ty2);
printf("Acuuracy (%.4f, X:%.4f, Y:%.4f)\n\n", sqrt((dx1-tx2)*(dx1-tx2)+(dy1-ty2)*(dy1-ty2)), (dx1-tx2), (dy1-ty2));

/////////////////////////////////////////////
//
//
// convert TWD97->TWD67
//
//
/////////////////////////////////////////////

dx = x2;
dy = y2;
dz = z2;

toTWD67(&dx, &dy, &dz);

dx1 = dx;
dy1 = dy;

toTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx1, &dy1);

dx1 += 250000; // TM2 in Taiwan should add 250000

printf("TWD97->TWD67\nTWD97 (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x2, y2, z2, tx2, ty2);
printf("Convert (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", dx, dy, dz, dx1, dy1);
printf("Origin (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x1, y1, z1, tx1, ty1);
printf("Acuuracy (%.4f, X:%.4f, Y:%.4f)\n\n", sqrt((dx1-tx1)*(dx1-tx1)+(dy1-ty1)*(dy1-ty1)), (dx1-tx1), (dy1-ty1));
}
(繼續閱讀)

, , ::意見(7)

台灣座標系統相關參數(TWD67, TWD97, WGS84, Hu-Tzu-Shan)

2004年6月1日 18:42

底下是目前所瞭解的台灣各座標系統的詳細參數. 不過各參數間似乎有著不小的差距.

 



1. WGS84
a: 6378137.0
b: 6356752.3142
f: 298.257223563
dx: 0
dy: 0
dz: 0
rx: 0
ry: 0
rz: 0
s: 0

2. Hu-Tzu-Shan(虎子山, International 1924橢球)
a: 6378388.0
b: 6356911.9461
f: 297
dx: -637
dy: -549
dz: -203
rx: 0
ry: 0
rz: 0
s: 0

3. TWD97(GRS80橢球)
a: 6378137.0
b: 6356752.3141
f: 298.257222101
dx: 0
dy: 0
dz: 0
rx: 0
ry: 0
rz: 0
s: 0
來源: 內政部

4. TWD67-1(International 1967 橢球)
a: 6378160.0
b: 6356774.7192
f: 298.25
dx: -764.558
dy: -359.515
dz: -180.510
rx: 0.00003863
ry: 0.00001721
rz: 0.00000197
s: 0.99998180
5. TWD67-2(International 1967 橢球)
a: 6378160.0
b: 6356774.7192
f: 298.25
dx: -750.739
dy: -361.229
dz: -178.374
rx: -0.0000011698
ry: 0.0000018398
rz: 0.0000009822
s: 0.00002329
6. TWD67-3(橢球未知)
a: 6378155.0
b: 6356769.7359
f: 298.257223481
dx: -767
dy: -358
dz: -176
rx: 0
ry: 0
rz: 0
s: 0
來源: Garmin, 見積丹尼網站
(繼續閱讀)

, ::意見(1)

座標轉換 - 三參數法與七參數法

2004年5月25日 22:00

座標轉換通常使用兩種方法. 第一種是平移, 第二種是平移加旋轉. 平移法要使用到三個參數, 通稱三參數法(3 parameter formula). 平移旋轉法要使用到七個參數, 通稱七參數法(7 parameter formula). 底下為兩種參數法的解釋說明:

Common Transformation Models

If the XYZ coordinate axes of two datums are known to be parallel and identically scaled, a three parameter transformation can be derived to represent their relationship (see Diagram 13 and Insert 5).

If the coordinate axes are not parallel and identically scaled, a seven parameter transformation can be derived (see Diagram 14 and Insert 6).

Diagram 13
3-Parameter Transformation

Insert 5
Three-Parameter Transformation Formulae

X1 = X2 + DX Y1 = Y2 + DY Z1 = Z2 + DZ

where

X1, Y1, Z1 = Cartesian Coordinates of Datum 1
X2, Y2, Z2 = Cartesian Coordinates of Datum 2
DX, DY, DZ = The difference between the centres of the two spheroids


Diagram 14
7-Parameter Transformation

Insert 6
Seven Parameter Transformation Formulae

(Bursa-Wolf Model)

where

X1, Y1, Z1 = Cartesian Coordinates of Datum 1
X2, Y2, Z2 = Cartesian Coordinates of Datum 2
DX, DY, DZ = The difference between the centres of the two spheroids
RX, RY, RZ = The rotations around the three coordinate axes
SC = The scale difference between the coordinate systems

Rotations are positive anticlockwise about the axes of Datum 2 coordinate system when viewing the origin from the positive axes.

全文參照: http://www.environment.sa.gov.au/mapland/sicom/sicom/tp_scs.html
(繼續閱讀)

::意見(0)

座標系轉換參數

2004年5月24日 23:28

底下是一些座標系與轉換參數.

出處為 http://www.colorado.edu/geography/gcraft/notes/datum/edlist.html



橢球:

EllipsoidSemi-major axis 1/flattening
Airy 1830, 6377563.396 299.3249646
Modified Airy6377340.189 299.3249646
Australian National 6378160298.25
Bessel 1841 (Namibia) 6377483.865299.1528128
Bessel 18416377397.155 299.1528128
Clarke 1866,6378206.4 294.9786982
Clarke 1880,6378249.145 293.465
Everest (India 1830)" 6377276.345300.8017
Everest (Sabah Sarawak) 6377298.556300.8017
Everest (India 1956)6377301.243300.8017
Everest (Malaysia 1969) 6377295.664300.8017
Everest (Malay. & Sing) 6377304.063300.8017
Everest (Pakistan)6377309.613300.8017
Modified Fischer 19606378155298.3
Helmert 19066378200 298.3
Hough 19606378270 297
Indonesian 19746378160 298.247
International 19246378388297
Krassovsky 19406378245 298.3
GRS 806378137 298.257222101
South American 19696378160298.25
WGS 726378135 298.26
WGS 846378137 298.257223563

座標系與參數(3 parameters):

d = delta in meters; e = error estimate in meters; #S = number of satellite measurement stations

Datum Ellipsoid dXdYdZ Region of use eX eYeZ #S
AdindanClarke 1880-118 -14218Burkina Faso 2525251
AdindanClarke 1880-134 -2210Cameroon 2525251
AdindanClarke 1880-165 -11206Ethiopia 3338
AdindanClarke 1880-123 -20220Mali 2525251
AdindanClarke 1880-166 -15204MEAN FOR Ethiopia; Sudan 55322
AdindanClarke 1880-128 -18224Senegal 2525252
AdindanClarke 1880-161 -14205Sudan 35314
AfgooyeKrassovsky 1940-43-16345Somalia2525251
Ain el Abd 1970International 1924 -150-250-1 Bahrain2525 252
Ain el Abd 1970International 1924 -143-2367 Saudi Arabia1010 109
American Samoa 1962Clarke 1866 -115118426 American Samoa Islands2525252
Anna 1 Astro 1965Australian National -491-22435 Cocos Islands2525 251
Antigua Island Astro 1943Clarke 1880 -2701362 Antigua (Leeward Islands)25 25251
Arc 1950Clarke 1880-138-105-289 Botswana35 39
Arc 1950Clarke 1880-153-5-292Burundi2020203
Arc 1950Clarke 1880-125-108-295 Lesotho33 85
Arc 1950Clarke 1880-161-73-317 Malawi924 86
Arc 1950Clarke 1880-143-90-294 MEAN FOR Botswana; Lesotho; Malawi; Swaziland; Zaire; Zambia; Zimbabwe 20332041
Arc 1950Clarke 1880-134-105-295 Swaziland1515 154
Arc 1950Clarke 1880-169-19-278 Zaire2525 252
Arc 1950Clarke 1880-147-74-283 Zambia2121 275
Arc 1950Clarke 1880-142-96-293 Zimbabwe58 1110
Arc 1960Clarke 1880-160-6-302MEAN FOR Kenya; Tanzania20202025
Arc 1960Clarke 1880-157-2-299Kenya43324
Arc 1960Clarke 1880-175-23-303 Taanzania69 1012
Ascension Island 1958International 1924 -20510753 Ascension Island2525 252
Astro Beacon E 1945International 1924 14575-272 Iwo Jima2525 251
Astro DOS 71/4International 1924 -320550-494 St Helena Island2525 251
Astro Tern Island (FRIG) 1961International 1924 114-116-333 Tern Island2525 251
Astronomical Station 1952International 1924 124-234-25 Marcus Island2525 251
Australian Geodetic 1966Australian National -133-48148 Australia; Tasmania333105
Australian Geodetic 1984Australian National -134-48149 Australia; Tasmania22290
Ayabelle LighthouseClarke 1880 -79-129145 Djibouti2525 251
Bellevue (IGN)International 1924 -127-769472 Efate & Erromango Islands20 20203
Bermuda 1957Clarke 1866-73213296Bermuda2020203
BissauInternational 1924 -17325327 Guinea-Bissau2525 252
Bogota ObservatoryInternational 1924 307304-318 Colombia65 67
Bukit RimpahBessel 1841-384664-48Indonesia (Bangka & Belitung Ids)-1 -1-10
Camp Area AstroInternational 1924 -104-129239 Antarctica (McMurdo Camp Area)-1 -1-10
Campo InchauspeInternational 1924 -14813690 Argentina55 520
Canton Astro 1966International 1924 298-304-375 Phoenix Islands1515 154
CapeClarke 1880-136 -108-292South Africa 3665
Cape CanaveralClarke 1866 -2151181 Bahamas; Florida33 319
CarthageClarke 1880-2636431Tunisia6985
Chatham Island Astro 1971International 1924 175-38113 New Zealand (Chatham Island)15 15154
Chua AstroInternational 1924 -134229-29 Paraguay69 56
Corrego AlegreInternational 1924 -206172-6 Brazil53 517
DabolaClarke 1880-83 37124Guinea 1515154
Deception IslandClarke 1880 26012-147 Deception Island; Antarctia20 20203
Djakarta (Batavia)Bessel 1841 -377681-50 Indonesia (Sumatra)3335
DOS 1968International 1924 230-199-752 New Georgia Islands (Gizo Island)25 25251
Easter Island 1967International 1924 211147111 Easter Island2525 251
Estonia; Coordinate System 1937Bessel 1841 374150588 Estonia23 319
European 1950International 1924 -104-101-140 Cyprus1515 154
European 1950International 1924 -130-117-151 Egypt68 814
European 1950International 1924 -86-96-120 England; Channel Islands; Scotland; Shetland Islands 33340
European 1950International 1924 -86-96-120 England; Ireland; Scotland; Shetland Islands 33347
European 1950International 1924 -87-95-120 Finland; Norway35 320
European 1950International 1924 -84-95-130 Greece2525 252
European 1950International 1924 -117-132-164 Iran912 1127
European 1950International 1924 -97-103-120 Italy (Sardinia)2525 252
European 1950International 1924 -97-88-135 Italy (Sicily)2020 203
European 1950International 1924 -107-88-149 Malta2525 251
European 1950International 1924 -87-98-121 MEAN FOR Austria; Belgium; Denmark; Finland; France; W Germany; Gibraltar; Greece; Italy; Luxembourg; Netherlands; Norway; Portugal; Spain; Sweden; Switzerland 38585
European 1950International 1924 -87-96-120 MEAN FOR Austria; Denmark; France; W Germany; Netherlands; Switzerland 33352
European 1950International 1924 -103-106-141 MEAN FOR Iraq; Israel; Jordan; Lebanon; Kuwait; Saudi Arabia; Syria -1-1-10
European 1950International 1924 -84-107-120 Portugal; Spain56 318
European 1950International 1924 -112-77-145 Tunisia2525 254
European 1979International 1924 -86-98-119 MEAN FOR Austria; Finland; Netherlands; Norway; Spain; Sweden; Switzerland 33322
Fort Thomas 1955Clarke 1880 -7215225 Nevis; St. Kitts (Leeward Islands)25 25252
Gan 1970International 1924 -133-32150 Republic of Maldives2525251
Geodetic Datum 1949International 1924 84-22209 New Zealand53 514
Graciosa Base SW 1948International 1924 -104167-38 Azores (Faial; Graciosa; Pico; Sao Jorge; Terceira) 3335
Guam 1963Clarke 1866-100-248259 Guam333 5
Gunung SegaraBessel 1841 -40368441 Indonesia (Kalimantan)-1-1-10
GUX 1 AstroInternational 1924 252-209-751 Guadalcanal Island2525251
Herat NorthInternational 1924 -333-222114 Afghanistan-1-1 -10
Hermannskogel DatumBessel 1841 (Namibia) 653-212449 Croatia -Serbia, Bosnia-Herzegovina -1 -1-1 0
Hjorsey 1955International 1924 -7346-86 Iceland33 66
Hong Kong 1963International 1924 -156-271-189 Hong Kong2525 252
Hu-Tzu-ShanInternational 1924 -637-549-203 Taiwan1515 154
IndianEverest (India 1830) 282726254 Bangladesh108 126
IndianEverest (India 1956) 295736257 India; Nepal1210 157
IndianEverest (Pakistan) 283682231 Pakistan-1-1 -10
Indian 1954Everest (India 1830) 217823299 Thailand156 1211
Indian 1960Everest (India 1830) 182915344 Vietnam (Con Son Island)25 25251
Indian 1960Everest (India 1830) 198881317 Vietnam (Near 16øN)25 25252
Indian 1975Everest (India 1830) 210814289 Thailand32 362
Indonesian 1974Indonesian 1974 -24-155 Indonesia2525 251
Ireland 1965Modified Airy 506-122611 Ireland33 37
ISTS 061 Astro 1968International 1924 -794119-298 South Georgia Islands2525251
ISTS 073 Astro 1969International 1924 208-435-229 Diego Garcia2525 252
Johnston Island 1961International 1924 189-79-202 Johnston Island2525 251
KandawalaEverest (India 1830) -9778786 Sri Lanka2020 203
Kerguelen Island 1949International 1924 145-187103 Kerguelen Island2525 251
Kertau 1948Everest (Malay. & Sing) -118515 West Malaysia & Singapore10 866
Kusaie Astro 1951International 1924 6471777-1124 Caroline Islands2525 251
Korean Geodetic SystemGRS 80 000South Korea22 212
L. C. 5 Astro 1961Clarke 1866 42124147 Cayman Brac Island2525251
LeigonClarke 1880-130 29364Ghana 2328
Liberia 1964Clarke 1880-904088Liberia1515154
LuzonClarke 1866-133 -77-51Philippines (Excluding Mindanao) 81196
LuzonClarke 1866-133 -79-72Philippines (Mindanao) 2525251
M'PoralokoClarke 1880-74-13042Gabon2525251
Mahe 1971Clarke 188041-220-134Mahe Island2525 251
MassawaBessel 1841639 40560Ethiopia (Eritrea) 2525251
MerchichClarke 18803114647Morocco5339
Midway Astro 1961International 1924 912-581227 Midway Islands2525 251
MinnaClarke 1880-81 -84115Cameroon 2525252
MinnaClarke 1880-92 -93122Nigeria 3656
(繼續閱讀)

::意見(0)