-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathannual2.m
159 lines (130 loc) · 6.03 KB
/
annual2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
% Constants
solar_constant = 1353; % W/m²
Cn = 0.75; % PV system efficiency
latitude = 42.984333; % Syracuse, NY
longitude = -76.142167;
tilt_angle = 30; % Assumed tilt angle (degrees)
surface_azimuth = 10; % Assumed collector azimuth angle (degrees)
albedo = 0.2; % Ground reflectance
household_consumption = 10000; % Yearly household electricity consumption in kWh
conversion_efficiency = 0.15; % 15% solar to electricity conversion
electricity_cost_per_kWh = 0.13; % Cost of electricity per kWh (USD)
% Dates for simulation (equinoxes and solstices)
dates = ["2024-03-21", "2024-06-21", "2024-09-21", "2024-12-21"];
labels = ["March 21", "June 21", "September 21", "December 21"];
% Times from 4 AM to 9 PM (in 1-hour increments)
times = arrayfun(@(x) sprintf('%02d:00', x), 4:21, 'UniformOutput', false);
% Prepare figure
figure;
hold on;
colors = lines(length(dates));
% Variables to store monthly energy and total yearly energy
monthly_energy = zeros(1, length(dates));
total_yearly_energy = 0;
sunrise_times = strings(1, length(dates));
sunset_times = strings(1, length(dates));
for i = 1:length(dates)
date = dates{i};
label = labels{i};
radiation_values = zeros(1, length(times));
non_zero_times = {};
for j = 1:length(times)
time = times{j};
Ic = solar_radiation(latitude, longitude, date, time, tilt_angle, surface_azimuth, albedo, Cn, 0.144, 0.06, solar_constant);
radiation_values(j) = Ic;
if Ic > 0
non_zero_times{end+1} = time;
end
end
% Calculate daily energy in kWh/m²
daily_energy = sum(radiation_values) / 1000; % Convert Wh/m² to kWh/m²
monthly_energy(i) = daily_energy * 30; % Approximate monthly energy (30 days)
% Plot the radiation values
plot(1:length(times), radiation_values, 'DisplayName', label, 'Color', colors(i,:));
% Calculate and store sunrise and sunset times
if ~isempty(non_zero_times)
first_non_zero = datetime(non_zero_times{1}, 'InputFormat', 'HH:mm');
last_non_zero = datetime(non_zero_times{end}, 'InputFormat', 'HH:mm');
sunrise_times(i) = datestr(first_non_zero, 'HH:MM AM');
sunset_times(i) = datestr(last_non_zero, 'HH:MM PM');
else
sunrise_times(i) = 'N/A';
sunset_times(i) = 'N/A';
end
end
% Customize the plot
xlabel('Time (Hour)');
ylabel('Solar Radiation (W/m²)');
title('Total Solar Radiation on the Panel vs. Local Time');
legend show;
grid on;
hold off;
% Sum up total annual energy
total_yearly_energy = sum(monthly_energy);
% Calculate electricity generated by the panel
electricity_generated = total_yearly_energy * conversion_efficiency; % kWh/m²
% Calculate the required PV area
required_PV_area = household_consumption / electricity_generated; % m²
% Calculate the monthly and annual electricity savings cost
monthly_savings_cost = monthly_energy * conversion_efficiency * required_PV_area * electricity_cost_per_kWh; % USD
annual_savings_cost = sum(monthly_savings_cost); % USD
% Display sunrise and sunset times
for i = 1:length(dates)
fprintf('%s:\n', labels{i});
fprintf(' Sunrise: %s\n', sunrise_times(i));
fprintf(' Sunset: %s\n', sunset_times(i));
end
% Display monthly energy, costs, and area
fprintf('\nMonthly Energy (kWh/m²) and Savings Cost (USD):\n');
for i = 1:length(dates)
fprintf('%s: %.2f kWh/m², $%.2f\n', labels{i}, monthly_energy(i), monthly_savings_cost(i));
end
fprintf('\nTotal Annual Energy: %.2f kWh/m²\n', total_yearly_energy);
fprintf('Electricity Generated by the Panel: %.2f kWh/m²\n', electricity_generated);
fprintf('Required PV Area to meet household consumption: %.2f m²\n', required_PV_area);
fprintf('Annual Electricity Savings Cost: $%.2f\n', annual_savings_cost);
% Define the function to calculate solar radiation
function Ic = solar_radiation(lat, lon, date, time, tilt_angle, surface_azimuth, albedo, Cn, k, C, solar_constant)
% Convert inputs to radians
lat = deg2rad(lat);
tilt_angle = deg2rad(tilt_angle);
surface_azimuth = deg2rad(surface_azimuth);
% Parse the date and time
datetime_obj = datetime(date + " " + time, 'InputFormat', 'yyyy-MM-dd HH:mm');
% Calculate day of the year (n)
day_of_year = day(datetime_obj, 'dayofyear');
% Equation of Time (ET)
B = deg2rad(360 * (day_of_year - 81) / 364); % Radians
ET = 9.87 * sin(2 * B) - 7.53 * cos(B) - 1.5 * sin(B); % Minutes
% Convert EST to Solar Time
local_standard_meridian = -75; % For EST (degrees)
solar_time_offset = ET + 4 * (lon - local_standard_meridian);
solar_time = datetime_obj + minutes(solar_time_offset);
solar_hour_angle = deg2rad(15 * (hour(solar_time) + minute(solar_time) / 60 - 12)); % Radians
% Solar Declination (delta)
declination = deg2rad(23.45 * sin(deg2rad(360 / 365 * (day_of_year - 81))));
% Solar Altitude Angle (alpha)
sin_alpha = sin(lat) * sin(declination) + cos(lat) * cos(declination) * cos(solar_hour_angle);
alpha = asin(sin_alpha);
% Direct Normal Irradiance (Ib,n)
if rad2deg(alpha) > 0
Ib_n = Cn * solar_constant * (1 + 0.034 * cos(deg2rad(360 * day_of_year / 365.25))) * exp(-k / sin(alpha));
else
Ib_n = 0;
end
% Diffuse Radiation (Id)
Id = C * Ib_n * cos(tilt_angle / 2) * cos(tilt_angle / 2);
% Ground Reflected Radiation (Ir)
Ir = albedo * Ib_n * sin(alpha + C) * sin(tilt_angle / 2) * sin(tilt_angle / 2);
% Angle of Incidence (theta)
cos_theta = (sin(declination) * sin(lat) * cos(tilt_angle)) - ...
(sin(declination) * cos(lat) * sin(tilt_angle) * cos(surface_azimuth)) + ...
(cos(declination) * cos(lat) * cos(tilt_angle) * cos(solar_hour_angle)) + ...
(cos(declination) * sin(lat) * sin(tilt_angle) * cos(surface_azimuth) * cos(solar_hour_angle)) + ...
(cos(declination) * sin(tilt_angle) * sin(surface_azimuth) * sin(solar_hour_angle));
cos_theta = max(cos_theta, 0);
% Beam Radiation on the Collector (Ib_n * cos(theta))
Ib_cos_theta = Ib_n * cos_theta;
% Total Radiation on the Collector (Ic)
Ic = Ib_cos_theta + Id + Ir;
end