FUNCTION RANDOMQ,Q,N,SEED=seed,MEAN=mu
;+
;z=RandomQ(Q,N) returns an array z of normal random vectors
;with N columns which have the covariances defined by Q.
;
;OPTIONAL PARAMETERS
;
;SEED is a keyword specifiying the seed to the random number generator.
;MEAN is a keyword specifying the desired mean vector. Default is 0.
;
;NOTES:
;Q must be a positive definite symmetric matrix.
;Uses RANDOMN and CHOLDC.
;
;Harvey Rhody
;November, 2002.
;-
sq=size(Q)
IF sq[0] NE 2 OR sq[1] NE sq[2] THEN $
	MESSAGE,'USAGE: Z=RANDOMQ(Q,N) where Q is 2D symmetric positive definite array.'
m=sq[1]

;Factor Q into Q=L##Transpose(L) using CHOLDC.
A=Q
CHOLDC,A,P

;Extract the lower triangular part of A, which contains the
;part of L below the diagonal.
k=LINDGEN(m,m)
kl=WHERE(k/m GT (k Mod m))
L=Fltarr(m,m)
L[kl]=A[kl]
;Put in the diagonal values from the vector P.
L[(m+1)*Lindgen(m)]=P

;Generate n random vectors and multiply by L
x=L##Randomn(seed,n,m)

;If a mean vector is provided then shift the means.
IF N_ELEMENTS(mu) EQ m THEN $
	x=x+replicate(1,n)#(mu[*])

RETURN,X
END