In [1]:

    

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.lines as mlines


    

In [47]:

    

def newline(p1, p2):
    ax = plt.gca()
    xmin, xmax = ax.get_xbound()

    if(p2[0] == p1[0]):
        xmin = xmax = p1[0]
        ymin, ymax = ax.get_ybound()
    else:
        ymax = p1[1]+(p2[1]-p1[1])/(p2[0]-p1[0])*(xmax-p1[0])
        ymin = p1[1]+(p2[1]-p1[1])/(p2[0]-p1[0])*(xmin-p1[0])

    l = mlines.Line2D([xmin,xmax], [ymin,ymax])
    ax.add_line(l)
    return l


    

In [2]:

    

#Declination of the sun, counted from Spring equinox
def sundec(day):
    return np.arcsin(np.sin(23.44/180*math.pi)*np.sin(day/365.24*2*math.pi))

def sundecdeg(day):
    return 180/math.pi*sundec(day)

#More accurate computation
def sundeq(day):
    dph = 3 #perihelion on Jan 3 (Jan 5, 2, 4, 4, 3, 4, 3, 3, 5, 2 in 2020, 2021, ...)
    w=2*math.pi/365.24
#    a=day*w+math.pi/2
    a=(day+89)*w
    b=a+2*0.0167*np.sin(w*(day+79-dph)) #solar ecliptic longitude counted from Winter solstice
    return -np.arcsin(np.sin(23.44/180*math.pi)*np.cos(b))


    

In [3]:

    

#Equation of time, counted from Spring equinox
def eqot(day):
    dph = 3 #perihelion on Jan 3 (Jan 5, 2, 4, 4, 3, 4, 3, 3, 5, 2 in 2020, 2021, ...)
    w=2*math.pi/365.24
#    a=day*w+math.pi/2
    a=(day+89)*w
    b=a+2*0.0167*np.sin(w*(day+79-dph)) #solar ecliptic longitude counted from Winter solstice
#    c=(a+np.arctan2(np.cos(b),np.cos(23.44/180*math.pi)*np.sin(b)))/math.pi+.5
    c=(a-np.arctan(np.tan(b)/np.cos(23.44/180*math.pi)))/math.pi
    return 12*(c-np.round(c))

#counting from Jan 1
def eqot1(d): #d=0 on Jan 1
    return eqot(d-79)

#counting from Jan 1
def eqot2(d): #d=0 on Jan 1
    dph = 3 #perihelion on Jan 3 (Jan 5, 2, 4, 4, 3, 4, 3, 3, 5, 2 in 2020, 2021, ...)
    w=2*math.pi/365.24
    a=(d+10)*w #counting from Winter solstice 
    b=a+2*0.0167*np.sin(w*(d+1-dph)) #solar ecliptic longitude counted from Winter solstice
    c=(a-np.arctan(np.tan(b)/np.cos(23.44/180*math.pi)))/math.pi
    return 12*(c-np.round(c))


    

In [4]:

    

#Equation of time
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi

x = np.linspace(0, 366, 500)
fig = plt.figure(figsize=(15,7))
ax = fig.gca()
plt.plot(x,60*eqot1(x))
plt.plot(x,60*eqot2(x),'r')
plt.margins(.01,.1)
#ax.set_xlim([-1,3])
#ax.set_ylim([-15,18])
ax.set_xticks(np.arange(0, 365, 20))
ax.set_yticks(np.arange(-14, 16.1, 2))
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
#plt.legend(loc=2)
#plt.savefig('scaphe.pdf',bbox_inches='tight')
plt.show()


    

In [5]:

    

#Declination curve for (unit) spherical sundial
x = np.linspace(0, 366, 500)
fig = plt.figure(figsize=(15,7))
ax = fig.gca()
plt.plot(x,sundecdeg(x))
plt.plot(x,180/math.pi*sundeq(x))
plt.margins(.01,.1)
#ax.set_xlim([-1,3])
ax.set_ylim([-1,10])
ax.set_xticks(np.arange(0, 365, 20))
ax.set_yticks(np.arange(-24, 25, 2))
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
#plt.legend(loc=2)
#plt.savefig('scaphe.pdf',bbox_inches='tight')
plt.show()


    

In [129]:

    

