Reason اعتبارك خبير في البرمجة بلغة الماتلاب ولديك خبرة كبيرة جدا في تصميم ومحاكاة النماذج الرياضية لتقدير الاشعاعات الشمسية المختلفة المباشرة والمتشرة والكلية وغيرها المستلمة على الأسطح في أماكن محددة المطلوب منك تحليل الكود التالي وفهمه بعمق مع التركيز أكثر على تحليل الصيغ الرياضية المقدمة وتحليلها بالتفصيل وفهمها بعمق لتقوم بتقديم تشريح مفصل له   مع التدقيق جيدا والتأكد من تطبيق الصيغ الرياضية للنموذج  التالي حرفيا بدون زيادة أونقصان أو تعديل اذا وجدت أي اختلاف فيي الصيغ قدمها لي     % كود MATLAB لحساب الإشعاع الشمسي على الأسطح الأفقية والمائلة في ورقلة، الجزائر
% باستخدام حسابات فلكية دقيقة لأوقات الشروق والغروب ومدى النهار

% مسح جميع المتغيرات وإغلاق الأشكال
clear all;
close all;
clc;

% ------------------------- تعريف المتغيرات والثوابت -------------------------

% الثابت الشمسي (I_sc): قيمة ثابتة تمثل الإشعاع الشمسي خارج الغلاف الجوي (W/m^2)
I_sc = 1367;

% موقع ورقلة، الجزائر
phi_deg = 31.95; % خط العرض بالدرجات شمالًا
phi_rad = deg2rad(phi_deg);

L_deg = 5.33; % خط الطول بالدرجات شرقًا

% الارتفاع عن سطح البحر (Z): الارتفاع بالأمتار
Z = 140; % ورقلة، الجزائر

% المنطقة الزمنية (Timezone): للجزائر UTC+1
Timezone = 1;

% التاريخ المطلوب (year, month, day)
year = 2025;
month = 1;
day = 4; % مثال: 3 يناير 2025

% زاوية ميل السطح المائل (beta_deg): زاوية الميل بالدرجات (°)
beta_deg_inclined = 30;
beta_rad_inclined = deg2rad(beta_deg_inclined);

% معامل الانعكاس السطحي (rho): قيمة بين 0 و 1 تمثل نسبة الإشعاع المنعكس من السطح
rho = 0.2;

% ------------------------- حساب رقم اليوم في السنة (N) -------------------------

% حساب رقم اليوم في السنة (N)
N = dayOfYear(year, month, day);

% ------------------------- حساب الميل الشمسي (delta) -------------------------

delta_deg = 23.45 * sin(deg2rad(360 * (284 + N) / 365));
delta_rad = deg2rad(delta_deg);

% ------------------------- حساب معادلة الزمن (E) -------------------------

B = deg2rad(360 * (N - 81) / 365);
E = 9.87 * sin(2 * B) - 7.53 * cos(B) - 1.5 * sin(B); % بالدقائق

% ------------------------- حساب أوقات الشروق والغروب بدقة عالية -------------------------

% تحويل التاريخ إلى اليوم الجولياني
JD = juliandate(year, month, day);

% حساب الزمن النجمي الأوسط في غرينتش
GMST = mod(18.697374558 + 24.06570982441908 * (JD - 2451545), 24);

% حساب زاوية الميل (declination)
delta_rad = deg2rad(23.45) * sin(deg2rad(360 * (284 + N) / 365));

% حساب زاوية ساعة غروب الشمس (omega_s)
cos_omega_s = -tan(phi_rad) * tan(delta_rad);
omega_s_rad = acos(cos_omega_s);

% مدة النهار (بالساعات)
day_length_hours = 2 * rad2deg(omega_s_rad) / 15;

% وقت الظهر الحقيقي (Solar Noon)
Eqtime = E; % معادلة الزمن بالدقائق
solnoon_frac = (720 - 4 * L_deg - Eqtime + Timezone * 60) / 1440;
solnoon = solnoon_frac * 24;

% وقت شروق الشمس
sunrise_time = solnoon - day_length_hours / 2;
sunrise_hour = floor(sunrise_time);
sunrise_minute = floor((sunrise_time - sunrise_hour) * 60);

