module VoigtKampff
    use accuracy
    implicit none

    public 	VoigtKampffCollection

    private

    integer(ik),parameter :: EXPAND_SIZE = 50

    real(rk),parameter    ::    REFERENCE_NU0 = 1.0
    real(rk),parameter    ::    DISTANCE_MAGIC_NUMBER = 2.0

    type ::        VoigtKampffT
        private
        integer(ik)    ::    m_Npoints
        real(rk)    ::    m_res
        real(rk)    ::    m_lorentz_cutoff
        integer(ik)    ::    m_middle_point
        real(rk)     ::    m_gammaD
        real(rk)     ::    m_gammaL
        real(rk)    ::    m_mag
        logical        ::    normalize
        logical        ::    constructed = .false.
        real(rk),pointer    ::    m_voigt_grid(:)


    contains
        procedure,public    ::    construct    =>    construct_voigt
    
        procedure,public    ::    compute        =>     compute_voigt
    
        procedure,public        ::    destroy

        procedure,public    ::    get_gammaL

    end type    
    
    
    
    type ::	 VoigtKampffCollection

        type(VoigtKampffT),allocatable :: fast_voigts(:)
        integer(ik) :: capacity
        integer(ik) :: n
        logical     :: normalize
        real(rk)    :: dpwcoeff
        real(rk)    :: offset
        real(rk)    :: dfreq

    contains
        procedure,public  :: construct
        procedure,public  :: expand_voigts
        procedure,public  :: generate_indices
        procedure,public  :: get_size
        procedure,private :: add_voigt
        procedure,private :: search_gamma
        procedure,public  :: compute => compute_fast_voigt
    end type

    
    

contains

    subroutine construct(this,dpwcoeff,offset,dfreq,normalize)
        class(VoigtKampffCollection)	::	this
        logical,intent(in)		::	normalize
        real(rk),intent(in)		::	dpwcoeff,offset,dfreq

        this%n = 0
        this%capacity = 0

        this%dpwcoeff =dpwcoeff
        this%offset = offset
        this%dfreq = dfreq
        this%normalize = normalize

        call this%expand_voigts()

    end subroutine

    subroutine expand_voigts(this)
        class(VoigtKampffCollection) :: this
        type(VoigtKampffT),allocatable :: temp(:)
        integer(ik) :: old_capacity

        if(this%capacity == 0) then
            this%capacity = this%capacity + EXPAND_SIZE
            allocate(this%fast_voigts(this%capacity))
            return
        endif

        allocate(temp(this%capacity))
        temp(:) = this%fast_voigts(:)

        deallocate(this%fast_voigts)

        old_capacity=this%capacity

        this%capacity = this%capacity + EXPAND_SIZE

        allocate(this%fast_voigts(this%capacity))

        this%fast_voigts(1:old_capacity) = temp(1:old_capacity)

        deallocate(temp)

    end subroutine

    integer(ik) function get_size(this)
        class(VoigtKampffCollection),intent(in)	::	this

        get_size = this%n
    end function


    subroutine generate_indices(this,gammaL,gamma_idx,Jmax)
        class(VoigtKampffCollection) :: this
        real(rk),intent(in) :: gammaL
        integer(ik),intent(inout) :: gamma_idx
        integer(ik),intent(in) :: Jmax

        integer(ik) :: ido,jdo
        real(rk)    :: val
        integer     :: idx

        !Search of the gamma value, if its similar then resue the index otherwise we add antoher voigt

        val = gammaL
        idx = this%search_gamma(val)
              !If we didnt find it then add a voigt
        if(idx==0) then
           gamma_idx = this%add_voigt(val)
        else
           gamma_idx = idx
        endif

    end subroutine

    subroutine compute_fast_voigt(this,freq,intens,abscoef,ib,ie,start_nu,nu0,index)
        class(VoigtKampffCollection),intent(in)	::	this
        real(rk),intent(in)     :: freq(:),abscoef,start_nu,nu0
        integer(ik),intent(in)  :: ib,ie
        real(rk),intent(inout)  :: intens(:)
        integer,intent(in)      :: index
        if(index > this%n) stop "FastVoigt::Out of bounds"

        call this%fast_voigts(index)%compute(freq,intens,abscoef,ib,ie,start_nu,nu0)

    end subroutine

    integer(ik) function add_voigt(this,gammaL)
        class(VoigtKampffCollection) :: this
        real(rk),intent(in) :: gammaL

        this%n = this%n+1

        if(this%n > this%capacity) call this%expand_voigts

        add_voigt = this%n

        call this%fast_voigts(this%n)%construct(this%dpwcoeff,gammaL,this%dfreq,this%offset,this%normalize)
   
    end function

    integer(ik) function search_gamma(this,gammaL)
        class(VoigtKampffCollection) :: this
        real(rk),intent(in) :: gammaL

        search_gamma = 0

        if(this%n == 0) return
        search_gamma = binarySearch_R(this%fast_voigts(1:this%n),gammaL)


    end function

