anchor.diff:      
       ! calculate the distance from the wake center to the line perpendicular to the 
       ! line between the two circle intersection points
  -    if (OVdYd >= 0.0_dp + tol) then ! check case to avoid division by zero
  +    !if (OVdYd >= 0.0_dp + tol) then ! check case to avoid division by zero
  +    if (OVdYd >= 0.0_dp) then ! check case to avoid division by zero
           OVL = (-OVr*OVr+OVRR*OVRR+OVdYd*OVdYd)/(2.0_dp*OVdYd)
       else
           OVL = 0.0_dp
  
positive.code:     OVz = OVRR*OVRR-OVL*OVL
  
      ! Finish calculating the distance from the intersection line to the outer edge of the wake
      if (OVz > 0.0_dp + tol) then
          OVz = sqrt(OVz)
      else
          OVz = 0.0_dp
      end if
  
      if (OVdYd < (OVr+OVRR)) then ! if the rotor overlaps the wake
          !print *, "OVL: ", OVL
          !if (OVL < OVRR .and. (OVdYd-OVL) < OVr) then
          if (OVdYd > 0.0_dp + tol) then
              wake_overlap = OVRR*OVRR*dacos(OVL/OVRR) + OVr*OVr*dacos((OVdYd-OVL)/OVr) - OVdYd*OVz
          else if (OVRR > OVr) then
              wake_overlap = pi*OVr*OVr
              !print *, "wake ovl: ", wake_overlap
  
negative.code:     
  end subroutine discontinuity_point_func
  
  subroutine smooth_max(x, y, g)
  
      implicit none
      
      ! define precision to be the standard for a double precision ! on local system
      integer, parameter :: dp = kind(0.d0)
  
      ! in
      real(dp), intent(in) :: x, y
      
      ! out
      real(dp), intent(out) :: g
      
      ! local
      real(dp) :: s
      
      intrinsic log, exp
      
      s = 1000.0_dp
      
      g = (log(exp(s*x) + exp(s*y)))/s
  
      
  end subroutine smooth_max
      
   !    yd, n = _checkIfFloat(yd)
  ! 
  !     y1 = (1-pct_offset)*ymax
  !     y2 = (1+pct_offset)*ymax
  ! 
  !     dy1 = (1-pct_offset)
  !     dy2 = (1+pct_offset)
  ! 
  !     if (maxmin == 1) then
  !         f1 = y1
  !         f2 = ymax
  !         g1 = 1.0_dp
  !         g2 = 0.0_dp
  !         if (yd .ge. y2) then
  !             idx_constant = False
  !         else
  !             idx_constant = True
  !         end if
  ! 
  !         df1 = dy1
  !         df2 = 1.0_dp
  ! 
  ! 
  !     else if (maxmin == 0) then
  !         f1 = ymax
  !         f2 = y2
  !         g1 = 0.0_dp
  !         g2 = 1.0_dp
  !         if (yd .ge. y1) then
  !             idx_constant = False
  !         else
  !             idx_constant = True
  !         end if
  ! 
  !         df1 = 1.0_dp
  !         df2 = dy2
  !         
  !     end if
  ! 
  !     f = CubicSplineSegment(y1, y2, f1, f2, g1, g2)
  ! 
  !     # main region
  !     ya = np.copy(yd)
  !     if dyd is None:
  !         dya_dyd = np.ones_like(yd)
  !     else:
  !         dya_dyd = np.copy(dyd)
  ! 
  !     dya_dymax = np.zeros_like(ya)
  ! 
  !     # cubic spline region
  !     idx = np.logical_and(yd > y1, yd < y2)
  !     ya[idx] = f.eval(yd[idx])
  !     dya_dyd[idx] = f.eval_deriv(yd[idx])
  !     dya_dymax[idx] = f.eval_deriv_params(yd[idx], dy1, dy2, df1, df2, 0.0, 0.0)
  ! 
  !     # constant region
  !     ya[idx_constant] = ymax
  !     dya_dyd[idx_constant] = 0.0
  !     dya_dymax[idx_constant] = 1.0
  ! 
  !     if n == 1:
  !         ya = ya[0]
  !         dya_dyd = dya_dyd[0]
  !         dya_dymax = dya_dymax[0]
  ! 
  ! 
  !     return ya, dya_dyd, dya_dymax
  ! 
  ! 
  ! def smooth_max(yd, ymax, pct_offset=0.01, dyd=None):
  !     """array max, uses cubic spline to smoothly transition.  derivatives with respect to array and max value.
  !     width of transition can be controlled, and chain rules for differentiation"""
  !     return _smooth_maxmin(yd, ymax, 'max', pct_offset, dyd)
  ! 
  ! 
  ! def smooth_min(yd, ymin, pct_offset=0.01, dyd=None):
  !     """array min, uses cubic spline to smoothly transition.  derivatives with respect to array and min value.
  !     width of transition can be controlled, and chain rules for differentiation"""
  !     return _smooth_maxmin(yd, ymin, 'min', pct_offset, dyd)
  ! 
  