pro trivariate_gaussian_f,x,y,z,n,rhoxy,rhoxz,rhoyz,sx=sx,sy=sy,sz=sz,$
                          plot=plot,seed0=seed0

if(n_params() le 0) then begin
    print,'--------------------------------------------------------'
    print,'trivariate_gaussian_f,x,y,z,N,rhoxy,rhoxz,rhoyz,sx=sx,sy=sy,sz=sz'
    print,'create N triples x,y,z following a trivariate Gaussian pdf'
    print,'with zero mean and unit variances'
    print,'keywords:'
    print,'  sx=sigmax, sy=sigmay, sz=sigmaz (defs=1)'
    print,'  plot=1,2,3     1-xy, 2-xz   3-yz
    print,'EXAMPLE:'
    print,'trivariate_gaussian_f,x,y,10000,0.5,0.5,0.25,plot=1'
    print,'--------------------------------------------------------'
    return
endif

;------------------------------------------------------------------------
;Trivariate distribution with correlation coefficient rho
;i) Set up the covariance matrix: cov[x,y]=rho * sigma_x * sigma_y
; sigma_x^2      cov[x,y]   cov[x,z]
; cov[y,x]       sigma_y^2  cov[y,z]
; cov[z,x]       cov[z,y]   sigma_z^2  

;ii) Find the eigenvalues and eigenvectors of covariance
;matrix. Combine eigenvectors (column vectors) to a matrix T

;iii) Create independent variable x',y',z' using variances
;equal to the three eigenvalues

;iv) Finally, (x,y,z) triples are obtained with
;    transfomation (x,y,z) =  T (x', y', z')
;where T is formed from the eigenvectors
;------------------------------------------------------------------------

;defaults:
  sigmax=1.
  sigmay=1.
  sigmaz=1.

  if(keyword_set(sx)) then sigmax=sx
  if(keyword_set(sy)) then sigmay=sy
  if(keyword_set(sz)) then sigmaz=sz

  sigmax2=sigmax^2
  sigmay2=sigmay^2
  sigmaz2=sigmaz^2

  cov_xy=rhoxy*sigmax*sigmay
  cov_xz=rhoxz*sigmax*sigmaz
  cov_yz=rhoyz*sigmay*sigmaz
  cov_yx=cov_xy
  cov_zx=cov_xz
  cov_zy=cov_yz


  cc=[[sigmax2, cov_xy,  cov_xz ],$
      [cov_yx,  sigmay2, cov_yz ],$
      [cov_zx,  cov_zy,  sigmaz2]]
 

;obtain eigenvalues & eigenvectors  j=0,1,2
;Using IDL function EIGENQL (for real symmetric matrix)

  lambda=eigenql(cc,eigenvectors=ee)

;print,'cc'
;print,cc

;print,'lambda'
;print,lambda

;print,'ee'
;print,ee

;calculate three independent variables xp,yp,zp

if(keyword_set(seed)) then seed=seed0
  xp=randomn(seed,n)*sqrt(lambda(0))
  yp=randomn(seed,n)*sqrt(lambda(1)) 
  zp=randomn(seed,n)*sqrt(lambda(2)) 


;transform these to x,y,z triples
;    e1         e2         e3
  x=ee(0,0)*xp+ee(0,1)*yp+ee(0,2)*zp
  y=ee(1,0)*xp+ee(1,1)*yp+ee(1,2)*zp
  z=ee(2,0)*xp+ee(2,1)*yp+ee(2,2)*zp

;--------------------------------------------------------------------
;
  if(keyword_set(plot)) then begin



;--------------------------------
;projections to xy, xz, yz planes

      if(n_elements(plot) eq 1) then begin
          if(plot(0) eq 1) then begin
              xp=x
              yp=y
              xtitle='x'
              ytitle='y'
              title='rho_xy='+string(rhoxy,'(f8.4)')
          endif         
          if(plot(0) eq 2) then begin
              xp=x
              yp=z
              xtitle='x'
              ytitle='z'
              title='rho_xz='+string(rhoxz,'(f8.4)')
          endif
          if(plot(0) eq 3) then begin
              xp=y
              yp=z
              xtitle='y'
              ytitle='z'
              title='rho_yz='+string(rhoyz,'(f8.4)')
          endif

          widx=5
          widy=5
          r=correlate(xp,yp)
          title=title+'   r='+string(r,'(f8.4)')
          nwin
          plot,xp,yp,xr=[-1.,1.]*widx,xs=1,yr=[-1.,1.]*widy,$
            ys=1,/iso,psym=3,title=title,xtitle=xtitle,ytitle=ytitle
      endif


;--------------------------------
;oblique projection
      
      if(n_elements(plot) gt 1) then begin

;turn with respect to x axis by ax degree (def=30)
;turn with respect to z axis by az degree (def=30)

          if(plot(0) ne 1) then  ax=plot(0)
          if(plot(1) ne 1) then  az=plot(1)
          
;just create the transformation (/save) with surface, do not plot
          xt=findgen(11)-5.
          shade_surf,xt#xt,xt,xt,/save,xr=[-5,5],yr=[-5,5],zr=[-5,5],$
            xs=1,ys=1,zs=1,/nod,chars=2,xtitle='x',ytitle='y',ztitle='z',$
            ax=ax,az=az
          
;3-d plotting with the created transformation (/t3d)
          plots,x,y,z,/t3d,psym=3
                    
;for the plotting of 3-sigma ellipsoid
          phi=findgen(201)/200.*2*!pi
          sinp=sin(phi)
          cosp=cos(phi)
        

          e1=ee(*,0)
          e2=ee(*,1)
          e3=ee(*,2)
          l1=sqrt(lambda(0))*3.
          l2=sqrt(lambda(1))*3.
          l3=sqrt(lambda(2))*3.

;plot largest eigenvector + 3-sigma circle perpendicular to it
          print,'e1:',e1,   'lambda1:',lambda(0)
          plots,[0,e1(0)*l1],[0,e1(1)*l1],[0,e1(2)*l1],/t3d,$
            col=2,lines=0,thick=5
          xt=l2*cosp*e2(0)+l3*sinp*e3(0)
          yt=l2*cosp*e2(1)+l3*sinp*e3(1)
          zt=l2*cosp*e2(2)+l3*sinp*e3(2)
          plots,xt,yt,zt,/t3d,lines=0,col=2,thick=5
          
;plot second largest 
          print,'e2:',e2,   'lambda2:',lambda(1)
          plots,[0,e2(0)*l2],[0,e2(1)*l2],[0,e2(2)*l2],/t3d,$
            col=3,thick=2,lines=0
          xt=l1*cosp*e1(0)+l3*sinp*e3(0)
          yt=l1*cosp*e1(1)+l3*sinp*e3(1)
          zt=l1*cosp*e1(2)+l3*sinp*e3(2)
          plots,xt,yt,zt,/t3d,thick=2,lines=0,col=3
                    
;and third (the smallest)
          print,'e3:',e3,   'lambda3:',lambda(2)
          plots,[0,e3(0)*l3],[0,e3(1)*l3],[0,e3(2)*l3],/t3d,$
            col=5,thick=2,lines=2
          xt=l1*cosp*e1(0)+l2*sinp*e2(0)
          yt=l1*cosp*e1(1)+l2*sinp*e2(1)
          zt=l1*cosp*e1(2)+l2*sinp*e2(2)
          plots,xt,yt,zt,/t3d,col=5,thick=2,lines=2
          
      endif

  endif      


end