;+
; NAME:
;       SWARM__DEFINE
;
; PURPOSE:
;
;	Object class definition for artificial self-driven particles that are
;	capable of swarming.  Dynamic implementation of this class can be
;	found in the widget program called XSWARM.PRO, appended to the end
;	of this class definition.  To run the widget program, compile this file
;	and type XSWARM at the IDL command line.
;
;	The motion is governed by a few rules as discussed in the following paper:
;	T.Vicsek et. al. Physical Review Letters, vol. 75, pp. 1226-1229 (1995).
;
;	The rules obeyed by the particles in this simulation are the following:
;
;	1.	Each particle moves with a constant speed.
;	2.	The nearest neighbors are defined as all of the other particles whose center
;		lies within a fixed distance (2 times the particle radius) of the particle
;		in question.
;	3.	At each time step a particle assumes the average heading of all of its
;		nearest neighbors.
;	4.	At each time step a random fluctuation in the heading of each particle
;		is applied.
;	5.	Periodic boundary conditions are imposed on the particles.
;
;	An order parameter defined as the magnitude of the average velocity of all
;	of the particles normalized to the individual speed is displayed showing
;	the onset of ordered motion for various values of the "NOISE" parameter.
;
; AUTHOR:
;
;       Robert M. Dimeo, Ph.D.
;	NIST Center for Neutron Research
;       100 Bureau Drive
;	Gaithersburg, MD 20899
;       Phone: (301) 975-8135
;       E-mail: robert.dimeo@nist.gov
;       http://www.ncnr.nist.gov/staff/dimeo
;
; CATEGORY:
;
;       Objects, physics
;
; CALLING SEQUENCE:
;
;       o = obj_new("SWARM",	VELOCITY = velocity,		$
;				POSITION = position,		$
;				RADIUS = radius)
;
;	VELOCITY:	two-element array specifying particles' velocity (float)
;	POSITION:	two-element array specifying particles' position (float)
;	RADIUS:		size of particles in pixels
;
;	To run the widget program XSWARM, compile this file and type
;
;	IDL> XSWARM
;
; PARAMETERS
;
;	None
;
; COMMON BLOCKS:
;
;	None
;
; DISCLAIMER
;
;	This software is provided as is without any warranty whatsoever.
;	Permission to use, copy, modify, and distribute modified or
;	unmodified copies is granted, provided this disclaimer
;	is included unchanged.
;
; MODIFICATION HISTORY:
;
;       Written by Rob Dimeo, December 29, 2002
;	01/03/02 (RMD):	Added option to display the order parameter defined
;			as the magnitude of the average total velocity.
;	01/04/02 (RMD):	Added option to display a histogram of the order
;			parameter so that an estimate of the order parameter
;			and its standard deviation can be obtained.
;	01/05/02 (RMD):	Added capability to run a sequence of "experiments"
;			for different values of NOISE and plot the resulting
;			mean of the distribution of order parameter.
;	01/07/02 (RMD):	Added continuously varying colors for the swarm objects
;			whose color index is based on their headings.  Also added
;			option to change color tables using XLOADCT.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro swarm::cleanup
ptr_free,self.cPtr
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro swarm::print
print,'Position:',self.position
print,'Velocity:',self.velocity
print,'Radius:',self.radius
print,'Index:',self.index
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro swarm::getProperty,	velocity = velocity, $
			position = position, $
			radius = radius, $
			index = index

index = self.index
radius = self.radius
velocity = self.velocity
position = self.position
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro swarm::draw
r = self.radius
numArg = self.velocity[1]
denArg = self.velocity[0]
tanArg = numArg/denArg
angle = atan(tanArg)
if denArg lt 0 then angle = angle + !pi
arg = 0.5*!pi-angle
xo = self.position[0]-self.radius*cos(angle)
x1 = self.position[0]+self.radius*cos(angle)
yo = self.position[1]-self.radius*sin(angle)
y1 = self.position[1]+self.radius*sin(angle)
nc = n_elements(*self.cPtr)
if angle lt 0. then cAngle = 2.*!pi+angle else cAngle = angle
index = fix(1.*nc*(cAngle)/(2.*!pi)-1)
color = (*self.cPtr)[index]
arrow, xo,yo,x1,y1,color = color,thick = 1.5
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro swarm::setProperty,	velocity = velocity, $
			position = position, $
			radius = radius, $
			index = index

