      SUBROUTINE R_pure_Dmatr(RMAT,MSIZE,ROUTIN)
C
C	########################################################################
C
C     Valerio Cappellini                                         February 2007
C
C     USES :: rgnf_lux
C     COMPATIBILITY :: Fortran 90
C     VERSION :: 1.0.1
C
C     Generator of a real density matrix representing a random pure state.
C   	It takes as input:
C
C	+ RMAT  :: a real matrix of size MSIZE
C	+ MSIZE :: size of the real matrix RMAT
C
C     and gives as output the same matrix RMAT, now a density matrix,
C	representing a random pure state as the kind of the generic outputs 
C	of the subroutine "R_pure_vect.f".
C     Of course RMAT is Symmetric and of trace 1.
C
C	This code uses an auxiliary (real) vector RVEC of size MSIZE
C
C     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     ROUTIN   : Name of the subroutine you are going to use for generating
C                the random numbers (suggested :: rgnf_lux.f)
C
C     IMPORTANT [1] !!! Remember to initialize the Subroutine rgnf_lux.f
C                       by means of "Call Ranset_()"
C     IMPORTANT [2] !!! Remember to put the file "i_seed_rng" into the 
C                       working directory.
C     IMPORTANT [3] !!! Remember to insert main_program: "EXTERNAL rgnf_lux"
C
C     ########################################################################
C
$	debug
      SAVE  S, T, A, B, R1, R2
	INTEGER MSIZE, IDEV, JDEV
      REAL U(2), RVEC(MSIZE), RMAT(MSIZE,MSIZE), S, T, A, B, R1, R2
	REAL V, X, Y, Q, DEVIAT, SOMMA
      EXTERNAL ROUTIN
      DATA  S, T, A, B / 0.449871, -0.386595, 0.19600, 0.25472/
      DATA  R1, R2 / 0.27597, 0.27846/
C         generate pair of uniform deviates

	SOMMA = 0.
      DO 200 IDEV = 1, MSIZE
   50 CALL ROUTIN(U,2)
      V = 1.7156 * (U(2) - 0.5)
      X = U(1) - S
      Y = ABS(V) - T
      Q = X**2 + Y*(A*Y - B*X)
C           accept P if inside inner ellipse
      IF (Q .LT. R1)  GO TO 100
C           reject P if outside outer ellipse
      IF (Q .GT. R2)  GO TO 50
C           reject P if outside acceptance region
      IF (V**2 .GT. -4.0 *ALOG(U(1)) *U(1)**2)  GO TO 50
C           ratio of P's coordinates is normal deviate
  100 DEVIAT = V/U(1)
	SOMMA = SOMMA + DEVIAT * DEVIAT
  200 RVEC(IDEV) = DEVIAT

      DO IDEV = 1, MSIZE
		RMAT(IDEV,IDEV) = RVEC(IDEV) * RVEC(IDEV) / SOMMA
	    DO JDEV = IDEV+1, MSIZE
			RMAT(IDEV,JDEV) = RVEC(IDEV) * RVEC(JDEV) / SOMMA
			RMAT(JDEV,IDEV) = RMAT(IDEV,JDEV)
		ENDDO
	ENDDO

      RETURN
      END
