#!/usr/bin/env python
# coding: utf-8

# ## Obliczanie optymalnego kąta nachylenia kolektora słonecznego
# 
# Przy projektowaniu i analizie funkcjonowania słonecznych instalacji służących do podgrzewania wody użytkowej oraz instalacji grzewczych zachodzi konieczność obliczenia natężenia promieniowania na powierzchni czołowej kolektora słonecznego zainstalowanego pod określonym kątem $\beta$ oraz w odniesieniu względem osi północ-południe $\gamma$.
# 
# Do obliczeń są niezbędne dane dotyczące natężenia promieniowania słonecznego dla każdego dnia roku. Mogą to być np.  wygenerowane losowo dane z programu SolarSym (http://fluid.itcmp.pwr.wroc.pl/~kasper/solarsym/).

# In[1]:


# wyciszenie komunikatów o błędach
import warnings;warnings.simplefilter('ignore')

# załadowanie modułów
get_ipython().run_line_magic('pylab', 'inline')
get_ipython().run_line_magic('config', "InlineBackend.figure_format = 'retina'")
rcParams['savefig.dpi'] = 120

import math
import operator


# ### Szerokość geograficzna oraz orientacja
# 
# * Szerokość geograficzna $\phi = 51.10$ stopni została przyjęta dla Wrocławia (podobnie jak zestaw danych)
# * Kolektora zorientowany jest idealnie w osi północ-południe $\gamma = 0$ 

# In[2]:


ϕ = math.radians(51.10) # [°]
γ = 0 # [°]


# * Dzień roku (do obliczenia deklinacji)
# * Różnica w minutach pomiędzy godziną 12:00, a godziną w której zarejestrowana maksymalne natężenie. Wartość ujemna oznacza, że wskazówka wychylona jest w lewo, przykładowo: godzina **11:50 = -10 minut**.

# In[3]:


n = 180 # numer dnia
minuty = -10


# ### Deklinacja słońca dla wskazanego dnia
# 
# Deklinacja jest to długość łuku chwilowego koła godzinnego Słońca zawartego pomiędzy płaszczyzną równika niebieskiego a Słońcem. Deklinacja dowolnego punktu sfery niebieskiej przyjmuje wartości od 0° (na równiku niebieskim), do +/-90° ("+" na N od równika i "-" na S od równika). Maksymalna wartość deklinacji Słońca (w momentach przesileń) wynosi obecnie +/-23°26,4'. Deklinację oblicza się wg nastepującego równania:
# 
# $$\delta = 23.45 \cdot \rm{sin}\left( 360 \cdot \left( \frac{284 + n}{365} \right)\right)$$

# In[4]:


δ = math.radians(23.45 * sin(math.radians(360 * (284 + n)/365)))

print ('Deklinacja słońca: %.2f°' % math.degrees(δ))


# ### Kąt godzinowy $\omega$

# In[5]:


ω = math.radians(15 * (minuty / 60.))
print ('Kąt godzinowy: %.2f°' % math.degrees(ω))


# ### Obliczanie kąta nachylenia kolektora $\beta$
# Celem jest ustalenie takiego kąta nachylenia kolektora $\beta$ aby promienie słońca padały na niego w sposób optymalny, tj. mozliwie zbliżony do kierunku normalnego do jego powierzchni.

# #### Definicje funkcji do obliczeń $\rm{cos}\left(\theta_{\beta}\right)$ oraz kąta $\theta_{\beta}$
# 
# Zależność opisująca $\rm{cos}\left(\theta_{\beta}\right)$ w funkcji kątów $\beta$ i $\gamma$ oraz szerokości geograficznej $\phi$, deklinacji słońca $\delta$ i kąta godzinnego $\omega$.
# 
# $$\rm{cos}\theta_{\beta} = \rm{sin}\left(\delta\right) \rm{sin}\left(\phi\right) \rm{cos}\left(\beta\right) - \rm{sin}\left(\delta\right) \rm{cos}\left(\phi\right) \rm{sin}\left(\beta\right) \rm{cos}\left(\gamma\right) + \rm{cos}\left(\delta\right) \rm{cos}\left(\phi\right) \rm{cos}\left(\beta\right) \rm{cos}\left(\omega\right) + \\ + \rm{cos}\left(\delta\right) \rm{sin}\left(\phi\right) \rm{sin}\left(\beta\right) \rm{cos}\left(\gamma\right) + \rm{cos}\left(\delta\right) \rm{sin}\left(\beta\right) \rm{sin}\left(\gamma\right) \rm{sin}\left(\omega\right)$$

# In[6]:


def CTB(aphi, adelta, agamma, abeta, aomega):
    a1 = math.sin(adelta)*math.sin(aphi)*math.cos(abeta) 
    a2 = - math.sin(adelta)*math.cos(aphi)*math.sin(abeta)*math.cos(agamma) 
    a3 = math.cos(adelta)*math.cos(aphi)*math.cos(abeta)*math.cos(aomega) 
    a4 = math.cos(adelta)*math.sin(aphi)*math.sin(abeta)*math.cos(agamma)*math.cos(aomega)
    a5 = math.cos(adelta)*math.sin(abeta)*math.sin(agamma)*math.sin(aomega)
    return a1+a2+a3+a4+a5

def TB(ctb):
    return math.degrees(math.acos(ctb))


# ### Wyznaczenie optymalnego kąta nachylenia kolektora $\beta$ w zakresie 0 do 90 stopni.

# In[7]:


# Definicja tablic
abetaD = np.linspace(0,90,100) # ... beta[] jest tylko po to aby wyświetlić na osi x wartości w stopniach
cos_TB = np.linspace(0,90,100) # ... tablica wyników

abetaR = np.linspace(0,90,100) # ... do obliczeń należy użyć tablicy wartości w radianach

# Obliczenie cos(θ_β)
x = 0
for (i) in abetaR:
    abetaR[x] = math.radians(i) # ... przeliczanie stopni na radiany
    cos_TB[x] = CTB(ϕ,δ,γ,abetaR[x],ω) # obliczanie wartości cos_theta_beta
    x += 1

# Poszukiwanie wartości maksymalnej cos(θ_β)
max_index, max_value = max(enumerate(cos_TB), key=operator.itemgetter(1))

# Optymalna wartość kąta β
β = math.radians(abetaD[max_index])
θ_β = math.acos(cos_TB[max_index])

# Wizualizacja graficzna wyniku obliczeń
plt.figure()
plt.plot(abetaD, cos_TB, 'k-', linewidth=0.75, ms=4, fillstyle='none',label="Prim. chiller")
plt.ylabel(r'$\rm{cos} \left( \theta_{\beta} \right)$',fontsize=9)
plt.xlabel(r'$\beta$',fontsize=9)
plt.annotate(r'$\beta = $' + str(round(abetaD[max_index],2)) + 
             r'$\degree$ $\theta_{\beta} = $' +str(round(math.degrees(math.acos(cos_TB[max_index])),2)), 
             xy=(abetaD[max_index],cos_TB[max_index]), xytext=(abetaD[max_index], cos_TB[max_index]-0.3),
             arrowprops=dict(facecolor='black', shrink=0.05), horizontalalignment='center', 
             verticalalignment='top', fontsize=8)
plt.ylim(0.5,1.1);