if n_elements(index) ne 0 then self.index = index
if n_elements(velocity) ne 0 then self.velocity = velocity
if n_elements(position) ne 0 then self.position = position
if n_elements(radius) ne 0 then self.radius = radius
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function swarm::init,	velocity = velocity, $
			position = position, $
			radius = radius

if n_elements(radius) eq 0 then radius = 1.0
if n_elements(velocity) eq 0 then velocity = [1.0,1.0]
if n_elements(position) eq 0 then position = [0.0,0.0]
cLo = 16
cHi = 255
nc_one = cHi-cLo+1
nc_two = (cHi-1)-(cLo+1)+1
cVec = [cLo+indgen(nc_one),reverse(cLo+1+indgen(nc_two))]
self.cPtr = ptr_new(cVec)
self.velocity = velocity
self.position = position
self.radius = radius
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro swarm__define
swarm =	{	swarm,				$
		radius:0.0,			$
		index:0,			$
		color:0,			$
		cPtr:ptr_new(/allocate_heap),	$
		position:fltarr(2),		$
		velocity:fltarr(2)		$
	}

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 		END of swarm class definition
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;		BEGINNING of widget program called XSWARM
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_cleanup,tlb
widget_control,tlb,get_uvalue = pState
obj_destroy,(*pState).swarmContainer
wdelete,(*pState).winPix
ptr_free,(*pState).avPtr,(*pState).noisePtr
ptr_free,(*pState).meanPtr,(*pState).sigmaPtr
ptr_free,pState
loadct,0,/silent
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_quit,event
widget_control,event.top,/destroy
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_fireTimer,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).timerId,timer = (*pState).duration
(*pState).elapsed = (*pState).elapsed + 1
if (*pState).sequence eq 1 then begin
    if (*pState).elapsed gt (*pState).timeLimit then begin
      (*pState).loop = 0
	  xswarm_estimate,event,sigma = sigma,cen = cen
	  if (*pState).runIndex gt 1 then begin
	    *(*pState).meanPtr = [*(*pState).meanPtr,cen]
	    *(*pState).sigmaPtr = [*(*pState).sigmaPtr,sigma]
	  endif else begin
	    *(*pState).meanPtr = cen
	    *(*pState).sigmaPtr = sigma
	  endelse

      (*pState).runIndex = (*pState).runIndex + 1

      if (*pState).runIndex gt n_elements(*(*pState).noisePtr) then begin
        xswarm_plotSequence,event
        return
      endif

      noiseId = widget_info(event.top,find_by_uname = 'NOISE')
      widget_control,noiseId,set_value = (*(*pState).noisePtr)[(*pState).runIndex-1]
      xswarm_init,event
      xswarm_start_seq,event
    endif
endif
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_histo,event
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).avPtr) lt 2 then return
nbins = 50
blo = 0.0 & bhi = 1.0
db = (bhi-blo)/(nbins)
bins = blo+db*findgen(nbins)
nAv = n_elements(*(*pState).avPtr)
av = (*(*pState).avPtr)[nav/2:nav-1]
hist = histogram(av,min = blo,max = bhi,nbins = nbins)
newhist = [0.0,hist[0],hist,hist[nbins-1]]
newbins = [bins[0],bins[0],bins+0.5*db,bins[nbins-1]+db]
newhist = newhist/(db*nbins*total(newhist))
yrange = [0.0,1.5*max(newHist)]
plot,newbins,newhist,psym = 10,color = (*pState).black, $
	background = (*pState).white,xtitle = '!3<V>', $
	ytitle = '!3P(<V>)',yrange = yrange,/ysty, $
	title = '!3Order Parameter Probability Density'
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_velo,event
widget_control,event.top,get_uvalue = pState
oall = (*pState).swarmContainer->get(/all)
vxArray = fltarr((*pState).nParticles)
vyArray = fltarr((*pState).nParticles)
xArray = fltarr((*pState).nParticles)
for i = 0,(*pState).nparticles-1 do begin
  oall[i]->getProperty,velocity = velocity,position = position
  xArray[i] = position[0]
  vo = (*pState).speed
  vmag = norm(velocity)
  vxArray[i] = velocity[0]*vo/vmag
  vyArray[i] = velocity[1]*vo/vmag
