2019年春学期のパソコン教室での実習でつかった教材をすこし編集しなおしたものです。

「大気上端」に入射する太陽放射

主題: 地球の公転・自転にともなって、 「大気上端」に入射する太陽放射 (単位時間・単位面積あたりのエネルギーの流れ)は どのように変化するか。

手段: 「R」をプログラム言語処理系として使い、数値計算と作図のプログラムを動かす。

(1) 日変化

緯度・経度をあたえられた地点での、 1日のうちの各瞬間の太陽の方位・高度角を計算し、 「太陽定数」を定数として、「大気上端」に入射する日射量を計算する。 現在の地球の軌道要素にもとづく 「均時差」(『図表』 49ページ [4 A] 参照) の計算がはいっている (ただし単位が、『図表』では時間の分、プログラムでは角度のラジアンである)。

計算式は、 立正大学 地球環境科学部の 中川 清隆 先生 [外][中川 清隆 研究室]の [外][太陽方位、高度、大気外日射量の計算]の ページの記述にもとづいた。

増田が、プログラムを R 言語で書きかえ、グラフの作図を追加した。 (ファイル名が nakagawa ではじまっているのは、中川先生の計算式にもとづいたという意味であり、 中川先生のプログラムという意味ではありません。)

プログラム例

• nakagawa_test_kinzisa.R ... 「均時差」を計算してファイル "eqtime.csv"に書き出す。例文作成作業の途中の動作確認用。
  • nakagawa_test_kinzisa_utf8.R 文字コード UTF-8
  • nakagawa_test_kinzisa_sjis.R 文字コード Shift-JIS
• nakagawa_work1.R ... 1つの時刻の数値を計算する。
  • nakagawa_work1_utf8.R 文字コード UTF-8
  • nakagawa_work1_sjis.R 文字コード Shift-JIS
• 1時間きざみで24回の数値を計算し、CSVで書きだす。
  • nakagawa_work2_utf8.R 文字コード UTF-8
  • nakagawa_work2_sjis.R 文字コード Shift-JIS
  • 出力CSV例 hour_insol_day140.csv
• 1時間きざみのCSVファイルを読み、グラフにする。そのままだと画面表示。「#G」をはずすと画像ファイル出力。
  • plothourly_utf8.R 文字コード UTF-8
  • plothourly_sjis.R 文字コード Shift-JIS
  • 出力画像例 hour_insol_day140.png
• 1時間きざみのCSVファイルを 2つ読み、同じ画面上にちがう記号でグラフにする。
  • plothourly2_utf8.R 文字コード UTF-8
  • plothourly2_sjis.R 文字コード Shift-JIS
  • 出力画像例 hour_insol_day140_173.png
• 1時間きざみのCSVファイルを 2つ読み、両者の差の値をグラフにする。
  • plothourlydiff_utf8.R 文字コード UTF-8
  • plothourlydiff_sjis.R 文字コード Shift-JIS
  • 出力画像例 hour_insol_day173_minus_day140.png

注