#Declination curve for polar sundial with (unit) cylindrical surface (gnomon at the centre)
x = np.linspace(0, 366, 500)
fig = plt.figure(figsize=(15,7))
ax = fig.gca()
plt.plot(x,np.tan(sundec(x)))
plt.plot(x,np.tan(sundeq(x)))
plt.margins(.01,.1)
#ax.set_xlim([-1,3])
#ax.set_ylim([-1,10])
ax.set_xticks(np.arange(0, 365, 20))
ax.set_yticks(np.arange(-0.5, 0.55, .1))
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
#plt.legend(loc=2)
#plt.savefig('sundecflat.pdf',bbox_inches='tight')
plt.show()


    

In [42]:

    

#Hour curve for polar sundial with flat surface (gnomon at unit distance)
x = np.linspace(7, 17, 100)
fig = plt.figure(figsize=(15,7))
ax = fig.gca()
plt.plot(x,np.tan((x-12)*math.pi/12))
plt.margins(.01,.1)
#ax.set_xlim([-1,3])
#ax.set_ylim([-1,10])
ax.set_xticks(np.arange(7, 17.1, 1))
#ax.set_yticks(np.arange(-0.4, 0.5, .1))
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
#plt.legend(loc=2)
#plt.savefig('sunhourflat.pdf',bbox_inches='tight')
plt.show()


    

In [77]:

    

