PHPで太陽黒点相対数の周期変化を描く

(1/1)
世界は本当に温暖化へ向かっているのだろうか――。
地球のエネルギー収支では、入ってくる方のほぼ 100%が太陽エネルギーである。これが二酸化炭素などの影響で宇宙空間へ出て行かず、温暖化が進むとされている。
では、入ってくる太陽エネルギーに変動はないのか。

観測により、太陽エネルギーの放射は周期的に変動していることが分かっている。太陽の表面にある黒点の数の変化が、これに連動している。
そこで今回は、天文台の黒点観測データを利用し、PHP を使って黒点の数の周期変化をグラフに描いてみることにする。

サンプル・プログラムの実行例

PHPで太陽黒点相対数の周期変化を描く
水色の折れ線グラフは黒点相対数を表す。後述するが、黒点相対数が多いときは太陽活動が盛んである。グラフの通り、黒点相対数は約 11 年周期で増減している。
だが、相対数のピークには長周期変動があることが見て取れる。そこで、区間を 11 年(132 ヶ月)にして移動平均をとってみたのが赤い折れ線グラフである。1810 年前後と直近に、活動の低下があることが見て取れる。
1810 年前後はダルトン極小期と呼ばれている。この時期、日本では天明の飢饉が起きており、1816 年は夏のない年とされている。しかし、この年にインドネシアのタンボラ山が大爆発を起こしており、これが寒冷化の主要因とされている。
とはいうものの、ダルトン極小期では、移動平均が急に下降して、急に上昇している点に注目したい。そして直近で、移動平均が再び急下降している。新たな極小期がやってくるかもしれない。

太陽黒点と黒点相対数

第23太陽周期で最大級の黒点(2003年)
2003 年の太陽(国立天文台)
太陽黒点とは、太陽表面にあって周囲より温度が低いために黒く見える部分のことである。
左の写真が黒点で、群れになって出現する。大きなものでは、地球の数倍の大きさに及ぶ。
太陽活動が盛んなときには黒点の数が増えることが分かっている。
1ヶ月続いた無黒点(2009年)
2009 年の太陽(国立天文台)
太陽活動が低下すると、黒点の数も減る。
左の写真は、2009 年、1 ヶ月間も黒点がなかったときの太陽である。

太陽黒点の群れや個々の総量を数値化するのに、18 世紀のスイスの天文学者ルドルフ・ウォルフが考案した黒点相対数(ウォルフ数)を用いる。黒点数を f、黒点群数を g、観測地点家観測方法による補正係数を k とすると、黒点相対数は mimetex で表される。

サンプル・プログラム

解説:準備

0021: //jqPlotのあるフォルダ
0022: define('JQPLOT', '../../../../common/jqplot/');
0023: 
0024: //データファイル
0025: //define('FILE_DATA', 'http://sidc.oma.be/silso/DATA/SN_m_tot_V2.0.txt');
0026: define('FILE_DATA', './SN_m_tot_V2.0.txt');
0027: 
0028: //移動平均区間;12ヶ月×11年
0029: define('INTERVAL', (12 * 11));
0030: 
0031: //グラフの幅、高さ
0032: define('GRAPH_WIDTH',  600);
0033: define('GRAPH_HEIGHT', 500);
0034: 
0035: //グラフの色
0036: define('COLOR_SSN',  '#88CCFF');
0037: define('COLOR_MEAN', '#CC0000');
0038: 
0039: //太陽黒点相対数 推移グラフの名前
0040: define('GRAPH_SSN', 'jqPlot_SSN');
0041: 
0042: //統計に関わるクラス:include_pathが通ったディレクトリに配置
0043: require_once('pahooStat.php');

黒点相対数の推移グラフを描くのに、jQuery プラグイン「jqPlot」を利用する。
jqPlot の入手方法、インストール方法については、「PHP で NHK 政治意識月例調査をグラフ表示」を参照してほしい。

移動平均を求めるのに、自前のクラスファイル "pahooStat.php" を利用する。利用方法については、「PHP で相関係数と回帰直線を表示」を参照してほしい。

