      SUBROUTINE R_pure_vect(RVEC,VLEN,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 components for a REAL random pure state.
C   	It takes as input:
C
C	+ RVEC :: a real vector of length VLEN
C	+ VLEN :: length of the real vector RVEC
C
C     and gives as output the same vector RVEC, whose entries now 
C     consist of VLEN (rescaled) independent Gaussian random R-numbers.
C     Rescaled means that the L_2-norm of RVEC is 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 VLEN, IDEV
      REAL U(2), RVEC(VLEN), 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, VLEN
   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

	SOMMA = SQRT(SOMMA)
      DO 300 IDEV = 1, VLEN
  300	RVEC(IDEV) = RVEC(IDEV) / SOMMA
      RETURN
      END
