2014年10月31日 (金)

三角関数の計算回数

■三角関数の計算回数
 Sunrise
 
1.計算式どおりに計算する場合
(赤緯δ,均時差e,北緯φ,東経θ, 時差td,補正定数Ds)
 時刻計算              //三角関数計算回数:6回
  a =  -sin(Ds)-sin(δ)*sin(φ);     //3
  b =  cos(δ)*cos(φ);              //2
  時角 t = 180/M_PI*acos(a/b);       //1
 角度計算               //三角関数計算回数:15回
  高度 h = asin(sin(φ)*sin(δ)+cos(φ)*cos(δ)*cos(t));//6
  方位 A
    sinA = cos(δ)*sin(t)/cos(h);         //3
    cosA = (sin(h)*sin(φ)-sin(d))/cos(h)/cos(φ); //5
    A = atan2(sinA, cosA) + M_PI;          //1
(合計計算回数): 36回
  1地点につき,時刻1回,方位2回(6+15*2=36)
  
2.1度計算した値を保持する場合
 時刻計算              //三角関数計算回数:5回
  (同一地点では固定値のもの:日出日入時刻の値)
  δのsin,cos  ==>  変数sind, cosd に保持  //2回
  φのsin,cos  ==>  変数sinp, cosp に保持  //2回
  時角 t = 180/M_PI*acos(a/b);             //1回
 角度計算              //三角関数計算回数:6回
  (任意時刻の値)
  時角 t のsin,cos  ==>  変数sind, cosd に保持  //2回
  高度 h = asin(sinp*sind+cosp*cosd*cost);      //1回
  高度 h のsin,cos  ==>  変数sinh, cosh に保持  //2回
  方位 A = atan2(sinA, cosA) + M_PI;          //1回
     (sinA = cosd*sint/cosh;      )
     (cosA = (sinh*sinp-sind)/cosh/cosp;)
(合計計算回数):17回
  1地点につき,時刻1回,方位2回(5+6*2=17)
  (南中高度は太陽赤緯と緯度の加減算で計算できる)
  
・一覧表画面の1ページ分の三角関数計算回数
 360回 ==> 170回
 
・グラフ計算の場合の三角関数計算回数
 グラフ計算は,25点の時刻における方位と高度を計算する.
  1の場合:6+15*25 = 381回/1グラフ
  2の場合:5+6*25  = 155回/1グラフ
 381回 ==> 155回
 
いずれの場合も,実行時間が,半分程度に短縮される.
 
[HP200LX/LXDIC]

|

2014年10月20日 (月)

太陽方位角の変化

■太陽方位角の変化
太陽の方位角は,時間の変化に対して一定の割合で変化するものではない.
北半球の夏場は,昼の12時前後に大きく変化し,日出日入前後では変化
が小さくなる.そして,真夜中の変化率も大きくなるが,昼には及ばない.
冬場には,真夜中の変化が大きくなり,昼の変化を上まわる.
南半球では,夏場と冬場が逆になる.
 
したがって,0時の方位角と12時の方位角の差が,180度になるとは限ら
ない.(南中時刻が12時の場合は180度になる)
 
例)東京 12/20
Sr6902
1時間ごとの変化
Sr6901
この例では,12時の方位は,真南(180度)から少し進んでいる.ここから
12時間たつと,夜中の方位角の変化が大きいため,0時の方位角は,昼間よ
少し余分に進む.結果として,0時と12時の間に方位角は180度以上進む.
 
  
東京 6/20
Sr6911
 
明石(南中時刻がほぼ12時) 方位角の差は約180度
Sr6912
 
沖縄  方位差180度からのずれが大きい
Sr6913
(太陽がほぼ真上を通過するため,方位は南中時刻の前後でほぼ180度変化する)
 
