#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=1		// Use modern global access method.

//reserved letters: e i j p q r s t x y z
// Last changed February 17, 2017 by IP

function PMF(F,N,DL,L0,G0,M0)
	variable F,n,L0,DL,G0,M0
	variable Lmax,kT,p,L,i,j,Dx,R,b,w, xb
	NVAR p_val= root:SIM:p_val
	NVAR sigmaG = root:SIM:sigmaG
	make/O/N=(N+1) pl=p_val
	make/O/N=(N+2,2) MinUxy=NaN
	//Morse parameters
	//M0=19.27	// in pN nm
	R=2.8//6//4				// size of a folded domain - changed to 2.49*sqrt(6)=6
	b=40			// Morse width
	//Gaussian parameters
	//G0=32.8  	// in pN nm
	w=sigmaG*sqrt(2)//.2 //this is sigma*sqrt(2), Gaussian width
	xb=0.44		//.3  this is delta X, distance to transition state
	kT=4.11		// thermal energy
	Dx=0.001		// resolution landscape in nm
	///////////////get U setup as a template///////
	L=L0+N*DL
	make/O/N=(L/Dx) U
	setscale/P x,0,Dx,U
	//get shortest WLC; E-1
	L=L0
	U=(kT/pl[0])*((L/4)*(1/(1-x/L)-1)-(x/4)+x^2/(2*L))-F*x
	FindAPeak 0.1,2,1, U
	MinUxy[1][0]=V_peakX
	Duplicate/O/R=(0,MinUxy[1][0]) U Cut_U
	string Ui="U"+num2str(1)
	string Si="S"+num2str(1)
	string Ti="T"+num2str(1)
	Duplicate/O Cut_U $Si
	Duplicate/O Cut_U $Ti
	Duplicate/O U $Ui

	////////////////do the rest of the WLC///////////////////
	for(i=2;i<N+2;i+=1)
		Ui="U"+num2str(i)
		L=L0+(i-1)*DL
		U=(kT/pl[i])*((L/4)*(1/(1-x/L)-1)-(x/4)+x^2/(2*L))-F*x
		FindAPeak 0.1,2,1, U
		MinUxy[i][0]=V_peakX
		Duplicate/O/R=(MinUxy[i-1][0],MinUxy[i][0]) U Cut_U
		Ui="U"+num2str(i)
		Si="S"+num2str(i)
		Ti="T"+num2str(i)
		Duplicate/O Cut_U $Si
		Duplicate/O Cut_U $Ti
		Duplicate/O U $Ui
	endfor
	////////////////Add the Enthalpy to each segment and shift///////////////////
	WAVE S1
	Deletepoints numpnts(S1)-1,1,S1
	for(i=2;i<N+2;i+=1)
		string S0="S"+num2str(i-1)
		Duplicate/O $S0 H0
		Si="S"+num2str(i)
		Duplicate/O $Si H
		H=H+M0*((1-exp((-2*b/R)*((x+R-pnt2x(H,0))-R)))^2-1)+G0*Exp(-((x-pnt2x(H,0)-xb)/w)^2)
		variable shift=H[0]-H0(numpnts(H0))
		H=H-shift
		Deletepoints numpnts(H)-1,1,H
		Duplicate/O H $Si
	endfor
	////////////////create final segment and shift///////////////////
	Si="S"+num2str(N+2)
	Ui="U"+num2str(N+1)
	Duplicate/O/R=(MinUxy[N+1][0],MinUxy[N+1][0]+100) $Ui Last
	shift=H(numpnts(H))-Last[0]
	Last=Last+shift
	Duplicate/O last $Si
	////////////////Concatenate for final wave///////////////////
	string full=""
	for(i=1;i<N+3;i+=1)
		full=full+"S"+num2str(i)+";"
	endfor
	concatenate/O full, PMF_complete

	//fill PMF_complete with NAN's above 500 pN nm
	findlevel/P/Q PMF_complete, 500
	for(j=round(V_LevelX); j<numpnts(PMF_complete);j+=1)
		PMF_complete[j]=nan
	endfor

	//save MinUy at the force used
	variable z
	for(i=0;i<N+2;i+=1)
		z=PMF_complete(MinUxy[i][0])
		MinUxy[i][1]=z
	endfor
	//remove all other Waves//
	for(i=1;i<N+2;i+=1)
		Si="S"+num2str(i)
		Ui="U"+num2str(i)
		Ti="T"+num2str(i)
		killwaves/Z $Si, $Ui,$Ti
	endfor
	Si="S"+num2str(N+2)
	killwaves/Z $Si
	killwaves/Z H,H0,Last,Cut_U,U,PL
End


Function E_curves(N,delta_F,nF,DL,L0,G0,M0)
	variable N,delta_F,nF,Dl,L0,G0,M0
	variable j,i,F
	variable nmax
	make/O PMF_complete//dummy wave
	make/O/N=(N+2,2) MinUxy=NaN //dummy wave
	make/O/N=(N+2,2,nF) root:MinXYF=NaN
	wave MinXYF = root:MinXYF
	make/O/N=(nF) root:F_table
	wave F_table = root:F_table
	F_Table=1+delta_F*(nF-p)//the minimum force is always 1 pN............
	//see how long is the longest PMF at the highest force!!
	PMF(F_table[0],N,DL,L0,G0,M0)
	nmax=numpnts(PMF_complete)
	make/O/N=(nmax,nF) root:PMFM=NaN
	wave PMFM = root:PMFM
	setscale/P x,0, deltax(PMF_complete), PMFM
	setscale/P y,F_table[0], -delta_F, PMFM//set the scale in descending mode
	//create #nF PMF's and their Minima//
	for(i=0;i<nF;i+=1)
		F=F_table[i]
		PMF(F,N,DL,L0,G0,M0)//returns both PMF_complete and MinUxy
		PMFM[][i]=PMF_complete[p]
		MinXYF[][][i]=MinUxy[p][q]
	endfor
	MinXYF[0][0][]=0
	Display_PMF(N,nF)
	killwaves/Z MinUxy, PMF_complete
	setscale/I z,(1+delta_F*nF),(1+delta_F), MinXYF

End


Function Display_PMF(N,nF)
	variable N, nF
	variable i
	Wave MinxyF = root:MinxyF
	Wave PMFM = root:PMFM
	DoWindow/K PMF0
	//display E_curves first
	Display/W=(500,10,1200,600)/N=PMF0/K=1 MinxyF[0][1][] vs MinxyF[0][0][]
	for(i=1;i<N+2;i+=1)
		Appendtograph MinxyF[i][1][] vs MinxyF[i][0][]
	endfor
	ModifyGraph lsize=2, RGB=(65535,0,0)
	//now display the PMF's
	for(i=0;i<nF+1;i+=10)//show only every 10 curves
		Appendtograph/C=(1,12815,52428) PMFM[][i]
	endfor
	Global_min()
	wave G_min = root:G_min
	AppendToGraph G_min[][1] vs G_min[][0]
	ModifyGraph lsize(G_min)=3,rgb(G_min)=(2,39321,1)
End

function Global_min()
	NVAR N_modules = root:sim:N_modules
	NVAR L0 = root:sim:L0
	NVAR Lc = root:sim:Lc
	NVAR U0 = root:sim:U0
	NVAR A0 = root:sim:A0
	NVAR deltaF = root:SIM:deltaF
	NVAR noF = root:SIM:noF

	wave MinXYF = root:MinXYF
	wave F_table = root:F_table
	make/O/N=(dimsize(MinXYF,2)-1,4) root:G_min
	wave G_min = root:G_min
	G_min[][2]=F_table[p]
	variable i
	make/O/N=(dimsize(MinXYF,0)) PMF_min
	for(i=0;i<(dimsize(MinXYF,2));i+=1)
		PMF_min= MinXYF[p][1][i]
		wavestats/Q PMF_min
		G_min[i][1]=MinXYF[V_minloc][1][i]
		G_min[i][0]=MinXYF[V_minloc][0][i]
		G_min[i][4]=V_minloc-1 // No. unfolded domains
	endfor
	Make/O/N=(numpnts(F_table)) DummyN
	DummyN=G_min[p][3]
	SetScale/P x, F_table[0], F_table[1]-F_table[0], DummyN
	Findlevel/Q DummyN,N_modules/2
	Printf "Half force is at %1.2f pN", V_LevelX
	Print " for U0:", U0/4.11, " KT, A0:", A0/4.11, "KT"
	killwaves/Z PMF_min

End


Function Display_F_X()
	wave MinXYF = root:MinXYF
	wave G_min = root:G_min
	wave F_table = root:F_table
	make/O/N=(dimsize(MinXYF,2),2) E_1
	E_1[][0]= MinXYF[1][0][p]
	E_1[][1]= F_table[p]
	DoWindow/K F_X
	//display E_1
	Display/W=(500,10,1200,600)/N=F_X/K=1 E_1[][1] vs E_1[][0]
	SetAxis left 0,25;DelayUpdate
	SetAxis bottom 0,120
	ModifyGraph lsize=2,rgb=(1,3,39321)
	AppendToGraph G_min[][2] vs G_min[][0]
	ModifyGraph lsize=2,rgb(G_min)=(65000,0,0)
End



/////////////////////New Code from here for GUI
Menu "SIM"
	"-"
	"Start Sim Panel", InitializeVars()

End

Function InitializeVars()
	NewDataFolder/O root:SIM
	Make/O/N=3 root:SIM:CommandF, root:SIM:Force_List={3,20,10},root:SIM:Duration_List={2,2,2}, root:SIM:Ramp_List={0,0,0}, root:SIM:SIM_PMF01=NAN
	Variable/G root:SIM:N_modules = 8
	Variable/G root:SIM:Lc = 18.6//19
	Variable /G root:SIM:L0 = 0.1
	Variable/G root:SIM:U0_KT = 1.7//3
	NVAR U0_KT=root:SIM:U0_KT
	Variable/G root:SIM:U0 = U0_KT*4.11
	//	SetFormula root:SIM:U0, "root:SIM:U0_KT*4.11"
	Variable/G root:SIM:A0_KT = 14.3//13
	NVAR A0_KT=root:SIM:A0_KT
	Variable/G root:SIM:A0 = A0_KT*4.11
	//	SetFormula root:SIM:A0, "root:Sim:A0_KT*4.11"
	Variable/G root:SIM:p_val = 0.58// persistence length
	Variable/G root:SIM:deltaF = 1	// min force step
	Variable/G root:SIM:noF = 150	// no steps in force
	Variable/G root:SIM:sigmaG = 0.1	// which comes from (2.51-2.49)*sqrt(6)=0.05 nm
	Variable/G root:SIM:dt=1e-7	// sampling time
	Variable/G	root:SIM:D_coef=15000


	BuildSimPanel()
	UpdatePMFMatrix("",0,"","")
	CreateSimFCCommand()

End

Function BuildSimPanel()
	Variable PanelXSize	=	1090
	Variable PanelYSize	=	550
	Make/O root:SIM:SIM_PMF01=NAN
	PauseUpdate
	String FldrSave= GetDataFolder(1)
	//	Define Position and Dimensions
	GetWindow kwFrameInner, wsizeDC				//	Get size of the Igor Main window
	Variable XOffset	=	10
	Variable	Left		=	V_right-PanelXSize-XOffset	//	The panel will be on the right side
	Variable	Top		=	0
	Variable	right		=	V_right-XOffset
	Variable bottom	=	PanelYSize+Top
	//	Create the panel
	DoWindow/K MainFrame
	NewPanel/N=MainFrame /W=(left,top,right,bottom) as "SIM Mainframe"
	SetWindow MainFrame,hook(testhook)=KeyboardPanelHook
	// The globals
	SetVariable NDisplaySet, pos={49,32},size={58,21},title="N",variable=root:SIM:N_modules,proc = UpdatePMFMatrix
	SetVariable NDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.0f"

	SetVariable U0DisplaySet,pos={123,32},size={110,21},title="U0, kT", variable=root:SIM:U0_kT, proc = UpdatePMFMatrix
	SetVariable U0DisplaySet,limits={1,inf,1}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"

	SetVariable A0DisplaySet,pos={232,32},size={115,21},title="A0, kT",variable=root:SIM:A0_kt, proc = UpdatePMFMatrix
	SetVariable A0DisplaySet,limits={1,inf,1}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"

	SetVariable L0DisplaySet,pos={12,77},size={100,21},title="L0, nm",variable=root:SIM:L0, proc = UpdatePMFMatrix
	SetVariable L0DisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"

	SetVariable LcDisplaySet,pos={121,77},size={103,21},title="Lc, nm",variable=root:SIM:Lc, proc = UpdatePMFMatrix
	SetVariable LcDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"

	SetVariable pDisplaySet,pos={237,77},size={89,2},title="p, nm",variable=root:SIM:p_val, proc = UpdatePMFMatrix
	SetVariable pDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"

	SetVariable DDisplaySet,pos={222,121},size={120,21},title="D,nm2/s",variable=root:SIM:D_coef
	SetVariable DDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.0f"

	SetVariable deltaFDisplaySet,pos={20,121},size={92,21},title="dF, pN",variable=root:SIM:deltaF, proc = UpdatePMFMatrix
	SetVariable deltaFDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"

	SetVariable noFDisplaySet,pos={150,121},size={60,21}, title="nF",variable=root:SIM:noF, proc = UpdatePMFMatrix
	SetVariable noFDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.0f"


	//The protocol and PMF figures
	SetActiveSubwindow MainFrame
	Display /HIDE=0 /W=(350,30,750,280)/HOST=MainFrame /N=ProtocolWavePlot root:SIM:CommandF
	Label/W=MainFrame#ProtocolWavePlot left "Force Protocol (pN)"
	Label/W=MainFrame#ProtocolWavePlot bottom "Time (s)"
	ModifyGraph/W=MainFrame#ProtocolWavePlot tick=2,mirror=1, lsize=2, zero(left)=3, tickUnit=1

	Display /HIDE=0 /W=(350,300,750,550)/HOST=MainFrame /N=PMFPlot root:SIM:SIM_PMF01
	Label/W=MainFrame#PMFPlot left "PMF(pN*nm)"
	Label/W=MainFrame#PMFPlot bottom "Extension (nm)"
	ModifyGraph/W=MainFrame#PMFPlot tick=2,mirror=1, lsize=2, zero(left)=3, tickUnit=1

	//	FC > Tables
	Edit /HIDE=0/HOST=MainFrame/W=(800,30,998,530)/N=ProtocolForceTable root:SIM:Force_List,root:SIM:Duration_List, root:SIM:Ramp_List
	SetWindow MainFrame#ProtocolForceTable , userdata(tabnum)=  "1"
	SetWindow MainFrame#ProtocolForceTable , userdata(tabcontrol)=  "tabequipment"
	ModifyTable/W=MainFrame#ProtocolForceTable alignment=1, width=65
	ModifyTable/W=MainFrame#ProtocolForceTable width(Point)=35, showParts=199
	ModifyTable/W=MainFrame#ProtocolForceTable title(root:SIM:Ramp_List)="Ramp on?",Title(root:SIM:Force_List)="Force",Title(root:SIM:Duration_List)="Time"
	SetActiveSubwindow MainFrame

	// FC > Button > Delete Rows
	Button RemoveFromForceList,pos={1012,30},size={65,36},disable=0,proc=FCSimProtocolDeleteLastRow,title="Delete\rLast Row"

	SetDataFolder FldrSave
	ResumeUpdate
