! module that contains some various useful constants and definitions
module constants

    ! import modern fortran kinds (real32/64, int8/16/32, etc.)
    use, intrinsic :: iso_fortran_env

    ! turn off f77 "feature" forcing counter variables (ijklmn) to be int
    implicit none

    ! define various variable types
    integer, parameter :: SP = real32  ! single-precision floating point
    integer, parameter :: DP = real64  ! double-precision floating point
    integer, parameter :: I1 = int8    !  8 bit (1 byte)  integer
    integer, parameter :: I2 = int16   ! 16 bit (2 bytes) integer
    integer, parameter :: I4 = int32   ! 32 bit (4 bytes) integer

    ! pi
    real(DP), parameter :: PI_DP = 4._DP * atan (1._DP)  ! double-precision

    ! define some constants for converting between degrees<-->radians    
    real(DP), parameter :: DEG_PER_RAD = 180._DP / PI_DP  ! degrees per radian
    real(DP), parameter :: RAD_PER_DEG = PI_DP / 180._DP  ! radians per degree


end module constants

! module that contains useful error checking code
module error_checking

    ! import the useful constants we've defined
    use constants

    ! turn off f77 "feature" forcing counter variables (ijklmn) to be int
    implicit none

    ! subroutine/function block
    contains

    ! error catching function
    subroutine catch(errorCode)

        ! define errorCode as an integer and is an input variable
        integer(I4), intent(in) :: errorCode

        ! declare a character array to hold our error message
        character(len=1024) :: errorMsg

        ! check if there is an error code (if it's not 0 or 1)
        if (errorCode /= 0 .AND. errorCode /= 1) then

            ! get the error message
            call slbm_get_error_message( errorMsg, sizeof(errorMsg) )

            ! print the error message
            print '(a,i0,a)', "ERROR ", errorCode, ":"
            print '(a)', errorMsg

            ! exit the program with error
            stop 1

        endif

    end subroutine catch
    
end module error_checking


! main progran
program main

    ! import statements
    use constants       ! some useful constants we've defined
    use error_checking  ! error checking functions

    ! turn off f77 "feature" forcing counter variables (ijklmn) to be int
    implicit none

    ! declare a whole bunch of variables
    integer(I4)        :: errorCode  ! integer error code
    character(len=80)  :: modelPath  ! /path/to/model_file
    integer(I4)        :: phase      ! phase id
    real(DP) :: srcLatDeg, srcLonDeg, srcDepKm,  &  ! source location
                rcvLatDeg, rcvLonDeg, rcvDepKm      ! receiver location
    real(DP) :: srcLatRad, srcLonRad,  &  ! source location (radians)
                rcvLatRad, rcvLonRad      ! source location (radians)
    real(DP) :: distDeg, distRad,      &  ! source <--> receiver distance
                travelTimeSec,         &  ! travel time (sec) from source --> receiver
                travelTimeUncertSec,   &  ! path-dependent travel time uncertainty (sec)
                travelTime1DUncertSec     ! old-style 1D travel time uncertainty (sec)

    ! declare the relative path to the model file we're going to be using
    modelPath = "../models/pdu202009Du.geotess"

    ! choose an arbitrary phase (Pn=0, Sn=1, Pg=2, or Lg=3)
    phase = 0

    ! arbitrary source (Meteor Crater, AZ)
    srcLatDeg =   35.0274_DP  ! latitude (degrees)
    srcLonDeg = -111.0228_DP  ! longitude (degrees)
    srcDepKm  =    1.2345_DP  ! depth (km)

    ! ! arbitrary receiver (Albuquerque, NM)
    rcvLatDeg =   35.1053_DP  ! latitude (degrees)
    rcvLonDeg = -106.6294_DP  ! longitude (degrees)
    rcvDepKm  =   -1.6000_DP  ! depth (km)

    ! convert lat/lon from degrees to radians
    srcLatRad = srcLatDeg * rad_per_deg
    srcLonRad = srcLonDeg * rad_per_deg
    rcvLatRad = rcvLatDeg * rad_per_deg
    rcvLonRad = rcvLonDeg * rad_per_deg

    ! instantiate an RSTT object
    call slbm_create(errorCode)
    call catch(errorCode)  ! check if the previous function returned an error

    ! load the velocity model
    call slbm_load_velocity_model(trim(modelPath), len_trim(modelPath), errorCode)
    call catch(errorCode)  ! check if the previous function returned an error

    ! create a great circle from source to the receiver
    call slbm_create_great_circle(phase,  &
        srcLatRad, srcLonRad, srcDepKm,   & 
        rcvLatRad, rcvLonRad, rcvDepKm,   &
        errorCode)
    call catch(errorCode)  ! check if the previous function returned an error

    ! get the distance and travel time from source --> receiver
    call slbm_get_distance(distRad, errorCode)           ! compute distance (rad) and store its value
    call catch(errorCode)  ! check if the previous function returned an error
    call slbm_get_travel_time(travelTimeSec, errorCode)  ! compute travel time (sec) and store its value
    call catch(errorCode)  ! check if the previous function returned an error
    distDeg = distRad * DEG_PER_RAD  ! convert radians --> degrees

    ! get the travel time uncertainty
    call slbm_get_tt_uncertainty(travelTimeUncertSec, errorCode)      ! compute uncertainty
    call catch(errorCode)  ! check if the previous function returned an error
    call slbm_get_tt_uncertainty1d(travelTime1DUncertSec, errorCode)  ! compute 1D uncertainty
    call catch(errorCode)  ! check if the previous function returned an error

    ! print our results
    print '(a)', ""
    print '(a)', "Fortran Example:"
    print '(a)', "--------------------------------"
    print '(a,f9.4,a)', "Distance:          ", distDeg, " deg"
    print '(a,f9.4,a)', "Travel Time:       ", travelTimeSec, " sec"
    print '(a,f9.4,a)', "Uncertainty:       ", travelTimeUncertSec, " sec"
    print '(a,f9.4,a)', "Uncertainty (1D):  ", travelTime1DUncertSec, " sec"
    print '(a)', ""

    ! return successfully
    return

end program main
