anchor.diff:      
       fxs = log_fwa(xs)
       amax = xs(maxloc(fxs,1,fxs < largest))
  -    log_fpxa_amax = log_fpxa(amax)
  -    log_fpa_amax = log_fpa(amax)
  +    log_fwa_amax = log_fwa(amax)
       
     end subroutine compute_integration_range
   
  
positive.code: 
      real(real64), intent(in) :: alpha
  
      log_fwa = (log_fpa(alpha) - log_fpa_amax) + (log_fpxa(alpha) - log_fpxa_amax)
  
    end function log_fwa
  
  
negative.code:     real(real64) :: a
  
      a = exp(x)
      m_func = exp(log_fwa(a) - log_fwa_amax) * hdir(a) * a
      
    end function m_func
  
  