undef ("WEIBULL")
function WEIBULL(x[*]:numeric, scale[1]:numeric, shape[1]:numeric, weiDist[1]:integer) 
; source: http://support.sas.com/documentation/onlinedoc/ets/ex_code/131/svrtwbul.html
local eps, z, eps, one, Fill
begin
   
  if (.not.(typeof(x).eq."double".or.typeof(x).eq."float")) then
      print("WARNING: WEIBULL: function expects type double or float")
  end if

  if (typeof(x).eq."double") then
      eps  = 1d-14               ;  2.220446e-16
      one  = 1.0d
      Fill = 1d20
  else
      eps = 1e-6  
      one = 1.0
      Fill = 1e20
  end if

  z = (x/scale)^shape

  if (weiDist.eq.1) then        ; Weibull=>Log_PDF
      y = log(shape * z / x) - z
  end if

  if (weiDist.eq.2) then        ; Weibull=>PDF
      y = (shape / x) * z * exp(-z)
  end if

  if (weiDist.eq.3) then        ; Weibull=>Log_CDF
      c = exp(-z)
      y = where(c.lt.one, log(one-c), Fill) 
  end if

  if (weiDist.eq.4) then        ; Weibull=>CDF
      y = 1 - exp(-z)
  end if

  if (weiDist.eq.5) then        ; Weibull=>SDF
      y = exp(-z)
  end if

  return(y)
end
;---
undef ("WEIBULL_MEAN")
function WEIBULL_MEAN(scale[1]:numeric, shape[1]:numeric) 
begin
  if (typeof(scale).eq."double") then
      Fill = 1d20
  else
      Fill = 1e20
  end if

   if (shape.gt.0) then
       ymean = scale*gamma(1+1/shape)
   else
       ymean = Fill
       ymean@_FillValue = Fill
   end if
   ymean@long_name = "Weibull Mean derived from shape and scale"
   return(ymean)
end