太陽の位置（方位角および仰角）を正確に算出するには、地理座標・日時情報に基づく天文学的計算が必要である。以下では、観測点の緯度・経度とUTC時刻を入力として、太陽の赤緯角（δ）および時角（T_s）を導出し、そこから方位角（A）と仰角（h）を求める一連の処理を示す。

基本公式

太陽の仰角 h および方位角 A は以下の式で与えられる：

h = arcsin( sin δ · sin φ + cos δ · cos φ · cos T<sub>s</sub> )
A = 
  { arccos( (sin h · sin φ − sin δ) / (cos h · cos φ) )      if T<sub>s</sub> ≤ 0
  { 360° − arccos( (sin h · sin φ − sin δ) / (cos h · cos φ) ) if T<sub>s</sub> > 0

ここで、φ は観測地の地理緯度（ラジアン）、δ は太陽赤緯角（ラジアン）、T_s は太陽時角（ラジアン）である。方位角 A は正北を基準（0°）とし、時計回りに増加する。

UTC時刻の取得

時刻情報はUTCで扱う必要がある。MATLAB環境では次のように現在時刻を取得できる：

now_utc = datetime('now', 'TimeZone', 'UTC');
year_val = year(now_utc);
month_val = month(now_utc);
day_val = day(now_utc);
hour_utc = hour(now_utc);
minute_utc = minute(now_utc);

days_in_month = [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
day_of_year = sum(days_in_month(1:month_val)) + day_val;

外部時間同期装置（時統端末）からタイムスタンプ（2000年1月1日0時からのUTC秒数）が得られる場合、次のように分解する：

timestamp = 766162398 + 12*3600; % 例：北京時間 2024/04/12 10:53:18

year_est = floor(timestamp / 31556736) + 2000;
day_of_year = floor(mod(timestamp, 31556736) / 86400) + 1;
hour_utc = floor(mod(timestamp / 3600, 24));
minute_utc = mod(timestamp / 60, 60);
decimal_hour = hour_utc + minute_utc / 60.0;

観測地の緯度設定

latitude_deg = 43.976850;  % 長春市の緯度（十進度）
phi = deg2rad(latitude_deg);  % ラジアンに変換
longitude_deg = 125.394173;  % 経度（東経）

太陽赤緯角 δ の計算

Bourges法による高精度な赤緯角計算を採用：

n0 = 79.676 + 0.2422 * (year_val - 1985) - round(0.25 * (year_val - 1985));
b = 2 * pi * (day_of_year - 1 - n0) / 365.2422;

delta_deg = 0.3723 + 23.2567 * sin(b) + 0.1149 * sin(2*b) ...
            - 0.1712 * sin(3*b) - 0.758 * cos(b) ...
            + 0.3656 * cos(2*b) + 0.0201 * cos(3*b);

delta = deg2rad(delta_deg);

真太陽時の計算（均時差補正）

標準時から真太陽時への変換には、均時差（EoT）と経度補正を適用：

eot = 0.0028 - 1.9587 * sin(b) + 9.9059 * sin(2*b) ...
      - 7.0924 * cos(b) - 0.6882 * cos(2*b);  % 単位：分

% 北京時間（UTC+8）を現地平均太陽時に変換
local_mean_time = (hour_utc + 8) + (minute_utc - (120 - longitude_deg) * 4) / 60.0;

% 真太陽時（分単位の均時差を加算）
solar_time = local_mean_time + eot / 60.0;

太陽時角 T_s の導出

ts_rad = (solar_time - 12.0) * pi / 12.0;  % ラジアン

仰角と方位角の最終計算

% 仰角（ラジアン → 度）
sun_elevation = asin(sin(delta) * sin(phi) + cos(delta) * cos(phi) * cos(ts_rad));
elev_deg = rad2deg(sun_elevation);

% 方位角（南基準 → 北基準に変換）
temp_az = acos((sin(sun_elevation) * sin(phi) - sin(delta)) / (cos(sun_elevation) * cos(phi)));
az_deg = 180 - rad2deg(temp_az);

% 午後（T<sub>s</sub> > 0）の場合は360°から減算
if ts_rad > 0
    az_deg = 360 - az_deg;
end

この方法により、実測値との誤差は0.5度以内となることが確認されている。