sggamma <- function(beta, sigma, lambda, timet)
{

s <- lambda^(-2)
r <- s * exp(( - beta * lambda)/sigma)

# consider timet below exp(-500) at exp(-500) to avoid division by zero
timet[timet<exp(-500)] <- exp(-500)

b <- lambda/sigma
surv <- (lambda>0)*(1-pgamma(timet^b, shape=s, rate=r))
surv <- surv+(lambda<0)*pgamma(timet^b, shape=s, rate=r)
surv[surv<exp(-500)] <- exp(-500)
surv
}