We track the sun from the point O of latitude α.
The position of the sun, P, is characterized by the rotation angle ψ in the horizontal east-west direction and the rotation angle φ rising from the horizon.
ψ=0 is directed toward south and ψ<0 is in the east. φ=0 is the horizon.

Actually the graph of ψ-φ has already been obtained in suntracker1.html (in Japanese, Modula-2). Here we adopt Python to obtain similar results.

• Two modules in layers
• Lower layser, Stglobal, defines the functions for rotation.
• Upper layer, SunTrackerEx (MAIN), includes STpsiphi(day,time).
  ψ[i,j],φ[i,j] are created within class getPsiPhi, followed by graphics display.

• Three modules in layers
• Lower layer,Stglobal, defines types, variables, functions, and arrays.
• Middle layer, STpplictions, defines the function STpsiphi(day,time).
• Upper layer, SunTracker (MAIN), determines arrays ψ[i,j] and φ[i,j], and displays graphics.
• Plotter output is shown above; screen output is of poor quality.

ψ and φ are calculated at red code.

[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()
		

ψ and φ are calculated at red code.

[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) *)

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;

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       *)