アレアレ

お役立ち情報、お悩み解決情報を発信!

xcode

Swiftで天体の位置計算

Swiftで天体の位置計算

今回は、天体の位置をプログラミングで計算してみたい人へのお役立ち情報として、Swiftを使ってその計算を行う方法をご紹介します。

天体の位置計算をわかりやすく解説した名著として、長沢工の「天体の位置計算」という本があります。で、その本で解説されている、「赤道座標で示された赤緯・赤経の恒星が、地平座標上のどの方位・高度に見えるか」という変換処理を「プログラミング言語Swiftを使って、こう実装しました」というのが、今回ご紹介したい内容です。

動作環境

はじめに、今回ご紹介する方法の動作確認が取れている環境をご紹介すると、次の通りです。

  • OS: OS X El Capitan (Version 10.11.5)
  • Xcode: 7.3.1 (7D1014)
  • swift: 2.2 (swiftlang-703.0.18.8 clang-703.0.31)

開発環境には、Xcodeを使います。今回ご紹介のプログラムは、そのXcodeのPlaygroundと呼ばれる、Swiftのプログラムを簡単に実行できる環境で実行できるものです。

「天体の位置計算」の古典的な手法をプログラミング化

では、本題に入ります。まず今回ご紹介する天体の位置計算方法ですが、上述した書籍「天体の位置計算」内の「2.6 赤道座標から地平座標への変換」までで解説された内容をプログラミング化したものです。

で、この手法は、この書籍内でも「どちらかと言えば古典的なやり方」と解説されているのですが、そこをあえてプログラミング化しています。

まず今回ご紹介するプログラムですが、「StarLocator」というクラスが中心です。このクラスに、観測点の緯度・経度、日時を渡して初期化し、観測したい恒星の赤緯・赤経を登録して利用します。

具体的なコードとしては次のようになります。

//1978年6月20日22時32分17秒(日本標準時間)
//PZTの緯度: 35度40'20".707N 経度: 9h18m9s.936E
let locator = StarLocator(
    time: (year: 1978, month: 6, day: 20, h: 22, m: 32, s: 20)
    ,latitude: (d:35, m:40, s:20.707, orientaion:Orientation.North)
    ,longitude: (h:9, m:18, s:9.936, orientaion:Orientation.East)
)

//Aur Capellaの赤緯:45度56'58".04 赤経:5h12m59s.466
locator.addStar("Aur Capella", declination: (d: 45, m: 56, s: 58.04), rightAscension: (h: 5, m: 12, s: 59.466))
locator.addStar("61 Cyg", declination: (d: 38, m: 29, s: 59.10), rightAscension: (h: 21, m: 4, s: 39.935))

var theta = locator.localSiderealTimeDecimal
let phai = locator.latitudeDecimal

//変換結果を出力する
for star in locator.stars{
    
    let alfa = star.rightAscension.timeToDegrees()
    let delta = star.declination.timeToDegrees()
    
    let hA = locator.calcHighangleAndAzimuth(phai, theta: theta, alfa: alfa, delta: delta)
    
    NSLog("\(star.name): h=\(hA.h) A=\(hA.A)")
}

で、以上のプログラムですが、真っさらな「Playground」に貼り付けても、動作しません。「StarLocator」クラスの定義がないから当然そうなります。

そこで、以上のコードのに、次のコードをコピペしてください。「StarLocator」クラスも含めて、必要なコードの定義となっています。

import Foundation
import UIKit

//赤緯の点
public class Declination{
    var degrees: Int
    var minites: Int
    var seconds: Double
    
    public init(d:Int, m:Int, s:Double){
        self.degrees = d
        self.minites = m
        self.seconds = s
    }
    
    //時間の赤緯表記から角度に変換する
    public func timeToDegrees() -> Double{
        return Double(degrees) + Double(minites)/60 + seconds/3600
    }
}

//赤経の点
public class RightAscension{
    var hours: Int
    var minites: Int
    var seconds: Double
    
    let anHourToDegree: Double = 15.0
    let anMinutesToDegree: Double = 1.0 / 4
    let anSecondToDegree: Double = 15.0 / 3600
    
    public init(h:Int, m:Int, s:Double){
        self.hours = h
        self.minites = m
        self.seconds = s
    }
    