endfor
dx = 0.25*!d.x_size
plot,[-vo,vo],[-vo,vo],/nodata,xrange = [-dx,!d.x_size+dx], $
  	yrange = 1.5*[-vo,vo],/xsty,/ysty,xtitle = '!3X', $
  	ytitle = '!3V!Dx!N',title = 'Phase space plot'
for i = 0,(*pState).nParticles-1 do begin
  plots,[xArray[i]],[vyArray[i]],psym = 8,/data
endfor
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_refresh,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).pGroup,get_value = val
case val[0] of
0:	begin	; Bubbles
		oall = (*pState).swarmContainer->get(/all)
		wset,(*pState).winPix
		erase
		for i = 0,(*pState).nparticles-1 do begin
  			oall[i]->draw
		endfor
	end
1:	begin	; Order parameter
		wset,(*pState).winPix
		num = 2
		if n_elements(*(*pState).avPtr) gt num then begin
			npts = n_elements(*(*pState).avPtr)
			plot,[*(*pState).avPtr],psym = 0,xtitle = '!3Time step', $
				ytitle = '!3<V>',charsize = 1.0,title = 'Order parameter', $
				nsum = num,yrange = [0.0,1.2],/ysty,background = (*pState).white, $
				color = (*pState).black,xrange = [0.0,1.2*npts],/xsty
			oplot,!x.crange,[1.0,1.0],linestyle = 2,color = (*pState).black,thick = 2.0
		endif else begin
			erase
		endelse
	end
2:	begin	; Histogram
		wset,(*pState).winPix
		xswarm_histo,event
	end
3:	begin	; Velocity space plot
		wset,(*pState).winPix
		xswarm_velo,event
	end
else:
endcase

wset,(*pState).winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_resize,event
widget_control,event.top,get_uvalue = pState
wdelete,(*pState).winPix
geom = widget_info((*pState).timerId,/geometry)
newxsize = event.x-geom.xsize
newysize = event.y
widget_control,(*pState).win,draw_xsize = newxsize,draw_ysize = newysize
widget_control,(*pState).win,get_value = winVis
(*pState).winVis = winVis
window,/free,/pixmap,xsize = newxsize,ysize = newysize
(*pState).winPix = !d.window
xswarm_refresh,event
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_init,event
widget_control,event.top,get_uvalue = pState
ptr_free,(*pState).avPtr
(*pState).avPtr = ptr_new(/allocate_heap)
(*pState).elapsed = 0
noiseId = widget_info(event.top,find_by_uname = 'NOISE')
widget_control,noiseId,get_value = noise
(*pState).pert = 0.01*float(noise[0])

numId = widget_info(event.top,find_by_uname = 'NUMSWARM')
widget_control,numId,get_value = numVal
(*pState).nparticles = fix(numVal[0])

heading = 2.0*!pi*randomu(s,(*pState).nparticles)
x = !d.x_size*randomu(s,(*pState).nparticles)
y = !d.y_size*randomu(s,(*pState).nparticles)

; If full of existing object then destroy and recreate the container.
obj_destroy,(*pState).swarmContainer
(*pState).swarmContainer = obj_new('IDL_CONTAINER')

id = widget_info(event.top,find_by_uname = 'RADIUS')
widget_control,id,get_value = val
(*pState).radius = float(val[0])

id = widget_info(event.top,find_by_uname = 'SPEED')
widget_control,id,get_value = val
(*pState).speed = float(val[0])

id = widget_info(event.top,find_by_uname = 'NUMSWARM')
widget_control,id,get_value = val
(*pState).nparticles = fix(val[0])

