１　周期を求める

　観測データを以下のフォーマットで用意する。1カラム目は観測時間で単位は日。ドキュメントではJD, MJDとなっているが、ここでは観測初日のUT0時をゼロとしている。2カラム目は等級、3カラム目は等級の誤差である。ここでは mag.dat というファイルに記録してある。

 0.64220         1.783    0.023
 0.64368         1.840    0.025
 0.64516         1.733    0.024
 0.64662         1.833    0.024
 0.64810         1.814    0.023
 0.64958         1.787    0.022

　cyclocode を起動し、1 を入力して Periodogarm を実行する。

　観測データのファイル名を入力する。ここでは mag.dat とする。

　周期解析の方法を選択する。ここでは省略値の Spectral Analysis を選ぶので、そのままリターンキーを押す。

　２つのパラメータ ofac と hifac を入力する。これらについてはこの後詳しく説明する。

　次に、ライトカーブがシングルピークかダブルピークか聞いてくるが、一般的な小惑星はダブルピークのため、そのままリターンキーを押す。

　計算結果として、周期（単位は日）、有意水準が出力される。

　※ ofac に関しては調べても詳しくわからなかった。以下の記述は筆者自身の考察である。

　ofac は oversampling parameter で "4" 以上を入力するとしており、算出する周波数の下限（短波長側）を決定するパラメータと説明されている。だが実態として、算出する周波数の間隔を決めるパラメータでもある。算出する周波数は、ゼロから最短周波数の値ごとのとびとびの値となるからである。具体的には *.spc ファイルの1カラム目の周波数値である。

　ofac = 4 の場合、求まる周波数の間隔（＝最短周波数）は 1/8T となる。ここで T は観測の期間（単位は日）である。

　例えば、観測が 2 夜と短い場合、仮に観測期間を T = 1.2 とすると、ofac = 4 では周波数の間隔は 0.104 なる。周期 7 時間と仮定すると周期の間隔が 10 分程度にもなる。そこで ofac = 40 とすると、周期が約 1 分間隔で得られる。

　ofac=4 と ofac=40 とで得られた Spectral Analysis の結果のグラフを以下に示す。ofac = 4 では周波数の算出間隔が荒いため、ofac = 40 では Power が最高値となった 0.3011 という周期が算出できていない。ofac = 4 ではその箇所のピークの形がいびつとなり 0.3052 という異なった値が得られている。また時間間隔が荒いことからグラフ全体が滑らかでない。

　このように、特に観測期間が短い場合や周期が長い場合は、ofac の値を適切に与えて周波数間隔が十分に小さくなるよう注意する必要がある。実際には、まず ofac = 4 で周期を求めてみて、*.spc ファイルの結果から求める周期の時間間隔を調べ、必要に応じ ofac を増やして再度実行するのが良いだろう。

　hifac = 1 が妥当である。

　hifac は求める周波数の上限（長周期側）を決定するパラメータである。FFT で求められる周波数は、原理的にサンプリング周波数の 1/2 までに限られる。この 1/2 の周波数をナイキスト周波数 (fc) という。求めたい周波数の上限を fh とした場合、hifac = fh / fc の関係となるため、hifac は 1 またはそれ以下でなければならない。ただし、極端に長い周期は小惑星の自転周期として採用されないため神経質になる必要はないだろう。

　また、結果として得られる周波数の個数 Np（*.spc ファイルに出力される周波数の数）は ofac と hifac の値に依存し、次のとおりである。

　　　　

　ここでは、gnuplot を用いて *.spc ファイル Spectral Analysis の結果のグラフを描く例を示す。

　まず、*.spc ファイルのヘッダ部分、先頭4行を削除する必要があるので、オリジナルファイルの複製を作り、それに加工を加える。ここでは、加工した複製ファイルを mag_gnuplot.spc とする。

　gnuplot を起動し、次のコマンドを入力すると下のグラフが得られる。

set key off
set log x
set xlabel "Period (d)"
set ylabel "Power"
plot [0.01:10][0:60] "mag_gnuplot.spc" using 2:3 with lines

　このグラフでは 0.2616 と 0.3011 にピークがあり、わずかに 0.3011 の方が Power が大きいが、単純に Power が大きい方を正確な周期とするわけではない。このあとライトカーブを描き判断する。