recursive function binarySearch_R (a, value) result (bsresult)
    type(VoigtKampffT), intent(in) :: a(:)
    real(rk)	::	value
    integer          :: bsresult, mid
 
    mid = size(a)/2 + 1
    if (size(a) == 0) then
        bsresult = 0        ! not found
    else if(is_close(a(mid)%get_gammaL(),value,1e-5_rk)) then
        bsresult = mid      ! SUCCESS!!
    else if (a(mid)%get_gammaL() > value) then
        bsresult= binarySearch_R(a(:mid-1), value)
    else if (a(mid)%get_gammaL() < value) then
        bsresult = binarySearch_R(a(mid+1:), value)
        if (bsresult /= 0) then
            bsresult = mid + bsresult
        end if
   endif
end function binarySearch_R



    !Check if the lorentzian is the same
    logical function is_close(a, b, rel_tol )
      real(rk),intent(in)	::	a,b,rel_tol
      	is_close = abs(a-b) <= rel_tol * max(abs(a), abs(b))
    end function



  subroutine construct_voigt(this,pGammaD,pGammaL,pRes,pLorentzCutoff,pNormalize)
      class(VoigtKampffT),intent(inout)    ::    this
      real(rk),intent(in)            ::    pGammaD,pGammaL
      real(rk),intent(in)            ::    pRes
      real(rk),intent(in)            ::    pLorentzCutoff
      logical,intent(in)            ::    pNormalize

      integer(ik)                ::    ierr,ido

      real(rk)                ::    nu,nu0

  
      this%m_res = pRes
      this%m_gammaD = pGammaD
      this%m_gammaL = pGammaL
      this%m_lorentz_cutoff = pLorentzCutoff
     this%normalize = pNormalize
      this%m_Npoints = int(2.0*this%m_lorentz_cutoff/this%m_res,ik);
      this%m_middle_point = (this%m_Npoints/2) + 1
      !Construct the voigt grid

      this%m_mag = 0.0

      allocate(this%m_voigt_grid(this%m_Npoints),stat=ierr)

      !call ArrayStart('Voigt-F-Grid',info,size(this%m_voigt_grid),kind(this%m_voigt_grid))
      nu0 = REFERENCE_NU0
      nu = nu0-this%m_lorentz_cutoff

      do ido=1,this%m_Npoints
      
          this%m_voigt_grid(ido) = voigt_humlicek(nu,nu0,this%m_gammaD,this%m_gammaL)
          this%m_mag = this%m_mag + this%m_voigt_grid(ido)
          nu = nu+ this%m_res
      enddo

      this%m_mag = this%m_mag*this%m_res

      this%constructed = .true.

  end subroutine

  subroutine destroy(this)
     class(VoigtKampffT),intent(inout)    ::    this
     
     if(.not. this%constructed)return
     
     if(associated(this%m_voigt_grid)) then
         deallocate(this%m_voigt_grid)
     endif
  end subroutine

  
  subroutine compute_voigt(this,freq,intens,abscoef,ib,ie,start_nu,nu0)
      class(VoigtKampffT),intent(in)    ::    this
      real(rk),intent(in)        ::    freq(:),abscoef,start_nu,nu0
      integer(ik),intent(in)        ::    ib,ie
      real(rk),intent(inout)        ::    intens(:)

  
      integer                ::    center_point,dist
      integer                ::    middle_shift
      integer                ::    start_dist,end_dist,left_start,left_end,right_start,right_end
      integer                ::    ib_rel,ie_rel,num_hum_points,ido
      integer                ::    Npoints
      real(rk)            ::    gammaD,hum_res,mag,nu
      real(rk),allocatable        ::    temp_humlicek(:)

  
      center_point = (nu0-start_nu)/this%m_res + 1
      middle_shift = (center_point - this%m_middle_point) + 1

  

      dist = max(DISTANCE_MAGIC_NUMBER/this%m_res,3.0)

      Npoints = ie - ib + 1

      ib_rel = ib - middle_shift
    ie_rel = ie - middle_shift

      gammaD = this%m_gammaD*nu0


    mag = this%m_mag

    start_dist = max(center_point - dist + 1, ib)
    end_dist = min(center_point + dist, ie)
    
    left_start = ib_rel
    left_end = min(this%m_middle_point - dist, ie_rel)
    right_start = max(this%m_middle_point + dist, ib_rel)
    right_end = ie_rel

    !print *,mag


    !If we have any humlicek points
    if(start_dist < end_dist) then


        if(this%normalize) then 
            num_hum_points = end_dist - start_dist + 1;
            allocate(temp_humlicek(num_hum_points))

            do ido=start_dist,end_dist
                nu=freq(ido)
                hum_res = voigt_humlicek(nu,nu0,gammaD,this%m_gammaL)
                temp_humlicek(ido-start_dist+1) = hum_res
                mag = mag - this%m_voigt_grid(ido -middle_shift+1)*this%m_res
                mag = mag + hum_res*this%m_res
            enddo
    
    
            intens(start_dist:end_dist) = intens(start_dist:end_dist) + temp_humlicek(:)*abscoef/mag
    
            !print *,mag
            deallocate(temp_humlicek)
        else
            do ido=start_dist,end_dist
                nu=freq(ido)
                intens(ido) = intens(ido) + voigt_humlicek(nu,nu0,gammaD,this%m_gammaL)*abscoef
            enddo    
        
            mag = 1.0
    
    
        endif
      endif
      if (left_start < left_end ) then
          call vectorized_voigt(intens(ib-ib_rel+left_start:ib-ib_rel+left_end),this%m_voigt_grid(left_start:left_end),(abscoef/mag))

      endif

      if (right_start < right_end ) then
          call vectorized_voigt(intens(ib-ib_rel+right_start:ib-ib_rel+right_end),this%m_voigt_grid(right_start:right_end),(abscoef/mag))
      endif

  end subroutine
          
  real(rk) function get_gammaL(this)
      class(VoigtKampffT),intent(in)    ::    this

      get_gammaL = this%m_gammaL

      return
  end function




