Top Level Namespace

Includes:
Math

Defined Under Namespace

Modules: CfSim

Instance Method Summary collapse

Instance Method Details

#bl2xy(phi, lamda) ⇒ Object

緯度経度を平面直角座標に変換する計算



41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
 
# File 'lib/cblxy.rb', line 41

def bl2xy(phi, lamda)
   $e = sqrt(2.0 * $F - 1.0) / $F

   phi1 = deg2rad(phi)
   lamda1 = deg2rad(lamda)

   s0 = kocyou($phi0)
   s1 = kocyou(phi1)

   ut = $a / sqrt(1.0 - $e ** 2.0 * sin(phi1) ** 2.0)
   conp = cos(phi1)
   t1 = tan(phi1)
   dlamda = lamda1 - $lamda0
   eta2 = ($e ** 2.0 / (1.0 - $e ** 2.0)) * conp ** 2.0

   #X座標値の算出
   v1 = 5.0 - t1 ** 2.0 + 9.0 * eta2 + 4.0 * eta2 ** 2.0
   v2 = -61.0 + 58.0 * t1 ** 2.0 - t1 ** 4.0 - 270.0 * eta2 + 330.0 * t1 ** 2.0 * eta2
   v3 = -1385.0 + 3111.0 * t1 ** 2.0 - 543.0 * t1 ** 4.0 + t1 ** 6.0

   x = ((s1 - s0) + ut * conp ** 2.0 * t1 * dlamda ** 2.0 / 2.0 + ut * conp ** 4.0 * t1 * v1 * dlamda ** 4.0 / 24.0 - ut * conp ** 6.0 * t1 * v2 * dlamda ** 6.0 / 720.0 - ut * conp ** 8.0 * t1 * v3 * dlamda ** 8.0 / 40320.0) * $sbyS

   #Y座標値の算出
   v1 = -1.0 + t1 ** 2.0 - eta2
   v2 = -5.0 + 18.0 * t1 ** 2.0 - t1 ** 4.0 - 14.0 * eta2 + 58.0 * t1 ** 2.0 * eta2
   v3 = -61.0 + 479.0 * t1 ** 2.0 - 179.0 * t1 ** 4.0 + t1 ** 6.0

   y = (ut * conp * dlamda - ut * conp ** 3.0 * v1 * dlamda ** 3.0 / 6.0 - ut * conp ** 5.0 * v2 * dlamda ** 5.0 / 120.0 - ut * conp ** 7.0 * v3 * dlamda ** 7.0 / 5040.0) * $sbyS

   return [x,y]

end

#blxy(phi1, lamda1, k) ⇒ Object

世界測地系で緯度経度を平面直角座標に変換するメソッド



22
23
24
25
26
27
28
29
 
# File 'lib/cblxy.rb', line 22

def blxy(phi1, lamda1,k)
   $a = $aG #GRS80楕円体の定数を代入
   $F = $FG
   genten(k)

   bl2xy(phi1, lamda1)

end

#deg2rad(deg) ⇒ Object

度分秒からラジアンへの変換



177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
 
# File 'lib/cblxy.rb', line 177

def deg2rad(deg)
   radbase = PI / 180.0
   if deg < 0
      fugou  = -1.0
   else
      fugou = 1.0
   end

   deg = deg.abs

   angle = (deg / 10000.0).to_i
   minute = ((deg / 100.0).to_i - (angle * 100.0)).to_f
   second =  deg - (angle * 10000.0) - (minute * 100.0)
   rad = fugou * (angle + (minute + (second / 60.0)) / 60.0)
   rad = rad * radbase

end

#genten(k) ⇒ Object

座標系原点の緯度経度



118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
 
# File 'lib/cblxy.rb', line 118