２　ライトカーブを描く

　cyclocode を起動し、 2 を入力して Synthetic Lightcurve を実行する。

　観測データのファイル名を入力する。ここでは mag.dat とする。

　Base-level of the mag. data is zero mag. [y/Y (default) or n/N]: には n を入力する。
　次に、1 の Periodogram で得られた周期を入力する。ここでは、先ほど得られた 0.2616 と 0.3011 のうち、0.3011 を入力している。

カーブフィッティングの次数の最小値と最大値を入力すると、それぞれの次数のベース値とχ 2   の値

(Chi-Square) が出力される。

　フィッティングの次数は、原理的には Chi-Square の値が最小となる次数が最適だが、現実には高次すぎると誤差も含めて無理にフィットしてしまうので、必要最小限の次数とすべきである。

　ベース値は、ライトカーブの振幅の中心値で単位は等級。ベース値を用いると、異なる比較星を用いた観測であっても、それぞれのデータで周期解析ができた場合、ベース値の差だけシフトすることで簡単にライトカーブを合成できる。

　実行後、以下のファイルが出力される。

*.pha　Phase の昇順にソートされている。JD, mag, err はオリジナルデータの値。normag はカーブフィッティングの推定値で、値はベース値からの相対等級である。 1 行目の数値は、各次数に対応したベース値。
　このファイルは、ライトカーブとともに観測データをプロットする場合に用いる。（オリジナルドキュメント P6 サンプル中段図参照）

                                        1.7083    1.7071    1.7084    1.7084
  Phase      JD       mag      err    3:normag  4:normag  5:normag  6:normag
  0.0078   0.7868    1.763    0.018     0.055     0.056     0.055     0.055 
  0.0135   0.7883    1.753    0.023     0.045     0.046     0.045     0.045 
  0.0191   0.7898    1.750    0.021     0.042     0.043     0.042     0.042 
  0.0248   0.7913    1.785    0.022     0.077     0.078     0.077     0.077 
  0.0305   0.7928    1.780    0.022     0.072     0.073     0.072     0.072 

*.roc　Phase 0 〜 1 を 分割して 0.0125 毎にカーブフィッティングの推定値を並べたファイル。mag, nomag どちらも推定値、nomag はベース値からの相対等級。
　このファイルは、1 周期分のライトカーブそのものを描くために用いる。（オリジナルドキュメント P6 サンプル中段図参照）

              3      1.7083     4      1.7071     5      1.7084     6      1.7084   
  Phase      mag    normag     mag    normag     mag    normag     mag    normag 
  0.0000    1.774    0.065    1.779    0.072    1.760    0.052    1.756    0.048 
  0.0125    1.771    0.062    1.777    0.070    1.759    0.050    1.758    0.049 
  0.0250    1.766    0.057    1.772    0.065    1.757    0.049    1.759    0.050 
  0.0375    1.759    0.050    1.764    0.057    1.756    0.047    1.759    0.051 
  0.0500    1.750    0.042    1.754    0.047    1.753    0.045    1.757    0.049 
 

*.fuc　Phase を 80 分割し 0.0125 毎に推定値を並べているのは *.roc と同じだが、こちらは pahse 0 〜 1 で折り返さず、最初の観測から最後の観測までの期間を通して計算している。このため、1 カラム目に JD があり、JD の昇順にソートされている。計算の期間は前後にそれぞれ 1/4 周期の余裕をもっている。
　このファイルは、観測期間全体のライトカーブそのものを描くために用いる。（オリジナルドキュメント P6 サンプル下図参照）

                       3      1.7083     4      1.7071     5      1.7084     6      1.7084   
    JD     Phase      mag    normag     mag    normag     mag    normag     mag    normag 
  0.5768   0.2049    1.650   -0.058    1.648   -0.059    1.618   -0.091    1.624   -0.084 
  0.5801   0.2174    1.652   -0.056    1.652   -0.055    1.621   -0.088    1.629   -0.080 
  0.5833   0.2299    1.656   -0.052    1.659   -0.048    1.629   -0.079    1.637   -0.071 
  0.5866   0.2424    1.663   -0.046    1.667   -0.040    1.643   -0.066    1.650   -0.058 
  0.5899   0.2549    1.671   -0.038    1.677   -0.030    1.661   -0.048    1.666   -0.043 