End

Function FCSimProtocolDeleteLastRow(ba) : ButtonControl
	STRUCT WMButtonAction &ba

	switch( ba.eventCode )
		case 2: // mouse up
			Wave w1 = root:SIM:Force_List
			Wave w2 = root:SIM:Duration_List
			Wave w3 = root:SIM:Ramp_List
			Variable n = Max(Max(NumPnts(w1),NumPnts(w2)),Max(NumPnts(w3),1))-1
			If(n>0)
				DeletePoints n,1, w1,w2,w3
			EndIF
			break
	endswitch

	return 0
End

Function UpdatePMFMatrix(ctrlName,varNum,varStr,varName) : SetVariableControl
	String ctrlName
	Variable varNum
	String varStr, varName

	NVAR N_modules = root:sim:N_modules
	NVAR L0 = root:sim:L0
	NVAR Lc = root:sim:Lc
	NVAR U0 = root:sim:U0
	NVAR A0 = root:sim:A0
	NVAR deltaF = root:SIM:deltaF
	NVAR noF = root:SIM:noF
	NVAR U0_KT=root:SIM:U0_KT
	NVAR A0_KT=root:SIM:A0_KT
	U0=U0_KT*4.11
	A0=A0_KT*4.11
	E_curves(N_modules,deltaF,noF,Lc,L0,A0,U0)
	CreateSimPMF()

End

static constant hook_keyboard=11
static constant hook_mouse=3
Function KeyboardPanelHook(s)
	STRUCT WMWinHookStruct &s
	NVAR X_Stage =root:Config:PiezoStage:X_Stage
	NVAR Y_Stage = root:Config:PiezoStage:Y_Stage
	NVAR  StepPiezoStage = root:Config:PiezoStage:StepPiezoStage

	switch(s.eventCode)
		case hook_keyboard:
			CreateSimFCCommand()
			CreateSimPMF()
			break
		case hook_mouse:
			CreateSimFCCommand()
			CreateSimPMF()
			break
	endswitch
	return 0
End

Function CreateSimFCCommand()
	WAVE	Force_List		=	root:SIM:Force_List
	WAVE	Duration_List		=	root:SIM:Duration_List
	WAVE	Ramp_List		=	root:SIM:Ramp_List
	NVAR deltaF = root:SIM:deltaF
	NVAR noF = root:SIM:noF
	NVAR	dt = root:SIM:dt
	Variable Duration	 =	Sum(Duration_List)
	Variable PointsPerTrace=round(Duration/(dt*1e3))
	Make/O/N=(PointsPerTrace)	root:SIM:CommandF=NaN	//	In units of pN
	Wave CF	=	root:SIM:CommandF
	//	Setscale/I x,0,Duration, "s", CF
	Setscale/P x, 0,dt*1e3, "s", CF
	Variable		Time1
	Variable		Time2
	Variable		Point1
	Variable 		Point2
	Variable 		i		=	0
	Variable 		nForces	=	NumPnts(Force_List)
	Variable		RampSlope
	Variable		F1, F2
	Do
		Time1	=	Time2
		Time2	+=	Duration_List[i]
		Point1	=	Point2
		Point2	=	x2Pnt(CF,Time2)
		F1		=	Force_List[i]
		if (F1>noF*deltaF)
			Print "Please increase the dF/nF to cover higher forces"
			Force_List[i]=noF*deltaF
			F1 = Force_List[i]
		endif
		if (F1<2)
			Print "Minimum force is 2 pN"
			Force_List[i]=2
			F1 = Force_List[i]
		endif
		if (Ramp_List[i]==1)
			//			F2					=	Force_List[i+1]
			//			RampSlope			=	(F2-F1)/(Point2 - Point1)
			//			CF[Point1, Point2]	=	F1+ RampSlope*(p-Point1)
			CF[Point1, Point2]	=	F1 //no ramps yet
		else
			CF[Point1, Point2]	=	F1
		endif
		i	+=	1
	While(i<nForces)

	WaveStats/Q CF
	DoWindow/F MainFrame
	if( V_Flag!=0 )
		SetAxis/W=MainFrame#ProtocolWavePlot left 0,1.1*V_max
	endif
End

Function CreateSimPMF()
	Wave	Force_List		=	root:SIM:Force_List
	Wave F_table = root:F_table
	Wave PMFM = root:PMFM
	NVAR N_modules = root:SIM:N_modules
	Variable i
	// remove old PMF waves from panel & kill them
	DoWindow/F MainFrame
	if( V_Flag!=0 )
		String FldrSave= GetDataFolder(1)
		SetDataFolder root:SIM
		do
			String list_string = TraceNameList("MainFrame#PMFPlot", ";",1) 	//This function makes a string with all the names of the waves in graph
			RemoveFromGraph /W=MainFrame#PMFPlot $StringFromList(0, list_string, ";")
			if(!StringMatch (StringFromList(0,list_string,";"), "MinXYF*"))
				KillWaves/Z $StringFromList(0, list_string, ";") // THIS IS WHY THE MinxyF & MinxyF SHOULD NOT BE IN THIS FOLDER
			endif
		while(ItemsInList	(list_string)!=1)
		SetDataFolder FldrSave
	endif
	//make new PMFs on which the simulation will run
	for(i=0;i<numpnts(Force_List);i+=1)
		String nameWave = "root:Sim:SIM_PMF" + num2str(round(Force_List[i]))
		Make/O/N=(dimsize(PMFM,0)) $nameWave
		WAVE workWave = $nameWave
		workWave = PMFM[x][binarysearch(F_table,Force_List[i])]
		setscale/P x,0, deltax(PMFM), workWave
		Appendtograph /C=(1,12815,52428)/W=MainFrame#PMFPlot $nameWave
		Label/W=MainFrame#PMFPlot left "PMF(pN*nm)"
		Label/W=MainFrame#PMFPlot bottom "Extension (nm)"
	endfor
	// Set axis before adding the E curves
	SetAxis/A/W=MainFrame#PMFPlot /N=1 left;
	SetAxis/A/W=MainFrame#PMFPlot /N=1 bottom
	DoUpdate
	GetAxis/Q/W=MainFrame#PMFPlot left
	SetAxis/W=MainFrame#PMFPlot left V_min, V_max
	for(i=0;i<=N_modules+1;i+=1)
		wave MinxyF = root:MinxyF
		Appendtograph/W=MainFrame#PMFPlot root:MinxyF[i][1][] vs root:MinxyF[i][0][]
	endfor
End

Function BuildSIMDisplay()
	Wave	F_wave		=	root:SIM:F_wave
	Wave	x_wave		=	root:SIM:x_wave
	Wave	CommandF		=	root:SIM:CommandF
	NVAR N_modules = root:sim:N_modules
	NVAR L0 = root:sim:L0
	NVAR Lc = root:sim:Lc

	DoWindow/K FcDualDisplay
	Display/K=1 /W=(0,40,550,450)/L=ForceAxis/N=FcDualDisplay root:SIM:CommandF, root:SIM:F_wave
	AppendToGraph/W=FcDualDisplay/L=left root:SIM:x_wave
	ModifyGraph rgb(F_wave)=(0,12800,52224)
	ModifyGraph/W=FcDualDisplay grid(left)=1,grid(ForceAxis)=1
	ModifyGraph/W=FcDualDisplay grid(bottom)=0,tick=2,mirror=1
	ModifyGraph/W=FcDualDisplay gridRGB(ForceAxis)=(34816,34816,34816)
	ModifyGraph/W=FcDualDisplay gridRGB(left)=(34816,34816,34816)
	ModifyGraph/W=FcDualDisplay nticks(left)=10
	ModifyGraph/W=FcDualDisplay lblPosMode(ForceAxis)=2,lblPosMode(left)=2
	ModifyGraph/W=FcDualDisplay lblPos(ForceAxis)=50,lblPos(left)=50
	ModifyGraph/W=FcDualDisplay axisEnab(ForceAxis)={0,0.35}
	ModifyGraph/W=FcDualDisplay axisEnab(left)={0.4,1}
	ModifyGraph/W=FcDualDisplay freePos(ForceAxis)=0
	ModifyGraph/W=FcDualDisplay lstyle=0,lsize(CommandF)=1.1 ,rgb(CommandF)=(0,0,0)
	ModifyGraph/W=FcDualDisplay live= 1
	Label/W=FcDualDisplay ForceAxis "Force (pN)"
	Label/W=FcDualDisplay bottom "Time (\\U)"
	Label/W=FcDualDisplay left "Length (nm)"
	WaveStats/Q CommandF
	SetAxis/W=FcDualDisplay ForceAxis, 0.8*V_min, 1.1*V_max
	SetAxis/W=FcDualDisplay left, 0, L0+(N_modules+1)*Lc
	DoUpdate/W=FcDualDisplay
End
Function Find_A_and_U()
	NVAR N_modules = root:sim:N_modules
	NVAR L0 = root:sim:L0
	NVAR Lc = root:sim:Lc
	NVAR U0 = root:sim:U0
	NVAR A0 = root:sim:A0
	NVAR deltaF = root:SIM:deltaF
	NVAR noF = root:SIM:noF
	NVAR U0_KT=root:SIM:U0_KT
	NVAR A0_KT=root:SIM:A0_KT

	Variable i, sumAU
	sumAU=16
	Make/O/N=(sumAU-1,3) UAN_matrix
	for(i=0;i<round(sumAU)-1;i+=1)
		U0_KT=i*0.1+1
		A0_KT=sumAU-U0_KT
		U0=U0_KT/4.11
		A0=A0_KT/4.11
		UpdatePMFMatrix("",0,"","")
		DoUpdate
		//Global_min()
		WAVE DummyN
		Findlevel/Q DummyN,4
		UAN_matrix[i][0]= U0_KT
		UAN_matrix[i][1] = A0_KT
		UAN_matrix[i][2] = V_LevelX
	endfor
	Display UAN_matrix[2][] vs UAN_matrix[0][]
End
Function plot_normal_distrib(sigma)
	Variable sigma
	Make/o/N=1000 normD
	Setscale/I x,-5*sigma, 5*sigma, normD
	normD=1/sqrt(2*(sigma^2)*pi)*exp(-(x)^2/(2*sigma^2))