これらの結果は,実生活上の経験と乖離しているように感じるかもしれないが,
これは,方位角が太陽の方向を地平面に投影した角度であり,太陽の位置情報
(高度と方位)の一側面でしかないことによる.
実際,太陽が真上近くにある時は,方位角にはあまり意味がない.
 
[HP200LX/LXDIC]

|

2014年10月17日 (金)

太陽高度方位計算プログラム

太陽高度方位計算プログラム
Sunrise Ver 7.2--7.9
 
■機能
世界各地の日出/日入/南中時刻,太陽高度および方位を計算する.
任意の位置,日付,時刻についての計算が可能で,太陽の高度と方位に
ついては,グラフ表示も行う.
 
(概要)
・日出日入時刻,太陽高度/方位の計算
・任意の地点に対応
  東経: -180°~ 180°
  北緯:  -90°~ 90°
  時差:  -12 ~ +12  (日本は+9:JST=UTC+9)
・浮動小数点数を用いた詳細な計算
  誤差:日出日入時刻:1分以内
     方位,高度 :1度以内
・表示項目
  場所(緯度,経度)
  日付,時刻
  日出時刻,日出方位
  南中時刻,南中高度
  日入時刻,日入方位
  太陽高度,太陽方位(指定した時刻の値)
  日出日入一覧表画面の表示
  グラフ表示(高度の時間変化,高度・方位関係)
・任意の位置,日付,時刻を指定しての計算が可能
  
■ファイル一覧
・SRISE.EXM       実行ファイル
・SRISE.HLP       ヘルプファイル
・SRISE.INI       設定ファイル

■セットアップ
・SRISE.EXM を任意のディレクトリに置いて,AppMgr または,MoreEXM
 に登録する.
・SRISE.INI を編集して,好みの地点(都市)の位置情報を登録する.
 10都市登録できるが,起動時には先頭の都市が表示される.
・SRISE.HLP, SRISE.INI を C:\_DAT にコピーする.

■使い方

●メイン画面
・起動すると,現在の日付・時刻における値が表示される.
・グラフ表示は,F6(グラフ)キー.
  方位図==>グラフ1==>グラフ2==>方位図
・表示情報は,時間が経過しても自動更新されない.
・現在時刻の情報に更新する場合は,Spaceキー.
 任意の時刻を指定することも可能 ==> 設定画面
 
 都市選択
・F3(都市)キーで,登録済み都市を選択できる.
 選択キーは,Fキー(F1~F10)または,数字キー(0~9).
・数字キー(0~9)による選択は,常時可能になっている.
・都市名一覧表示(表示のみ)は,'c'キー.
 消去するときは,F2(再描画)キーで再描画する.
 
 ヘルプ表示(F1キー)
・F1でヘルプを表示する.(ウィンドウは作成されない)
・F1キーを押すたびに,表示ページが切替わる.(2ページ)
・消去は,F2(再描画).
 
 小数表示
・時刻・角度の小数表示が煩わしい場合は,TABキーで変更できる.
 TABキー:整数表示 <==> 小数表示(小数第1位まで表示)
 
●設定画面
・世界の任意の位置,日付,時刻を設定して計算できる.
・それぞれ,F8(位置),F10(日付),F9(時刻) キーで個別に設定する.
・不正な数値の入力は受けつけない.
・設定中止は,ESCキー.
 
 
●データ一覧画面
・メイン画面で,F7(一覧表)キーにより表示する.
・上下の矢印キー(↓↑)で,1行ずつスクロール
 左右の矢印キー(→←)で,10行ずつスクロール
・表示日付の間隔を変更できる.F9(日間隔)(10=>1=>7=>10)
・Shift+↓↑キーで,表示日付を1日ずつ調整できる.
・都市の変更は,数字キー(0~9)で常時可能.
・F10(Main)キー,または,Enter キーで,メイン画面に戻る.
 
[HP200LX/LXDIC]

|

2014年10月11日 (土)

整数演算化

