Python vs Modula-2

(4) 太陽はどこを動いているか?
Where is the Sun traversing?

Workbench Python に戻る    (5)に進む

＜Task＞

地球上の任意の地点 O （緯度α）から天空上の太陽の位置 P（ψ(東西方向),φ(仰角)）を追いかけます。ψ=0 は真南でψ<0は東の方向、φ=0 は水平面上の高さとします。もし地球が球体であり、地球が太陽の周りを円運動すると仮定れば P の動きは幾何学で求められ、id=1〜365, jr=0〜Nphis-1 の二重配列で表せます。

ψ-φ のグラフは既に「太陽はどこに見える?」で求めてあります (Modula-2)。 P は球殻上を東から西に向かいますが、この球殻から遠く離れた点 Q から P の動きを眺めれば季節感の感じ方が変わってくるかもしれません。

実は、この様子は太陽の動きの「見える化」で既に動画にしてありますが、Modula-2のグラフィックスが元々プロッター対応のため残念ながらクオリティがよくありません。そこでPythonで同等のことをやってみようというわけです。まずは ψ-φ 図を作成します。

太陽の位置 (ψ, φ) を求める（北緯35°）
 
30日刻み
7日刻み
1日刻み
＜Python (Design A)＞

  • モジュールは3層
  • 最下層 Stglobal では変数の設定、関数の定義(回転行列の生成とベクトルの回転)を行う。
  • 次の層 STpplictionsでは、関数 STpsiphi(day,time)を定義。
  • 最上層 SunTracker (MAIN) では(ψ,φ)をSTpsiphiによって決定。そのあとグラフィックス表示。
  • 配列ψ[i,j],φ[i,j]の大きさを設定してから実行。
＜Python (Design B)＞

  • モジュールは2層
  • 下層 Stglobal は関数の定義(回転行列の生成とベクトルの回転)を行う。
  • 上層 SunTrackerEx (MAIN) は STpsiphi(day,time)を含む。 class getPsiPhi の中で配列ψ[i,j],φ[i,j]を生成。そのあとグラフィックス表示。
＜Modula-2＞

  • モジュールは3層
  • 最下層 Stglobal は　主に型の定義、変数の宣言、関数の定義(回転行列の生成とベクトルの回転、ψ(t)とφ(t)の計算)を行う。配列[ψ,φ]は余裕をもって定義する。
  • 次の層 STpplictionsでは、関数 STpsiphi(day,time)を定義。。
  • 最上層 SunTracker (MAIN) ではパラメータの設定、配列ψ[i,j],φ[i,j]の決定とグラフィックスが主な作業。
＜Python＞朱でψとφを計算

[Main] sunTracker.py 

import matplotlib.pyplot as plt
import numpy as np
import STglobal as STg
import STapplications as STa

Omegad0 = STg.Omegad0
Omegad0new = STa.sunTrack(Omegad0)

Nx = np.zeros(STg.Nis)
pd = np.pi/180.0
vrQ = np.zeros(3)
alpha  = STg.alphad*np.pi/180.0

def plotPsiPhi():            
    for i in range(STg.Nis):
        x = []
        y = []
        
        k = 0
        for j in range(STg.Nphis):
            if STg.phi2[i,j] > 0.0:
                x.append(STg.psi2[i,j])
                y.append(STg.phi2[i,j])
                k += 1
        Nx[i] = k
	.......
[Module] STapplications.py
def STpsiphi(id, jr, Omegad0): #, delta, alpha):
    alpha = STg.alphad*np.pi/180.0
    DTheta = 360.0/float(STg.Nphis)
    Omegad = float(id) + DTheta*float(jr)/360.0 
    Omegad = Omegad*360.0/STg.days_in_1yr + Omegad0
    Omega = Omegad*np.pi/180.0
    Ry = STg.makeRoty(delta)
    if STg.isClockwise:
        Rz = STg.makeRotz(-Omega)
        Rzinv = STg.makeRotz(Omega)
    else:
        Rz = STg.makeRotz(Omega)
        Rzinv = STg.makeRotz(-Omega)
            
    Thetad = DTheta*float(jr) - 180.0
    Theta = Thetad*np.pi/180.0
        
    rQ = STg.initVector(alpha, Theta)
    rQ = STg.rotVector(rQ, Rz)
    rQ = STg.rotVector(rQ, Ry)
    rQ = STg.rotVector(rQ, Rzinv) 
     
    zz = STg.polarisVector()

	......
    sQ = STg.vectorProduct(rQ, zz)

       
    rA = STg.incidentLightVector(0.0) 
    cc = STg.scalarProduct(rA, rQ) 

	......
    psi = STg.atan2(ca, cb)
    psid = psi*180.0/np.pi
    phi = np.arcsin(cc)
    phid = phi*180.0/np.pi
    return psi, phi, psid, phid, Omega, Omegad  
......

def sunTrack(Omegad0):
    alpha = STg.alphad*np.pi/180.0

    i = 0
    ii = 0
    while(True):
        phdprev = -1.0
        for j in range(STg.Nphis):
            ps,ph,psd,phd,Om,Omd = STpsiphi(i,j,Omegad0)
	.......
＜Python (Object Oriented Programming)＞朱でψとφを計算