#Hour marks for polar sundial with flat surface (gnomon at unit distance)
fig = plt.figure(figsize=(15,5))
ax = fig.gca()
d = .05
for i in range(7,18):
    x = np.tan((i-12)*math.pi/12)
    plt.plot([x,x],[0,1],'b')
    plt.text(x, 1+d, i, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
plt.margins(.01,.1)
#ax.set_xlim([-1,3])
ax.set_ylim([-0.1,1.1])
ax.set_xticks(np.arange(-4, 4.1, 1))
#ax.set_yticks(np.arange(-0.4, 0.5, .1))
ax.set_yticks([])
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
#plt.legend(loc=2)
#plt.savefig('sunhourflat.pdf',bbox_inches='tight')
plt.show()


    

In [143]:

    

#Seasonal and hourly marks for polar sundial with (unit) cylindrical surface (gnomon at the centre)
fig = plt.figure(figsize=(15,5))
ax = fig.gca()
d = .03
for i in range(-3,4):
    y = np.tan(sundec(i*365.5/12))
    plt.plot([-1.5,1.5],[y,y],'g')
for i in range(7,18):
    x = (i-12)*math.pi/12
    plt.plot([x,x],[-.5,.5],'b')
    plt.text(x, .5+d, i, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
plt.margins(.01,.1)
ax.set_xlim([-1.5,1.5])
ax.set_ylim([-0.5,0.5])
ax.set_xticks(np.arange(-1.5, 1.1, .5))
ax.set_yticks(np.arange(-0.5, 0.51, .1))
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
ax.set_aspect(1)
#plt.legend(loc=2)
#plt.savefig('sunhourflat.pdf',bbox_inches='tight')
plt.show()


    

In [130]:

    

#Seasonal and hourly marks for polar sundial with flat surface (gnomon at unit distance)
def sdpolar(x):
    return np.sqrt(x*x+1)*np.tan(delta)
def sspolar(x):
    return cta
    
fig = plt.figure(figsize=(15,8))
ax = fig.gca()
d = .1
lat = 47.8864
cta = 1/np.tan(lat*math.pi/180)
x = np.linspace(-4, 4, 100)
for i in range(-3,4):
    delta = sundeq(i*365.5/12)
    plt.plot(x,sdpolar(x),'g')
for i in range(7,18):
    x = np.tan((i-12)*math.pi/12)
    plt.plot([x,x],[-2,1],'b')
    plt.text(x, 1+d, i, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
plt.plot([-4,4],[cta,cta],'r')
cc = plt.Circle((0,0),.03,color='r')
ax.add_artist(cc)

plt.margins(.01,.1)
ax.set_xlim([-2,2])
ax.set_ylim([-2,1])
ax.set_xticks(np.arange(-4, 4.1, 1))
ax.set_yticks(np.arange(-2, 1.1, 1))
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
ax.set_aspect(1)
ax.xaxis.set_major_formatter(plt.NullFormatter())
ax.xaxis.set_ticks_position('none')
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.yaxis.set_ticks_position('none')
#plt.legend(loc=2)
#plt.savefig('polarsundial.pdf',bbox_inches='tight')
plt.show()


    

In [131]:

    

#Seasonal and hourly marks for horizontal sundial (gnomon at unit height at the origin)
def sdhoriz(x):
    td=np.tan(delta)
    t2d=td*td
    a = 1-t2a*t2d
    b = ta*(t2d+1)
    c = x*x*sc2a*t2d+t2d-t2a
    d = a*c+b*b
    v = (np.sqrt(d)*np.sign(delta)+b)/a
    return v

fig = plt.figure(figsize=(15,9))
ax = fig.gca()
d = .03
xm = 10
ym = 7
lat = 47.8864
h = np.tan(lat*math.pi/180)
h0 = 1/h
h1 = ym
ua = (h0+h1)*np.sin(lat*math.pi/180)

ta=h
t2a=ta*ta
ca=np.cos(lat*math.pi/180)
sc2a=1/(ca*ca)


x = np.linspace(-xm, xm, 100)
for i in range(-3,4):
    delta = sundeq(i*365.5/12)
    plt.plot(x,sdhoriz(x),'g')

for i in range(7,18):
    x = ua*np.tan((i-12)*math.pi/12)
    plt.plot([0,x],[-h0,h1],'b')
cc = plt.Circle((0,0),.08,color='r')
ax.add_artist(cc)

for i in range(9,16):
    x = ua*np.tan((i-12)*math.pi/12)
    plt.text(x, h1+.2, i, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')


plt.margins(.01,.1)
ax.set_xlim([-xm,xm])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(-xm, xm+.1, 1))
ax.set_yticks(np.arange(-1, ym+.1, 1))
plt.rcParams['grid.linestyle'] = "--"
ax.set_aspect(1)
plt.grid()
ax.xaxis.set_major_formatter(plt.NullFormatter())
ax.xaxis.set_ticks_position('none')
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.yaxis.set_ticks_position('none')
#plt.legend(loc=2)
#plt.savefig('horizontalsundial.pdf',bbox_inches='tight')
plt.show()


    

In [290]:

    

#Seasonal and hourly marks for horizontal sundial (gnomon at unit height at the origin)
def sdhx(t):
    td=np.tan(delta)
    return np.sin(t)/(csa*np.cos(t)-sna*td)
def sdhv(t):
    td=np.tan(delta)
    return (sna*np.cos(t)+csa*td)/(csa*np.cos(t)-sna*td)

fig = plt.figure(figsize=(15,9))
ax = fig.gca()
d = .03
xm = 10
ym = 7
tm = 4*math.pi/12
lat = 47.8864
h = np.tan(lat*math.pi/180)
h0 = 1/h
h1 = ym
ua = (h0+h1)*np.sin(lat*math.pi/180)

sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)


for i in range(-3,4):
    delta = sundec(i*365.5/12)
    tm = np.arccos(h*np.tan(delta))-.5*math.pi/12
    t = np.linspace(-tm, tm, 100)
    plt.plot(sdhx(t),sdhv(t),'g')

for i in range(7,18):
    x = ua*np.tan((i-12)*math.pi/12)
    plt.plot([0,x],[-h0,h1],'b')
cc = plt.Circle((0,0),.08,color='r')
ax.add_artist(cc)

for i in range(9,16):
    x = ua*np.tan((i-12)*math.pi/12)
    plt.text(x, h1+.2, i, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')


plt.margins(.01,.1)
ax.set_xlim([-xm,xm])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(-xm, xm+.1, 1))
ax.set_yticks(np.arange(-1, ym+.1, 1))
plt.rcParams['grid.linestyle'] = "--"
ax.set_aspect(1)
plt.grid()
ax.xaxis.set_major_formatter(plt.NullFormatter())
ax.xaxis.set_ticks_position('none')
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.yaxis.set_ticks_position('none')
#plt.legend(loc=2)
#plt.savefig('horizontalsundial1.pdf',bbox_inches='tight')
plt.show()


    

In [132]:

    

#Seasonal and hourly marks for south facing vertical sundial (gnomon at unit distance at the origin)
def sdvert(x):
    td=np.tan(delta)
    t2d=td*td
    a = 1-t2a*t2d
    b = ta*(t2d+1)
    c = x*x*sc2a*t2d+t2d-t2a
    d = a*c+b*b
    v = (np.sqrt(d)*np.sign(delta)+b)/a
    return v

fig = plt.figure(figsize=(15,9))
ax = fig.gca()
d = .03
xm = 10
ym = 7
lat = 47.8864
colat = 90-lat
h = np.tan(colat*math.pi/180)
h0 = 1/h
h1 = ym
ua = (h0+h1)*np.sin(colat*math.pi/180)

ta=h
t2a=ta*ta
ca=np.cos(lat*math.pi/180)
sc2a=1/(ca*ca)


x = np.linspace(-xm, xm, 100)
for i in range(-3,4):
    delta = sundeq(i*365.5/12)
    plt.plot(x,-sdvert(x),'g')

for i in range(7,18):
    x = ua*np.tan((i-12)*math.pi/12)
    plt.plot([0,x],[h0,-h1],'b')

plt.plot([-xm,xm],[0,0],'r')

cc = plt.Circle((0,0),.08,color='r')
ax.add_artist(cc)

for i in range(8,17):
    x = ua*np.tan((i-12)*math.pi/12)
    plt.text(x, -h1-.2, i, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')


plt.margins(.01,.1)
ax.set_xlim([-xm,xm])
ax.set_ylim([-ym,1.3])
ax.set_xticks(np.arange(-xm, xm+.1, 1))
ax.set_yticks(np.arange(-ym, 1.3, 1))
plt.rcParams['grid.linestyle'] = "--"
ax.set_aspect(1)
plt.grid()
ax.xaxis.set_major_formatter(plt.NullFormatter())
ax.xaxis.set_ticks_position('none')
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.yaxis.set_ticks_position('none')
#plt.legend(loc=2)
#plt.savefig('verticalsundial.pdf',bbox_inches='tight')
plt.show()


    

In [133]:

    

#Solar altitude
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi
def setrise(d):
    t=np.arccos(tna*np.tan(d))
    return a2h(t)
def solalt(tt):
    td=np.tan(delta)
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solaz(tt):
    td=np.tan(delta)
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)
def solaltd(d):
    td=np.tan(d)
    t=h2a(ttt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solazd(d):
    td=np.tan(d)
    t=h2a(ttt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)

d = .03
xm = 10
ym = 70
t1=4
t2=20

lat = 47.780622
lon = 106.644736

sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)
tna=np.tan(lat*math.pi/180)

fig = plt.figure(figsize=(15,9))
ax = fig.gca()

t = np.linspace(t1, t2, 100)

for i in range(-3,4):
    delta = -sundeq(i*365.24/12)
    plt.plot(t,solalt(t),'g')

plt.margins(.01,.1)
ax.set_xlim([t1,t2])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(t1, t2+.1, 1))
ax.set_yticks(np.arange(0, ym+.1, 10))
plt.rcParams['grid.linestyle'] = "--"
#ax.set_aspect(1)
plt.grid()
#plt.legend(loc=2)
#plt.savefig('solaraltitude.pdf',bbox_inches='tight')
plt.show()


    

In [134]:

    

#Solar azimuth
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi
def setrise(d):
    t=np.arccos(tna*np.tan(d))
    return a2h(t)
def solalt(tt):
    td=np.tan(delta)
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solaz(tt):
    td=np.tan(delta)
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)