• 日付は1月1日を1とした年ごとのとおしの日付を使っている。 「月・日」とのあいだの相互変換は、Rのライブラリを使う方法も、 自まえでプログラムする方法もあるが、例文にはまだくみこんでいない。
• 日本語文字をあつかうところで、Windows 7 上の R では、 文字コードを 「Shift_JIS」にする必要があるようだ。 ほかのOSでは「UTF-8」にする必要があるかもしれない。 プログラム例はなるべく両方用意するが、文字化けしたらコード変化をこころみてほしい。
• 例文のプログラムを画面に表示させて保存すると、 ブラウザによっては、ファイル名の「拡張子」部分がが「.txt」になったり、 文字コードが変わったりする。 プログラムを表示させず、ファイルへのリンクが書かれているページで、 リンクを(たいていのブラウザでは)右クリックして、 「(Firefoxのばあいでいえば）名まえを付けてリンク先を保存」の機能で保存するのがよい。

演習課題

• nakagawa_work2_sjis.R をもとに、 次の変数に代入される数値をひとつずつ、(プラスマイナス10日以上・10度以上の差で)かえて、結果を比較し、結果の特徴を簡単にのべてください。
  • daytoyear 元旦からの通し日数 [日] (例文では5月20日にあたる数値を 31+28+31+30+20 としてあたえている)
  • geolatdeg 緯度 [度、北緯が正] (例文の値は板倉の緯度)
  • geolondeg 経度 [度、東経が正] (例文の値は板倉の経度)
  • 結果をファイルに書き出すばあいは、値を変えたものを区別できるようにファイル名をつけること。 プログラム nakagawa_work2_sjis.R のテキストファイルに書き出すところでは、 dayofyearの値がファイル名にくみこまれるようにしたが、緯度・経度はくみこんでいない (自動化してうまくいくか不安なので)。 必要に応じて、nakagawa_work_sjis.R のおわり付近の「fname <- ...」の文の 文字列 "hour_insol_day" をエディタ作業で変更してほしい。
  • solconst いわゆる太陽定数 [W/m2] をかえてみてもよいが、定数倍されるだけなので、レポート課題にはふくめない。

結果の例

緯度・経度・日付を指定し、1日ぶんの1時間きざみの値を、CSVファイルに書きだす。

• 緯度・経度は板倉のもの、日付は5月20日。
  • nakagawa_work2_sjis.R
  • 出力CSV hour_insol_day140.csv
• 日付を変える
  • 板倉の緯度経度、6月22日 (ほぼ夏至)
    • calc_hourly_day173_sjis.R
    • 出力CSV hour_insol_day173.csv
  • 板倉の緯度経度、12月22日 (ほぼ冬至)
    • calc_hourly_day356_sjis.R
    • 出力CSV hour_insol_day356.csv
• 経度を変える
  • 緯度は板倉と同じ、東経120度、5月20日
    • calc_hourly_lon120_sjis.R
    • 出力CSV hour_insol_lon120_day140.csv
  • 緯度は板倉と同じ、東経150度、5月20日
    • calc_hourly_lon150_sjis.R
    • 出力CSV hour_insol_lon150_day140.csv
• 緯度を変える
  • 北緯 23.4度(北回帰線)、東経 135度、6月22日
    • calc_hourly_lat234_day173_sjis.R
    • 出力CSV hour_insol_lat234_day173.csv
  • 北緯 66.6度(北極圏の緯線)、東経 135度、6月22日
    • calc_hourly_lat666_day173_sjis.R
    • 出力CSV hour_insol_lat666_day173.csv

1時間きざみのCSVファイルを2つ読み、グラフにする。PNG画像ファイルに書きだす。

• 緯度・経度は板倉のもの、5月20日(140)と6月22日(173)。
  • plothourly2_sjis.R
  • 出力画像 hour_insol_day140_173.png
• 日付を(大きく)変える。6月22日(173, ほぼ夏至)と12月22日(356, ほぼ冬至)。緯度・経度は板倉のもの。
  • plothourly2_day173_356_sjis.R
  • 出力画像 hour_insol_day173_356.png
  • 結果の特徴: 振幅が夏至に大きく冬至に小さい。
• 経度を変える。東経120度と東経150度。緯度は板倉のもの、日付は5月20日(140)。
  • plothourly2_lon120_150_sjis.R
  • 出力画像 hour_insol_day120_150.png
  • 結果の特徴: 位相がちがう。東経120度では日本標準時(東経135度の時刻)で13時、東経150度では11時に極大となる。
• 緯度を変える。北緯23.4度と北緯66.6度。緯度は東経135度、日付は6月22日(173, ほぼ夏至)。
  • plothourly2_lat234_666_sjis.R
  • 出力画像 hour_insol_lat234_666.png
  • 結果の特徴: 極大値は低緯度のほうが大きいが、日射のある時間は高緯度のほうが長い。