End
//#pragma rtGlobals=	3		// Use modern global access method and strict wave access.
Function make_lines(x,y,z,conc,savepath)
	Variable x,y,z,conc
	string savepath
	if (x==0 || y==0 || z==0)
		abort "x, y, or z cannot be one. That's not even gel."
	endif
	NewDataFolder/O root:RESULTS
	Variable/G N_lines=x*y*z
	Variable/G xby=x
	Variable/G yby=y
	Variable/G zby=z
	Variable/G Length = 4 //nm
	Variable/G modules=8
	Variable/G conc_molec=conc//molec/nm^3 (from 1000 molec per 100 nm^3)
	Variable/G box_size
	Variable/G NumberPoints=1000
	Variable/G displayw=1
	Variable/G showconnect=0
	Variable/G LinkSize=4//twice a radius domain (2*2)
	string/G CrosslinkingPath,UnfoldingPath,OptimizationPath,mainpath
	mainpath=savepath
	NewPath/C/O/Q CrosslinkingPath savepath+":Crosslinking"
	NewPath/C/O/Q UnfoldingPath savepath+":Unfolding:"
	NewPath/C/O/Q OptimizationPath savepath+":Optimization"
	Variable i,j,k,l,m, CMdif,theta, phi,a,b,c,timecount
	String nameW, colorW, InfoNote, exec, NameSW
	Variable SetX, SetY, deltaStep, step, SetZ
	//Make the waves that define everything
	make/o/n=(n_lines,n_lines) root:results:ConnectMat=0
	make/o/n=(1,n_lines) root:results:a_list=0
	make/o/n=(1,n_lines) root:results:b_list=0
	make/o/n=(1,n_lines) root:results:c_list=0
	make/o/n=(1,n_lines) root:results:theta_list=0
	make/o/n=(1,n_lines) root:results:phi_list=0
	make/o/n=(n_lines) root:results:deltawave=0
	make/o/n=(modules,3,n_lines,0) root:results:LineMat
	wave linemat=root:results:linemat
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:results:theta_list
	wave phi_list=root:results:phi_list
	wave deltawave = root:results:deltawave
	make/o/n=(modules,n_lines) root:results:lengthmat=0
	make/o/n=(modules,n_lines) root:results:UnfoldMat=0
	variable xbox,ybox,zbox,deltaStepx,deltaStepy,deltaStepz
	variable/G connectPerClust=3
	variable/G Dtran=10^(-6)*4.11*(ln(modules)+0.312+0.565/modules-0.1/modules^2)/(3*pi*0.000911*10^(-6))/(length*modules)//nm^2
	variable/G Drot=10^(-6)*3*4.11*(ln(modules)-0.662+0.917/modules-0.050/modules^2)/(pi*0.000911*10^(-6))/((length*modules)^3)//nm^2
	variable/G Dein=10^(-6)*4.11/(6*pi*0.000911*10^(-6))//without rg
	//const is dt*kbT/(3*pi*eta), with kbT=4.11, dt=10^-6, and eta=0.000911 kg/(m*s) <- dynamic viscocity of water at 24 C
	N_lines=xby*yby*zby
	xbox=xby/(conc_molec)^(1/3)
	//xbox=(xby^2/yby/zby)*(N_lines/conc_molec)^(1/3)
	deltaStepx = xbox/(xby)
	//ybox=(yby^2/xby/zby)*(N_lines/conc_molec)^(1/3)
	ybox=yby/(conc_molec)^(1/3)
	deltaStepy = ybox/(yby)
	//zbox=(zby^2/xby/yby)*(N_lines/conc_molec)^(1/3)
	zbox=zby/(conc_molec)^(1/3)
	deltaStepz = zbox/(zby)
	ColorTab2Wave rainbow
	WAVE M_Colors
	if(displayW)
		DoWindow/K DisplayG
		exec = "NewGizmo /W = (10,0,1000,1000)/N=DisplayG; ModifyGizmo ShowInfo"
		Execute exec
	endif
	make/o/n=(modules,n_lines) Root:Results:UnfoldMat=0////rows define domain, columns define line
	i=0
	for(j=1;j<xby+1;j+=1)
		setX=j*deltaStepx//deltaStep
		for(k=1;k<yby+1;k+=1)
			setY=k*deltaStepy//deltaStep
			for(l=1;l<zby+1;l+=1)
				setZ=l*deltaStepz//deltaStep
				nameW="Line_"+num2str(i)
				Make/O/N=(modules,3) root:RESULTS:$nameW
				WAVE LineW = root:RESULTS:$nameW
				LineW=NAN
				InfoNote=""
				a=SetX+enoise(deltaStep/2)
				b=SetY+enoise(deltaStep/2)
				c=SetZ+enoise(deltaStep/2)
				theta=enoise(pi/2)+pi/2
				phi=enoise(pi)+pi
				LineW[][0][]=(p-modules/2)*Length*sin(theta)*cos(phi)+a
				LineW[][1][]=(p-modules/2)*Length*sin(theta)*sin(phi)+b
				LineW[][2][]=(p-modules/2)*Length*cos(theta)+c
				InfoNote += "Force=" + num2str(abs(cos(theta)))+ "; theta=" + num2str(theta) + "; phi=" + num2str(phi) + "; a=" + num2str(a) + "; b=" + num2str(b) + "; c=" + num2str(c)+ "; connections=0; startX="+num2str(setx)+"; startY="+num2str(sety)+"; startZ="+num2str(setz)+";"
				Note LineW, InfoNote
				if(displayW)
					exec = "ModifyGizmo setDisplayList=" + num2str(2*i)+ ", object=path"+ num2str(i)
					Execute exec

					exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:"+ nameofwave(LineW)
					Execute exec

					exec = "ModifyGizmo setDisplayList=" + num2str(2*i+1)+ ", object=scatter"+ num2str(i)
					Execute exec

					exec = "AppendToGizmo/N=DisplayG  DefaultScatter=root:RESULTS:"+ nameofwave(LineW)
					Execute exec

					exec ="ModifyGizmo ModifyObject=scatter"+ num2str(i)+", property={ size,0.2}"
					Execute exec
				endif
				i+=1
			endfor
		endfor
	endfor
	if(displayw)
		exec="ModifyGizmo scalingOption=0"
		execute exec
		exec="ModifyGizmo setOuterBox={"+num2str(-3*xbox)+","+num2str(3*xbox)+","+num2str(-3*ybox)+","+num2str(3*ybox)+","+num2str(-6*zbox)+","+num2str(6*zbox)+"}"
		Execute exec
		exec="modifygizmo home={115,-110,23};modifygizmo gohome"
		execute exec
	endif
	variable connects
	variable limits
	limits=1
	connections()
	variable savetime=0
	nvar totalnot
	timecount=0
	if (displayw)
		Saveit("Crosslinking","time_"+num2str(savetime))
	endif
	savetime+=1
	getaggregates()
	wave agg=root:results:aggwave
	do
		moveit(sqrt(2*Dtran),sqrt(2*Drot))///
		for (i=0;i<numpnts(agg);i+=1)
			if (agg[i]>1)
				movecluster(i)
			endif
		endfor
		connections()
		getaggregates()
		timecount+=1
		if (mod(timecount,10)==0)
			Saveit("Crosslinking","time_"+num2str(savetime))
			sleep/s 0.1
			savetime+=1
		endif
	while (numpnts(agg)!=1)
	if(displayw)
		exec="ModifyGizmo scalingOption=0"
		execute exec
		exec="ModifyGizmo setOuterBox={"+num2str(-xbox-modules*length)+","+num2str(xbox+modules*length)+","+num2str(-ybox-modules*length)+","+num2str(ybox+modules*length)+","+num2str(-zbox-modules*length)+","+num2str(zbox+modules*length)+"}"
		Execute exec
		exec="modifygizmo home={115,-110,23};modifygizmo gohome"
		execute exec
	endif
	if(showconnect)
		exec = "ModifyGizmo setDisplayList=" + num2str(2*N_lines-1)+ ", object=path"+ num2str(N_lines)
		Execute exec

		exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:Crosslinks"
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={pathColor,0,0,0,1}"///65535
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={lineWidthType,1}"
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={linewidth,5}"
		Execute exec
	endif
	connections()
	Saveit("Crosslinking","time_"+num2str(savetime))
	savetime+=1
	for(i=0;i<n_lines;i+=1)
		namew="Line_"+num2str(i)
		wave linew=root:results:$namew
		infonote=note(linew)
		a=NumberByKey(" a", infonote,"=")
		a_list[0][i]=a
		b=NumberByKey(" b", infonote,"=")
		b_list[0][i]=b
		c=NumberByKey(" c", infonote,"=")
		c_list[0][i]=c
		theta=NumberByKey(" theta", InfoNote, "=")
		theta_list[0][i]=theta
		phi=NumberByKey(" phi", InfoNote, "=")
		phi_list[0][i]=phi
	endfor
	//doproteinshift()
	CreateSizeW()
	updatecolors()
	variable xcm,ycm,zcm,r
	for(r=0;r<n_lines;r+=1)
		wave linew=root:results:$("Line_"+num2str(r))
		linew[][0]=(p-modules/2)*(Length)*sin(theta_list[dimsize(theta_list,0)-1][r])*cos(phi_list[dimsize(phi_list,0)-1][r])+a_list[dimsize(a_list,0)-1][r]
		linew[][1]=(p-modules/2)*(Length)*sin(theta_list[dimsize(theta_list,0)-1][r])*sin(phi_list[dimsize(phi_list,0)-1][r])+b_list[dimsize(b_list,0)-1][r]
		linew[][2]=(p-modules/2)*(Length)*cos(theta_list[dimsize(theta_list,0)-1][r])+c_list[dimsize(c_list,0)-1][r]
		for(i=0;i<modules;i+=1)
			xcm+=linew[i][0]/(n_lines*modules)
			ycm+=linew[i][1]/(n_lines*modules)
			zcm+=linew[i][2]/(n_lines*modules)
		endfor
	endfor
	for(r=0;r<n_lines;r+=1)
		wave linew=root:results:$("Line_"+num2str(r))
		linew[][0]-=xcm;linew[][1]-=ycm;linew[][2]-=zcm
	endfor
	variable/G timetotal=timecount*10^-4
	Saveit("Crosslinking","time_"+num2str(savetime))
	print num2str(getaggregates())+" clusters formed"
	//print "Gel polymerized in "+num2str(timecount*10^-4)+" seconds"
End
function moveit(stept,stepr)
	variable stept,stepr
	variable i
	variable connects,a,b,c,theta,phi,movea,moveb,movec,flag,movetheta,movephi,xbox,ybox,zbox
	string namew,infonote,exec
	nvar n_lines,modules,length,box_size,xby,yby,zby,conc_molec,Displayw
	variable pmx,pmy,pmz
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	xbox=(xby^2/yby/zby)*(N_lines/conc_molec)^(1/3)
	ybox=(yby^2/xby/zby)*(N_lines/conc_molec)^(1/3)
	zbox=(zby^2/xby/yby)*(N_lines/conc_molec)^(1/3)
	variable volumeFactor=3 // so the volume of the box is 1.5 times its intial volume
	pmx=xbox*(volumeFactor^(1/3)-1)/2
	pmy=ybox*(volumeFactor^(1/3)-1)/2
	pmz=zbox*(volumeFactor^(1/3)-1)/2
	for(i=0;i<n_lines;i+=1)
		namew="root:results:Line_"+num2str(i)
		wave linew=$namew
		infonote=note(linew)
		connects = NumberByKey(" connections", InfoNote, "=")
		if (connects==0)
			a=NumberByKey(" a", InfoNote, "=")
			b=NumberByKey(" b", InfoNote, "=")
			c=NumberByKey(" c", InfoNote, "=")
			theta=NumberByKey(" theta", InfoNote, "=")
			phi=NumberByKey(" phi", InfoNote, "=")
			do
				flag=0
				movea=gnoise(stept)
				moveb=gnoise(stept)
				movec=gnoise(stept)
				if ((movea+a)<(-pmx) || (movea+a)>(pmx+xbox))
					flag+=1
				endif
				if (flag==0)
					if ((moveb+b)<(-pmy) || (moveb+b)>(pmy+ybox))
						flag+=1
					endif
				endif
				if (flag==0)
					if ((movec+c)<(-pmz) || (movec+c)>(pmz+zbox))
						flag+=1
					endif
				endif
			while (flag)
			movetheta=gnoise(stepr*pi/180)
			movephi=gnoise(stepr*pi/180)
			LineW[][0][]=(p-modules/2)*Length*sin(theta+movetheta)*cos(phi+movephi)+a+movea
			LineW[][1][]=(p-modules/2)*Length*sin(theta+movetheta)*sin(phi+movephi)+b+moveb
			LineW[][2][]=(p-modules/2)*Length*cos(theta+movetheta)+c+movec
			infonote=ReplaceNumberByKey(" a",infonote,a+movea,"=")
			infonote=ReplaceNumberByKey(" b",infonote,b+moveb,"=")
			infonote=ReplaceNumberByKey(" c",infonote,c+movec,"=")
			infonote=ReplaceNumberByKey(" theta",infonote,theta+movetheta,"=")
			infonote=ReplaceNumberByKey(" phi",infonote,phi+movephi,"=")
			note/K linew,infonote
		endif
	endfor
end
function movecluster(clusternum)
	variable clusternum
	variable i,j
	variable connects,a,b,c,theta,phi,movea,moveb,movec,flag,movetheta,movephi,xbox,ybox,zbox,points
	string namew,infonote,exec
	variable cmx,cmy,cmz,pmx,pmy,pmz
	nvar n_lines,modules,length,box_size,xby,yby,zby,conc_molec,Displayw,dein
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	xbox=(xby^2/yby/zby)*(N_lines/conc_molec)^(1/3)
	ybox=(yby^2/xby/zby)*(N_lines/conc_molec)^(1/3)
	zbox=(zby^2/xby/yby)*(N_lines/conc_molec)^(1/3)
	pmx=xbox/xby
	pmy=ybox/yby
	pmz=zbox/zby
	wave clust=root:results:$("Cluster_"+num2str(clusternum))
	points=numpnts(clust)*modules
	make/o/n=(1,3) totalLines=0
	for(i=0;i<numpnts(clust);i+=1)
		wave linew=root:results:$("Line_"+num2str(clust[i]))
		Matrixop/o totalLines=catrows(totallines,linew)
		infonote=note(linew)
		a=NumberByKey(" a", InfoNote, "=")
		b=NumberByKey(" b", InfoNote, "=")
		c=NumberByKey(" c", InfoNote, "=")
		cmx+=a/numpnts(clust)
		cmy+=b/numpnts(clust)
		cmz+=c/numpnts(clust)
	endfor
	DeletePoints 0,1,totallines
	matrixop/o rg2=sumsqr(subtractmean(totallines,0))
	variable rg=sqrt(rg2[0]/points)
	variable step=sqrt(2*dein/rg)
	do
		flag=0
		movea=gnoise(step)
		moveb=gnoise(step)
		movec=gnoise(step)
		if ((movea+cmx)<(-pmx) || (movea+cmx)>(pmx+xbox))
			flag+=1
		endif
		if (flag==0)
			if ((moveb+cmy)<(-pmy) || (moveb+cmy)>(pmy+ybox))
				flag+=1
			endif
		endif
		if (flag==0)
			if ((movec+cmz)<(-pmz) || (movec+cmz)>(pmz+zbox))
				flag+=1
			endif
		endif
	while (flag)
	make/o/n=(numpnts(clust)*modules,3) root:Results:cluster_wave=0
	wave clust_wave=root:results:cluster_wave
	for(i=0;i<numpnts(clust);i+=1)
		wave linew=root:results:$("Line_"+num2str(clust[i]))
		infonote=note(linew)
		connects = NumberByKey(" connections", InfoNote, "=")
		a=NumberByKey(" a", InfoNote, "=")
		b=NumberByKey(" b", InfoNote, "=")
		c=NumberByKey(" c", InfoNote, "=")
		theta=NumberByKey(" theta", InfoNote, "=")
		phi=NumberByKey(" phi", InfoNote, "=")
		LineW[][0]=(p-modules/2)*Length*sin(theta)*cos(phi)+a+movea
		LineW[][1]=(p-modules/2)*Length*sin(theta)*sin(phi)+b+moveb
		LineW[][2]=(p-modules/2)*Length*cos(theta)+c+movec
		infonote=ReplaceNumberByKey(" a",infonote,a+movea,"=")
		infonote=ReplaceNumberByKey(" b",infonote,b+moveb,"=")
		infonote=ReplaceNumberByKey(" c",infonote,c+movec,"=")
		note/K linew,infonote
		for(j=i*modules;j<i*modules+modules;j+=1)
			clust_wave[j][]=linew[j-i*modules][q]
		endfor
	endfor
	//makerotationmat(gnoise(pi/2),gnoise(pi/2))//change this to be physical
	//wave rotationmat
	//matrixop/o clust_wave=(clust_wave-rowrepeat(sumcols(clust_wave)/points,points)) x rotationmat^t + rowrepeat(sumcols(clust_wave)/points,points)
	//make/o/n=(modules,3) tempwave
	//for(i=0;i<numpnts(clust);i+=1)
	//	wave linew=root:results:$("Line_"+num2str(clust[i]))
	//	for(j=i*modules;j<i*modules+modules;j+=1)
	//		linew[j-i*modules][]=clust_wave[j][q]
	//	endfor
	//	infonote=note(linew)
	//	theta=NumberByKey(" theta", InfoNote, "=")
	//	phi=NumberByKey(" phi", InfoNote, "=")
	//	tempwave[][0]=(p-modules/2)*Length*sin(theta)*cos(phi)
	//	tempwave[][1]=(p-modules/2)*Length*sin(theta)*sin(phi)
	//	tempwave[][2]=(p-modules/2)*Length*cos(theta)
	//	matrixop/o coords=sumcols(linew-tempwave)/modules
	//	infonote=ReplaceNumberByKey(" a",infonote,coords[0],"=")
	//	infonote=ReplaceNumberByKey(" b",infonote,coords[1],"=")
	//	infonote=ReplaceNumberByKey(" c",infonote,coords[2],"=")
	//	note/K linew,infonote
	//endfor
end
function makeRotationMat(b,c)
	variable b,c
	variable a=0
	wave RotationMat
	RotationMat[0][0]=cos(a)*cos(b);RotationMat[0][1]=cos(c)*sin(a)+cos(a)*sin(b)*sin(c);RotationMat[0][2]=cos(a)*cos(c)*sin(b)+sin(a)*sin(c)
	RotationMat[1][0]=cos(b)*sin(a);RotationMat[1][1]=cos(a)*cos(c)+sin(a)*sin(b)*sin(c);RotationMat[1][2]=cos(c)*sin(a)*sin(b)-cos(a)*sin(c)
	RotationMat[2][0]=-sin(b);RotationMat[2][1]=cos(b)*sin(c);RotationMat[2][2]=cos(b)*cos(c)