1749 年から現在までの毎月の黒点相対数の平均値を公開しているベルギー王立天文台のデータを利用することにした。Sunspot Numberからダウンロードしたテキスト・ファイルを、定数 FILE_DATA で示す名前で保存しておく。URL の方を有効にすれば、直接、ベルギー王立天文台のデータを参照する。

解説:単純移動平均

0190: /**
0191:  * 単純移動平均
0192:  * @param int   $m 区間
0193:  * @param array $x データ配列
0194:  * @param array $y 移動平均を格納する配列
0195:  * @return bool TRUE/FALSE
0196: */
0197: function simple_moving_average($m$x, &$y) {
0198:     if (! $this->is_data($x)) return FALSE;
0199: 
0200:     //区間が奇数
0201:     if ($m % 2 == 1) {
0202:         $n = floor($m / 2);
0203:         for ($i = $n$i < count($x) - $n$i++) {
0204:             $z = 0;
0205:             for ($j = $i - $n$j <= $i + $n$j++)   $z += $x[$j];
0206:             $y[$i] = $z / $m;
0207:         }
0208:     //区間が偶数
0209:     } else {
0210:         $n = floor($m / 2);
0211:         for ($i = $n$i < count($x) - $n$i++) {
0212:             $z1 = 0;
0213:             for ($j = $i - $n$j < $i + $n$j++)            $z1 += $x[$j];
0214:             $z2 = 0;
0215:             for ($j = $i - $n + 1; $j <= $i + $n$j++)    $z2 += $x[$j];
0216:             $y[$i] = (($z1 / $m) + ($z2 / $m)) / 2;
0217:         }
0218:     }
0219:     return TRUE;
0220: }

変動のあるデータをなだらかにし、全体の傾向を見るときに移動平均という統計手法を用いる。ここではデータに重み付けを行わない単純移動平均を用いることにする。

移動平均の計算方法は、移動平均を求めたい区間が奇数か偶数かによって異なる。

  • 区間が奇数( mimetex )のとき:
 mimetex 


  • 区間が偶数( mimetex )のとき:
 mimetex 

解説:データ・ファイルを読み込む

