daily pastebin goal
91%
SHARE
TWEET

Untitled

a guest May 16th, 2018 119 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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top