end
function connections()
	nvar n_lines,modules,box_size,numberpoints,connectPerClust,LinkSize
	variable i,j,k,l,cmdif
	string namew,infonote,exec
	wave connectmat = root:results:ConnectMat
	wave deltawave=root:results:deltawave
	nvar displayw,showconnect
	make/O/n=0 root:results:Conditions
	wave conditions = root:results:Conditions
	Make/O/N=(0,3) root:RESULTS:Crosslinks
	WAVE Crosslinks=root:RESULTS:Crosslinks
	Make/O/n=(0,5) root:RESULTS:Movement
	WAVE Movement=root:RESULTS:Movement
	for(i=0;i<N_lines;i+=1)
		for(j=i+1;j<N_lines;j+=1)
			WAVE L1=root:RESULTS:$("Line_"+num2str(i))
			WAVE L2=root:RESULTS:$("Line_"+num2str(j))
			for(k=0;k<modules;k+=1)
				variable counter=0
				if (counter>connectPerClust)
					break
				endif
				for(l=0;l<modules;l+=1)
					CMdif=sqrt((L1[k][0]-L2[l][0])^2+(L1[k][1]-L2[l][1])^2+(L1[k][2]-L2[l][2])^2)
					if(CMdif<LinkSize)//4 is the size of twice the radius of a domain
						InsertPoints/M=0 dimsize(Crosslinks,0),2,Crosslinks
						Crosslinks[dimsize(Crosslinks,0)-2][0]=L1[k][0]
						Crosslinks[dimsize(Crosslinks,0)-2][1]=L1[k][1]
						Crosslinks[dimsize(Crosslinks,0)-2][2]=L1[k][2]
						Crosslinks[dimsize(Crosslinks,0)-1][0]=L2[l][0]
						Crosslinks[dimsize(Crosslinks,0)-1][1]=L2[l][1]
						Crosslinks[dimsize(Crosslinks,0)-1][2]=L2[l][2]
						InsertPoints/M=0 dimsize(Crosslinks,0),1,Crosslinks
						Crosslinks[dimsize(Crosslinks,0)-1][]=NAN
						InsertPoints/M=0 dimsize(Movement,0),1,Movement
						connectmat[i][j]=1;connectmat[j][i]=1
						Movement[dimsize(Movement,0)-1][0]=i	//molecule number
						Movement[dimsize(Movement,0)-1][1]=k //domain number
						Movement[dimsize(Movement,0)-1][2]=j //molecule number
						Movement[dimsize(Movement,0)-1][3]=l //domain number
						Movement[dimsize(Movement,0)-1][4]=Cmdif //Distance
						counter+=1
					endif
				endfor
			endfor
		endfor
	endfor
	variable connect
	variable/G totalnot
	totalnot=0
	for(i=0;i<n_lines;i+=1)
		connect=0
		for(j=0;j<dimsize(movement,0);j+=1)
			if(movement[j][0]==i)
				connect+=1
			endif
			if(movement[j][2]==i)
				connect+=1
			endif
		endfor
		nameW="Line_"+num2str(i)
		WAVE LineW = root:RESULTS:$nameW
		InfoNote = note(root:RESULTS:$nameW)
		infonote=ReplaceNumberByKey(" connections", infonote, connect,"=")
		if (connect==0)
			totalnot+=1
		endif
		Note/K LineW, InfoNote
		variable interm
		if (displayw)
			ColorTab2Wave Rainbow
			WAVE M_Colors
			interm=(connect/modules)
			if (interm>1)
				interm=0.99
			endif
			Variable pick=((dimsize(M_Colors,0)-1)*(interm))
			exec ="ModifyGizmo ModifyObject=scatter"+ num2str(i)+", property={color,"+ num2str(M_Colors[pick][0]/65535) + ","+ num2str(M_Colors[pick][1]/65535) + ","+num2str(M_Colors[pick][2]/65535) + ",0}"
			Execute exec
		endif
	endfor
	if(showconnect)
		exec = "ModifyGizmo setDisplayList=" + num2str(2*N_lines-1)+ ", object=path"+ num2str(N_lines)
		Execute exec

		exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:Crosslinks"
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={pathColor,0,0,0,1}"//////65535
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={lineWidthType,1}"
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={linewidth,5}"
		Execute exec
	endif
end
function fixtemplooklist()
	wave w=root:templooklist
	variable i
	for(i=0;i<numpnts(w);i+=1)
		if (w[i]!=w[i])
			w[i]=-1
		endif
	endfor
end
function getaggregates()
	wave connectmat=root:Results:connectmat
	nvar n_lines
	variable i,j,k,aggcount,r,fnan,initnonnan,finalnonnan
	aggcount=0
	make/o/n=(N_lines) looklist
	make/o/n=0 root:results:aggwave
	wave aggwave=root:results:aggwave
	looklist[]=p
	aggcount=0
	fnan=0
	do
		duplicate/o looklist,templooklist
		fixtemplooklist()
		initnonnan=numbernonnan(looklist)
		fnan=firstnonnan(looklist)
		r=looklist[fnan]
		looklist[fnan]=nan
		make/o/n=0 dolist
		for(i=0;i<n_lines;i+=1)
			if (connectmat[r][i]==1)
				insertpoints numpnts(dolist),1,dolist
				dolist[numpnts(dolist)-1]=i
				looklist[i]=nan
			endif
		endfor
		k=0
		do
			k=dodolist(dolist,looklist)
		while (k==1)
		finalnonnan=numbernonnan(looklist)
		aggcount+=1
		compare(templooklist,looklist,numpnts(aggwave))
		insertpoints numpnts(aggwave),1,aggwave
		aggwave[numpnts(aggwave)-1]=initnonnan-finalnonnan
	while (isallnan(looklist)==0)
	return aggcount
end
function compare(w1,w2,clusterno)
	wave w1,w2
	variable clusterno
	variable len,i
	make/o/n=0 $("root:results:Cluster_"+num2str(clusterno))
	wave clust=$("root:results:Cluster_"+num2str(clusterno))
	len=numpnts(w1)
	for(i=0;i<len;i+=1)
		if (w2[i]!=w2[i])
			if (w1[i]!=-1)
				insertpoints numpnts(clust),1,clust
				clust[numpnts(clust)-1]=i
			endif
		endif
	endfor
end
function numbernonnan(w)
	wave w
	variable i,j
	for(i=0;i<(numpnts(w));i+=1)
		if(numtype(w[i])==0)
			j+=1
		endif
	endfor
	return j
end
function firstnonnan(w)
	wave w
	variable i
	for(i=0;i<numpnts(w);i+=1)
		if(numtype(w[i])==0)
			return i
		endif
	endfor
end
function isallnan(w)
	wave w
	variable i,r
	for(i=0;i<numpnts(w);i+=1)
		if (numtype(w[i])==0)
			return 0
		endif
	endfor
	return 1
end
function dodolist(dolist,looklist)
	wave dolist,looklist
	wave connectmat=root:results:connectmat
	nvar n_lines
	make/o/n=0 tempdolist
	variable i,r,j,t,find
	find=0
	for(i=0;i<numpnts(dolist);i+=1)
		r=dolist[i]
		for(j=0;j<n_lines;j+=1)
			if (connectmat[r][j]==1)
				if (isinwave(looklist,j))
					insertpoints numpnts(tempdolist),1,tempdolist
					tempdolist[numpnts(tempdolist)-1]=j
					looklist[j]=nan
					find=1
				endif
			endif
		endfor
	endfor
	if (numpnts(tempdolist)==0)
		return Find
	else
		redimension/N=(numpnts(tempdolist)) dolist
		dolist=tempdolist
		return find
	endif
end
function isinwave(w,pnt)
	wave w
	variable pnt
	variable j
	for(j=0;j<numpnts(w);j+=1)
		if(w[j]==pnt)
			return 1
		endif
	endfor
	return 0
end
function moveandconnect(stept,stepr)
	variable stept,stepr
	moveit(stept,stepr)
	connections()
end
Function CreateSizeW()
	NVAR N_lines
	Variable i
	String NameSW
	for(i=0;i<N_lines;i+=1)
		NameSW = "SizeW_"+num2str(i)
		Make/D/O/N=(8,3) root:RESULTS:$NameSW=0.2
	endfor
End
function linelength(w)
	// returns the line length its geometry
	wave w
	nvar modules
	return distance(w,0,w,modules-1)
end
function BuildLineCopies()
	// Builds a Copy of the current "Line_#" as "LineCopy_#"
	// Line copies are the ones modulated during the optimzaiton
	string namew,infonote,exec
	variable a,b,c,theta,phi,setx,sety,setz,deltastep,i,LengthLine
	nvar modules,length,n_lines,numberpoints
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:results:theta_list
	wave phi_list=root:results:phi_list
	wave movement=root:results:movement
	wave lengthmat=root:results:lengthmat
	wave unfoldmat=root:results:unfoldmat
	make/o/n=(dimsize(movement,0)) root:results:$("CriticalDis")
	wave criticaldis=root:results:$("CriticalDis")
	duplicate/o lengthmat,root:results:$("CompareMat")
	//criticaldis[]=movement[p][4]
	InsertPoints/M=0 dimsize(a_list,0), 1, a_list
	InsertPoints/M=0 dimsize(b_list,0), 1, b_list
	InsertPoints/M=0 dimsize(c_list,0), 1, c_list
	InsertPoints/M=0 dimsize(theta_list,0), 1, theta_list
	InsertPoints/M=0 dimsize(phi_list,0), 1, phi_list
	for(i=0;i<n_lines;i+=1)
		a_list[dimsize(a_list,0)-1][i]=a_list[dimsize(a_list,0)-2][i]
		b_list[dimsize(b_list,0)-1][i]=b_list[dimsize(b_list,0)-2][i]
		c_list[dimsize(c_list,0)-1][i]=c_list[dimsize(c_list,0)-2][i]
		theta_list[dimsize(theta_list,0)-1][i]=theta_list[dimsize(theta_list,0)-2][i]
		phi_list[dimsize(phi_list,0)-1][i]=phi_list[dimsize(phi_list,0)-2][i]
		namew="Line_"+num2str(i)
		wave linew=root:results:$namew
		infonote=note(linew)
		namew="LineCopy_"+num2str(i)
		make/o/n=(modules,3) root:results:$namew
		wave linetot=root:results:$namew
		LengthLine=linelength(root:results:$("Line_"+num2str(i)))
		Linetot[][0][]=(p-modules/2)*(Length)*sin(theta_list[dimsize(theta_list,0)-1][i])*cos(phi_list[dimsize(phi_list,0)-1][i])+a_list[dimsize(a_list,0)-1][i]+Lengthmat[p][i]*sin(theta_list[dimsize(theta_list,0)-1][i])*cos(phi_list[dimsize(phi_list,0)-1][i])
		Linetot[][1][]=(p-modules/2)*(Length)*sin(theta_list[dimsize(theta_list,0)-1][i])*sin(phi_list[dimsize(phi_list,0)-1][i])+b_list[dimsize(b_list,0)-1][i]+Lengthmat[p][i]*sin(theta_list[dimsize(theta_list,0)-1][i])*sin(phi_list[dimsize(phi_list,0)-1][i])
		Linetot[][2][]=(p-modules/2)*(Length)*cos(theta_list[dimsize(theta_list,0)-1][i])+c_list[dimsize(c_list,0)-1][i]+Lengthmat[p][i]*cos(theta_list[dimsize(theta_list,0)-1][i])
	endfor
	for(i=0;i<dimsize(movement,0);i+=1)
		wave w1=root:results:$("line_"+num2str(movement[i][0]))
		wave w2=root:results:$("line_"+num2str(movement[i][2]))
		movement[i][4]=distance(w1,movement[i][1],w2,movement[i][3])
		wave w1=root:results:$("linecopy_"+num2str(movement[i][0]))
		wave w2=root:results:$("linecopy_"+num2str(movement[i][2]))
		criticaldis[i]=distance(w1,movement[i][1],w2,movement[i][3])
	endfor
end
function distance(w1,a,w2,b)
	// returns the ditance between "a" indexed tuple on wave w1 vs the "b" indexed tuple on wave w2
	wave w1,w2
	variable a,b
	return Sqrt(((w1[a][0]-w2[b][0])^2+(w1[a][1]-w2[b][1])^2+(w1[a][2]-w2[b][2])^2))
end
Function DistanceOptimization(w,passwave)
	// functiont that returns an associated error with a given set of dynamical paramers
	// Only modifies angular and positional parameters, does not touch length
	// i.e considers constant lengthwave with variable a,b,c,theta,phi
	wave w
	wave passwave
	nvar n_lines,modules,linksize
	wave criticaldis=root:results:CriticalDis
	wave movement=root:results:movement
	make/o/n=(dimsize(movement,0)) root:results:CompareWave
	wave Comparewave=root:results:CompareWave
	wave deltawave=root:results:deltawave
	CompareWave[]=Movement[p][4]
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:results:theta_list
	wave phi_list=root:results:phi_list
	wave lengthmat=root:results:lengthmat
	variable i,j
	make/o/n=(n_lines) awavecurr
	make/o/n=(n_lines) bwavecurr
	make/o/n=(n_lines) cwavecurr
	make/o/n=(n_lines) thetawavecurr
	make/o/n=(n_lines) phiwavecurr
	make/o/n=(modules,n_lines) lengthmatcurr
	for(i=0;i<N_lines;i+=1)
		awavecurr[i]=passwave[i]
		bwavecurr[i]=passwave[i+n_lines]
		cwavecurr[i]=passwave[i+n_lines*2]
		thetawavecurr[i]=passwave[i+n_lines*3]
		phiwavecurr[i]=passwave[i+n_lines*4]
	endfor
	for(i=0;i<modules;i+=1)
		for(j=0;j<n_lines;j+=1)
			lengthmatcurr[i][j]=passwave[n_lines*(5+i)+j]
		endfor
	endfor
	updateLinesandDis(awavecurr,bwavecurr,cwavecurr,thetawavecurr,phiwavecurr,lengthmatcurr)
	variable energy=0
	variable LengthWeight=1,RepulseWeight=1
	duplicate/o criticaldis,tempwave
	//Term that fixes the crosslinks between timesteps
	tempwave[]=(criticaldis[p]-linksize)^2
	energy=sum(tempwave)
	//for(i=0;i<numpnts(criticaldis);i+=1)
	//	energy+=(criticaldis[i]-LinkSize)^2
	//endfor
	//Term that weights the lengthchanges
	//for(i=0;i<modules;i+=1)
	//	for(j=0;j<n_lines;j+=1)
	//		energy+=LengthWeight*(lengthmat[i][j]-lengthmatcurr[i][j])^2
	//	endfor
	//endfor
	//Term that fixes the repulsion between timesteps
	//for(i=0;i<numpnts(criticaldis);i+=1)
	//	energy+=1/(criticaldis[i])^2
	//endfor
	return energy
end
function buildpasswave()
	// generates the One-D wave from the 5-D set of dynamic parameters a,b,c,theta,phi
	nvar n_lines,modules
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:results:theta_list
	wave phi_list=root:results:phi_list
	wave lengthmat=root:results:lengthmat
	wave criticaldis=root:results:CriticalDis
	wave movement=root:results:movement
	make/o/n=(n_lines*(5+modules)) root:results:Passwave
	wave passwave=root:results:passwave
	variable i,j
	for(i=0;i<n_lines;i+=1)
		passwave[i]=a_list[dimsize(a_list,0)-1][i]
		passwave[i+n_lines]=b_list[dimsize(b_list,0)-1][i]
		passwave[i+n_lines*2]=c_list[dimsize(c_list,0)-1][i]
		passwave[i+n_lines*3]=theta_list[dimsize(theta_list,0)-1][i]+enoise(Pi/8)
		passwave[i+n_lines*4]=phi_list[dimsize(phi_list,0)-1][i]+enoise(Pi/4)
		//passwave[i+n_lines*3]=enoise(Pi/2)+Pi/2
		//passwave[i+n_lines*4]=enoise(Pi)+Pi
	endfor
	for(i=0;i<modules;i+=1)
		for(j=0;j<n_lines;j+=1)
			passwave[n_lines*(5+i)+j]=lengthmat[i][j]
		endfor
	endfor