for i = 0,(*pState).nparticles-1 do begin
  ; Create a random velocity vector
  r = 2.0*!pi*randomu(s,1)
  vx = (*pState).speed*cos(r)
  vy = (*pState).speed*sin(r)
  radius = (*pState).radius
  o = obj_new('swarm',	position = [x[i],y[i]],	$
  			velocity = [vx,vy], 	$
  			radius = radius)
  o->setProperty,index = i
  (*pState).swarmContainer->add,o
endfor
xswarm_refresh,event
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_start_seq,event
widget_control,event.top,get_uvalue = pState
xswarm_init,event
(*pState).loop = 1
xswarm_fireTimer,event
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_start_manual,event
widget_control,event.top,get_uvalue = pState
(*pState).sequence = 0
xswarm_init,event
(*pState).loop = 1
xswarm_fireTimer,event
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_pause,event
widget_control,event.top,get_uvalue = pState
if (*pState).loop eq 1 then begin
  (*pState).loop = 0
endif else begin
  (*pState).loop = 1
  xswarm_fireTimer,event
endelse
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_update,event
; This is where we implement the intelligence!
widget_control,event.top,get_uvalue = pState
oall = (*pState).swarmContainer->get(/all)
xArray = fltarr((*pState).nparticles)
yArray = fltarr((*pState).nparticles)
rArray = fltarr((*pState).nparticles)
vxArray = fltarr((*pState).nparticles)
vyArray = fltarr((*pState).nparticles)

; First get all of the current information out...
for i = 0,(*pState).nparticles-1 do begin
  oall[i]->getProperty,position = position, $
  			radius = radius,velocity = velocity
  xArray[i] = position[0]
  yArray[i] = position[1]
  rArray[i] = radius
  vo = (*pState).speed
  vmag = norm(velocity)
  vxArray[i] = velocity[0]*vo/vmag
  vyArray[i] = velocity[1]*vo/vmag
endfor

;av = (moment(vxArray))[0]
av = norm([total(vxArray),total(vyArray)])/(((*pState).nParticles)*(*pState).speed)
sig = sqrt((moment(vxArray))[1])
if n_elements(*(*pState).avPtr) eq 0 then begin
	*(*pState).avPtr = av
endif else begin
	*(*pState).avPtr = [*(*pState).avPtr,av]
endelse

thRand = ((*pState).pert)*0.5*(-!pi+2.0*!pi*randomu(s,1,(*pState).nparticles))
for i = 0,(*pState).nparticles-1 do begin

  	distance = sqrt((xArray-xArray[i])^2+(yArray-yArray[i])^2)

	myNeighbors = where(distance le 2.0*rArray[i],countNeighbors)
    if countNeighbors gt 1 then begin
    ; Determine the average heading of all nearest neighbors
		notMe = where(myNeighbors ne i)
    	if countNeighbors eq 2 then begin
    		vxAve = vxArray[notMe]
    		vyAve = vyArray[notMe]
    	endif else begin
     		notMe = where(myNeighbors ne i)
    		vxAve = (moment(vxArray[notMe]))[0]
    		vyAve = (moment(vyArray[notMe]))[0]
    	endelse
    	vxArray[i] = vxAve
    	vyArray[i] = vyAve
    endif
    vxNew = vxArray[i]*cos(thRand[i])-vyArray[i]*sin(thRand[i])
    vyNew = vxArray[i]*sin(thRand[i])+vyArray[i]*cos(thRand[i])
	velocity = [vxNew,vyNew]

	; Is this swarm out of bounds?  If so then change the direction
	; of the velocity vector
  	if xArray[i] gt !d.x_size then xArray[i] = 0
  	if yArray[i] gt !d.y_size then yArray[i] = 0
  	if xArray[i] lt 0 then xArray[i] = !d.x_size
	if yArray[i] lt 0 then yArray[i] = !d.y_size

	; Now finally update the positions
  	position = [xArray[i],yArray[i]]
  	position = position + velocity

  	oall[i]->setProperty,velocity = velocity,position = position

