#pragma rtGlobals=1		// Use modern global access method.

// Date: 11/21/2012 //Version 1.2
// This procedure creates a simulation of a particle experiencing a series of random molecular collisions (each of which steps the molecule a mean free path length).
// After a defined number of collisions, the particle is at a new, random location. The movement between the starting locations and final location is termed a "displacement".
//In the laboratory setting, this is equivalent to the displacement a particle undergoes in between two consecutive frames when imaged with a CCD camera.
// A series of random displacements can then be used to estimate Brownian Motion for a particle in 2-D space.
//
// Written by Josef G. Trapani (jtrapani@amherst.edu) and Ashley R. Carter (acarter@amherst.edu)

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////

Menu "Macros"

	"Simulate Brownian Motion", Start_BM_Panel("start")
	
End

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
Function Start_BM_Panel(ctrlName)
	String ctrlName
	String savedDF= GetDataFolder(1)	// Remember CDF in a string.

	DoWindow BM_Panel0
	If (V_flag!=0)
		KillWindow BM_Panel0
		SetDataFolder root:BM
		DoWindow /K BM_histgraph1
		DoWindow /K BM_histgraph2
		DoWindow /K BM_graph1
		DoWindow /K BM_graph2
		DoWindow /K BM_graph3
		DoWindow /K D_table1
		KillWaves /Z /A
		SetDataFolder savedDF			// and restore
		KillDataFolder /Z  root:BM
	Endif

	If (cmpstr(ctrlName,"start")==0)
		NewDataFolder/O/S root:BM
		Variable /G collisionNum = 99
		Variable /G displacementNum = 4500
		Variable /G trackNum = 1
		Variable /G pathNum = 3	//based on ~3 angstrom = the mean free path length for a bead in water, before it collides with another molecule
		SetDataFolder savedDF			// and restore
		BM_Panel(ctrlName)
	Endif

End

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
Function BM_Panel(ctrlName)
	String ctrlName
	
	String info = IgorInfo(0)
	String screen1RectStr = StringByKey("SCREEN1", info)		//e.g., "DEPTH=23,RECT=0,0,1280,1024"
	Variable depth, left, top, right, bottom
	sscanf screen1RectStr, "DEPTH=%d,RECT=%d,%d,%d,%d", depth, left,top, right, bottom
	
	String platform= UpperStr(igorinfo(2))
	If (cmpstr(platform,"WINDOWS")==0)
		NewPanel /N=BM_Panel0 /W=(right-300,top,right,200) /K=2
	Else	
		NewPanel /N=BM_Panel0 /W=(right-400,top,right-100,200) /K=2
	Endif
	ModifyPanel /W=BM_Panel0 cbRGB=(62535,62535,62535), frameStyle=3, frameInset=5, noEdit=1, fixedSize=1


	SetDrawLayer /W=BM_Panel0 ProgFront
	DrawRRect /W=BM_Panel0 10,35,290,120
	SetDrawEnv /W=BM_Panel0 fsize= 16, fillbgc=(49151,53155,65535)
	DrawText /W=BM_Panel0 30,65,"2-D Brownian Motion Simulator"
	SetDrawEnv /W=BM_Panel0 fsize= 10, fillbgc=(49151,53155,65535)

	Button closeBMPanelbutton,pos={210,10},size={60,20},proc=Start_BM_Panel,title="Close", win=BM_Panel0

	SetVariable BMsetvar0,pos={20,160},size={100,15},title="Collisions:", win=BM_Panel0
	SetVariable BMsetvar0,fSize=10,limits={1,100000,0},value= root:BM:collisionNum, win=BM_Panel0

	SetVariable BMsetvar1,pos={130,160},size={140,15},title="Displacements:", win=BM_Panel0
	SetVariable BMsetvar1,fSize=10,limits={1,10000,0},value= root:BM:displacementNum, win=BM_Panel0
	
	SetVariable BMsetvar2,pos={200,130},size={80,15},title="Tracks:", win=BM_Panel0
	SetVariable BMsetvar2,fSize=10,limits={1,100,0},value= root:BM:trackNum, win=BM_Panel0

	SetVariable BMsetvar3,pos={20,130},size={140,15},title="Mean Free Path:", win=BM_Panel0
	SetVariable BMsetvar3,fSize=10,limits={0,100,0},value= root:BM:pathNum, win=BM_Panel0


	Button BMbutton0, pos={90,70}, size={120,40}, title="Take a\rrandom walk", proc=runBM, fColor=(65535,54607,32768), win=BM_Panel0
	
