How to poperly sample truncated distirbutions

I am trying to learn how to sample truncated distributions. To begin with I decided to try a simple example I found here example

I didn't really understand the division by the CDF, therefore I decided to tweak the algorithm a bit. Being sampled is an exponential distribution for values x>0 Here is an example python code:

# Sample exponential distribution for the case x>0
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def pdf(x):
        return x*np.exp(-x)

xvec=np.zeros(1000000)
x=1.
for i in range(1000000):
      a=x+np.random.normal()
      xs=x
      if a > 0. :
        xs=a
      A=pdf(xs)/pdf(x)
      if np.random.uniform()<A :
        x=xs
        xvec[i]=x

x=np.linspace(0,15,1000)
plt.plot(x,pdf(x))
plt.hist([x for x in xvec if x != 0],bins=150,normed=True)
plt.show()

Ant the output is:

The code above seems to work fine only for when using the condition if a > 0. :, i.e. positive x, choosing another condition (e.g. if a > 0.5 :) produces wrong results.

Since my final goal was to sample a 2D-Gaussian - pdf on a truncated interval I tried extending the simple example using the exponential distribution (see the code below). Unfortunately, since the simple case didn't work, I assume that the code given below would yield wrong results.

I assume that all this can be done using the advanced tools of python. However, since my primary idea was to understand the principle behind, I would greatly appreciate your help to understand my mistake. Thank you for your help.

# sample truncated Gaussian 
import numpy as np
from scipy.stats import multivariate_normal
RANGE=100000

a=3.01467E-02
b=1.32786E+00
a_range=[0.001,0.5]
b_range=[0.01, 2.5]
cov=[[8.8127990E-06,  5.3212327E-04],[  5.3212327E-04,  3.2322653E-02]]
cov_2=[[8.8127990E-06,  0.],[ 0.,  3.2322653E-02]]
x=a
y=b
for i in range(RANGE):
    dx,dy=np.random.multivariate_normal([a,b],cov_2)
    a_t=x+dx
    b_t=y+dy
    xs=x
    ys=y
# accept if within bounds
    if a_range[0]<a_t and a_t<a_range[1] and b_range[0]<b_t and b_t<b_range[1]:
        xs=a_t
        ys=b_t
    fr=multivariate_normal([a,b],cov)
    A=fr.pdf([xs,ys])/fr.pdf([x,y])
    if  np.random.uniform()<A :
        x=xs
        y=ys