■整数演算化
前回の結果により,基準年のデータ表参照により,すべての年の日出日入
データを導出できることが分かった.(整数演算のみ使用)
この表から整数演算のみで,隣接する他の地点のデータを導出でるかを検
討する.
 
(整数演算化)
 想定する範囲:東経 -180°~ 180°
        北緯 20°~ 54°
 対象となる項目:日出時刻,日入時刻,南中時刻,
         日出方位,日入方位,南中高度
         
・日出時刻,日入時刻
 経度の変化(M)は,経度15度=1時間で補正する.
 同一の経度のもとで,緯度を変化させると,日出日入時刻が変化する.
 理科年表の補正法(北緯20~50)では,この変化量を2つの因子に分解し
 ている.
   日出 T = T0 + M - n*N (N,n は表から求める)
   日入 T = T0 + M + n*N 
 ここで,n は昼長の1/2(東京)のみに依存し,N は緯度にのみ依存する.
 ただし,ある程度以上の精度を求めると,nもNと同様に,緯度の両方に依
 存してしまい,2つの因子に分解したことが無意味になる.
 
 コンピュータプログラムとしては,最初から単純に1つの補正値 nN を計
 算する方が効率的で精度が高い.(nNは緯度と昼長に依存する)
 北緯(20~54°)は1度刻み,昼長/2(290~440分)5分刻みに,nN-値を
 計算し参照テーブルを作成する.(35×31×2=2170バイト)
 予備的な計算では,誤差は 0.2 分程度に抑えられた.
 
・南中時刻
 経度の補正のみ.経度15度=1時間.
 
・南中高度
 東京における南中高度を,目的地点の緯度と東京の緯度の差により補正
 する.
 
・日出方位,日入方位,
 昼長が長くなると,その分,日入の方位角は増加し,日出の方位角は減
 少する.(この変化量は,緯度に依存する.)
 北緯(20~54°)は1度刻み,昼長/2(290~440分)5分刻みに,角度
 補正値を計算し,参照テーブルを作成する.(35×31×2=2170バイト)
 予備的な計算では,誤差は 0.2°程度に抑えられた.
 (PCでの計算:浮動小数点を用いた計算数値との差)
 
※整数演算化できた項目
 日出時刻,日入時刻,南中時刻,昼の長さ,
 日出方位,日入方位,南中高度
※整数演算化できていない項目
 任意時刻の太陽方位,太陽高度
 
[HP200LX/LXDIC]

|

2014年9月30日 (火)

Sunrise::データ表参照による計算

■データ表参照による計算
日出日入時刻等のデータの日々の変化は,毎年同じパターンの繰返しになる.
場所を固定すれば,基準年1年分の計算データから,補正により全ての年の
計算結果を導出できる.
事前計算データ
 基準年:2016年(通算日数の補正値が0になる)
 ファイルサイズ:2940バイト(4+367*8)
 データ内容:
   位置情報 2+2バイト
   日出時刻,日入時刻,日出方位,日入方位,南中高度 8バイト
 ファイル名:c:\_dat\srise$?$.zzd  (?に都市番号が入る)
 
処理内容
・良く利用する地点のデータを1年分保存する.
 データ保存機能を追加(10地点まで保存可能)
・利用の際は,都市番号と位置情報を照合し,一致した場合だけ使用する.
 不一致の場合,事前計算データがない場合は,通常の計算を行う.
・年ごとの補正は,整数計算で行い,浮動小数点は使わない.(高速)
・この計算は,一覧表画面の計算で用いる.
・メイン画面では,浮動小数点による近似式計算を用いる.
 
計算の比較
1.近似式による計算(==>部分がHP200LXでリアルタイムに行う計算)
  赤緯/均時差データ --> フーリエ係数計算 --> 赤緯/均時差計算 ==> 日出日入計算
2.データ表参照による計算
  赤緯/均時差データ --> 日出日入計算 ==> 年次変化補正計算
