#From: Columbia Yen <tech547@yahoo.com> Nov 12 2001
#Subject: 台電座標轉2度TM, 及經緯度
#
#   丹尼兄： 小弟在 2001.10.22 於 PDA 板向您討教過的 tech，現在我已經 把程
#   式寫了個大概，提供給您，你可以把程式放在您的網頁上。 我的程式是版權沒
#   有，歡迎引用。 這個 UTM2.bas 是 palm 上的 basic 程式，需要用到 Hotpaw
#   BASIC (ybasic.prc) 及 math.lib，目前可以正確地做台灣本島的台電 座
#   標<->2 度 TM 之間的轉換，而經緯度部份在經度與 proj 的計算 結果是一致的
#   ，但緯度就有些誤差了，約在十秒以下。
#
#   Dan> 好呀，暫掛我這兒但建議您劈開自己的網站， 以便於更新。我本人也不
#   會 palm... 至於誤差：大概可參 proj source code ?
#
#   以下是 UTM2 的操作說明，而 Hotpaw BASIC 的安裝及操作就請恕小弟偷懶了。
#    1. 首先以 G8150, HD78 為例子。
#    2. 在畫面中的【大區座標】右方輸入「G8150」(大小寫均可)。
#    3. 在畫面中的【小區座標】右方輸入「HD78」(大小寫均可)。
#    4. 點【Calc】，則開始計算，然後在【TM2-X】、【TM2-Y】 、【東經】及【
#       北緯】右方出現計算結果。
#          TM2-X  235570
#          TM2-Y 2675380
#          東經 120 度 51 分 28.69 秒
#          北緯  24 度 10 分 54.98 秒
#    5. 如果【大區座標】及【小區座標】是空白的，而【TM2-X】及 【TM2-Y】有
#       資料，點【Calc】，則可以反算【大區座標】及 【小區座標】，同時經緯
#       度也會算出來。
#    6. 另外【Calc2】的功用是清除資料，以便重新輸入，【Done】 是結束程式，
#       【OK】目前沒有功用，未來將作為把計算結果 寫入到memo中。
#
#UTM2.bas
# 台電座標轉2度TM座標
n$="感謝積丹尼提供資料"
sc=1805.1
dim s$(17), xb(23), yb(23)
dim xa$(22)
for i=0 to 17: s$(i)="": next i
for i=0 to 22: xb(i)=0: yb(i)=0:next i
for i=0 to 22: xa$(i)="": next i
s$(0)=n$
s$(2)="大區座標"
s$(4)="小區座標"
s$(6)="TM2-X"
s$(8)="TM2-Y"
s$(10)="東經"
s$(14)="北緯"
a=-1:j=0
while(j&lt;3)
for i=0 to 2
            a=a+1
            xa$(a)=chr$(a+65)
             xb(a)=170000+i*80000
             yb(a)=2750000-j*50000
 next i
j=j+1
wend
j=0
while(j&lt;4)
for i=0 to 2
        a=a+1
        xa$(a)=chr$(a+65)
        xb(a)=90000+i*80000
        yb(a)=2600000-j*50000
 next i
j=j+1
wend
a=21
xa$(a)=chr$(a+65)
xb(a)=170000
yb(a)=2400000

while (1)
 j = form(12,0,"台電座標轉2度TM座標")
if j=0 then end
if j=1 then gosub 1000
#clear variable
 if j=3 then gosub 100
's$(3)="G8150"
's$(5)="HD78"

#==================
#====input 台電座標
#==================
if ((len(s$(3))&gt;0) and (j=2)) then
gosub 500
j=1
endif
#==================
#====END   input 台電座標
#==================

#==================
#====input TM-2
#==================
if ((len(s$(7))&gt;0)  and (j=2)) then
gosub 600
j=1
endif
#==================
#====END  input TM-2
#==================
wend
end

#==================
#====clear 
#==================
100 for i=3 to 17 step 2
s$(i)=""
next i
return
#==================
#====clear 
#==================

#==================
#====input 台電座標
#==================
500 s$(3)=ucase$(s$(3))
a=asc(mid$(s$(3),1,1))-asc("A")
b=val(mid$(s$(3),2,2))
c=val(mid$(s$(3),4,2))

