program multilayer_insulation_cryogenic
    implicit none
    integer :: iostat_val, N_shields
    double precision :: Th_hot_K, Tc_cold_K, eps_emiss, P_vac_Torr, thick_mm, A_tank_m2, hfg_kJkg
    double precision :: N_density_lcm, T_mean_K, DeltaT_K, thick_m
    double precision :: q_rad_Wm2, q_solid_Wm2, q_gas_Wm2, q_total_Wm2, keff_WmK, boiloff_kgday
    double precision, parameter :: SIGMA = 5.670374d-8

    ! Read inputs
    read(*,*,iostat=iostat_val) N_shields       ! Number of Shield Layers N (e.g. 30)
    read(*,*,iostat=iostat_val) thick_mm        ! Total Blanket Thickness [mm] (e.g. 15.0)
    read(*,*,iostat=iostat_val) Th_hot_K        ! Warm Boundary Temp [K] (e.g. 293.0)
    read(*,*,iostat=iostat_val) Tc_cold_K       ! Cold Cryo Boundary Temp [K] (e.g. 77.0 for LN2 or 20.0 for LH2)
    read(*,*,iostat=iostat_val) eps_emiss       ! Shield Emissivity (e.g. 0.03)
    read(*,*,iostat=iostat_val) P_vac_Torr      ! Residual Vacuum Pressure [Torr] (e.g. 1.0e-5)
    read(*,*,iostat=iostat_val) A_tank_m2       ! Tank Surface Area [m2] (e.g. 5.0)
    read(*,*,iostat=iostat_val) hfg_kJkg        ! Latent Heat of Cryogen [kJ/kg] (e.g. 199.0 for LN2, 446.0 for LH2)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for MLI calculation.'
        stop
    end if

    if (N_shields < 1 .or. thick_mm <= 0.0d0 .or. Th_hot_K <= Tc_cold_K .or. eps_emiss <= 0.0d0) then
        write(*,*) 'ERROR: Shield count, thickness, temperatures (Th > Tc), and emissivity must be valid.'
        stop
    end if

    thick_m = thick_mm * 1.0d-3
    N_density_lcm = dble(N_shields) / (thick_mm / 10.0d0) ! layers per cm
    DeltaT_K = Th_hot_K - Tc_cold_K
    T_mean_K = 0.5d0 * (Th_hot_K + Tc_cold_K)

    ! Modified Lockheed / McIntosh Semi-Empirical MLI Model
    ! 1. Radiation Component
    q_rad_Wm2 = (SIGMA * (Th_hot_K**4 - Tc_cold_K**4)) / (dble(N_shields + 1) * (2.0d0/eps_emiss - 1.0d0))

    ! 2. Solid Spacer Conduction Component (Dacron mesh)
    q_solid_Wm2 = (0.000073d0 * (N_density_lcm**2.63d0) * (T_mean_K) * DeltaT_K) / dble(N_shields)

    ! 3. Residual Gas Molecular Conduction Component
    q_gas_Wm2 = (14600.0d0 * max(1.0d-9, P_vac_Torr) * (Th_hot_K**0.52d0 - Tc_cold_K**0.52d0)) / dble(N_shields)

    q_total_Wm2 = q_rad_Wm2 + q_solid_Wm2 + q_gas_Wm2
    keff_WmK = (q_total_Wm2 * thick_m) / DeltaT_K

    ! Cryogen Boil-Off Rate [kg / day]
    boiloff_kgday = (q_total_Wm2 * A_tank_m2 * 86400.0d0) / (hfg_kJkg * 1000.0d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CRYOGENIC MULTILAYER INSULATION (MLI)'
    write(*,'(A)') '============================================================'
    write(*,'(A,I5,A,F10.2,A)')    'Shield Layers N / Thickness= ', N_shields, ' layers / ', thick_mm, ' mm'
    write(*,'(A,F10.1,A,F10.1,A)') 'Hot Temp Th / Cold Temp Tc = ', Th_hot_K, ' K / ', Tc_cold_K, ' K'
    write(*,'(A,F10.4,A,E12.3,A)') 'Emissivity / Vacuum Press  = ', eps_emiss, ' / ', P_vac_Torr, ' Torr'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4,A)')  'TOTAL MLI HEAT FLUX q"    = ', q_total_Wm2, ' W/m2'
    write(*,'(A,F10.4,A,F10.4,A)') 'Rad / Solid / Gas Flux     = ', q_rad_Wm2, ' / ', q_solid_Wm2, ' / ', q_gas_Wm2, ' W/m2'
    write(*,'(A,E12.4,A)')  'EFFECTIVE CONDUCTIVITY keff= ', keff_WmK, ' W/(m.K)'
    write(*,'(A,F10.2,A)')  'Cryogen Boil-Off Loss Rate = ', boiloff_kgday, ' kg/day'
    write(*,'(A)') '============================================================'

end program multilayer_insulation_cryogenic