% وقت غروب الشمس
sunset_time = solnoon + day_length_hours / 2;
sunset_hour = floor(sunset_time);
sunset_minute = floor((sunset_time - sunset_hour) * 60);

% ------------------------- تحضير المتجه الزمني -------------------------

% تحديد زمن البداية والنهاية للحسابات (فترة النهار)
start_time = sunrise_time;
end_time = sunset_time;

% إنشاء متجه الزمن خلال فترة النهار فقط بفواصل 10 ثوانٍ
TL = start_time:10/3600:end_time;

% التحقق من أن TL ليس فارغًا
if isempty(TL)
    error('لا يوجد إشعاع شمسي في هذا اليوم (الشمس لا تشرق).');
end

num_steps = length(TL);

% تهيئة المصفوفات
I_bh = zeros(1, num_steps);
I_dh = zeros(1, num_steps);
I_gh = zeros(1, num_steps);
I_b_inclined = zeros(1, num_steps);
I_d_inclined = zeros(1, num_steps);
I_g_inclined = zeros(1, num_steps);

% ------------------------- حلقة الحسابات -------------------------

for i = 1:num_steps
    % التوقيت المحلي الحالي (TL)
    LT = TL(i);
    
    % الزمن الشمسي الحقيقي (TST)
    TC = E + 4 * (L_deg - 15 * Timezone); % بالدقائق
    TST = LT + TC / 60; % تحويل الدقائق إلى ساعات
    
    % زاوية الساعة (omega)
    omega_deg = 15 * (TST - 12);
    omega_rad = deg2rad(omega_deg);
    
    % زاوية ارتفاع الشمس (h)
    sin_h = sin(phi_rad) * sin(delta_rad) + cos(phi_rad) * cos(delta_rad) * cos(omega_rad);
    h_rad = asin(sin_h);
    h_deg = rad2deg(h_rad);
    
    % تحقق من أن الشمس فوق الأفق
    if sin_h > 0
        % كتلة الهواء النسبية (m)
        m = 1 / (sin_h + 0.50572 * (6.07995 + h_deg)^-1.6364);
        
        % معامل التوهين (tau)
        tau_b = exp(-0.128 * m); % قيمة تقريبية لمعامل التوهين المباشر
        tau_d = 0.271 - 0.294 * tau_b; % معامل التوهين المنتشر
        
        % الإشعاع خارج الغلاف الجوي (I_0)
        I_0 = I_sc * (1 + 0.033 * cos(deg2rad(360 * N / 365)));
        
        % الإشعاع المباشر على السطح الأفقي
        I_bh(i) = I_0 * tau_b * sin_h;
        
        % الإشعاع المنتشر على السطح الأفقي
        I_dh(i) = I_0 * tau_d;
        
        % الإشعاع الكلي على السطح الأفقي
        I_gh(i) = I_bh(i) + I_dh(i);
        
        % زاوية السقوط على السطح المائل
        cos_theta_tilted = sin(delta_rad) * sin(phi_rad - beta_rad_inclined) + cos(delta_rad) * cos(phi_rad - beta_rad_inclined) * cos(omega_rad);
        cos_theta_tilted = max(cos_theta_tilted, 0); % التأكد من أن القيمة غير سالبة
        
        % الإشعاع المباشر على السطح المائل
        I_b_inclined(i) = I_0 * tau_b * cos_theta_tilted;
        
        % الإشعاع المنتشر على السطح المائل
        I_d_inclined(i) = I_dh(i) * ((1 + cos(beta_rad_inclined)) / 2) + I_gh(i) * rho * ((1 - cos(beta_rad_inclined)) / 2);
        
        % الإشعاع الكلي على السطح المائل
        I_g_inclined(i) = I_b_inclined(i) + I_d_inclined(i);
    else
        % الشمس تحت الأفق
        I_bh(i) = 0;
        I_dh(i) = 0;
        I_gh(i) = 0;
        I_b_inclined(i) = 0;
        I_d_inclined(i) = 0;
        I_g_inclined(i) = 0;
    end
end

% ------------------------- حساب الإشعاع الشمسي اليومي الكلي -------------------------