*.fud　JD はオリジナルデータで、JD の昇順にソートされている。Phase はそれぞれの JD に対応した値。mag, nomag は推定値。
　このファイルは、ライトカーブとともに推定値をプロットする場合に用いる。（オリジナルドキュメント P6 サンプル下図参照）

                       3      1.7083     4      1.7071     5      1.7084     6      1.7084   
    JD     Phase      mag    normag     mag    normag     mag    normag     mag    normag 
  0.6422   0.4549    1.803    0.095    1.797    0.090    1.787    0.079    1.788    0.079 
  0.6437   0.4606    1.802    0.094    1.797    0.090    1.789    0.081    1.789    0.080 
  0.6452   0.4662    1.801    0.092    1.797    0.090    1.791    0.083    1.790    0.082 
  0.6466   0.4718    1.799    0.091    1.796    0.089    1.793    0.085    1.791    0.083 
  0.6481   0.4774    1.797    0.089    1.796    0.089    1.795    0.086    1.792    0.084 

　ここでは、gnuplot を用いて 1 周期分のライトカーブを描く例を示す。観測日毎にプロットの色を変えたり、Phase の原点を変えたり、すべての観測期間に渡るライトカーブを描くなど、他のケースは各自工夫されたい。

　まず、*.pha ファイルと *.roc ファイルのヘッダ部分、それぞれ先頭2行を削除する必要があるので、オリジナルファイルの複製を作り、それに加工を加える。ここでは、加工した複製ファイルをそれぞれ mag_gnuplot.pha、mag_gnuplot.roc とする。

　gnuplot を起動し、次のコマンドを入力すると下のグラフが得られる。

set key off
set xlabel "Phase"
set ylabel "Relative Magnitude (R)"
plot [0:1][1.9:1.5] "mag_gnuplot.pha" using 1:3:4 with err, "mag_gnuplot.roc" using 1:2 with lines

３　周期の誤差を求める

　cyclocode には周期の誤差を求める機能はない。周期の誤差は以下の式から求める。

(1)

　ここで f は周期の周波数。cyclocode は周期 (P) を日の単位で出力するので、f = 1 / P から求める。

　 T は観測の期間（単位は日）、N は観測数、A は変光幅(amplitude)で *.roc ファイルの mag の値の最大値と最小値との差から求める。

σ 2   はデータの誤差の分散である。

　ここでデータの誤差（Error）とは、カーブフィッティングの誤差（Δfit）と測光誤差（Δmag）との合成である。

　カーブフィッティングの誤差（Δfit）は、観測値（ *.pha ファイルの mag の値）とカーブフィッティングによる推定値（ *.fud ファイルの mag の値）との差である。

　測光誤差（Δmag）は、*.pha ファイルの err の値である。

　これらを次式により合成しデータの誤差（Error）を得る。

　　　　　　　　　　　　　(2)

得られたすべての Error の分散が σ 2   となる。

　式（１）の左辺は、周期の誤差の割合である。この値に周期(P)を乗じると、周期に対する誤差の値が得られる。

＜計算例＞

　３夜連続で測光観測して cyclocode を用いて周期解析を行った。誤差の計算に必要な各値は次のとおり。

　周期　5.2033時間
　観測の期間　2.144日
　観測数　145

　変光幅 (amplitude) の計算には、*.roc の 3 次の mag の値を採用し、最大値が 1.755、最小値が 1.532 、その差は 0.223等級。

　周波数 f は、ここでは周期を時間の単位としているので、f = 1 / ( 5.2033 / 24 ) = 4.6125 となる。

　データの誤差の分散は、式（２）によりすべての観測値について Error を求め、得られたすべての Error の分散を求めると 0.0069 となった。（これらの計算は表計算ソフトを利用すると便利）

　以上の値を式（１）に代入すると、誤差の割合として 0.0029 を得た。周期の精度は約0.3%となる。

　これに周期 5.2033 を乗じて、周期に対する誤差を求めると 0.0029 × 5.2033 = 0.015

　よって、P = 5.2033 ± 0.015 hr が得られた。

参考文献