def genten(k)

   err = "座標系が正しく指定されていません"

   case k
   when 1
      degen = [330000.0,1293000.0]
   when 2
      degen = [330000.0,1310000.0]
   when 3
      degen = [360000.0,1321000.0]
   when 4
      degen = [330000.0,1333000.0]
   when 5
      degen = [360000.0,1342000.0]
   when 6
      degen = [360000.0,1360000.0]
   when 7
      degen = [360000.0,1371000.0]
   when 8
      degen = [360000.0,1383000.0]
   when 9
      degen = [360000.0,1395000.0]
   when 10
      degen = [400000.0,1405000.0]
   when 11
      degen = [440000.0,1401500.0]
   when 12
      degen = [440000.0,1421500.0]
   when 13
      degen = [440000.0,1441500.0]
   when 14
      degen = [260000.0,1420000.0]
   when 15
      degen = [260000.0,1273000.0]
   when 16
      degen = [260000.0,1240000.0]
   when 17
      degen = [260000.0,1310000.0]
   when 18
      degen = [200000.0,1360000.0]
   when 19
      degen = [260000.0,1540000.0]
   else
      degen = err
   end

   unless degen == err
      $phi0 = (deg2rad degen[0])
      $lamda0 = (deg2rad degen[1])
      return [$phi0,$lamda0]
   else
      return err
   end
end

#kocyou(ido) ⇒ Object

緯度から赤道への子午線弧長を計算する



218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
 
# File 'lib/cblxy.rb', line 218

def kocyou(ido)
   e2 = $e ** 2.0
   e4 = $e ** 4.0
   e6 = $e ** 6.0
   e8 = $e ** 8.0
   e10 = $e ** 10.0
   e12 = $e ** 12.0
   e14 = $e ** 14.0
   e16 = $e ** 16.0

   a = 1.0 + 3.0 / 4.0 * e2 + 45.0 / 64.0 * e4 + 175.0 / 256.0 * e6 + 11025.0 / 16384.0 * e8 + 43659.0 / 65536.0 * e10 + 693693.0 / 1048576.0 * e12 + 19324305.0 / 29360128.0 * e14 + 4927697775.0 / 7516192768.0 * e16
   b = 3.0 / 4.0 * e2 + 15.0 / 16.0 * e4 + 525.0 / 512.0 * e6 + 2205.0 / 2048.0 * e8 + 72765.0 / 65536.0 * e10 + 297297.0 / 262144.0 * e12 + 135270135.0 / 117440512.0 * e14 + 547521975.0 / 469762048.0 * e16
   c = 15.0 / 64.0 * e4 + 105.0 / 256.0 * e6 + 2205.0 / 4096.0 * e8 + 10395.0 / 16384.0 * e10 + 1486485.0 / 2097152.0 * e12 + 45090045.0 / 58720256.0 * e14 + 766530765.0 / 939524096.0 * e16
   d = 35.0 / 512.0 * e6 + 315.0 / 2048.0 * e8 + 31185.0 / 131072.0 * e10 + 165165.0 / 524288.0 * e12 + 45090045.0 / 117440512.0 * e14 + 209053845.0 / 469762048.0 * e16
   e = 315.0 / 16384.0 * e8 + 3465.0 / 65536.0 * e10 + 99099.0 / 1048576.0 * e12 + 4099095.0 / 29360128.0 * e14 + 348423075.0 / 1879048192.0 * e16
   f = 693.0 / 131072 * e10 + 9009.0 / 524288.0 * e12 + 4099095.0 / 117440512.0 * e14 + 26801775.0 / 469762048.0 * e16
   g = 3003 / 2097152.0 * e12 + 315315.0 / 58720256.0 * e14 + 11486475.0 / 939524096.0 * e16
   h = 45045.0 / 117440512.0 * e14 + 765765.0 / 469762048.0 * e16
   i = 765765.0 / 7516192768.0 * e16

   sigosen = $a * (1.0 - e2) * (a * ido - b * sin(ido * 2.0) / 2.0 + c * sin(ido * 4.0) / 4.0 - d * sin(ido * 6.0) / 6.0 + e * sin(ido * 8.0) / 8.0 - f * sin(ido * 10.0) / 10.0 + g * sin(ido * 12.0) / 12.0 - h * sin(ido * 14.0) / 14.0 + i * sin(ido * 16.0) / 16.0)
end

#rad2deg(rad) ⇒ Object

ラジアンから度分秒への変換



196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
 
# File 'lib/cblxy.rb', line 196

def rad2deg(rad)

   degbase = 180.0 / PI

   if rad < 0
      fugou = -1.0
   else
      fugou = 1.0
   end

   rad = rad.abs

   rad = rad * degbase
   angle = rad.to_i
   rad = (rad - angle).to_f * 60.0
   minute = rad.to_i
   second = (rad - minute).to_f * 60.0

   deg = fugou * (angle * 10000.0 + minute * 100.0 + second)
