Solar Radiation & Position

program solar_radiation
    implicit none

    ! Inputs
    double precision :: phi, L_loc, H, TZ, beta, gamma, rho_g
    integer :: n

    ! Calculated solar times
    double precision :: B, E, L_st, AST

    ! Calculated solar angles
    double precision :: delta, omega, alpha, gamma_s, theta

    ! Radiative rates
    double precision :: I_on, I_dn, I_beam, I_diffuse, I_albedo, I_total

    ! Helpers
    double precision :: pi, phi_rad, delta_rad, omega_rad, alpha_rad
    double precision :: beta_rad, gamma_rad, gamma_s_rad
    double precision :: sin_alpha, cos_gamma_s, cos_theta
    integer :: read_status

    pi = 3.141592653589793d0

    ! Read inputs from stdin
    read(*, *, iostat=read_status) phi
    if (read_status /= 0) then
        print *, "ERROR: Failed to read Latitude"
        stop
    endif
    read(*, *) L_loc
    read(*, *) n
    read(*, *) H
    read(*, *) TZ
    read(*, *) beta
    read(*, *) gamma
    read(*, *) rho_g

    ! Inputs validations
    if (phi < -90.0d0 .or. phi > 90.0d0) then
        print *, "ERROR: Latitude must be between -90 and 90 degrees"
        stop
    endif
    if (L_loc < -180.0d0 .or. L_loc > 180.0d0) then
        print *, "ERROR: Longitude must be between -180 and 180 degrees"
        stop
    endif
    if (n < 1 .or. n > 366) then
        print *, "ERROR: Day of year must be between 1 and 366"
        stop
    endif
    if (H < 0.0d0 .or. H > 24.0d0) then
        print *, "ERROR: Local Hour must be between 0 and 24"
        stop
    endif
    if (beta < 0.0d0 .or. beta > 90.0d0) then
        print *, "ERROR: Panel Tilt must be between 0 and 90 degrees"
        stop
    endif
    if (gamma < -180.0d0 .or. gamma > 180.0d0) then
        print *, "ERROR: Panel Azimuth must be between -180 and 180 degrees"
        stop
    endif
    if (rho_g < 0.0d0 .or. rho_g > 1.0d0) then
        print *, "ERROR: Ground Albedo must be between 0.0 and 1.0"
        stop
    endif

    ! Step 1: Equation of Time (E) and Solar Time (AST)
    B = dble(n - 1) * 360.0d0 / 365.0d0
    E = 229.2d0 * (0.000075d0 + 0.001868d0 * cos(B * pi / 180.0d0) &
        - 0.032077d0 * sin(B * pi / 180.0d0) &
        - 0.014615d0 * cos(2.0d0 * B * pi / 180.0d0) &
        - 0.040849d0 * sin(2.0d0 * B * pi / 180.0d0))

    L_st = 15.0d0 * TZ
    AST = H + (4.0d0 * (L_st - L_loc) + E) / 60.0d0

    ! Wrap AST
    do while (AST < 0.0d0)
        AST = AST + 24.0d0
    enddo
    do while (AST >= 24.0d0)
        AST = AST - 24.0d0
    enddo

    ! Step 2: Position Angles
    delta = 23.45d0 * sin(360.0d0 * dble(284 + n) / 365.0d0 * pi / 180.0d0)
    omega = 15.0d0 * (AST - 12.0d0)

    phi_rad = phi * pi / 180.0d0
    delta_rad = delta * pi / 180.0d0
    omega_rad = omega * pi / 180.0d0

    sin_alpha = sin(phi_rad) * sin(delta_rad) + cos(phi_rad) * cos(delta_rad) * cos(omega_rad)
    if (sin_alpha > 1.0d0) sin_alpha = 1.0d0
    if (sin_alpha < -1.0d0) sin_alpha = -1.0d0

    alpha_rad = asin(sin_alpha)
    alpha = alpha_rad * 180.0d0 / pi

    if (alpha > 0.0d0) then
        cos_gamma_s = (sin(alpha_rad) * sin(phi_rad) - sin(delta_rad)) / (cos(alpha_rad) * cos(phi_rad))
        if (cos_gamma_s > 1.0d0) cos_gamma_s = 1.0d0
        if (cos_gamma_s < -1.0d0) cos_gamma_s = -1.0d0
        gamma_s_rad = acos(cos_gamma_s)

        ! Azimuth is negative in morning, positive in afternoon (relative to South)
        if (omega < 0.0d0) then
            gamma_s = -gamma_s_rad * 180.0d0 / pi
        else
            gamma_s = gamma_s_rad * 180.0d0 / pi
        endif
    else
        gamma_s = 0.0d0
    endif

    ! Angle of incidence
    beta_rad = beta * pi / 180.0d0
    gamma_rad = gamma * pi / 180.0d0
    gamma_s_rad = gamma_s * pi / 180.0d0

    cos_theta = sin(alpha_rad) * cos(beta_rad) + cos(alpha_rad) * sin(beta_rad) * cos(gamma_s_rad - gamma_rad)
    if (cos_theta > 1.0d0) cos_theta = 1.0d0
    if (cos_theta < -1.0d0) cos_theta = -1.0d0
    theta = acos(cos_theta) * 180.0d0 / pi

    ! Step 3: Irradiation calculations
    if (alpha > 0.0d0) then
        I_on = 1367.0d0 * (1.0d0 + 0.033d0 * cos(360.0d0 * dble(n) / 365.0d0 * pi / 180.0d0))
        I_dn = I_on * exp(-0.128d0 / max(0.001d0, sin(alpha_rad)))

        if (cos_theta > 0.0d0) then
            I_beam = I_dn * cos_theta
        else
            I_beam = 0.0d0
        endif

        I_diffuse = I_dn * sin(alpha_rad) * 0.3d0 * (1.0d0 + cos(beta_rad)) / 2.0d0
        I_albedo = I_dn * sin(alpha_rad) * rho_g * (1.0d0 - cos(beta_rad)) / 2.0d0
    else
        I_beam = 0.0d0
        I_diffuse = 0.0d0
        I_albedo = 0.0d0
    endif

    I_total = I_beam + I_diffuse + I_albedo

    ! Write outputs
    print *, "AST=", AST
    print *, "E_MIN=", E
    print *, "DEC=", delta
    print *, "HR_ANGLE=", omega
    print *, "ALT=", alpha
    print *, "AZIMUTH=", gamma_s
    print *, "INCIDENCE=", theta
    print *, "I_BEAM=", I_beam
    print *, "I_DIFFUSE=", I_diffuse
    print *, "I_ALBEDO=", I_albedo
    print *, "I_TOTAL=", I_total

end program solar_radiation