End

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
Function runBM(ctrlName):ButtonControl
	String ctrlName
	NVAR numtracks= root:BM:trackNum

	If (numtracks==1)
		runBM1(numtracks)
	Else
		runBM2(numtracks)
	Endif
End
	
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
Function runBM1(tempNum)

	Variable tempNum
	Variable a,z,f,j
	NVAR numCollisions= root:BM:collisionNum
	NVAR numDisplacements= root:BM:displacementNum
	NVAR numtracks= root:BM:trackNum
	NVAR numpath= root:BM:pathNum
	Variable elapsetime, ticktime, i
	String savedDF= GetDataFolder(1)	// Remember CDF in a string.


	String BM_distances = "root:BM:BM_distances"
	String BM_XWave = "root:BM:BM_X"
	String BM_YWave = "root:BM:BM_Y"
	String BM_XYwave = "root:BM:BM_XYvals"
	String BMhistwave = "root:BM:BM_stdHist"
	String histwave = "root:BM:BM_distancesHist"
	
	//Make waves to store tally of displacements as consecutive X,Y points
	String runningBM_XWave = "root:BM:BM_X_tally"
	String runningBM_YWave = "root:BM:BM_Y_tally"
	Make /O /N=(numDisplacements) $runningBM_XWave
	Wave runningBMX = $runningBM_XWave
	Make /O /N=(numDisplacements) $runningBM_YWave
	Wave runningBMY = $runningBM_YWave
	
	//Make a wave to store each displacement
	Make /O /N=0 $BM_distances
	Wave BM_distancesWave = $BM_distances

	Make /O /N=(numDisplacements) $BM_XWave
	Wave BMX = $BM_XWave
	Make /O /N=(numDisplacements) $BM_YWave
	Wave BMY = $BM_YWave

	String X_startsum = "root:BM:BM_startXvalues"	
	String X_endsum = "root:BM:BM_endXvalues"	
	Make /O /N=(numDisplacements) $X_startsum
	Wave xstartsum =  $X_startsum
	Make /O /N=(numDisplacements) $X_endsum
	Wave xendsum =  $X_endsum

	String Y_startsum = "root:BM:BM_startYvalues"
	String Y_endsum = "root:BM:BM_endYvalues"
	Make /O /N=(numDisplacements) $y_startsum
	Wave ystartsum =  $y_startsum
	Make /O /N=(numDisplacements) $Y_endsum
	Wave yendsum =  $y_endsum

	//Kill all the graphs
	DoWindow /K BM_graph1
	DoWindow /K BM_graph2
	DoWindow /K BM_graph3
	DoWindow /K BM_histgraph1
	DoWindow /K BM_histgraph2

	
	For (f=0; f<numDisplacements; f+=1)

		String X_steps = "root:BM:BM_frameXvals"
		String X_start = "root:BM:BM_startXval"
		String X_end = "root:BM:BM_endXval"
		String Y_steps = "root:BM:BM_frameYvals"
		String Y_start = "root:BM:BM_startYval"
		String Y_end = "root:BM:BM_endYval"


		//For making multiple graphs of random walks
		//**in the future make a check box and warning that this is time intensive for many frames!
		//	String X_steps = "root:BM:BM_frameXvals"+"_"+num2str(f)
		//	String X_start = "root:BM:BM_startXval"+"_"+num2str(f)
		//	String X_end = "root:BM:BM_endXval"+"_"+num2str(f)
		//	String Y_steps = "root:BM:BM_frameYvals"+"_"+num2str(f)
		//	String Y_start = "root:BM:BM_startYval"+"_"+num2str(f)
		//	String Y_end = "root:BM:BM_endYval"+"_"+num2str(f)



		Make /O /N=(numCollisions) $(X_steps)
		Make /O /N=(numCollisions) $(Y_steps) 
		Wave framexsteps =  $(X_steps) 
		Wave frameysteps =  $(Y_steps) 

		Make /O /N=(1) $(X_start)
		Wave xstart =  $(X_start)
		
		Make /O /N=(1) $(X_end)
		Wave xend =  $(X_end)
		
		Make /O /N=(1) $(Y_start)
		Wave ystart =  $(y_start)

		Make /O /N=(1) $(Y_end)
		Wave yend =  $(y_end)

		//Initialize the waves and variables for each displacement
		framexsteps=0
		frameysteps=0
		Variable horiz=0
		Variable vert=0
	
		For (a=0; a<numCollisions; a+=1)

			//Generate random number for length of free path (around the mean)
			j=sqrt((gnoise(numPath))^2)

			//Generate random number for a horizontal, vertical, or diagonal collision
			z=sqrt((8*enoise(1))^2)
			z=Trunc(z)
 
			If (z==0)
				vert=vert+j
			Elseif (z==1)
				vert=vert-j
			Elseif (z==2)
				vert=vert+j/(sqrt(2))  //Makes a diagonal "collision" = j units
				horiz=horiz+j/(sqrt(2)) //Makes a diagonal "collision" = j units
			ElseIf (z==3)
				horiz=horiz+j
			Elseif (z==4)
				horiz=horiz-j
			Elseif (z==5)
				vert=vert-j/(sqrt(2)) //Makes a diagonal "collision" = j unit
				horiz=horiz+j/(sqrt(2)) //Makes a diagonal "collision" = 1 unit
			Elseif (z==6)
				vert=vert+j/(sqrt(2)) //Makes a diagonal "collision" = j unit
				horiz=horiz-j/(sqrt(2)) //Makes a diagonal "collision" = 1 unit
			Elseif (z==7)
				vert=vert-j/(sqrt(2)) //Makes a diagonal "collision" = j unit
				horiz=horiz-j/(sqrt(2)) //Makes a diagonal "collision" = j unit
			Endif


			framexsteps[a]=horiz
			frameysteps[a]=vert
			
			//Calculate start position
			If (a==0)
				Variable startXpos = framexsteps[a]
				xstartsum[f]=startXpos //accumulate for each frame
				Variable startYpos = frameysteps[a]
				ystartsum[f]=startypos //accumalte for each frame
			Endif

			//Calculate end X,Y position after all collisions
			If (a==numCollisions-1)
				Variable endXpos = framexsteps[a]
				xendsum[f]=endXpos
		
				Variable endYpos = frameysteps[a]
				yendsum[f]=endYpos

			Endif
			if(f==(numDisplacements-1))
				xstart[0]=startXpos //store in a wave for plotting
				xend[0]=endXpos
				ystart[0]=startypos //store in a wave for plotting
				yend[0]=endYpos
			endif

		Endfor //end For statement for collisions

		//calculate distance from start to end position after all collisions
		Variable BMlengthX = (endXpos-startXpos)
		BMX[f]=BMlengthX
		Variable BMlengthY = (endYpos-startYpos)
		BMY[f]=BMlengthY

		//Running tally of X and Y for Track display
		If (f==0)
			runningBMX[f] = BMX[f]
			runningBMY[f] = BMY[f]
		Else
			runningBMX[f] = runningBMX[f-1]+BMX[f]
			runningBMY[f] = runningBMY[f-1]+BMY[f]
		Endif		

		Variable BMlength = sqrt((endXpos-startXpos)^2+(endYpos-startYpos)^2)
		Variable d= DimSize($BM_distances,0)
		InsertPoints d,1, $BM_distances
		BM_distancesWave[d+1]=BMlength

		//		//Graph display inside <For/Endfor> for plotting each frame in a separate graph
		//		Display /N=BM_graph1 /K=1 /W=(200,50,600,450 ) $Y_steps vs $X_steps
		//		ModifyGraph rgb=(0,0,0)
		//		ModifyGraph mirror=1
		//		SetAxis/A/N=1/E=2 left
		//		SetAxis/A/N=1/E=2 bottom
		////		
		////		AppendToGraph $Y_start vs $X_start
		////		ModifyGraph mode($NameOfWave(Ystart))=3,marker($NameOfWave(Ystart))=19, msize($NameOfWave(Ystart))=8
		////		ModifyGraph useMrkStrokeRGB($NameOfWave(Ystart))=1
		////		ModifyGraph rgb($NameOfWave(Ystart))=(2,39321,1)
		//		
		//		AppendToGraph $Y_end vs $X_end
		//		ModifyGraph mode($NameOfWave(Yend))=3,marker($NameOfWave(Yend))=55, msize($NameOfWave(Yend))=8
		//		ModifyGraph useMrkStrokeRGB($NameOfWave(Yend))=1
		
		If ((f==0)&&(numtracks==1))
			//Display panel to indicate data is loading.
			NewPanel /N=frameWin1 /FLT=2 /W=(800,400,1000,450) as "Walking..."
			ModifyPanel cbRGB=(10583,50583,10583), frameStyle=1, NoEdit=1
			DrawText 21,20,"Please wait."
			SetActiveSubwindow _endfloat_
			DoUpdate
		Endif

	Endfor //EndFor statement for each displacement
	
	If ((f==numDisplacements)&&(numtracks==1))
		KillWindow frameWin1			// Kill panel
	Endif

	//Calculate average length of net BM
	Variable BM_avgLength = mean(BM_distancesWave)

	If (numDisplacements!=1)


		//Concatenate the x and y path lengths (displacements) and then build histogram; mean of lengths shoulbe be equal to zero
		Concatenate /NP /O {$BM_XWave,$BM_YWave}, $BM_XYwave
		Make /O /N=100  $BMhistwave
		Histogram /C /B=1 $BM_XYwave,$BMhistwave
		Wave newhistwave = $BMhistwave
	
		Display /N=BM_histgraph1 /K=1 /W=(600,685,1000,880 ) $BMhistwave
		ModifyGraph offset($NameOfWave(newhistwave))={-deltaX(newhistwave),0}
		ModifyGraph mode=5,hbFill=2,rgb=(26214,26214,26214),useBarStrokeRGB=1
		Label left "N"
		Label bottom "Displacement"
		SetDataFolder root:BM
		If (numTracks!=1)
			CurveFit/M=2 /Q=1 /W=2 gauss, $BMhistwave /D
		Else
			CurveFit/M=2/W=0 gauss, $BMhistwave /D
		Endif

		WAVE W_coef
		Variable width=W_coef[3]
		
		//convert width of histogram to sigma squared (variance)
		Variable fitVariance = (width/sqrt(2))^2
		
		
		//Calculate diffusion from MSD (variance) = 2Dt (t = # collisions)
		//Note nominal D is meanpath^2/2  from k^2 = (2Dt), where k= mean unit of displacement and t =1 collision
		
		Variable meanpath = (0+0+numpath^2+numpath^2+(numpath^2)/2+(numpath^2)/2+(numpath^2)/2+(numpath^2)/2)/ 8
		
		Variable Diff = fitVariance/(2*numCollisions)
		
		SetDataFolder savedDF

		TextBox/C/N=textbox2/A=LT /X=0/Y=0 "D = "+num2str(Diff)+"\r\Z08Nominal D = "+num2str(((meanpath)/2))


		If (tempNum==numtracks)
			//Histogram of all displacements
			Make/N=100/O  $histwave
			Variable binNum =sqrt(numCollisions)
			Histogram  /C /B=1 $BM_distances,$histwave

			//Display graph of displacements from each random walk taken (total # of displacements; if multiple tracks)
			Display /N=BM_histgraph2 /K=1 /W=(600,480,1000,675 ) $histwave
			ModifyGraph mode=5,hbFill=2,rgb=(26214,26214,26214),useBarStrokeRGB=1
			TextBox/C/N=textbox2/A=LT /X=0/Y=0 "Histogram of "+num2str(numDisplacements)+" displacements"+"\r\Z10mean displacement = "+num2str(Round(BM_avgLength))
			Label left "N"
			Label bottom "rms displacement"

			Display /N=BM_graph2 /K=1 /W=(200,480,600,880 ) $Y_endsum vs $X_endsum
			ModifyGraph rgb=(0,0,0), Mode=2, mirror=1
			SetAxis/A/N=1/E=2 left
			SetAxis/A/N=1/E=2 bottom
			TextBox/C/N=textbox1/A=RB /X=0/Y=0 num2str(numDisplacements)+" displacements"
			Label left "ÆY"
			Label bottom "ÆX"
			DoUpdate
			GetAxis left
			Variable tempMax1 = V_max
			GetAxis bottom
			Variable tempMax2 = V_max
			If (tempMax1<tempMax2)
				SetAxis left, -tempMax2,tempMax2
			Elseif (tempMax2<tempMax1)
				SetAxis bottom, -tempMax1,tempMax1
			Endif
			
			AppendXYZContour :BM:BM_distances vs {:BM:BM_X,:BM:BM_Y}; DelayUpdate
			ModifyContour BM_distances labelBkg=3, labelFSize=8, labelFStyle=0, labelHV=4,  autoLevels={*,*,6}
			DoUpdate
		
			
			//Display Track
			Display /N=BM_graph3 /K=1 /W=(200,50,600,450 ) $runningBM_YWave vs $runningBM_XWave
			ModifyGraph rgb=(0,0,0)
			ModifyGraph mirror=1
			SetAxis/A/N=1/E=2 left
			SetAxis/A/N=1/E=2 bottom
			TextBox/C/N=textbox1/A=LT /X=0/Y=0 "\Z12Track #"+num2str(numtracks)+"\r\Z10("+num2str(numDisplacements)+" displacements)"
			Label left "Y"
			Label bottom "X"
			Tag/C/N=text0_1/A=LB /F=0 /L=2 /TL={lineRGB=(2,39321,1 )} BM_Y_tally, 0,"\K(2,39321,1)Start"
			Tag/C/N=text0_2/A=LB /F=0 /L=2 /TL={lineRGB=(65535,0,0 )} BM_Y_tally, numpnts($runningBM_YWave),"\K(65535,0,0)End"
		Endif
	Endif

	//Display the graph of the last displacement
	If (tempNum==numtracks)
		Display /N=BM_graph1 /K=1 /W=(600,50,1000,450) $Y_steps vs $X_steps m
		ModifyGraph rgb=(0,0,0)
		ModifyGraph mirror=1
		SetAxis/A/N=1/E=2 left
		SetAxis/A/N=1/E=2 bottom
		Label left "Y"
		Label bottom "X"
		DoUpdate
		GetAxis left
		variable tempMax3 = V_max
		GetAxis bottom
		variable tempMax4 = V_max
		If (tempMax3<tempMax4)
			SetAxis left, -tempMax4,tempMax4
		Elseif (tempMax4<tempMax3)
			SetAxis bottom, -tempMax3,tempMax3
		Endif		
	
		TextBox/C/N=textbox1/A=LT /X=0/Y=0 "\Z12Path of "+num2str(numcollisions)+" collisions \rfor the last displacement of Track #"+num2str(numtracks)+ "\r\Z10(displacement = "+num2str(Round(BMlength))+")"
						
		AppendToGraph $Y_start vs $X_start
		ModifyGraph mode($NameOfWave(Ystart))=3,marker($NameOfWave(Ystart))=16, msize($NameOfWave(Ystart))=8
		ModifyGraph useMrkStrokeRGB($NameOfWave(Ystart))=1
		ModifyGraph rgb($NameOfWave(Ystart))=(0,0,0)
		Tag/C/N=text0_3/A=MC /F=0 /L=0 BM_startYval, 0,"Start"

		
		AppendToGraph $Y_end vs $X_end
		ModifyGraph mode($NameOfWave(Yend))=3,marker($NameOfWave(Yend))=17, msize($NameOfWave(Yend))=8
		ModifyGraph useMrkStrokeRGB($NameOfWave(Yend))=1
		ModifyGraph rgb($NameOfWave(Yend))=(0,0,0)
		Tag/C/N=text0_4/A=MC /F=0 /L=0 BM_endYval, 0,"End"
			
	Endif



	Return Diff
	
