      SUBROUTINE C_mixed_Dmatr(N,K,CMAT,ROUTIN)

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 complex 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 COMPLEX 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 C_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 we get \sigma distributed according to the Hilbert
C	      Schmidt measure of COMPLEX density matrices...
C
C	It takes as input N and K and gives as output complex density matrix CMAT,
C	of size NxN, representing a random mixed state distributed according to
C	the induced measure P_(N,K). Of course CMAT 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
	INTEGER N, K, emme, enne, mu
      COMPLEX SW(0:1), CVEC(N,K), CMAT(N,N)
      REAL U(2), S, T, A, B, R1, R2, V, X, Y, Q, DEVIAT, SOMMA
      SAVE  S, T, A, B, R1, R2, SW
      EXTERNAL ROUTIN
      DATA  S, T, A, B / 0.449871, -0.386595, 0.19600, 0.25472/
      DATA  R1, R2 ,SW / 0.27597, 0.27846, (1. , 0.), (0. , 1.)/
C         generate pair of uniform deviates

	SOMMA = 0.
      DO 200 emme = 1, N
	    DO 200 mu = 2, 2*K+1
	IF(mu.LE.(K+1)) CVEC(emme,mu-1) = (0.,0.)
   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 CVEC(emme,mu/2) = CVEC(emme,mu/2) + DEVIAT * SW(MOD(mu,2))

     	Do emme = 1, N

			CMAT(emme,emme) = (0.,0.)
			Do mu=1,K	
  	CMAT(emme,emme)=CMAT(emme,emme)+CVEC(emme,mu)*CONJG(CVEC(emme,mu))
			enddo
			CMAT(emme,emme) = CMAT(emme,emme) / SOMMA
	
		Do enne = emme + 1, N

			CMAT(emme,enne) = (0.,0.)
			Do mu=1,K	
  	CMAT(emme,enne)=CMAT(emme,enne)+CVEC(emme,mu)*CONJG(CVEC(enne,mu))
			enddo
			CMAT(emme,enne) = CMAT(emme,enne) / SOMMA
			CMAT(enne,emme) = CONJG(CMAT(emme,enne))

		enddo
	enddo

      RETURN
      END