元になるデータが同じなので,両者の計算値には,ほとんど差がない.
差は全ての項目で,最大 0.1~0.2(端数処理の誤差程度).
 
[HP200LX/LXDIC]

|

2014年9月24日 (水)

Sunrise::近似式再考

■近似式再考
前回の結論は,通算日数に補正を加えれば,ある基準年の係数を用いて,
すべての年の赤緯/均時差の近似計算ができるということであった.
しかし,良く考えてみれば,わざわざ近似式を計算する必要がないことに
気付く.すなわち,赤緯/均時差の実データを直接利用し,補正を加えな
がら,他の年にも使い回わせば良い.
これなら,1件につき合計で 16回になる三角関数の計算が不要になり,
CPUの負荷が軽減される.特に,HP200LX のような非力なCPUを用いる場合
には好ましい.
 
ある年の赤緯/均時差の 366日分のデータをとり,これを元に補正を加え
ながら,他の全ての年に適用する.(次年の1/1を含むデータ)
(閏年を基準年とする場合の通算日数の補正)
 y%4==0 (閏年):  0    :J = days;
 y%4==1 の場合: +0.75 :J = days + 0.75;
 y%4==2 の場合: +0.5  :J = days + 0.5;
 y%4==3 の場合: +0.25 :J = days + 0.25;
J が整数ではなくなるので,求める値は,整数点のデータから,内挿によ
り計算する.
 
●データ量:2202バイト(6×367)(比較的小さなデータになる)
 赤緯 :long(4バイト)秒数で表し整数化
 均時差:int(2バイト) 秒数で表し整数化
(理科年表あるいは高精度計算サイトのデータを使うことが可能)
 
●計算結果
・4年分の赤緯/均時差データ表を用いた計算
・1年分の赤緯/均時差データ表(2016)を補正しながら使う計算
上記2つの計算の誤差:4秒程度(近似式を使うより誤差が小さい)
 2014: 日出 (-00.05, 00.02)分   日入 (-00.06, 00.03)分
 2015: 日出 (-00.03, 00.01)分   日入 (-00.04, 00.03)分
 2016: 日出 ( 00.00, 00.00)分   日入 ( 00.00, 00.00)分 <--基準年
 2017: 日出 (-00.03, 00.03)分   日入 (-00.02, 00.03)分
数年のスパンで考えれば,365.25日ごとに全く同一の変動パターンを繰り返す
と見なしても,誤差はほとんど生じない.
 
●三角関数の計算回数
 一覧表1日分: 67回 ==>  51回 ( -16回)
 グラフ1日分:1125回 ==> 1125回 (変らず)
 一覧表の都市切替時間 1.5秒 ==> 1.2秒 (7日分の計算時間)
 
●計算高速化の可能性
数年のスパンで考えれば,365.25日ごとに全く同一の変動パターンを繰り返す
のであるから,良く利用する都市の日出日入時刻等のデータについては,予め
1年分を計算しておけば,すべての年に対応できる.
浮動小数点計算の必要がなく,大幅な高速化が可能になる.
  
[HP200LX/LXDIC]

|

2014年9月22日 (月)

Sunrise 開発メモ

■Sunrise 開発メモ

●近似式の改良
赤緯/均時差の計算は近似式によっているが,これが全ての計算結果を左右
する.これまでの近似式はあまり精度が高くなく,日出日入時刻の計算誤差
が1~2分程度になっていた.
 
今回は,実データを与えて,より精度の高い係数を求めた.
この近似式は,赤緯/均時差を表す関数をフーリエ級数展開したもので,係
数は以下の計算で求められる.
 an = 1/π∫f(x)cos(nx)dx   (x=0 --> 2π)
 bn = 1/π∫f(x)sin(nx)dx   (x=0 --> 2π)
 f(x) = a0/2+Σ(an*cos(nx)+bn*sin(nx))  (n=1,2,...)
 