endfor
; Redisplay the swarm with a new position
xswarm_refresh,event
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro gaussfunct, x, a, f
area = 1.d
cen = a[0]
sigma = a[1]
arg = 0.5*((x-cen)/sigma)^2
f = (area/sqrt(2.0*!dpi*sigma^2))*exp(-arg)
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_estimate,event,sigma = sigma,cen = cen
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).pGroup,get_value = val
if (*pState).loop eq 0 then begin
  if n_elements(*(*pState).avPtr) lt 2 then return

  nbins = 50
  blo = 0.0 & bhi = 1.0
  db = (bhi-blo)/(nbins)
  bins = blo+db*findgen(nbins)
  nAv = n_elements(*(*pState).avPtr)
  av = (*(*pState).avPtr)[nav/2:nav-1]
  hist = histogram(av,min = blo,max = bhi,nbins = nbins)
  newhist = [hist[0],hist,hist[nbins-1]]
  newbins = [bins[0],bins+0.5*db,bins[nbins-1]+db]
  err = sqrt(newbins)
  wherezero = where(newbins eq 0.0,countZero)
  if countZero gt 0 then err[whereZero] = 1.0
  factor = (db*total(newhist))
  newhist = newhist/factor
  newerr = err/factor
  yrange = [0.0,1.5*max(newHist)]
  area = 1.0 & cen = 0.5 & fwhm = 0.5
  sigma = fwhm/2.354
  parms = [cen,sigma]
  w = 1.d/newerr^2
  yfit = CURVEFIT(newbins, newhist, w, $
		parms, function_name='gaussfunct', $
		/noderivative,itmax = 50,chisq = chisq)
  sigma = parms[1]
  cen = parms[0]
  if val[0] eq 2 then begin
    wset,(*pState).winPix
    plot,newbins,newhist,psym = 10,color = (*pState).black, $
	  background = (*pState).white,xtitle = '!3<V>', $
	  ytitle = '!3P(<V>)',yrange = yrange,/ysty, $
	  title = '!3Order Parameter Probability Density'
    errplot,newbins,newhist-newerr,newhist+newerr,width = 0.0
    oplot,newbins,yfit,psym = 0,thick = 2.0,color = (*pState).black
    s1 = 'Mean='+strtrim(string(cen),2)
    s2 = 'Sigma='+strtrim(string(sigma),2)
    s3 = 'Chi-squared='+strtrim(string(chisq),2)
    xyouts,0.2,0.9,s1,/normal,color = (*pState).black
    xyouts,0.2,0.85,s2,/normal,color = (*pState).black
    xyouts,0.2,0.8,s3,/normal,color = (*pState).black
    wset,(*pState).winVis
    device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
  endif
endif

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro gsi_event,event
widget_control,event.top,get_uvalue = pState
uname = widget_info(event.id,/uname)
if uname eq 'EXECUTE' then (*pState).cancel = 0
if uname eq 'CANCEL' then (*pState).cancel = 1
widget_control,(*pState).loField,get_value = lo
(*pState).lo = float(lo[0])
widget_control,(*pState).hiField,get_value = hi
(*pState).hi = float(hi[0])
widget_control,(*pState).numField,get_value = num
(*pState).num = fix(num[0])
widget_control,(*pState).timeField,get_value = time
(*pState).time = fix(time[0])
widget_control,event.top,/destroy
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function getSeqInfo,group_leader = group_leader,preferences = preferences
if n_elements(preferences) ne 0 then begin
  lo = preferences.lo
  hi = preferences.hi
  num = preferences.num
  time = preferences.time
endif else begin
  lo = 0.1
  hi = 5.0
  num = 5
  value = 100
endelse
tlb = widget_base(group_leader = group_leader,/col, $
	title = 'Sequence Information',/modal,/base_align_right)
loField = cw_field(tlb,title = 'Lower noise value', $
	value = strtrim(string(lo),2),/row)
hiField = cw_field(tlb,title = 'Upper noise value', $
	value = strtrim(string(hi),2),/row)
numField = cw_field(tlb,title = '# noise values', $
	value = strtrim(string(num),2),/row)
timeField = cw_field(tlb,title = '# time steps', $
	value = strtrim(string(time),2),/row)