    //時間の赤経から角度に変換する
    public func timeToDegrees() -> Double{
        return Double(hours)*anHourToDegree + Double(minites)*anMinutesToDegree + seconds*anSecondToDegree
    }
}

//星の座標
public class Star{
    public var name: String
    public var declination: Declination
    public var rightAscension : RightAscension
    
    public init(name:String, declination: (d:Int, m:Int, s:Double), rightAscension : (h:Int, m:Int, s:Double))
    {
        self.name = name
        self.declination = Declination(d: declination.d, m: declination.m, s: declination.s)
        self.rightAscension = RightAscension(h: rightAscension.h, m: rightAscension.m, s: rightAscension.s)
    }
}

//GMT恒星時計算クラス
public class SiderealTime{
    public let baseYear:Int = 1899
    public let baseMonth:Int = 12
    public let baseDay:Int = 31
    public let baseHour:Int = 12
    
    
    public var year:Int
    public var month:Int
    public var day:Int
    public var hour:Int
    public var minute:Int
    public var second:Double
    
    public init(year:Int, month:Int, day:Int, hour:Int, minute:Int, second:Double){
        self.year = year
        self.month = month
        self.day = day
        self.hour = hour
        self.minute = minute
        self.second = second
    }
    
    public func diffYears()->Int{
        let years = year - baseYear - 1
        return years
    }
    
    public func uruuYearsCount()->Int{
        var count = 0
        
        for y in baseYear ..< year+1{ if isUruuYear(y) { count = count + 1 } } return count } private func isUruuYear(year:Int)->Bool{
        return (year % 4 ) == 0 &amp;&amp; ( year % 100 ) != 0 || ( year % 400 ) == 0
    }
    
    public func daysOfTheYear()->Int{
        let daysCount:[Int] = [31,28,31,30,31,30,31,31,30,31,30,31]
        
        var days = 0
        
        for m in 1 ..< self.month{ days = days + daysCount[m] if m==2 &amp;&amp; self.isUruuYear(self.year) &amp;&amp; self.month > 2 {
                days = days + 1
            }
        }
        
        return days + self.day
    }
    
    public func diffDays()->Double{
        return Double(diffYears()*365 + uruuYearsCount() + daysOfTheYear()) + Double(self.baseHour) / 24
    }
    
    public func Tu()->Double{
        return diffDays() / 36525
    }
    
    //世界時0時のグリニジ平均恒星時Θ(0)の第2項
    public func theta2ndSeconds()->Double{
        return 8640184.542 * Tu()
    }
    
    //世界時0時のグリニジ平均恒星時Θ(0)の第3項
    public func theta3ndSeconds()->Double{
        return 0.0929 * pow(Tu(),2)
    }
    
    //世界時0時のグリニジ平均恒星時Θ(0)
    public func thetaSeconds()->Double{
        return  timeToSeconds(0, hour: 6, minute: 38, second: 45.836) + theta2ndSeconds() + theta3ndSeconds()
    }
    
    //秒を日/時/分/秒形式に変換する
    public func secondsToTime(seconds:Double)->(day: Int, hour: Int, minute: Int, second: Double){
        let day = seconds / (60 * 60 * 24)
        let hour = (day - Double(Int(day))) * 24
        let minute = (hour - Double(Int(hour))) * 60
        let second = (minute - Double(Int(minute))) * 60
        
        return (Int(day), Int(hour), Int(minute), second)
    }
    
    //日/時/分/秒形式を秒に変換する
    public func timeToSeconds(day: Int, hour: Int, minute: Int, second: Double)->Double{
        return Double(day * 24 * 60 * 60 + hour * 60 * 60 + minute * 60) + second
    }
    
    //参考:経過日数計算の簡易版
    public func calcK(year:Int, month:Int, day:Int) -> Double{
        var y = Double(year)
        var m = Double(month)
        let d = Double(day)
        
        if m==1{
            y = y - 1
            m = m + 1
            
        }else if m==2{
            y = y - 1
            m = m + 2
        }
        
        var temp = 365*y
        temp = temp + 30*m
        temp = temp + d
        temp = temp - 33.5
        temp = temp + Double(Int(3.0/5.0 * (m+1)))
        temp = temp + Double(Int(y/4.0))
        return temp
    }
}

//方角を表す列挙
public enum Orientation{
    case South
    case East
    case North
    case West
}

