two sources of occlusion self occlusion, earth occlusion
earth occlusion when norm.solar<0 self occlusion when point.solar<0
how does how does the panel normal rotate. on the pole the normal is constant but a panel will be rotated throughout the course of a day. the tangent space is getting rotated, sorta on the equator the solar panel normal is unchanging
stationary planet, x,y,z axes
then any local vector pointing off of surface is defined by $$ \hat{r}\quad \hat{θ} \quad \hat{φ} $$ physical location: $$ θ_l $$ -pi/2 -z pi/2 +z panel orientation: $$ θ_p \quad φ $$ stationary planet panel normal: $$ \hat{z}sin(θ_l+θ_p) + cos(θ_l+θ_p)(\hat{x}cos(φ)+\hat{y}sin(φ)) $$ rotating normal $$ \hat{x}→ \hat{x}cos(ω_dt)+\hat{y}sin(ω_dt) \quad \hat{y}→ \hat{y}cos(ω_dt)-\hat{x}sin(ω_dt) $$ solar frame $$ \hat{y}→ \hat{y} \quad\hat{z}→ \hat{z}cos(23^ˆ)+\hat{x}sin(23^ˆ) \newline $$ $$ \hat{x}→ \hat{x}cos(23^ˆ)-\hat{z}sin(23^ˆ) $$
can then model sun vector as $$ \hat{x}cos(ω_yt)+\hat{y}sin(ω_yt) $$ or just $$ \hat{x}cos(θ_s)+\hat{y}sin(θ_s) $$
rotating frame panel normal: zsin(stheta)+cos(stheta)*(xcos(phi+wdt)+ysin(phi+wdt)) solar frame panel normal: $$ \hat{y}cos(θ_l+θ_p)sin(φ+ω_dt) + \hat{z}\left[sin(θ_l+θ_p)cos(23^ˆ) - sin(23^ˆ)cos(φ+ω_dt)cos(θ_l+θ_p)\right] $$ $$ \hat{x}\left[ cos(23^ˆ)cos(θ_l+θ_p)cos(φ+ω_dt) + sin(23^ˆ)sin(θ_l +θ_p) \right] $$
then for surface normal power we have $$ φ=0\quad θ_p=0 $$ $$ S⋅ n= sin(23^ˆ)sin(θ_l)cos(θ_s) + cos(23^ˆ)cos(θ_l)cos(ω_dt)cos(θ_s)+cos(θ_l)sin(ω_dt)sin(θ_s) $$
we can see all occourences of $$ θ_l $$ in the above expression were originally $$ θ_l + θ_p $$
want to calculate $$ ∫_02π/ω_ddt max(0,a+bsin(ω_dt)) $$ where
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
def panel_power_time(tday,tl,tp,ts):
# tday is going to measured in rads
a=np.sin(np.pi*23.0/180)*np.sin(tl)*np.cos(ts)
b1=np.cos(np.pi*23.0/180)*np.cos(tl)*np.cos(ts)
b2=np.cos(tl)*np.sin(ts)
b=np.sqrt(b1**2+b2**2)
solar_occlude=(a+b*np.sin(tday))<0
tsum=tl+tp
a=np.sin(np.pi*23.0/180)*np.sin(tsum)*np.cos(ts)
b1=np.cos(np.pi*23.0/180)*np.cos(tsum)*np.cos(ts)
b2=np.cos(tsum)*np.sin(ts)
b=np.sqrt(b1**2+b2**2)
v=a + b*np.sin(tday)
v[solar_occlude]=0
v[v<0]=0
return v
def total_power_day(tl,tp,ts):
n=400
xs=np.linspace(0,2*np.pi,n)
v=np.sum(panel_power_time(xs[1:],tl,tp,ts))/float(n-1)
return v*24
def annual_daily_power(tl,tp):
days=np.linspace(0,2*np.pi,100)
vals=[]
for d in days:
vals.append(total_power_day(tl,tp,d))
return vals
def deg(r):
return r*180.0/np.pi
def daily_optimal_angle(tl,ts):
fun=lambda x: -total_power_day(tl,x,ts)
sv=deg(tl)-23*np.cos(ts)
# sv=-45
return minimize(fun,-rad(sv),method='BFGS',options={"eps":1e-5,"gtol":1e-3})
def multi_angle_annual_power(tl,tp_list):
vs=np.asarray([annual_daily_power(tl,tp) for tp in tp_list])
return np.max(vs,axis=0)
def annual_mult_opt(tl,n):
start_angs=-np.linspace(tl-rad(23),tl+rad(23),n)
fun=lambda x: -np.sum(multi_angle_annual_power(tl,x))
return minimize(fun,start_angs,method="BFGS")
def rad(x):
return x/180.0*np.pi
tday=np.linspace(0,6.3,1000)
tl=45.0/180.0*np.pi
# tpl=[0.0,-tl/2,-tl,-3*tl/2]
tpl=[rad(10),0,-rad(20),-rad(45),-rad(70),]
ts=np.pi*0
labels=[]
for v in tpl:
lstr="{} degrees".format(int(-v*180/np.pi))
line,=plt.plot(tday,panel_power_time(tday,tl,v,ts))
labels.append(line.set_label(lstr))
plt.legend()
plt.show()
tl=45/180.0*np.pi
# tp=-67/180.0*np.pi
angs=[22,45,68,0,-5]
days=np.linspace(0,360,100)
# print(annual_daily_power(tl,tp))
labels=[]
for ang in angs:
ls="{} degrees".format(ang)
line, =plt.plot(days,annual_daily_power(tl,-rad(ang)))
labels.append(line.set_label(ls))
plt.title("Daily Output by Panel Angle")
plt.ylabel("Watt Hours")
plt.xlabel("Day in year")
plt.legend()
res=daily_optimal_angle(rad(45),0.0)
print (res)
print (deg(res['x']))
print (total_power_day(rad(45),res['x'],0.0))
reses=[daily_optimal_angle(rad(45),x) for x in np.linspace(0,2*np.pi,100)]
angles=[-deg(v["x"][0]) for v in reses]
optimal_power=[-v["fun"] for v in reses]
line,=plt.plot(np.linspace(0,365,100),angles)
line.set_label("optimal angle")
line,=plt.plot(np.linspace(0,365,100),45-23*np.cos(np.linspace(0,2*np.pi,100)))
line.set_label("45-23*cos(theta)")
plt.legend()
plt.figure(figsize=(12,8))
tl=45/180.0*np.pi
# tp=-67/180.0*np.pi
angs=[22,45,68,0,-5]
days=np.linspace(0,360,100)
# print(annual_daily_power(tl,tp))
labels=[]
for ang in angs:
ls="{} degrees".format(ang)
line, =plt.plot(days,annual_daily_power(tl,-rad(ang)))
labels.append(line.set_label(ls))
line,=plt.plot(days,optimal_power)
line.set_label("daily optimal power")
plt.title("Daily Output by Panel Angle")
plt.ylabel("Watt Hours")
plt.xlabel("Day in year")
plt.legend()
v=( multi_angle_annual_power(rad(45),[-rad(18),-rad(62)]))
plt.plot(v)
plt.plot(annual_daily_power(rad(45),-rad(45)))
plt.ylabel("Daily Solar Hours")
plt.xlabel("Time of Year")
plt.title("Annual Power Output for Best of 18 and 62 Degrees")
v= annual_mult_opt(rad(45),1)
print (v)
print (-deg(v["x"]))
v= annual_mult_opt(rad(45),2)
print (v)
print (-deg(v["x"]))
v= annual_mult_opt(rad(45),3)
print (v)
print (-deg(v["x"]))
plt.figure(figsize=(12,8))
line,=plt.plot(optimal_power)
line.set_label("method theoretical max power")
v=( multi_angle_annual_power(rad(45),[-rad(13),-rad(40),-rad(66)]))
t,=plt.plot(v)
t.set_label("3 positions optimal power")
t,=plt.plot(multi_angle_annual_power(rad(45),[-rad(18),-rad(62)]))
t.set_label("2 positions optimal power")
t,=plt.plot(annual_daily_power(rad(45),-rad(42)))
t.set_label("completely fixed panel")
plt.ylabel("Daily Solar Hours")
plt.xlabel("Time of Year")
plt.title("Annual Power Output for Best of 13, 40, and 66 Degrees")
plt.legend()
v1=multi_angle_annual_power(rad(45),[-rad(42)])
efficiency_1=np.sum(v1)/np.sum(optimal_power)
print("1 angles achieves {:.1f}% efficiency".format(efficiency_1*100))
v3=multi_angle_annual_power(rad(45),[-rad(13),-rad(40),-rad(66)])
efficiency_3=np.sum(v3)/np.sum(optimal_power)
inc3=(efficiency_3/efficiency_1)-1
print("3 angles achieves {:.1f}% efficiency: {:.1f}% more than one position"
.format(efficiency_3*100,inc3*100))
v2=multi_angle_annual_power(rad(45),[-rad(18),-rad(62)])
efficiency_2=np.sum(v2)/np.sum(optimal_power)
inc2=(efficiency_2/efficiency_1)-1
print("2 angles achieves {:.1f}% efficiency: {:.1f}% more than one position".format(efficiency_2*100,inc2*100))
plt.title("")
t,=plt.plot(v3/v1)
# t.set_label("One fixed position vs Daily Max")
t,=plt.plot(v2/v1)
# t.set_label("Two fixed positions vs Daily Max")
# plt.legend()