!-----------Humlicek methods----------------------!

 ! This should be heavily vectorized by the compiler
  subroutine vectorized_voigt(intens,voigt,abscoeff)
      real(rk),intent(inout)    ::    intens(:)
      real(rk),intent(in)    ::    voigt(:)
      real(rk),intent(in)    ::    abscoeff

      intens(:) = intens(:)+voigt(:)*abscoeff
  end subroutine


  function voigt_humlicek(nu,nu0,alpha,gamma) result(f)
   !
   real(rk),intent(in) :: nu,nu0,alpha,gamma
   real(rk) :: f,sigma,x,y

   real(rk),parameter ::  RT2LN2 = 1.1774100225154747_rk     ! sqrt(2ln2)
   real(rk),parameter ::  RTPI   = 1.7724538509055159_rk     ! sqrt(pi)
   real(rk),parameter ::  RT2PI  = 2.5066282746310002_rk     ! sqrt(2pi)
   real(rk),parameter ::  RT2    = 1.4142135623730951_rk     ! sqrt(pi)
   !
   sigma = alpha / RT2LN2
   x = (nu0 - nu) / sigma / RT2
   y = gamma / sigma / RT2
   !
   f = humlic(x, y) / sigma / RT2PI
   !
  end function voigt_humlicek

  !
  function humlic(x,y) result (f)
    real(rk),intent(in) :: x,y
    real(rk) :: f
    complex(rk) :: t,humlic1,u
    real(rk) :: s
    !
    T = cmplx(y,-x)
    S = abs(x) + y
    !
    if (S >= 15) then
        ! Region I
        humlic1 = T*0.5641896_rk/(0.5_rk+T*T)
        !fprintf(stdout, "I")
        !
        f = real(humlic1,rk)
        return
    endif
    !
    if (S >= 5.5) then
        ! Region II
        U = T * T;
        humlic1 = T * (1.410474_rk + U*0.5641896_rk)/(0.75_rk + U*(3.0_rk+U))
        !
        !fprintf(stdout, "II");
        !
        f = real(humlic1,rk)
        return
        !
    endif
    !
    if (y >= 0.195_rk * abs(x) - 0.176_rk) then
        ! Region III
        humlic1 = (16.4955_rk+T*(20.20933_rk+T*(11.96482_rk &
                +T*(3.778987_rk+T*.5642236_rk)))) / (16.4955_rk+T*(38.82363_rk &
                +T*(39.27121_rk+T*(21.69274_rk+T*(6.699398_rk+T)))))
        !fprintf(stdout, "III");
        f = real(humlic1,rk)
        return
        !
    endif
    !
    ! Region IV
    U = T * T
    humlic1 = exp(U)-T*(36183.31_rk-U*(3321.9905_rk-U*(1540.787_rk-U*(219.0313_rk-U* &
       (35.76683_rk-U*(1.320522_rk-U*.56419_rk))))))/(32066.6_rk-U*(24322.84_rk-U* &
       (9022.228_rk-U*(2186.181_rk-U*(364.2191_rk-U*(61.57037_rk-U*(1.841439_rk-U)))))))
    !fprintf(stdout, "IV");
    !
    f = real(humlic1,rk)
    !
    end function humlic


end module VoigtKampff