End

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////

Function runBM2(tempNum)
	Variable tempNum
	Variable i

	String BM_Dwave = "root:BM:BM_Dvalues"
	//Make a wave to store each displacement
	Make /O /N=(tempNum) $BM_Dwave
	Wave Dwave = $BM_Dwave

	For (i=1; i<=tempNum;i+=1)
		Variable D = runBM1(i)
		Dwave[i-1]=D
	
		If (i==1)
			//Display panel to indicate data is loading.
			NewPanel /N=loopWin1 /FLT=2 /W=(800,400,1000,450) as "Walking..."
			ModifyPanel cbRGB=(10583,50583,10583), frameStyle=1, NoEdit=1
			DrawText 21,20,"Please wait."
			SetActiveSubwindow _endfloat_
			DoUpdate
		Elseif (i==tempNum)
			KillWindow loopWin1			// Kill panel
		Endif
		
	Endfor

	DoWindow /K D_table1
	Edit /N=D_table1 /K=1 Dwave

End

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
//function to test random number generatory
Function makeRandom_Numbers(tempNum)

	Variable tempNum
	Variable i,z
	String z_wave = "random_zWave"
	String h_wave = "random_zWave_hist"
	Make /O /N=(tempNum) $z_wave
	Wave zWave= $z_wave


	For (i=0;i<(tempNum);i+=1)
		z=sqrt((8*enoise(1))^2)
		z=Trunc(z)
			
		zWave[i]=z
	Endfor
	
	Make/N=100/O  $h_wave
	Histogram  /C /B=1 $z_wave,$h_wave
	Wave newhistwave = $h_wave

	
	//Display a graph of the random values
	DoWindow /K random_graph1
	Display /N=random_graph1 /K=1 /W=(200,50,800,450 ) $h_wave
	ModifyGraph offset($NameOfWave(newhistwave))={-deltaX(newhistwave),0}
	ModifyGraph mode=5,hbFill=2,rgb=(1,39321,19939),useBarStrokeRGB=1
	SetAxis/A/N=1/E=0 left
	SetAxis bottom -1,8
			
			
End

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////