anchor.diff:    implicit none
   
     type modelpar
  -    real(fp) :: skybg = 0
  -    real(fp) :: alpha = 0.4
  -    real(fp) :: relabs = 0.75
  -    real(fp) :: hscale = 8000.0
  +    real(fp) :: skybg  = 1.5d-10
  +    real(fp) :: alpha  = 0.4
  +    real(fp) :: relabs = 0.5
  +    real(fp) :: hscale = 5000.0
  +    real(fp) :: beta   = 0.2
     end type
   
   contains
  
positive.code: 
      elemental real(fp) function E(cosg)
        real(fp), intent(in) :: cosg
        E = E0(cosg)
      end function
  
    end subroutine
  
negative.code: 
      elemental real(fp) function g(costh)
        real(fp), intent(in) :: costh
        g = 0.75_fp * (1 + costh**2)
      end function
  
      elemental real(fp) function E0(cosg) result(E)
  