2014年12月18日 (木)

Sunrise::標高補正の誤差評価

■Sunrise::標高補正の誤差評価
 
国立天文台のこよみの計算
http://eco.mtk.nao.ac.jp/cgi-bin/koyomi/koyomix.cgi
と比較では,概ね良好な結果になった.
 
※以下の設定数値は必ずしも正確な値ではないが,同一の数値で計算して
 いるので,比較上の問題はない.
 
  
●松本: 東経137.97,北緯36.25,標高610m
Matu2b
 
Srmatu2
 
●富士山頂:東経138.73,北緯35.36,標高3772m
Huji2c
  
Srhuji2  
 
富士山頂:補正がない場合(0m)
Srhuji0
 
[HP200LX/LXDIC]

|

Sunrise::公開版

■Sunrise Ver 10.2 公開版
 
太陽高度方位計算プログラム
Sunrise Ver 10.2R ==>
http://hp200lx.cocolog-nifty.com/blog/download.html 
 
注意
・このプログラムは,LXDICの認証機能を利用する.
・起動時に,A:\_DAT\LXDIC.KEY をチェックし認証を行う.
・著作権は kiz に帰属する.
・問題が発生した場合は,予告なく配布を中止し,プログラムを
 回収する.

  
変更点
・標高による時刻補正および方位角補正を導入した.
・位置情報の設定項目に,標高を追加した.
・設定ファイルの仕様を変更した.
 以前のバージョンとの互換性はない.
  
[HP200LX/LXDIC]

|

2014年12月16日 (火)

日出入時刻の標高による補正

■日出入時刻の標高による補正
 
標高が高くなると,日出時刻は早くなり,日入時刻は遅くなる.
その程度は以外に大きく,標高 0m との差は,標高 418m で4分程
度,標高 1300m で7分程度になる.
標高 6m の東京でも,0.4分程度の差が出てしまうため,標高による
補正を回避することはできない.
 
●補正方法
 補正計算は,国立天文台の資料に基づいて行う.
 
 補正式:標高 h[m]
  標高補正 = 2.09*sqrt(h) [分] = 0.034833*sqrt(h) [度]
  日出日入時の実際の太陽高度 = 0-(太陽視半径+大気差+標高補正)  ※
    太陽視半径:16′0″ (平均)
    大気差  :35′8″
  ※この高度になる時刻を計算すれば,日出日入時刻が得られる.
         
 補正範囲:標高 0~9999m
 
 
(参考資料)
国立天文台報第5 巻 91 95 (2001)
『日出入時刻計算における標高の効果について』
http://www.nao.ac.jp/contents/about-naoj/reports/report-naoj/p91.pdf
 
 
 
東京(標高6m)の補正値:0.4~0.5分,0.1度
Sralti1
標高補正なし(0m)
Sralti2_2 
 
[HP200LX/LXDIC]

|

2014年12月12日 (金)

LXDIC Ver 16.0f

■LXDIC Ver 16.0f
 
変更点
・LXDIC のライトスリープ制御を中止した.
・起動キーを更新した.(12か月分添付)
 
[HP200LX/LXDIC]

|

2014年12月11日 (木)

Sunrise::テスト版

■Sunrise::テスト版
 
太陽高度方位計算プログラム
Sunrise Ver 9.9 テスト版
----ダウンロード終了----- 
 
テスト版の制約
・想定外のバグや不具合が残っている可能性がある.
・高緯度地域については,誤差の評価が困難なため,対応地域から除外
 しているが,テスト版では計算可能になっている.
・テスト版では,設定値あるいは入力値の一部について,制限を外して
 いるため,ユーザー自身が数値をチェックする必要がある.
・テスト版では,標高を0メートルに固定している. 
   
■機能
各地の日出/日入/南中時刻,太陽高度および方位を計算する.
 
■仕様概要
・日出日入時刻,太陽高度/方位を固定小数で高速に計算
・計算は標高0メートル基準,現行暦準拠
・対応地域
  東経: -180°~ 180°
  北緯:  -65°~ 65°
・計算誤差
  日出日入時刻:1分以下 
  方位,高度 :1度以下 
