daily pastebin goal
90%
SHARE
TWEET

Untitled

a guest May 16th, 2018 117 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
  218. Menu "SIM"
  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)
  572.     string/G CrosslinkingPath,UnfoldingPath,OptimizationPath,mainpath
  573.     mainpath=savepath
  574.     NewPath/C/O/Q CrosslinkingPath savepath+":Crosslinking"
  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)
  680.         Saveit("Crosslinking","time_"+num2str(savetime))
  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)
  696.             Saveit("Crosslinking","time_"+num2str(savetime))
  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()
  726.     Saveit("Crosslinking","time_"+num2str(savetime))
  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
  763.     Saveit("Crosslinking","time_"+num2str(savetime))
  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()
  934.     nvar n_lines,modules,box_size,numberpoints,connectPerClust,LinkSize
  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
  942.     Make/O/N=(0,3) root:RESULTS:Crosslinks
  943.     WAVE Crosslinks=root:RESULTS:Crosslinks
  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
  958.                         InsertPoints/M=0 dimsize(Crosslinks,0),2,Crosslinks
  959.                         Crosslinks[dimsize(Crosslinks,0)-2][0]=L1[k][0]
  960.                         Crosslinks[dimsize(Crosslinks,0)-2][1]=L1[k][1]
  961.                         Crosslinks[dimsize(Crosslinks,0)-2][2]=L1[k][2]
  962.                         Crosslinks[dimsize(Crosslinks,0)-1][0]=L2[l][0]
  963.                         Crosslinks[dimsize(Crosslinks,0)-1][1]=L2[l][1]
  964.                         Crosslinks[dimsize(Crosslinks,0)-1][2]=L2[l][2]
  965.                         InsertPoints/M=0 dimsize(Crosslinks,0),1,Crosslinks
  966.                         Crosslinks[dimsize(Crosslinks,0)-1][]=NAN                  
  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
  1243.     nvar n_lines,modules,linksize
  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
  1280.     tempwave[]=(criticaldis[p]-linksize)^2
  1281.     energy=sum(tempwave)
  1282.     //for(i=0;i<numpnts(criticaldis);i+=1)
  1283.     //  energy+=(criticaldis[i]-LinkSize)^2
  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
  1421.     Make/O/N=(0,3) root:RESULTS:Crosslinks
  1422.     WAVE Crosslinks=root:RESULTS:Crosslinks
  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]))
  1431.         InsertPoints/M=0 dimsize(Crosslinks,0),2,Crosslinks
  1432.         Crosslinks[dimsize(Crosslinks,0)-2][0]=w1[movement[i][1]][0]
  1433.         Crosslinks[dimsize(Crosslinks,0)-2][1]=w1[movement[i][1]][1]
  1434.         Crosslinks[dimsize(Crosslinks,0)-2][2]=w1[movement[i][1]][2]
  1435.         Crosslinks[dimsize(Crosslinks,0)-1][0]=w2[movement[i][3]][0]
  1436.         Crosslinks[dimsize(Crosslinks,0)-1][1]=w2[movement[i][3]][1]
  1437.         Crosslinks[dimsize(Crosslinks,0)-1][2]=w2[movement[i][3]][2]
  1438.         InsertPoints/M=0 dimsize(Crosslinks,0),1,Crosslinks
  1439.         Crosslinks[dimsize(Crosslinks,0)-1][]=NAN          
  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
  1524.     svar crosslinkingpath,unfoldingpath
  1525.     if (displayw)
  1526.         if (stringmatch(path,"Crosslinking"))
  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
  1703.     variable i,j,k,r,flag,mask,avg,low
  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)
  1726.         mask=1
  1727.         low=1e10
  1728.         for(j=startp;j<endp;j+=1)
  1729.             if(w[j]>w[startp])
  1730.                 mask=0
  1731.             endif
  1732.             if(w[j]<low)
  1733.                 low=w[j]
  1734.             endif
  1735.         endfor
  1736.         if (mask==1)
  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))
  1838.             matrixop/o LinkEnergy=(DeltaU.UnitOld)
  1839.             energy+=(StrechPrefactor/oldveclength[0])*(LinkEnergy[0])^2
  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]))
  1861.         matrixop/o LinkEnergy=(DeltaU.UnitOld)
  1862.         energy+=(StrechPrefactor/oldveclength[0])*(LinkEnergy[0])^2
  1863.         //Claculate from link distances
  1864.         Matrixop/o LengthVector=row(passmat,clindex[i][0])-row(passmat,clindex[i][1]))
  1865.         matrixop/o LinkLength=powr(sumsqr(LengthVector),0.5)
  1866.         energy+=CrossLinkPrefactor*(LinkLength[0]-LinkSize)^2
  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
  1910.                 matrixop/o linknet=row(networkmat,j)-row(networkmat,j+1)
  1911.                 matrixop/o linkpass=row(passmat,j)-row(passmat,j+1)
  1912.                 deltheta+=(atan(sqrt((linkpass[0][0])^2 + (linkpass[0][1])^2)/linkpass[0][2])-atan(sqrt((linknet[0][0])^2 + (linknet[0][1])^2)/linknet[0][2]))/(domaincols[0][i]-1)
  1913.                 delphi+=(atan(linkpass[0][1]/linkpass[0][0])-atan(linknet[0][1]/linknet[0][0]))/(domaincols[0][i]-1)
  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
  2100.     //removecrosslinks()
  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
  2338.             LoadWave/G/O/J/n=status/p=OptimizationPath/M/Q/Z "status.txt"
  2339.             AbortonRTE
  2340.         catch
  2341.             variable err=getRTError(1)
  2342.         endtry
  2343.         sleep/s 0.2
  2344.     while (status0[0]==0)
  2345.     LoadWave/G/O/J/n=test/p=OptimizationPath/M/Q "optimParams.txt"
  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
Top