void = widget_button(tlb,value = 'CANCEL',uname = 'CANCEL')
void = widget_button(tlb,value = 'EXECUTE SEQUENCE',uname = 'EXECUTE')
widget_control,tlb,/realize
state = {	loField:loField,	$
			hiField:hiField,	$
			numField:numField,	$
			timeField:timeField,	$
			cancel:0,			$
			lo:0.1, hi:5.0, num:5,time:100}
pState = ptr_new(state)
widget_control,tlb,set_uvalue = pState
xmanager,'gsi',tlb,/no_block
lo = (*pState).lo & hi = (*pState).hi
num = (*pState).num & time = (*pState).time
cancel = (*pState).cancel
ptr_free,pState
return,{lo:lo,hi:hi,num:num,time:time,cancel:cancel}
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_seqSetup,event
; Set up sequence
widget_control,event.top,get_uvalue = pState
result = getSeqInfo(group_leader = event.top,preferences = (*pState).preferences)
(*pState).preferences = result
if result.cancel eq 1 then return
(*pState).sequence = 1
d = (result.hi-result.lo)/(result.num-1.)
eta = result.lo+d*findgen(result.num)
*(*pState).noisePtr = eta
(*pState).timeLimit = result.time
xswarm_seqRun,event
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_seqRun,event
; Run sequence
widget_control,event.top,get_uvalue = pState
nruns = n_elements(*(*pState).noisePtr)
if nruns eq 0 then return
ptr_free,(*pState).sigmaPtr,(*pState).meanPtr
(*pState).sigmaPtr = ptr_new(/allocate_heap)
(*pState).meanPtr = ptr_new(/allocate_heap)
noiseId = widget_info(event.top,find_by_uname = 'NOISE')
(*pState).runIndex = 1
widget_control,noiseId,set_value = (*(*pState).noisePtr)[(*pState).runIndex-1]
xswarm_init,event
xswarm_start_seq,event
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_plotSequence,event
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).noisePtr) eq 0 then return
x = *(*pState).noisePtr
dx = 0.2*(max(x)-min(x))
xrange = [min(x)-dx,max(x)+dx]
y = *(*pState).meanPtr
yerr = abs(*(*pState).sigmaPtr)

dy = 0.2*(max(y+yerr))
yrange = [0.0,max(y+yerr)+dy]
wset,(*pState).winPix
plot,x,y,psym = 4,color = (*pState).white,background = (*pState).black, $
	xrange = xrange,/xsty,xtitle = '!3Noise (%)',ytitle = 'Order Parameter', $
	yrange = yrange,/ysty
;plot,x,y,psym = 4,color = (*pState).white,background = (*pState).black, $
;	xtitle = '!3Noise (%)',ytitle = 'Order Parameter', $
;	/xlog,/ylog
errplot,x,y-yerr,y+yerr,width = 0.0
wset,(*pState).winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_refreshcb,data = data
xswarm_refresh,data
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm_handler,event
; Main event dispatcher
widget_control,event.top,get_uvalue = pState
if widget_info(event.id,/uname) eq 'RADIUS' then begin
  widget_control,(event.id),get_value = val
  (*pState).radius = float(val[0])
  oall = (*pState).swarmContainer->get(/all)
  for i = 0,n_elements(oall)-1 do begin
    oall[i]->setProperty,radius = (*pState).radius
  endfor
endif

if widget_info(event.id,/uname) eq 'SPEED' then begin
  widget_control,(event.id),get_value = val
  (*pState).speed = float(val[0])
endif

if widget_info(event.id,/uname) eq 'NOISE' then begin
  widget_control,(event.id),get_value = val
  (*pState).pert = 0.01*float(val[0])
endif

if widget_info(event.id,/uname) eq 'NUMSWARM' then begin
  (*pState).loop = 0
  widget_control,(event.id),get_value = val
  (*pState).nparticles = fix(val[0])
  obj_destroy,(*pState).swarmContainer
  (*pState).swarmContainer = obj_new('IDL_CONTAINER')
  xswarm_init,event
endif

if widget_info(event.id,/uname) eq 'PLOTTYPE' then begin
  xswarm_refresh,event
endif