//指定した観測地点で計測したい恒星がどの方位高度に見えるかを計算するクラス
public class StarLocator{
    
    //観測時間(現地時間)
    public var time: (year:Int, month:Int, day:Int, h:Int, m:Int, s:Double)
    //観測点の緯度(度 分 秒)
    public var latitude: (d:Int, m:Int, s:Double, orientaion:Orientation)
    //観測点の緯度(小数)
    public var latitudeDecimal:Double{
        get{
            var temp: Double = Double(latitude.d)
            temp = temp + Double(latitude.m)/60.0
            temp = temp + (latitude.s / (60 * 60))
            
            if longitude.orientaion == Orientation.South {
                temp = -1 * temp
            }
            
            return temp
        }
    }
    //観測点の経度(時 分 秒)
    public var longitude: (h:Int, m:Int, s:Double, orientaion:Orientation)
    //観測点の経度(小数)
    public var longitudeDecimal:Double{
        get{
            var temp: Double = Double(longitude.h) * 15.0
            temp = temp + Double(longitude.m)/60.0 * 15.0
            temp = temp + (longitude.s / (60 * 60)) * 15.0
            
            if longitude.orientaion == Orientation.South {
                temp = -1 * temp
            }
            
            return temp
        }
    }
    //観測時間のGMTとの差(単位:時)
    public var gmtDiff: Double{
        get{
            var diff = Double(longitude.h)/24.0 * 360 / 15.0
            
            if longitude.orientaion==Orientation.West{
                diff = -1*diff
            }
            
            return diff
        }
    }
    //観測対象の恒星
    public var stars:[Star] = []
    //観測対象の恒星を追加する
    public func addStar(name:String, declination:(d:Int, m:Int, s:Double), rightAscension: (h:Int, m:Int, s:Double)){
        let star = Star(name: name, declination:declination, rightAscension: rightAscension)
        
        self.stars.append(star)
    }
    //世界0時のグリニジ平均恒星時(秒)
    public var gmtSiderealTimeSeconds:Double{
        get{
            let sd = SiderealTime(year: time.year, month: time.month, day: time.day, hour: 0, minute: 0, second: 0)
            return sd.thetaSeconds()
        }
    }
    //世界0時のグリニジ平均恒星時(日時分秒)
    public var gmtSiderealTime:(day: Int, hour: Int, minute: Int, second: Double){
        get{
            let sd = SiderealTime(year: time.year, month: time.month, day: time.day, hour: 0, minute: 0, second: 0)
            return sd.secondsToTime(sd.thetaSeconds())
        }
    }
    //観測点の恒星時(日時分秒)
    public var localSiderealTime:(hour: Int, minute: Int, second: Double){
        get{
            
            let gmt = gmtSiderealTime
            
            //経過時間=観測時刻-GMTとの差の時間
            let keikaSeconds:Double = Double(Double(time.h) - gmtDiff)*60*60 + Double(time.m)*60 + Double(time.s)
            
            //補正値
            let hoseiichiSeconds = keikaSeconds * 0.00273791
            
            //経度の時間変換(秒)
            var pztLatLanSeconds = Double(longitude.h)*60*60 + Double(longitude.m)*60 + longitude.s
            
            if longitude.orientaion == Orientation.West {
                pztLatLanSeconds = -1 * pztLatLanSeconds
            }
            
            //グリニジ視恒星時の秒
            let gKouseijiSeconds =  Double(gmt.hour*60*60 + gmt.minute*60) + gmt.second
            
            //合計(秒)
            let sum = keikaSeconds + hoseiichiSeconds + pztLatLanSeconds + gKouseijiSeconds
            let hours = sum / (60*60)
            let hoursInteger = Int(hours) % 24
            let minites = (hours - Double(Int(hours))) * 60
            let seconds = (minites - Double(Int(minites))) * 60
            
            return (Int(hoursInteger), Int(minites), seconds)
        }
    }
    //観測点の恒星時(秒)
    public var localSiderealTimeDecimal:Double{
        get{
            let lst = localSiderealTime
            
            var temp: Double = Double(lst.hour) * 15.0
            temp = temp + Double(lst.minute)/60.0 * 15.0
            temp = temp + (lst.second / (60 * 60)) * 15.0
            
            return temp
            
        }
    }
    //コンストラクタ
    public init(time: (year:Int, month:Int, day:Int, h:Int, m:Int, s:Double)
        , latitude: (d:Int, m:Int, s:Double, orientaion:Orientation)
        , longitude: (h:Int, m:Int, s:Double, orientaion:Orientation))
    {
        self.time = time
        self.latitude = latitude
        self.longitude = longitude
    }
    