d = .03
xm = 10
ym = 100
t1=4
t2=20

lat = 47.780622
lon = 106.644736

sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)
tna=np.tan(lat*math.pi/180)

fig = plt.figure(figsize=(15,9))
ax = fig.gca()

t = np.linspace(t1, t2, 100)

for i in range(-3,4):
    delta = -sundeq(i*365.24/12)
    plt.plot(t,solaz(t),'g')

plt.margins(.01,.1)
ax.set_xlim([t1,t2])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(t1, t2+.1, 1))
ax.set_yticks(np.arange(-ym, ym+.1, 10))
plt.rcParams['grid.linestyle'] = "--"
#ax.set_aspect(1)
plt.grid()
#plt.legend(loc=2)
#plt.savefig('solarazimuth.pdf',bbox_inches='tight')
plt.show()


    

In [135]:

    

#Solar altitude-azimuth curve
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi
def setrise(d):
    t=np.arccos(tna*np.tan(d))
    return a2h(t)
def solalt(tt):
    td=np.tan(delta)
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solaz(tt):
    td=np.tan(delta)
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)
def solaltd(d):
    td=np.tan(d)
    t=h2a(ttt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solazd(d):
    td=np.tan(d)
    t=h2a(ttt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)

d = .03
xm = 130
ym = 72
t1=4
t2=20

lat = 47.780622
lon = 106.644736

sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)
tna=np.tan(lat*math.pi/180)