計算結果から,4次の係数までとれば,十分な精度が得られることが判明.
赤緯(計算結果)[単位:度]
 a0/2=  0.394472
  a1=-22.899438
  a2= -0.388787
  a3= -0.153999
  a4= -0.009749
  a5= -0.003751
  b1= +4.027593
  b2= +0.045014
  b3= +0.084584
  b4= +0.005171
  b5= +0.003806
均時差 (計算結果)[単位:時間]
 a0/2=  0.000334
  a1= +0.008998
  a2= -0.056444
  a3= -0.001506
  a4= -0.002352
  a5= -0.000155
  b1= -0.122402
  b2= -0.155564
  b3= -0.005427
  b4= -0.003044
  b5= -0.000344
 
●年変化の補正
赤緯/均時差の年変化は,実用上の精度に関しては,無いに等しい.
すなわち,365.25日を周期として同じ変動パターンを繰り返すと見るこ
とができる.
ただし,年ごとの変動パターンと暦の間には,毎年約 0.25日ずつのズレ
が生じるため,これを補正する必要がある.
 
1/1からの通算日数の計算に補正を加える.(yは西暦)
 y%4==0 (閏年): -0.75
 y%4==1 の場合:  0
 y%4==2 の場合: -0.25
 y%4==3 の場合: -0.50
毎年,変動パターンの同一位置より,約0.25日分前にずれて,1/1を基点
とする通算日数は大きくカウントされるため,その分を差引いて補正する.
閏年の翌年は,前年の閏日の挿入により,暦と変動パターンのずれがリ
セットされる.(正確には,0.03日分残る)
 
(例)閏年の1/1を基点にする場合の通算日数補正
 y%4==0 (閏年):  0    
 y%4==1 の場合: +0.75 
 y%4==2 の場合: +0.5  
 y%4==3 の場合: +0.25 
フーリエ係数は,基点となる年のデータで計算する.
 
(補足)
上記計算の誤差(1年365.2425日で計算)
 基準年:   0
 4年後:  +0.03(1.0-0.2425*4)
 8年後:  +0.06
 12年後:  +0.09
 16年後:  +0.12
 20年後:  +0.15

[HP200LX/LXDIC]

|

2014年9月17日 (水)

Sunrise 開発メモ

■Sunrise 開発メモ
追加メモ
・計算精度向上(時刻1分以内,角度1度以内)
・閏年対応(精度向上に伴い年変化が無視できない) 
・日付設定に"西暦"を追加
・”都市設定モード”を追加
 ノーマルフォントでファンクションガイドに都市名表示
 
年設定
Sr50lp1
閏年 
Sr50lp2
 
太陽方位の時間変化
Sr50b24h  
  
ノーマルフォントのファンクションガイド
Sr51bmn
  
[HP200LX/LXDIC]

|

2014年9月13日 (土)

Sunrise 開発メモ

■Sunrise 開発メモ
Ver 4.6c
 
現在の仕様 
 
●ウィンドウ構成
・メイン画面(1)
  数値表+グラフ1,又は,グラフ2
・設定ダイアログ(3)
  位置設定
  日付設定
  時刻設定
・データ一覧画面(1)
  7週分の日出日入等の数値表
・メニュー
  &Location, &Date, &Time, &Help, &Quit
  
  
●詳細内容
・グラフの種類
 グラフ1:太陽高度の時間変化+太陽方位の簡易表示
 グラフ2:太陽高度と方位の関係図
 (グラフ1)太陽高度の時間変化
  Sn46c02
 (グラフ2)太陽高度と方位の関係
  Sn46c01
 
・メイン画面のファンクションキー割当て
 F1  Help         'h'
 F2  再描画       '/'
 F3  Cities       '\'
 F4  夏時間       Del
 F5  現時刻       space
 F6  グラフ       right, left
 F7  一覧表       up,down
 F8   位置        'l'
 F9   時刻        't'
 F10  日付        'd'
 
 