    //度をラジアンに変換する
    public func toRadian(degree:Double)->Double{
        return degree * M_PI / 180
    }
    
    //ラジアンを度に変換する
    public func toDegree(radian:Double)->Double{
        return radian / (M_PI / 180)
    }
    
    //観測地点の緯度、恒星時、赤経、赤緯を高度と方位角に変換する関数
    public func calcHighangleAndAzimuth(latitude:Double, theta:Double, alfa:Double, delta:Double)->(h:Double, A:Double){
        let f1 = theta - alfa                       //Θ-α           (1)
        let f2 = cos(f1 * M_PI / 180)               //cos(Θ-α)      (2)
        let f3 = sin(f1 * M_PI / 180)               //sin(Θ-α)      (3)
        let f4 = cos(latitude * M_PI / 180)         //cos(φ)        (4)
        let f5 = sin(latitude * M_PI / 180)         //sin(φ)        (5)
        let f6 = cos(delta * M_PI / 180)            //cos(δ)        (6)
        let f7 = sin(delta * M_PI / 180)            //sin(δ)        (7)
        let f8 = f5 * f7                            //sin(φ)sin(δ)  (8)=(5)x(7)
        let f9 = f4*f6*f2                           //cos(φ)cos(δ)cos(Θ-α)  (9)=(4)x(6)x(2)
        let f10 = f8+f9                             //sin(h)        (10)=(8)+(9)
        let f11 = asin(f10) / (M_PI / 180)          //h             (11)
        let f12 = -1 * f4 * f7                      //-cos(φ)sin(δ) (12)=-(4)x(7)
        let f13 = f5 * f6 * f2                      //sin(φ)cos(δ)cos(Θ-α)  (13)=(5)x(6)x(2)
        let f14 = f12+f13                           //cos(h)cos(A)  (14)=(12)+(13) ※ここの符号で象限が決まる
        let f15 = f6 * f3                           //cos(h)sin(A)  (15)=(6)x(3)
        let f16 = f15 / f14                         //tan(A)        (16)=(15)÷(14)
        let f17 = atan(f16) / (M_PI / 180)          //A'            (17)
        
        var f18 = f17                               //A=A'
        
        //cos(h)cos(A)の符号で象限の調整
        if(f14<0) {
            f18 = f18 + 180.0                       //(18)=(17)+180度
        }
        
        //Aは正符号で表記する
        if f18 < 0 { f18 = f18 + 360 } return (f11, f18) } //時間表記の角度を度表記に変換する public func hourToDegree(hour:(h:Int, m:Int, s:Double))->(d:Int, m:Int, s:Double){
        
        let totalSeconds:Double = (Double(hour.h * 60 * 60 + hour.m * 60) + hour.s) * 15
        
        let d = totalSeconds / (60 * 60)
        let m = (d - Double(Int(d))) * 60
        let s = (m - Double(Int(m))) * 60
        
        return (Int(d), Int(m), s)
    }
    //度表記の角度をに時間表記に変換する
    public func degreeToHour(degree:(d:Int, m:Int, s:Double))->(h:Int, m:Int, s:Double){
        let h: Double = Double(degree.d)/15
        let m: Double = (h - Double(Int(h))) * 60 + (Double(degree.m) / 15)
        let s: Double = (m - Double(Int(m))) * 60 + degree.s/15.0
        
        return (Int(h), Int(m), s)
    }
}	

以上のコードがPlayground上でうまく動作すると、次のような標準出力がPlaygroundの下の方に得られるはずです。

Aur Capella: h=-8.11912573218584 A=174.938785277314
61 Cyg: h=39.0758820971646 A=246.157497408152

ss

この出力のhが高度、Aが方位を表しています。つまり、Capellaは、1978年6月20日22時32分17秒(日本標準時間)の指定の緯度経度の観測点では、高度hがマイナスなので、見られないことがわかります。

一方、61 Cygは、高度hが約39度、方位が約246度の方向に見られることがわかるわけです。

Return Top