In the preceding part 4 the angles (ψ, φ) have been obtained, which characterizes the position of the sun in the sky as seen from the observer at a certain point O on the globe.

The situation may be understood as if the sun, a bright light source, were moving on the ceiling of a big dome. (You are standing in the center of the ground.)

Now let's assume that the ceiling is transparent and that you are lifted far beyond the ceiling. You can still observe the movement of the light source that simulates the sun.

We characterize the geometry with the origin located at O. Let the x-axis be directed toward south and y-axis toward east. z-axis is perpendicular to these two, forming a right-handed system..

Next you place a sheet of paper in front of you perpendicular to the (α, β) vector, and record the trajectory of the simulated sun, thus giving the drawing in two-dimensional (u, v) space. If possible, the movement is recorded in movie.

In fact the task has been done in suntracker2 (in Japanese, Modula-2), but we adopt Python here.

[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); (* outline *)
  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;