! Version 1.0 released by Yi-Chuan Lu on May 20, 2022.
! 
! When using this code, please cite:
! 
! @article{20heatindex,
!   Title   = {Extending the Heat Index},
!   Author  = {Yi-Chuan Lu and David M. Romps},
!   Journal = {Journal of Applied Meteorology and Climatology},
!   Year    = {2022},
!   Volume  = {61},
!   Number  = {10},
!   Pages   = {1367--1383},
!   Year    = {2022},
! }
!
! This headindex function returns the Heat Index in Kelvin. The inputs are:
! - T, the temperature in Kelvin
! - RH, the relative humidity, which is a value from 0 to 1

program test
    use heatindex_mod
    implicit none
    integer           :: i, j
    double precision  :: error = 1d-8
    double precision, dimension(18)    :: Ta    = (/200.,210.,220.,230.,240.,250.,260.,270.,280.,&
                                                    290.,300.,310.,320.,330.,340.,350.,360.,370./)
    double precision, dimension(3)     :: RH    = (/0.,0.5,1./) 
    double precision, dimension(3,18)  :: HI    = reshape((/&
                                                199.9994020652,199.9997010342,200.0000000021,&
                                                209.9975943902,209.9987971085,209.9999998068,&
                                                219.9915822029,219.9957912306,219.9999999912,&
                                                229.9739691979,229.9869861009,230.0000001850,&
                                                239.9253828022,239.9626700074,240.0000000003,&
                                                249.7676757244,249.8837049107,250.0000000037,&
                                                259.3735990024,259.6864068902,259.9999999944,&
                                                268.5453870455,269.2745889562,270.0000002224,&
                                                277.2234200026,278.6369451963,280.0000000091,&
                                                285.7510545370,288.2813660100,290.7860610129,&
                                                297.5737503539,300.2922595865,305.3947127590,&
                                                305.5549530893,318.6225524695,359.9063248191,&
                                                313.0298872791,359.0538750602,407.5345212438,&
                                                320.5088548469,398.5759733823,464.9949352940,&
                                                328.0358006469,445.8599463105,530.5524786708,&
                                                333.2806160592,500.0421800191,601.9518435268,&
                                                343.6312984164,559.6640227151,677.2462089759,&
                                                354.1825692377,623.1960299857,755.0832658147/), (/3,18/))
    logical, dimension(18,3)  :: truth_table
                                            
    do i = 1, 18
    do j = 1, 3
        truth_table(i,j) = ( abs(heatindex(Ta(i),RH(j))/HI(j,i) - 1.) < error )
    end do
    end do
    
    if (all(truth_table)) then
        print *,'Success'
    else
        print *,'Failure'
    end if

end program test