How to generate random numbers normally distributed in R or any software with limitations (bounds)? I am working on a project where I need to generate random numbers for a given task time which is normally distributed with mean = 40, and standard deviation = 150.
Because of the high SD, I will get some negative values and low values when I generated numbers directly which is unrealistic.
Is there any way where I can generate random numbers normally distributed with limitations (i.e. bounds)?
 A: This sounds like you want to sample from truncated normal distribution. If you only want to truncate the tails of the distribution (regions with low probability), than the approach suggested by Dave is probably enough. In other cases it might however quickly get inefficient. Better approach was suggested by Christian P. Robert in

Robert, C.P. (1995). Simulation of truncated normal variables.
Statistics and Computing 5(2): 121-125.

The algorithm is a bit more complicated, so I suggest you check the paper. Alternatively, if you prefer the code here you can find a C++ implementation from  R package extraDistr (disclaimer: it's written by me).
A: I could imagine something where you use some if/else logic to screen for unrealistic values. There would be some kind of recursion where you keep drawing random numbers u til you’ve gotten 1000 (or whatever) realistic values. Some pseudocode:
i=0
while i < 1000:
    x = make your draw here 
    # (np.random.normal or rnorm, for instance)
    if x is realistic:
        sample[i] = x
        increase i by 1

This will keep drawing random numbers for observation i until it gets a realistic value.
Note that you are not simulating normal data if you do this, as any real number is technically possible for any normal distribution, and you eliminate some values.
A: While your question is not entirely clear about what you are trying to achieve (how do you wish to go from a Gaussian distribution to a distribution that is truncated at 0?)...
... I thought that it was interesting to show something about the limit of the ratio between the mean and standard deviation of a Gaussian distribution that is truncated at $x=0$. (this issue has been mentioned in some of the comments)

Below is a piece of code and a graph that shows the computation of Gaussian distributions that have been truncated at different z-values (and shifted and re-scaled appropriately in order to have the truncating occur at $x=0$ and have the population mean equal to $\bar{x} = 40$).
What we can notice is that by changing the point where we truncate the distribution we can shift from a curve that looks like a Gaussian distribution (when we cut at a low z-value) to a curve that approaches an exponential distribution (when we cut at a high z-value and only have the right tail, which approximates an exponential function).
From this display I guess, intuitively, that the ratio of the standard deviation and the mean for this truncated distribution, is not able to become larger than this ratio for an exponential distribution (for an exponential distribution this ratio is 1).
Therefore: By truncating a normal distribution such that no negative values appear, we can not get a distribution whose standard deviation is larger than it's mean. (and you are looking for sd = 150 and mean = 40, which means that truncating a normal distribution won't be able to do it)

library(truncnorm)

x = seq(-10^3,10^3,0.1)

### empty canvas/plot
plot(-100,-100, 
     ylim = c(0,0.025), xlim = c(0,200),
     xlab = "x", ylab = "density")

d = 20 ### number of curves
i = 0  ### counter used in for-loop

varst = rep(0,d-1)

for (trunc in qnorm(seq(1/d,1-1/d,1/d))) {
  
  ### compute truncated standard normal
  ### and it's mean and variance
  y <- dtruncnorm(x, mean = 0, sd = 1, a = trunc)
  mean = dnorm(trunc)/(1-pnorm(trunc))
  var  = (1+trunc*dnorm(trunc)/(1-pnorm(trunc)) - mean^2)
  
  ### transform such that the mean is equal to 40
  xtrans <- (x-trunc)*40/(mean-trunc)
  ytrans <- y/(40/(mean-trunc))
  
  
  ### storing variance of transformed trucated standard normal (multipliying with square of scalefactor)
  varst[i+1] = var*(40/(mean-trunc))^2
  
  ### plot
  lines(xtrans[xtrans>=0],ytrans[xtrans>=0], 
        col = hsv(0.15+i/2/d,1-i/2/d,1-(d-i)/4/d,1))
  i = i+1
}

### exponential distribution
lines(x[x>=0],dexp(x,rate=1/40)[x>=0], lty = 2)

i = 1:(d-1)
legend(200,0.025,xjust = 1, cex = 0.7,
       legend = c("exponential distribution", "normal distribution cut at 5%", "normal distribution cut at 95%"),
       lty = c(2,1,1), col = c(1, hsv(0.15+i/2/d,1-i/2/d,1-(d-i)/4/d,1)[c(1,d-1)]))

For the equations used to compute the mean and variance of the truncated normal distribution see: https://en.wikipedia.org/wiki/Truncated_normal_distribution
