% optimizing the vertical tracking angle beta at 25 degC

day_s = [172, 355]; %june and dec 21
sunny_OCI = 0;
lat_pec = 30.2605; %latitude of palmer events center, degrees
temp_panel = 25; %degC, given parameter
gamma_o = 46; %surface azimuth angle, given

t_res = 5; %resolution of time, mins
mins_day = 60*24;
time_of_day = 0:t_res:mins_day;

days_all = 0:1:365;

beta_res = 1; %resolution of inclination, degrees
beta_all = -90:beta_res:90;

%irradiance and power vs time of day, december 21

%allocating memory to store day's values
time_solar = zeros(1,size(time_of_day, 2));
omega_min = time_solar; %hour angle
theta_z = time_solar; %solar zenith
alpha_z = time_solar; %solar altitude
gamma_s = time_solar; %solar azimuth
theta_i = time_solar; %angle of incidence
tau_b = time_solar; %atmospheric transmittance
i_beam = time_solar; %beam irradiance
i_diffuse = time_solar; %diffuse irradiance

irradiance_total = zeros(size(time_of_day,2), size(beta_all,2) );
power_instant = time_solar;

% irradiance_max = zeros(size(day_s,1), size(time_of_day, 2));
irradiance_max = time_solar;
beta_optimal = irradiance_max;
power_max = irradiance_max;
% total_energy_day = irradiance_max;
total_energy_day = zeros(1,365);

for n = 1:365
% n = day_s(2)
    delta = declination(n);
    for k = 1:size(time_of_day, 2)
        for m = 1:size(beta_all,2)-1
            time_solar(k) = solar_time(time_of_day(k), day_s(2)); %calculate solar time
            omega_min(k) = hour_angle(time_of_day(k));
            theta_z(k) = solar_zenith(delta, lat_pec, omega_min(k));
            alpha_z(k) = asind(cosd(theta_z(k)));
            gamma_s(k) = solar_azimuth(delta, omega_min(k), alpha_z(k));
            theta_i(k) = incidence(alpha_z(k), beta_all(m), gamma_o, gamma_s(k));
            tau_b(k) = atmospheric_transmittance(alpha_z(k));
            i_beam(k) = beam_irradiance(n, theta_i(k), tau_b(k));
            i_diffuse(k) = diffuse_irradiance(n, theta_z(k), tau_b(k), beta_all(m));
            irradiance_total(k,m) = i_beam(k) + i_diffuse(k);
            % power_instant(k) = calculate_power(irradiance_total(k), temp_panel); %power produced at 25 degC
        end
    end

    [irradiance_max, beta_optimal] = max(irradiance_total, [], 2);
    power_max = calculate_power(irradiance_max, temp_panel);
    total_energy_day(n) = trapz(time_of_day, power_max) / (60*10000000);
end

total_energy_yearly = sum(total_energy_day);

% figure(1) % n = 2
% yyaxis left
% plot(time_of_day/60, power_max/10000)
% title('Optimized Vertical Tracking, Dec 21 @ 25 degC')
% xlabel('Time of Day (hours)')
% ylabel('Maximum System Power Delivery (kW)')

% yyaxis right
% plot(time_of_day/60, irradiance_max)
% ylabel('Maximum Irradiance (W/m^2)')

% figure(2)
% yyaxis left
% plot(time_of_day/60, power_max/10000)
% title('Optimized Vertical Tracking, Jun 21 @ 25 degC')
% xlabel('Time of Day (hours)')
% ylabel('Maximum System Power Delivery (kW)')

% yyaxis right
% plot(time_of_day/60, irradiance_max)
% ylabel('Maximum Irradiance (W/m^2)')

% figure(3)
% yyaxis left
% plot(time_of_day/60, beta_optimal - 90)
% title('Optimized Panel Angle with Horizontal, Jun 21')
% xlabel('Time of Day (hours)')
% ylabel('Tilt Angle \beta (deg)')

% yyaxis right
% plot(time_of_day/60, power_max/10000)
% ylabel('Maximum System Power Delivery (kW)')

% figure(4)
% plot(1:365, total_energy_day)
% title('Yearly Maximum Energy Production with Vertical Tracking')
% xlabel('Days of the Year')
% ylabel('Total Daily Energy Production (MWh)')
% % ylim([0 1])
% xlim([0 365])