しばらく続いていた多忙な日々も収束してきて、ぼちぼち諸々の活動を再開しようかと思っています。

今回手始めのネタは多忙中に時間を盗んでいろいろ仕込んでいたもののひとつで、自身の備忘録を兼ねています。
一部未解決で、通りがっかった方でこのあたり詳しい人がおられれば是非教えてください。

まずは至極当然な前提から。

撮影にとって最も重要な要素は何か？

言うまでもなくそれは「光」であって、特に自然界での撮影だとその源は太陽であり、あるいは場合によって月になる。

というわけで、太陽や月の方向というのは屋外で撮影する限りは常に注意を払わなければならない…当然ですが。

というわけで、そのあたりを机上でも定量的に検討できるよう天文学の素人でも計算できるようなライブラリを以前から物色しておりました。

本格的なものだと、その筋のデファクト（らしい）SOFAとか

http://www.iausofa.org/

NASAご謹製のやつとか

https://naif.jpl.nasa.gov/naif/index.html

しかしながらこのあたり、どう見ても素人が手軽に使うには敷居が高すぎる。

そこで目をつけたのが以下の2つのライブラリ。

PyEphem

libnova

前者はPython系、後者はC++系でこのあたりだと素人でもなんとか扱えそう。
ただし使いこなすには最低限の天文学の基礎知識は必須だが、以下が日本語記述かつ必要十分。

https://www.astroarts.co.jp/alacarte/kiso/index-j.shtml

ということで双方のライブラリを用いて「その日・その時・その場所」の太陽・月の地平座標及び月の満ち欠けを計算するスニペットを作ってみました。

まずはPyEphemから。実行にはPyEphemライブラリが必要。

# sun_moon.py

import ephem
import datetime
import sys
from math import radians as rad,degrees as deg  
from pytz import timezone

args = sys.argv

if len(args) < 9:
	print("Usage: sun_moon.py year mon date hour min sec longitude latitude")
	sys.exit()
	
# arg to parameter
year=int(args[1])
month=int(args[2])
day=int(args[3])
hour=int(args[4])
min=int(args[5])
sec=int(args[6])
longitude = float(args[7])
latitude = float(args[8])

# UTC conversion
jst = datetime.datetime(year, month, day, hour, min, sec)
local_tz = timezone('Asia/Tokyo')
jst = local_tz.localize(jst)
print("JST ",jst)
utc = jst.astimezone(timezone('UTC'))
print("UTC ",utc)

ephem_tm = ephem.Date(utc)
print ("Ephem Time ", ephem_tm)

gatech = ephem.Observer()
gatech.lon = rad(longitude)
gatech.lat = rad(latitude)
gatech.date = utc

sun = ephem.Sun()
moon = ephem.Moon()
sun.compute(gatech)
moon.compute(gatech)

print("Solar AZ: %s ALT: %s" % (deg(sun.az),deg(sun.alt)))
print("Lunar AZ: %s ALT: %s" % (deg(moon.az),deg(moon.alt)))

nnm = ephem.next_new_moon(gatech.date)
pnm = ephem.previous_new_moon(gatech.date)

lunation=(gatech.date-pnm)/(nnm-pnm)  
print("Lunar lunation: ",lunation)

起動方法は

python sun_moon.py year month date hour min sec 経度緯度

次にlibnovaバージョン。コンパイル・実行にはlibnovaライブラリが必要。

// sun_moon.cpp

#include <iostream>
#include <string>
#include <string.h>

#include <libnova/utility.h>
#include <libnova/julian_day.h>
#include <libnova/angular_separation.h>
#include <libnova/lunar.h>
#include <libnova/solar.h>
#include <libnova/transform.h>

// S:0 W:90 N:180 E:270 -> N:0 E:90 S:180 W:270
double ConvAz(double pAZorg) {
    double ret = pAZorg - 180.; 
    return ((ret > 0.) ? ret : ret + 360.);
}