[Main] sunTrackerEx.py 

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

def STpsiphi(id, jr, Omegad0, alphad, deltad, Nphis): 
	......

class getPsiPhi:
    def __init__(self,istep,daysyr,Omegad0,alphad,deltad,Nphis):
        Nis = daysyr // istep + 1
        if Nphis==0:
            Nphis = 12
        alpha = alphad*np.pi/180.0
        delta = deltad*np.pi/180.0

        DThetad = 360.0/float(Nphis)
        DTheta = DThetad*np.pi/180.0

        i = 0
        ii = 0
        
        psida = np.arange(Nis)
        psidb = np.arange(Nis)
        psi2 = np.zeros((Nis, Nphis))
        phi2 = np.zeros((Nis, Nphis))

        while(True):
            phdprev = -1.0
            for j in range(Nphis):
                ps, ph, psd, phd, Om, Omd 
		= STpsiphi(i,j,Omegad0,alphad,deltad,Nphis)
                if ii<Nis:
                    psi2[ii,j] = psd
                    phi2[ii,j] = phd
                if phd*phdprev < 0.0:
                    if phd > 0.0:
                        psida[ii] = psd
                    else:
                        psidb[ii] = psd
                phdprev = phd 
                
            ii += 1   
            if (i + istep) >= daysyr:
	            break
            else:
                i += istep 
                
        if Omd>1. :
            Omd -= 360.0
        self.Omegad = Omd
        self.Nis = Nis
        self.psida = psida
        self.psidb = psidb
        self.psi2 = psi2
        self.phi2 = phi2

def plotPsiPhi():
	......
	
w = getPsiPhi(daystep, daysyr, Omegad0, alphad, deltad, Nphis)
psi2 = w.psi2
phi2 = w.phi2 
Nis = w.Nis
  
plotPsiPhi()
＜Modula-2 (Structured Programming)＞朱でψとφを計算

[MAIN] SunTracker.MOD 

  MODULE SunTracker
    ･･････
    RotationParams(mode);
    i:= 0;    (* for day number, INC = iskip or istep*)
    ii:= 0;   (* for arrays, INC = 1             *)
    LOOP      (* to scan days *)
      phidprev:= -1.0;
      FOR j:= 0 TO Nphis DO　(* Dtheta = 360/Nphis *)
        STpsiphi(psid,phid,psi,phi,i,j,delta,alpha,DTheta);
        phione[j]:= phid; psione[j]:= psid;
        phi2[ii,j]:= phid; psi2[ii,j]:= psid;
        phidprev:= phid;
      END;   (* of FOR j:=0 TO Nphis-1   *)
 
      IF (i MOD imon = 0) THEN 
        IF toPlt3D THEN
          ChooseDataSequence(psione,phione,ndata,ii); 
           (* psione and phione redefined *)
          Make3Ddata(psione,phione,ndata,u,v,w); 
          ReverseArray(v,ndata);
          ReverseArray(w,ndata);
          One3Dplot(v,w,ndata,uu,vv,njj,0);
           (* Draw uu & vv in black, v & w in color *)
        END;
      END;   (* of IF (i MOD imon = 0) *)
	
      INC(ii); (* next dataset *)
	
      IF isLeap AND(istep=2)
	      AND (daysyr-1>i)AND((daysyr-i)<=istep) THEN
        i:= daysyr-1;       (* daysyr = 365 or 366 *)
      ELSIF (i+istep)>=daysyr THEN EXIT;
      ELSE INC(i,istep);
      END;
    END;    (* of LOOP, same as WHILE (i<=365) *)
[DEFINITION MODULE STglobal] 
VAR
  psi2,phi2,time2: ARRAY[0..Ahigh] OF 
		ARRAY[0..Angles] OF LONGREAL;

PROCEDURE MakeRotx(VAR xmatrix: ARRAY OF LongVector;
                           phi: LONGREAL);
PROCEDURE RotVector(VAR v: LongVector;
	                R: ARRAY OF LongVector);
PROCEDURE VectorProduct(VAR C: LongVector;
			 A, B: LongVector);

PROCEDURE ScalarProduct(A, B: LongVector): LONGREAL;
[DEFINITION MODULE STapplications] 
PROCEDURE STpsiphi(VAR psid: LONGREAL;
		   VAR phid: LONGREAL;
		    VAR psi: LONGREAL;
		    VAR phi: LONGREAL;
             id: INTEGER;   (* days 0..365      *)
             jr: INTEGER;   (* earth's rotation *)
    	  delta: LONGREAL;  (* inclination angle*)
     	  alpha: LONGREAL;  (* latitude         *)
         DTheta: LONGREAL); (* angle step       *)
コメント
この方式はModula-2の発想を受け継いでいます。オブジェクトは使っていません。 STglobal.py でいくつかの変数の値を定めましたが、Modula-2と違ってMain層から値を変更することはできません。 class getPsiPhiを使ってグラフィックス用のデータψ[i,j],φ[i,j]を作成します。制御パラメータはgetPsiPhiの入力とします。 最下層に制御パラメータと配列変数ψ[i,j],φ[i,j]を定義しています。ψとφはSTpsiphiからのVAR渡しで更新します。

5-22-2023, S. Hayashi