s$(5)=ucase$(s$(5))
d=asc(mid$(s$(5),1,1))-asc("A")
e=asc(mid$(s$(5),2,1))-asc("A")
f=val(mid$(s$(5),3,1))
g=val(mid$(s$(5),4,1))
k=val(mid$(s$(5),5,1))
l=val(mid$(s$(5),6,1))
x=xb(a)+b*800+d*100+f*10+k
y=yb(a)+c*500+e*100+g*10+l
s$(7)=str$(x)
s$(9)=str$(y)
gosub 2000
return
#==================
#====END   input 台電座標
#==================

#==================
#====input TM-2
#==================
600 a=-1: z=0
while (z=0)
a=a+1
x=val(s$(7))
y=val(s$(9))
t=x-xb(a)
u=y-yb(a)
if ((t&gt;=0) and (t&lt;80000)) then
    if ((u&gt;=0) and (u&lt;50000))  then
        z=1
  endif
endif
wend
k=t mod 10
f=int((t mod 100)/10)
b=int(t/800)
cb$=str$(b)
if (len(cb$)&lt;2) then cb$="0"+cb$
d=int((t-b*800)/100)
l=u mod 10
g=int((u mod 100)/10)
c=int(u/500)
cc$=str$(c)
if (len(cc$)&lt;2) then cc$="0"+cc$
e=int((u-c*500)/100)
s$(3)=chr$(a+65)+cb$
s$(3)=s$(3)+cc$
s$(5)=chr$(d+65)+chr$(e+65)
s$(5)=s$(5)+str$(f)+str$(g)
if ((k&lt;&gt;0) or (l&lt;&gt;0)) then
s$(5)=s$(5)+str$(k)+str$(l)
endif
s$(7)=str$(x)
s$(9)=str$(y)
gosub 2000
return
#==================
#====END   input TM-2
#================== 
#===save DATA
#==================
1000 print "DATA SAVING..."
return
#==================
#===convert longtitude---ARCVIEW
2000 q=(x-250000)
q=q*0.000001/0.9999
nd=23+sc/3600
nd=nd-(2600015.4208-y)/110745.08
Nrad=nd/180*4*atn(1)
u=1/3600/180*4*atn(1)
u=sin(u)
Rn=(sin(Nrad)^2)
Rn=((1-0.00669460532856*Rn)^0.5)
Rn=6378160/Rn
Rm=(sin(Nrad)^2)
Rm=(1-0.00669460532856*Rm)^1.5
Rm=(1-0.00669460532856)/Rm
Rm=6378160*Rm
e1=1/(Rn*u*cos(Nrad))
e2=(Rn/Rm+2*(tan(Nrad)^2))
e2=e2/(3*(Rn^3)*u*cos(Nrad))
e3=(1+3*(tan(Nrad)^2))
e3=(tan(Nrad)^2)*e3
e3=e3/(15*(Rn^5)*u*cos(Nrad))
n1=tan(Nrad)/(2*Rm*Rn*u)
n2=5+3*(tan(Nrad)^2)
n2=n2*tan(Nrad)
n2=n2/(24*Rm*(Rn^3)*u)
Edegree=e1*q-e2*(q^3)+e3*(q^5)
Edegree=Edegree*1000000/3600
Edegree=121+Edegree
x1=int(Edegree)
x2=int((Edegree-x1)*60)
x3=((Edegree-x1)*60-x2)*60
x3=int(x3*100)/100
p1=(n1*(q^2)-n2*(q^4))*1000000
p1=p1/3600
p2=45+(121-x1)*36-x2
p3=(x1-120)*1.32
p3=0.00004106790125*p3
Ndegree=nd-p1+p2*p3
y1=int(Ndegree)
y2=int((Ndegree-y1)*60)
y3=((Ndegree-y1)*60-y2)*60
y3=int(y3*100)/100
s$(10)="東經"
s$(11)=str$(x1)+"度"+str$(x2)+"分"
s$(13)= str$(x3)+"秒"
s$(14)="北緯"
s$(15)=str$(y1)+"度"+str$(y2)+"分"
s$(17)=str$(y3)+"秒"
return
#====END convert longtitude