program blackbody_radiation
    implicit none

    ! Inputs
    integer :: mode
    double precision :: T1_C, T2_C
    double precision :: Area, F12

    ! Constants
    double precision, parameter :: SIGMA = 5.670374d-8
    double precision, parameter :: WIEN_CONST = 2897.8d0
    double precision, parameter :: C1 = 3.74177d8   ! W.um^4 / m^2
    double precision, parameter :: C2 = 14387.8d0   ! um.K

    ! Calculated properties
    double precision :: T1_K, T2_K
    double precision :: Eb1, Eb2
    double precision :: lambda_max1, lambda_max2
    double precision :: Q12, F_eff

    ! Wavelength profile variables
    integer :: i, n_points
    double precision :: lambda, d_lambda, max_lambda
    double precision :: Eb_lam1, Eb_lam2
    double precision :: exp_term1, exp_term2
    integer :: iostat_val

    ! Read inputs
    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid calculation mode.'
        stop
    end if

    read(*,*,iostat=iostat_val) T1_C
    read(*,*,iostat=iostat_val) T2_C
    read(*,*,iostat=iostat_val) Area
    read(*,*,iostat=iostat_val) F12
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read calculation parameters.'
        stop
    end if

    ! Valider temperatures
    T1_K = T1_C + 273.15d0
    T2_K = T2_C + 273.15d0

    if (T1_K <= 0.0d0 .or. T2_K <= 0.0d0) then
        write(*,*) 'ERROR: Temperatures must be above absolute zero.'
        stop
    end if

    if (Area <= 0.0d0) then
        write(*,*) 'ERROR: Area must be positive.'
        stop
    end if

    if (F12 < 0.0d0 .or. F12 > 1.0d0) then
        write(*,*) 'ERROR: View factor must be between 0.0 and 1.0.'
        stop
    end if

    ! Determine effective view factor based on geometry mode
    ! 1 = Basic Blackbody (no exchange)
    ! 2 = Infinite Parallel Plates (F12 = 1.0)
    ! 3 = Finite Parallel Plates (F12 custom)
    ! 4 = Concentric Cylinders (F12 = 1.0)
    ! 5 = Concentric Spheres (F12 = 1.0)
    select case (mode)
    case (1)
        F_eff = 0.0d0
    case (2, 4, 5)
        F_eff = 1.0d0
    case (3)
        F_eff = F12
    case default
        write(*,*) 'ERROR: Invalid geometry mode selected.'
        stop
    end select

    ! Calculate total emissive power (Stefan-Boltzmann)
    Eb1 = SIGMA * T1_K**4
    Eb2 = SIGMA * T2_K**4

    ! Calculate peak wavelength (Wien's Law)
    lambda_max1 = WIEN_CONST / T1_K
    lambda_max2 = WIEN_CONST / T2_K

    ! Calculate radiative heat exchange
    Q12 = Area * F_eff * SIGMA * (T1_K**4 - T2_K**4)

    ! Output results
    write(*,'(A)') '============================================================'
    write(*,'(A)') '      BLACKBODY RADIATION ANALYSIS ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,I2)')         '  Mode Code               = ', mode
    write(*,'(A,F12.2,A)')  '  Temperature 1 (T1)      = ', T1_C, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Temperature 2 (T2)      = ', T2_C, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Temperature 1 (T1_K)    = ', T1_K, ' K'
    write(*,'(A,F12.2,A)')  '  Temperature 2 (T2_K)    = ', T2_K, ' K'
    write(*,'(A,F12.4,A)')  '  Surface Area (A1)       = ', Area, ' m2'
    write(*,'(A,F12.4)')    '  Effective View Factor   = ', F_eff
    write(*,*)
    write(*,'(A)') '--- STEFAN-BOLTZMANN & WIEN RESULTS ------------------------'
    write(*,'(A,ES12.4,A)') '  T1 Emissive Power (Eb1) = ', Eb1, ' W/m2'
    write(*,'(A,ES12.4,A)') '  T2 Emissive Power (Eb2) = ', Eb2, ' W/m2'
    write(*,'(A,F12.4,A)')  '  T1 Peak Wavelength      = ', lambda_max1, ' um'
    write(*,'(A,F12.4,A)')  '  T2 Peak Wavelength      = ', lambda_max2, ' um'
    write(*,'(A,ES12.4,A)') '  Net Heat Transfer (Q12) = ', Q12, ' W'
    write(*,*)

    ! =======================================================
    ! PLANCK SPECTRAL DISTRIBUTION DATA
    ! =======================================================
    write(*,'(A)') '--- PLANCK SPECTRAL DISTRIBUTION ---------------------------'
    write(*,'(A)') '  lambda [um]   Eb_lambda1 [W/m2.um]   Eb_lambda2 [W/m2.um]'
    write(*,'(A)') '  ----------------------------------------------------------'

    ! Adapt max wavelength to temperatures to show a nice curve
    ! For cold temperatures, peak is at longer wavelengths.
    max_lambda = 3.0d0 * max(lambda_max1, lambda_max2)
    if (max_lambda < 15.0d0) max_lambda = 15.0d0
    if (max_lambda > 50.0d0) max_lambda = 50.0d0

    n_points = 60
    d_lambda = max_lambda / dble(n_points)

    do i = 1, n_points
        lambda = dble(i) * d_lambda
        
        ! Curve 1
        exp_term1 = C2 / (lambda * T1_K)
        if (exp_term1 > 80.0d0) then
            Eb_lam1 = 0.0d0
        else
            Eb_lam1 = C1 / (lambda**5 * (exp(exp_term1) - 1.0d0))
        end if

        ! Curve 2
        exp_term2 = C2 / (lambda * T2_K)
        if (exp_term2 > 80.0d0) then
            Eb_lam2 = 0.0d0
        else
            Eb_lam2 = C1 / (lambda**5 * (exp(exp_term2) - 1.0d0))
        end if

        write(*,'(F10.4,4X,ES15.6,4X,ES15.6)') lambda, Eb_lam1, Eb_lam2
    end do

    write(*,*)
    write(*,'(A)') '--- THEORY REFERENCE ----------------------------------------'
    write(*,'(A)') '  Stefan-Boltzmann constant sigma = 5.670374e-8 W/m2.K4'
    write(*,'(A)') '  Planck constants: C1 = 3.74177e8 W.um4/m2, C2 = 14387.8 um.K'
    write(*,'(A)') '  Peak wavelength matches Wien Displacement constant = 2897.8 um.K'

end program blackbody_radiation