module uncer_subs

  use cmdfit_subs

  implicit none

contains

  subroutine create_cum_prob(tau2_grid, tau2_cum, cum_prob, one_sig_tau2)

    ! Creates the plot of tau^2 against cumulative probability.

    ! The input tau2 grid.
    real, dimension(:,:), intent(in) :: tau2_grid
    ! The output plot.
    real, allocatable, dimension(:), intent(out) :: tau2_cum, cum_prob
    ! And the tau2 corresponding to one sigma.
    real :: one_sig_tau2

    ! Local variables.
    integer, dimension(2) :: naxis
    integer :: icount, i1, i2
    real :: tau2_min, tot_prob
    character(len=10) :: flag

    naxis(1)=size(tau2_grid,1)
    naxis(2)=size(tau2_grid,2)
    allocate(tau2_cum(naxis(1)*naxis(2)), cum_prob(naxis(1)*naxis(2)))

    icount=0
    do i2=1, naxis(2)
      do i1=1, naxis(1)
        icount=icount+1
        tau2_cum(icount)=tau2_grid(i1,i2)
      end do
    end do
    call tau_sort(tau2_cum)

    ! The numerics can fail if tau**2 is too high.
    tau2_min=tau2_cum(1)

    tot_prob=sum(exp(-0.5*(tau2_cum-tau2_min)))

    do icount=1, naxis(1)*naxis(2)
      cum_prob(icount)=&
      100.0*sum(exp(-0.5*(tau2_cum(1:icount)-tau2_min)))/tot_prob
    end do

    ! Find the 1 sigma contour by quadratic interpolation.
    one_sig_tau2=quadint(cum_prob, tau2_cum, 68.26, flag)

  end subroutine create_cum_prob

  subroutine conf(one_sig_tau2, tau2_grid, axgrid, naxis)

    ! The one sigma level for tau^2
    real, intent(in) :: one_sig_tau2
    ! The tau^2 grid and associated axis.
    real, intent(in), dimension(:) :: tau2_grid, axgrid
    integer, intent(in) :: naxis

    real :: tau_level, test, work
    real, parameter :: conf_level=68.26
    integer :: i, ilow, ihigh
    real :: low, high, best

    from_bottom: do i=1, naxis
      if (one_sig_tau2 > tau2_grid(i)) then
        low=axgrid(i)
        ilow=i
        exit from_bottom
      end if
    end do from_bottom
    from_top: do i=naxis, 1, -1
      if (one_sig_tau2 > tau2_grid(i)) then
        high=axgrid(i)
        ihigh=i
        best=axgrid(minloc(tau2_grid,1))
        print*, best, '+', high-best, low-best
        print*, 'i.e. ', low, ' -> ', high, '(=', ihigh-ilow, ' grid pixels.)'
        best=(low+high)/2.0
        print*, best, '+/-', (high-low)/2.0 ,'(symmetrised)'
        exit from_top
      end if
    end do from_top
  
  end subroutine conf

end module uncer_subs

program uncer

  ! Calculates the uncertainties using a new experimental method.

  use ark_file_io
  use uncer_subs
  use quad

  implicit none

  integer :: iflag
  real, allocatable, dimension(:,:) :: axgrid, tau2_grid
  real, allocatable, dimension(:) :: cum_prob, tau2_cum
  real, allocatable, dimension(:,:) :: oneDgrid
  integer, dimension(2) :: naxis

  integer :: i, i1, i2, iaxis
  real :: tot_prob
  ! The minimum tau**2, removed to make numerics work.
  real :: tau2_min

  real :: one_sig_tau2
  character(len=80) :: axisnam

  call nxtark_in('grid.fit')
  iflag=inpark(naxis, tau2_grid, axgrid)

  call create_cum_prob(tau2_grid, tau2_cum, cum_prob, one_sig_tau2)

  print*, 'The 68 percent confidence contour is at tau^2 =', one_sig_tau2

  call nxtark_out('uncer.out')
  call typark(2)
  iflag=makark(naxis(1)*naxis(2), cum_prob, tau2_cum)

  deallocate(cum_prob, tau2_cum)

  do iaxis=1, 2
    allocate(oneDgrid(naxis(iaxis),1))
    tau2_min=minval(tau2_grid)
    do i=1, naxis(iaxis)
      if (iaxis == 1) then
        oneDgrid(i,1)=&
        -2.0*log(sum(exp(-0.5*(tau2_grid(i,1:naxis(2))-tau2_min))))
      else if (iaxis == 2) then
        oneDgrid(i,1)=&
        -2.0*log(sum(exp(-0.5*(tau2_grid(1:naxis(1),i)-tau2_min))))
      end if
    end do
    call create_cum_prob(oneDgrid, tau2_cum, cum_prob, one_sig_tau2)
    print*, ' '
    if (iaxis == 1) then 
      iflag=get_header_s('CTYPE1', axisnam)
      if (iflag < 0) axisnam='AXIS 1'
    else if (iaxis == 2) then
      iflag=get_header_s('CTYPE2', axisnam)
      if (iflag < 0) axisnam='AXIS 2'
    end if
    print*, 'For ', trim(axisnam)
    call conf(one_sig_tau2, oneDgrid(:,1), axgrid(:,iaxis), naxis(iaxis))
    call nxtark_out('tau2.'//axisnam(1:3))
    call typark(2)
    iflag=makark(naxis(iaxis), oneDgrid(:,1), axgrid(1:naxis(iaxis),iaxis))
    deallocate(cum_prob, tau2_cum, oneDgrid)
  end do

end program uncer