end

#suisen(x) ⇒ Object 



241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
 
# File 'lib/cblxy.rb', line 241

def suisen(x)
   s0 = kocyou($phi0)
   m = s0 + x / $sbyS
   cnt = 0
   phin = $phi0
   e2 = $e ** 2.0
   phi0 = phin

   while
      cnt = cnt + 1
      phi0 = phin
      sn = kocyou(phin)
      v1 = 2.0 * (sn - m) * ((1.0 - e2 * sin(phin) ** 2.0) ** 1.5)
      v2 = 3.0 * e2 * (sn - m) * sin(phin) * cos(phin) * sqrt(1.0 - e2 * sin(phin) ** 2.0) - 2.0 * $a * (1.0 - e2)
      phin = phin +v1 / v2
      if ((phin - phi0).abs < 0.00000000000001) or cnt > 100
         break
      end
   end
   return phin

end

#xy2bl(x, y) ⇒ Object

平面直角座標から緯度経度に変換する計算



75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
 
# File 'lib/cblxy.rb', line 75

def xy2bl(x,y)
   $e = sqrt(2.0 * $F - 1.0) / $F

   phi1 = suisen(x)

   ut = $a / sqrt( 1.0 - $e ** 2.0 * sin(phi1) ** 2.0)
   conp = cos(phi1)
   t1 = tan(phi1)
   eta2 = ($e ** 2.0 / (1.0 - $e ** 2.0)) * conp ** 2.0

   #緯度算出
   yy = y / $sbyS
   v1 = 1.0 + eta2
   v2 = 5.0 + 3.0 * t1 ** 2.0 + 6.0 * eta2 - 6.0 * t1 ** 2.0 * eta2 - 3.0 * eta2 ** 2.0 - 9.0 * t1 ** 2.0 * eta2 ** 2.0
   v3 = 61.0 + 90.0 * t1 ** 2.0 + 45.0 * t1 ** 4.0 + 107.0 * eta2 - 162.0 * t1 ** 2.0 * eta2 - 45.0 * t1 ** 4.0 * eta2
   v4 = 1385.0 + 3633.0 * t1 ** 2.0 + 4095.0 * t1 ** 4.0 + 1575.0 * t1 ** 6.0

   phir = -(v1 / (2.0 * ut ** 2.0)) * yy ** 2.0
   phir = phir + (v2 / (24.0 * ut ** 4.0)) * yy ** 4.0
   phir = phir - (v3 / (720.0 * ut ** 6.0)) * yy ** 6.0
   phir = phir + (v4 / (40320.0 * ut ** 8.0)) * yy ** 8.0
   phir = phir * t1
   phir = phir + phi1
   phir = rad2deg(phir)

   #経度算出
   v1 = ut * conp
   v2 = 1.0 + 2.0 * t1 ** 2.0 + eta2
   v3 = 5.0 + 28.0 * t1 ** 2.0 + 24.0 * t1 ** 4.0 + 6.0 * eta2 + 8.0 * t1 ** 2.0 * eta2
   v4 = 61.0 + 662.0 * t1 ** 2.0 + 1320.0 * t1 ** 4.0 + 720.0 * t1 ** 6.0

   lamdar = (1.0 / v1) * yy
   lamdar = lamdar - (v2 / (6.0 * ut ** 2.0 * v1)) * yy ** 3.0
   lamdar = lamdar + (v3 / (120.0 * ut ** 4.0 * v1)) * yy ** 5.0
   lamdar = lamdar - (v4 / (5040.0 * ut ** 6.0 * v1)) * yy ** 7.0
   lamdar = lamdar + $lamda0

   lamdar = rad2deg(lamdar)

   return [phir,lamdar]
end

#xybl(phi1, lamda1, k) ⇒ Object

世界測地系で平面直角座標を緯度経度に変換するメソッド



32
33
34
35
36
37
38
 
# File 'lib/cblxy.rb', line 32

def xybl(phi1,lamda1,k)
   $a = $aG #GRS80楕円体の定数を代入
   $F = $FG
   genten(k)

   xy2bl(phi1, lamda1)
end