・データ一覧画面のファンクションキー割当て
 F2  再描画        '/'
 F3  Cities        '\'
 F10 Cancel        ESC
(データ一覧画面)                           ※方位表示の不具合修正
 Sr46cerr
  
 
・登録できる都市数:10
 都市選択は,メイン画面と一覧表画面で数字キーにより行う.
 
・ヘルプ表示(F1キー)
・都市名一覧表示(F3キー)
 メイン画面への単なる上書表示で,簡易的なもの.
 消去するときは,F2で再描画
 (ヘルプ画面)表示中もすべての操作が有効
 Sn46c05
 (都市名表示)表示中もすべての操作が有効
 Sn46c03 
 
・データキャッシュによる高速化
 浮動小数点演算の負荷が大きく,計算に時間がかかるため,
 表示の遅延を回避する目的で以下のデータに実施する.
  日出日入一覧表用計算データ(50日分)
  グラフ用計算データ(10日分)
 
・設定ファイル
#(設定順序)
#   東経,北緯,時差,地名
#(設定項目)
#   東経:実数または整数  -180°<= lo <= 180°(20バイト以内)
#   北緯:実数または整数   -90°<= la <= 90°(20バイト以内)
#   時差:整数で指定  -12 <= h <= +12  (日本JST=UTC+9)
#         夏時間を考慮しない場合の時差を設定する.
#   地名:8 バイト以内.それ以上は無視される.
#(規則)
#   ・各項目は半角コンマ [,] で区切る.
#   ・各項目の間には,半角スペースを挿入して整形できる.
#   ・最大で10地点を設定できる.
#   ・11個目以降の設定は無視される.
#   ・#で始まる行は,コメント行
#   
139.7414, 35.6581,  +9, 東京
136.9167, 35.1667,  +9, 名古屋
135.0000, 35.0000,   9, 明石
127.6667, 26.2167,  +9, 那覇
141.3500, 43.0667,  +9, 札幌
-0.12806, 51.50778,  0, London
-73.54,   40.46,    -5, NewYork
 
・エラー表示
 設定ファイルエラーの表示(位置情報が全く取得できない場合)
 Sn46c06
 
[HP200LX/LXDIC]

|

2014年9月 7日 (日)

Sunrise Ver 3.8c

■Sunrise Ver 3.8c
開発メモ
 
追加
・緯度経度の制限を外す(全地域をカバー)
 東経: -180°~ +180°, 北緯: -90°~ +90°
・日出・日入時刻が存在しない場合が発生する.(高緯度の場合)
・日出日入で計算エラーが出ても,計算を中断せず高度と方位を計算
 1日中昼の場合は,日出00:00,日入24:00,昼長24:00 と表示
 1日中夜の場合は,日出--:--,日入--:--,昼長 0:00 と表示
・グラフ枠を3ドット下へ拡大
 高度0度の線が見えるようにする.
・グラフの現在時刻(又は指定時刻)の位置に縦線を入れる.
 
 
浮動小数点使用の影響(Borland C++ 3.1)
・三角関数を自由に使えるので,自在に計算できる.
・数学ライブラリをリンクするため,プログラムサイズが大きい.
・関数計算処理の負担が大きく,グラフ表示と一覧表示が遅くなる.
 
 
今後の課題
・場所選択方式の変更
 各都市を常時ファンクションキーに割当てると他の多くの機能の割当てができない.
 必要な時に選択可能な場所を表示し,選択する方式にする必要がある.
・計算速度の改善
 浮動小数点ではなく,整数による計算が可能かどうか検討する.
 計算が高速化すれば,純粋に表示に起因する部分の速度向上も必要になる.
 
 
実行例
高緯度の場合(北緯85度)
Sn38c1  
 
北極(北緯90度)
Sn38c6  
 
南半球(北緯 -35度)
Sn38c5  
 
[HP200LX/LXDIC]

|

«Sunrise Ver 3.7