end
function DoProteinShift()
	// Finds the set of a,b,c,theta,phi parameters that minimize the error
	nvar n_lines,modules
	variable i,j
	string namew
	//updateLengthMat()
	BuildLineCopies()
	buildpasswave()
	wave passwave=root:results:passwave
	DistanceOptimization(passwave,passwave)
	duplicate/O passwave,root:results:x_param_wave
	wave x_Param_wave=root:results:x_param_wave
	Optimize/A=0/X=X_param_wave/Y=1/I=100/M={0,0}/R=X_param_wave/Q DistanceOptimization,passwave
	print "error = "+num2str(V_min)
	//x_param_wave gives set of a,b,c,theta,phi that minimize the error
	DistanceOptimization(passwave,x_param_wave)
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:results:theta_list
	wave phi_list=root:results:phi_list
	wave lengthmat=root:results:lengthmat
	wave deltawave=root:results:deltawave
	for(i=0;i<N_lines;i+=1)
		a_list[dimsize(a_list,0)-1][i]=x_param_wave[i]
		b_list[dimsize(b_list,0)-1][i]=x_param_wave[i+n_lines]
		c_list[dimsize(c_list,0)-1][i]=x_param_wave[i+n_lines*2]
		theta_list[dimsize(theta_list,0)-1][i]=x_param_wave[i+n_lines*3]
		phi_list[dimsize(phi_list,0)-1][i]=x_param_wave[i+n_lines*4]
	endfor
	for(i=0;i<modules;i+=1)
		for(j=0;j<n_lines;j+=1)
			lengthmat[i][j]=x_param_wave[n_lines*(5+i)+j]
		endfor
	endfor
	updatelines()
end
function updateLengthMat()
	// Describes the change in length for each domain on the protein
	//Note that the lengthmat is actually the change between two timepoints
	wave lengthmat=root:results:lengthmat
	wave display_lengthmat=root:results:display_lengthmat
	wave unfoldmat=root:results:unfoldmat
	variable shift,j,k,domain_no,i,nchange
	nvar n_lines,modules,displayw
	wave deltaLwave=root:results:deltawave
	wave networkmat=root:results:networkmat
	string exec
	lengthmat=display_lengthmat
	display_lengthmat=0
	//matricies are labeled nameMat[domain_no][module_no]([j][i])
	for(i=0;i<n_lines;i+=1)//i is domain number
		wave sizew=root:results:$("Sizew_"+num2str(i))
		wave nwave=root:results:$("N_"+num2str(i))
		if ((nwave[numpnts(nwave)-1]>nwave[numpnts(nwave)-2]))
			nchange=nwave[numpnts(nwave)-1]-nwave[numpnts(nwave)-2]
			for(k=0;k<nchange;k+=1)
				Domain_no=getRandomPoint(i,0)
				unfoldmat[domain_no][i]=deltaLwave[i]/nchange
			endfor
		elseif ((nwave[numpnts(nwave)-1]<nwave[numpnts(nwave)-2]))
			nchange=nwave[numpnts(nwave)-2]-nwave[numpnts(nwave)-1]
			for(k=0;k<nchange;k+=1)
				Domain_no=getRandomPoint(i,1)
				unfoldmat[domain_no][i]=0
			endfor
		endif
		for(j=0;j<modules;j+=1)
			shift=unfoldmat[j][i]
			for(k=0;k<modules;k+=1)
				if (k<j)
					display_lengthmat[k][i]+=-shift/2*(1-krondelta(j,k))
				else
					display_lengthmat[k][i]+=shift/2*(1-krondelta(j,k))
				endif
			endfor
			sizew[j][]=0.01*(1-krondelta(0,shift))+0.2*krondelta(0,shift)
		endfor
		if (displayw)
			exec ="ModifyGizmo ModifyObject=scatter"+num2str(i)+", property={sizetype,1}"
			Execute exec
			exec="ModifyGizmo ModifyObject=scatter"+num2str(i)+", property={sizewave,root:RESULTS:" + nameofwave(SizeW) + "}"
			Execute exec
		endif
	endfor
	matrixop/O lengthmat=display_lengthmat-lengthmat
end
function UpdateLines()
	// transfers over line copies to actual line wave
	variable i,j
	nvar n_lines,showconnect,modules
	string exec
	wave movement=root:results:movement
	wave lengthmat=root:results:lengthmat
	wave comparemat=root:results:comparemat
	Make/O/N=(0,3) root:RESULTS:Crosslinks
	WAVE Crosslinks=root:RESULTS:Crosslinks
	for(i=0;i<n_lines;i+=1)
		wave w1=root:results:$("Linecopy_"+num2str(i))
		wave w2=root:Results:$("Line_"+num2str(i))
		w2=w1
	endfor
	for(i=0;i<dimsize(movement,0);i+=1)
		wave w1=root:results:$("line_"+num2str(movement[i][0]))
		wave w2=root:results:$("line_"+num2str(movement[i][2]))
		InsertPoints/M=0 dimsize(Crosslinks,0),2,Crosslinks
		Crosslinks[dimsize(Crosslinks,0)-2][0]=w1[movement[i][1]][0]
		Crosslinks[dimsize(Crosslinks,0)-2][1]=w1[movement[i][1]][1]
		Crosslinks[dimsize(Crosslinks,0)-2][2]=w1[movement[i][1]][2]
		Crosslinks[dimsize(Crosslinks,0)-1][0]=w2[movement[i][3]][0]
		Crosslinks[dimsize(Crosslinks,0)-1][1]=w2[movement[i][3]][1]
		Crosslinks[dimsize(Crosslinks,0)-1][2]=w2[movement[i][3]][2]
		InsertPoints/M=0 dimsize(Crosslinks,0),1,Crosslinks
		Crosslinks[dimsize(Crosslinks,0)-1][]=NAN
		movement[i][4]=distance(w1,movement[i][1],w2,movement[i][3])
	endfor
	if(showconnect)
		exec = "ModifyGizmo setDisplayList=" + num2str(2*N_lines-1)+ ", object=path"+ num2str(N_lines)
		Execute exec

		exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:Crosslinks"
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={pathColor,0,0,0,1}"//////65535
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={lineWidthType,1}"
		Execute exec

		exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={linewidth,5}"
		Execute exec
	endif
end
function getrandompoint(line_no,bool)//bool is zero for folded->unfolded and 1 otherwise
	// function that returns a random number for a domain to fold/unfold
	variable line_no,bool
	variable j,rand
	nvar modules
	wave unfoldmat=root:Results:unfoldmat
	make/o/n=0 dummywave
	if (bool==0)
		for(j=0;j<modules;j+=1)
			if (unfoldmat[j][line_no]==0)
				insertpoints numpnts(dummywave),1,dummywave
				dummywave[numpnts(dummywave)-1]=j
			endif
		endfor
	else
		for(j=0;j<modules;j+=1)
			if (unfoldmat[j][line_no]!=0)
				insertpoints numpnts(dummywave),1,dummywave
				dummywave[numpnts(dummywave)-1]=j
			endif
		endfor
	endif
	if (numpnts(dummywave)==0)
		print "Wtf"
	endif
	rand= abs(round((enoise(0.5)+0.5)*(numpnts(dummywave)-1)))
	return dummywave[rand]
end
function KronDelta(i,j)
	// simply a kronecker delta
	variable i,j
	if(i==j)
		return 1
	else
		return 0
	endif
end
function updateLinesandDis(a_list,b_list,c_list,theta_list,phi_list,lengthmat)
	// updates the line copies given the current LengthMat and set of dynamic parameters
	wave ,,a_list,b_list,c_list,theta_list,phi_list,lengthmat
	variable i
	nvar n_lines,modules,length
	string namew
	wave criticaldis=root:results:$("CriticalDis")
	wave movement=root:results:movement
	//wave lengthmat=root:results:lengthmat
	for(i=0;i<n_lines;i+=1)
		namew="LineCopy_"+num2str(i)
		make/o/n=(modules,3) root:results:$namew
		wave linetot=root:results:$namew
		Linetot[][0][]=(p-modules/2)*(Length)*sin(theta_list[i])*cos(phi_list[i])+a_list[i]+Lengthmat[p][i]*sin(theta_list[i])*cos(phi_list[i])
		Linetot[][1][]=(p-modules/2)*(Length)*sin(theta_list[i])*sin(phi_list[i])+b_list[i]+Lengthmat[p][i]*sin(theta_list[i])*sin(phi_list[i])
		Linetot[][2][]=(p-modules/2)*(Length)*cos(theta_list[i])+c_list[i]+Lengthmat[p][i]*cos(theta_list[i])
	endfor
	for(i=0;i<dimsize(movement,0);i+=1)
		wave w1=root:results:$("linecopy_"+num2str(movement[i][0]))
		wave w2=root:results:$("linecopy_"+num2str(movement[i][2]))
		criticaldis[i]=distance(w1,movement[i][1],w2,movement[i][3])
	endfor
end
function Saveit(path,name)
	// function to save a snap shot of the gizmo
	// need to change the path destination if not working on Kirill's MacBook Pro
	string name,path
	nvar displayw
	svar crosslinkingpath,unfoldingpath
	if (displayw)
		if (stringmatch(path,"Crosslinking"))
			SavePICT/O/E=-6/RES=2400/p=Crosslinkingpath as name+".png"
		elseif (stringmatch(path,"Unfolding"))
			SavePICT/O/E=-6/RES=2400/p=unfoldingpath as name+".png"
		endif
	endif
end
function updateColors()
	// updates the colors of the lines given their angle relative to the z-axis
	wave theta_list=root:results:theta_list
	string namew,exec,infonote
	variable i,j,force
	nvar n_lines,modules,displayw
	if (displayw)
		for(i=0;i<N_lines;i+=1)
			force=abs(cos(theta_list[dimsize(theta_list,0)-1][i]))
			ColorTab2Wave Rainbow
			WAVE M_Colors
			Variable pick=((dimsize(M_Colors,0)-1)*(force))
			exec ="ModifyGizmo ModifyObject=scatter"+ num2str(i)+", property={color,"+ num2str(M_Colors[pick][0]/65535) + ","+ num2str(M_Colors[pick][1]/65535) + ","+num2str(M_Colors[pick][2]/65535) + ",0}"
			Execute exec
		endfor
	endif
end
function ztothists()
	wave zwave=root:results:zdistmat
	variable i,j,k,l
	nvar n_lines, modules
	for(i=0;i<(dimsize(zwave,0));i+=100)
		make/o/n=(n_lines*modules) $("ZDist_"+num2str(i))
		wave curr=$("ZDist_"+num2str(i))
		for(j=0;j<(n_lines*modules);j+=1)
			curr[j]=zwave[i][j]
		endfor
	endfor
end
function GelLengthHist()
	wave zwave=root:results:zdistmat
	variable i,j,k,l,r,q,A,s,t,binsize
	nvar n_lines, modules,zby,conc_molec
	variable zbox=zby/(conc_molec)^(1/3)
	make/o/n=(dimsize(zwave,0)-1) root:results:Gel_Length
	make/o/n=(dimsize(zwave,0)-1) root:results:s_Gel_Length
	make/o/n=(dimsize(zwave,0)-1) A_wave
	wave len=root:results:Gel_Length
	wave s_len=root:results:s_Gel_Length
	r=zbox
	q=2
	s=0
	t=N_lines/5
	Make/O/T/N=1 T_Constraints
	binsize=4//11*exp(-0.015*N_lines)//Value-Tested
	T_Constraints[0] = {"K3 > 0","K4 < 100","K4 > 1"}
	for(i=0;i<(dimsize(zwave,0)-1);i+=1)
		make/o/n=(dimsize(zwave,1)) newz
		newz[]=zwave[i][p]
		Make/N=100/O newz_hist
		wavestats/Q newz
		l=v_min-20
		Histogram/C/N/B={l,binsize,100} newz,newz_hist
		wavestats/q newz_hist
		T_Constraints[0] = {"K3 > "+num2str(r-zbox/10),"K3 < "+num2str(r+zbox/10),"K4 < 100","K4 > 1"}//,"K3 < "+num2str(v_maxloc-v_minloc)}
		Make/D/N=5/O W_coef
		W_coef[0] = {0,t,s,r,q}
		FuncFit/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
		//curvefit/x=1/Q gauss newz_hist /D
		wave w_coef,w_sigma
		len[i]=abs(w_coef[3])
		s_len[i]=abs(w_sigma[3])
		A_wave[i]=w_coef[4]
		r=abs(w_coef[3])
		q=abs(w_coef[4])
		if (q>20)
			q=5
		endif
		t=abs(w_coef[1])
		s=w_coef[2]
		print num2str((i/(dimsize(zwave,0))*100))+"% Done"
	endfor
	//drops(len)
	avgsFixes()
end
function Calc_GelLength()
	//Need to create Z-wave still
	variable i,j,k,l,r,q,A,s,t,binsize,u
	nvar n_lines, modules,zby,conc_molec,const_0,const_1,const_2,const_3
	wave len = root:results:Gel_Length
	wave s_len = root:results:s_Gel_Length
	wave A_wave
	//Inc is the the increment of the projection segmentation per line
	variable Inc=0.01
	variable zbox=zby/(conc_molec)^(1/3)
	insertpoints numpnts(len),1,len,s_len,A_wave
	make/o/n=0 newz
	for(i=0;i<n_lines;i+=1)
		wave linew=$("Root:results:Line_"+num2str(i))
		for(j=min(linew[0][2],linew[modules-1][2]);j<max(linew[0][2],linew[modules-1][2]);j+=inc)
			insertpoints numpnts(newz),1,newz
			newz[numpnts(newz)-1]=j
		endfor
	endfor
	Make/O/T/N=1 T_Constraints
	Make/N=100/O newz_hist
	wavestats/Q newz
	l=v_min-40
	//binsize=4//11*exp(-0.015*N_lines)//Value-Tested
	binsize=3.49*V_sdev*(numpnts(newz))^(-1/3)
	u=(v_max+40-l)/binsize
	Histogram/C/N/B={l,binsize,u} newz,newz_hist
	wavestats/Q newz_hist
	T_Constraints[0] = {"K3 > "+num2str(const_0-zbox/10),"K3 < "+num2str(const_0+zbox/10),"K4 < 100","K4 > 1"}//,"K3 < "+num2str(v_maxloc-v_minloc)}
	Make/D/N=5/O W_coef
	W_coef[0] = {0,const_3,const_2,const_0,const_1}
	FuncFit/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
	//curvefit/x=1/Q gauss newz_hist /D
	wave w_coef,w_sigma
	len[numpnts(len)-1]=abs(w_coef[3])
	s_len[numpnts(len)-1]=abs(w_sigma[3])
	A_wave[numpnts(len)-1]=w_coef[4]
	const_0=abs(w_coef[3])
	const_1=abs(w_coef[4])
	if (const_1>20)
		const_1=5
	endif
	const_3=abs(w_coef[1])
	const_2=w_coef[2]
end
function Calc_GelLength_ionel()
	//Need to create Z-wave still
	variable i,j,k,l,r,q,A,s,t,binsize,u
	nvar n_lines, modules,zby,conc_molec,xby,yby
	wave len = root:results:Gel_Length
	variable avgnum=xby*yby
	insertpoints numpnts(len),1,len
	//Inc is the the increment of the projection segmentation per line
	make/o/n=(modules*n_lines) zwave
	k=0
	for(i=0;i<N_lines;i+=1)
		wave linew=$("root:results:line_"+num2str(i))
		for(j=0;j<modules;j+=1)
			zwave[k]+=linew[j][2]
			k+=1
		endfor
	endfor
	sort zwave,zwave
	variable maxavg=0,minavg=0
	for(i=0;i<avgnum;i+=1)
		maxavg+=zwave[numpnts(zwave)-i]
		minavg+=zwave[i]
	endfor
	print (maxavg-minavg)/avgnum
	len[numpnts(len)-1]=(maxavg-minavg)/avgnum