if tag_names(event,/structure_name) eq 'WIDGET_BASE' then begin
  xswarm_resize,event
  xswarm_fireTimer,event
endif
if widget_info(event.id,/uname) eq 'CHANGE COLORS' then begin
	xloadct,group = event.top,updatecallback = 'xswarm_refreshcb', $
		updatecbdata = event,/silent,/modal
endif

if (*pState).loop eq 0 then return
if tag_names(event,/structure_name) eq 'WIDGET_TIMER' then begin
  xswarm_update,event
  xswarm_fireTimer,event
endif
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro xswarm
; widget definition module
!except = 0
device,decomposed = 0
colorTable = 2
loadct,colorTable,/silent
tlb = widget_base(/row,title = 'xswarm',/tlb_size_events,mbar = bar)
filemenu = widget_button(bar,value = 'FILE',/menu)
void = widget_button(filemenu,value = 'SET UP SEQUENCE', $
		event_pro = 'xswarm_seqSetup')
void = widget_button(filemenu,value = 'CHANGE COLOR TABLE', $
		uname = 'CHANGE COLORS')
void = widget_button(filemenu,value = 'QUIT',event_pro = 'xswarm_quit')

colBase = widget_base(tlb,/col)
void = cw_field(colBase,/col,title = '# swarm',value = '60', $
		uname = 'NUMSWARM',/return_events)
void = cw_field(colBase,/col,title = 'Speed',value = '6.0', $
		uname = 'SPEED',/return_events)
void = cw_field(colBase,/col,title = 'swarm size',value = '8.0', $
		uname = 'RADIUS',/return_events)
void = cw_field(colBase,/col,title = 'NOISE (% of 360 degrees)',value = '15.0', $
		uname = 'NOISE',/return_events)
void = widget_button(colBase,value = 'START',event_pro = 'xswarm_start_manual', $
		uname = 'START')
void = widget_button(colBase,value = 'PAUSE/RESUME',event_pro = 'xswarm_pause')
void = widget_button(colBase,value = 'ESTIMATE <V>',event_pro = 'xswarm_estimate')

ptypes = ['Swarm','Order parameter','Histogram','Phase space']
pgroup = cw_bgroup(colBase, ptypes, /col, /exclusive,$
         label_top='Plot Type',/no_release,$
         frame=1,set_value=0,/return_index, $
         uname = 'PLOTTYPE')

xsize = 400 & ysize = 10
win = widget_draw(tlb,xsize = xsize,ysize = ysize)

widget_control,tlb,/realize
geom = widget_info(colBase,/geometry)

widget_control,win,draw_ysize = geom.ysize,get_value = winVis
window,/free,/pixmap,xsize = xsize,ysize = geom.ysize
winPix = !d.window
preferences = { lo:0.1, hi:5.0, num:5,time:100 ,cancel:0}
state = {	winVis:winVis,			$
		winPix:winPix,			$
		win:win,				$
		pGroup:pGroup,			$
		preferences:preferences,	$
		speed:8.0,				$
		green:28,				$
		red:73,					$
		black:0,				$
		white:255,				$
		avPtr:ptr_new(/allocate_heap),		$
		noisePtr:ptr_new(/allocate_heap),	$
		runIndex:0,				$
		timeLimit:1,			$
		sequence:0,				$
		elapsed:0,				$
		loop:0,					$
		nParticles:60,			$
		pert:0.03,				$
		radius:1.0,				$
		duration:0.0001,		$
		timerID:colBase,		$
		meanPtr:ptr_new(/allocate_heap),	$
		sigmaPtr:ptr_new(/allocate_heap),	$
		swarmContainer:obj_new('IDL_CONTAINER')}

nth = 20
th = 2.0*!pi*findgen(nth)/(nth-1.)
xc = cos(th) & yc = sin(th)
usersym,xc,yc,/fill

pState = ptr_new(state,/no_copy)
widget_control,tlb,set_uvalue = pState
ev = {pEvent,id:void,top:tlb,handler:0L}
xswarm_init,ev
xmanager,'xswarm',tlb,/no_block,cleanup = 'xswarm_cleanup', $
	event_handler = 'xswarm_handler'
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;