0112: /**
0113:  * データ・ファイルを読み込む
0114:  * @param string $fname  入力ファイル名
0115:  * @param array  $items  データを格納する配列
0116:  * @return int データ件数/FALSE
0117: */
0118: function readDataFile($fname, &$items) {
0119:     $infp = fopen($fname, 'r');
0120:     if (! $infpreturn FALSE;
0121: 
0122:     $cnt = 0;
0123:     $arr = array();
0124:     while (! feof($infp)) {
0125:         $str = fgets($infp);
0126:         $arr = preg_split('/[ ]+/i', $str);
0127:         if (($arr != FALSE) && (count($arr) >= 3)) {
0128:             $items[$cnt]['year'] = $arr[2];
0129:             $items[$cnt]['SSN']  = $arr[3];
0130:         }
0131:         $cnt++;
0132:     }
0133:     fclose($infp);
0134: 
0135:     return ($cnt - 1);
0136: }

先述のベルギー王立天文台のテキスト・データの 1行は、下記のような構造になっている。
内容
1int西暦年
2int
3float西暦年月(月は小数)
4float黒点相対数(月平均)
5float標準偏差
6int観測数
ここでは列 3 を X軸に、列 4 を Y軸にプロットしていくことにする。
行は固定長ではあるが、スペースで区切られていることから、組み込み関数  preg_split  で個々の要素に分離し、配列へ格納してゆくことにする。

解説:推移グラフ

0138: /**
0139:  * jqPlot用のスクリプト:年次推移グラフ
0140:  * @param array  $items データ配列
0141:  * @param string $name オブジェクト名
0142:  * @param string $title グラフ・タイトル
0143:  * @return string スクリプト
0144: */
0145: function plotProgress($items$name$title) {
0146:     //グラフの色
0147:     $color_ssn  = COLOR_SSN;
0148:     $color_mean = COLOR_MEAN;
0149:     //移動平均を求める区間(月)
0150:     $interval = INTERVAL;
0151: 
0152:     //移動平均
0153:     $x = array();
0154:     $y = array();
0155:     foreach ($items as $key=>$val)  $x[$key] = $val['SSN'];
0156:     //統計オブジェクト
0157:     $pst = new pahooStat();
0158:     $pst->simple_moving_average(INTERVAL, $x$y);        //移動平均
0159:     $pst = NULL;
0160:     foreach ($y as $key=>$val)  $items[$key]['mean'] = $val;
0161: 
0162:     //X軸の最小・最大
0163:     $xmin = floor($items[0]['year'] / 10) * 10;
0164:     $xmax = (floor(date('Y') / 10) + 1) * 10;
0165: 
0166:     //系列の生成
0167:     $s1 = '';
0168:     foreach ($items as $key=>$val) {
0169:         if (isset($val['SSN'])) {
0170:             $s1 .= sprintf('[%f, %f], ', $val['year'], $val['SSN']);
0171:         }
0172:     }
0173:     $s2 = '';
0174:     foreach ($items as $key=>$val) {
0175:         if (isset($val['mean'])) {
0176:             $s2 .= sprintf('[%f, %f], ', $val['year'], $val['mean']);
0177:         }
0178:     }
0179: 
0180:     $width  = GRAPH_WIDTH;
0181:     $height = GRAPH_HEIGHT;
0182: $js =<<< EOD
0183: <script>
0184: jQuery(function() {
0185:     jQuery.jqplot('{$name}',
0186:     [
0187:         [ {$s1} ],
0188:         [ {$s2} ]
0189:     ],
0190:     {
0191:         //ラベル
0192:         series: [
0193:         {
0194:             label: '月平均',
0195:             color: '{$color_ssn}'
0196:         },
0197:         {
0198:             label: '{$interval}ヵ月移動平均',
0199:             color: '{$color_mean}'
0200:         },
0201:         ],
0202:         legend: {
0203:             show: true,
0204:             placement: 'inside',
0205:             location: 'ne',
0206:             renderer: $.jqplot.EnhancedLegendRenderer,
0207:             rendererOptions: { numberRows: 1 }
0208:         },
0209:         //グラフ
0210:         seriesDefaults: {
0211:             showLine: true,
0212:             rendererOptions: { smooth: false },
0213:             markerOptions: { size: 0 },
0214:         },
0215:         //軸
0216:         axes: {
0217:             xaxis: {
0218:                 label: '西暦年',
0219:                 tickOptions: { formatString: '%d' },
0220:                 min: {$xmin},
0221:                 max: {$xmax},
0222:             },
0223:             yaxis: {
0224:                 label: '黒点相対数',
0225:                 tickOptions: { formatString: '%.1f' },
0226:                 min: 0.0
0227:             }
0228:         }
0229:     }
0230:     );
0231: });
0232: </script>
0233: <div id="{$name}" style="width:{$width}px; height:{$height}px;"></div>
0234: 
0235: EOD;
0236: 
0237:     return $js;
0238: }

配列に格納したデータから移動平均を求め、黒点相対数と移動平均の推移を jqPlot スクリプトへ展開し、折れ線グラフとする処理は、ユーザー関数 plotProgress で行う。

マウンダー極小期

過去400年間の黒点数
過去 400 年間の黒点数
観測記録が少ないものの、1645 年から 1715 年にかけ、黒点数が著しく減少した期間があり、マウンダー極小期と呼ばれている。
イギリスではテムズ川が凍結し、アメリカではニューヨーク湾が凍結した。また、アイスランドは海氷に取り囲まれ、食糧不足で人口が半減した。日本でも、この時代に飢饉が頻発し、百姓一揆が起きるようになった。

太陽活動と小氷期の関係については、専門家による研究を待たねばならないが、地球のエネルギー収支があるのは事実だから、温暖化ガスの影響だけで温暖化や寒冷化を論じるべきではないだろう。

参考サイト

(この項おわり)
header