・表示項目
  位置(緯度,経度)
  日付,時刻
  日出時刻,日出方位
  南中時刻,南中高度
  日入時刻,日入方位
  太陽高度,太陽方位
・画面構成
  メイン画面:基本データ+太陽方位図/グラフ
  一覧表画面:データ一覧(10日分)
・設定ファイルに位置情報を登録して使用(10地点)
・位置/日付/時刻の一時的な設定も可能
 
■ファイル一覧
・SRISE.EXM       実行ファイル
・SRISE.HLP       ヘルプファイル
・SRISE.INI       設定ファイル

■セットアップ
・SRISE.EXM を任意のディレクトリに置いて,AppMgr または,MoreEXM
 に登録する.
・SRISE.INI を編集して,好みの地点(都市)の位置情報を登録する.
・SRISE.HLP, SRISE.INI を C:\_DAT にコピーする.
 
[HP200LX/LXDIC]

|

2014年11月29日 (土)

LX用64ビット電卓

■LX用64ビット電卓 
 
整備した64 ビット整数を用いて,LX用の64ビット電卓プログラムを作成
 
●64ビット電卓プログラム
・基本機能のみのシンプルな試作プログラム(動作確認用)
・数値は,10進,16進,2進で入力
(主な演算)
 右シフト,左シフト,NOT, AND, OR, XOR
 加算,減算,乗算,除算,剰余算
 
表示例:除算(商+剰余)
64calc02  
 
乗算
64calc03  
 
XOR
64calc04
 
ゼロ除算エラー 1/0 
Zerodiv  
[HP200LX/LXDIC]

|

2014年11月27日 (木)

LX用64ビット整数の整備

■LX用64ビット整数の整備
通常では32ビットまでしか扱えないHP200LXのプログラミング環境で,64
ビット整数を扱えるようにするため,簡易的な64ビット整数を定義し,自作
の演算関数を整備した.
 
 int64:64ビットデータを格納する構造体
 
主な関数(仮の名称)
 int64 lxshr(int64 a, int n);    //右シフト a>>n
 int64 lxshl(int64 b, int n);    //左シフト b<<n
 int  lxcmp(int64 x, int64 y);  //大小判定
 int64 lxadd(int64 x, int64 y);  //加算 x+y
 int64 lxsub(int64 x, int64 y);  //減算 x-y
 int64 lxmul(int64 x, int64 y);  //乗算 x*y
 int64 lxdivs(int64 n, unsigned short d, int64 *r); //除算 n/d, r=余り
 int64 lxdiv(int64 n, int64 d, int64 *r); //除算 n/d, r=余り
 int64 lxdiv0(int64 n, int64 d, int64 *r); //除算 n/d, r=余り
 
実行時間
 lxmul()  :32ビット乗算の約10倍
 lxdiv()  :64ビット除算の約 4倍(VCの_int64との比較)[改良版] 
 lxdiv0() :64ビット除算の約20倍 回復型除算(〃)[改良版]
 ※速度を要求される計算には不向き 
  数値を改良後のものに変更12/5(商と剰余両方の計算時間で比較)

(参考)
64ビット整数構築手順
・64ビット整数を格納する構造体を定義する.
・大小比較演算を定義する.
・論理演算を定義する.(xor, or, and, notなど)
・シフト演算を定義する.
・加算を定義する.
・減算を定義する.
・乗算を定義する.
 32ビットの各種演算を利用
・unsigned shortによる除算を定義する.
 32ビットの各種演算と
 定義済の64ビット演算(シフト,加,減,乗)を利用
・64ビット整数による除算を定義する.
・除算を定義する.
 32ビットの各種演算と定義済の64ビット演算を利用
 
64ビット除算の定義(正数の場合の概要)
0.ゼロ除算エラー処理.
1.分母が16ビットになるまで,分子分母とも同じ数だけ右シフトする.
2.これに上記unsigned shortによる除算を適用して,概算値を求める.
  概算値=(分子>>n)/(分母>>n)  (概算値は真値より大きい)
  暫定剰余=分子-概算値×分母
  補正値=概算値+(暫定剰余>>n)/(分母>>n)
  剰余=分子-補正値×分母
3.剰余が依然として,分母より大きい場合は,再度補正を行う.
  2度の補正で,誤差は1以下に収束する.
