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

// Date: 011/12/2012 //Version 1.2
// Procedure  for simulating a random walk to model Brownian Motion. A single random walk (frame) consists of some number of steps and results in a random mean displacement.
// A series of these random displacements (frames) 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 Carter (acarter@amherst.edu) at Amherst College

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

Menu "Macros"

	"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 D_table1
		KillWaves /Z /A
		KillDataFolder /Z  root:BM
	Endif

	If (cmpstr(ctrlName,"start")==0)
		NewDataFolder/O/S root:BM
		Variable /G stepNum = 1000
		Variable /G frameNum = 100
		Variable /G trackNum = 1
		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 ///FLT=2
	Else	
		NewPanel /N=BM_Panel0 /W=(right-400,top,right-100,200) /K=2 ///FLT=2
	Endif
	ModifyPanel /W=BM_Panel0 cbRGB=(62535,62535,62535), frameStyle=3, frameInset=5, noEdit=1, fixedSize=1


	//Info and controls that should always be displayed; i.e. controls above the tabs
	SetDrawLayer /W=BM_Panel0 ProgFront
	DrawRRect /W=BM_Panel0 10,35,290,140
	SetDrawEnv /W=BM_Panel0 fsize= 18, fillbgc=(49151,53155,65535)
	DrawText /W=BM_Panel0 40,70,"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={10,160},size={90,15},title="Steps:", win=BM_Panel0
	SetVariable BMsetvar0,fSize=10,limits={1,100000,0},value= root:BM:stepNum, win=BM_Panel0

	SetVariable BMsetvar1,pos={110,160},size={90,15},title="Frames:", win=BM_Panel0
	SetVariable BMsetvar1,fSize=10,limits={1,10000,0},value= root:BM:frameNum, win=BM_Panel0
	
	SetVariable BMsetvar2,pos={210,160},size={80,15},title="Tracks:", win=BM_Panel0
	SetVariable BMsetvar2,fSize=10,limits={1,100,0},value= root:BM:trackNum, win=BM_Panel0

	Button BMbutton0, pos={100,80}, size={120,40}, title="Start\rRandom Walk(s)", 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
	NVAR numSteps= root:BM:stepNum
	NVAR numframes= root:BM:frameNum
	NVAR numtracks= root:BM:trackNum
	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 a wave to store each frame's displacement
	Make /O /N=0 $BM_distances
	Wave BM_distancesWave = $BM_distances

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

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

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

	
	For (f=0; f<numframes; 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=(numSteps) $(X_steps)
		Make /O /N=(numSteps) $(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 frame
		framexsteps=0
		frameysteps=0
		Variable horiz=0
		Variable vert=0
	
		For (a=0; a<numSteps; a+=1)

			//Generate random number for vertical step
			z=sqrt((8*enoise(1))^2)
			z=Trunc(z)
 
			If (z==0)
				vert=vert+1
			Elseif (z==1)
				vert=vert-1
			Elseif (z==2)
				vert=vert+sqrt(0.5)  //Makes every "step" = 1 unit
				horiz=horiz+sqrt(0.5)
			ElseIf (z==3)
				horiz=horiz+1
			Elseif (z==4)
				horiz=horiz-1
			Elseif (z==5)
				vert=vert-sqrt(0.5)
				horiz=horiz+sqrt(0.5)
			Elseif (z==6)
				vert=vert+sqrt(0.5)
				horiz=horiz-sqrt(0.5)
			Elseif (z==7)
				vert=vert-sqrt(0.5)
				horiz=horiz-sqrt(0.5)
			Endif


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

		Endfor //end For statement for stepping

		//Calculate end X,Y position
		Variable endXpos = framexsteps[a]
		xend[0]=endXpos
		xendsum[f]=endXpos
		
		Variable endYpos = frameysteps[a]
		yend[0]=endYpos
		yendsum[f]=endYpos

		//calculate distance from start
		Variable BMlengthX = (endXpos-startXpos)
		BMX[f]=BMlengthX
		Variable BMlengthY = (endYpos-startYpos)
		BMY[f]=BMlengthY
		

		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=(30583,30583,30583), frameStyle=1, NoEdit=1
			DrawText 21,20,"Please wait."
			SetActiveSubwindow _endfloat_
			DoUpdate
		Endif
	Endfor //end For statement for each frame
	If ((f==numframes)&&(numtracks==1))
		KillWindow frameWin1			// Kill panel
	Endif

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

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

	If (numframes!=1)


		//Concatenate the x and y path lengths 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,480,1000,700 ) $BMhistwave
		ModifyGraph offset($NameOfWave(newhistwave))={-deltaX(newhistwave),0}
		ModifyGraph mode=5,hbFill=2,rgb=(26214,26214,26214),useBarStrokeRGB=1
		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 to sigma squared (variance)
		Variable fitVariance = (width/sqrt(2))^2 
		//Calculate diffusion from MSD (variance) = 2Dt (t = # steps)
		//Note nominal D is (1/sqrt2)^2 = 1/2 = 0.5 from k = sqrt(2Dt), where k=1 displacement and t =1 step
		//Average of X^2 = 0.5 ==> 0+0+1+1+0.5+0.5+0.5+0.5 / 8
		Variable Diff = fitVariance/(2*numSteps)
		SetDataFolder savedDF

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


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

			//Display graph of displacements from each random walk taken (frame)
			Display /N=BM_histgraph2 /K=1 /W=(200,480,600,700 ) $histwave
			ModifyGraph mode=5,hbFill=2,rgb=(26214,26214,26214),useBarStrokeRGB=1
			TextBox/C/N=textbox2/A=LT /X=0/Y=0 "Mean displacement length = "+num2str(Round(BM_avgLength))+"\r\Z08Histogram of displacements for "+num2str(numframes)+" frames."

			Display /N=BM_graph2 /K=1 /W=(600,50,1000,450 ) $Y_endsum vs $X_endsum
			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=RB /X=0/Y=0 "Consecutive displacements for the "+num2str(numframes)+" frames"
			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		
			
		Endif
	Endif

	If (tempNum==numtracks)
		//Display the graph of the last random walk taken
		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
			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		
	
		If(numframes!=1)	
			TextBox/C/N=textbox1/A=RB /X=0/Y=0 "\Z12Displacement length = "+num2str(Round(BMlength))+"\r\Z08This graph is the walk for frame #"+num2str(numframes)
		Else
			TextBox/C/N=textbox1/A=RB /X=0/Y=0 "\Z12Displacement length = "+num2str(Round(BMlength))
		Endif
			
		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
	Endif
	
	Return Diff
	
End

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

Function runBM2(tempNum)
	Variable tempNum
	Variable i

	String BM_Dwave = "root:BM:BM_Dvalues"
	//Make a wave to store each frame's 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=(30583,30583,30583), 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

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