module Radiation

export AbstractRadiation

    abstract type AbstractRadiation end

    using GSL
    using Interpolations
    # const GSL_CONST = GSL.Const

    # const kpc =             1e3*GSL_CONST.CGSM_PARSEC
    # const cee =             GSL_CONST.CGSM_SPEED_OF_LIGHT
    # const emgm =            GSL_CONST.CGSM_MASS_ELECTRON
    # const pmgm =            GSL_CONST.CGSM_MASS_PROTON
    # const kboltz =	        GSL_CONST.CGSM_BOLTZMANN
    # const kboltz_kev2erg =  1.6022e-9 		  # Boltzman constant in keV/erg 
    # const emerg =           (GSL_CONST.CGSM_MASS_ELECTRON*GSL_CONST.CGSM_SPEED_OF_LIGHT^2)
    # const charg =           4.8e-10
    # const sigtom =          GSL_CONST.CGSM_THOMSON_CROSS_SECTION
    # const herg =            GSL_CONST.CGSM_PLANCKS_CONSTANT_H
    # const hkev =            (GSL_CONST.CGSM_PLANCKS_CONSTANT_H*6.2415e8)
    # const mjy =             1.e-26
    # const re0 =             2.81794e-13
    # const gconst =          GSL_CONST.CGSM_GRAVITATIONAL_CONSTANT
    # const sbconst =         GSL_CONST.CGSM_STEFAN_BOLTZMANN_CONSTANT
    # const aconst =          7.56e-15
    # const msun =            GSL_CONST.CGSM_SOLAR_MASS

    struct cyclosyn_emis_params
        nu::Float64
        b::Float64
        syn::Interpolations.Extrapolation
        derivs::Interpolations.Extrapolation
        # acc_syn::GSL.gsl_interp_accel
        # acc_derivs::GSL.gsl_interp_accel
    end

    struct cyclosyn_abs_params
        nu::Float64
        b::Float64
        syn::Interpolations.Extrapolation
        derivs::Interpolations.Extrapolation
        # acc_syn::gsl_interp_accel
        # acc_derivs::gsl_interp_accel
    end

    struct comint_params
        eph::Float64
        ephmin::Float64
        ephmax::Float64
        eldis::Interpolations.Extrapolation
        phodis::Interpolations.Extrapolation
        # acc_eldis::gsl_interp_accel
        # acc_phodis::gsl_interp_accel
    end

    struct comfnc_params
        game::Float64 
        e1::Float64 
        phodis::Interpolations.Extrapolation 
        # acc_phodis::gsl_interp_accel
    end

    struct disk_obs_params
        tin::Float64
        rin::Float64
        nu::Float64
    end

    struct disk_ic_params
        gamma::Float64
        beta::Float64
        tin::Float64
        rin::Float64
        rout::Float64
        h::Float64
        z::Float64	
        nu::Float64
    end

    function get_energy(rad::AbstractRadiation)     return rad.en_phot end
    function get_nphot(rad::AbstractRadiation)      return num_phot end
    function get_energy_obs(rad::AbstractRadiation) return en_phot_obs end
    function get_nphot_obs(rad::AbstractRadiation)  return num_phot_obs end
    function get_size(rad::AbstractRadiation)       return size end


    function set_beaming!(rad::AbstractRadiation, θ, speed, doppler)
        rad.angle = deg2rad(theta)
        rad.beta = speed
        rad.dopfac = doppler
    end

    function set_inclination!(rad::AbstractRadiation, θ)
        rad.angle = deg2rad(θ)
    end

    function set_geometry!(rad::AbstractRadiation, geom, l1, l2)

        if isequal(geom, "cylinder")
            r = l1
            z = l2
            vol = π*r^2*z
            dopnum = 2.0
        elseif isequal(geom, "sphere")
            r = l1
            z = r
            vol = (4/3)*π*r^3
            dopnum = 3.0
        else
            @warn "Input wrong, assuming sphere of radius 1 cm. \n Choose either sphere or cylinder"
            r = 1.0
            z = 1.0
            vol = 4π/3
            dopnum = 3.0
        end

        rad.r = r
        rad.z = z
        rad.vol = vol
        rad.dopnum = dopnum
        rad.geometry = geom
    end

    function integrated_lumunosity!(rad::AbstractRadiation, numin, numax)
        temp = 0.0
        for i ∈ 1:rad.size - 1
            if ((rad.en_phot_obs[i]  /rad.herg > numin) && 
                (rad.en_phot_obs[i+1]/rad.herg < numax))

                temp = temp + 0.5*(rad.en_phot_obs[i+1]/rad.herg - rad.en_phot_obs[i]/rad.herg)*
                                (rad.num_phot_obs[i+1] + rad.num_phot_obs[i])
            end
        end
        return temp
    end

    function set_counterjet!(rad::AbstractRadiation, flag)
        rad.counterjet = flag
    end


end