4.最終段階で,+1,-1の調整を行う.
 ※この関数ではループを使わない.
 
[HP200LX/LXDIC]

|

2014年11月25日 (火)

Sunrise::仕様確定

■Sunrise::仕様確定
Ver 9.9
固定小数を用い,一定の計算精度を保ちながら,同時に計算も高速化する.
 
変更点
・浮動小数点数による計算を中止する.
・事前計算データ表参照による補正計算方式を廃止する.
・計算には,小数部19ビットの固定小数点数を用いる.
・高緯度地域は対応地域から除外する.
 
仕様概要
・対応地域
  東経: -180°~ 180°
  北緯:  -65°~ 65°
  時差:  -12 ~ +12  (日本は+9:JST=UTC+9)
・計算誤差:日出日入時刻:1分以内(緯度-65~65°)
      方位,高度 :1度以内(緯度-65~65°)
・表示項目
  場所(緯度,経度)
  日付,時刻
  日出時刻,日出方位
  南中時刻,南中高度
  日入時刻,日入方位
  太陽高度,太陽方位(指定した時刻の値)
  日出日入一覧表画面の表示
  グラフ表示(高度の時間変化,高度・方位関係)
・設定ファイルの位置情報登録数:10地点
 
計算
・プログラム内蔵の太陽赤緯/均時差テーブルから得たデータに,暦と
 季節のずれに基づく補正を加える.
・この赤緯/均時差,及び,位置情報を計算式に代入し,時角を求める.
 時角から,日出日入南中時刻を計算する.
・得られた時刻における方位角および高度を計算する.
・三角関数の計算は,小数部19ビットの固定小数点数で行う.
 三角関数値は,内蔵の三角関数値表から,内挿計算により算出する.
・固定小数点数の計算は,全て32ビット整数の演算で行う.
・グラフ用データの計算は,日出から日入までの時間を50等分して,
 各時刻における方位と高度を計算する.
・計算は現行の暦に基づいて行う.
 
コンパイラー
・Borland C++ 3.1(PC上でコンパイル)
 
表示速度の実測データ
(メイン画面)
 ・データ+グラフ表示時間:0.45秒(20件連続で9.0秒)
 ・データ+方位図表示時間:0.30秒(20件連続で6.0秒)
(データ一覧画面)
 ・ページデータの更新時間:0.38秒(20ページ更新:7.5秒)
 
[HP200LX/LXDIC]

|

2014年11月19日 (水)

Sunrise::固定小数点数の乗算

■固定小数点数の乗算
Sunrise
 
・固定小数点数の乗算も,小数部19ビットの場合は,計算の途中でオーバーフ
 ローが発生する.
   a*b ==> (a*b)>>19
 a=0x80000 , b=0x80000 の場合(実数では a=1.0, b=1.0)
  a*b = 0x40 0000 0000  : 39ビット(32ビットを越える)
  (a*b)>>19  : 20ビット(結果は20ビット以下になる)
 三角関数値の演算は,1.0以下の数値の演算なので,必ず 1.0以下の結果が得
 られる.すなわち,20ビット以下になる.
 
 したがって,64ビット整数をつかわなくても,工夫すれば,32ビット整数で
 高速に計算することができる.
 
 実際の計算は,32ビット整数を16ビット整数に分割して行う.
 //符号なし32ビット整数 a ==> (a1,a2): a1, a2 は32ビット整数型
  a1 = a>>16;         //a を16+16ビットに分割
  a2 = a&0x0000ffff; 
  b1 = b>>16;         //b を16+16ビットに分割
  b2 = b&0x0000ffff; 
 //積(a1, a2)*(b1,b2) = (a1*b1, a2*b1+a1*b2, a2*b2)
  c1 = (a1*b1);       //32--64ビット:c1,c2,c3は32ビット整数型
  c2 = (a2*b1+a1*b2); //16--48ビット:c1,c2,c3はオーバーフローしない※
  c3 = (a2*b2);       // 0--32ビット
 //c1,c2,c3を19ビット右シフト(c1,c2,c3はベースが16ビットずつ異なる)
  c  = (c1<<(32-19))+(c2>>(19-16))+(c3>>19);
 (演算回数) 乗算 4回,加算 3回,シフト 5回
 