int main(int argc,char *argv[]) {
    if(argc < 9) {
        std::cout << "sun_moon year mon date hour min sec longitude latitude" << std::endl;
        return -1;
    }

    // arg conversion
    std::string s_year(argv[1]);
    std::string s_mon (argv[2]);
    std::string s_day (argv[3]);
    std::string s_hour(argv[4]);
    std::string s_min (argv[5]);
    std::string s_sec (argv[6]);
    std::string s_lon (argv[7]);
    std::string s_lat (argv[8]);
    
    // Get Julian day
    time_t time_of_the_day;
    struct tm ptm;
    double jd; 
    
    memset(&ptm, 0, sizeof(struct tm));
    ptm.tm_year = std::stoi(s_year) - 1900;
    ptm.tm_mon  = std::stoi(s_mon) - 1;
    ptm.tm_mday = std::stoi(s_day);
    ptm.tm_hour = std::stoi(s_hour);
    ptm.tm_min  = std::stoi(s_min);
    ptm.tm_sec  = std::stoi(s_sec);
    time_of_the_day = mktime(&ptm);
    
    std::cout << "Local time " << (ptm.tm_year+1900) << "/" << (ptm.tm_mon + 1) << "/" << ptm.tm_mday << 
        " " << ptm.tm_hour << ":" << ptm.tm_min << ":" << ptm.tm_sec << std::endl;
    jd = ln_get_julian_from_timet(&time_of_the_day);
    
    std::cout.setf(std::ios::fixed);    
    std::cout << "Julian day  " << jd << std::endl;

    // Get equ position of sun and moon.
    struct ln_equ_posn lunar_equ, solar_equ;
    ln_get_solar_equ_coords(jd, &solar_equ);
    ln_get_lunar_equ_coords(jd, &lunar_equ);
    
    ln_lnlat_posn observer;
    ln_hrz_posn position;
    
    observer.lng = std::stof(s_lon);
    observer.lat = std::stof(s_lat);
    
    // Get Horizontal coord
    ln_get_hrz_from_equ(&solar_equ, &observer,jd,&position);
    std::cout << "Solar AZ:" << ConvAz(position.az) << " ALT:" << position.alt << std::endl;

    ln_get_hrz_from_equ(&lunar_equ, &observer,jd,&position);
    std::cout << "Lunar AZ:" << ConvAz(position.az) << " ALT:" << position.alt << std::endl;
    
    // Get Lunar Phase
    double phase = ln_get_lunar_phase(jd);
    std::cout << "Lunar phase " << phase << std::endl;

    return 0;
}

コンパイル後、実行は

sun_moon.exe year month date hour min sec 経度緯度

さて、これらの実地検証を…

こんなときに役に立つのが対象の方角と仰角を写せるスマホ用の測量アプリ。
自分はAndroid用のSurveyCompassというのが精度がよいので使ってますが、荒川土手で太陽を撮ったのがこんな感じ。

撮ったのが2017/2/4 15:27なので東京の緯度・経度と共に両者を起動してみる。AZ=Azimuth(方位角) ALT=Altitude(高度)

# PyEphem

>python sun_moon.py  2017 2 4 15 27 0 139 35
JST  2017-02-04 15:27:00+09:00
UTC  2017-02-04 06:27:00+00:00
Ephem Time  2017/2/4 06:27:00
Solar AZ: 233.3656526924783 ALT: 18.823892069528423
Lunar AZ: 116.80623420416293 ALT: 49.93173592735221
Lunar lunation:  0.2452435712394331

// libnova

>./sun_moon.exe 2017 2 4 15 27 0 139 35
Local time 2017/2/4 15:27:0
Julian day  2457788.768750
Solar AZ:233.488559 ALT:18.575095
Lunar AZ:117.142069 ALT:50.684293
Lunar phase 88.705032

若干の誤差はあるが、まあいい感じです。

ほぼ同じ時間に撮った月がこんな感じ。

何故か月の場合は経度方向に10度以上も誤差が…この時間に限らず常に同じくらいのズレが継続して出ます。原因はまだ分かってません。

何か勘違いをしているのか…通りがっかった方でこのあたり詳しい人がおられれば是非教えてください。

ちなみに、月の満ち欠けの値ですが、PyEphemは0～1の間を遷移、0と1が新月、0.5が満月となります。

libnovaは0～180を遷移、0が満月、180が新月となります。
こちらは正しい値が出るようです。

光あれ～PyEphem vs libnova