      SUBROUTINE R_mixed_Dmatr(N,K,RMAT,ROUTIN)
C
C	########################################################################
C
C	Valerio Cappellini                                            March 2007
C
C     USES :: rgnf_lux
C     COMPATIBILITY :: Fortran 90
C     VERSION :: 1.0.0
C
C     Generator of a real density matrix representing mixed quantum states,
C	distributed according to the induced measure P_(N,K).
C
C	Such a states \sigma are generated by this algorithmic procedures:
C
C	- we generate a REAL random pure state \psi in the composite 
C	  Hilbert Space of dimension (N*K), according to the natural 
C	  Fubini-Study measure (e.g.: by using the subroutine R_pure_vect); 
C
C	- we represent \psi as a (N*K)x(N*K) density matrix \rho = |\psi> <\psi| ;
C
C     - we obtain the desidered NxN mixed states \sigma by partial tracing the 
C	  pure state \rho over the ancillary space.
C	
C	NOTE: when K=(N+1) we get \sigma distributed according to the Hilbert
C	      Schmidt measure od REAL density matrices...
C
C	It takes as input N and K and gives as output real density matrix RMAT, of
C	size NxN, representing a random mixed state distributed according to the
C	induced measure P_(N,K). Of course RMAT is Hermitian and of trace 1. 
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 N, K, emme, enne, mu
      REAL U(2), RVEC(N,K), RMAT(N,N), 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 emme = 1, N
	    DO 200 mu = 1, K
   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(emme,mu) = DEVIAT

     	Do emme = 1, N

			RMAT(emme,emme) = 0.
			Do mu=1,K	
  	RMAT(emme,emme) = RMAT(emme,emme) + RVEC(emme,mu) * RVEC(emme,mu) 
			enddo
			RMAT(emme,emme) = RMAT(emme,emme) / SOMMA
	
		Do enne = emme+1, N

			RMAT(emme,enne) = 0.
			Do mu=1,K	
  	RMAT(emme,enne) = RMAT(emme,enne) + RVEC(emme,mu) * RVEC(enne,mu) 
			enddo
			RMAT(emme,enne) = RMAT(emme,enne) / SOMMA
			RMAT(enne,emme) = RMAT(emme,enne)

		enddo
	enddo

      RETURN
      END