end
function avgsFixes()
	nvar n_lines
	variable i,j,k
	wave duration = root:sim:duration_list
	wave gel_length=root:results:Gel_Length
	make/o/n=(numpnts(root:results:x_0)) root:results:x_avg
	wave x_avg= root:results:x_avg
	x_avg=0
	for(i=0;i<n_lines;i+=1)
		wave xwave=$("root:results:x_"+num2str(i))
		wave f=$("root:results:F_"+num2str(i))
		if (i==0)
			duplicate/o F,Root:results:F_avg
			wave favg=Root:results:F_avg
		else
			favg+=F
		endif
		setscale/I x,0,(sum(duration)),"s",xwave
		x_avg[]+=xwave[p]
	endfor
	setscale/I x,0,(sum(duration)),"s",gel_length,x_avg
end
function drops(w)
	wave w
	variable i,j,k,r,flag,mask,avg,low
	nvar n_lines, modules,zby,conc_molec
	variable zbox=zby/(conc_molec)^(1/3)
	variable cut=zbox/5
	variable pointnum=100
	variable threshhold=0.5
	variable numflag,startp,endp
	i=0
	flag=1
	do
		if (i==8760)
			print ""
		endif
		startp=i
		r=0
		do
			i+=1
			endp=i
			r+=1
			if (r==200)
				break
			endif
		while ((round(w[startp])-round(w[endp]))>threshhold)
		mask=1
		low=1e10
		for(j=startp;j<endp;j+=1)
			if(w[j]>w[startp])
				mask=0
			endif
			if(w[j]<low)
				low=w[j]
			endif
		endfor
		if (mask==1)
			avg=(w[endp]+w[startp])/2
			for(i=startp;i<endp;i+=1)
				w[i]=avg
			endfor
		endif
		if(i==(numpnts(w)-2))
			flag=0
		endif
	while (flag)
end
function zdists()
	wave zdistwave=root:results:zdistmat
	variable i,j,k,r,l
	nvar n_lines,modules
	k=0
	for(i=0;i<n_lines;i+=1)
		wave linew=$("Root:results:Line_"+num2str(i))
		for(j=0;j<modules;j+=1)
			zdistwave[dimsize(zdistwave,0)-1][k]=linew[j][2]
			k+=1
		endfor
	endfor
	insertpoints dimsize(zdistwave,0),1,zdistwave
end
Function SuperGauss(w,x) : FitFunc
	WAVE w;Variable x
	//CurveFitDialog/
	//CurveFitDialog/ Coefficients 5
	//CurveFitDialog/ w[0] = y0
	//CurveFitDialog/ w[1] = A
	//CurveFitDialog/ w[2] = x0
	//CurveFitDialog/ w[3] = Width
	//CurveFitDialog/ w[4] = P
	return w[0]+w[1]*exp(-(((x-w[2])/w[3])^2)^(w[4]))
End
//Function that corrects the crosslink distances
function AppendLinesAndNetwork()
	// Appends an elemt to Linematrix and builds the the passmat, network mat representing the past step
	// Line copies are the ones modulated during the optimzaiton
	updateLengthMat()
	string namew,exec
	variable i,j,k,m,modnum,linenum
	nvar modules,length,n_lines,numberpoints
	wave LineMatrix=root:results:LineMatrix
	wave passmat=root:results:passmat
	wave networkmat=root:results:networkmat
	wave domainconnections=root:Results:domainconnections
	wave pivot=root:results:pivot
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:Results:theta_list
	wave phi_list=root:results:phi_list
	wave lengthmat=root:results:lengthmat
	//shift passmat points according to lengthmat
	//passmat is past iteration of lengthmat plus the shift according to angles set on past interation
	// Shifts are due to slight extensions and unfoldings
	k=0
	redimension/N=(modules*n_lines) domainconnections
	for(i=0;i<numpnts(domainconnections);i+=1)
		if (domainconnections[i]==1)
			modnum=mod(i,modules);linenum=(i-mod(i,modules))/modules
			passmat[k][0]=networkmat[k][0]+Lengthmat[modnum][linenum]*sin(theta_list[dimsize(theta_list,0)-1][linenum])*cos(phi_list[dimsize(phi_list,0)-1][linenum])
			passmat[k][1]=networkmat[k][1]+Lengthmat[modnum][linenum]*sin(theta_list[dimsize(theta_list,0)-1][linenum])*sin(phi_list[dimsize(phi_list,0)-1][linenum])
			passmat[k][2]=networkmat[k][2]+Lengthmat[modnum][linenum]*cos(theta_list[dimsize(theta_list,0)-1][linenum])
			k+=1
		endif
	endfor
	redimension/N=(modules,n_lines) domainconnections
end
function EnergyCalc(w,PassMat)
	//Calculates the strech energy between the new set of coordinates
	//Inputs: PriorMat is the Matrix of crosslinks from the PAST iteration
	// PassMat is the matrix of connections being energy calculated
	//Potential imporvements: Combine the inter- and intra-line connections into
	// a single for-loop
	wave w,PassMat
	wave priorMat=root:results:networkmat
	NVAR N_Lines, modules,linksize,length
	NVAR p_val = root:sim:P_val
	wave domainconnections=root:results:domainconnections
	wave clindex=root:results:clindex
	wave bendindex=root:results:bendindex
	matrixop/o deltamat=PassMat-PriorMat
	variable i,j,energy=0,k

	// Energy PreFactors
	variable StrechPrefactor=2*p_val*4.11/pi/(length/2)^4	//from article and wiki (does not include length divider (1/prior))
	variable BendPrefactor=p_val*4.11/2 // from article and wiki (does not include length divider (1/prior^3))
	variable CrossLinkPrefactor=372030//calculated from OPLS force-field

	//Determine Difference in Streching between crosslinks
	k=0
	for(i=0;i<n_lines;i+=1)
		//CLsum is # of crosslink locations on that line
		matrixop/o CLSum=sum(col(domainconnections,i))
		//Compute Strech Energy between Line Crosslinks
		for(j=k;j<k+CLSum[0]-1;j+=1)
			Matrixop/o DeltaU=row(deltamat,j)-row(deltamat,j+1)
			matrixop/o OldVecLength=powr(sum(powr(row(priormat,j)-row(priormat,j+1),2)),0.5)
			Matrixop/o UnitOld=normalize(row(PriorMat,j)-row(priormat,j+1))
			matrixop/o LinkEnergy=(DeltaU.UnitOld)
			energy+=(StrechPrefactor/oldveclength[0])*(LinkEnergy[0])^2
		endfor
		//Compute Bending Between Crosslinks
		if(CLSUm[0]>2)
			for(j=k;j<k+CLSum[0]-2;j+=1)
				matrixop/O DeltaU=row(deltaMat,j)+row(deltamat,j+1)-2*row(deltamat,j+2)
				matrixop/o OldVecLength=powr(sum(powr(row(priormat,j)-row(priormat,j+1),2)),0.5)
				Matrixop/o UnitOld=normalize(row(PriorMat,j)-row(priormat,j+1))
				redimension/N=3 deltaU,unitold
				cross DeltaU,unitOld
				wave w_cross
				energy+=(BendPrefactor/oldveclength[0]^3)*(w_cross[0]^2+w_cross[1]^2+w_cross[2]^2)
			endfor
		endif
		k+=CLSum[0]
	endfor
	//Calculate crosslink difference in strech and linksize
	for(i=0;i<dimsize(clindex,0);i+=1)
		//Calculate Strech energy between crosslinks
		Matrixop/o DeltaU=row(deltamat,clindex[i][0])-row(deltamat,clindex[i][1])
		matrixop/o OldVecLength=powr(sum(powr(row(priormat,clindex[i][0])-row(priormat,clindex[i][1]),2)),0.5)
		Matrixop/o UnitOld=normalize(row(PriorMat,clindex[i][0])-row(priormat,clindex[i][1]))
		matrixop/o LinkEnergy=(DeltaU.UnitOld)
		energy+=(StrechPrefactor/oldveclength[0])*(LinkEnergy[0])^2
		//Claculate from link distances
		Matrixop/o LengthVector=row(passmat,clindex[i][0])-row(passmat,clindex[i][1]))
		matrixop/o LinkLength=powr(sumsqr(LengthVector),0.5)
		energy+=CrossLinkPrefactor*(LinkLength[0]-LinkSize)^2
	endfor
	for(i=0;i<dimsize(bendindex,0);i+=1)
		//Calculate Bending between crosslinks
		matrixop/O DeltaU=row(deltaMat,bendindex[i][2])+row(deltamat,bendindex[i][0])-2*row(deltamat,bendindex[i][1])
		matrixop/o OldVecLength=powr(sum(powr(row(priormat,bendindex[i][0])-row(priormat,bendindex[i][1]),2)),0.5)
		Matrixop/o UnitOld=normalize(row(PriorMat,bendindex[i][0])-row(priormat,bendindex[i][1]))
		redimension/N=3 deltaU,unitold
		cross DeltaU,unitOld
		wave w_cross
		energy+=(BendPrefactor/oldveclength[0]^3)*(w_cross[0]^2+w_cross[1]^2+w_cross[2]^2)
	endfor
	return energy
end
function UpdateVariables()
	//This function updates all line parameter waves and updates the networkmat
	//STILL NEED TO DO THE CHANGE IN LENGTH
	nvar modules,length,n_lines,numberpoints
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:results:theta_list
	wave phi_list=root:results:phi_list
	wave passmat=root:results:passmat
	wave networkMat=root:results:networkMat
	wave domainconnections=root:results:domainconnections
	wave linematrix=root:results:linematrix
	wave lengthmat=root:Results:lengthmat
	wave pivot=root:results:pivot
	variable i,j,k,deltheta,delphi,dela,delb,delc
	//pivot passmat to zero
	matrixop/o ChangeMat=PassMat-NetWorkMat
	matrixop/o avgcol=averagecols(ChangeMat)
	matrixop/o domaincols=sumcols(domainconnections)
	k=0
	//update all the variabnle lists
	InsertPoints/M=0 dimsize(a_list,0), 1, a_list, b_list,c_list,theta_list,phi_list
	for(i=0;i<n_lines;i+=1)
		deltheta=0;delphi=0;dela=0;delb=0;delc=0;
		for(j=k;j<k+domaincols[0][i];j+=1)
			// 'del{x}' accounts for center of mass changes to that protien
			dela+=Changemat[j][0];delb+=changemat[j][1];delc+=changemat[j][2]
			if ((domaincols[0][i]>1) && (j!=(k+domaincols[0][i]-1)))
				// A single connection protein should not rotate
				matrixop/o linknet=row(networkmat,j)-row(networkmat,j+1)
				matrixop/o linkpass=row(passmat,j)-row(passmat,j+1)
				deltheta+=(atan(sqrt((linkpass[0][0])^2 + (linkpass[0][1])^2)/linkpass[0][2])-atan(sqrt((linknet[0][0])^2 + (linknet[0][1])^2)/linknet[0][2]))/(domaincols[0][i]-1)
				delphi+=(atan(linkpass[0][1]/linkpass[0][0])-atan(linknet[0][1]/linknet[0][0]))/(domaincols[0][i]-1)
			endif
		endfor
		a_list[dimsize(a_list,0)-1][i]+=a_list[dimsize(a_list,0)-2][i]+dela
		b_list[dimsize(b_list,0)-1][i]+=b_list[dimsize(a_list,0)-2][i]+delb
		c_list[dimsize(c_list,0)-1][i]+=c_list[dimsize(a_list,0)-2][i]+delc
		theta_list[dimsize(theta_list,0)-1][i]=theta_list[dimsize(theta_list,0)-2][i]+deltheta
		phi_list[dimsize(phi_list,0)-1][i]=phi_list[dimsize(phi_list,0)-2][i]+delphi
		k+=domaincols[0][i]
	endfor
	//Current Iteration becomes optimized iteration
	networkmat=passmat
end
function updateGraph()
	nvar N_lines,modules,length
	variable r,xcm=0,ycm=0,zcm=0,i
	wave linematrix=root:results:linematrix
	wave lengthmat=root:Results:display_lengthmat
	wave a_list=root:results:a_list
	wave b_list=root:results:b_list
	wave c_list=root:results:c_list
	wave theta_list=root:results:theta_list
	wave phi_list=root:results:phi_list
	//use variables from past iteration to construct lines for updates
	for(r=0;r<n_lines;r+=1)
		wave linew=root:results:$("Line_"+num2str(r))
		linew[][0]=(p-modules/2)*(Length)*sin(theta_list[dimsize(theta_list,0)-1][r])*cos(phi_list[dimsize(phi_list,0)-1][r])+a_list[dimsize(a_list,0)-1][r]+Lengthmat[p][r]*sin(theta_list[dimsize(theta_list,0)-1][r])*cos(phi_list[dimsize(phi_list,0)-1][r])
		linew[][1]=(p-modules/2)*(Length)*sin(theta_list[dimsize(theta_list,0)-1][r])*sin(phi_list[dimsize(phi_list,0)-1][r])+b_list[dimsize(b_list,0)-1][r]+Lengthmat[p][r]*sin(theta_list[dimsize(theta_list,0)-1][r])*sin(phi_list[dimsize(phi_list,0)-1][r])
		linew[][2]=(p-modules/2)*(Length)*cos(theta_list[dimsize(theta_list,0)-1][r])+c_list[dimsize(c_list,0)-1][r]+Lengthmat[p][r]*cos(theta_list[dimsize(theta_list,0)-1][r])
		for(i=0;i<modules;i+=1)
			xcm+=linew[i][0]/(n_lines*modules)
			ycm+=linew[i][1]/(n_lines*modules)
			zcm+=linew[i][2]/(n_lines*modules)
		endfor
	endfor
	for(r=0;r<n_lines;r+=1)
		wave linew=root:results:$("Line_"+num2str(r))
		wave linehis=root:results:$("LineHistory_"+num2str(r))
		linew[][0]-=xcm;linew[][1]-=ycm;linew[][2]-=zcm
		linehis[][][dimsize(linehis,2)-1]=linew[p][q]
		insertpoints/M=2 dimsize(linehis,2),1,linehis
	endfor