fig = plt.figure(figsize=(15,7))
ax = fig.gca()

t = np.linspace(t1, t2, 100)

for i in range(-3,4):
    delta = -sundeq(i*365.24/12)
    plt.plot(solaz(t),solalt(t),'g')

k=2.5
for ttt in range(5,20):
    d = np.linspace(sundeq(-3*365.5/12),sundeq(3*365.5/12),100)
    plt.plot(solazd(d),solaltd(d),'r')
    d = -sundeq(2*365.24/12)
    dx = solazd(d)
    dy = solaltd(d)
    d = -sundeq(3*365.24/12)
    x = solazd(d)
    y = solaltd(d)
    dx = x-dx
    dy = y-dy
    d = np.sqrt(dx*dx+dy*dy)
    dx = dx/d
    dy = dy/d
    plt.text(x+k*dx,y+k*dy, ttt, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
    
plt.margins(.01,.1)
ax.set_xlim([t1,t2])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(-xm, xm+.1, 10))
ax.set_yticks(np.arange(0, ym+.1, 10))
plt.rcParams['grid.linestyle'] = "--"
ax.set_aspect(1)
plt.grid()
#plt.legend(loc=2)
#plt.savefig('azalt.pdf',bbox_inches='tight')
plt.show()


    

In [37]:

    

#Solar height-azimuth curve
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi
def setrise(d):
    t=np.arccos(tna*np.tan(d))
    return a2h(t)
def solalttg(tt):
    td=np.tan(delta)
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.tan(alt)
def solaz(tt):
    td=np.tan(delta)
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return az
def solalttgd(d):
    td=np.tan(d)
    t=h2a(ttt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.tan(alt)
def solazd(d):
    td=np.tan(d)
    t=h2a(ttt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return az

d = .03
xm = 2.5
ym = 2.3
t1=4
t2=20

lat = 47.780622
lon = 106.644736

sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)
tna=np.tan(lat*math.pi/180)

fig = plt.figure(figsize=(20,15))
ax = fig.gca()

t = np.linspace(t1, t2, 100)

for i in range(-3,4):
    delta = -sundeq(i*365.24/12)
    plt.plot(solaz(t),solalttg(t),'g')

kx=.05
ky=.05
for ttt in range(5,20):
    d = np.linspace(sundeq(-3*365.5/12),sundeq(3*365.5/12),200)
    plt.plot(solazd(d),solalttgd(d),'r')
    d = -sundeq(2*365.24/12)
    dx = solazd(d)
    dy = solalttgd(d)
    d = -sundeq(3*365.24/12)
    x = solazd(d)
    y = solalttgd(d)
    dx = x-dx
    dy = y-dy
    d = np.sqrt(dx*dx+dy*dy)
    dx = np.sign(dx)
    plt.text(x+kx*dx,y+ky, ttt, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
    
plt.margins(.01,.1)
#ax.set_xlim([t1,t2])
ax.set_ylim([0,ym])
ax.set_xticks(np.arange(-xm, xm+.1, .5))
ax.set_yticks(np.arange(0, ym+.1, .5))
plt.rcParams['grid.linestyle'] = "--"
ax.set_aspect(1)
plt.grid()
#plt.legend(loc=2)
#plt.savefig('azalttg.pdf',bbox_inches='tight')
plt.show()


    

In [112]:

    

#Equation of time
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi

x = np.linspace(0, 366, 500)
fig = plt.figure(figsize=(15,7))
ax = fig.gca()
plt.plot(x,60*eqot1(x))
plt.plot(x,60*eqot2(x),'r')
plt.margins(.01,.1)
#ax.set_xlim([-1,3])
ax.set_ylim([-15,18])
ax.set_xticks(np.arange(0, 365, 20))
ax.set_yticks(np.arange(-14, 16.1, 2))
plt.rcParams['grid.linestyle'] = "--"
plt.grid()
#plt.legend(loc=2)
#plt.savefig('scaphe.pdf',bbox_inches='tight')
plt.show()


    

In [9]:

    

#Solar altitude-azimuth curve with equation of time accounted for
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi
def setrise(d):
    t=np.arccos(tna*np.tan(d))
    return a2h(t)
def solalt(tt):
    d=-sundeq(ddd)
    td=np.tan(d)
    t=h2a(tt+eqot(ddd)+lon/15-timezone-dst)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solaz(tt):
    d=-sundeq(ddd)
    td=np.tan(d)
    t=h2a(tt+eqot(ddd)+lon/15-timezone-dst)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)
def solaltd(day):
    d=-sundeq(day)
    td=np.tan(d)
    tt=ttt+eqot(day)+lon/15-timezone-dst
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solazd(day):
    d=-sundeq(day)
    td=np.tan(d)
    tt=ttt+eqot(day)+lon/15-timezone-dst
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)

