When::Ephemeris

     # File when/ephemeris.rb, line 204
204:   def _rot(x, y, t)
205:     c, s = cosc(t), sinc(t)
206:     return [x*c - y*s, x*s + y*c]
207:   end

     # File when/ephemeris.rb, line 188
188:   def _to_p2(x, y)
189:     return [0.0, 0.0] if x==0 && y==0
190:     return [atan2(y,x)/CIRCLE, sqrt(x*x+y*y)]
191:   end

     # File when/ephemeris.rb, line 172
172:   def _to_p3(x, y, z)
173:     ll, rr = _to_p2(x,  y)
174:     tb, rr = _to_p2(rr, z)
175:     return [ll, tb, rr]
176:   end

     # File when/ephemeris.rb, line 155
155:   def _to_r3(ll, tb, rr)
156:     c, s = cosc(tb), sinc(tb)
157:     return [rr*c*cosc(ll), rr*c*sinc(ll), rr*s]
158:   end

     # File when/ephemeris.rb, line 109
109:   def acos(x)
110:     atan2(sqrt(1-x*x), x)
111:   end

     # File when/ephemeris.rb, line 104
104:   def asin(x)
105:     atan2(x, sqrt(1-x*x))
106:   end

     # File when/ephemeris.rb, line 129
129:   def cosc(x)
130:     cos(x * CIRCLE)
131:   end

     # File when/ephemeris.rb, line 114
114:   def cosd(x)
115:     cos(x * DEG)
116:   end

     # File when/ephemeris.rb, line 233
233:   def delta_e(c)
234:     trigonometric(c, 
235:       [[COS   ,  125.04      ,   1934.136    , 0.00256  ],
236:        [COS   ,  200.93      ,  72001.539    , 0.00016  ]]) / 360
237:   end

     # File when/ephemeris.rb, line 246
246:   def delta_p(c)
247:     trigonometric(c, 
248:       [[SIN   ,  125.04      ,   1934.136    , 0.00478  ],
249:        [SIN   ,  200.93      ,  72001.539    , 0.00037  ]]) / 360
250:   end

     # File when/ephemeris.rb, line 223
223:   def julian_century_from_2000(tt)
224:     return (tt - EPOCH4) / JCENT
225:   end

     # File when/ephemeris.rb, line 215
215:   def julian_year_from_1975(tt)
216:     return (tt - EPOCH3) / JYEAR
217:   end

     # File when/ephemeris.rb, line 259
259:   def obl(c)
260:     return (23.43929 + 0.013004*c + 1.0*delta_e(c)) / 360
261:   end

    # File when/ephemeris.rb, line 70
70:   def polynomial(t, equ)
71:     equ.reverse.inject(0) {|sum, v| sum * t + v}
72:   end

     # File when/ephemeris.rb, line 271
271:   def root(t0, y0=nil, &func)
272: 
273:     # 近似値0,1
274:     # printf("y0=%20.7f\n",y0)
275:     t = [t0,              t0+0.1        ]
276:     y = [func.call(t[0]), func.call(t[1])]
277:     # printf("t=%20.7f,L=%20.7f\n",t[1],y[1])
278: 
279:     # 近似値2(1次関数による近似)
280:     t << (y0 ? (t[1]-t[0])/(y[1]-y[0])*(y0-y[0])+t[0] : t0-0.1)
281: 
282:     # 繰り返し
283:     i = 10
284:     while ((t[2]-t[1]).abs > 1E-6) && (i > 0)
285:       # 予備計算
286:       y << func.call(t[2])
287:       break if y0 && (y[2]-y0).abs <= 1E-7
288: 
289:       # printf("t=%20.7f,L=%20.7f\n",t[2],y[2])
290:       t01     = t[0]-t[1]
291:       t02,y02 = t[0]-t[2], y[0]-y[2]
292:       t12,y12 = t[1]-t[2], y[1]-y[2]
293: 
294:       # 2次関数の係数
295:       a = ( y02     / t02 - y12     / t12) / t01
296:       b = (-y02*t12 / t02 + y12*t02 / t12) / t01
297:       c = y[2]
298: 
299:       if y0
300:         # 判別式
301:         if (d = b*b-4*a*(c-y0)) < 0
302:           i = 1
303:           break
304:         end
305: 
306:         # 近似値(2次関数による近似)
307:         sqrtd = Math.sqrt(d)
308:         sqrtd = -sqrtd if b < 0
309:         t << (t[2] + 2*(y0-c)/(b+sqrtd)) # <-桁落ち回避
310:       else
311:         t << (t[2] - b / (2*a))
312:       end
313: 
314:       t.shift
315:       y.shift
316:     end
317:     $stderr.puts "The result is not reliable" if i<=0
318:     return t[2]
319:   end

     # File when/ephemeris.rb, line 134
134:   def sinc(x)
135:     sin(x * CIRCLE)
136:   end

     # File when/ephemeris.rb, line 119
119:   def sind(x)
120:     sin(x * DEG)
121:   end

     # File when/ephemeris.rb, line 139
139:   def tanc(x)
140:     tan(x * CIRCLE)
141:   end

     # File when/ephemeris.rb, line 124
124:   def tand(x)
125:     tan(x * DEG)
126:   end

     # File when/ephemeris.rb, line 86
 86:   def trigonometric(t, equ, dl=0.0, count=0)
 87:     t2 = t * t
 88:     equ[0..(count-1)].inject(0) do |sum, v|
 89:       op, epoch, freq, amp, sqr = v
 90:       ds = epoch + t * freq
 91:       if (op < 0) # 直線
 92:         ds += t2 * amp
 93:       else        # 三角関数
 94:         ds += t2 * sqr if sqr
 95:         ds += dl  if (op[2]==1) # delta L is exist
 96:         ds  = amp * ((op[0]==1) ? cosd(ds) : sind(ds))
 97:         ds *= t   if (op[1]==1) # time proportional
 98:       end
 99:       sum += ds
100:     end
101:   end