二つの角度(ψ, φ)を用いて太陽の位置を表す 太陽の位置の「見える化」をします。 南をx軸、東をy軸、天頂をz軸とし、目線の方向を偏角 α と方位角 β とします。座標軸を2回回転し、 一つの軸を目線に合わせて(ψ, φ)に対応した点 Q(u,v) を得ます。前回用いたプログラムに若干のQ の動きが太陽の動きとなります。

既にModula-2を用いた太陽の動きの「見える化」でこの作業ををやりましたが、画質がよくなかったのでPythonで同等のことをやってみようという訳です。

• モジュールは2層
• 下層 Stglobal は関数の定義(回転行列の生成とベクトルの回転)を行う。
• 上層 SunTrackerExx (MAIN) ではclass getPsiPhi の中で配列ψ[i,j],φ[i,j]を生成したあと座標回転によって軌跡 (xx,yy)を得る。時間ステップ毎に (xx,yy)を描く。 旧版の sunTrackerEx にclass rotMatrices を追加。このオブジェクトを座標軸（フレーム）と軌跡の座標回転で共有。 そのあとグラフィックス表示。

• モジュールは3層
• 最下層 STglobalEx は STglobal に座標軸の回転を行う Make3Ddata, 枠線を描く Prep3Dplot, データ曲線を描く One3Dplot を追加。この層の回転行列データをすべての座標回転で共有する。
• 最上層 SunTracker (MAIN) はパラメータの設定とグラフィックス表示を行う。
• グラフィックス画面の初期化を一度行う。古い配列(u,v)を無色で上書きして消去したあと、新しい配列(u,v)を描いてリフレッシュ。

[Main] sunTrackerExx.py

import matplotlib.pyplot as plt
import numpy as np
import STglobalEx as STg

class rotMatrices:
	......
class getPsiPhi:
	......
def rotAxes(ww, vrQ):
	......
def drawSingleTraj(valpha, vbeta, vextra, isel, Nphis, delay): 
    ww = rotMatrices(valpha, vbeta, vextra)   
    vrQ = np.zeros(3)
    def drawFrame():
        Nsec = 8
        Nseg = 36
        for j in range(2):            # axes
            u = []
            v = []
            k = 0
            for i in range(2):
                if j==0:
                    xi = i; yi = 0.0
                else:
                    xi = 0.0; yi = float(2*i-1)            
                zi = 0.0
                vrQ = STg.initVectorXYZ(xi, yi, zi)
                vrQ = rotAxes(ww,vrQ)
                u.append(vrQ[compx])
                v.append(vrQ[compy])

            plt.axis('off')    
            plt.plot(u, v, color="black", lw=0.5)
			......

          
    i = isel
    xx = []
    yy = []
    for j in range(Nphis):

        if phi2[i,j] > 0.0:
            drawFrame()
            
            xi = np.cos(psi2[i,j]*pd)*np.cos(phi2[i,j]*pd)
            yi = -np.sin(psi2[i,j]*pd)*np.cos(phi2[i,j]*pd)
            zi = np.sin(phi2[i,j]*pd)
            vrQ = STg.initVectorXYZ(xi, yi, zi)
            vrQ = rotAxes(ww,vrQ)

            xx.append(vrQ[compx])
            yy.append(vrQ[compy])
            plt.axis('off')    
            plt.plot(xx, yy, color="red", lw=2)
            if delay>0:
                plt.pause(delay)           

Nphis = 180
delay = 0.2
w = getPsiPhi(daystep, daysyr, Omegad0, alphad, deltad, Nphis)

psi2 = w.psi2
phi2 = w.phi2 
Nis = w.Nis

ax = plt.axes()
for i in range(Nis):
	SingleTraj(valpha, vbeta, vextra, i, Nphis, delay)
plt.show() 

[MAIN] SunTracker.MOD

MODULE SunTracker
    ･･････
PROCEDURE Plot3D(Ndevice: INTEGER);　(* 概要 *)
  BEGIN
    Prep3Dplot(alpha,beta, X,Y,101, "x", "y", Ndevice);
    ChooseDataSequence(psione,phione,ndata,ii);
    Make3Ddata(psione,phione,ndata,u,v,w);
    One3Dplot(v,w,ndata,uu,vv,njj,Ndevice);

FROM HgrWindows IMPORT iniplt3D, ...; 

PROCEDURE One3Dplot( VAR x,y: ARRAY OF LONGREAL;
		        nphi: INTEGER;
		   VAR xx,yy: ARRAY OF LONGREAL;
		   VAR nphii: INTEGER;
		     Ndevice: INTEGER);
 
  PROCEDURE Plot;
    VAR
      u, v: ARRAY [0..1] OF LONGREAL;
      i, j: INTEGER;
    BEGIN
      IF Ndevice=0 THEN
        IF nphii>0 THEN   (* erase previous curve *)
          polyline(xx, yy, nphii, nullcolor, 0, Size2, 0, 0);
        END;
        color:= col14;
        FOR j:= 2 TO nphi-1 DO
          polyline(x, y, j, color, 0, Size2, 0, 0);
          Lib.Delay(idling3D);
        END;	
        FOR i:=0 TO nphi-1 DO
          xx[i]:= x[i]; yy[i]:= y[i];
        END;
        nphii:= nphi;
      ELSE
        color:= 1;
        polyline(x, y, nphi, color, 0, Size2, 0, 0);
      END;			
    END Plot;

  BEGIN
    Plot;
    IF Ndevice=2 THEN endplt; END;
  END One3Dplot;