xm = 130
ym = 72
t1=4
t2=20

lat = 47.780622
lon = 106.644736
timezone = 8
dst = 0

sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)
tna=np.tan(lat*math.pi/180)

fig = plt.figure(figsize=(20,8))
ax = fig.gca()

t = np.linspace(t1, t2, 100)

for i in range(-3,4):
    ddd=i*365.24/12
    plt.plot(solaz(t),solalt(t),'b')

k=2.5
for ttt in range(6,21):
    d = np.linspace(-3*365.24/12,3*365.24/12,100)
    plt.plot(solazd(d),solaltd(d),'r')
    d = np.linspace(3*365.24/12,9*365.24/12,100)
    plt.plot(solazd(d),solaltd(d),'g')
    d = -3*365.24/12
    dx = solazd(d)
    dy = solaltd(d)
    d = 3*365.24/12
    x = solazd(d)
    y = solaltd(d)
    dx = x-dx
    dy = y-dy
    d = np.sqrt(dx*dx+dy*dy)
    dx = dx/d
    dy = dy/d
    plt.text(x+k*dx,y+k*dy, ttt, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
    
plt.margins(.01,.1)
ax.set_xlim([t1,t2])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(-xm, xm+.1, 10))
ax.set_yticks(np.arange(0, ym+.1, 10))
plt.rcParams['grid.linestyle'] = "--"
#ax.set_aspect(1)
plt.grid()
#plt.legend(loc=2)
#plt.savefig('azalteqot.pdf',bbox_inches='tight')
plt.show()


    

In [14]:

    

60-(eqot(90)+lon/15+1-timezone)*60


    

    
Out[14]:





54.492702130592434

In [15]:

    

60-(eqot(-90)+lon/15+1-timezone)*60


    

    
Out[15]:





51.07666859966933

In [7]:

    

#Solar altitude-azimuth curve with equation of time accounted for
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi
def setrise(d):
    t=np.arccos(tna*np.tan(d))
    return a2h(t)
def solalt(tt):
    td=np.tan(delta)
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solaz(tt):
    td=np.tan(delta)
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)
def solaltd(day):
    d=-sundeq(day)
    td=np.tan(d)
    tt=ttt+eqot(day)+lon/15-timezone-dst
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solazd(day):
    d=-sundeq(day)
    td=np.tan(d)
    tt=ttt+eqot(day)+lon/15-timezone-dst
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)

xm = 150
ym = 59
t1=2.5
t2=21.5

#Moscow
#lat = 55.755833
#lon = 37.617222
#timezone = 3
#dst = 0

#Helsinki
lat = 60.170833
lon = 24.9375
timezone = 2
dst = 0

sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)
tna=np.tan(lat*math.pi/180)

fig = plt.figure(figsize=(20,8))
ax = fig.gca()

t = np.linspace(t1, t2, 100)