end
Function TheGel()
	// The function that runs the gel
	///CHECK IF CLINDEX FUCKING WORKS
	NVAR dt = root:SIM:dt
	NVAR D	=	root:SIM:D_coef
	WAVE	Force_List		=	root:SIM:Force_List
	WAVE	Duration_List		=	root:SIM:Duration_List
	WAVE	Ramp_List		=	root:SIM:Ramp_List
	WAVE 	PMFM = root:PMFM
	Wave 	CommandF = root:SIM:CommandF //workd with dt = 1e-7 and jmax = 1e3
	Wave 	MinXYF = root:MinXYF
	Wave 	F_table = root:F_table
	WAVE deltawave = root:results:deltawave
	WAVE theta_list = root:results:theta_list
	wave movement=root:results:movement
	wave lengthmat=root:Results:lengthmat
	duplicate/o lengthmat,root:results:display_lengthmat
	lengthmat=0
	NVAR N_lines = root:N_lines
	svar UnfoldingPath,OptimizationPath,mainpath
	nvar displayw
	nvar modules,zby,conc_molec
	make/o/n=(modules,n_lines) root:results:unfoldmat_history=0
	Variable  kT=4.11
	Variable Duration=dt*1e3*numpnts(CommandF)
	Variable i, j,k,l ,xc=3, tc=0,Fr=0, F_b=0, current_E_var, currentF, oldF
	variable updateflag, numupdate
	NVAR 	N_modules = root:SIM:N_modules
	Duplicate/O CommandF root:Sim:x_wave, root:Sim:F_wave,  root:Sim:N_wave
	Wave x_wave = root:Sim:x_wave
	Wave F_wave = root:Sim:F_wave
	Wave N_wave = root:Sim:N_wave
	x_wave=NaN; F_wave=NaN; N_wave=NaN
	Variable D0=15000,f_factor,xcold
	String nameWx, nameWf, nameWn,namehis
	// some gel length variables
	variable/G Const_0=zby/(conc_molec)^(1/3)
	variable/G Const_1=2
	variable/G Const_2=0
	variable/G Const_3=N_lines/5
	//convert all lines into into a matrix
	//Note, Linematrix is not the same points that are used in display purposes.
	make/o/n=(modules,3,n_lines) root:results:linematrix=0
	make/o/n=3 root:results:pivot
	wave linematrix=root:results:linematrix
	wave pivot=root:results:pivot
	for(i=0;i<n_lines;i+=1)
		wave linew=root:results:$("Line_"+num2str(i))
		lineMatrix[][][i]=linew[p][q]
	endfor
	pivot[]=linematrix[0][p][0]
	//Convert movement matrix into a connection matrix for operations in optimization
	//Matrix value is 1 is connection exists at that domain, otherwise is zero
	//Define the Domainconnections mat
	make/o/n=(modules,n_lines) root:results:DomainConnections=0 //rows label modules, columns label lines
	wave domainconnections=root:results:domainconnections
	make/o/n=(0,2) root:results:CLindex_iter
	wave clindex_iter = root:results:CLindex_iter
	for(i=0;i<dimsize(movement,0);i+=1)
		if (DomainConnections[movement[i][1]][movement[i][0]]==0)
			DomainConnections[movement[i][1]][movement[i][0]]=1
		endif
		if (DomainConnections[movement[i][3]][movement[i][2]]==0)
			DomainConnections[movement[i][3]][movement[i][2]]=1
		endif
		insertpoints/M=0 dimsize(clindex_iter,0),1,clindex_iter
		clindex_iter[dimsize(clindex_iter,0)-1][0]=movement[i][1]*(movement[i][0]+1)
		clindex_iter[dimsize(clindex_iter,0)-1][1]=movement[i][3]*(movement[i][2]+1)
	endfor
	//Create Network mat, this is the one the specifies the dynamics of the system
	//The function Make_NetworkMat does this for arbitrary time
	//Networkmat is the previous iteration of the network (priormat)
	make/o/n=(sum(domainconnections),3) root:results:NetworkMat,root:results:PassMat
	wave networkMat=root:results:NetworkMat
	wave passmat=root:results:passmat
	k=0
	for(j=0;j<n_lines;j+=1)
		for(i=0;i<modules;i+=1)
			if(domainconnections[i][j]==1)
				networkMat[k][]=linematrix[i][q][j]
				k+=1
			endif
		endfor
	endfor
	//Make the crosslink_iter mat, values correspond to where connections are
	//unpack domainconnections to allow for indexing
	redimension/N=(modules*n_lines) domainconnections
	make/o/n=(0,2) root:Results:clindex
	wave clindex=root:results:clindex
	//create temp_movement wave for searching
	duplicate/o movement,temp_movement
	for(j=0;j<n_lines;j+=1)
		for(i=0;i<modules;i+=1)
			for(k=0;k<dimsize(temp_movement,0);k+=1)
				//domain must be in one of the two places
				if (((temp_movement[k][0]==j) && (temp_movement[k][1]==i)) || ((temp_movement[k][2]==j) && (temp_movement[k][3]==i)))
					insertpoints/M=0 dimsize(clindex,0),1,clindex
					// the 'l-value' is: (Line_#)*modules+(Module_#)
					l=modules*temp_movement[k][0]+temp_movement[k][1]
					//Sum over domain connections specifies index placement, subtract one for 0-indexing
					clindex[dimsize(clindex,0)-1][0]=sum(domainconnections,0,l)-1
					l=modules*temp_movement[k][2]+temp_movement[k][3]
					clindex[dimsize(clindex,0)-1][1]=sum(domainconnections,0,l)-1
					//delete used point in temp_movement to prevent double-counting
					deletepoints k,1,temp_movement
				endif
			endfor
		endfor
	endfor
	// Make Bending wave
	make/o/n=(0,3) root:results:bendindex
	wave bendindex=root:Results:bendindex
	variable smallindex,largeindex
	make/o/n=3 tempwave
	for(i=0;i<dimsize(clindex,0);i+=1)
		smallindex=min(clindex[i][0],clindex[i][1]);largeindex=max(clindex[i][0],clindex[i][1])
		tempwave[0]=smallindex;tempwave[1]=largeindex
		for(j=i+1;j<dimsize(clindex,0);j+=1)
			if (smallindex==clindex[j][0])
				tempwave[2]=clindex[j][1];sort tempwave,tempwave
				insertpoints/M=0 dimsize(bendindex,0),1,bendindex
				bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
			elseif (smallindex==clindex[j][1])
				tempwave[2]=clindex[j][0];sort tempwave,tempwave
				insertpoints/M=0 dimsize(bendindex,0),1,bendindex
				bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
			elseif (largeindex==clindex[j][0])
				tempwave[2]=clindex[j][1];sort tempwave,tempwave
				insertpoints/M=0 dimsize(bendindex,0),1,bendindex
				bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
			elseif (largeindex==clindex[j][1])
				tempwave[2]=clindex[j][0];sort tempwave,tempwave
				insertpoints/M=0 dimsize(bendindex,0),1,bendindex
				bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
			endif
		endfor
	endfor
	//Make some of the gel length waves
	make/o/n=0 root:results:Gel_Length
	make/o/n=0 root:results:s_Gel_Length
	make/o/n=0 A_wave
	//kill unused wave and redimension domainconnections
	killwaves/Z temp_movement,tempwave
	redimension/N=(modules,n_lines) domainconnections
	//removecrosslinks()
	Differentiate PMFM /D=D_PMFM
	Make/O/N=(dimSize(MINXYF,0)) OneD
	string nameflist
	for(l=0;l<N_lines;l+=1)
		nameWx="x_"+num2str(l)
		nameWf="F_"+num2str(l)
		nameWn="N_"+num2str(l)
		namehis="LineHistory_"+num2str(l)
		nameflist="F_command_"+num2str(l)
		Make/O/N=1 root:results:$nameWx, root:results:$nameWf, root:results:$nameWn
		make/o/n=((modules),3,1) root:results:$namehis
		wave linehis=$("root:results:LineHistory_"+num2str(l))
		wave linew=$("root:results:line_"+num2str(l))
		WAVE Wx=$("root:results:x_"+num2str(l))
		WAVE Wf=$("root:results:F_"+num2str(l))
		WAVE Wn=$("root:results:N_"+num2str(l))
		wx=0
		wf=0
		wn=0
		linehis[][][0]=linew[p][q]
		insertpoints/M=2 dimsize(linehis,2),1,linehis
		Setscale/P x,0,dt, "s", Wx, Wf, Wn
	endfor
	Make/O/N=(N_lines) theta_curr
	Make/O/N=(N_lines) theta_past
	make/o/n=(1,n_lines*modules) root:results:zdistmat
	//first column counts iteration in current state, second denotes last state change
	// 1 means last change was an unfolding, 0 means last change was a refolding
	make/o/n=(N_lines,2) state=0;state[][1]=1
	zdists()
	make/o/n=1 status0=0
	NVAR p_val = root:sim:P_val,length
	make/o/n=3 prefactors
	prefactors[0]=2*p_val*4.11/pi/(length/2)^4
	prefactors[1]=5*p_val*4.11/2
	prefactors[2]=372030
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:Results:clIndex as "clindex.txt"
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:Results:bendindex as "bendindex.txt"
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:Results:domainconnections as "domcon.txt"
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath status0 as "status.txt"
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath prefactors as "prefactors.txt"
	Print "Working on Equilibration..."
	Appendlinesandnetwork()
	runOptimization()
	wave test0
	passmat=test0
	UpdateVariables()
	updateGraph()
	updateColors()
	Saveit("Unfolding","Time_0")
	for(i=0;i<Duration/(dt*1e3)-1;i+=1)
		updateflag=0
		theta_curr[]=theta_list[dimsize(theta_list,0)-1][p]
		if (i>0)
			theta_past[]=theta_list[dimsize(theta_list,0)-2][p]
		endif
		for(l=0;l<N_lines;l+=1)
			WAVE Wx=$("root:results:x_"+num2str(l))
			WAVE Wf=$("root:results:F_"+num2str(l))
			WAVE Wn=$("root:results:N_"+num2str(l))
			wave vlevel=$("root:results:vlevel_"+num2str(l))
			if (i>0)
				currentF=round(commandF[i]*abs(cos(theta_curr[l])))
				oldF=round(commandF[i-1]*abs(cos(theta_past[l])))
			else
				currentF=round(commandF[i]*abs(cos(theta_curr[l])))
			endif
			if (currentF<2)
				currentF=2
			endif
			if (oldF<2)
				OlDF=2
			endif
			if(i>0)
				if (oldF!=CurrentF)
					OneD[]=MinXYF[p][0](oldF)
					xc=wx[i]
					findlevel/Q/P/R=[0,dimSize(MINXYF,0)] OneD,xc// which E-curve were you on at i-1?
					if(V_flag)	//	sometimes at full unfolding/folding findlevel gives NAN
						Wavestats/Q OneD
						if(xc>V_max)
							V_LevelX = N_modules +1
						endif
						if (xc<V_min)
							V_LevelX = 1
						endif
					endif
					OneD[]=MinXYF[p][0](currentF)
					xc=OneD[round(V_LevelX)]// move to that E-curve at the new force
				else
					OneD[]=MinXYF[p][0](currentF)
					xc=wx[i]
				endif
			else
				OneD[]=MinXYF[p][0](currentF)
				xc=MinXYF[1][0](currentF)
			endif
			if(currentF<20)
				D=18450
			else
				D=D0
			endif
			if(i>0)
				xcold=Wx[numpnts(Wx)-1]
				if (xcold<deltax(PMFM))
					xcold=deltax(PMFM)
				endif
				for(j=0;j<1e3;j+=1)
					do
						xc=xcold
						tc+=dt
						F_b = D_PMFM(xc)(currentF)
						Fr = sqrt(2*D*dt)*gnoise(1)
						xc +=  -dt*D * F_b/kT 	+ Fr
					while ((xc<deltax(PMFM)))
					//Correcting for tunneling at large x
					if (xc!=xc)
						xc=xcold
					endif
					xcold=xc
				endfor
			endif
			Findlevel/Q/P/R=[0,dimSize(MINXYF,0)] OneD,xc
			if(V_flag)	//	sometimes at full unfolding/folding findlevel gives NAN
				V_LevelX = Wn[i]+1
			endif
			current_E_var = round(V_LevelX)-1
			if (current_E_var<0)
				current_E_var=0
			endif
			if (i>0)
				state[l][0]+=1
				//if state changes reset time and note direction in state[l][1]
				if (current_e_var>wn[numpnts(wn)-1])
					//this is trying to unfold scenario
					if ((state[l][1]==0) && (state[l][0]>=999))
						// trying to unfold, post refolding, after at least 10 ms
						state[l][0]=0;state[l][0]=1
					elseif ((state[l][1]==0) && (state[l][0]<999))
						// if not enough time has passed for an unfolding, reset state back to
						// last local minimum and adjust current_e_var
						current_e_var=wn[numpnts(wn)-1]
						OneD[]=MinXYF[p][0](currentF)
						xc=OneD[current_e_var+1]
					elseif (state[l][1]==1)
						//if it previously unfolding and trying to unfold again, let it. and reset
						//state[l][0] to zero
						state[l][0]=0
					endif
				endif
				if (current_e_var<wn[numpnts(wn)-1])
					//this is trying to refold now
					if ((state[l][1]==1) && state[l][0]>=999)
						//refolding after at least 10 ms of being unfolded
						state[l][0]=0;state[l][1]=0
					elseif ((state[l][1]==1) && (state[l][0]<999))
						// if not enough time has passed for an refolding, reset state back to
						// last local minimum and adjust current_e_var
						current_e_var=wn[numpnts(wn)-1]
						OneD[]=MinXYF[p][0](currentF)
						xc=OneD[current_e_var+1]
					elseif (state[l][1]==0)
						//refolding after past refolding is fine, no need to change state[l][1]
						state[l][0]=0
					endif
				endif
			endif
			InsertPoints numpnts(Wx),1, Wx, Wf, Wn
			Wx[numpnts(Wx)-1]=xc
			Wf[numpnts(Wf)-1]=currentF+Fr
			Wn[numpnts(Wn)-1]=round(current_E_var)
			//If there is an unfolding, run optimization
			if (i>0)
				if (Wn[numpnts(Wn)-1]!=Wn[numpnts(Wn)-2])
					updateflag=1
				endif
			endif
			if(numpnts(Wx)>1)
				deltawave[l]=Wx[numpnts(Wx)-1]-Wx[numpnts(Wx)-2]
			else
				deltawave[l]=Wx[numpnts(Wx)-1]
			endif
		endfor
		AppendLinesAndNetwork()
		if (updateflag)
			Print "Working on Optimization..."
			sleep/s 0.01
			runOptimization()
			wave test0
			passmat=test0
			//Optimize/A=0/X=passmat/I=15/M={0,0}/T={0.01,0.001}/Q EnergyCalc,passmat
		endif
		UpdateVariables()
		updateGraph()
		updateColors()
		updateUnfoldHistory()
		print num2str(100*i/(Duration/(dt*1e3)-1))+" % Done"
		if ((mod(i,100)==0) && Displayw)
			Saveit("Unfolding","Time_"+num2str(i+1))
			sleep/s 0.01
		endif
		//calc_gellength_ionel()
	endfor
	avgsFixes()
	print "Finished"
End
Function/S UnixShellCommand()
	// Paths must be POSIX paths (using /).
	// Paths containing spaces or other nonstandard characters
	// must be single-quoted. See Apple Techical Note TN2065 for
	// more on shell scripting via AppleScript.
	String unixCmd
	String igorCmd
	svar mainpath
	string folder=parsefilepath(5,mainpath,"/",0,0)+"/Optimization"
	unixCmd = "cd '"+folder+"'; sh 'doOptim.sh' > /dev/null 2>&1"
	sprintf igorCmd, "do shell script \"%s\"", unixCmd
	ExecuteScriptText/UNQ igorCmd
	Print S_value		// For debugging only.
	return S_value
