16QAM BER with Gray Mapping

Derive the theoretical 16QAM bit error rate (BER) with Gray coded constellation mapping in AWGN

The Matlab/Octave script performs the following:
(a) Generation of random binary sequence
(b) Assigning group of 4 bits to each 16-QAM constellation symbol per the Gray mapping
(c) Addition of white Gaussian Noise
(d) Demodulation of 16-QAM symbols and
(e) De-mapping per decimal to Gray conversion
(f) Counting the number of bit errors
(g) Running this for each value of Eb/No in steps of 1dB.
For the discussion on the theoretical BER for 16QAM, please refer to the post
http://www.dsplog.com/2008/06/05/16qam-bit-error-gray-mapping/