• API
• FAQ
• Tools
• Archive
daily pastebin goal
84%
SHARE
TWEET

# Untitled

a guest May 16th, 2018 125 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #pragma TextEncoding = "UTF-8"
2. #pragma rtGlobals=1     // Use modern global access method.
3.
4. //reserved letters: e i j p q r s t x y z
5. // Last changed February 17, 2017 by IP
6.
7. function PMF(F,N,DL,L0,G0,M0)
8.     variable F,n,L0,DL,G0,M0
9.     variable Lmax,kT,p,L,i,j,Dx,R,b,w, xb
10.     NVAR p_val= root:SIM:p_val
11.     NVAR sigmaG = root:SIM:sigmaG
12.     make/O/N=(N+1) pl=p_val
13.     make/O/N=(N+2,2) MinUxy=NaN
14.     //Morse parameters
15.     //M0=19.27  // in pN nm
16.     R=2.8//6//4             // size of a folded domain - changed to 2.49*sqrt(6)=6
17.     b=40            // Morse width
18.     //Gaussian parameters
19.     //G0=32.8   // in pN nm
20.     w=sigmaG*sqrt(2)//.2 //this is sigma*sqrt(2), Gaussian width
21.     xb=0.44     //.3  this is delta X, distance to transition state
22.     kT=4.11     // thermal energy
23.     Dx=0.001        // resolution landscape in nm
24.     ///////////////get U setup as a template///////
25.     L=L0+N*DL
26.     make/O/N=(L/Dx) U
27.     setscale/P x,0,Dx,U
28.     //get shortest WLC; E-1
29.     L=L0
30.     U=(kT/pl[0])*((L/4)*(1/(1-x/L)-1)-(x/4)+x^2/(2*L))-F*x
31.     FindAPeak 0.1,2,1, U
32.     MinUxy[1][0]=V_peakX
33.     Duplicate/O/R=(0,MinUxy[1][0]) U Cut_U
34.     string Ui="U"+num2str(1)
35.     string Si="S"+num2str(1)
36.     string Ti="T"+num2str(1)
37.     Duplicate/O Cut_U \$Si
38.     Duplicate/O Cut_U \$Ti
39.     Duplicate/O U \$Ui
40.
41.     ////////////////do the rest of the WLC///////////////////
42.     for(i=2;i<N+2;i+=1)
43.         Ui="U"+num2str(i)
44.         L=L0+(i-1)*DL
45.         U=(kT/pl[i])*((L/4)*(1/(1-x/L)-1)-(x/4)+x^2/(2*L))-F*x
46.         FindAPeak 0.1,2,1, U
47.         MinUxy[i][0]=V_peakX
48.         Duplicate/O/R=(MinUxy[i-1][0],MinUxy[i][0]) U Cut_U
49.         Ui="U"+num2str(i)
50.         Si="S"+num2str(i)
51.         Ti="T"+num2str(i)
52.         Duplicate/O Cut_U \$Si
53.         Duplicate/O Cut_U \$Ti
54.         Duplicate/O U \$Ui
55.     endfor
56.     ////////////////Add the Enthalpy to each segment and shift///////////////////
57.     WAVE S1
58.     Deletepoints numpnts(S1)-1,1,S1
59.     for(i=2;i<N+2;i+=1)
60.         string S0="S"+num2str(i-1)
61.         Duplicate/O \$S0 H0
62.         Si="S"+num2str(i)
63.         Duplicate/O \$Si H
64.         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)
65.         variable shift=H[0]-H0(numpnts(H0))
66.         H=H-shift
67.         Deletepoints numpnts(H)-1,1,H
68.         Duplicate/O H \$Si
69.     endfor
70.     ////////////////create final segment and shift///////////////////
71.     Si="S"+num2str(N+2)
72.     Ui="U"+num2str(N+1)
73.     Duplicate/O/R=(MinUxy[N+1][0],MinUxy[N+1][0]+100) \$Ui Last
74.     shift=H(numpnts(H))-Last[0]
75.     Last=Last+shift
76.     Duplicate/O last \$Si
77.     ////////////////Concatenate for final wave///////////////////
78.     string full=""
79.     for(i=1;i<N+3;i+=1)
80.         full=full+"S"+num2str(i)+";"
81.     endfor
82.     concatenate/O full, PMF_complete
83.
84.     //fill PMF_complete with NAN's above 500 pN nm
85.     findlevel/P/Q PMF_complete, 500
86.     for(j=round(V_LevelX); j<numpnts(PMF_complete);j+=1)
87.         PMF_complete[j]=nan
88.     endfor
89.
90.     //save MinUy at the force used
91.     variable z
92.     for(i=0;i<N+2;i+=1)
93.         z=PMF_complete(MinUxy[i][0])
94.         MinUxy[i][1]=z
95.     endfor
96.     //remove all other Waves//
97.     for(i=1;i<N+2;i+=1)
98.         Si="S"+num2str(i)
99.         Ui="U"+num2str(i)
100.         Ti="T"+num2str(i)
101.         killwaves/Z \$Si, \$Ui,\$Ti
102.     endfor
103.     Si="S"+num2str(N+2)
104.     killwaves/Z \$Si
105.     killwaves/Z H,H0,Last,Cut_U,U,PL
106. End
107.
108.
109. Function E_curves(N,delta_F,nF,DL,L0,G0,M0)
110.     variable N,delta_F,nF,Dl,L0,G0,M0
111.     variable j,i,F
112.     variable nmax
113.     make/O PMF_complete//dummy wave
114.     make/O/N=(N+2,2) MinUxy=NaN //dummy wave
115.     make/O/N=(N+2,2,nF) root:MinXYF=NaN
116.     wave MinXYF = root:MinXYF
117.     make/O/N=(nF) root:F_table
118.     wave F_table = root:F_table
119.     F_Table=1+delta_F*(nF-p)//the minimum force is always 1 pN............
120.     //see how long is the longest PMF at the highest force!!
121.     PMF(F_table[0],N,DL,L0,G0,M0)
122.     nmax=numpnts(PMF_complete)
123.     make/O/N=(nmax,nF) root:PMFM=NaN
124.     wave PMFM = root:PMFM
125.     setscale/P x,0, deltax(PMF_complete), PMFM
126.     setscale/P y,F_table[0], -delta_F, PMFM//set the scale in descending mode
127.     //create #nF PMF's and their Minima//
128.     for(i=0;i<nF;i+=1)
129.         F=F_table[i]
130.         PMF(F,N,DL,L0,G0,M0)//returns both PMF_complete and MinUxy
131.         PMFM[][i]=PMF_complete[p]
132.         MinXYF[][][i]=MinUxy[p][q]
133.     endfor
134.     MinXYF[0][0][]=0
135.     Display_PMF(N,nF)
136.     killwaves/Z MinUxy, PMF_complete
137.     setscale/I z,(1+delta_F*nF),(1+delta_F), MinXYF
138.
139. End
140.
141.
142. Function Display_PMF(N,nF)
143.     variable N, nF
144.     variable i
145.     Wave MinxyF = root:MinxyF
146.     Wave PMFM = root:PMFM
147.     DoWindow/K PMF0
148.     //display E_curves first
149.     Display/W=(500,10,1200,600)/N=PMF0/K=1 MinxyF[0][1][] vs MinxyF[0][0][]
150.     for(i=1;i<N+2;i+=1)
151.         Appendtograph MinxyF[i][1][] vs MinxyF[i][0][]
152.     endfor
153.     ModifyGraph lsize=2, RGB=(65535,0,0)
154.     //now display the PMF's
155.     for(i=0;i<nF+1;i+=10)//show only every 10 curves
156.         Appendtograph/C=(1,12815,52428) PMFM[][i]
157.     endfor
158.     Global_min()
159.     wave G_min = root:G_min
160.     AppendToGraph G_min[][1] vs G_min[][0]
161.     ModifyGraph lsize(G_min)=3,rgb(G_min)=(2,39321,1)
162. End
163.
164. function Global_min()
165.     NVAR N_modules = root:sim:N_modules
166.     NVAR L0 = root:sim:L0
167.     NVAR Lc = root:sim:Lc
168.     NVAR U0 = root:sim:U0
169.     NVAR A0 = root:sim:A0
170.     NVAR deltaF = root:SIM:deltaF
171.     NVAR noF = root:SIM:noF
172.
173.     wave MinXYF = root:MinXYF
174.     wave F_table = root:F_table
175.     make/O/N=(dimsize(MinXYF,2)-1,4) root:G_min
176.     wave G_min = root:G_min
177.     G_min[][2]=F_table[p]
178.     variable i
179.     make/O/N=(dimsize(MinXYF,0)) PMF_min
180.     for(i=0;i<(dimsize(MinXYF,2));i+=1)
181.         PMF_min= MinXYF[p][1][i]
182.         wavestats/Q PMF_min
183.         G_min[i][1]=MinXYF[V_minloc][1][i]
184.         G_min[i][0]=MinXYF[V_minloc][0][i]
185.         G_min[i][4]=V_minloc-1 // No. unfolded domains
186.     endfor
187.     Make/O/N=(numpnts(F_table)) DummyN
188.     DummyN=G_min[p][3]
189.     SetScale/P x, F_table[0], F_table[1]-F_table[0], DummyN
190.     Findlevel/Q DummyN,N_modules/2
191.     Printf "Half force is at %1.2f pN", V_LevelX
192.     Print " for U0:", U0/4.11, " KT, A0:", A0/4.11, "KT"
193.     killwaves/Z PMF_min
194.
195. End
196.
197.
198. Function Display_F_X()
199.     wave MinXYF = root:MinXYF
200.     wave G_min = root:G_min
201.     wave F_table = root:F_table
202.     make/O/N=(dimsize(MinXYF,2),2) E_1
203.     E_1[][0]= MinXYF[1][0][p]
204.     E_1[][1]= F_table[p]
205.     DoWindow/K F_X
206.     //display E_1
207.     Display/W=(500,10,1200,600)/N=F_X/K=1 E_1[][1] vs E_1[][0]
208.     SetAxis left 0,25;DelayUpdate
209.     SetAxis bottom 0,120
210.     ModifyGraph lsize=2,rgb=(1,3,39321)
211.     AppendToGraph G_min[][2] vs G_min[][0]
212.     ModifyGraph lsize=2,rgb(G_min)=(65000,0,0)
213. End
214.
215.
216.
217. /////////////////////New Code from here for GUI
219.     "-"
220.     "Start Sim Panel", InitializeVars()
221.
222. End
223.
224. Function InitializeVars()
225.     NewDataFolder/O root:SIM
226.     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
227.     Variable/G root:SIM:N_modules = 8
228.     Variable/G root:SIM:Lc = 18.6//19
229.     Variable /G root:SIM:L0 = 0.1
230.     Variable/G root:SIM:U0_KT = 1.7//3
231.     NVAR U0_KT=root:SIM:U0_KT
232.     Variable/G root:SIM:U0 = U0_KT*4.11
233.     //  SetFormula root:SIM:U0, "root:SIM:U0_KT*4.11"
234.     Variable/G root:SIM:A0_KT = 14.3//13
235.     NVAR A0_KT=root:SIM:A0_KT
236.     Variable/G root:SIM:A0 = A0_KT*4.11
237.     //  SetFormula root:SIM:A0, "root:Sim:A0_KT*4.11"
238.     Variable/G root:SIM:p_val = 0.58// persistence length
239.     Variable/G root:SIM:deltaF = 1  // min force step
240.     Variable/G root:SIM:noF = 150   // no steps in force
241.     Variable/G root:SIM:sigmaG = 0.1    // which comes from (2.51-2.49)*sqrt(6)=0.05 nm
242.     Variable/G root:SIM:dt=1e-7 // sampling time
243.     Variable/G  root:SIM:D_coef=15000
244.
245.
246.     BuildSimPanel()
247.     UpdatePMFMatrix("",0,"","")
248.     CreateSimFCCommand()
249.
250. End
251.
252. Function BuildSimPanel()
253.     Variable PanelXSize =   1090
254.     Variable PanelYSize =   550
255.     Make/O root:SIM:SIM_PMF01=NAN
256.     PauseUpdate
257.     String FldrSave= GetDataFolder(1)
258.     //  Define Position and Dimensions
259.     GetWindow kwFrameInner, wsizeDC             //  Get size of the Igor Main window
260.     Variable XOffset    =   10
261.     Variable    Left        =   V_right-PanelXSize-XOffset  //  The panel will be on the right side
262.     Variable    Top     =   0
263.     Variable    right       =   V_right-XOffset
264.     Variable bottom =   PanelYSize+Top
265.     //  Create the panel
266.     DoWindow/K MainFrame
267.     NewPanel/N=MainFrame /W=(left,top,right,bottom) as "SIM Mainframe"
268.     SetWindow MainFrame,hook(testhook)=KeyboardPanelHook
269.     // The globals
270.     SetVariable NDisplaySet, pos={49,32},size={58,21},title="N",variable=root:SIM:N_modules,proc = UpdatePMFMatrix
271.     SetVariable NDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.0f"
272.
273.     SetVariable U0DisplaySet,pos={123,32},size={110,21},title="U0, kT", variable=root:SIM:U0_kT, proc = UpdatePMFMatrix
274.     SetVariable U0DisplaySet,limits={1,inf,1}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"
275.
276.     SetVariable A0DisplaySet,pos={232,32},size={115,21},title="A0, kT",variable=root:SIM:A0_kt, proc = UpdatePMFMatrix
277.     SetVariable A0DisplaySet,limits={1,inf,1}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"
278.
279.     SetVariable L0DisplaySet,pos={12,77},size={100,21},title="L0, nm",variable=root:SIM:L0, proc = UpdatePMFMatrix
280.     SetVariable L0DisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"
281.
282.     SetVariable LcDisplaySet,pos={121,77},size={103,21},title="Lc, nm",variable=root:SIM:Lc, proc = UpdatePMFMatrix
283.     SetVariable LcDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"
284.
285.     SetVariable pDisplaySet,pos={237,77},size={89,2},title="p, nm",variable=root:SIM:p_val, proc = UpdatePMFMatrix
286.     SetVariable pDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"
287.
288.     SetVariable DDisplaySet,pos={222,121},size={120,21},title="D,nm2/s",variable=root:SIM:D_coef
289.     SetVariable DDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.0f"
290.
291.     SetVariable deltaFDisplaySet,pos={20,121},size={92,21},title="dF, pN",variable=root:SIM:deltaF, proc = UpdatePMFMatrix
292.     SetVariable deltaFDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.2f"
293.
294.     SetVariable noFDisplaySet,pos={150,121},size={60,21}, title="nF",variable=root:SIM:noF, proc = UpdatePMFMatrix
295.     SetVariable noFDisplaySet,limits={0,inf,0}, font="Arial",fSize=16,frame=5,fStyle=0, format="%1.0f"
296.
297.
298.     //The protocol and PMF figures
299.     SetActiveSubwindow MainFrame
300.     Display /HIDE=0 /W=(350,30,750,280)/HOST=MainFrame /N=ProtocolWavePlot root:SIM:CommandF
301.     Label/W=MainFrame#ProtocolWavePlot left "Force Protocol (pN)"
302.     Label/W=MainFrame#ProtocolWavePlot bottom "Time (s)"
303.     ModifyGraph/W=MainFrame#ProtocolWavePlot tick=2,mirror=1, lsize=2, zero(left)=3, tickUnit=1
304.
305.     Display /HIDE=0 /W=(350,300,750,550)/HOST=MainFrame /N=PMFPlot root:SIM:SIM_PMF01
306.     Label/W=MainFrame#PMFPlot left "PMF(pN*nm)"
307.     Label/W=MainFrame#PMFPlot bottom "Extension (nm)"
308.     ModifyGraph/W=MainFrame#PMFPlot tick=2,mirror=1, lsize=2, zero(left)=3, tickUnit=1
309.
310.     //  FC > Tables
311.     Edit /HIDE=0/HOST=MainFrame/W=(800,30,998,530)/N=ProtocolForceTable root:SIM:Force_List,root:SIM:Duration_List, root:SIM:Ramp_List
312.     SetWindow MainFrame#ProtocolForceTable , userdata(tabnum)=  "1"
313.     SetWindow MainFrame#ProtocolForceTable , userdata(tabcontrol)=  "tabequipment"
314.     ModifyTable/W=MainFrame#ProtocolForceTable alignment=1, width=65
315.     ModifyTable/W=MainFrame#ProtocolForceTable width(Point)=35, showParts=199
316.     ModifyTable/W=MainFrame#ProtocolForceTable title(root:SIM:Ramp_List)="Ramp on?",Title(root:SIM:Force_List)="Force",Title(root:SIM:Duration_List)="Time"
317.     SetActiveSubwindow MainFrame
318.
319.     // FC > Button > Delete Rows
320.     Button RemoveFromForceList,pos={1012,30},size={65,36},disable=0,proc=FCSimProtocolDeleteLastRow,title="Delete\rLast Row"
321.
322.     SetDataFolder FldrSave
323.     ResumeUpdate
324. End
325.
326. Function FCSimProtocolDeleteLastRow(ba) : ButtonControl
327.     STRUCT WMButtonAction &ba
328.
329.     switch( ba.eventCode )
330.         case 2: // mouse up
331.             Wave w1 = root:SIM:Force_List
332.             Wave w2 = root:SIM:Duration_List
333.             Wave w3 = root:SIM:Ramp_List
334.             Variable n = Max(Max(NumPnts(w1),NumPnts(w2)),Max(NumPnts(w3),1))-1
335.             If(n>0)
336.                 DeletePoints n,1, w1,w2,w3
337.             EndIF
338.             break
339.     endswitch
340.
341.     return 0
342. End
343.
344. Function UpdatePMFMatrix(ctrlName,varNum,varStr,varName) : SetVariableControl
345.     String ctrlName
346.     Variable varNum
347.     String varStr, varName
348.
349.     NVAR N_modules = root:sim:N_modules
350.     NVAR L0 = root:sim:L0
351.     NVAR Lc = root:sim:Lc
352.     NVAR U0 = root:sim:U0
353.     NVAR A0 = root:sim:A0
354.     NVAR deltaF = root:SIM:deltaF
355.     NVAR noF = root:SIM:noF
356.     NVAR U0_KT=root:SIM:U0_KT
357.     NVAR A0_KT=root:SIM:A0_KT
358.     U0=U0_KT*4.11
359.     A0=A0_KT*4.11
360.     E_curves(N_modules,deltaF,noF,Lc,L0,A0,U0)
361.     CreateSimPMF()
362.
363. End
364.
365. static constant hook_keyboard=11
366. static constant hook_mouse=3
367. Function KeyboardPanelHook(s)
368.     STRUCT WMWinHookStruct &s
369.     NVAR X_Stage =root:Config:PiezoStage:X_Stage
370.     NVAR Y_Stage = root:Config:PiezoStage:Y_Stage
371.     NVAR  StepPiezoStage = root:Config:PiezoStage:StepPiezoStage
372.
373.     switch(s.eventCode)
374.         case hook_keyboard:
375.             CreateSimFCCommand()
376.             CreateSimPMF()
377.             break
378.         case hook_mouse:
379.             CreateSimFCCommand()
380.             CreateSimPMF()
381.             break
382.     endswitch
383.     return 0
384. End
385.
386. Function CreateSimFCCommand()
387.     WAVE    Force_List      =   root:SIM:Force_List
388.     WAVE    Duration_List       =   root:SIM:Duration_List
389.     WAVE    Ramp_List       =   root:SIM:Ramp_List
390.     NVAR deltaF = root:SIM:deltaF
391.     NVAR noF = root:SIM:noF
392.     NVAR    dt = root:SIM:dt
393.     Variable Duration    =  Sum(Duration_List)
394.     Variable PointsPerTrace=round(Duration/(dt*1e3))
395.     Make/O/N=(PointsPerTrace)   root:SIM:CommandF=NaN   //  In units of pN
396.     Wave CF =   root:SIM:CommandF
397.     //  Setscale/I x,0,Duration, "s", CF
398.     Setscale/P x, 0,dt*1e3, "s", CF
399.     Variable        Time1
400.     Variable        Time2
401.     Variable        Point1
402.     Variable        Point2
403.     Variable        i       =   0
404.     Variable        nForces =   NumPnts(Force_List)
405.     Variable        RampSlope
406.     Variable        F1, F2
407.     Do
408.         Time1   =   Time2
409.         Time2   +=  Duration_List[i]
410.         Point1  =   Point2
411.         Point2  =   x2Pnt(CF,Time2)
412.         F1      =   Force_List[i]
413.         if (F1>noF*deltaF)
414.             Print "Please increase the dF/nF to cover higher forces"
415.             Force_List[i]=noF*deltaF
416.             F1 = Force_List[i]
417.         endif
418.         if (F1<2)
419.             Print "Minimum force is 2 pN"
420.             Force_List[i]=2
421.             F1 = Force_List[i]
422.         endif
423.         if (Ramp_List[i]==1)
424.             //          F2                  =   Force_List[i+1]
425.             //          RampSlope           =   (F2-F1)/(Point2 - Point1)
426.             //          CF[Point1, Point2]  =   F1+ RampSlope*(p-Point1)
427.             CF[Point1, Point2]  =   F1 //no ramps yet
428.         else
429.             CF[Point1, Point2]  =   F1
430.         endif
431.         i   +=  1
432.     While(i<nForces)
433.
434.     WaveStats/Q CF
435.     DoWindow/F MainFrame
436.     if( V_Flag!=0 )
437.         SetAxis/W=MainFrame#ProtocolWavePlot left 0,1.1*V_max
438.     endif
439. End
440.
441. Function CreateSimPMF()
442.     Wave    Force_List      =   root:SIM:Force_List
443.     Wave F_table = root:F_table
444.     Wave PMFM = root:PMFM
445.     NVAR N_modules = root:SIM:N_modules
446.     Variable i
447.     // remove old PMF waves from panel & kill them
448.     DoWindow/F MainFrame
449.     if( V_Flag!=0 )
450.         String FldrSave= GetDataFolder(1)
451.         SetDataFolder root:SIM
452.         do
453.             String list_string = TraceNameList("MainFrame#PMFPlot", ";",1)  //This function makes a string with all the names of the waves in graph
454.             RemoveFromGraph /W=MainFrame#PMFPlot \$StringFromList(0, list_string, ";")
455.             if(!StringMatch (StringFromList(0,list_string,";"), "MinXYF*"))
456.                 KillWaves/Z \$StringFromList(0, list_string, ";") // THIS IS WHY THE MinxyF & MinxyF SHOULD NOT BE IN THIS FOLDER
457.             endif
458.         while(ItemsInList   (list_string)!=1)
459.         SetDataFolder FldrSave
460.     endif
461.     //make new PMFs on which the simulation will run
462.     for(i=0;i<numpnts(Force_List);i+=1)
463.         String nameWave = "root:Sim:SIM_PMF" + num2str(round(Force_List[i]))
464.         Make/O/N=(dimsize(PMFM,0)) \$nameWave
465.         WAVE workWave = \$nameWave
466.         workWave = PMFM[x][binarysearch(F_table,Force_List[i])]
467.         setscale/P x,0, deltax(PMFM), workWave
468.         Appendtograph /C=(1,12815,52428)/W=MainFrame#PMFPlot \$nameWave
469.         Label/W=MainFrame#PMFPlot left "PMF(pN*nm)"
470.         Label/W=MainFrame#PMFPlot bottom "Extension (nm)"
471.     endfor
472.     // Set axis before adding the E curves
473.     SetAxis/A/W=MainFrame#PMFPlot /N=1 left;
474.     SetAxis/A/W=MainFrame#PMFPlot /N=1 bottom
475.     DoUpdate
476.     GetAxis/Q/W=MainFrame#PMFPlot left
477.     SetAxis/W=MainFrame#PMFPlot left V_min, V_max
478.     for(i=0;i<=N_modules+1;i+=1)
479.         wave MinxyF = root:MinxyF
480.         Appendtograph/W=MainFrame#PMFPlot root:MinxyF[i][1][] vs root:MinxyF[i][0][]
481.     endfor
482. End
483.
484. Function BuildSIMDisplay()
485.     Wave    F_wave      =   root:SIM:F_wave
486.     Wave    x_wave      =   root:SIM:x_wave
487.     Wave    CommandF        =   root:SIM:CommandF
488.     NVAR N_modules = root:sim:N_modules
489.     NVAR L0 = root:sim:L0
490.     NVAR Lc = root:sim:Lc
491.
492.     DoWindow/K FcDualDisplay
493.     Display/K=1 /W=(0,40,550,450)/L=ForceAxis/N=FcDualDisplay root:SIM:CommandF, root:SIM:F_wave
494.     AppendToGraph/W=FcDualDisplay/L=left root:SIM:x_wave
495.     ModifyGraph rgb(F_wave)=(0,12800,52224)
496.     ModifyGraph/W=FcDualDisplay grid(left)=1,grid(ForceAxis)=1
497.     ModifyGraph/W=FcDualDisplay grid(bottom)=0,tick=2,mirror=1
498.     ModifyGraph/W=FcDualDisplay gridRGB(ForceAxis)=(34816,34816,34816)
499.     ModifyGraph/W=FcDualDisplay gridRGB(left)=(34816,34816,34816)
500.     ModifyGraph/W=FcDualDisplay nticks(left)=10
501.     ModifyGraph/W=FcDualDisplay lblPosMode(ForceAxis)=2,lblPosMode(left)=2
502.     ModifyGraph/W=FcDualDisplay lblPos(ForceAxis)=50,lblPos(left)=50
503.     ModifyGraph/W=FcDualDisplay axisEnab(ForceAxis)={0,0.35}
504.     ModifyGraph/W=FcDualDisplay axisEnab(left)={0.4,1}
505.     ModifyGraph/W=FcDualDisplay freePos(ForceAxis)=0
506.     ModifyGraph/W=FcDualDisplay lstyle=0,lsize(CommandF)=1.1 ,rgb(CommandF)=(0,0,0)
507.     ModifyGraph/W=FcDualDisplay live= 1
508.     Label/W=FcDualDisplay ForceAxis "Force (pN)"
509.     Label/W=FcDualDisplay bottom "Time (\\U)"
510.     Label/W=FcDualDisplay left "Length (nm)"
511.     WaveStats/Q CommandF
512.     SetAxis/W=FcDualDisplay ForceAxis, 0.8*V_min, 1.1*V_max
513.     SetAxis/W=FcDualDisplay left, 0, L0+(N_modules+1)*Lc
514.     DoUpdate/W=FcDualDisplay
515. End
516. Function Find_A_and_U()
517.     NVAR N_modules = root:sim:N_modules
518.     NVAR L0 = root:sim:L0
519.     NVAR Lc = root:sim:Lc
520.     NVAR U0 = root:sim:U0
521.     NVAR A0 = root:sim:A0
522.     NVAR deltaF = root:SIM:deltaF
523.     NVAR noF = root:SIM:noF
524.     NVAR U0_KT=root:SIM:U0_KT
525.     NVAR A0_KT=root:SIM:A0_KT
526.
527.     Variable i, sumAU
528.     sumAU=16
529.     Make/O/N=(sumAU-1,3) UAN_matrix
530.     for(i=0;i<round(sumAU)-1;i+=1)
531.         U0_KT=i*0.1+1
532.         A0_KT=sumAU-U0_KT
533.         U0=U0_KT/4.11
534.         A0=A0_KT/4.11
535.         UpdatePMFMatrix("",0,"","")
536.         DoUpdate
537.         //Global_min()
538.         WAVE DummyN
539.         Findlevel/Q DummyN,4
540.         UAN_matrix[i][0]= U0_KT
541.         UAN_matrix[i][1] = A0_KT
542.         UAN_matrix[i][2] = V_LevelX
543.     endfor
544.     Display UAN_matrix[2][] vs UAN_matrix[0][]
545. End
546. Function plot_normal_distrib(sigma)
547.     Variable sigma
548.     Make/o/N=1000 normD
549.     Setscale/I x,-5*sigma, 5*sigma, normD
550.     normD=1/sqrt(2*(sigma^2)*pi)*exp(-(x)^2/(2*sigma^2))
551. End
552. //#pragma rtGlobals=    3       // Use modern global access method and strict wave access.
553. Function make_lines(x,y,z,conc,savepath)
554.     Variable x,y,z,conc
555.     string savepath
556.     if (x==0 || y==0 || z==0)
557.         abort "x, y, or z cannot be one. That's not even gel."
558.     endif
559.     NewDataFolder/O root:RESULTS
560.     Variable/G N_lines=x*y*z
561.     Variable/G xby=x
562.     Variable/G yby=y
563.     Variable/G zby=z
564.     Variable/G Length = 4 //nm
565.     Variable/G modules=8
566.     Variable/G conc_molec=conc//molec/nm^3 (from 1000 molec per 100 nm^3)
567.     Variable/G box_size
568.     Variable/G NumberPoints=1000
569.     Variable/G displayw=1
570.     Variable/G showconnect=0
571.     Variable/G LinkSize=4//twice a radius domain (2*2)
573.     mainpath=savepath
575.     NewPath/C/O/Q UnfoldingPath savepath+":Unfolding:"
576.     NewPath/C/O/Q OptimizationPath savepath+":Optimization"
577.     Variable i,j,k,l,m, CMdif,theta, phi,a,b,c,timecount
578.     String nameW, colorW, InfoNote, exec, NameSW
579.     Variable SetX, SetY, deltaStep, step, SetZ
580.     //Make the waves that define everything
581.     make/o/n=(n_lines,n_lines) root:results:ConnectMat=0
582.     make/o/n=(1,n_lines) root:results:a_list=0
583.     make/o/n=(1,n_lines) root:results:b_list=0
584.     make/o/n=(1,n_lines) root:results:c_list=0
585.     make/o/n=(1,n_lines) root:results:theta_list=0
586.     make/o/n=(1,n_lines) root:results:phi_list=0
587.     make/o/n=(n_lines) root:results:deltawave=0
588.     make/o/n=(modules,3,n_lines,0) root:results:LineMat
589.     wave linemat=root:results:linemat
590.     wave a_list=root:results:a_list
591.     wave b_list=root:results:b_list
592.     wave c_list=root:results:c_list
593.     wave theta_list=root:results:theta_list
594.     wave phi_list=root:results:phi_list
595.     wave deltawave = root:results:deltawave
596.     make/o/n=(modules,n_lines) root:results:lengthmat=0
597.     make/o/n=(modules,n_lines) root:results:UnfoldMat=0
598.     variable xbox,ybox,zbox,deltaStepx,deltaStepy,deltaStepz
599.     variable/G connectPerClust=3
600.     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
601.     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
602.     variable/G Dein=10^(-6)*4.11/(6*pi*0.000911*10^(-6))//without rg
603.     //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
604.     N_lines=xby*yby*zby
605.     xbox=xby/(conc_molec)^(1/3)
606.     //xbox=(xby^2/yby/zby)*(N_lines/conc_molec)^(1/3)
607.     deltaStepx = xbox/(xby)
608.     //ybox=(yby^2/xby/zby)*(N_lines/conc_molec)^(1/3)
609.     ybox=yby/(conc_molec)^(1/3)
610.     deltaStepy = ybox/(yby)
611.     //zbox=(zby^2/xby/yby)*(N_lines/conc_molec)^(1/3)
612.     zbox=zby/(conc_molec)^(1/3)
613.     deltaStepz = zbox/(zby)
614.     ColorTab2Wave rainbow
615.     WAVE M_Colors
616.     if(displayW)
617.         DoWindow/K DisplayG
618.         exec = "NewGizmo /W = (10,0,1000,1000)/N=DisplayG; ModifyGizmo ShowInfo"
619.         Execute exec
620.     endif
621.     make/o/n=(modules,n_lines) Root:Results:UnfoldMat=0////rows define domain, columns define line
622.     i=0
623.     for(j=1;j<xby+1;j+=1)
624.         setX=j*deltaStepx//deltaStep
625.         for(k=1;k<yby+1;k+=1)
626.             setY=k*deltaStepy//deltaStep
627.             for(l=1;l<zby+1;l+=1)
628.                 setZ=l*deltaStepz//deltaStep
629.                 nameW="Line_"+num2str(i)
630.                 Make/O/N=(modules,3) root:RESULTS:\$nameW
631.                 WAVE LineW = root:RESULTS:\$nameW
632.                 LineW=NAN
633.                 InfoNote=""
634.                 a=SetX+enoise(deltaStep/2)
635.                 b=SetY+enoise(deltaStep/2)
636.                 c=SetZ+enoise(deltaStep/2)
637.                 theta=enoise(pi/2)+pi/2
638.                 phi=enoise(pi)+pi
639.                 LineW[][0][]=(p-modules/2)*Length*sin(theta)*cos(phi)+a
640.                 LineW[][1][]=(p-modules/2)*Length*sin(theta)*sin(phi)+b
641.                 LineW[][2][]=(p-modules/2)*Length*cos(theta)+c
642.                 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)+";"
643.                 Note LineW, InfoNote
644.                 if(displayW)
645.                     exec = "ModifyGizmo setDisplayList=" + num2str(2*i)+ ", object=path"+ num2str(i)
646.                     Execute exec
647.
648.                     exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:"+ nameofwave(LineW)
649.                     Execute exec
650.
651.                     exec = "ModifyGizmo setDisplayList=" + num2str(2*i+1)+ ", object=scatter"+ num2str(i)
652.                     Execute exec
653.
654.                     exec = "AppendToGizmo/N=DisplayG  DefaultScatter=root:RESULTS:"+ nameofwave(LineW)
655.                     Execute exec
656.
657.                     exec ="ModifyGizmo ModifyObject=scatter"+ num2str(i)+", property={ size,0.2}"
658.                     Execute exec
659.                 endif
660.                 i+=1
661.             endfor
662.         endfor
663.     endfor
664.     if(displayw)
665.         exec="ModifyGizmo scalingOption=0"
666.         execute exec
667.         exec="ModifyGizmo setOuterBox={"+num2str(-3*xbox)+","+num2str(3*xbox)+","+num2str(-3*ybox)+","+num2str(3*ybox)+","+num2str(-6*zbox)+","+num2str(6*zbox)+"}"
668.         Execute exec
669.         exec="modifygizmo home={115,-110,23};modifygizmo gohome"
670.         execute exec
671.     endif
672.     variable connects
673.     variable limits
674.     limits=1
675.     connections()
676.     variable savetime=0
677.     nvar totalnot
678.     timecount=0
679.     if (displayw)
681.     endif
682.     savetime+=1
683.     getaggregates()
684.     wave agg=root:results:aggwave
685.     do
686.         moveit(sqrt(2*Dtran),sqrt(2*Drot))///
687.         for (i=0;i<numpnts(agg);i+=1)
688.             if (agg[i]>1)
689.                 movecluster(i)
690.             endif
691.         endfor
692.         connections()
693.         getaggregates()
694.         timecount+=1
695.         if (mod(timecount,10)==0)
697.             sleep/s 0.1
698.             savetime+=1
699.         endif
700.     while (numpnts(agg)!=1)
701.     if(displayw)
702.         exec="ModifyGizmo scalingOption=0"
703.         execute exec
704.         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)+"}"
705.         Execute exec
706.         exec="modifygizmo home={115,-110,23};modifygizmo gohome"
707.         execute exec
708.     endif
709.     if(showconnect)
710.         exec = "ModifyGizmo setDisplayList=" + num2str(2*N_lines-1)+ ", object=path"+ num2str(N_lines)
711.         Execute exec
712.
713.         exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:Crosslinks"
714.         Execute exec
715.
716.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={pathColor,0,0,0,1}"///65535
717.         Execute exec
718.
719.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={lineWidthType,1}"
720.         Execute exec
721.
722.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={linewidth,5}"
723.         Execute exec
724.     endif
725.     connections()
727.     savetime+=1
728.     for(i=0;i<n_lines;i+=1)
729.         namew="Line_"+num2str(i)
730.         wave linew=root:results:\$namew
731.         infonote=note(linew)
732.         a=NumberByKey(" a", infonote,"=")
733.         a_list[0][i]=a
734.         b=NumberByKey(" b", infonote,"=")
735.         b_list[0][i]=b
736.         c=NumberByKey(" c", infonote,"=")
737.         c_list[0][i]=c
738.         theta=NumberByKey(" theta", InfoNote, "=")
739.         theta_list[0][i]=theta
740.         phi=NumberByKey(" phi", InfoNote, "=")
741.         phi_list[0][i]=phi
742.     endfor
743.     //doproteinshift()
744.     CreateSizeW()
745.     updatecolors()
746.     variable xcm,ycm,zcm,r
747.     for(r=0;r<n_lines;r+=1)
748.         wave linew=root:results:\$("Line_"+num2str(r))
749.         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]
750.         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]
751.         linew[][2]=(p-modules/2)*(Length)*cos(theta_list[dimsize(theta_list,0)-1][r])+c_list[dimsize(c_list,0)-1][r]
752.         for(i=0;i<modules;i+=1)
753.             xcm+=linew[i][0]/(n_lines*modules)
754.             ycm+=linew[i][1]/(n_lines*modules)
755.             zcm+=linew[i][2]/(n_lines*modules)
756.         endfor
757.     endfor
758.     for(r=0;r<n_lines;r+=1)
759.         wave linew=root:results:\$("Line_"+num2str(r))
760.         linew[][0]-=xcm;linew[][1]-=ycm;linew[][2]-=zcm
761.     endfor
762.     variable/G timetotal=timecount*10^-4
764.     print num2str(getaggregates())+" clusters formed"
765.     //print "Gel polymerized in "+num2str(timecount*10^-4)+" seconds"
766. End
767. function moveit(stept,stepr)
768.     variable stept,stepr
769.     variable i
770.     variable connects,a,b,c,theta,phi,movea,moveb,movec,flag,movetheta,movephi,xbox,ybox,zbox
771.     string namew,infonote,exec
772.     nvar n_lines,modules,length,box_size,xby,yby,zby,conc_molec,Displayw
773.     variable pmx,pmy,pmz
774.     wave a_list=root:results:a_list
775.     wave b_list=root:results:b_list
776.     wave c_list=root:results:c_list
777.     xbox=(xby^2/yby/zby)*(N_lines/conc_molec)^(1/3)
778.     ybox=(yby^2/xby/zby)*(N_lines/conc_molec)^(1/3)
779.     zbox=(zby^2/xby/yby)*(N_lines/conc_molec)^(1/3)
780.     variable volumeFactor=3 // so the volume of the box is 1.5 times its intial volume
781.     pmx=xbox*(volumeFactor^(1/3)-1)/2
782.     pmy=ybox*(volumeFactor^(1/3)-1)/2
783.     pmz=zbox*(volumeFactor^(1/3)-1)/2
784.     for(i=0;i<n_lines;i+=1)
785.         namew="root:results:Line_"+num2str(i)
786.         wave linew=\$namew
787.         infonote=note(linew)
788.         connects = NumberByKey(" connections", InfoNote, "=")
789.         if (connects==0)
790.             a=NumberByKey(" a", InfoNote, "=")
791.             b=NumberByKey(" b", InfoNote, "=")
792.             c=NumberByKey(" c", InfoNote, "=")
793.             theta=NumberByKey(" theta", InfoNote, "=")
794.             phi=NumberByKey(" phi", InfoNote, "=")
795.             do
796.                 flag=0
797.                 movea=gnoise(stept)
798.                 moveb=gnoise(stept)
799.                 movec=gnoise(stept)
800.                 if ((movea+a)<(-pmx) || (movea+a)>(pmx+xbox))
801.                     flag+=1
802.                 endif
803.                 if (flag==0)
804.                     if ((moveb+b)<(-pmy) || (moveb+b)>(pmy+ybox))
805.                         flag+=1
806.                     endif
807.                 endif
808.                 if (flag==0)
809.                     if ((movec+c)<(-pmz) || (movec+c)>(pmz+zbox))
810.                         flag+=1
811.                     endif
812.                 endif
813.             while (flag)
814.             movetheta=gnoise(stepr*pi/180)
815.             movephi=gnoise(stepr*pi/180)
816.             LineW[][0][]=(p-modules/2)*Length*sin(theta+movetheta)*cos(phi+movephi)+a+movea
817.             LineW[][1][]=(p-modules/2)*Length*sin(theta+movetheta)*sin(phi+movephi)+b+moveb
818.             LineW[][2][]=(p-modules/2)*Length*cos(theta+movetheta)+c+movec
819.             infonote=ReplaceNumberByKey(" a",infonote,a+movea,"=")
820.             infonote=ReplaceNumberByKey(" b",infonote,b+moveb,"=")
821.             infonote=ReplaceNumberByKey(" c",infonote,c+movec,"=")
822.             infonote=ReplaceNumberByKey(" theta",infonote,theta+movetheta,"=")
823.             infonote=ReplaceNumberByKey(" phi",infonote,phi+movephi,"=")
824.             note/K linew,infonote
825.         endif
826.     endfor
827. end
828. function movecluster(clusternum)
829.     variable clusternum
830.     variable i,j
831.     variable connects,a,b,c,theta,phi,movea,moveb,movec,flag,movetheta,movephi,xbox,ybox,zbox,points
832.     string namew,infonote,exec
833.     variable cmx,cmy,cmz,pmx,pmy,pmz
834.     nvar n_lines,modules,length,box_size,xby,yby,zby,conc_molec,Displayw,dein
835.     wave a_list=root:results:a_list
836.     wave b_list=root:results:b_list
837.     wave c_list=root:results:c_list
838.     xbox=(xby^2/yby/zby)*(N_lines/conc_molec)^(1/3)
839.     ybox=(yby^2/xby/zby)*(N_lines/conc_molec)^(1/3)
840.     zbox=(zby^2/xby/yby)*(N_lines/conc_molec)^(1/3)
841.     pmx=xbox/xby
842.     pmy=ybox/yby
843.     pmz=zbox/zby
844.     wave clust=root:results:\$("Cluster_"+num2str(clusternum))
845.     points=numpnts(clust)*modules
846.     make/o/n=(1,3) totalLines=0
847.     for(i=0;i<numpnts(clust);i+=1)
848.         wave linew=root:results:\$("Line_"+num2str(clust[i]))
849.         Matrixop/o totalLines=catrows(totallines,linew)
850.         infonote=note(linew)
851.         a=NumberByKey(" a", InfoNote, "=")
852.         b=NumberByKey(" b", InfoNote, "=")
853.         c=NumberByKey(" c", InfoNote, "=")
854.         cmx+=a/numpnts(clust)
855.         cmy+=b/numpnts(clust)
856.         cmz+=c/numpnts(clust)
857.     endfor
858.     DeletePoints 0,1,totallines
859.     matrixop/o rg2=sumsqr(subtractmean(totallines,0))
860.     variable rg=sqrt(rg2[0]/points)
861.     variable step=sqrt(2*dein/rg)
862.     do
863.         flag=0
864.         movea=gnoise(step)
865.         moveb=gnoise(step)
866.         movec=gnoise(step)
867.         if ((movea+cmx)<(-pmx) || (movea+cmx)>(pmx+xbox))
868.             flag+=1
869.         endif
870.         if (flag==0)
871.             if ((moveb+cmy)<(-pmy) || (moveb+cmy)>(pmy+ybox))
872.                 flag+=1
873.             endif
874.         endif
875.         if (flag==0)
876.             if ((movec+cmz)<(-pmz) || (movec+cmz)>(pmz+zbox))
877.                 flag+=1
878.             endif
879.         endif
880.     while (flag)
881.     make/o/n=(numpnts(clust)*modules,3) root:Results:cluster_wave=0
882.     wave clust_wave=root:results:cluster_wave
883.     for(i=0;i<numpnts(clust);i+=1)
884.         wave linew=root:results:\$("Line_"+num2str(clust[i]))
885.         infonote=note(linew)
886.         connects = NumberByKey(" connections", InfoNote, "=")
887.         a=NumberByKey(" a", InfoNote, "=")
888.         b=NumberByKey(" b", InfoNote, "=")
889.         c=NumberByKey(" c", InfoNote, "=")
890.         theta=NumberByKey(" theta", InfoNote, "=")
891.         phi=NumberByKey(" phi", InfoNote, "=")
892.         LineW[][0]=(p-modules/2)*Length*sin(theta)*cos(phi)+a+movea
893.         LineW[][1]=(p-modules/2)*Length*sin(theta)*sin(phi)+b+moveb
894.         LineW[][2]=(p-modules/2)*Length*cos(theta)+c+movec
895.         infonote=ReplaceNumberByKey(" a",infonote,a+movea,"=")
896.         infonote=ReplaceNumberByKey(" b",infonote,b+moveb,"=")
897.         infonote=ReplaceNumberByKey(" c",infonote,c+movec,"=")
898.         note/K linew,infonote
899.         for(j=i*modules;j<i*modules+modules;j+=1)
900.             clust_wave[j][]=linew[j-i*modules][q]
901.         endfor
902.     endfor
903.     //makerotationmat(gnoise(pi/2),gnoise(pi/2))//change this to be physical
904.     //wave rotationmat
905.     //matrixop/o clust_wave=(clust_wave-rowrepeat(sumcols(clust_wave)/points,points)) x rotationmat^t + rowrepeat(sumcols(clust_wave)/points,points)
906.     //make/o/n=(modules,3) tempwave
907.     //for(i=0;i<numpnts(clust);i+=1)
908.     //  wave linew=root:results:\$("Line_"+num2str(clust[i]))
909.     //  for(j=i*modules;j<i*modules+modules;j+=1)
910.     //      linew[j-i*modules][]=clust_wave[j][q]
911.     //  endfor
912.     //  infonote=note(linew)
913.     //  theta=NumberByKey(" theta", InfoNote, "=")
914.     //  phi=NumberByKey(" phi", InfoNote, "=")
915.     //  tempwave[][0]=(p-modules/2)*Length*sin(theta)*cos(phi)
916.     //  tempwave[][1]=(p-modules/2)*Length*sin(theta)*sin(phi)
917.     //  tempwave[][2]=(p-modules/2)*Length*cos(theta)
918.     //  matrixop/o coords=sumcols(linew-tempwave)/modules
919.     //  infonote=ReplaceNumberByKey(" a",infonote,coords[0],"=")
920.     //  infonote=ReplaceNumberByKey(" b",infonote,coords[1],"=")
921.     //  infonote=ReplaceNumberByKey(" c",infonote,coords[2],"=")
922.     //  note/K linew,infonote
923.     //endfor
924. end
925. function makeRotationMat(b,c)
926.     variable b,c
927.     variable a=0
928.     wave RotationMat
929.     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)
930.     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)
931.     RotationMat[2][0]=-sin(b);RotationMat[2][1]=cos(b)*sin(c);RotationMat[2][2]=cos(b)*cos(c)
932. end
933. function connections()
935.     variable i,j,k,l,cmdif
936.     string namew,infonote,exec
937.     wave connectmat = root:results:ConnectMat
938.     wave deltawave=root:results:deltawave
939.     nvar displayw,showconnect
940.     make/O/n=0 root:results:Conditions
941.     wave conditions = root:results:Conditions
944.     Make/O/n=(0,5) root:RESULTS:Movement
945.     WAVE Movement=root:RESULTS:Movement
946.     for(i=0;i<N_lines;i+=1)
947.         for(j=i+1;j<N_lines;j+=1)
948.             WAVE L1=root:RESULTS:\$("Line_"+num2str(i))
949.             WAVE L2=root:RESULTS:\$("Line_"+num2str(j))
950.             for(k=0;k<modules;k+=1)
951.                 variable counter=0
952.                 if (counter>connectPerClust)
953.                     break
954.                 endif
955.                 for(l=0;l<modules;l+=1)
956.                     CMdif=sqrt((L1[k][0]-L2[l][0])^2+(L1[k][1]-L2[l][1])^2+(L1[k][2]-L2[l][2])^2)
957.                     if(CMdif<LinkSize)//4 is the size of twice the radius of a domain
967.                         InsertPoints/M=0 dimsize(Movement,0),1,Movement
968.                         connectmat[i][j]=1;connectmat[j][i]=1
969.                         Movement[dimsize(Movement,0)-1][0]=i    //molecule number
970.                         Movement[dimsize(Movement,0)-1][1]=k //domain number
971.                         Movement[dimsize(Movement,0)-1][2]=j //molecule number
972.                         Movement[dimsize(Movement,0)-1][3]=l //domain number
973.                         Movement[dimsize(Movement,0)-1][4]=Cmdif //Distance
974.                         counter+=1
975.                     endif
976.                 endfor
977.             endfor
978.         endfor
979.     endfor
980.     variable connect
981.     variable/G totalnot
982.     totalnot=0
983.     for(i=0;i<n_lines;i+=1)
984.         connect=0
985.         for(j=0;j<dimsize(movement,0);j+=1)
986.             if(movement[j][0]==i)
987.                 connect+=1
988.             endif
989.             if(movement[j][2]==i)
990.                 connect+=1
991.             endif
992.         endfor
993.         nameW="Line_"+num2str(i)
994.         WAVE LineW = root:RESULTS:\$nameW
995.         InfoNote = note(root:RESULTS:\$nameW)
996.         infonote=ReplaceNumberByKey(" connections", infonote, connect,"=")
997.         if (connect==0)
998.             totalnot+=1
999.         endif
1000.         Note/K LineW, InfoNote
1001.         variable interm
1002.         if (displayw)
1003.             ColorTab2Wave Rainbow
1004.             WAVE M_Colors
1005.             interm=(connect/modules)
1006.             if (interm>1)
1007.                 interm=0.99
1008.             endif
1009.             Variable pick=((dimsize(M_Colors,0)-1)*(interm))
1010.             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}"
1011.             Execute exec
1012.         endif
1013.     endfor
1014.     if(showconnect)
1015.         exec = "ModifyGizmo setDisplayList=" + num2str(2*N_lines-1)+ ", object=path"+ num2str(N_lines)
1016.         Execute exec
1017.
1018.         exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:Crosslinks"
1019.         Execute exec
1020.
1021.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={pathColor,0,0,0,1}"//////65535
1022.         Execute exec
1023.
1024.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={lineWidthType,1}"
1025.         Execute exec
1026.
1027.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={linewidth,5}"
1028.         Execute exec
1029.     endif
1030. end
1031. function fixtemplooklist()
1032.     wave w=root:templooklist
1033.     variable i
1034.     for(i=0;i<numpnts(w);i+=1)
1035.         if (w[i]!=w[i])
1036.             w[i]=-1
1037.         endif
1038.     endfor
1039. end
1040. function getaggregates()
1041.     wave connectmat=root:Results:connectmat
1042.     nvar n_lines
1043.     variable i,j,k,aggcount,r,fnan,initnonnan,finalnonnan
1044.     aggcount=0
1045.     make/o/n=(N_lines) looklist
1046.     make/o/n=0 root:results:aggwave
1047.     wave aggwave=root:results:aggwave
1048.     looklist[]=p
1049.     aggcount=0
1050.     fnan=0
1051.     do
1052.         duplicate/o looklist,templooklist
1053.         fixtemplooklist()
1054.         initnonnan=numbernonnan(looklist)
1055.         fnan=firstnonnan(looklist)
1056.         r=looklist[fnan]
1057.         looklist[fnan]=nan
1058.         make/o/n=0 dolist
1059.         for(i=0;i<n_lines;i+=1)
1060.             if (connectmat[r][i]==1)
1061.                 insertpoints numpnts(dolist),1,dolist
1062.                 dolist[numpnts(dolist)-1]=i
1063.                 looklist[i]=nan
1064.             endif
1065.         endfor
1066.         k=0
1067.         do
1068.             k=dodolist(dolist,looklist)
1069.         while (k==1)
1070.         finalnonnan=numbernonnan(looklist)
1071.         aggcount+=1
1072.         compare(templooklist,looklist,numpnts(aggwave))
1073.         insertpoints numpnts(aggwave),1,aggwave
1074.         aggwave[numpnts(aggwave)-1]=initnonnan-finalnonnan
1075.     while (isallnan(looklist)==0)
1076.     return aggcount
1077. end
1078. function compare(w1,w2,clusterno)
1079.     wave w1,w2
1080.     variable clusterno
1081.     variable len,i
1082.     make/o/n=0 \$("root:results:Cluster_"+num2str(clusterno))
1083.     wave clust=\$("root:results:Cluster_"+num2str(clusterno))
1084.     len=numpnts(w1)
1085.     for(i=0;i<len;i+=1)
1086.         if (w2[i]!=w2[i])
1087.             if (w1[i]!=-1)
1088.                 insertpoints numpnts(clust),1,clust
1089.                 clust[numpnts(clust)-1]=i
1090.             endif
1091.         endif
1092.     endfor
1093. end
1094. function numbernonnan(w)
1095.     wave w
1096.     variable i,j
1097.     for(i=0;i<(numpnts(w));i+=1)
1098.         if(numtype(w[i])==0)
1099.             j+=1
1100.         endif
1101.     endfor
1102.     return j
1103. end
1104. function firstnonnan(w)
1105.     wave w
1106.     variable i
1107.     for(i=0;i<numpnts(w);i+=1)
1108.         if(numtype(w[i])==0)
1109.             return i
1110.         endif
1111.     endfor
1112. end
1113. function isallnan(w)
1114.     wave w
1115.     variable i,r
1116.     for(i=0;i<numpnts(w);i+=1)
1117.         if (numtype(w[i])==0)
1118.             return 0
1119.         endif
1120.     endfor
1121.     return 1
1122. end
1123. function dodolist(dolist,looklist)
1124.     wave dolist,looklist
1125.     wave connectmat=root:results:connectmat
1126.     nvar n_lines
1127.     make/o/n=0 tempdolist
1128.     variable i,r,j,t,find
1129.     find=0
1130.     for(i=0;i<numpnts(dolist);i+=1)
1131.         r=dolist[i]
1132.         for(j=0;j<n_lines;j+=1)
1133.             if (connectmat[r][j]==1)
1134.                 if (isinwave(looklist,j))
1135.                     insertpoints numpnts(tempdolist),1,tempdolist
1136.                     tempdolist[numpnts(tempdolist)-1]=j
1137.                     looklist[j]=nan
1138.                     find=1
1139.                 endif
1140.             endif
1141.         endfor
1142.     endfor
1143.     if (numpnts(tempdolist)==0)
1144.         return Find
1145.     else
1146.         redimension/N=(numpnts(tempdolist)) dolist
1147.         dolist=tempdolist
1148.         return find
1149.     endif
1150. end
1151. function isinwave(w,pnt)
1152.     wave w
1153.     variable pnt
1154.     variable j
1155.     for(j=0;j<numpnts(w);j+=1)
1156.         if(w[j]==pnt)
1157.             return 1
1158.         endif
1159.     endfor
1160.     return 0
1161. end
1162. function moveandconnect(stept,stepr)
1163.     variable stept,stepr
1164.     moveit(stept,stepr)
1165.     connections()
1166. end
1167. Function CreateSizeW()
1168.     NVAR N_lines
1169.     Variable i
1170.     String NameSW
1171.     for(i=0;i<N_lines;i+=1)
1172.         NameSW = "SizeW_"+num2str(i)
1173.         Make/D/O/N=(8,3) root:RESULTS:\$NameSW=0.2
1174.     endfor
1175. End
1176. function linelength(w)
1177.     // returns the line length its geometry
1178.     wave w
1179.     nvar modules
1180.     return distance(w,0,w,modules-1)
1181. end
1182. function BuildLineCopies()
1183.     // Builds a Copy of the current "Line_#" as "LineCopy_#"
1184.     // Line copies are the ones modulated during the optimzaiton
1185.     string namew,infonote,exec
1186.     variable a,b,c,theta,phi,setx,sety,setz,deltastep,i,LengthLine
1187.     nvar modules,length,n_lines,numberpoints
1188.     wave a_list=root:results:a_list
1189.     wave b_list=root:results:b_list
1190.     wave c_list=root:results:c_list
1191.     wave theta_list=root:results:theta_list
1192.     wave phi_list=root:results:phi_list
1193.     wave movement=root:results:movement
1194.     wave lengthmat=root:results:lengthmat
1195.     wave unfoldmat=root:results:unfoldmat
1196.     make/o/n=(dimsize(movement,0)) root:results:\$("CriticalDis")
1197.     wave criticaldis=root:results:\$("CriticalDis")
1198.     duplicate/o lengthmat,root:results:\$("CompareMat")
1199.     //criticaldis[]=movement[p][4]
1200.     InsertPoints/M=0 dimsize(a_list,0), 1, a_list
1201.     InsertPoints/M=0 dimsize(b_list,0), 1, b_list
1202.     InsertPoints/M=0 dimsize(c_list,0), 1, c_list
1203.     InsertPoints/M=0 dimsize(theta_list,0), 1, theta_list
1204.     InsertPoints/M=0 dimsize(phi_list,0), 1, phi_list
1205.     for(i=0;i<n_lines;i+=1)
1206.         a_list[dimsize(a_list,0)-1][i]=a_list[dimsize(a_list,0)-2][i]
1207.         b_list[dimsize(b_list,0)-1][i]=b_list[dimsize(b_list,0)-2][i]
1208.         c_list[dimsize(c_list,0)-1][i]=c_list[dimsize(c_list,0)-2][i]
1209.         theta_list[dimsize(theta_list,0)-1][i]=theta_list[dimsize(theta_list,0)-2][i]
1210.         phi_list[dimsize(phi_list,0)-1][i]=phi_list[dimsize(phi_list,0)-2][i]
1211.         namew="Line_"+num2str(i)
1212.         wave linew=root:results:\$namew
1213.         infonote=note(linew)
1214.         namew="LineCopy_"+num2str(i)
1215.         make/o/n=(modules,3) root:results:\$namew
1216.         wave linetot=root:results:\$namew
1217.         LengthLine=linelength(root:results:\$("Line_"+num2str(i)))
1218.         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])
1219.         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])
1220.         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])
1221.     endfor
1222.     for(i=0;i<dimsize(movement,0);i+=1)
1223.         wave w1=root:results:\$("line_"+num2str(movement[i][0]))
1224.         wave w2=root:results:\$("line_"+num2str(movement[i][2]))
1225.         movement[i][4]=distance(w1,movement[i][1],w2,movement[i][3])
1226.         wave w1=root:results:\$("linecopy_"+num2str(movement[i][0]))
1227.         wave w2=root:results:\$("linecopy_"+num2str(movement[i][2]))
1228.         criticaldis[i]=distance(w1,movement[i][1],w2,movement[i][3])
1229.     endfor
1230. end
1231. function distance(w1,a,w2,b)
1232.     // returns the ditance between "a" indexed tuple on wave w1 vs the "b" indexed tuple on wave w2
1233.     wave w1,w2
1234.     variable a,b
1235.     return Sqrt(((w1[a][0]-w2[b][0])^2+(w1[a][1]-w2[b][1])^2+(w1[a][2]-w2[b][2])^2))
1236. end
1237. Function DistanceOptimization(w,passwave)
1238.     // functiont that returns an associated error with a given set of dynamical paramers
1239.     // Only modifies angular and positional parameters, does not touch length
1240.     // i.e considers constant lengthwave with variable a,b,c,theta,phi
1241.     wave w
1242.     wave passwave
1244.     wave criticaldis=root:results:CriticalDis
1245.     wave movement=root:results:movement
1246.     make/o/n=(dimsize(movement,0)) root:results:CompareWave
1247.     wave Comparewave=root:results:CompareWave
1248.     wave deltawave=root:results:deltawave
1249.     CompareWave[]=Movement[p][4]
1250.     wave a_list=root:results:a_list
1251.     wave b_list=root:results:b_list
1252.     wave c_list=root:results:c_list
1253.     wave theta_list=root:results:theta_list
1254.     wave phi_list=root:results:phi_list
1255.     wave lengthmat=root:results:lengthmat
1256.     variable i,j
1257.     make/o/n=(n_lines) awavecurr
1258.     make/o/n=(n_lines) bwavecurr
1259.     make/o/n=(n_lines) cwavecurr
1260.     make/o/n=(n_lines) thetawavecurr
1261.     make/o/n=(n_lines) phiwavecurr
1262.     make/o/n=(modules,n_lines) lengthmatcurr
1263.     for(i=0;i<N_lines;i+=1)
1264.         awavecurr[i]=passwave[i]
1265.         bwavecurr[i]=passwave[i+n_lines]
1266.         cwavecurr[i]=passwave[i+n_lines*2]
1267.         thetawavecurr[i]=passwave[i+n_lines*3]
1268.         phiwavecurr[i]=passwave[i+n_lines*4]
1269.     endfor
1270.     for(i=0;i<modules;i+=1)
1271.         for(j=0;j<n_lines;j+=1)
1272.             lengthmatcurr[i][j]=passwave[n_lines*(5+i)+j]
1273.         endfor
1274.     endfor
1275.     updateLinesandDis(awavecurr,bwavecurr,cwavecurr,thetawavecurr,phiwavecurr,lengthmatcurr)
1276.     variable energy=0
1277.     variable LengthWeight=1,RepulseWeight=1
1278.     duplicate/o criticaldis,tempwave
1279.     //Term that fixes the crosslinks between timesteps
1281.     energy=sum(tempwave)
1282.     //for(i=0;i<numpnts(criticaldis);i+=1)
1284.     //endfor
1285.     //Term that weights the lengthchanges
1286.     //for(i=0;i<modules;i+=1)
1287.     //  for(j=0;j<n_lines;j+=1)
1288.     //      energy+=LengthWeight*(lengthmat[i][j]-lengthmatcurr[i][j])^2
1289.     //  endfor
1290.     //endfor
1291.     //Term that fixes the repulsion between timesteps
1292.     //for(i=0;i<numpnts(criticaldis);i+=1)
1293.     //  energy+=1/(criticaldis[i])^2
1294.     //endfor
1295.     return energy
1296. end
1297. function buildpasswave()
1298.     // generates the One-D wave from the 5-D set of dynamic parameters a,b,c,theta,phi
1299.     nvar n_lines,modules
1300.     wave a_list=root:results:a_list
1301.     wave b_list=root:results:b_list
1302.     wave c_list=root:results:c_list
1303.     wave theta_list=root:results:theta_list
1304.     wave phi_list=root:results:phi_list
1305.     wave lengthmat=root:results:lengthmat
1306.     wave criticaldis=root:results:CriticalDis
1307.     wave movement=root:results:movement
1308.     make/o/n=(n_lines*(5+modules)) root:results:Passwave
1309.     wave passwave=root:results:passwave
1310.     variable i,j
1311.     for(i=0;i<n_lines;i+=1)
1312.         passwave[i]=a_list[dimsize(a_list,0)-1][i]
1313.         passwave[i+n_lines]=b_list[dimsize(b_list,0)-1][i]
1314.         passwave[i+n_lines*2]=c_list[dimsize(c_list,0)-1][i]
1315.         passwave[i+n_lines*3]=theta_list[dimsize(theta_list,0)-1][i]+enoise(Pi/8)
1316.         passwave[i+n_lines*4]=phi_list[dimsize(phi_list,0)-1][i]+enoise(Pi/4)
1317.         //passwave[i+n_lines*3]=enoise(Pi/2)+Pi/2
1318.         //passwave[i+n_lines*4]=enoise(Pi)+Pi
1319.     endfor
1320.     for(i=0;i<modules;i+=1)
1321.         for(j=0;j<n_lines;j+=1)
1322.             passwave[n_lines*(5+i)+j]=lengthmat[i][j]
1323.         endfor
1324.     endfor
1325. end
1326. function DoProteinShift()
1327.     // Finds the set of a,b,c,theta,phi parameters that minimize the error
1328.     nvar n_lines,modules
1329.     variable i,j
1330.     string namew
1331.     //updateLengthMat()
1332.     BuildLineCopies()
1333.     buildpasswave()
1334.     wave passwave=root:results:passwave
1335.     DistanceOptimization(passwave,passwave)
1336.     duplicate/O passwave,root:results:x_param_wave
1337.     wave x_Param_wave=root:results:x_param_wave
1338.     Optimize/A=0/X=X_param_wave/Y=1/I=100/M={0,0}/R=X_param_wave/Q DistanceOptimization,passwave
1339.     print "error = "+num2str(V_min)
1340.     //x_param_wave gives set of a,b,c,theta,phi that minimize the error
1341.     DistanceOptimization(passwave,x_param_wave)
1342.     wave a_list=root:results:a_list
1343.     wave b_list=root:results:b_list
1344.     wave c_list=root:results:c_list
1345.     wave theta_list=root:results:theta_list
1346.     wave phi_list=root:results:phi_list
1347.     wave lengthmat=root:results:lengthmat
1348.     wave deltawave=root:results:deltawave
1349.     for(i=0;i<N_lines;i+=1)
1350.         a_list[dimsize(a_list,0)-1][i]=x_param_wave[i]
1351.         b_list[dimsize(b_list,0)-1][i]=x_param_wave[i+n_lines]
1352.         c_list[dimsize(c_list,0)-1][i]=x_param_wave[i+n_lines*2]
1353.         theta_list[dimsize(theta_list,0)-1][i]=x_param_wave[i+n_lines*3]
1354.         phi_list[dimsize(phi_list,0)-1][i]=x_param_wave[i+n_lines*4]
1355.     endfor
1356.     for(i=0;i<modules;i+=1)
1357.         for(j=0;j<n_lines;j+=1)
1358.             lengthmat[i][j]=x_param_wave[n_lines*(5+i)+j]
1359.         endfor
1360.     endfor
1361.     updatelines()
1362. end
1363. function updateLengthMat()
1364.     // Describes the change in length for each domain on the protein
1365.     //Note that the lengthmat is actually the change between two timepoints
1366.     wave lengthmat=root:results:lengthmat
1367.     wave display_lengthmat=root:results:display_lengthmat
1368.     wave unfoldmat=root:results:unfoldmat
1369.     variable shift,j,k,domain_no,i,nchange
1370.     nvar n_lines,modules,displayw
1371.     wave deltaLwave=root:results:deltawave
1372.     wave networkmat=root:results:networkmat
1373.     string exec
1374.     lengthmat=display_lengthmat
1375.     display_lengthmat=0
1376.     //matricies are labeled nameMat[domain_no][module_no]([j][i])
1377.     for(i=0;i<n_lines;i+=1)//i is domain number
1378.         wave sizew=root:results:\$("Sizew_"+num2str(i))
1379.         wave nwave=root:results:\$("N_"+num2str(i))
1380.         if ((nwave[numpnts(nwave)-1]>nwave[numpnts(nwave)-2]))
1381.             nchange=nwave[numpnts(nwave)-1]-nwave[numpnts(nwave)-2]
1382.             for(k=0;k<nchange;k+=1)
1383.                 Domain_no=getRandomPoint(i,0)
1384.                 unfoldmat[domain_no][i]=deltaLwave[i]/nchange
1385.             endfor
1386.         elseif ((nwave[numpnts(nwave)-1]<nwave[numpnts(nwave)-2]))
1387.             nchange=nwave[numpnts(nwave)-2]-nwave[numpnts(nwave)-1]
1388.             for(k=0;k<nchange;k+=1)
1389.                 Domain_no=getRandomPoint(i,1)
1390.                 unfoldmat[domain_no][i]=0
1391.             endfor
1392.         endif
1393.         for(j=0;j<modules;j+=1)
1394.             shift=unfoldmat[j][i]
1395.             for(k=0;k<modules;k+=1)
1396.                 if (k<j)
1397.                     display_lengthmat[k][i]+=-shift/2*(1-krondelta(j,k))
1398.                 else
1399.                     display_lengthmat[k][i]+=shift/2*(1-krondelta(j,k))
1400.                 endif
1401.             endfor
1402.             sizew[j][]=0.01*(1-krondelta(0,shift))+0.2*krondelta(0,shift)
1403.         endfor
1404.         if (displayw)
1405.             exec ="ModifyGizmo ModifyObject=scatter"+num2str(i)+", property={sizetype,1}"
1406.             Execute exec
1407.             exec="ModifyGizmo ModifyObject=scatter"+num2str(i)+", property={sizewave,root:RESULTS:" + nameofwave(SizeW) + "}"
1408.             Execute exec
1409.         endif
1410.     endfor
1411.     matrixop/O lengthmat=display_lengthmat-lengthmat
1412. end
1413. function UpdateLines()
1414.     // transfers over line copies to actual line wave
1415.     variable i,j
1416.     nvar n_lines,showconnect,modules
1417.     string exec
1418.     wave movement=root:results:movement
1419.     wave lengthmat=root:results:lengthmat
1420.     wave comparemat=root:results:comparemat
1423.     for(i=0;i<n_lines;i+=1)
1424.         wave w1=root:results:\$("Linecopy_"+num2str(i))
1425.         wave w2=root:Results:\$("Line_"+num2str(i))
1426.         w2=w1
1427.     endfor
1428.     for(i=0;i<dimsize(movement,0);i+=1)
1429.         wave w1=root:results:\$("line_"+num2str(movement[i][0]))
1430.         wave w2=root:results:\$("line_"+num2str(movement[i][2]))
1440.         movement[i][4]=distance(w1,movement[i][1],w2,movement[i][3])
1441.     endfor
1442.     if(showconnect)
1443.         exec = "ModifyGizmo setDisplayList=" + num2str(2*N_lines-1)+ ", object=path"+ num2str(N_lines)
1444.         Execute exec
1445.
1446.         exec = "AppendToGizmo/N=DisplayG Path=root:RESULTS:Crosslinks"
1447.         Execute exec
1448.
1449.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={pathColor,0,0,0,1}"//////65535
1450.         Execute exec
1451.
1452.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={lineWidthType,1}"
1453.         Execute exec
1454.
1455.         exec ="ModifyGizmo ModifyObject=path"+ num2str(N_lines)+", property={linewidth,5}"
1456.         Execute exec
1457.     endif
1458. end
1459. function getrandompoint(line_no,bool)//bool is zero for folded->unfolded and 1 otherwise
1460.     // function that returns a random number for a domain to fold/unfold
1461.     variable line_no,bool
1462.     variable j,rand
1463.     nvar modules
1464.     wave unfoldmat=root:Results:unfoldmat
1465.     make/o/n=0 dummywave
1466.     if (bool==0)
1467.         for(j=0;j<modules;j+=1)
1468.             if (unfoldmat[j][line_no]==0)
1469.                 insertpoints numpnts(dummywave),1,dummywave
1470.                 dummywave[numpnts(dummywave)-1]=j
1471.             endif
1472.         endfor
1473.     else
1474.         for(j=0;j<modules;j+=1)
1475.             if (unfoldmat[j][line_no]!=0)
1476.                 insertpoints numpnts(dummywave),1,dummywave
1477.                 dummywave[numpnts(dummywave)-1]=j
1478.             endif
1479.         endfor
1480.     endif
1481.     if (numpnts(dummywave)==0)
1482.         print "Wtf"
1483.     endif
1484.     rand= abs(round((enoise(0.5)+0.5)*(numpnts(dummywave)-1)))
1485.     return dummywave[rand]
1486. end
1487. function KronDelta(i,j)
1488.     // simply a kronecker delta
1489.     variable i,j
1490.     if(i==j)
1491.         return 1
1492.     else
1493.         return 0
1494.     endif
1495. end
1496. function updateLinesandDis(a_list,b_list,c_list,theta_list,phi_list,lengthmat)
1497.     // updates the line copies given the current LengthMat and set of dynamic parameters
1498.     wave ,,a_list,b_list,c_list,theta_list,phi_list,lengthmat
1499.     variable i
1500.     nvar n_lines,modules,length
1501.     string namew
1502.     wave criticaldis=root:results:\$("CriticalDis")
1503.     wave movement=root:results:movement
1504.     //wave lengthmat=root:results:lengthmat
1505.     for(i=0;i<n_lines;i+=1)
1506.         namew="LineCopy_"+num2str(i)
1507.         make/o/n=(modules,3) root:results:\$namew
1508.         wave linetot=root:results:\$namew
1509.         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])
1510.         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])
1511.         Linetot[][2][]=(p-modules/2)*(Length)*cos(theta_list[i])+c_list[i]+Lengthmat[p][i]*cos(theta_list[i])
1512.     endfor
1513.     for(i=0;i<dimsize(movement,0);i+=1)
1514.         wave w1=root:results:\$("linecopy_"+num2str(movement[i][0]))
1515.         wave w2=root:results:\$("linecopy_"+num2str(movement[i][2]))
1516.         criticaldis[i]=distance(w1,movement[i][1],w2,movement[i][3])
1517.     endfor
1518. end
1519. function Saveit(path,name)
1520.     // function to save a snap shot of the gizmo
1521.     // need to change the path destination if not working on Kirill's MacBook Pro
1522.     string name,path
1523.     nvar displayw
1525.     if (displayw)
1527.             SavePICT/O/E=-6/RES=2400/p=Crosslinkingpath as name+".png"
1528.         elseif (stringmatch(path,"Unfolding"))
1529.             SavePICT/O/E=-6/RES=2400/p=unfoldingpath as name+".png"
1530.         endif
1531.     endif
1532. end
1533. function updateColors()
1534.     // updates the colors of the lines given their angle relative to the z-axis
1535.     wave theta_list=root:results:theta_list
1536.     string namew,exec,infonote
1537.     variable i,j,force
1538.     nvar n_lines,modules,displayw
1539.     if (displayw)
1540.         for(i=0;i<N_lines;i+=1)
1541.             force=abs(cos(theta_list[dimsize(theta_list,0)-1][i]))
1542.             ColorTab2Wave Rainbow
1543.             WAVE M_Colors
1544.             Variable pick=((dimsize(M_Colors,0)-1)*(force))
1545.             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}"
1546.             Execute exec
1547.         endfor
1548.     endif
1549. end
1550. function ztothists()
1551.     wave zwave=root:results:zdistmat
1552.     variable i,j,k,l
1553.     nvar n_lines, modules
1554.     for(i=0;i<(dimsize(zwave,0));i+=100)
1555.         make/o/n=(n_lines*modules) \$("ZDist_"+num2str(i))
1556.         wave curr=\$("ZDist_"+num2str(i))
1557.         for(j=0;j<(n_lines*modules);j+=1)
1558.             curr[j]=zwave[i][j]
1559.         endfor
1560.     endfor
1561. end
1562. function GelLengthHist()
1563.     wave zwave=root:results:zdistmat
1564.     variable i,j,k,l,r,q,A,s,t,binsize
1565.     nvar n_lines, modules,zby,conc_molec
1566.     variable zbox=zby/(conc_molec)^(1/3)
1567.     make/o/n=(dimsize(zwave,0)-1) root:results:Gel_Length
1568.     make/o/n=(dimsize(zwave,0)-1) root:results:s_Gel_Length
1569.     make/o/n=(dimsize(zwave,0)-1) A_wave
1570.     wave len=root:results:Gel_Length
1571.     wave s_len=root:results:s_Gel_Length
1572.     r=zbox
1573.     q=2
1574.     s=0
1575.     t=N_lines/5
1576.     Make/O/T/N=1 T_Constraints
1577.     binsize=4//11*exp(-0.015*N_lines)//Value-Tested
1578.     T_Constraints[0] = {"K3 > 0","K4 < 100","K4 > 1"}
1579.     for(i=0;i<(dimsize(zwave,0)-1);i+=1)
1580.         make/o/n=(dimsize(zwave,1)) newz
1581.         newz[]=zwave[i][p]
1582.         Make/N=100/O newz_hist
1583.         wavestats/Q newz
1584.         l=v_min-20
1585.         Histogram/C/N/B={l,binsize,100} newz,newz_hist
1586.         wavestats/q newz_hist
1587.         T_Constraints[0] = {"K3 > "+num2str(r-zbox/10),"K3 < "+num2str(r+zbox/10),"K4 < 100","K4 > 1"}//,"K3 < "+num2str(v_maxloc-v_minloc)}
1588.         Make/D/N=5/O W_coef
1589.         W_coef[0] = {0,t,s,r,q}
1590.         FuncFit/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
1591.         //curvefit/x=1/Q gauss newz_hist /D
1592.         wave w_coef,w_sigma
1593.         len[i]=abs(w_coef[3])
1594.         s_len[i]=abs(w_sigma[3])
1595.         A_wave[i]=w_coef[4]
1596.         r=abs(w_coef[3])
1597.         q=abs(w_coef[4])
1598.         if (q>20)
1599.             q=5
1600.         endif
1601.         t=abs(w_coef[1])
1602.         s=w_coef[2]
1603.         print num2str((i/(dimsize(zwave,0))*100))+"% Done"
1604.     endfor
1605.     //drops(len)
1606.     avgsFixes()
1607. end
1608. function Calc_GelLength()
1609.     //Need to create Z-wave still
1610.     variable i,j,k,l,r,q,A,s,t,binsize,u
1611.     nvar n_lines, modules,zby,conc_molec,const_0,const_1,const_2,const_3
1612.     wave len = root:results:Gel_Length
1613.     wave s_len = root:results:s_Gel_Length
1614.     wave A_wave
1615.     //Inc is the the increment of the projection segmentation per line
1616.     variable Inc=0.01
1617.     variable zbox=zby/(conc_molec)^(1/3)
1618.     insertpoints numpnts(len),1,len,s_len,A_wave
1619.     make/o/n=0 newz
1620.     for(i=0;i<n_lines;i+=1)
1621.         wave linew=\$("Root:results:Line_"+num2str(i))
1622.         for(j=min(linew[0][2],linew[modules-1][2]);j<max(linew[0][2],linew[modules-1][2]);j+=inc)
1623.             insertpoints numpnts(newz),1,newz
1624.             newz[numpnts(newz)-1]=j
1625.         endfor
1626.     endfor
1627.     Make/O/T/N=1 T_Constraints
1628.     Make/N=100/O newz_hist
1629.     wavestats/Q newz
1630.     l=v_min-40
1631.     //binsize=4//11*exp(-0.015*N_lines)//Value-Tested
1632.     binsize=3.49*V_sdev*(numpnts(newz))^(-1/3)
1633.     u=(v_max+40-l)/binsize
1634.     Histogram/C/N/B={l,binsize,u} newz,newz_hist
1635.     wavestats/Q newz_hist
1636.     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)}
1637.     Make/D/N=5/O W_coef
1638.     W_coef[0] = {0,const_3,const_2,const_0,const_1}
1639.     FuncFit/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
1640.     //curvefit/x=1/Q gauss newz_hist /D
1641.     wave w_coef,w_sigma
1642.     len[numpnts(len)-1]=abs(w_coef[3])
1643.     s_len[numpnts(len)-1]=abs(w_sigma[3])
1644.     A_wave[numpnts(len)-1]=w_coef[4]
1645.     const_0=abs(w_coef[3])
1646.     const_1=abs(w_coef[4])
1647.     if (const_1>20)
1648.         const_1=5
1649.     endif
1650.     const_3=abs(w_coef[1])
1651.     const_2=w_coef[2]
1652. end
1653. function Calc_GelLength_ionel()
1654.     //Need to create Z-wave still
1655.     variable i,j,k,l,r,q,A,s,t,binsize,u
1656.     nvar n_lines, modules,zby,conc_molec,xby,yby
1657.     wave len = root:results:Gel_Length
1658.     variable avgnum=xby*yby
1659.     insertpoints numpnts(len),1,len
1660.     //Inc is the the increment of the projection segmentation per line
1661.     make/o/n=(modules*n_lines) zwave
1662.     k=0
1663.     for(i=0;i<N_lines;i+=1)
1664.         wave linew=\$("root:results:line_"+num2str(i))
1665.         for(j=0;j<modules;j+=1)
1666.             zwave[k]+=linew[j][2]
1667.             k+=1
1668.         endfor
1669.     endfor
1670.     sort zwave,zwave
1671.     variable maxavg=0,minavg=0
1672.     for(i=0;i<avgnum;i+=1)
1673.         maxavg+=zwave[numpnts(zwave)-i]
1674.         minavg+=zwave[i]
1675.     endfor
1676.     print (maxavg-minavg)/avgnum
1677.     len[numpnts(len)-1]=(maxavg-minavg)/avgnum
1678. end
1679. function avgsFixes()
1680.     nvar n_lines
1681.     variable i,j,k
1682.     wave duration = root:sim:duration_list
1683.     wave gel_length=root:results:Gel_Length
1684.     make/o/n=(numpnts(root:results:x_0)) root:results:x_avg
1685.     wave x_avg= root:results:x_avg
1686.     x_avg=0
1687.     for(i=0;i<n_lines;i+=1)
1688.         wave xwave=\$("root:results:x_"+num2str(i))
1689.         wave f=\$("root:results:F_"+num2str(i))
1690.         if (i==0)
1691.             duplicate/o F,Root:results:F_avg
1692.             wave favg=Root:results:F_avg
1693.         else
1694.             favg+=F
1695.         endif
1696.         setscale/I x,0,(sum(duration)),"s",xwave
1697.         x_avg[]+=xwave[p]
1698.     endfor
1699.     setscale/I x,0,(sum(duration)),"s",gel_length,x_avg
1700. end
1701. function drops(w)
1702.     wave w
1704.     nvar n_lines, modules,zby,conc_molec
1705.     variable zbox=zby/(conc_molec)^(1/3)
1706.     variable cut=zbox/5
1707.     variable pointnum=100
1708.     variable threshhold=0.5
1709.     variable numflag,startp,endp
1710.     i=0
1711.     flag=1
1712.     do
1713.         if (i==8760)
1714.             print ""
1715.         endif
1716.         startp=i
1717.         r=0
1718.         do
1719.             i+=1
1720.             endp=i
1721.             r+=1
1722.             if (r==200)
1723.                 break
1724.             endif
1725.         while ((round(w[startp])-round(w[endp]))>threshhold)
1727.         low=1e10
1728.         for(j=startp;j<endp;j+=1)
1729.             if(w[j]>w[startp])
1731.             endif
1732.             if(w[j]<low)
1733.                 low=w[j]
1734.             endif
1735.         endfor
1737.             avg=(w[endp]+w[startp])/2
1738.             for(i=startp;i<endp;i+=1)
1739.                 w[i]=avg
1740.             endfor
1741.         endif
1742.         if(i==(numpnts(w)-2))
1743.             flag=0
1744.         endif
1745.     while (flag)
1746. end
1747. function zdists()
1748.     wave zdistwave=root:results:zdistmat
1749.     variable i,j,k,r,l
1750.     nvar n_lines,modules
1751.     k=0
1752.     for(i=0;i<n_lines;i+=1)
1753.         wave linew=\$("Root:results:Line_"+num2str(i))
1754.         for(j=0;j<modules;j+=1)
1755.             zdistwave[dimsize(zdistwave,0)-1][k]=linew[j][2]
1756.             k+=1
1757.         endfor
1758.     endfor
1759.     insertpoints dimsize(zdistwave,0),1,zdistwave
1760. end
1761. Function SuperGauss(w,x) : FitFunc
1762.     WAVE w;Variable x
1763.     //CurveFitDialog/
1764.     //CurveFitDialog/ Coefficients 5
1765.     //CurveFitDialog/ w[0] = y0
1766.     //CurveFitDialog/ w[1] = A
1767.     //CurveFitDialog/ w[2] = x0
1768.     //CurveFitDialog/ w[3] = Width
1769.     //CurveFitDialog/ w[4] = P
1770.     return w[0]+w[1]*exp(-(((x-w[2])/w[3])^2)^(w[4]))
1771. End
1772. //Function that corrects the crosslink distances
1773. function AppendLinesAndNetwork()
1774.     // Appends an elemt to Linematrix and builds the the passmat, network mat representing the past step
1775.     // Line copies are the ones modulated during the optimzaiton
1776.     updateLengthMat()
1777.     string namew,exec
1778.     variable i,j,k,m,modnum,linenum
1779.     nvar modules,length,n_lines,numberpoints
1780.     wave LineMatrix=root:results:LineMatrix
1781.     wave passmat=root:results:passmat
1782.     wave networkmat=root:results:networkmat
1783.     wave domainconnections=root:Results:domainconnections
1784.     wave pivot=root:results:pivot
1785.     wave a_list=root:results:a_list
1786.     wave b_list=root:results:b_list
1787.     wave c_list=root:results:c_list
1788.     wave theta_list=root:Results:theta_list
1789.     wave phi_list=root:results:phi_list
1790.     wave lengthmat=root:results:lengthmat
1791.     //shift passmat points according to lengthmat
1792.     //passmat is past iteration of lengthmat plus the shift according to angles set on past interation
1793.     // Shifts are due to slight extensions and unfoldings
1794.     k=0
1795.     redimension/N=(modules*n_lines) domainconnections
1796.     for(i=0;i<numpnts(domainconnections);i+=1)
1797.         if (domainconnections[i]==1)
1798.             modnum=mod(i,modules);linenum=(i-mod(i,modules))/modules
1799.             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])
1800.             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])
1801.             passmat[k][2]=networkmat[k][2]+Lengthmat[modnum][linenum]*cos(theta_list[dimsize(theta_list,0)-1][linenum])
1802.             k+=1
1803.         endif
1804.     endfor
1805.     redimension/N=(modules,n_lines) domainconnections
1806. end
1807. function EnergyCalc(w,PassMat)
1808.     //Calculates the strech energy between the new set of coordinates
1809.     //Inputs: PriorMat is the Matrix of crosslinks from the PAST iteration
1810.     // PassMat is the matrix of connections being energy calculated
1811.     //Potential imporvements: Combine the inter- and intra-line connections into
1812.     // a single for-loop
1813.     wave w,PassMat
1814.     wave priorMat=root:results:networkmat
1815.     NVAR N_Lines, modules,linksize,length
1816.     NVAR p_val = root:sim:P_val
1817.     wave domainconnections=root:results:domainconnections
1818.     wave clindex=root:results:clindex
1819.     wave bendindex=root:results:bendindex
1820.     matrixop/o deltamat=PassMat-PriorMat
1821.     variable i,j,energy=0,k
1822.
1823.     // Energy PreFactors
1824.     variable StrechPrefactor=2*p_val*4.11/pi/(length/2)^4   //from article and wiki (does not include length divider (1/prior))
1825.     variable BendPrefactor=p_val*4.11/2 // from article and wiki (does not include length divider (1/prior^3))
1826.     variable CrossLinkPrefactor=372030//calculated from OPLS force-field
1827.
1828.     //Determine Difference in Streching between crosslinks
1829.     k=0
1830.     for(i=0;i<n_lines;i+=1)
1831.         //CLsum is # of crosslink locations on that line
1832.         matrixop/o CLSum=sum(col(domainconnections,i))
1833.         //Compute Strech Energy between Line Crosslinks
1834.         for(j=k;j<k+CLSum[0]-1;j+=1)
1835.             Matrixop/o DeltaU=row(deltamat,j)-row(deltamat,j+1)
1836.             matrixop/o OldVecLength=powr(sum(powr(row(priormat,j)-row(priormat,j+1),2)),0.5)
1837.             Matrixop/o UnitOld=normalize(row(PriorMat,j)-row(priormat,j+1))
1840.         endfor
1841.         //Compute Bending Between Crosslinks
1842.         if(CLSUm[0]>2)
1843.             for(j=k;j<k+CLSum[0]-2;j+=1)
1844.                 matrixop/O DeltaU=row(deltaMat,j)+row(deltamat,j+1)-2*row(deltamat,j+2)
1845.                 matrixop/o OldVecLength=powr(sum(powr(row(priormat,j)-row(priormat,j+1),2)),0.5)
1846.                 Matrixop/o UnitOld=normalize(row(PriorMat,j)-row(priormat,j+1))
1847.                 redimension/N=3 deltaU,unitold
1848.                 cross DeltaU,unitOld
1849.                 wave w_cross
1850.                 energy+=(BendPrefactor/oldveclength[0]^3)*(w_cross[0]^2+w_cross[1]^2+w_cross[2]^2)
1851.             endfor
1852.         endif
1853.         k+=CLSum[0]
1854.     endfor
1855.     //Calculate crosslink difference in strech and linksize
1856.     for(i=0;i<dimsize(clindex,0);i+=1)
1857.         //Calculate Strech energy between crosslinks
1858.         Matrixop/o DeltaU=row(deltamat,clindex[i][0])-row(deltamat,clindex[i][1])
1859.         matrixop/o OldVecLength=powr(sum(powr(row(priormat,clindex[i][0])-row(priormat,clindex[i][1]),2)),0.5)
1860.         Matrixop/o UnitOld=normalize(row(PriorMat,clindex[i][0])-row(priormat,clindex[i][1]))
1863.         //Claculate from link distances
1864.         Matrixop/o LengthVector=row(passmat,clindex[i][0])-row(passmat,clindex[i][1]))
1867.     endfor
1868.     for(i=0;i<dimsize(bendindex,0);i+=1)
1869.         //Calculate Bending between crosslinks
1870.         matrixop/O DeltaU=row(deltaMat,bendindex[i][2])+row(deltamat,bendindex[i][0])-2*row(deltamat,bendindex[i][1])
1871.         matrixop/o OldVecLength=powr(sum(powr(row(priormat,bendindex[i][0])-row(priormat,bendindex[i][1]),2)),0.5)
1872.         Matrixop/o UnitOld=normalize(row(PriorMat,bendindex[i][0])-row(priormat,bendindex[i][1]))
1873.         redimension/N=3 deltaU,unitold
1874.         cross DeltaU,unitOld
1875.         wave w_cross
1876.         energy+=(BendPrefactor/oldveclength[0]^3)*(w_cross[0]^2+w_cross[1]^2+w_cross[2]^2)
1877.     endfor
1878.     return energy
1879. end
1880. function UpdateVariables()
1881.     //This function updates all line parameter waves and updates the networkmat
1882.     //STILL NEED TO DO THE CHANGE IN LENGTH
1883.     nvar modules,length,n_lines,numberpoints
1884.     wave a_list=root:results:a_list
1885.     wave b_list=root:results:b_list
1886.     wave c_list=root:results:c_list
1887.     wave theta_list=root:results:theta_list
1888.     wave phi_list=root:results:phi_list
1889.     wave passmat=root:results:passmat
1890.     wave networkMat=root:results:networkMat
1891.     wave domainconnections=root:results:domainconnections
1892.     wave linematrix=root:results:linematrix
1893.     wave lengthmat=root:Results:lengthmat
1894.     wave pivot=root:results:pivot
1895.     variable i,j,k,deltheta,delphi,dela,delb,delc
1896.     //pivot passmat to zero
1897.     matrixop/o ChangeMat=PassMat-NetWorkMat
1898.     matrixop/o avgcol=averagecols(ChangeMat)
1899.     matrixop/o domaincols=sumcols(domainconnections)
1900.     k=0
1901.     //update all the variabnle lists
1902.     InsertPoints/M=0 dimsize(a_list,0), 1, a_list, b_list,c_list,theta_list,phi_list
1903.     for(i=0;i<n_lines;i+=1)
1904.         deltheta=0;delphi=0;dela=0;delb=0;delc=0;
1905.         for(j=k;j<k+domaincols[0][i];j+=1)
1906.             // 'del{x}' accounts for center of mass changes to that protien
1907.             dela+=Changemat[j][0];delb+=changemat[j][1];delc+=changemat[j][2]
1908.             if ((domaincols[0][i]>1) && (j!=(k+domaincols[0][i]-1)))
1909.                 // A single connection protein should not rotate
1914.             endif
1915.         endfor
1916.         a_list[dimsize(a_list,0)-1][i]+=a_list[dimsize(a_list,0)-2][i]+dela
1917.         b_list[dimsize(b_list,0)-1][i]+=b_list[dimsize(a_list,0)-2][i]+delb
1918.         c_list[dimsize(c_list,0)-1][i]+=c_list[dimsize(a_list,0)-2][i]+delc
1919.         theta_list[dimsize(theta_list,0)-1][i]=theta_list[dimsize(theta_list,0)-2][i]+deltheta
1920.         phi_list[dimsize(phi_list,0)-1][i]=phi_list[dimsize(phi_list,0)-2][i]+delphi
1921.         k+=domaincols[0][i]
1922.     endfor
1923.     //Current Iteration becomes optimized iteration
1924.     networkmat=passmat
1925. end
1926. function updateGraph()
1927.     nvar N_lines,modules,length
1928.     variable r,xcm=0,ycm=0,zcm=0,i
1929.     wave linematrix=root:results:linematrix
1930.     wave lengthmat=root:Results:display_lengthmat
1931.     wave a_list=root:results:a_list
1932.     wave b_list=root:results:b_list
1933.     wave c_list=root:results:c_list
1934.     wave theta_list=root:results:theta_list
1935.     wave phi_list=root:results:phi_list
1936.     //use variables from past iteration to construct lines for updates
1937.     for(r=0;r<n_lines;r+=1)
1938.         wave linew=root:results:\$("Line_"+num2str(r))
1939.         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])
1940.         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])
1941.         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])
1942.         for(i=0;i<modules;i+=1)
1943.             xcm+=linew[i][0]/(n_lines*modules)
1944.             ycm+=linew[i][1]/(n_lines*modules)
1945.             zcm+=linew[i][2]/(n_lines*modules)
1946.         endfor
1947.     endfor
1948.     for(r=0;r<n_lines;r+=1)
1949.         wave linew=root:results:\$("Line_"+num2str(r))
1950.         wave linehis=root:results:\$("LineHistory_"+num2str(r))
1951.         linew[][0]-=xcm;linew[][1]-=ycm;linew[][2]-=zcm
1952.         linehis[][][dimsize(linehis,2)-1]=linew[p][q]
1953.         insertpoints/M=2 dimsize(linehis,2),1,linehis
1954.     endfor
1955. end
1956. Function TheGel()
1957.     // The function that runs the gel
1958.     ///CHECK IF CLINDEX FUCKING WORKS
1959.     NVAR dt = root:SIM:dt
1960.     NVAR D  =   root:SIM:D_coef
1961.     WAVE    Force_List      =   root:SIM:Force_List
1962.     WAVE    Duration_List       =   root:SIM:Duration_List
1963.     WAVE    Ramp_List       =   root:SIM:Ramp_List
1964.     WAVE    PMFM = root:PMFM
1965.     Wave    CommandF = root:SIM:CommandF //workd with dt = 1e-7 and jmax = 1e3
1966.     Wave    MinXYF = root:MinXYF
1967.     Wave    F_table = root:F_table
1968.     WAVE deltawave = root:results:deltawave
1969.     WAVE theta_list = root:results:theta_list
1970.     wave movement=root:results:movement
1971.     wave lengthmat=root:Results:lengthmat
1972.     duplicate/o lengthmat,root:results:display_lengthmat
1973.     lengthmat=0
1974.     NVAR N_lines = root:N_lines
1975.     svar UnfoldingPath,OptimizationPath,mainpath
1976.     nvar displayw
1977.     nvar modules,zby,conc_molec
1978.     make/o/n=(modules,n_lines) root:results:unfoldmat_history=0
1979.     Variable  kT=4.11
1980.     Variable Duration=dt*1e3*numpnts(CommandF)
1981.     Variable i, j,k,l ,xc=3, tc=0,Fr=0, F_b=0, current_E_var, currentF, oldF
1982.     variable updateflag, numupdate
1983.     NVAR    N_modules = root:SIM:N_modules
1984.     Duplicate/O CommandF root:Sim:x_wave, root:Sim:F_wave,  root:Sim:N_wave
1985.     Wave x_wave = root:Sim:x_wave
1986.     Wave F_wave = root:Sim:F_wave
1987.     Wave N_wave = root:Sim:N_wave
1988.     x_wave=NaN; F_wave=NaN; N_wave=NaN
1989.     Variable D0=15000,f_factor,xcold
1990.     String nameWx, nameWf, nameWn,namehis
1991.     // some gel length variables
1992.     variable/G Const_0=zby/(conc_molec)^(1/3)
1993.     variable/G Const_1=2
1994.     variable/G Const_2=0
1995.     variable/G Const_3=N_lines/5
1996.     //convert all lines into into a matrix
1997.     //Note, Linematrix is not the same points that are used in display purposes.
1998.     make/o/n=(modules,3,n_lines) root:results:linematrix=0
1999.     make/o/n=3 root:results:pivot
2000.     wave linematrix=root:results:linematrix
2001.     wave pivot=root:results:pivot
2002.     for(i=0;i<n_lines;i+=1)
2003.         wave linew=root:results:\$("Line_"+num2str(i))
2004.         lineMatrix[][][i]=linew[p][q]
2005.     endfor
2006.     pivot[]=linematrix[0][p][0]
2007.     //Convert movement matrix into a connection matrix for operations in optimization
2008.     //Matrix value is 1 is connection exists at that domain, otherwise is zero
2009.     //Define the Domainconnections mat
2010.     make/o/n=(modules,n_lines) root:results:DomainConnections=0 //rows label modules, columns label lines
2011.     wave domainconnections=root:results:domainconnections
2012.     make/o/n=(0,2) root:results:CLindex_iter
2013.     wave clindex_iter = root:results:CLindex_iter
2014.     for(i=0;i<dimsize(movement,0);i+=1)
2015.         if (DomainConnections[movement[i][1]][movement[i][0]]==0)
2016.             DomainConnections[movement[i][1]][movement[i][0]]=1
2017.         endif
2018.         if (DomainConnections[movement[i][3]][movement[i][2]]==0)
2019.             DomainConnections[movement[i][3]][movement[i][2]]=1
2020.         endif
2021.         insertpoints/M=0 dimsize(clindex_iter,0),1,clindex_iter
2022.         clindex_iter[dimsize(clindex_iter,0)-1][0]=movement[i][1]*(movement[i][0]+1)
2023.         clindex_iter[dimsize(clindex_iter,0)-1][1]=movement[i][3]*(movement[i][2]+1)
2024.     endfor
2025.     //Create Network mat, this is the one the specifies the dynamics of the system
2026.     //The function Make_NetworkMat does this for arbitrary time
2027.     //Networkmat is the previous iteration of the network (priormat)
2028.     make/o/n=(sum(domainconnections),3) root:results:NetworkMat,root:results:PassMat
2029.     wave networkMat=root:results:NetworkMat
2030.     wave passmat=root:results:passmat
2031.     k=0
2032.     for(j=0;j<n_lines;j+=1)
2033.         for(i=0;i<modules;i+=1)
2034.             if(domainconnections[i][j]==1)
2035.                 networkMat[k][]=linematrix[i][q][j]
2036.                 k+=1
2037.             endif
2038.         endfor
2039.     endfor
2040.     //Make the crosslink_iter mat, values correspond to where connections are
2041.     //unpack domainconnections to allow for indexing
2042.     redimension/N=(modules*n_lines) domainconnections
2043.     make/o/n=(0,2) root:Results:clindex
2044.     wave clindex=root:results:clindex
2045.     //create temp_movement wave for searching
2046.     duplicate/o movement,temp_movement
2047.     for(j=0;j<n_lines;j+=1)
2048.         for(i=0;i<modules;i+=1)
2049.             for(k=0;k<dimsize(temp_movement,0);k+=1)
2050.                 //domain must be in one of the two places
2051.                 if (((temp_movement[k][0]==j) && (temp_movement[k][1]==i)) || ((temp_movement[k][2]==j) && (temp_movement[k][3]==i)))
2052.                     insertpoints/M=0 dimsize(clindex,0),1,clindex
2053.                     // the 'l-value' is: (Line_#)*modules+(Module_#)
2054.                     l=modules*temp_movement[k][0]+temp_movement[k][1]
2055.                     //Sum over domain connections specifies index placement, subtract one for 0-indexing
2056.                     clindex[dimsize(clindex,0)-1][0]=sum(domainconnections,0,l)-1
2057.                     l=modules*temp_movement[k][2]+temp_movement[k][3]
2058.                     clindex[dimsize(clindex,0)-1][1]=sum(domainconnections,0,l)-1
2059.                     //delete used point in temp_movement to prevent double-counting
2060.                     deletepoints k,1,temp_movement
2061.                 endif
2062.             endfor
2063.         endfor
2064.     endfor
2065.     // Make Bending wave
2066.     make/o/n=(0,3) root:results:bendindex
2067.     wave bendindex=root:Results:bendindex
2068.     variable smallindex,largeindex
2069.     make/o/n=3 tempwave
2070.     for(i=0;i<dimsize(clindex,0);i+=1)
2071.         smallindex=min(clindex[i][0],clindex[i][1]);largeindex=max(clindex[i][0],clindex[i][1])
2072.         tempwave[0]=smallindex;tempwave[1]=largeindex
2073.         for(j=i+1;j<dimsize(clindex,0);j+=1)
2074.             if (smallindex==clindex[j][0])
2075.                 tempwave[2]=clindex[j][1];sort tempwave,tempwave
2076.                 insertpoints/M=0 dimsize(bendindex,0),1,bendindex
2077.                 bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
2078.             elseif (smallindex==clindex[j][1])
2079.                 tempwave[2]=clindex[j][0];sort tempwave,tempwave
2080.                 insertpoints/M=0 dimsize(bendindex,0),1,bendindex
2081.                 bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
2082.             elseif (largeindex==clindex[j][0])
2083.                 tempwave[2]=clindex[j][1];sort tempwave,tempwave
2084.                 insertpoints/M=0 dimsize(bendindex,0),1,bendindex
2085.                 bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
2086.             elseif (largeindex==clindex[j][1])
2087.                 tempwave[2]=clindex[j][0];sort tempwave,tempwave
2088.                 insertpoints/M=0 dimsize(bendindex,0),1,bendindex
2089.                 bendindex[dimsize(bendindex,0)-1][]=tempwave[q]
2090.             endif
2091.         endfor
2092.     endfor
2093.     //Make some of the gel length waves
2094.     make/o/n=0 root:results:Gel_Length
2095.     make/o/n=0 root:results:s_Gel_Length
2096.     make/o/n=0 A_wave
2097.     //kill unused wave and redimension domainconnections
2098.     killwaves/Z temp_movement,tempwave
2099.     redimension/N=(modules,n_lines) domainconnections
2101.     Differentiate PMFM /D=D_PMFM
2102.     Make/O/N=(dimSize(MINXYF,0)) OneD
2103.     string nameflist
2104.     for(l=0;l<N_lines;l+=1)
2105.         nameWx="x_"+num2str(l)
2106.         nameWf="F_"+num2str(l)
2107.         nameWn="N_"+num2str(l)
2108.         namehis="LineHistory_"+num2str(l)
2109.         nameflist="F_command_"+num2str(l)
2110.         Make/O/N=1 root:results:\$nameWx, root:results:\$nameWf, root:results:\$nameWn
2111.         make/o/n=((modules),3,1) root:results:\$namehis
2112.         wave linehis=\$("root:results:LineHistory_"+num2str(l))
2113.         wave linew=\$("root:results:line_"+num2str(l))
2114.         WAVE Wx=\$("root:results:x_"+num2str(l))
2115.         WAVE Wf=\$("root:results:F_"+num2str(l))
2116.         WAVE Wn=\$("root:results:N_"+num2str(l))
2117.         wx=0
2118.         wf=0
2119.         wn=0
2120.         linehis[][][0]=linew[p][q]
2121.         insertpoints/M=2 dimsize(linehis,2),1,linehis
2122.         Setscale/P x,0,dt, "s", Wx, Wf, Wn
2123.     endfor
2124.     Make/O/N=(N_lines) theta_curr
2125.     Make/O/N=(N_lines) theta_past
2126.     make/o/n=(1,n_lines*modules) root:results:zdistmat
2127.     //first column counts iteration in current state, second denotes last state change
2128.     // 1 means last change was an unfolding, 0 means last change was a refolding
2129.     make/o/n=(N_lines,2) state=0;state[][1]=1
2130.     zdists()
2131.     make/o/n=1 status0=0
2132.     NVAR p_val = root:sim:P_val,length
2133.     make/o/n=3 prefactors
2134.     prefactors[0]=2*p_val*4.11/pi/(length/2)^4
2135.     prefactors[1]=5*p_val*4.11/2
2136.     prefactors[2]=372030
2137.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:Results:clIndex as "clindex.txt"
2138.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:Results:bendindex as "bendindex.txt"
2139.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:Results:domainconnections as "domcon.txt"
2140.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath status0 as "status.txt"
2141.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath prefactors as "prefactors.txt"
2142.     Print "Working on Equilibration..."
2143.     Appendlinesandnetwork()
2144.     runOptimization()
2145.     wave test0
2146.     passmat=test0
2147.     UpdateVariables()
2148.     updateGraph()
2149.     updateColors()
2150.     Saveit("Unfolding","Time_0")
2151.     for(i=0;i<Duration/(dt*1e3)-1;i+=1)
2152.         updateflag=0
2153.         theta_curr[]=theta_list[dimsize(theta_list,0)-1][p]
2154.         if (i>0)
2155.             theta_past[]=theta_list[dimsize(theta_list,0)-2][p]
2156.         endif
2157.         for(l=0;l<N_lines;l+=1)
2158.             WAVE Wx=\$("root:results:x_"+num2str(l))
2159.             WAVE Wf=\$("root:results:F_"+num2str(l))
2160.             WAVE Wn=\$("root:results:N_"+num2str(l))
2161.             wave vlevel=\$("root:results:vlevel_"+num2str(l))
2162.             if (i>0)
2163.                 currentF=round(commandF[i]*abs(cos(theta_curr[l])))
2164.                 oldF=round(commandF[i-1]*abs(cos(theta_past[l])))
2165.             else
2166.                 currentF=round(commandF[i]*abs(cos(theta_curr[l])))
2167.             endif
2168.             if (currentF<2)
2169.                 currentF=2
2170.             endif
2171.             if (oldF<2)
2172.                 OlDF=2
2173.             endif
2174.             if(i>0)
2175.                 if (oldF!=CurrentF)
2176.                     OneD[]=MinXYF[p][0](oldF)
2177.                     xc=wx[i]
2178.                     findlevel/Q/P/R=[0,dimSize(MINXYF,0)] OneD,xc// which E-curve were you on at i-1?
2179.                     if(V_flag)  //  sometimes at full unfolding/folding findlevel gives NAN
2180.                         Wavestats/Q OneD
2181.                         if(xc>V_max)
2182.                             V_LevelX = N_modules +1
2183.                         endif
2184.                         if (xc<V_min)
2185.                             V_LevelX = 1
2186.                         endif
2187.                     endif
2188.                     OneD[]=MinXYF[p][0](currentF)
2189.                     xc=OneD[round(V_LevelX)]// move to that E-curve at the new force
2190.                 else
2191.                     OneD[]=MinXYF[p][0](currentF)
2192.                     xc=wx[i]
2193.                 endif
2194.             else
2195.                 OneD[]=MinXYF[p][0](currentF)
2196.                 xc=MinXYF[1][0](currentF)
2197.             endif
2198.             if(currentF<20)
2199.                 D=18450
2200.             else
2201.                 D=D0
2202.             endif
2203.             if(i>0)
2204.                 xcold=Wx[numpnts(Wx)-1]
2205.                 if (xcold<deltax(PMFM))
2206.                     xcold=deltax(PMFM)
2207.                 endif
2208.                 for(j=0;j<1e3;j+=1)
2209.                     do
2210.                         xc=xcold
2211.                         tc+=dt
2212.                         F_b = D_PMFM(xc)(currentF)
2213.                         Fr = sqrt(2*D*dt)*gnoise(1)
2214.                         xc +=  -dt*D * F_b/kT   + Fr
2215.                     while ((xc<deltax(PMFM)))
2216.                     //Correcting for tunneling at large x
2217.                     if (xc!=xc)
2218.                         xc=xcold
2219.                     endif
2220.                     xcold=xc
2221.                 endfor
2222.             endif
2223.             Findlevel/Q/P/R=[0,dimSize(MINXYF,0)] OneD,xc
2224.             if(V_flag)  //  sometimes at full unfolding/folding findlevel gives NAN
2225.                 V_LevelX = Wn[i]+1
2226.             endif
2227.             current_E_var = round(V_LevelX)-1
2228.             if (current_E_var<0)
2229.                 current_E_var=0
2230.             endif
2231.             if (i>0)
2232.                 state[l][0]+=1
2233.                 //if state changes reset time and note direction in state[l][1]
2234.                 if (current_e_var>wn[numpnts(wn)-1])
2235.                     //this is trying to unfold scenario
2236.                     if ((state[l][1]==0) && (state[l][0]>=999))
2237.                         // trying to unfold, post refolding, after at least 10 ms
2238.                         state[l][0]=0;state[l][0]=1
2239.                     elseif ((state[l][1]==0) && (state[l][0]<999))
2240.                         // if not enough time has passed for an unfolding, reset state back to
2241.                         // last local minimum and adjust current_e_var
2242.                         current_e_var=wn[numpnts(wn)-1]
2243.                         OneD[]=MinXYF[p][0](currentF)
2244.                         xc=OneD[current_e_var+1]
2245.                     elseif (state[l][1]==1)
2246.                         //if it previously unfolding and trying to unfold again, let it. and reset
2247.                         //state[l][0] to zero
2248.                         state[l][0]=0
2249.                     endif
2250.                 endif
2251.                 if (current_e_var<wn[numpnts(wn)-1])
2252.                     //this is trying to refold now
2253.                     if ((state[l][1]==1) && state[l][0]>=999)
2254.                         //refolding after at least 10 ms of being unfolded
2255.                         state[l][0]=0;state[l][1]=0
2256.                     elseif ((state[l][1]==1) && (state[l][0]<999))
2257.                         // if not enough time has passed for an refolding, reset state back to
2258.                         // last local minimum and adjust current_e_var
2259.                         current_e_var=wn[numpnts(wn)-1]
2260.                         OneD[]=MinXYF[p][0](currentF)
2261.                         xc=OneD[current_e_var+1]
2262.                     elseif (state[l][1]==0)
2263.                         //refolding after past refolding is fine, no need to change state[l][1]
2264.                         state[l][0]=0
2265.                     endif
2266.                 endif
2267.             endif
2268.             InsertPoints numpnts(Wx),1, Wx, Wf, Wn
2269.             Wx[numpnts(Wx)-1]=xc
2270.             Wf[numpnts(Wf)-1]=currentF+Fr
2271.             Wn[numpnts(Wn)-1]=round(current_E_var)
2272.             //If there is an unfolding, run optimization
2273.             if (i>0)
2274.                 if (Wn[numpnts(Wn)-1]!=Wn[numpnts(Wn)-2])
2275.                     updateflag=1
2276.                 endif
2277.             endif
2278.             if(numpnts(Wx)>1)
2279.                 deltawave[l]=Wx[numpnts(Wx)-1]-Wx[numpnts(Wx)-2]
2280.             else
2281.                 deltawave[l]=Wx[numpnts(Wx)-1]
2282.             endif
2283.         endfor
2284.         AppendLinesAndNetwork()
2285.         if (updateflag)
2286.             Print "Working on Optimization..."
2287.             sleep/s 0.01
2288.             runOptimization()
2289.             wave test0
2290.             passmat=test0
2291.             //Optimize/A=0/X=passmat/I=15/M={0,0}/T={0.01,0.001}/Q EnergyCalc,passmat
2292.         endif
2293.         UpdateVariables()
2294.         updateGraph()
2295.         updateColors()
2296.         updateUnfoldHistory()
2297.         print num2str(100*i/(Duration/(dt*1e3)-1))+" % Done"
2298.         if ((mod(i,100)==0) && Displayw)
2299.             Saveit("Unfolding","Time_"+num2str(i+1))
2300.             sleep/s 0.01
2301.         endif
2302.         //calc_gellength_ionel()
2303.     endfor
2304.     avgsFixes()
2305.     print "Finished"
2306. End
2307. Function/S UnixShellCommand()
2308.     // Paths must be POSIX paths (using /).
2309.     // Paths containing spaces or other nonstandard characters
2310.     // must be single-quoted. See Apple Techical Note TN2065 for
2311.     // more on shell scripting via AppleScript.
2312.     String unixCmd
2313.     String igorCmd
2314.     svar mainpath
2315.     string folder=parsefilepath(5,mainpath,"/",0,0)+"/Optimization"
2316.     unixCmd = "cd '"+folder+"'; sh 'doOptim.sh' > /dev/null 2>&1"
2317.     sprintf igorCmd, "do shell script \"%s\"", unixCmd
2318.     ExecuteScriptText/UNQ igorCmd
2319.     Print S_value       // For debugging only.
2320.     return S_value
2321. End
2322. function runOptimization()
2323.     wave test0
2324.     nvar isWindows
2325.     svar OptimizationPath
2326.     wave status0
2327.     status0[0]=0
2328.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath status0 as "status.txt"
2329.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:results:networkmat as "priorMat.txt"
2330.     Save/J/M="\n"/DLIM=","/O/P=OptimizationPath root:results:passmat as "passMat.txt"
2331.     if (isWindows)
2332.         windowsShell()
2333.     else
2334.         UnixShellCommand()
2335.     endif
2336.     do
2337.         try
2339.             AbortonRTE
2340.         catch
2341.             variable err=getRTError(1)
2342.         endtry
2343.         sleep/s 0.2
2344.     while (status0[0]==0)
2346. end
2347. function windowsShell()
2348.     svar mainpath
2349.     string folder=parsefilepath(5,mainpath,"\\",0,0)+"\\Optimization"
2350.     ExecuteScriptText "cmd.exe /K \""+folder+"\\optim.bat "+folder+"\\\""
2351. end
2352. function runSim(x,y,z,conc,basepath,newFoldername)
2353.     // Basepath folder must have 'errorFunc.m', 'findParams.m', 'optim.m', and 'setUp.m' inside of it
2354.     variable x,y,z,conc
2355.     string basepath,newFoldername
2356.     variable/G isWindows=1
2357.     string mainfolder=parsefilepath(5,basepath,"\\",0,0)
2358.     ExecuteScriptText "cmd.exe /K \""+mainfolder+"\\setUp.bat "+mainfolder+" "+newFoldername+"\\\""
2359.     string simfolder=basepath+newFoldername
2360.     make_lines(x,y,z,conc,simfolder)
2361.     thegel()
2362.     newpath/o savepath simfolder
2363.     saveexperiment/p=savepath/F={1,"",2} as newfoldername+".pxp"
2364. end
2365. function runSimMac(x,y,z,conc,basepath,newFoldername)
2366.     // Basepath folder must have 'errorFunc.m', 'findParams.m', 'optim.m', and 'setUp.m' inside of it
2367.     variable x,y,z,conc
2368.     string basepath,newFoldername
2369.     variable/G isWindows=0
2370.     String unixCmd
2371.     String igorCmd
2372.     string mainfolder=parsefilepath(5,basepath,"/",0,0)
2373.     unixCmd = "cd '"+mainfolder+"'; sh 'setUp.sh' "+mainfolder+" "+newfoldername
2374.     sprintf igorCmd, "do shell script \"%s\"", unixCmd
2375.     ExecuteScriptText/UNQ igorCmd
2376.     string simfolder=basepath+newFoldername
2377.     make_lines(x,y,z,conc,simfolder)
2378.     thegel()
2379.     newpath/o savepath simfolder
2380.     saveexperiment/p=savepath/F={1,"",2} as newfoldername+".pxp"
2381. end
2382. function updateUnfoldHistory()
2383.     wave unfoldmat_history=root:results:unfoldmat_history
2384.     wave unfoldmat=root:results:unfoldmat
2385.     variable i,j
2386.     nvar modules,n_lines
2387.     insertpoints/M=2 dimsize(unfoldmat_history,2),1,unfoldmat_history
2388.     for(i=0;i<n_lines;i+=1)
2389.         wave sizew=\$("root:results:sizew_"+num2str(i))
2390.         for(j=0;j<modules;j+=1)
2391.             if (sizew[j][0]==0.2)
2392.                 unfoldmat_history[j][i]=0
2393.             elseif (sizew[j][0]==0.01)
2394.                 unfoldmat_history[j][i][dimsize(unfoldmat_history,2)-1]=1
2395.             endif
2396.         endfor
2397.     endfor
2398. end
2399. function dotwo()
2400.     variable i,j,f,c
2401.     runSim(4,4,2,0.0009,"C:Users:phy-popalab-g:Desktop:Kirill:Important:Tests:","oldSim")
2402. end
2403. #pragma TextEncoding = "UTF-8"
2404. #pragma rtGlobals=3     // Use modern global access method and strict wave access.
2405. function run_all_length_analysis()
2406.     iter()
2407.     GelLength_hists()
2408.     //get_x_avg()
2409. end
2410. function test_lengths(name,n_lines)
2411.     string name
2412.     variable N_lines
2413.     setdatafolder name
2414.     variable i,k
2415.     variable modules=8
2416.     make/o/n=(N_lines*modules,dimsize(\$("LineHistory_0"),2)) z_Mat
2417.     for(i=0;i<N_lines;i+=1)
2418.         wave LineHistory=\$("LineHistory_"+num2str(i))
2419.         for(k=i*modules;k<(i+1)*modules;k+=1)
2420.                 z_Mat[k][]=linehistory[k-i*modules][2][q]
2421.         endfor
2422.         make/o/n=(dimsize(LineHistory,2)-1) \$("X_Line_"+num2str(i))
2423.         wave x_wave=\$("X_Line_"+num2str(i))
2424.         wave killwave=\$("X_"+num2str(i))
2425.         killwaves killwave
2426.         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)
2427.     endfor
2428. end
2429. function iter()
2430.     variable F,N,i,minF,maxF,minN,maxN
2431.     minF=30;MaxF=60;minN=2;maxN=6
2432.     for(F=MinF;F<(MaxF+1);F+=10)
2433.         if (F==30)
2434.             for(N=5;N<(maxN+1);N+=1)
2435.                 for(i=0;i<5;i+=1)
2436.                     test_lengths("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i),N*N*2)
2437.                 endfor
2438.             endfor
2439.         else
2440.             for(N=minN;N<(MaxN+1);N+=1)
2441.                 for(i=0;i<5;i+=1)
2442.                     test_lengths("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i),N*N*2)
2443.                 endfor
2444.             endfor
2445.         endif
2446.     endfor
2447. end
2448. function top_n()
2449.     variable F,N,i,j,k
2450.     variable modules=8
2451.     variable n_avg=3,max_pnts,min_pnts
2452.     string namew="Top_5_gel_length"
2453.     for(F=50;F<60;F+=10)
2454.         for(N=5;N<6;N+=1)
2455.             for(i=0;i<1;i+=1)
2456.                 setdatafolder "root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2457.                 make/o/n=(N*N*2*modules,50002) z_Mat
2458.                 n_avg=5
2459.                 for(j=0;j<(N*N*2);j+=1)
2460.                     wave lineHis=\$("LineHistory_"+num2str(j))
2461.                     for(k=j*modules;k<(j+1)*modules;k+=1)
2462.                         z_Mat[k][]=lineHis[k-j*modules][2][q]
2463.                     endfor
2464.                 endfor
2465.                 make/o/n=(dimsize(z_mat,1)) \$(namew)
2466.                 wave curr=\$(namew)
2467.                 make/o/n=(N*N*2*modules) curr_z
2468.                 for(j=0;j<dimsize(z_mat,1);j+=1)
2469.                     curr_z[]=z_mat[p][j]
2470.                     sort/R curr_z,curr_z
2471.                     max_pnts=0;min_pnts=0
2472.                     for(k=0;k<n_avg;k+=1)
2473.                         max_pnts+=curr_z[k]/n_avg
2474.                         min_pnts+=curr_z[N*N*2*modules-1-k]/n_avg
2475.                     endfor
2476.                     curr[j]=max_pnts-min_pnts
2477.                 endfor
2478.                 killwaves/Z curr_z
2479.                 print "Done with F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2480.             endfor
2481.         endfor
2482.     endfor
2483.     get_avgs(namew)
2484. end
2485. function GelLength_hists()
2486.     variable F,N,i,j,k,minF,maxF,minN,maxN
2487.     variable modules=8
2488.     variable n_avg=3,max_pnts,min_pnts
2489.     minF=30;MaxF=60;minN=2;maxN=6
2490.     for(F=MinF;F<(MaxF+1);F+=10)
2491.         if (F==30)
2492.             for(N=5;N<(MaxN+1);N+=1)
2493.                 for(i=0;i<5;i+=1)
2494.                     setdatafolder "root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2495.                     wave z_mat
2496.                     print "Starting F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2497.                     get_hist_length(z_mat,n*N*2)
2498.                     print "Done with F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2499.                 endfor
2500.             endfor
2501.         else
2502.             for(N=MinN;N<(MaxN+1);N+=1)
2503.                 for(i=0;i<5;i+=1)
2504.                     setdatafolder "root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2505.                     wave z_mat
2506.                     print "Starting F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2507.                     get_hist_length(z_mat,n*N*2)
2508.                     print "Done with F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)
2509.                 endfor
2510.             endfor
2511.         endif
2512.     endfor
2513.     get_avgs("hist_Gel_Length")
2514. end
2515. function get_hist_length(w,n_lines)
2516.     wave w
2517.     variable n_lines
2518.     variable i,j,k,l,r,q,A,s,t,binsize
2519.     variable zbox=2/(0.0009)^(1/3)
2520.     make/o/n=(dimsize(w,1)) hist_Gel_Length=0
2521.     wave len=hist_Gel_Length
2522.     r=zbox
2523.     q=2
2524.     s=0
2525.     t=N_lines/5
2526.     Make/O/T/N=1 T_Constraints
2527.     binsize=4//11*exp(-0.015*N_lines)//Value-Tested
2528.     T_Constraints[0] = {"K3 > 0","K4 < 100","K4 > 1"}
2529.     for(i=0;i<(dimsize(w,1));i+=1)
2530.         make/o/n=(dimsize(w,0)) newz
2531.         newz[]=w[p][i]
2532.         wavestats/Q newz
2533.         Make/O newz_hist
2534.         Histogram/C/N/B={v_min-20,binsize,(v_max-v_min+40)/binsize} newz,newz_hist
2535.         wavestats/q newz_hist
2536.         T_Constraints[0] = {"K3 > "+num2str(r-zbox/10),"K3 < "+num2str(r+zbox/10),"K4 < 100","K4 > 1"}//,"K3 < "+num2str(v_maxloc-v_minloc)}
2537.         Make/D/N=5/O W_coef
2538.         W_coef[0] = {0,t,s,r,q}
2539.         FuncFit/H="10000"/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
2540.         //curvefit/x=1/Q gauss newz_hist /D
2541.         wave w_coef,w_sigma
2542.         len[i]=abs(w_coef[3])
2543.         r=abs(w_coef[3])
2544.         q=abs(w_coef[4])
2545.         if (q>20)
2546.             q=5
2547.         endif
2548.         t=abs(w_coef[1])
2549.         s=w_coef[2]
2550.         //print num2str(i/(dimsize(w,1)-1)*100)+" % Done"
2551.     endfor
2552.     killwaves/Z T_constraints, w_coef, W_sqrtN, newz,newz_hist,fit_newz_hist
2553. end
2554. function Extended_Gel_Length(w,n_lines)
2555.     wave w
2556.     variable n_lines
2557.     variable i,j,k,l,r,q,A,s,t,binsize,maxz,minz
2558.     variable zbox=2/(0.0009)^(1/3)
2559.     variable modules=8,interval=1
2560.     make/o/n=(dimsize(w,1)-1) hist_Extended_Gel_Length=0
2561.     wave len=hist_Extended_Gel_Length
2562.     r=zbox
2563.     q=2
2564.     s=0
2565.     t=N_lines/5
2566.     for(i=0;i<(dimsize(w,1)-1);i+=1)
2567.         make/o/n=0 newz
2568.         for(j=0;j<(sqrt(N_lines/2));j+=1)
2569.             wave LineHis=\$("LineHistory_"+num2str(j))
2570.             maxz=max(LineHis[0][2][i],LineHis[modules-1][2][i]);minz=min(LineHis[0][2][i],LineHis[modules-1][2][i])
2571.             Make/O/N=2 A_wave; a_wave={minz,maxz}
2572.             variable N_pnts=round((maxz-minz)/interval)
2573.             if (N_pnts<2)
2574.                 N_pnts=2
2575.             endif
2576.             Interpolate2/T=1/N=(N_pnts)/Y=z_pnts a_wave
2577.             duplicate/o newz,temp
2578.             wave temp
2579.             concatenate/NP/O {temp,z_pnts}, newz
2580.         endfor
2581.         Make/O/T/N=1 T_Constraints
2582.         T_Constraints[0] = {"K3 > 0","K4 < 100","K4 > 1"}
2583.         wavestats/Q newz
2584.         binsize=3.49*v_sdev/(numpnts(newz)^(1/3))
2585.         Make/O newz_hist
2586.         Histogram/C/N/B={v_min-21,binsize,(v_max-v_min+38)/binsize} newz,newz_hist
2587.         wavestats/q newz_hist
2588.         T_Constraints[0] = {"K3 > "+num2str(r-zbox/10),"K3 < "+num2str(r+zbox/10),"K4 < 100","K4 > 1"}//,"K3 < "+num2str(v_maxloc-v_minloc)}
2589.         Make/D/N=5/O W_coef
2590.         W_coef[0] = {0,t,s,r,q}
2591.         FuncFit/H="10000"/X=1/Q SuperGauss W_coef newz_hist /D /C=T_Constraints
2592.         //curvefit/x=1/Q gauss newz_hist /D
2593.         wave w_coef,w_sigma
2594.         len[i]=abs(w_coef[3])
2595.         r=abs(w_coef[3])
2596.         q=abs(w_coef[4])
2597.         if (q>20)
2598.             q=5
2599.         endif
2600.         t=abs(w_coef[1])
2601.         s=w_coef[2]
2602.         //print num2str(i/(dimsize(w,1)-1)*100)+" % Done"
2603.     endfor
2604.     killwaves/Z T_constraints, w_coef, W_sqrtN, newz,newz_hist,fit_newz_hist,temp,z_pnts,A_wave
2605. end
2606. function get_x_avg()
2607.     variable F,N,i,j,k,minF,maxF,minN,maxN
2608.     minF=30;MaxF=60;minN=2;maxN=6
2609.     for(F=minF;F<(MaxF+1);F+=10)
2610.         if (F==30)
2611.             for(N=5;N<(MaxN+1);N+=1)
2612.                 for(i=0;i<5;i+=1)
2613.                     wave x_0=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_0")
2614.                     duplicate/o x_0,\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
2615.                     wave xavg=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
2616.                     for(j=1;j<(N*N*2);j+=1)
2617.                         wave xwave=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":X_"+num2str(j))
2618.                         xavg+=xwave
2619.                     endfor
2620.                     xavg/=(n*N*2)
2621.                 endfor
2622.             endfor
2623.         else
2624.             for(N=minN;N<(MaxN+1);N+=1)
2625.                 for(i=0;i<5;i+=1)
2626.                     wave x_0=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_0")
2627.                     duplicate/o x_0,\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
2628.                     wave xavg=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":x_avg")
2629.                     for(j=1;j<(N*N*2);j+=1)
2630.                         wave xwave=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":X_"+num2str(j))
2631.                         xavg+=xwave
2632.                     endfor
2633.                     xavg/=(n*N*2)
2634.                 endfor
2635.             endfor
2636.         endif
2637.     endfor
2638. end
2639. function get_avgs(namew)
2640.     string namew
2641.     variable i,F,N,j,minF,maxF,minN,maxN
2642.     minF=30;maxf=60;minN=2;maxN=6
2643.     for(F=minF;F<(MaxF+1);F+=10)
2644.         if (F==30)
2645.             for(N=5;N<(maxN+1);N+=1)
2646.                 for(i=0;i<5;i+=1)
2647.                     wave inwave=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":"+namew)
2648.                     if (i==0)
2649.                         duplicate/o inwave, \$("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
2650.                         wave avg=\$("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
2651.                     else
2652.                         avg+=inwave
2653.                     endif
2654.                 endfor
2655.                 avg/=5
2656.             endfor
2657.         else
2658.             for(N=minN;N<(MaxN+1);N+=1)
2659.                 for(i=0;i<5;i+=1)
2660.                     wave inwave=\$("root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:Gel_"+num2str(i)+":"+namew)
2661.                     if (i==0)
2662.                         duplicate/o inwave, \$("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
2663.                         wave avg=\$("Root:F_"+num2str(F)+":N_"+num2str(N)+"x"+num2str(N)+"x2:avg_"+namew)
2664.                     else
2665.                         avg+=inwave
2666.                     endif
2667.                 endfor
2668.                 avg/=5
2669.             endfor
2670.         endif
2671.     endfor
2672. end
2673. function normalize()
2674.     variable start,last,i,maxN,minN,j
2675.     minN=2;maxN=6
2676.     for(i=minN;i<=maxN;i+=1)
2677.         wave wavg=\$("N_"+num2str(i)+"_avg")
2678.         start=wavemin(wavg)
2679.         wavg-=start
2680.         last=wavg[inf]
2681.         wavg/=last
2682.         setscale/i x,0,5,"s",wavg
2683.         for(j=0;j<5;j+=1)
2684.             wave w=\$("N_"+num2str(i)+"_"+num2str(j))
2685.             start=wavemin(w)
2686.             w-=start
2687.             last=w[inf]
2688.             w/=last
2689.             setscale/i x,0,5,"s",w
2690.         endfor
2691.     endfor
2692.
2693.
2694. end
2695. function calc_resid()
2696.     variable N,i
2697.     make/o/n=5 resid_wave
2698.     make/o/n=5 N_wave
2699.     for(N=2;N<=6;N+=1)
2700.         duplicate/o \$("N_"+num2str(N)+"_avg") \$("resid_"+num2str(N))
2701.         wave avg=\$("N_"+num2str(N)+"_avg")
2702.         wave resid=\$("resid_"+num2str(N))
2703.         for(i=0;i<5;i+=1)
2704.             wave curr=\$("N_"+num2str(N)+"_"+num2str(i))
2705.             resid[]+=abs(avg[p]-curr[p])/5
2706.         endfor
2707.         N_wave[N-2]=n*n*2;resid_wave[N-2]=sum(resid)/numpnts(resid)
2708.     endfor
2709. end
2710.
2711. function quick_resid(num)
2712.     variable num
2713.     variable i
2714.     wave avg=\$("N_"+num2str(num)+"_avg")
2715.     for(i=0;i<5;i+=1)
2716.         wave orig=\$("N_"+num2str(num)+"_"+num2str(i))
2717.         duplicate/o orig, \$("res_"+"N_"+num2str(num)+"_"+num2str(i))
2718.         wave res=\$("res_"+"N_"+num2str(num)+"_"+num2str(i))
2719.         res=avg-res
2720.         if (i==0)
2721.             display res
2722.         else
2723.             appendtograph res
2724.         endif
2725.     endfor
2726. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top