※通常は,c2で発生する桁上りの処理が必要になるが,Sunriseの計算では必
 要ない.
 最大値をとる場合 a=0x80000, b=0x80000(==>1.0,1.0)を考える.
   a1=0x0008, a2=0x0000
   b1=0x0008, b2=0x0000
 a1, b1が常に小さいので,c2=a1*b2+a2*b1 が32ビット以上になる事はない.
 a2, b2が最大値0xffffをとる場合でも,a1*b2, a2*b1は0x7fff8:19ビット
 
[HP200LX/LXDIC]
 

|

2014年11月18日 (火)

Sunrise::固定小数点数の除算

■固定小数点数の除算
Sunrise
固定小数点数の計算では,小数部のビット数を大きくするほど,計算誤差を
小さくできる.(Sunriseの計算では19ビット)
しかし,32ビット整数の計算では,小数部を16ビット以上にすると,計算の
途中でオーバーフローが発生する.
 小数部19ビットの場合の除算:
   a/b ==> (a<<19)/b
 a=0x80000 の場合(実数では a=1.0 に相当)
  a ==> a<<19   20ビット ==> 39ビット(32ビットを越える)
  1000 0000 0000 0000 0000 = 0x80000
  0100 0000 0000 0000 0000 0000 0000 0000 0000 0000 = 0x4000000000
したがって,正しい計算結果を得るためには,特別の工夫が必要になる.
 
試した主な方法
1.分子と分母を同じ数で割って,小さくしてから計算する.
  (a<<19)/b ===>  ((a<<19)/128)/(b/128) = (a<<12)/(b>>7)
  計算は高速だが,分母で切捨てたビット数分だけ誤差が拡大する.
  
2.64ビット整数で計算する.
  HP200LXの環境では,64ビット整数は扱えないので,独自の演算関数を自
  作する必要があるが,単純な低速アルゴリズムを用いると,演算速度が
  浮動小数点数よりも数倍遅くなってしまい,まったく使いものにならない.
 
3.1の計算結果を補正する.(採用した方法)
  商の概算値を q0 とする.q0=(a<<12)/(b>>7)
  概算値と真の値の差を b倍して考える.正しい商 q=(a<<19)/b
     d = b*q0-(a<<19) (d>0)
  ここで,b*q0 ,(a<<19) は共に32ビットを越えるが,d が32ビットを越
  えることはない.これは,概算値と真の値の差が小さいことによる.
  したがって,d の値は途中のオーバーフローにもかかわらず,正しい値
  になっている.(計算例参照)
  商の真の値 q は,
    q = q0-d/b    となる.(q<q0)
  (これは,正確な値で,64ビット整数による計算値と一致する)
  
 (まとめ)固定小数点数の除算  a/b (小数部19ビット)
    q0=(a<<12)/(b>>7)
    d = b*q0-(a<<19)  (d>0)
    q = q0-d/b
 (32ビット整数の演算回数)乗算1回,除算2回,減算2回,シフト3回
  
●計算例 0x80000=524288
  a=524012=0x7FEEC,  b=524280=0x7FFF8
  真値  a/b ==>   0x7FEEC*0x80000/0x7FFF8
            = 0x7FEF3 = 524019 ※
  概算値 a/b ==> q0 = (0x7FEEC<<12)/(0X7FFF8>>7)
            = 0x7FEEC000/0x0FFF
            = 0x7FF6B = 524139 
            (120の誤差は無視できない)
  計算値 d=q0*b-(a<<19)
            = 0x7FF6B*0x7FFF8-(0x7FEEC<<19)
            = 0x3FFB1804AB-0x3FF7600000(上位3桁が同じ)
            (オーバーフロー部分は,両方とも3Fで等しい)
            = 0x03B804AB
       d/b    = 0x03B804AB/0x7FFF8
            = 0x77 = 119
      q = q0-d/b = 0x7FF6B-0x0077
            = 0x7FEF4 = 524020 ※
            (誤差:120/524288 ==> 1/524288)
  ※実数の計算では,524019.99591 になるので,どちらも正しい.
   (切上げか切捨ての違い)
         
[HP200LX/LXDIC]

|

«Sunrise::三角関数の整数演算化