% تحويل الإشعاع بوحدة W/m^2 إلى Wh/m^2
time_interval_hours = (10/3600); % الفاصل الزمني بالساعة (10 ثوانٍ)

% حساب الإشعاع الشمسي اليومي الكلي على السطح الأفقي
daily_irradiation_horizontal = sum(I_gh * time_interval_hours); % بوحدة Wh/m^2

% حساب الإشعاع الشمسي اليومي الكلي على السطح المائل
daily_irradiation_inclined = sum(I_g_inclined * time_interval_hours); % بوحدة Wh/m^2

% ------------------------- رسم المنحنيات -------------------------

% زمن الرسم
time_vector = TL;

% رسم منحنيات السطح الأفقي
figure(1);
plot(time_vector, I_bh, 'r', time_vector, I_dh, 'b', time_vector, I_gh, 'g', 'LineWidth', 1.5);
xlabel('الوقت (بالساعات)');
ylabel('الإشعاع الشمسي (W/m^2)');
title(['الإشعاع الشمسي على سطح أفقي - ورقلة - ', num2str(day), '/', num2str(month), '/', num2str(year)]);
legend('I_{bh} (مباشر أفقي)', 'I_{dh} (منتشر أفقي)', 'I_{gh} (كلي أفقي)');
grid on;

% رسم منحنيات السطح المائل
figure(2);
plot(time_vector, I_b_inclined, 'r', time_vector, I_d_inclined, 'b', time_vector, I_g_inclined, 'g', 'LineWidth', 1.5);
xlabel('الوقت (بالساعات)');
ylabel('الإشعاع الشمسي (W/m^2)');
title(['الإشعاع الشمسي على سطح مائل (\beta = ', num2str(beta_deg_inclined), '°) - ورقلة - ', num2str(day), '/', num2str(month), '/', num2str(year)]);
legend('I_b (مباشر مائل)', 'I_d (منتشر مائل)', 'I_g (كلي مائل)');
grid on;

% ------------------------- عرض النتائج -------------------------

% إيجاد القيمة القصوى للإشعاع الكلي على السطح الأفقي والمائل
[max_Igh, idx_max_horiz] = max(I_gh);
[max_Ig_inclined, idx_max_inclined] = max(I_g_inclined);

% تحويل وقت الذروة إلى ساعات ودقائق
[max_time_hour, max_time_minute] = decimalToHourMinute(TL(idx_max_horiz));

fprintf('تاريخ: %02d/%02d/%04d\n', day, month, year);
fprintf('مدة النهار: %.2f ساعة\n', day_length_hours);
fprintf('وقت شروق الشمس: %02d:%02d\n', sunrise_hour, sunrise_minute);
fprintf('وقت غروب الشمس: %02d:%02d\n', sunset_hour, sunset_minute);
fprintf('أقصى إشعاع كلي على السطح الأفقي: %.2f W/m^2 عند %02d:%02d (LT)\n', max_Igh, max_time_hour, max_time_minute);
fprintf('الإشعاع الشمسي اليومي الكلي على السطح الأفقي: %.2f Wh/m^2\n', daily_irradiation_horizontal);
fprintf('الإشعاع الشمسي اليومي الكلي على السطح المائل (\beta = %d°): %.2f Wh/m^2\n', beta_deg_inclined, daily_irradiation_inclined);

% ------------------------- تعريف الدوال المساعدة -------------------------

% دالة لحساب رقم اليوم في السنة
function N = dayOfYear(year, month, day)
    N = datenum(year, month, day) - datenum(year, 1, 0);
end

% دالة لتحويل الوقت العشري إلى ساعات ودقائق
function [hour, minute] = decimalToHourMinute(time_decimal)
    hour = floor(time_decimal);
    minute = floor((time_decimal - hour) * 60);
    if minute >= 60
        minute = minute - 60;
        hour = hour + 1;
    end
    if hour >= 24
        hour = hour - 24;
    end
end

% دالة لحساب اليوم الجولياني
function JD = juliandate(year, month, day)
    if month <= 2
        year = year - 1;
        month = month + 12;
    end
    A = floor(year / 100);
    B = 2 - A + floor(A / 4);
    JD = floor(365.25 * (year + 4716)) + floor(30.6001 * (month + 1)) + day + B - 1524.5;
end