End
function runOptimization()
	wave test0
	nvar isWindows
	svar OptimizationPath
	wave status0
	status0[0]=0
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath status0 as "status.txt"
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:results:networkmat as "priorMat.txt"
	Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:results:passmat as "passMat.txt"
	if (isWindows)
		windowsShell()
	else
		UnixShellCommand()
	endif
	do
		try
			LoadWave/G/O/J/n=status/p=OptimizationPath/M/Q/Z "status.txt"
			AbortonRTE
		catch
			variable err=getRTError(1)
		endtry
		sleep/s 0.2
	while (status0[0]==0)
	LoadWave/G/O/J/n=test/p=OptimizationPath/M/Q "optimParams.txt"
end
function windowsShell()
	svar mainpath
	string folder=parsefilepath(5,mainpath,"\\",0,0)+"\\Optimization"
	ExecuteScriptText "cmd.exe /K \""+folder+"\\optim.bat "+folder+"\\\""
end
function runSim(x,y,z,conc,basepath,newFoldername)
	// Basepath folder must have 'errorFunc.m', 'findParams.m', 'optim.m', and 'setUp.m' inside of it
	variable x,y,z,conc
	string basepath,newFoldername
	variable/G isWindows=1
	string mainfolder=parsefilepath(5,basepath,"\\",0,0)
	ExecuteScriptText "cmd.exe /K \""+mainfolder+"\\setUp.bat "+mainfolder+" "+newFoldername+"\\\""
	string simfolder=basepath+newFoldername
	make_lines(x,y,z,conc,simfolder)
	thegel()
	newpath/o savepath simfolder
	saveexperiment/p=savepath/F={1,"",2} as newfoldername+".pxp"
end
function runSimMac(x,y,z,conc,basepath,newFoldername)
	// Basepath folder must have 'errorFunc.m', 'findParams.m', 'optim.m', and 'setUp.m' inside of it
	variable x,y,z,conc
	string basepath,newFoldername
	variable/G isWindows=0
	String unixCmd
	String igorCmd
	string mainfolder=parsefilepath(5,basepath,"/",0,0)
	unixCmd = "cd '"+mainfolder+"'; sh 'setUp.sh' "+mainfolder+" "+newfoldername
	sprintf igorCmd, "do shell script \"%s\"", unixCmd
	ExecuteScriptText/UNQ igorCmd
	string simfolder=basepath+newFoldername
	make_lines(x,y,z,conc,simfolder)
	thegel()
	newpath/o savepath simfolder
	saveexperiment/p=savepath/F={1,"",2} as newfoldername+".pxp"
end
function updateUnfoldHistory()
	wave unfoldmat_history=root:results:unfoldmat_history
	wave unfoldmat=root:results:unfoldmat
	variable i,j
	nvar modules,n_lines
	insertpoints/M=2 dimsize(unfoldmat_history,2),1,unfoldmat_history
	for(i=0;i<n_lines;i+=1)
		wave sizew=$("root:results:sizew_"+num2str(i))
		for(j=0;j<modules;j+=1)
			if (sizew[j][0]==0.2)
				unfoldmat_history[j][i]=0
			elseif (sizew[j][0]==0.01)
				unfoldmat_history[j][i][dimsize(unfoldmat_history,2)-1]=1
			endif
		endfor
	endfor
end
function dotwo()
	variable i,j,f,c
	runSim(4,4,2,0.0009,"C:Users:phy-popalab-g:Desktop:Kirill:Important:Tests:","oldSim")
end
#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.
function run_all_length_analysis()
	iter()
	GelLength_hists()
	//get_x_avg()
end
function test_lengths(name,n_lines)
	string name
	variable N_lines
	setdatafolder name
	variable i,k
	variable modules=8
	make/o/n=(N_lines*modules,dimsize($("LineHistory_0"),2)) z_Mat
	for(i=0;i<N_lines;i+=1)
		wave LineHistory=$("LineHistory_"+num2str(i))
		for(k=i*modules;k<(i+1)*modules;k+=1)
				z_Mat[k][]=linehistory[k-i*modules][2][q]
		endfor
		make/o/n=(dimsize(LineHistory,2)-1) $("X_Line_"+num2str(i))
		wave x_wave=$("X_Line_"+num2str(i))
		wave killwave=$("X_"+num2str(i))
		killwaves killwave
		x_wave[]=sqrt((LineHistory[0][0][p]-LineHistory[modules-1][0][p])^2+(LineHistory[0][1][p]-LineHistory[modules-1][1][p])^2+(LineHistory[0][2][p]-LineHistory[modules-1][2][p])^2)
	endfor
end
function iter()
	variable F,N,i,minF,maxF,minN,maxN
	minF=30;MaxF=60;minN=2;maxN=6
	for(F=MinF;F<(MaxF+1);F+=10)
		if (F==30)
			for(N=5;N<(maxN+1);N+=1)
				for(i=0;i<5;i+=1)
					test_lengths("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i),N*N*2)
				endfor
			endfor
		else
			for(N=minN;N<(MaxN+1);N+=1)
				for(i=0;i<5;i+=1)
					test_lengths("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i),N*N*2)
				endfor
			endfor
		endif
	endfor
end
function top_n()
	variable F,N,i,j,k
	variable modules=8
	variable n_avg=3,max_pnts,min_pnts
	string namew="Top_5_gel_length"
	for(F=50;F<60;F+=10)
		for(N=5;N<6;N+=1)
			for(i=0;i<1;i+=1)
				setdatafolder "root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
				make/o/n=(N*N*2*modules,50002) z_Mat
				n_avg=5
				for(j=0;j<(N*N*2);j+=1)
					wave lineHis=$("LineHistory_"+num2str(j))
					for(k=j*modules;k<(j+1)*modules;k+=1)
						z_Mat[k][]=lineHis[k-j*modules][2][q]
					endfor
				endfor
				make/o/n=(dimsize(z_mat,1)) $(namew)
				wave curr=$(namew)
				make/o/n=(N*N*2*modules) curr_z
				for(j=0;j<dimsize(z_mat,1);j+=1)
					curr_z[]=z_mat[p][j]
					sort/R curr_z,curr_z
					max_pnts=0;min_pnts=0
					for(k=0;k<n_avg;k+=1)
						max_pnts+=curr_z[k]/n_avg
						min_pnts+=curr_z[N*N*2*modules-1-k]/n_avg
					endfor
					curr[j]=max_pnts-min_pnts
				endfor
				killwaves/Z curr_z
				print "Done with F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
			endfor
		endfor
	endfor
	get_avgs(namew)
end
function GelLength_hists()
	variable F,N,i,j,k,minF,maxF,minN,maxN
	variable modules=8
	variable n_avg=3,max_pnts,min_pnts
	minF=30;MaxF=60;minN=2;maxN=6
	for(F=MinF;F<(MaxF+1);F+=10)
		if (F==30)
			for(N=5;N<(MaxN+1);N+=1)
				for(i=0;i<5;i+=1)
					setdatafolder "root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
					wave z_mat
					print "Starting F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
					get_hist_length(z_mat,n*N*2)
					print "Done with F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
				endfor
			endfor
		else
			for(N=MinN;N<(MaxN+1);N+=1)
				for(i=0;i<5;i+=1)
					setdatafolder "root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
					wave z_mat
					print "Starting F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
					get_hist_length(z_mat,n*N*2)
					print "Done with F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
				endfor
			endfor
		endif
	endfor
	get_avgs("hist_Gel_Length")
end
function get_hist_length(w,n_lines)
	wave w
	variable n_lines
	variable i,j,k,l,r,q,A,s,t,binsize
	variable zbox=2/(0.0009)^(1/3)
	make/o/n=(dimsize(w,1)) hist_Gel_Length=0
	wave len=hist_Gel_Length
	r=zbox
	q=2
	s=0
	t=N_lines/5
	Make/O/T/N=1 T_Constraints
	binsize=4//11*exp(-0.015*N_lines)//Value-Tested
	T_Constraints[0] = {"K3 > 0","K4 < 100","K4 > 1"}
	for(i=0;i<(dimsize(w,1));i+=1)
		make/o/n=(dimsize(w,0)) newz
		newz[]=w[p][i]
		wavestats/Q newz
		Make/O newz_hist
		Histogram/C/N/B={v_min-20,binsize,(v_max-v_min+40)/binsize} newz,newz_hist
		wavestats/q newz_hist
		T_Constraints[0] = {"K3 > "+num2str(r-zbox/10),"K3 < "+num2str(r+zbox/10),"K4 < 100","K4 > 1"}//,"K3 < "+num2str(v_maxloc-v_minloc)}
		Make/D/N=5/O W_coef
		W_coef[0] = {0,t,s,r,q}
		FuncFit/H="10000"/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
		//curvefit/x=1/Q gauss newz_hist /D
		wave w_coef,w_sigma
		len[i]=abs(w_coef[3])
		r=abs(w_coef[3])
		q=abs(w_coef[4])
		if (q>20)
			q=5
		endif
		t=abs(w_coef[1])
		s=w_coef[2]
		//print num2str(i/(dimsize(w,1)-1)*100)+" % Done"
	endfor
	killwaves/Z T_constraints, w_coef, W_sqrtN, newz,newz_hist,fit_newz_hist
end
function Extended_Gel_Length(w,n_lines)
	wave w
	variable n_lines
	variable i,j,k,l,r,q,A,s,t,binsize,maxz,minz
	variable zbox=2/(0.0009)^(1/3)
	variable modules=8,interval=1
	make/o/n=(dimsize(w,1)-1) hist_Extended_Gel_Length=0
	wave len=hist_Extended_Gel_Length
	r=zbox
	q=2
	s=0
	t=N_lines/5
	for(i=0;i<(dimsize(w,1)-1);i+=1)
		make/o/n=0 newz
		for(j=0;j<(sqrt(N_lines/2));j+=1)
			wave LineHis=$("LineHistory_"+num2str(j))
			maxz=max(LineHis[0][2][i],LineHis[modules-1][2][i]);minz=min(LineHis[0][2][i],LineHis[modules-1][2][i])
			Make/O/N=2 A_wave; a_wave={minz,maxz}
			variable N_pnts=round((maxz-minz)/interval)
			if (N_pnts<2)
				N_pnts=2
			endif
			Interpolate2/T=1/N=(N_pnts)/Y=z_pnts a_wave
			duplicate/o newz,temp
			wave temp
			concatenate/NP/O {temp,z_pnts}, newz
		endfor
		Make/O/T/N=1 T_Constraints
		T_Constraints[0] = {"K3 > 0","K4 < 100","K4 > 1"}
		wavestats/Q newz
		binsize=3.49*v_sdev/(numpnts(newz)^(1/3))
		Make/O newz_hist
		Histogram/C/N/B={v_min-21,binsize,(v_max-v_min+38)/binsize} newz,newz_hist
		wavestats/q newz_hist
		T_Constraints[0] = {"K3 > "+num2str(r-zbox/10),"K3 < "+num2str(r+zbox/10),"K4 < 100","K4 > 1"}//,"K3 < "+num2str(v_maxloc-v_minloc)}
		Make/D/N=5/O W_coef
		W_coef[0] = {0,t,s,r,q}
		FuncFit/H="10000"/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
		//curvefit/x=1/Q gauss newz_hist /D
		wave w_coef,w_sigma
		len[i]=abs(w_coef[3])
		r=abs(w_coef[3])
		q=abs(w_coef[4])
		if (q>20)
			q=5
		endif
		t=abs(w_coef[1])
		s=w_coef[2]
		//print num2str(i/(dimsize(w,1)-1)*100)+" % Done"
	endfor
	killwaves/Z T_constraints, w_coef, W_sqrtN, newz,newz_hist,fit_newz_hist,temp,z_pnts,A_wave
end
function get_x_avg()
	variable F,N,i,j,k,minF,maxF,minN,maxN
	minF=30;MaxF=60;minN=2;maxN=6
	for(F=minF;F<(MaxF+1);F+=10)
		if (F==30)
			for(N=5;N<(MaxN+1);N+=1)
				for(i=0;i<5;i+=1)
					wave x_0=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_0")
					duplicate/o x_0,$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
					wave xavg=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
					for(j=1;j<(N*N*2);j+=1)
						wave xwave=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":X_"+num2str(j))
						xavg+=xwave
					endfor
					xavg/=(n*N*2)
				endfor
			endfor
		else
			for(N=minN;N<(MaxN+1);N+=1)
				for(i=0;i<5;i+=1)
					wave x_0=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_0")
					duplicate/o x_0,$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
					wave xavg=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
					for(j=1;j<(N*N*2);j+=1)
						wave xwave=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":X_"+num2str(j))
						xavg+=xwave
					endfor
					xavg/=(n*N*2)
				endfor
			endfor
		endif
	endfor
end
function get_avgs(namew)
	string namew
	variable i,F,N,j,minF,maxF,minN,maxN
	minF=30;maxf=60;minN=2;maxN=6
	for(F=minF;F<(MaxF+1);F+=10)
		if (F==30)
			for(N=5;N<(maxN+1);N+=1)
				for(i=0;i<5;i+=1)
					wave inwave=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":"+namew)
					if (i==0)
						duplicate/o inwave, $("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
						wave avg=$("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
					else
						avg+=inwave
					endif
				endfor
				avg/=5
			endfor
		else
			for(N=minN;N<(MaxN+1);N+=1)
				for(i=0;i<5;i+=1)
					wave inwave=$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":"+namew)
					if (i==0)
						duplicate/o inwave, $("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
						wave avg=$("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
					else
						avg+=inwave
					endif
				endfor
				avg/=5
			endfor
		endif
	endfor
end
function normalize()
	variable start,last,i,maxN,minN,j
	minN=2;maxN=6
	for(i=minN;i<=maxN;i+=1)
		wave wavg=$("N_"+num2str(i)+"_avg")
		start=wavemin(wavg)
		wavg-=start
		last=wavg[inf]
		wavg/=last
		setscale/i x,0,5,"s",wavg
		for(j=0;j<5;j+=1)
			wave w=$("N_"+num2str(i)+"_"+num2str(j))
			start=wavemin(w)
			w-=start
			last=w[inf]
			w/=last
			setscale/i x,0,5,"s",w
		endfor
	endfor


end
function calc_resid()
	variable N,i
	make/o/n=5 resid_wave
	make/o/n=5 N_wave
	for(N=2;N<=6;N+=1)
		duplicate/o $("N_"+num2str(N)+"_avg") $("resid_"+num2str(N))
		wave avg=$("N_"+num2str(N)+"_avg")
		wave resid=$("resid_"+num2str(N))
		for(i=0;i<5;i+=1)
			wave curr=$("N_"+num2str(N)+"_"+num2str(i))
			resid[]+=abs(avg[p]-curr[p])/5
		endfor
		N_wave[N-2]=n*n*2;resid_wave[N-2]=sum(resid)/numpnts(resid)
	endfor
end

function quick_resid(num)
	variable num
	variable i
	wave avg=$("N_"+num2str(num)+"_avg")
	for(i=0;i<5;i+=1)
		wave orig=$("N_"+num2str(num)+"_"+num2str(i))
		duplicate/o orig, $("res_"+"N_"+num2str(num)+"_"+num2str(i))
		wave res=$("res_"+"N_"+num2str(num)+"_"+num2str(i))
		res=avg-res
		if (i==0)
			display res
		else
			appendtograph res
		endif
	endfor
end