for i in range(-3,4):
    delta = -sundeq(i*365.24/12)
    plt.plot(solaz(t),solalt(t),'b')

k=2.5
for ttt in range(4,22):
    d = np.linspace(-3*365.24/12,3*365.24/12,100)
    plt.plot(solazd(d),solaltd(d),'r')
    d = np.linspace(3*365.24/12,9*365.24/12,100)
    plt.plot(solazd(d),solaltd(d),'g')
    d = -3*365.24/12
    dx = solazd(d)
    dy = solaltd(d)
    d = 3*365.24/12
    x = solazd(d)
    y = solaltd(d)
    dx = x-dx
    dy = y-dy
    d = np.sqrt(dx*dx+dy*dy)
    dx = dx/d
    dy = dy/d
    plt.text(x+k*dx,y+k*dy, ttt, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
    
plt.margins(.01,.1)
ax.set_xlim([t1,t2])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(-xm, xm+.1, 10))
ax.set_yticks(np.arange(0, ym+.1, 10))
plt.rcParams['grid.linestyle'] = "--"
#ax.set_aspect(1)
plt.grid()
#plt.legend(loc=2)
#plt.savefig('azalteqot.pdf',bbox_inches='tight')
plt.show()


    

In [9]:

    

#Solar altitude-azimuth curve with equation of time accounted for
def h2a(h):
    return (h-12)*math.pi/12
def a2h(a):
    return 12+a*12/math.pi
def setrise(d):
    t=np.arccos(tna*np.tan(d))
    return a2h(t)
def solalt(tt):
    td=np.tan(delta)
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solaz(tt):
    td=np.tan(delta)
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)
def solaltd(day):
    d=-sundeq(day)
    td=np.tan(d)
    tt=ttt+eqot(day)+lon/15-timezone-dst
    t=h2a(tt)
    xx=np.sin(t)
    vv=sna*np.cos(t)+csa*td
    alt=np.arctan2(csa*np.cos(t)-sna*td,np.sqrt(xx*xx+vv*vv))
    return np.degrees(alt)
def solazd(day):
    d=-sundeq(day)
    td=np.tan(d)
    tt=ttt+eqot(day)+lon/15-timezone-dst
    t=h2a(tt)
    az=np.arctan2(np.sin(t),sna*np.cos(t)+csa*td)
    return np.degrees(az)

xm = 120
ym = 92
t1=5
t2=19

#Moscow
#lat = 55.755833
#lon = 37.617222
#timezone = 3
#dst = 0

#Helsinki
#lat = 60.170833
#lon = 24.9375
#timezone = 2
#dst = 0

#Taipei
lat = 25.066667
lon = 121.516667
timezone = 8
dst = 0


sna=np.sin(lat*math.pi/180)
csa=np.cos(lat*math.pi/180)
tna=np.tan(lat*math.pi/180)

fig = plt.figure(figsize=(20,8))
ax = fig.gca()

t = np.linspace(t1, t2, 300)

for i in range(-3,4):
    delta = -sundeq(i*365.24/12)
    plt.plot(solaz(t),solalt(t),'b')

k=2.5
for ttt in range(6,19):
    d = np.linspace(-3*365.24/12,3*365.24/12,100)
    plt.plot(solazd(d),solaltd(d),'r')
    d = np.linspace(3*365.24/12,9*365.24/12,100)
    plt.plot(solazd(d),solaltd(d),'g')
    d = -3*365.24/12
    dx = solazd(d)
    dy = solaltd(d)
    d = 3*365.24/12
    x = solazd(d)
    y = solaltd(d)
    dx = x-dx
    dy = y-dy
    d = np.sqrt(dx*dx+dy*dy)
    dx = dx/d
    dy = dy/d
    plt.text(x+k*dx,y+k*dy, ttt, horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')
    
plt.margins(.01,.1)
ax.set_xlim([t1,t2])
ax.set_ylim([-1,ym])
ax.set_xticks(np.arange(-xm, xm+.1, 10))
ax.set_yticks(np.arange(0, ym+.1, 10))
plt.rcParams['grid.linestyle'] = "--"
ax.set_aspect(1)
plt.grid()
#plt.legend(loc=2)
#plt.savefig('azalteqot.pdf',bbox_inches='tight')
plt.show()