Моделирование солнечного трекера (ч.1)
Author
In [ ]:
using Plots
using Printf
plotlyjs()
# ============================================================
# ПАРАМЕТРЫ МЕСТОПОЛОЖЕНИЯ
# Координаты объекта установки солнечных панелей
# Можно изменить на любые координаты в градусах/минутах/секундах
# ============================================================
phi_deg = 48; phi_min = 1; phi_sec = 54.42 # широта: 48°01'54.42"N
lam_deg = 37; lam_min = 47; lam_sec = 14.90 # долгота: 37°47'14.90"E
# Перевод в десятичные градусы (стандартный формат для расчётов)
phi = phi_deg + phi_min/60 + phi_sec/3600 # географическая широта, °
lambda = lam_deg + lam_min/60 + lam_sec/3600 # географическая долгота, °
timezone = 3 # часовой пояс UTC+N (Москва/Донецк = UTC+3)
E0 = 1360.8 # солнечная постоянная, Вт/м² — поток на границе атмосферы
A = 0.2 # альбедо подстилающей поверхности:
# 0.2 — типичный грунт/трава
# 0.5 — бетон/асфальт
# 0.8 — снежный покров
dO = 0.32 # толщина озонового слоя, см (типовое значение для средних широт)
alpha_aer = 0.0314 # коэффициент аэрозольного помутнения атмосферы
# (Ангстрем): 0.05 — чистый воздух, 0.1..0.2 — городской
# ============================================================
# ГОД РАСЧЁТА
# Меняйте year для расчёта любого года
# Код автоматически учитывает високосный год (366 дней)
# ============================================================
year = 2026
# Функция проверки високосного года по григорианскому календарю:
# делится на 4 И (не делится на 100 ИЛИ делится на 400)
is_leap(y) = (y % 4 == 0 && (y % 100 != 0 || y % 400 == 0))
# Количество дней в каждом месяце с учётом високосного года
function get_days_in_month(y)
d = [31,28,31,30,31,30,31,31,30,31,30,31]
is_leap(y) && (d[2] = 29) # февраль: 28 или 29 дней
return d
end
dim = get_days_in_month(year) # массив длин месяцев
days_in_year = is_leap(year) ? 366.0 : 365.0 # длина года в днях
@printf "Местоположение: %.4f°N %.4f°E\n" phi lambda
@printf "Год: %d %s\n\n" year (is_leap(year) ? "(високосный)" : "")
# ============================================================
# СТАТИСТИЧЕСКИЕ ДАННЫЕ ДЛЯ ВЕРИФИКАЦИИ МОДЕЛИ
#
# Пиковые солнечные часы (ПСЧ) для ГОРИЗОНТАЛЬНОЙ поверхности
# Единица: кВт·ч/м²/день = количество часов с мощностью 1 кВт/м²
# Источник: реальные многолетние измерения для данного региона
#
# ВАЖНО: эти данные используются для калибровки модели —
# модель считает "ясное небо", статистика учитывает облачность.
# Коэффициент k_cloud = статистика / модель автоматически
# приводит расчёт к реальным условиям.
#
# Для другого региона замените значения на актуальные данные
# (источники: NASA POWER, PVGIS, Solargis, местные метеостанции)
# ============================================================
stat_monthly = [1.21, 1.99, 2.94, 4.04, 5.48, 5.55,
5.66, 5.09, 3.67, 2.24, 1.23, 0.96]
# Янв Фев Мар Апр Май Июн
# Июл Авг Сен Окт Ноя Дек
stat_annual = 3.34 # среднегодовое значение, кВт·ч/м²/день
# ============================================================
# ТАБЛИЦА ОПТИМАЛЬНЫХ УГЛОВ НАКЛОНА ФИКСИРОВАННОЙ ПАНЕЛИ
#
# Получена методом полного перебора углов 0..90° для каждого
# ключевого дня года. Для каждого значения солнечного
# склонения δ определён угол γ, максимизирующий суточную
# выработку фиксированной панели, ориентированной на юг.
#
# Физический смысл: зимой (δ < 0) Солнце низко → панель
# ставится круче (большой угол). Летом (δ > 0) Солнце высоко
# → панель почти горизонтальна (малый угол).
#
# Для другой широты таблицу нужно пересчитать (запустить
# блок "Проверка оптимальных углов по ключевым дням")
#
# Хранится как именованный кортеж — без глобальных констант,
# без предупреждений при повторном запуске скрипта
# ============================================================
OPT_LOOKUP = (
delta = [-23.4, -23.0, -17.3, -7.8, 0.0, 8.5, 14.9, 18.2, 22.0, 23.4],
# зима лето
angle = [ 74.0, 74.0, 69.0, 59.0, 49.0, 37.0, 26.0, 21.0, 14.0, 11.0]
# круто полого
)
# Линейная интерполяция между узлами таблицы
# Позволяет получить оптимальный угол для любого значения δ
function optimal_angle_from_table(delta_deg, tbl)
# Граничные значения: за пределами таблицы — крайнее значение
delta_deg <= tbl.delta[1] && return tbl.angle[1]
delta_deg >= tbl.delta[end] && return tbl.angle[end]
for i in 1:length(tbl.delta)-1
if tbl.delta[i] <= delta_deg <= tbl.delta[i+1]
# Линейная интерполяция: t = 0 в левом узле, t = 1 в правом
t = (delta_deg - tbl.delta[i]) / (tbl.delta[i+1] - tbl.delta[i])
return tbl.angle[i] + t*(tbl.angle[i+1] - tbl.angle[i])
end
end
return 45.0 # fallback — не должен достигаться
end
# ============================================================
# ФУНКЦИЯ: положение Солнца + коэффициенты пропускания атмосферы
#
# Реализует математическую модель:
# Шаг 4 — реальное местное время (РМВ)
# Шаг 5 — часовой угол ω
# Шаг 6 — высота γ_C и азимут α_C Солнца
# Атмосфера — модель Берда–Атвотера (Bird & Atwater, 1985)
#
# Входные параметры:
# MV — местное время, ч (0..24)
# phi — широта, °
# lambda — долгота, °
# timezone — часовой пояс UTC+N
# delta — солнечное склонение, рад
# TE — уравнение времени, мин
# E0 — солнечная постоянная, Вт/м²
# dO — толщина озонового слоя, см
# alpha_aer — коэффициент аэрозольного помутнения
# day — номер дня в году (1..365/366)
# days_in_year — длина года в днях
#
# Возвращает:
# gamma_C — высота Солнца над горизонтом, °
# alpha_C — азимут Солнца от севера (0=С, 90=В, 180=Ю, 270=З), °
# Edir_hor — прямое излучение на горизонталь, Вт/м²
# Ediff_hor — диффузное излучение на горизонталь, Вт/м²
# EG_hor — суммарное излучение на горизонталь, Вт/м²
# cos_th — косинус зенитного угла Солнца
# ============================================================
function sun_and_atmo(MV, phi, lambda, timezone, delta, TE,
E0, dO, alpha_aer, day, days_in_year)
d2r = pi/180 # коэффициент перевода градусов в радианы
# --- Шаг 4. Реальное местное время ---
# СМВ — среднее местное время: учитывает положение внутри часового пояса
# (4 мин на каждый градус долготы от центрального меридиана пояса)
SMV = MV - timezone + 4.0*lambda/60.0 # среднее местное время, ч
# РМВ — реальное (истинное солнечное) время: добавляем уравнение времени
RMV = SMV + TE/60.0 # реальное местное время, ч
# --- Шаг 5. Часовой угол ω ---
# Угловое смещение Солнца от меридиана наблюдателя
# 15°/ч = 360°/24ч — угловая скорость вращения Земли
# ω > 0 утром (Солнце восточнее), ω = 0 в полдень, ω < 0 вечером
omega = (12.0 - RMV)*15.0 # градусы
# --- Шаг 6а. Высота Солнца γ_C ---
# Угол между направлением на Солнце и горизонтальной плоскостью
# 0° — на горизонте (восход/закат), максимум — в полдень
sin_gC = clamp(cos(omega*d2r)*cos(delta)*cos(phi*d2r) +
sin(delta)*sin(phi*d2r), -1.0, 1.0)
gamma_C = asin(sin_gC)/d2r # градусы
# Ночное время — Солнце под горизонтом, всё излучение = 0
gamma_C <= 0.0 && return gamma_C, 0.0, 0.0, 0.0, 0.0, 0.0
# --- Шаг 6б. Азимут Солнца α_C ---
# Используем atan2 — непрерывная функция без разрыва в полдень
cos_gC = max(cos(gamma_C*d2r), 1e-9)
sin_aC = cos(delta)*sin(omega*d2r)/cos_gC
cos_aC = (sin(gamma_C*d2r)*sin(phi*d2r) - sin(delta)) /
max(cos_gC*cos(phi*d2r), 1e-9)
alpha_C = mod(180.0 + atan(sin_aC, cos_aC)/d2r, 360.0) # 0..360°
# --- Воздушная масса Mi (формула Кастена) ---
# Относительная длина пути лучей через атмосферу
# Mi = 1 при Солнце в зените, резко растёт у горизонта
# Учитывает кривизну атмосферы при малых высотах Солнца
theta_hor = 90.0 - gamma_C # зенитный угол Солнца, °
cos_th = max(cos(theta_hor*d2r), 1e-6)
Mi = 1.0/(cos_th + 0.15*(93.885 - theta_hor)^(-1.253))
# --- Коэффициенты пропускания атмосферы (модель Берда–Атвотера) ---
# Каждый τ описывает отдельный физический механизм ослабления излучения.
# Суммарное ослабление = произведение всех τ.
# τ_R — рэлеевское рассеяние молекулами воздуха (N₂, O₂, Ar)
# Сильнее для коротких волн (синий свет), отсюда голубое небо
tR = exp(-0.0903*Mi^0.84*(1.0+Mi-Mi^1.01))
# τ_O — поглощение ультрафиолета озоновым слоем (O₃)
# dO — толщина слоя в сантиметрах при нормальных условиях
tO = 1.0 - 0.1611*dO*Mi*(1.0+139.48*dO*Mi)^(-0.3035) -
0.002715*dO*Mi*(1.0+0.044*dO*Mi+0.0003*(dO*Mi)^2)^(-1)
# τ_G — поглощение смесью газов атмосферы (CO₂, O₂)
# Влияет в основном на инфракрасную часть спектра
tG = exp(-0.0127*Mi^0.26)
# Сезонная поправка на влажность атмосферы:
# season = 0 зимой (δ = −23.5°), season = 1 летом (δ = +23.5°)
# Летом воздух теплее и суше → меньше водяного пара и слабее поглощение
delta_deg = delta*180/pi
season = clamp((delta_deg + 23.5)/47.0, 0.0, 1.0)
dH2O = 2.0 - 0.7 *season # содержание водяного пара: 2.0 зимой → 1.3 летом, см
Khor = 1.14 - 0.23*season # сезонная поправка потока: 1.14 зимой → 0.91 летом
# τ_H2O — поглощение водяными парами (инфракрасный диапазон)
tH2O = 1.0 - 2.4959*dH2O*Mi*
((1.0+79.034*dH2O*Mi)^0.6828 + 6.385*dH2O*Mi)^(-1)
# τ_A — аэрозольное рассеяние и поглощение (пыль, сажа, туман, дымка)
# TA — оптическая толщина аэрозоля (зависит от alpha_aer)
TA = 1.832*alpha_aer
tA = exp(-TA^0.873*(1.0+TA-TA^0.7088)*Mi^0.9108)
# tAA — пропускание без обратного рассеяния
tAA = 1.0 - 0.1*(1.0-Mi+Mi^1.06)*(1.0-tA)
# tAS — пропускание только за счёт прямого рассеяния вперёд
tAS = tA/tAA
# --- Прямое излучение на горизонтальную поверхность ---
# E0 × cos(зенит) × произведение всех τ × сезонная поправка
# Чем ниже Солнце → больше Mi → меньше каждый τ → меньше поток
Edir_hor = E0*cos_th*tR*tO*tG*tH2O*tA*Khor
# --- Диффузное излучение на горизонтальную поверхность ---
# Свет, рассеянный атмосферой и приходящий со всего небосвода
# Kdiff = 1.05 — эмпирический коэффициент модели Берда–Атвотера
Ediff_hor = E0*cos_th*tO*tG*tH2O*tAA*(1.0-tR*tAS)/(1.0-Mi+Mi^1.02)*1.05
# --- Суммарное (глобальное) излучение на горизонталь ---
EG_hor = Edir_hor + Ediff_hor
return gamma_C, alpha_C, Edir_hor, Ediff_hor, EG_hor, cos_th
end
# ============================================================
# ФУНКЦИЯ: облучённость наклонной панели
#
# Пересчитывает горизонтальное излучение на наклонную панель
# с произвольной ориентацией (модель Клюхера для диффузного).
#
# Входные параметры:
# gamma_C, alpha_C — положение Солнца (высота и азимут), °
# Edir_h, Ediff_h, EG_h — составляющие на горизонталь, Вт/м²
# cos_th — косинус зенитного угла Солнца
# gamma_E — угол наклона панели к горизонту, °
# 0° = горизонтальная, 90° = вертикальная
# alpha_E — азимут панели от севера, °
# 0° = север, 90° = восток, 180° = юг, 270° = запад
# A — альбедо подстилающей поверхности (доля отражения)
#
# Возвращает суммарную облучённость панели, Вт/м²
# ============================================================
function panel_irr(gamma_C, alpha_C, Edir_h, Ediff_h, EG_h, cos_th,
gamma_E, alpha_E, A)
d2r = pi/180
gamma_C <= 0.0 && return 0.0 # ночь — нет излучения
# --- Угол падения прямых лучей на панель ---
# cos_inc = cos угла между нормалью к панели и направлением на Солнце
# = 1 при перпендикулярном падении (максимум мощности)
# = 0 при параллельном падении (нет мощности)
cos_inc = clamp(sin(gamma_C*d2r)*cos(gamma_E*d2r) +
cos(gamma_C*d2r)*sin(gamma_E*d2r)*cos((alpha_C-alpha_E)*d2r),
0.0, 1.0)
# --- Прямое излучение на панель ---
# Пересчёт через соотношение косинусов углов падения:
# E_пан = E_гор × cos_inc / cos(θ_гор)
Edir_gen = Edir_h*cos_inc/cos_th
# --- Диффузное излучение на панель (модель Клюхера) ---
# F — коэффициент анизотропии (ясности) неба:
# F → 1 при ясном небе: небо ярче у Солнца и у горизонта
# F → 0 при облачности: небо изотропно (равномерное)
# (1+cos γ_E)/2 — геометрический фактор: доля видимого небосвода
# (1 + F·sin³(γ_E/2)) — усиление от яркого горизонтального пояса
F = 1.0 - (Ediff_h/max(EG_h,1e-6))^2
Ediff_gen = Ediff_h*0.5*(1.0+cos(gamma_E*d2r))*
(1.0 + F*sin(gamma_E*d2r/2.0)^3)
# --- Отражённое излучение от земли на панель ---
# Пропорционально альбедо A и доле видимой земли (1-cos γ_E)/2
# При γ_E = 0° (горизонталь) земли не видно → E_ref = 0
# При γ_E = 90° (вертикаль) земля видна максимально → E_ref = max
Eref_gen = EG_h*A*0.5*(1.0-cos(gamma_E*d2r))
# --- Суммарная облучённость панели ---
return Edir_gen + Ediff_gen + Eref_gen
end
# ============================================================
# ФУНКЦИЯ: суточная выработка четырёх режимов для одного дня
#
# Интегрирует мощность по времени методом трапеций (шаг 1 мин).
# Возвращает суточную энергию в Вт·ч/м²/день для каждого режима.
#
# Четыре режима:
# W_h — горизонтальная панель (γ=0°)
# Используется для калибровки модели по статистике
#
# W_f — фиксированная панель с оптимальным углом из таблицы
# Угол меняется плавно в течение года (≈74° зимой, ≈11° летом)
# Соответствует панели с ручной сезонной регулировкой
# Ориентация: строго на юг (α = 180°)
#
# W_a — одноосевой трекер по азимуту
# Наклон γ = φ (фиксирован весь год, равен широте места)
# Ось вращения горизонтальная, параллельна земной оси
# Азимут следит за Солнцем в течение дня
#
# W_t — двухосевой трекер
# Нормаль панели всегда направлена точно на Солнце
# γ_E = 90° − γ_C (зенитный угол), α_E = α_C (азимут Солнца)
# Теоретический максимум — верхняя граница выработки
# ============================================================
function calc_day(day, days_in_year, phi, lambda, timezone,
E0, dO, alpha_aer, A, opt_lookup)
dt = 1.0/60.0 # шаг интегрирования = 1 минута = 1/60 часа
MV_range = 0.0:dt:24.0
# Параметр дня J (градусы) — для формул склонения и уравнения времени
J = 360.0*day/days_in_year
# Солнечное склонение δ (рад) — формула DIN 5034
# Три слагаемых: наклон оси Земли + эллиптичность орбиты + возмущения
delta = (0.3948 - 23.2559*cos((J+9.1)*pi/180)
- 0.3915*cos((2*J+5.4)*pi/180)
- 0.1764*cos((3*J+26.0)*pi/180))*pi/180
# Уравнение времени TE (мин) — поправка на эллиптичность орбиты
# Меняется от −16 до +14 минут в течение года
TE = 0.0066 + 7.3525*cos((J+85.9)*pi/180) +
9.9359*cos((2*J+108.9)*pi/180) +
0.3387*cos((3*J+105.2)*pi/180)
delta_deg = delta * 180/pi
# Оптимальный угол фикс. панели для данного дня (из таблицы)
gamma_fixed = optimal_angle_from_table(delta_deg, opt_lookup)
# Угол азимутального трекера = широта места (постоянный весь год)
# При γ_E = φ ось вращения горизонтальна и параллельна земной оси
gamma_az = phi
# Накопители энергии и предыдущие значения мощности (для трапеций)
W_f=0.0; W_a=0.0; W_t=0.0; W_h=0.0
pf=0.0; pa=0.0; pt=0.0; ph=0.0
for (k,MV) in enumerate(MV_range)
gC,aC,Eh,Dh,EGh,cth =
sun_and_atmo(MV,phi,lambda,timezone,delta,TE,
E0,dO,alpha_aer,day,days_in_year)
cf = panel_irr(gC,aC,Eh,Dh,EGh,cth, gamma_fixed, 180.0, A)
ca = panel_irr(gC,aC,Eh,Dh,EGh,cth, gamma_az, aC, A)
ct = panel_irr(gC,aC,Eh,Dh,EGh,cth, 90.0-gC, aC, A)
ch = panel_irr(gC,aC,Eh,Dh,EGh,cth, 0.0, 180.0, A)
# Метод трапеций: W += (E_пред + E_тек)/2 × Δt
if k > 1
W_f += (pf+cf)/2*dt; W_a += (pa+ca)/2*dt
W_t += (pt+ct)/2*dt; W_h += (ph+ch)/2*dt
end
pf=cf; pa=ca; pt=ct; ph=ch
end
return W_f, W_a, W_t, W_h, gamma_fixed
end
# ============================================================
# ГЛАВНЫЙ ЦИКЛ: перебор всех дней года
# ============================================================
month_names = ["Янв","Фев","Мар","Апр","Май","Июн",
"Июл","Авг","Сен","Окт","Ноя","Дек"]
# Массивы суточных результатов (Вт·ч/м²/день)
W_f_daily=Float64[]; W_a_daily=Float64[]
W_t_daily=Float64[]; W_h_daily=Float64[]
ang_daily=Float64[] # оптимальный угол каждого дня, °
days_axis=Int[] # номер дня в году (1..365)
months_arr=Int[] # номер месяца каждого дня (1..12)
println("Расчёт суточной выработки за $year год...")
day_num = 0
for mo in 1:12
for d in 1:dim[mo]
day_num += 1
W_f,W_a,W_t,W_h,gf = calc_day(day_num, days_in_year, phi, lambda,
timezone, E0, dO, alpha_aer, A,
OPT_LOOKUP)
push!(W_f_daily,W_f); push!(W_a_daily,W_a)
push!(W_t_daily,W_t); push!(W_h_daily,W_h)
push!(ang_daily,gf); push!(days_axis,day_num)
push!(months_arr,mo)
end
@printf " %s: готово\n" month_names[mo]
end
# ============================================================
# АГРЕГАЦИЯ РЕЗУЛЬТАТОВ
# ============================================================
# Суммы по месяцам (Вт·ч/м²/мес)
W_f_month = [sum(W_f_daily[months_arr .== m]) for m in 1:12]
W_a_month = [sum(W_a_daily[months_arr .== m]) for m in 1:12]
W_t_month = [sum(W_t_daily[months_arr .== m]) for m in 1:12]
W_h_month = [sum(W_h_daily[months_arr .== m]) for m in 1:12]
# Годовые суммы (кВт·ч/м²/год)
W_f_year = sum(W_f_daily) / 1000.0
W_a_year = sum(W_a_daily) / 1000.0
W_t_year = sum(W_t_daily) / 1000.0
W_h_year = sum(W_h_daily) / 1000.0
# Среднесуточное горизонталь (кВт·ч/м²/день) — для расчёта k_cloud
W_h_daily_avg = W_h_year / days_in_year
# Среднесуточные по каждому месяцу (кВт·ч/м²/день)
W_h_avg_day = (W_h_month ./ Float64.(dim)) ./ 1000.0
# ============================================================
# КОЭФФИЦИЕНТ ОБЛАЧНОСТИ k_cloud
#
# Модель Берда–Атвотера рассчитывает облучённость для условий
# ЯСНОГО НЕБА — без облаков, тумана и осадков.
# Реальная выработка всегда меньше из-за этих факторов.
#
# k_cloud = статистика_горизонталь / модель_горизонталь
#
# Физический смысл:
# k_cloud = 0.58 означает что в реальности панель получает
# 58% от теоретического максимума ясного неба
# 42% теряется из-за облачности, тумана, осадков, дымки
#
# Важно: k_cloud применяется ОДИНАКОВО ко всем режимам,
# поэтому ОТНОСИТЕЛЬНЫЙ выигрыш трекеров от облачности не меняется.
# Меняются только абсолютные значения выработки.
#
# Для уточнения k_cloud используйте месячные данные из stat_monthly —
# тогда каждый месяц будет иметь свой коэффициент облачности.
# ============================================================
k_cloud = stat_annual / W_h_daily_avg
# Скорректированные значения с учётом реальной облачности
W_f_year_k = W_f_year * k_cloud # кВт·ч/м²/год, фикс. панель
W_a_year_k = W_a_year * k_cloud # кВт·ч/м²/год, аз. трекер
W_t_year_k = W_t_year * k_cloud # кВт·ч/м²/год, 2ос. трекер
W_h_year_k = W_h_year * k_cloud # кВт·ч/м²/год, горизонталь
W_h_avg_day_k = W_h_avg_day .* k_cloud # кВт·ч/м²/день по месяцам
W_f_month_k = W_f_month .* k_cloud # Вт·ч/м²/мес, фикс. панель
W_a_month_k = W_a_month .* k_cloud # Вт·ч/м²/мес, аз. трекер
W_t_month_k = W_t_month .* k_cloud # Вт·ч/м²/мес, 2ос. трекер
W_h_month_k = W_h_month .* k_cloud # Вт·ч/м²/мес, горизонталь
# ============================================================
# ВЫВОД ИТОГОВОЙ ТАБЛИЦЫ В КОНСОЛЬ
# ============================================================
println()
println("="^70)
@printf " Год %d | φ=%.4f° | k_cloud=%.4f (потери от облачности %.0f%%)\n" year phi k_cloud (1-k_cloud)*100
println("="^70)
@printf " Горизонталь — ясное небо : %7.1f кВт·ч/м²/год (%.2f кВт·ч/м²/день)\n" W_h_year W_h_daily_avg
@printf " Горизонталь — с облачностью : %7.1f кВт·ч/м²/год (%.2f кВт·ч/м²/день)\n" W_h_year_k stat_annual
@printf " Фикс. опт. угол — с облачностью : %7.1f кВт·ч/м²/год\n" W_f_year_k
@printf " Трекер по азимуту — с облачностью : %7.1f кВт·ч/м²/год\n" W_a_year_k
@printf " Двухосевой — с облачностью : %7.1f кВт·ч/м²/год\n" W_t_year_k
println("-"^70)
# Приросты считаются из модели ясного неба (k не влияет на %)
@printf " Прирост аз. vs фикс. : %+.1f %% (%+.0f кВт·ч/м²/год)\n" (W_a_year/W_f_year-1)*100 (W_a_year_k-W_f_year_k)
@printf " Прирост 2ос. vs фикс. : %+.1f %% (%+.0f кВт·ч/м²/год)\n" (W_t_year/W_f_year-1)*100 (W_t_year_k-W_f_year_k)
@printf " Прирост 2ос. vs аз. : %+.1f %% (%+.0f кВт·ч/м²/год)\n" (W_t_year/W_a_year-1)*100 (W_t_year_k-W_a_year_k)
println("="^70)
# Детальная таблица сравнения с статистикой по месяцам
println()
println("="^62)
println(" Сравнение расчёта с реальными измерениями")
println(" (горизонтальная поверхность, кВт·ч/м²/день)")
println("="^62)
@printf " %-8s %10s %10s %10s\n" "Месяц" "Расч.+обл." "Статист." "Откл.%"
println("="^62)
for m in 1:12
cv = W_h_avg_day_k[m] # расчёт с поправкой на облачность
sv = stat_monthly[m] # реальные измерения
dev = (cv/sv - 1)*100 # отклонение в %
@printf " %-8s %8.2f %8.2f %+6.1f %%\n" month_names[m] cv sv dev
end
println("-"^62)
@printf " %-8s %8.2f %8.2f %+6.1f %%\n" "Среднее" (W_h_year_k/days_in_year) stat_annual 0.0
println("="^62)
println()
# ============================================================
# ГРАФИКИ
# ============================================================
loc_str = @sprintf "%.4f°N %.4f°E (%d г.)" phi lambda year
month_starts = cumsum([1; dim[1:11]]) # первые дни каждого месяца
# --- График 1: модель горизонталь vs статистика по месяцам ---
# Показывает качество калибровки модели:
# светлые столбцы — модель ясного неба (завышена)
# тёмные столбцы — модель × k_cloud (должна совпасть с красными точками)
# красные точки — реальные измерения (статистика)
p1 = plot(
title = "Горизонталь: модель vs статистика\n$(loc_str)",
xlabel = "Месяц",
ylabel = "кВт·ч/м²/день",
grid=true, legend=:outertop,
xticks=(1:12, month_names), size=(1000,520))
bar!(p1, (1:12).-0.2, W_h_avg_day,
bar_width=0.35,
label="Модель ясное небо (ср.год=$(round(W_h_daily_avg,digits=2)))",
color=:steelblue, alpha=0.4)
bar!(p1, (1:12).+0.2, W_h_avg_day_k,
bar_width=0.35,
label="Модель × k=$(round(k_cloud,digits=3)) (ср.год=$(round(stat_annual,digits=2)))",
color=:steelblue, alpha=0.9)
plot!(p1, 1:12, stat_monthly,
label="Статистика (ср.год=$(stat_annual))",
lw=2, color=:red, marker=:circle, markersize=7)
hline!(p1, [stat_annual], color=:red, lw=1, ls=:dash, label="")
hline!(p1, [W_h_daily_avg], color=:steelblue, lw=1, ls=:dash, label="")
display(p1)
# --- График 2: суточная выработка трёх режимов за год ---
# С поправкой на облачность (умножение на k_cloud)
# Показывает сезонный ход выработки и разницу между режимами
p2 = plot(
title = "Суточная выработка (k=$(round(k_cloud,digits=3)))\n$(loc_str)",
xlabel = "День года",
ylabel = "Энергия, Вт·ч/м²/день",
grid=true, legend=:outertop, size=(1000,520))
plot!(p2, days_axis, W_f_daily.*k_cloud,
label="Фикс. опт. → $(round(W_f_year_k,digits=1)) кВт·ч/м²/год",
lw=2, color=:blue)
plot!(p2, days_axis, W_a_daily.*k_cloud,
label="Трекер по азимуту (γ=$(round(phi,digits=1))°) → $(round(W_a_year_k,digits=1)) кВт·ч/м²/год",
lw=2, color=:orange)
plot!(p2, days_axis, W_t_daily.*k_cloud,
label="Двухосевой → $(round(W_t_year_k,digits=1)) кВт·ч/м²/год",
lw=2, color=:red)
vline!(p2, month_starts, color=:lightgray, lw=1, label="")
for (i,ms) in enumerate(month_starts)
annotate!(p2, ms+3, maximum(W_t_daily.*k_cloud)*0.97,
text(month_names[i],7,:gray))
end
display(p2)
# --- График 3: помесячная выработка — столбчатая диаграмма ---
# Удобно для сравнения сезонного распределения энергии
# Позволяет оценить выигрыш каждого режима по месяцам
bw = 0.2 # ширина одного столбца
p3 = plot(
title = "Помесячная выработка (с поправкой на облачность)\n$(loc_str)",
xlabel = "Месяц",
ylabel = "кВт·ч/м²/мес",
grid=true, legend=:outertop,
xticks=(1:12, month_names), size=(1000,500))
bar!(p3, (1:12).-1.5bw, W_h_month_k./1000, bar_width=bw,
label="Горизонталь (γ=0°)", color=:gray, alpha=0.7)
bar!(p3, (1:12).-0.5bw, W_f_month_k./1000, bar_width=bw,
label="Фикс. опт. угол", color=:blue, alpha=0.8)
bar!(p3, (1:12).+0.5bw, W_a_month_k./1000, bar_width=bw,
label="Трекер по азимуту", color=:orange, alpha=0.8)
bar!(p3, (1:12).+1.5bw, W_t_month_k./1000, bar_width=bw,
label="Двухосевой трекер", color=:red, alpha=0.8)
display(p3)
# --- График 4: накопленная энергия нарастающим итогом ---
# Показывает суммарную выработку от начала года до каждого дня
# Наклон кривой = текущая суточная мощность
# Разрыв между кривыми = накопленный выигрыш от трекера
p4 = plot(
title = "Накопленная энергия за год (с поправкой)\n$(loc_str)",
xlabel = "День года",
ylabel = "кВт·ч/м²",
grid=true, legend=:outertop, size=(1000,500))
plot!(p4, days_axis, cumsum(W_f_daily.*k_cloud)./1000,
label="Фикс. → $(round(W_f_year_k,digits=1)) кВт·ч/м²/год",
lw=2, color=:blue)
plot!(p4, days_axis, cumsum(W_a_daily.*k_cloud)./1000,
label="Трекер по азимуту → $(round(W_a_year_k,digits=1)) кВт·ч/м²/год",
lw=2, color=:orange)
plot!(p4, days_axis, cumsum(W_t_daily.*k_cloud)./1000,
label="Двухосевой → $(round(W_t_year_k,digits=1)) кВт·ч/м²/год",
lw=2, color=:red)
vline!(p4, month_starts, color=:lightgray, lw=1, label="")
for (i,ms) in enumerate(month_starts)
annotate!(p4, ms+3, W_t_year_k*0.97, text(month_names[i],7,:gray))
end
display(p4)