rtrun <- function (mu, variance, a, b) 
{  
   sigma <- sqrt(variance)
   FA = pnorm(((a - mu)/sigma))
   FB = pnorm(((b - mu)/sigma))
   q<-runif(length(mu)) * (FB - FA) + FA
#   p<-qnorm(runif(length(mu)) * (FB - FA) + FA)

#   for(j in 1:length(mu)){while(q[j]==1){q[j]<-q[j]-(10^-16)
#		print(c(j,q[j]))
#		}
#   		while(q[j]==0){q[j]<-q[j]+(10^-200)
#	print(c(j,q[j]))
#	}}
q<-q-floor(q)*(10^-16)
q<-q+(1-ceiling(q))*(10^-200)
   		   p<-qnorm(q)
   return(mu + sigma * p)
}