anchor.diff:  
       ! initialize amax
       amax = 1.0_real64
  -    log_fpxa_amax = log_fpxa(amax) ! log p(n|a_max)
  -    log_fpa_amax = log_fpa(amax) ! propto p(a)
  +    log_fwa_amax = log_fwa(amax) ! log p(n|a_max)
   
       ! set intervals equally spaced on log scale
       dx = (log_alpha2 - log_alpha1) / (nx * 1.0_real64)
  
positive.code:     real(real64) :: a
  
      a = exp(x)
      m_func = exp(log_fwa(a)) * hdir(a) * a
      
    end function m_func
  
  
negative.code:     ! find amax such that the weight in the integral - fwa(alpha) - is maximal
      i = maxloc(fxs,1,fxs < largest)
      amax = xs(i)
      log_fwa_amax = log_fwa(amax)
  
      ! recompute fxs
      fxs = exp(log_fwa(xs) - log_fwa_amax)
  
      ! re-compute a reasonable integration range
      log_alpha1 = log(minval(xs, fxs > small_number))
  