; +
; NAME:
;  CREATE_SIMPLEX
;
; PURPOSE:
;  This function returns the coordinates of the vertices of
;  an n-dimensional regular.  The algorithm is based on a
;  suggestion given by James Kuyper on the Google IDL-PVWAVE
;  NG.
;
; PARAMETERS (Required):
;  NDIM: integer dimensionality of the resulting simplex
;
; KEYWORDS: (Optional)
;  CENTER:  float or double vector of length NDIM specifying
;           the coordinates of the center of the NDIM-dimensional
;           hypersphere. (Default: [0.,0.,0.,....])
;  RADIUS:  float value specifying the radius of the
;           NDIM-dimensional hypersphere.
;  TEST:    If set then the length of each simplex and the
;           center of all of the simplices are returned.  This
;           is for diagnostic purposes.
;
; RETURN VALUE:
;  The coordinates of the (NDIM+1) simplex nodes.  Float or double
;  array with dimensions NDIM by NDIM+1.
;
; AUTHOR:
;  Robert Dimeo
;  NIST Center for Neutron Research
;  National Institute of Standards and Technology
;  100 Bureau Drive-Stop 8562
;  Gaithersburg, MD 20899
;
; CATEGORY:
;  OPTIMIZATION, AMOEBA, SIMPLEX
;
; CALLING SEQUENCE:
;  P = CREATE_SIMPLEX(NDIM, CENTER = center, RADIUS = radius, /TEST)
;
; EXAMPLE USAGE:
;  IDL> P = CREATE_SIMPLEX(3,CENTER = [0.0,0.5,3.2],RADIUS = 2.25, /TEST)
;
; REQUIREMENTS:
;  NONE
;
; REQUIRED PROGRAMS:
;  NONE
;
; SIDE EFFECTS:
;  NONE
;
; COMMON BLOCKS:
;  NONE
;
; MODIFICATION HISTORY:
;  RMD -- Wrote 10/31/04
; -
; ******************** ;
function create_simplex,ndim,center = center,radius = radius,test = test
compile_opt idl2,hidden
; Uses a technique suggested by James Kuyper to
; create a regular simplex in n-dimensions.  The return
; value is an array ndim by ndim+1.  The coordinates of
; the i-th vertex of the simplex are obtained by
; the following call:
; p =  create_simplex(5)
; vertex_3 = p[*,2]  ; The 3rd vertex in the simplex
;
p = dblarr(ndim,ndim+1)
if n_elements(center) eq 0 then center = dblarr(ndim)
if n_elements(radius) eq 0 then radius = 1d
; Fill in the values for n = 1
p[0,0] = -1D & p[0,1] = 1D

if ndim gt 1 then begin
   for j = 2,ndim do begin
      p[j-1,j] = sqrt(total((p[0:j-1,1]-p[0:j-1,0])^2)-total(p[0:j-1,0]^2))
      p[j-1,0:j] = p[j-1,0:j] - p[j-1,j]/(j+1)
   endfor
endif
; Normalize the hypersphere to unit radius
p = p/norm(p[0:ndim-1,0])
; Scale the simplex to the appropriate sphere radius
p = p*radius
; Shift the center
p = p+rebin(center,ndim,ndim+1,/sample)
if keyword_set(test) then begin
   print,'Center'
   for i = 0,ndim-1 do print,(moment(p[i,*]))[0]
   print,'Radius'
   pp = p
   for i = 0,ndim-1 do pp[i,*]=pp[i,*]-center[i]
   for i = 0,ndim do print,norm(pp[*,i])
endif
return,p
end