source: trunk/User Procedures/Irena/IR3_3DModels.ipf @ 1169

Last change on this file since 1169 was 1169, checked in by ilavsky, 10 months ago

Remove debug code, never used. Fix Nika minor bugs

  • Property svn:eol-style set to native
File size: 180.0 KB
Line 
1#pragma TextEncoding = "UTF-8"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3#pragma version=1.01
4
5Constant IR3AMassFrAggVersionNumber     = 1.02
6Constant IR3TPOVPDBVersionNumber                = 1.00
7Constant IR3TTwoPhaseVersionNumber      = 1.00
8
9
10//*************************************************************************\
11//* Copyright (c) 2005 - 2021, Argonne National Laboratory
12//* This file is distributed subject to a Software License Agreement found
13//* in the file LICENSE that is included with this distribution.
14//*************************************************************************/
15
16
17//1.02 added MultiParticleAttraction parameter.
18//1.01 3dAggregate added ability ot grow N aggregates and Compare Stored graph
19//1.00 first version, added code for 3DMassFractalAggregate from Alex McGlassson
20//                      note: this ipf file also contains tools for import of pdb (ATSAS, GNOM produced output files) and for POV files produced by SAXSMorph.
21//                      and some Gizmo tools to be used to visualize these.
22
23
24
25//   Parameters description for referecne... This belongs to manual.
26//      MultiParticleAttraction = "Multi Part attr" controls how particle attached when approaching existing aggregate.
27// Sticking method controls how particle attached when approaching existing aggregate
28// MultiParticleAttraction = "Neutral;Attractive;Repulsive;Not Allowed;"
29// when : neutral, probablity of attaching does not depend on number of particles in neaest neighbor sphere aroudn the new position.
30// when : Attractive more particles increase the probability of attaching
31// when : Repulsive more particles decrease the probability of attaching.
32// consequence - Repulsive creases larger, more open particles, Attractive creates more compact particles
33// Sticking method = which neighbors are counted as "nearest neigbor"
34//      1 : only in x, y, znd z direction on lattice, their center distance is < 1.1
35//      2 : in x,z,y and in xy, yz, xz planes, their distacne is < 1.05*sqrt(2)
36//      3 : also in xyz body diagonal, their distacne is < 1.05*sqrt(3)
37// to grow compact particle set sticking method 1, low sticking probability and positive attraction, I got df up to 2.55   
38// to grow open particle, set sticking method 3, high sticking probability and negative attraction, I got df below 2 this way. 
39//
40//      SVAR MultiParticleAttraction = root:Packages:AggregateModeling:MultiParticleAttraction
41//      variable StickingProbability1, StickingProbabilityM1,StickingProbabilityM2, StickingProbabilityLoc
42//      if(StringMatch(MultiParticleAttraction, "Neutral"))
43//              StickingProbability1= StickingProbability
44//              StickingProbabilityM1= StickingProbability
45//              StickingProbabilityM2= StickingProbability
46//      elseif(StringMatch(MultiParticleAttraction, "Attractive"))
47//              StickingProbability1= StickingProbability
48//              StickingProbabilityM1= (StickingProbability+100)/2
49//              StickingProbabilityM2= (StickingProbability+300)/4
50//      elseif(StringMatch(MultiParticleAttraction, "Repulsive"))
51//              StickingProbability1 = StickingProbability
52//              StickingProbabilityM1 = (StickingProbability+10)/2
53//              StickingProbabilityM2 = (StickingProbability+30)/4
54//      elseif(StringMatch(MultiParticleAttraction, "Not allowed"))
55//              StickingProbability1 = StickingProbability
56//              StickingProbabilityM1 = 1
57//              StickingProbabilityM2 = 0
58//      else
59//              StickingProbability1 = StickingProbability
60//              StickingProbabilityM1 = StickingProbability
61//              StickingProbabilityM2 = StickingProbability
62//      endif
63
64
65
66//******************************************************************************************************************************************************
67//******************************************************************************************************************************************************
68//                      3D packages, 2018-12-26
69//******************************************************************************************************************************************************
70//******************************************************************************************************************************************************
71//******************************************************************************************************************************************************
72//                      Main packages as they are called from main menu.
73//******************************************************************************************************************************************************
74//******************************************************************************************************************************************************
75Function IR3A_MassFractalAggregate()
76                //this calls Fractal aggregate controls.
77        DoWIndow FractalAggregatePanel
78        if(V_Flag)
79                DoWIndow/K FractalAggregatePanel
80        endif
81        IN2G_CheckScreenSize("height",670)
82        IR3A_InitializeMassFractAgg()
83        IR3A_FractalAggregatePanel()
84        ING2_AddScrollControl()
85        IR1_UpdatePanelVersionNumber("FractalAggregatePanel", IR3AMassFrAggVersionNumber,1)
86        IR3A_UpdatePanelValues()       
87
88end
89//******************************************************************************************************************************************************
90//******************************************************************************************************************************************************
91//******************************************************************************************************************************************************
92//******************************************************************************************************************************************************
93Function IR3P_ImportPOVPDB()
94                //this calls GUI for import of various 3D formats - POV from SAXSMorph and PDB files from GNOM.
95                //this reads POV files: IR3T_ReadPOVFile(dimx,dimy,dimz)
96                //for pdb there is code in Subversion which needs to be added here...
97                //DoAlert /T="Unfinished" 0, "IR3A_ImportPOVPDB() is not finished yet"
98                //this calls POV/PDB import controls.
99        DoWIndow/K/Z POVPDBPanel
100        DoWIndow/K/Z POV3DData
101
102        IN2G_CheckScreenSize("height",670)
103        IR3P_InitializePOVPDB()
104        IR3P_POVPDBPanel()
105        ING2_AddScrollControl()
106        IR1_UpdatePanelVersionNumber("POVPDBPanel", IR3TPOVPDBVersionNumber,1)
107end
108//******************************************************************************************************************************************************
109//******************************************************************************************************************************************************
110Function IR3A_Display3DData()
111                //this calls controls which can display 3D data generated by other 3D code or imported by POV/PDB importer
112                //use Gizmo, add some way of selecting data and do some simple tricks...
113                DoAlert /T="Unfinished" 0, "IR3A_Display3DData() is not finished yet"
114                //this may be useful: IR3T_Convert3DMatrixToList(My3DVoxelWave, threshVal)
115end
116//******************************************************************************************************************************************************
117//                      Utility functions
118//******************************************************************************************************************************************************
119//*****************************************************************************************************************
120//*****************************************************************************************************************
121//*****************************************************************************************************************
122Function IR3T_MainCheckVersion()       
123        //this needs to get more of these lines for each tool/panel...
124        DoWindow TwoPhaseSystems
125        if(V_Flag)
126                if(!IR1_CheckPanelVersionNumber("TwoPhaseSystems", IR3TTwoPhaseVersionNumber))
127                        DoAlert /T="The Two Phase 3D modeling panel was created by incorrect version of Irena " 1, "Two Phase 3D modeling tool may need to be restarted to work properly. Restart now?"
128                        if(V_flag==1)
129                                DoWindow/K TwoPhaseSystems
130                                IR3T_TwoPhaseSystem()
131                        else            //at least reinitialize the variables so we avoid major crashes...
132                                IR3T_InitializeTwoPhaseSys()
133                        endif
134                endif
135        endif
136end
137//*****************************************************************************************************************
138//*****************************************************************************************************************
139Function IR3P_MainCheckVersion()       
140        //this needs to get more of these lines for each tool/panel...
141        DoWindow POVPDBPanel
142        if(V_Flag)
143                if(!IR1_CheckPanelVersionNumber("POVPDBPanel", IR3TPOVPDBVersionNumber))
144                        DoAlert /T="The POV/PDB panel was created by incorrect version of Irena " 1, "POV/PDB panel may need to be restarted to work properly. Restart now?"
145                        if(V_flag==1)
146                                DoWindow/K POVPDBPanel
147                                IR3P_ImportPOVPDB()
148                        else            //at least reinitialize the variables so we avoid major crashes...
149                                IR3P_InitializePOVPDB()
150                        endif
151                endif
152        endif
153end
154//******************************************************************************************************************************************************
155//******************************************************************************************************************************************************
156Function IR3A_MainCheckVersion()       
157        //this needs to get more of these lines for each tool/panel...
158        DoWindow FractalAggregatePanel
159        if(V_Flag)
160                if(!IR1_CheckPanelVersionNumber("FractalAggregatePanel", IR3AMassFrAggVersionNumber))
161                        DoAlert /T="The Mass Fractal Aggregate panel was created by incorrect version of Irena " 1, "Mass Fractal Aggregate may need to be restarted to work properly. Restart now?"
162                        if(V_flag==1)
163                                DoWindow/K FractalAggregatePanel
164                                IR3A_MassFractalAggregate()
165                        else            //at least reinitialize the variables so we avoid major crashes...
166                                IR3A_InitializeMassFractAgg()
167                        endif
168                endif
169        endif
170end
171//******************************************************************************************************************************************************
172//******************************************************************************************************************************************************
173//                      3D aggregate code, modified from Alex 2018-12-26
174//******************************************************************************************************************************************************
175//******************************************************************************************************************************************************
176
177Function IR3A_FractalAggregatePanel()
178        PauseUpdate             // building window...
179        NewPanel /K=1 /W=(5,20,395,680) as "Fractal Aggregate Model"
180        DoWindow/C FractalAggregatePanel
181        DefaultGUIControls /W=FractalAggregatePanel ///Mac os9
182        TitleBox MainTitle title="\Zr200Mass Fractal Aggregate model",pos={20,0},frame=0,fstyle=3, fixedSize=1,font= "Times New Roman", size={350,24},anchor=MC,fColor=(0,0,52224)
183        Button GetHelp,pos={305,50},size={80,15},fColor=(65535,32768,32768), proc=IR3A_PanelButtonProc,title="Get Help", help={"Open www manual page for this tool"}    //<<< fix button to help!!!
184        //COPY FROM IR2U_UnifiedEvaPanelFnct()
185        Checkbox  CurrentResults, pos={10,52}, size={50,15}, variable =root:packages:AggregateModeling:CurrentResults
186        Checkbox  CurrentResults, title="Current Unified Fit",mode=1,proc=IR3A_CheckProc
187        Checkbox  CurrentResults, help={"Select of you want to analyze current results in Unified Fit tool"}
188        Checkbox StoredResults, pos={150,52}, size={50,15}, variable =root:packages:AggregateModeling:StoredResults
189        Checkbox  StoredResults, title="Stored Unified Fit results",mode=1, proc=IR3A_CheckProc
190        Checkbox  StoredResults, help={"Select of you want to analyze Stored Unified fit data"}
191       
192        string UserDataTypes=""
193        string UserNameString=""
194        string XUserLookup=""
195        string EUserLookup=""
196        IR2C_AddDataControls("AggregateModeling","FractalAggregatePanel","","UnifiedFitIntensity",UserDataTypes,UserNameString,XUserLookup,EUserLookup, 1,1)
197        NVAR UseResults = root:Packages:AggregateModeling:UseResults
198        UseResults=1
199        STRUCT WMCheckboxAction CB_Struct
200        CB_Struct.ctrlName="UseResults"
201        CB_Struct.checked=1
202        CB_Struct.win="FractalAggregatePanel"
203        CB_Struct.eventcode=2
204       
205        IR2C_InputPanelCheckboxProc(CB_Struct)
206        Checkbox  UseResults, disable=1
207        KillControl UseModelData
208        KillCOntrol UseQRSData
209        PopupMenu ErrorDataName, disable=1
210        PopupMenu QvecDataName, disable=1
211        PopupMenu SelectDataFolder,pos={10,75}
212        PopupMenu SelectDataFolder proc=IR3A_PopMenuProc
213        PopupMenu IntensityDataName,pos={10,100}
214        PopupMenu IntensityDataName proc=IR3A_PopMenuProc       
215        Setvariable FolderMatchStr, pos={288,75}
216        KillControl Qmin
217        KillControl Qmax
218        KillControl QNumPoints
219
220        PopupMenu AvailableLevels,pos={288,100},size={109,20},proc=IR3A_PopMenuProc,title="Level:"
221        PopupMenu AvailableLevels,help={"Select level to use for data analysis"}
222        PopupMenu AvailableLevels,mode=1,popvalue="---", value=#"root:Packages:AggregateModeling:AvailableLevels"
223        SetVariable BrFract_ErrorMessage, title=" ",value=root:Packages:AggregateModeling:BrFract_ErrorMessage, noedit=1
224        SetVariable BrFract_ErrorMessage, pos={5,122}, size={375,20}, frame=0, help={"Error message, if any"}   
225
226        SetVariable BrFract_z, pos={10,140}, size={200,20}, title="Deg. of aggregation z     = ", help={"Degree of aggregation, 1+G2/G1"}, format="%.4g"
227        SetVariable BrFract_z, variable=root:Packages:AggregateModeling:BrFract_z, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0, bodyWidth=50
228        SetVariable BrFract_Rg1, pos={240,140}, size={150,20}, title="Rg primary part [A] =  ", help={"Rg of primary particle"}, format="%.4g",bodyWidth=50
229        SetVariable BrFract_Rg1, variable=root:Packages:AggregateModeling:BrFract_Rg1, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0
230
231        SetVariable BrFract_dmin, pos={10,160}, size={200,20}, title="Min. dimension  dmin = ", help={"Minimum dimension of the aggregate"}, format="%.4g"
232        SetVariable BrFract_dmin, variable=root:Packages:AggregateModeling:BrFract_dmin, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0, bodyWidth=50
233        SetVariable BrFract_df, pos={240,160}, size={150,20}, title="Fractal dimension df = ", help={"Mass fractal dimension of the aggregate"}, format="%.4g"
234        SetVariable BrFract_df, variable=root:Packages:AggregateModeling:BrFract_df, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0, bodyWidth=50
235//      SetVariable BrFract_fM, pos={25,190}, size={120,20}, title="fM =   ", help={"Parameter as defined in the references"}, format="%.4g"
236//      SetVariable BrFract_fM, variable=root:Packages:AggregateModeling:BrFract_fM, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0
237
238        SetVariable BrFract_c, pos={10,180}, size={200,20}, title="Connect. dimension c     = ", help={"Connectivity dimension of the aggregate"}, format="%.4g"
239        SetVariable BrFract_c, variable=root:Packages:AggregateModeling:BrFract_c, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0, bodyWidth=50
240        SetVariable BrFract_Rg2, pos={240,180}, size={150,20}, title="Rg aggregate [A]  =  ", help={"Rg of aggregate"}, format="%.4g",bodyWidth=50
241        SetVariable BrFract_Rg2, variable=root:Packages:AggregateModeling:BrFract_Rg2, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0
242       
243        //      R - aggregate size
244        //      df - Mass fractal dimension of the aggregate
245        //      p - short circuit path length
246        //      s - connective path length
247        //      dmin - minimum dimension of the aggregate
248        //      c - connecGvity dimension of the aggregate
249        //      s - connecGve path length of the aggregate
250       
251        //and this is new code
252        TitleBox Info1 title="\Zr160Unified results",pos={40,26},frame=0,fstyle=1, fixedSize=1,size={280,18},fColor=(0,0,52224)
253        TitleBox FakeLine1 title=" ",fixedSize=1,size={350,3},pos={16,215},frame=0,fColor=(0,0,52224), labelBack=(0,0,52224)
254        TitleBox Info2 title="\Zr160Model input parameters",pos={40,230},frame=0,fstyle=1, fixedSize=1,size={280,18},fColor=(0,0,52224)
255        SetVariable AggregateModeling,pos={10,260},size={190,16},noproc,title="Deg. of Agg. \"z\" (250) ", help={"Size of aggregate, degree of aggregation"}
256        SetVariable AggregateModeling,limits={10,Inf,0},variable= root:Packages:AggregateModeling:DegreeOfAggregation, bodyWidth=50
257        SetVariable StickingProbability,pos={10,285},size={190,16},noproc,title="Sticking prob. (10-100) ", help={"Sticking probablility, 100 for DLA, less for RLA"}
258        SetVariable StickingProbability,limits={10,100,0},variable= root:Packages:AggregateModeling:StickingProbability, bodyWidth=50
259        SetVariable NumberOfTestPaths,pos={10,310},size={190,16},noproc,title="Max paths/end (1k-10k) ", help={"Max measured paths per end point for parameter evaluation. Larger = possibly longer times but is more precise"}
260        SetVariable NumberOfTestPaths,limits={1000,100000,0},variable= root:Packages:AggregateModeling:NumberOfTestPaths, bodyWidth=50
261
262        SetVariable RgPrimary,pos={210,260},size={170,16},noproc,title="Primary Rg[A] (10) ", help={"Size of primary particle from which Aggregate is created"}
263        SetVariable RgPrimary,limits={10,Inf,0},variable= root:Packages:AggregateModeling:RgPrimary, bodyWidth=50
264
265        PopupMenu AllowedNearDistance,pos={220,290},size={150,20},proc=IR3A_PopMenuProc,title="Sticking method:"
266        PopupMenu AllowedNearDistance,help={"Which neighbors are allowed to stick"}
267        NVAR AllowedNearDistance=root:Packages:AggregateModeling:AllowedNearDistance
268        PopupMenu AllowedNearDistance,mode=1,popvalue=num2str(AllowedNearDistance), value="1;2;3;"
269
270        PopupMenu MParticleAttraction,pos={220,310},size={150,20},proc=IR3A_PopMenuProc,title="Multi Part. attr:"
271        PopupMenu MParticleAttraction,help={"If there are more particles neaby, is chance of attaching? "}
272        SVAR MultiParticleAttraction=root:Packages:AggregateModeling:MultiParticleAttraction
273        PopupMenu MParticleAttraction,mode=1,popvalue=MultiParticleAttraction, value="Neutral;Attractive;Repulsive;Not allowed;"
274       
275        Button Grow1AggAll,pos={5,338},size={150,20}, proc=IR3A_PanelButtonProc,title="Grow 1 Agg, graph", help={"Perform all steps and generate 3D graph"}
276        Button GrowNAggAll,pos={165,338},size={150,20}, proc=IR3A_PanelButtonProc,title="Grow N Agg.", help={"Generate N aggregates randomly"}
277
278        PopupMenu NUmberOfTestAggregates,pos={320,338},size={150,20},proc=IR3A_PopMenuProc,title="N ="
279        PopupMenu NUmberOfTestAggregates,help={"How many test aggregates to grow"}
280        NVAR NUmberOfTestAggregates=root:Packages:AggregateModeling:NUmberOfTestAggregates
281        PopupMenu NUmberOfTestAggregates,mode=1,popvalue=num2str(NUmberOfTestAggregates), value="5;10;20;30;50;"
282
283
284
285        //repeat the target values
286        SetVariable Target_dmin, pos={10,360}, size={100,20}, title="Target dmin = ", help={"Minimum dimension of the aggregate"}, format="%.2g",fColor=(1,16019,65535),valueColor=(16385,16388,65535)
287        SetVariable Target_dmin, variable=root:Packages:AggregateModeling:BrFract_dmin, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0, bodyWidth=50
288        SetVariable Target_c, pos={150,360}, size={100,20}, title="Target c =  ", help={"Connectivity dimension of the aggregate"}, format="%.2g", fColor=(1,16019,65535),valueColor=(16385,16388,65535)
289        SetVariable Target_c, variable=root:Packages:AggregateModeling:BrFract_c, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0, bodyWidth=50
290        SetVariable Target_df, pos={280,360}, size={100,20}, title="Target df = ", help={"Mass fractal dimension of the aggregate"}, format="%.2g", fColor=(1,16019,65535),valueColor=(16385,16388,65535)
291        SetVariable Target_df, variable=root:Packages:AggregateModeling:BrFract_df, noedit=1,limits={-inf,inf,0},  disable=0,noedit=1,frame=0, bodyWidth=50
292
293        // here are resulting values
294        SetVariable dminValue,pos={40,380},size={100,20},noproc,title="dmin = ", help={"Minimum dimension of the aggregate"},limits={10,100,0}
295        SetVariable dminValue,variable= root:Packages:AggregateModeling:dminValue, bodyWidth=80, disable=0,noedit=1,frame=0, format="%.4f"
296        SetVariable cValue,pos={150,380},size={100,20},noproc,title="c = ", help={"Connectivity dimension of the aggregate"},limits={10,100,0}
297        SetVariable cValue,variable= root:Packages:AggregateModeling:cValue, bodyWidth=80, disable=0,noedit=1,frame=0, format="%.4f"
298        SetVariable dfValue,pos={280,380},size={100,20},noproc,title="df = ", help={"Mass fractal dimension of the aggregate"},limits={10,100,0}, format="%.4f"
299        SetVariable dfValue,variable= root:Packages:AggregateModeling:dfValue, bodyWidth=80, disable=0,noedit=1,frame=0
300
301
302        SetVariable RValue,pos={40,400},size={100,16},noproc,title="R [A] = ", help={"R of the aggregate, Rg/dp"},limits={10,100,0}, format="%6.2f"
303        SetVariable RValue,variable= root:Packages:AggregateModeling:RxRgPrimaryValue, bodyWidth=80, disable=0,noedit=1,frame=0
304        SetVariable pValue,pos={150,400},size={100,16},noproc,title="p = ", help={"Short circuit path length"},limits={10,100,0}, format="%6.0f"
305        SetVariable pValue,variable= root:Packages:AggregateModeling:pValue, bodyWidth=80, disable=0,noedit=1,frame=0
306
307        SetVariable sValue,pos={280,400},size={100,16},noproc,title="s = ", help={"Connective path length of the aggregate"},limits={10,100,0}
308        SetVariable sValue,variable= root:Packages:AggregateModeling:sValue, bodyWidth=80, disable=0,noedit=1,frame=0, format="%6.0f"
309
310        SetVariable TrueStickingProbability pos={20,430},size={200,16},noproc,title="True Sticking Prob. [%] = ", help={"True Sticking Probability during aggregate growth"}
311        SetVariable TrueStickingProbability variable=root:Packages:AggregateModeling:TrueStickingProbability, noedit=1,frame=0,limits={10,100,0}, format="%2.2d"
312
313        Button Display3DMFASummary,pos={220,425},size={150,20}, proc=IR3A_PanelButtonProc,title="Summary Table", help={"Display summary for 3D Mass Fractal Aggregate"}
314        Button SaveAggregateData,pos={220,445},size={150,20}, proc=IR3A_PanelButtonProc,title="Store Current Aggregate", help={"Copy this aggregate with parameters in a folder"}
315        TitleBox FakeLine2 title=" ",fixedSize=1,size={330,3},pos={16,469},frame=0,fColor=(0,0,52224), labelBack=(0,0,52224)
316
317        TitleBox ListBoxTitle title="\Zr130Saved 3D Mass Fract aggregates",size={250,15},pos={5,480},frame=0,fColor=(0,0,52224)
318        wave/T Stored3DAggregates= root:Packages:AggregateModeling:Stored3DAggregates
319        wave Stored3DAggSelections= root:Packages:AggregateModeling:Stored3DAggSelections
320        ListBox StoredAggregates pos={5,500}, size={250,145}
321        ListBox StoredAggregates listWave=Stored3DAggregates,mode=1,selRow=-1,selWave=Stored3DAggSelections
322       
323        Button Display1DData,pos={262,480},size={120,20}, proc=IR3A_PanelButtonProc,title="Display 1D Graph", help={"Display 1 D graph with Intensity vs Q for these data"}
324        Button Display3DMassFracGizmo,pos={262,505},size={120,20}, proc=IR3A_PanelButtonProc,title="Display 3D Graph", help={"Display Gizmo with 3D Mass Fractal Aggregate"}
325        Button Calculate1DIntensity,pos={262,530},size={120,20}, proc=IR3A_PanelButtonProc,title="Calculate 1D Int.", help={"Calculate using UF 1D intensity and append to graph"}
326        Button Model1DIntensity,pos={262,555},size={120,20}, proc=IR3A_PanelButtonProc,title="Monte Carlo 1D Int.", help={"Calculate using Monte Carlo 1D intensity and append to graph"}
327        Button CompareStoredResults,pos={262,580},size={120,20}, proc=IR3A_PanelButtonProc,title="Compare Stored", help={"Present statistcs on stored results"}
328        Button DeleteStoredResults,pos={262,635},size={100,15}, proc=IR3A_PanelButtonProc,title="Delete all Stored", help={"Delete all stored results"}
329        //the above button works, but the results seem to take forever and do not look too good. 
330        IR3A_SetControlsInPanel()
331        IR3A_Create3DAggListForListbox()
332end
333
334//******************************************************************************************************************************************************
335//******************************************************************************************************************************************************
336//******************************************************************************************************************************************************
337//******************************************************************************************************************************************************
338//******************************************************************************************************************************************************
339//******************************************************************************************************************************************************
340//******************************************************************************************************************************************************
341static Function IR3A_Create3DAggListForListbox()
342        DFref oldDf= GetDataFolderDFR()
343
344
345        DOWindow FractalAggregatePanel
346        if(!V_Flag)
347                return 0
348        endif
349        wave/Z Stored3DAggSelections= root:Packages:AggregateModeling:Stored3DAggSelections
350        if(!WaveExists(Stored3DAggSelections))
351                abort
352        endif
353        wave/T Stored3DAggregates= root:Packages:AggregateModeling:Stored3DAggregates
354        wave/T Stored3DAggregatesPaths= root:Packages:AggregateModeling:Stored3DAggregatesPaths
355        variable NumOfFolders
356        string CurrentList, tempStr
357        if(DataFolderExists("root:MassFractalAggregates"))
358                setDataFolder root:MassFractalAggregates
359                CurrentList=stringByKey("FOLDERS",DataFolderDir(1),":")
360                NumOfFolders = ItemsInList(CurrentList,",")+1
361        else
362                CurrentList=""
363                NumOfFolders = 1
364        endif
365        redimension/N=(NumOfFolders) Stored3DAggregates, Stored3DAggSelections, Stored3DAggregatesPaths
366        Stored3DAggSelections = 0
367        Stored3DAggregates[0] = "Current model"
368        Stored3DAggregatesPaths[0] = "Current model"
369        variable i
370        if(NumOfFolders>1)
371                For(i=0;i<NumOfFolders-1;i+=1)
372                        tempStr = StringFromList(i, CurrentList, ",")
373                        Stored3DAggregatesPaths [i+1] ="root:MassFractalAggregates:"+tempStr
374                        Stored3DAggregates      [i+1] = num2str(i)+" : "+IR3A_BuildUser3DAggNames("root:MassFractalAggregates:"+tempStr)
375                endfor 
376        endif
377        DoUpdate  /W=FractalAggregatePanel
378        setDataFOlder OldDf
379end
380//******************************************************************************************************************************************************
381//******************************************************************************************************************************************************
382static Function/T IR3A_BuildUser3DAggNames(PathToFolder)
383        string PathToFOlder
384        DFref oldDf= GetDataFolderDFR()
385
386        string UserFriendlyString=""
387        SetDataFolder PathToFolder
388        NVAR DOA=DegreeOfAggregation
389        NVAR Stick=StickingProbability
390        NVAR SMeth=StickingMethod
391        NVAR dMin=DminValue
392        NVAR df = dfValue
393        NVAR cval= cValue
394        UserFriendlyString="z="+num2str(DOA)+",dmin="+num2str(IN2G_roundDecimalPlaces(dmin,2))+",c="+num2str(IN2G_roundDecimalPlaces(cval,2))+",df="+num2str(IN2G_roundDecimalPlaces(df,2))
395        setDataFOlder OldDf
396        return UserFriendlyString
397end
398
399//******************************************************************************************************************************************************
400//******************************************************************************************************************************************************
401static Function IR3A_Display1DIntensity()
402
403        //DoAlert /T="Need calcuation here" 0, "IR3A_Display1DIntensity is not finished yet"
404        NVAR useCurrentResults = root:packages:AggregateModeling:CurrentResults
405        NVAR useStoredResults = root:packages:AggregateModeling:StoredResults
406        if(useCurrentResults+useStoredResults!=1)
407                useStoredResults = 0
408                useCurrentResults = 1
409        endif
410        if(useCurrentResults)
411                Wave/Z IntWaveOriginal = root:Packages:Irena_UnifFit:OriginalIntensity
412                Wave/Z QwaveOriginal = root:Packages:Irena_UnifFit:OriginalQvector
413                Wave/Z ErrorOriginal = root:Packages:Irena_UnifFit:OriginalError
414                if(!WaveExists(IntWaveOriginal))
415                        abort
416                endif
417        else //use stored results, in this case the strings shoudl be set...
418                SVAR DataFolderName = root:Packages:AggregateModeling:DataFolderName
419                SVAR IntensityWaveName = root:Packages:AggregateModeling:IntensityWaveName
420                SVAR QWavename = root:Packages:AggregateModeling:QWavename
421                SVAR ErrorWaveName = root:Packages:AggregateModeling:ErrorWaveName
422                Wave/Z IntWaveOriginal = $(DataFolderName+IntensityWaveName)
423                Wave/Z QwaveOriginal = $(DataFolderName+QWavename)
424                Wave/Z ErrorOriginal = $(DataFolderName+ErrorWaveName)
425                if(!WaveExists(IntWaveOriginal))
426                        abort
427                endif
428        endif
429        string IntWvName=nameofWave(IntWaveOriginal)
430        DoWIndow MassFractalAggDataPlot
431        if(V_Flag)
432                DoWIndow/F MassFractalAggDataPlot
433        else
434                Display /W=(282.75,37.25,900,400)/K=1  IntWaveOriginal vs QwaveOriginal as "Mass Fractal Aggregate 1D Data Plot"
435                DoWindow/C MassFractalAggDataPlot
436                ModifyGraph mode($(nameofwave(IntWaveOriginal)))=3
437                ModifyGraph msize($(nameofwave(IntWaveOriginal)))=0
438                ModifyGraph log=1
439                ModifyGraph mirror=1
440                ShowInfo
441                String LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Intensity"
442                Label left LabelStr
443                LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Q [A\\S-1\\M\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"]"
444                Label bottom LabelStr
445                string LegendStr="\\F"+IN2G_LkUpDfltStr("FontType")+"\\Z"+IN2G_LkUpDfltVar("LegendSize")+"\\s("+IntWvName+") Data intensity"
446                Legend/W=MassFractalAggDataPlot/N=text0/J/F=0/A=MC/X=32.03/Y=38.79 LegendStr
447                TextBox/C/N=DateTimeTag/F=0/A=RB/E=2/X=2.00/Y=1.00 "\\Z07"+date()+", "+time()   
448                TextBox/C/N=SampleNameTag/F=0/A=LB/E=2/X=2.00/Y=1.00 "\\Z07"+GetWavesDataFolder(IntWaveOriginal,1)     
449        endif
450        AutoPositionWindow/R=FractalAggregatePanel MassFractalAggDataPlot
451       
452end
453
454
455//******************************************************************************************************************************************************
456//******************************************************************************************************************************************************
457static Function IR3A_Append1DInMassFracAgg(IntensityWV,QvectorWV)
458        wave IntensityWV,QvectorWV
459
460        DoWIndow MassFractalAggDataPlot
461        if(V_Flag)
462                DoWIndow/F MassFractalAggDataPlot
463        else
464                IR3A_Display1DIntensity()
465        endif
466
467        NVAR useCurrentResults = root:packages:AggregateModeling:CurrentResults
468        NVAR useStoredResults = root:packages:AggregateModeling:StoredResults
469        if(useCurrentResults+useStoredResults!=1)
470                useStoredResults = 0
471                useCurrentResults = 1
472        endif
473        //fix q rave of overlap
474        //proper scaling of model... this is not so easy here, but may be needed in teh future.
475        //variable QRg2 = 0.5*pi/Level2Rg
476        //variable QRg1 = 1.5*pi/Level1Rg
477        //print "Scaling Model to data using integral intensities in Q range from "+num2str(QRg2)+"  to "+num2str(QRg1)
478        //variable InvarModel=areaXY(A3DAgg1DQwave, Model3DAggIntensity, QRg2, QRg1)
479        variable InvarModel=areaXY(QvectorWV, IntensityWV )
480        variable InvarData
481        //Model3DAggIntensity
482        //need to also work, when no data are present. In this case the root:Packages:AggregateModeling:DataFolderName is set to non-sensical value...
483        SVAR DataFolderName = root:Packages:AggregateModeling:DataFolderName
484        if(strlen(DataFolderName)<4)
485                //this is non sensical case. Use Model3DAggIntensity to get proper normalization
486                Wave/Z IntWaveOriginal = root:Aggregate1DModel:Model3DAggIntensity
487                Wave/Z QwaveOriginal = root:Aggregate1DModel:A3DAgg1DQwave
488                if(!WaveExists(IntWaveOriginal))
489                        abort
490                endif
491                InvarData=areaXY(QwaveOriginal, IntWaveOriginal, QvectorWV[0], QvectorWV[numpnts(QvectorWV)-1])
492        else
493                if(useCurrentResults)
494                        Wave/Z IntWaveOriginal = root:Packages:Irena_UnifFit:OriginalIntensity
495                        Wave/Z QwaveOriginal = root:Packages:Irena_UnifFit:OriginalQvector
496                        Wave/Z ErrorOriginal = root:Packages:Irena_UnifFit:OriginalError
497                        if(!WaveExists(IntWaveOriginal))
498                                abort
499                        endif
500                else //use stored results, in this case the strings should be set...
501                        SVAR DataFolderName = root:Packages:AggregateModeling:DataFolderName
502                        SVAR IntensityWaveName = root:Packages:AggregateModeling:IntensityWaveName
503                        SVAR QWavename = root:Packages:AggregateModeling:QWavename
504                        SVAR ErrorWaveName = root:Packages:AggregateModeling:ErrorWaveName
505                        Wave/Z IntWaveOriginal = $(DataFolderName+IntensityWaveName)
506                        Wave/Z QwaveOriginal = $(DataFolderName+QWavename)
507                        Wave/Z ErrorOriginal = $(DataFolderName+ErrorWaveName)
508                        if(!WaveExists(IntWaveOriginal))
509                                abort
510                        endif
511                endif
512                InvarData=areaXY(QwaveOriginal, IntWaveOriginal, QvectorWV[0], QvectorWV[numpnts(QvectorWV)-1])
513        endif
514        IntensityWV*=InvarData/InvarModel
515        CheckDisplayed /W=MassFractalAggDataPlot IntensityWV
516        if(V_flag==0)
517                        AppendToGraph/W=MassFractalAggDataPlot IntensityWV vs QvectorWV
518        endif
519                 
520        ModifyGraph/W=MassFractalAggDataPlot mode($(nameofwave(IntensityWV)))=0,rgb($(nameofwave(IntensityWV)))=(0,0,0)         
521        ModifyGraph/W=MassFractalAggDataPlot log=1
522        ModifyGraph/W=MassFractalAggDataPlot mirror=1
523        ModifyGraph/W=MassFractalAggDataPlot lsize(PDFIntensityWv)=3,rgb(PDFIntensityWv)=(1,16019,65535)       
524        ShowInfo
525        String LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Intensity"
526        Label/W=MassFractalAggDataPlot left LabelStr
527        LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Q [A\\S-1\\M\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"]"
528        Label/W=MassFractalAggDataPlot bottom LabelStr
529//              string LegendStr="\\F"+IN2G_LkUpDfltStr("FontType")+"\\Z"+IN2G_LkUpDfltVar("LegendSize")+"\\s(PDFIntensityWv) Mass Fractal Model intensity"
530//              Legend/W=MassFractalAggDataPlot/N=text0/J/F=0/A=MC/X=32.03/Y=38.79 LegendStr
531
532//      string IntWvName=nameofWave(IntWaveOriginal)
533//      DoWIndow MassFractalAggDataPlot
534//      if(V_Flag)
535//              DoWIndow/F MassFractalAggDataPlot
536//      else
537//              Display /W=(282.75,37.25,900,400)/K=1  IntWaveOriginal vs QwaveOriginal as "Mass Fractal Aggregate 1D Data Plot"
538//              DoWindow/C MassFractalAggDataPlot
539//              ModifyGraph mode($(nameofwave(IntWaveOriginal)))=3
540//              ModifyGraph msize($(nameofwave(IntWaveOriginal)))=0
541//              ModifyGraph log=1
542//              ModifyGraph mirror=1
543//              ShowInfo
544//              String LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Intensity"
545//              Label left LabelStr
546//              LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Q [A\\S-1\\M\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"]"
547//              Label bottom LabelStr
548//              string LegendStr="\\F"+IN2G_LkUpDfltStr("FontType")+"\\Z"+IN2G_LkUpDfltVar("LegendSize")+"\\s("+IntWvName+") Data intensity"
549//              Legend/W=MassFractalAggDataPlot/N=text0/J/F=0/A=MC/X=32.03/Y=38.79 LegendStr
550//              TextBox/C/N=DateTimeTag/F=0/A=RB/E=2/X=2.00/Y=1.00 "\\Z07"+date()+", "+time()   
551//              TextBox/C/N=SampleNameTag/F=0/A=LB/E=2/X=2.00/Y=1.00 "\\Z07"+GetWavesDataFolder(IntWaveOriginal,1)     
552//      endif
553        AutoPositionWindow/R=FractalAggregatePanel MassFractalAggDataPlot
554end
555
556
557//******************************************************************************************************************************************************
558//******************************************************************************************************************************************************
559static Function IR3A_Calculate1DIntensity()
560       
561        DFref oldDf= GetDataFolderDFR()
562
563        SetDataFolder root:Packages:AggregateModeling
564        //decide which data - if table has selected data, use that, else current data in the tool
565        string selection = IR3A_FindSelectedAggData()
566        variable isInWorkFolder = 0
567        if(strlen(selection)<1 || stringmatch(selection, "root:Packages:*"))
568                isInWorkFolder =1
569        endif
570        Wave/Z MassFractalAggregate=$(selection)
571        if(!WaveExists(MassFractalAggregate))
572                        Print "MassFractalAggregate 3D wave does not exist, cannot do anything"
573                        return 0                        //end on user.
574        endif
575        setDataFOlder $(GetWavesDataFolder(MassFractalAggregate, 1 ))
576        string OldNote=note(MassFractalAggregate)
577
578        //figure out where we are working...
579        variable measuredDataExists = 0
580        string CurrentFolder
581        NVAR useCurrentResults = root:packages:AggregateModeling:CurrentResults
582        NVAR useStoredResults = root:packages:AggregateModeling:StoredResults
583        if(useCurrentResults+useStoredResults!=1)
584                useStoredResults = 0
585                useCurrentResults = 1
586        endif
587        if(useCurrentResults)
588                Wave/Z IntWaveOriginal = root:Packages:Irena_UnifFit:OriginalIntensity
589                Wave/Z QwaveOriginal = root:Packages:Irena_UnifFit:OriginalQvector
590                Wave/Z ErrorOriginal = root:Packages:Irena_UnifFit:OriginalError
591                if(WaveExists(IntWaveOriginal))
592                        measuredDataExists = 1
593                endif
594                CurrentFolder = "root:Packages:Irena_UnifFit"
595        else //use stored results, in this case the strings shoudl be set...
596                SVAR DataFolderName = root:Packages:AggregateModeling:DataFolderName
597                SVAR IntensityWaveName = root:Packages:AggregateModeling:IntensityWaveName
598                SVAR QWavename = root:Packages:AggregateModeling:QWavename
599                SVAR ErrorWaveName = root:Packages:AggregateModeling:ErrorWaveName
600                Wave/Z IntWaveOriginal = $(DataFolderName+IntensityWaveName)
601                Wave/Z QwaveOriginal = $(DataFolderName+QWavename)
602                Wave/Z ErrorOriginal = $(DataFolderName+ErrorWaveName)
603                if(WaveExists(IntWaveOriginal))
604                        measuredDataExists = 1
605                endif
606                CurrentFolder = DataFolderName
607        endif
608        if(!measuredDataExists)
609                NewDataFOlder/O root:Aggregate1DModel
610                CurrentFolder = "root:Aggregate1DModel"
611        endif
612       
613        //OK, now lets go where we have the data..
614        SetDataFolder $(CurrentFolder)
615        if(measuredDataExists)
616                Duplicate/O IntWaveOriginal, Model3DAggIntensity
617                Duplicate/O QwaveOriginal, A3DAgg1DQwave
618        else
619                Make/O/N=200 IntWaveOriginal, Model3DAggIntensity, QwaveOriginal, A3DAgg1DQwave
620                A3DAgg1DQwave = 10^(log(0.001)+p*((log(1)-log(0.001))/(199)))                           //sets k scaling on log-scale...
621                QwaveOriginal = 10^(log(0.001)+p*((log(1)-log(0.001))/(199)))                           //sets k scaling on log-scale...
622                IntWaveOriginal = 1
623        endif
624        //DoAlert /T="This is not working right yet" 0, "IR3A_Calculate1DIntensity() is not finished yet, using Guinier-Porod does nto work right... "
625        //SetDataFolder OldDf
626        //abort
627        //now, use info from Alex McGlasson who back calculated Unified parameters from Andrew Mulderig J. Aerosol Sci. 109 (2017), 28-37 manuscript
628        //NoteText="Mass Fractal Aggregate created="+date()+", "+time()+";z="+num2str(DegreeOfAggregation)+";StickingProbability="+num2str(StickingProbability)+
629        //";R="+num2str(RValue)+";Rprimary="+num2str(RgPrimary)+";p="+num2str(pValue)
630        //NoteText+=";RxRgPrimaryValue="+num2str(RValue*RgPrimary)+";s="+num2str(sValue)+
631        //";df="+num2str(IN2G_roundSignificant(dfValue,3))+";dmin="+num2str(IN2G_roundSignificant(dminValue,3))+
632        //";c="+num2str(IN2G_roundSignificant(cValue,3))+";True Sticking Probability="+num2str(100*DegreeOfAggregation/AttemptValue)+";"
633        variable Level1Rg = str2num(StringByKey("Rprimary", OldNote, "=", ";"))                                 //recorded is small dimension Rg in A
634        variable level1Radius =         sqrt(5/3)*Level1Rg                                                                                                                      //Rg^2 = 3 * R^2 / 5
635        variable Level2P = str2num(StringByKey("df", OldNote, "=", ";"))                                                        //this is mass fractal dimension
636        variable zval = str2num(StringByKey("z", OldNote, "=", ";"))
637        variable dfval = str2num(StringByKey("df", OldNote, "=", ";"))
638        variable dminval = str2num(StringByKey("dmin", OldNote, "=", ";"))
639        variable Level2G = zval - 1                                                                                                                                                                             //this is assuming level1G = 1 (which we use as default here)
640        //variable Level2RgAggregate = str2num(StringByKey("RxRgPrimaryValue", OldNote, "=", ";"))      //recorded is large dimension Rg in A, R = sqrt(5/3) * Rg
641        variable R2val = str2num(StringByKey("R", OldNote, "=", ";"))
642        variable cval = str2num(StringByKey("c", OldNote, "=", ";"))
643        //variable Level2Rg = (R2val^2 / 4 ) *  zval^(2/dfval)                                                                                                  //per Alex's original formula. way too large
644        variable Level2Rg =     Level1Rg * zval^((1/cval - 1)/(dminval - dfval))                                                                //per Alex's second formula.
645        //print "Aggregate Radius of Gyration is [A] : "+num2str(Level2Rg)
646        //this above is suppose to be R^2/4 * z^(2/df) 
647        //this seems to be unitless, do we need to do following?
648        //Level2Rg = Level2Rg * Level1Rg                                                                                                                                                //all calculates seem relative to size of the primary particle, so do we need to scale this by primary size somehow?
649        //this is Unified fit calculation...
650        //      if (MassFractal)
651        //              B=(G*P/Rg^P)*exp(gammln(P/2))
652        //      endif
653        //this below is df * Gamma(df/2) * G2 / (c * Rg2^df) 
654        //variable Level2B = dfval * gamma(dfval/2) * Level2G /(cval * Level2Rg^dfval) 
655        //use alternativer from mass fractal...
656        //variable Level2B = (Level2G*Level2P/Level2Rg^Level2P)*exp(gammln(Level2P/2))
657        variable Level2B = (Level2G*Level2P/Level2Rg^Level2P)*gamma(Level2P/2)
658       
659
660        variable Level2RgCO = Level1Rg                                                                                                                                          //this is logical, needs to terminate mass fractal curve here...
661        variable Level1G = 1                                                                                                                                                                                            //see above, this is model basic volume fraction term. Needs to be scaled to measured dat athrough invariant.
662        variable Level1B = 4*pi/(Level1Rg^4)                                                                                                                                    //does this need to be converted to 1/cm ???
663        // Alex has Rg1 calculation here: Rg1 = Rg2/(z^((1/c - 1)/(dmin-df)))                                                           //it seems round calculations, I am not sure I understand this whole thing philosophically.
664        variable Level1P = 4                                                                                                                                                                                            //see above, this is model basic volume fraction term. Needs to be scaled to measured dat athrough invariant.
665        print "Calculating 1D intensity, here are parameters for Unified fit model used to generate 1D curve"
666        print "Rg primary = "+num2str(Level1Rg)
667        print "G primary = "+num2str(Level1G)
668        print "B primary = "+num2str(Level1B)
669        print "P primary = "+num2str(Level1P)
670        print "Rg aggregate = "+num2str(Level2Rg)
671        print "G aggregate = "+num2str(Level2G)
672        print "B aggregate = "+num2str(Level2B)
673        print "P aggregate = "+num2str(Level2P)
674
675
676        Duplicate/Free Model3DAggIntensity, Model3DAggIntensityL1, Model3DAggIntensityL2
677        //level 1 - primaries...
678        IR3A_UnifiedFitIntOne(A3DAgg1DQwave,Model3DAggIntensityL1, Level1G, Level1Rg, Level1B, Level1P,  0)
679        //Level 2 - aggregate
680        IR3A_UnifiedFitIntOne(A3DAgg1DQwave,Model3DAggIntensityL2, Level2G, Level2Rg, Level2B, Level2P, Level2RgCO)
681        //summ together
682        Model3DAggIntensity = Model3DAggIntensityL1 + Model3DAggIntensityL2
683        //make graph if needed...
684        String LabelStr
685        DoWIndow MassFractalAggDataPlot
686        if(!V_Flag)
687                if(measuredDataExists)
688                        IR3A_Display1DIntensity()
689                else
690                        Display /W=(282.75,37.25,900,400)/K=1  IntWaveOriginal vs QwaveOriginal as "Mass Fractal Aggregate 1D Data Plot"
691                        DoWindow/C MassFractalAggDataPlot
692                        ModifyGraph mode($(nameofwave(IntWaveOriginal)))=3
693                        ModifyGraph msize($(nameofwave(IntWaveOriginal)))=0
694                        ModifyGraph log=1
695                        ModifyGraph mirror=1
696                        ShowInfo
697                        LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Intensity"
698                        Label left LabelStr
699                        LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Q [A\\S-1\\M\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"]"
700                        Label bottom LabelStr
701                        string LegendStr="\\F"+IN2G_LkUpDfltStr("FontType")+"\\Z"+IN2G_LkUpDfltVar("LegendSize")+"\\s(IntWaveOriginal) Data intensity"
702                        Legend/W=MassFractalAggDataPlot/N=text0/J/F=0/A=MC/X=32.03/Y=38.79 LegendStr
703                        TextBox/C/N=DateTimeTag/F=0/A=RB/E=2/X=2.00/Y=1.00 "\\Z07"+date()+", "+time()   
704                        TextBox/C/N=SampleNameTag/F=0/A=LB/E=2/X=2.00/Y=1.00 "\\Z07"+GetWavesDataFolder(IntWaveOriginal,1)     
705                endif
706        else
707                DoWIndow/F MassFractalAggDataPlot
708        endif
709        //proper scaling of model...
710        variable QRg2 = 0.5*pi/Level2Rg
711        variable QRg1 = 1.5*pi/Level1Rg
712        print "Scaling Model to data using integral intensities in Q range from "+num2str(QRg2)+"  to "+num2str(QRg1)
713        variable InvarModel=areaXY(A3DAgg1DQwave, Model3DAggIntensity, QRg2, QRg1)
714        variable InvarData=areaXY(QwaveOriginal, IntWaveOriginal, QRg2, QRg1)
715        Model3DAggIntensity*=InvarData/InvarModel
716        CheckDisplayed /W=MassFractalAggDataPlot Model3DAggIntensity
717        if(V_flag==0)
718                        AppendToGraph/W=MassFractalAggDataPlot Model3DAggIntensity vs A3DAgg1DQwave
719        endif
720                 
721        ModifyGraph/W=MassFractalAggDataPlot mode($(nameofwave(Model3DAggIntensity)))=0,rgb($(nameofwave(Model3DAggIntensity)))=(0,0,0)         
722        ModifyGraph/W=MassFractalAggDataPlot lsize(Model3DAggIntensity)=3
723        ModifyGraph/W=MassFractalAggDataPlot log=1
724        ModifyGraph/W=MassFractalAggDataPlot mirror=1
725        ShowInfo
726        LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Intensity"
727        Label/W=MassFractalAggDataPlot left LabelStr
728        LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Q [A\\S-1\\M\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"]"
729        Label/W=MassFractalAggDataPlot bottom LabelStr
730       
731        AutoPositionWindow/M=0 /R=FractalAggregatePanel MassFractalAggDataPlot
732        SetDataFolder OldDf
733       
734end
735///*********************************************************************************************************************
736///*********************************************************************************************************************
737
738
739Function IR3A_UnifiedFitIntOne(QvectorWave,IntensityWave, G, Rg, B, P, RGCO)
740        variable G, Rg, B, P, RGCO
741        Wave QvectorWave, IntensityWave
742       
743
744        Duplicate/Free QvectorWave, TempUnifiedIntensity, QstarVector
745       
746        variable Kval=1
747        QstarVector=QvectorWave/(erf(Kval * QvectorWave*Rg/sqrt(6)))^3
748//      if (MassFractal)
749//              B=(G*P/Rg^P)*exp(gammln(P/2))
750//      endif
751       
752        IntensityWave=G*exp(-QvectorWave^2*Rg^2/3)+(B/QstarVector^P) * exp(-RGCO^2 * QvectorWave^2/3)
753end
754
755
756///*********************************************************************************************************************
757///*********************************************************************************************************************
758//**************************************************************************************************************************************************************
759Function IR3A_CalcRgof3DAgg(MassFractalAggregate)
760        wave MassFractalAggregate
761        //note, this is 3 column listing of positions from Aggregate growth
762        //we will calculate Rg by using Rg^2 = sum(distances^2) / NumOfPoints
763        //this is in units of "simple cubic side" steps used in growth of the particles.
764        //this needs to be multiplied by center-to-center distance between particles in side of this simple cubic strucure
765        // if we know Rgp - primary particle Rg, then RadiusPrimary = sqrt(5/3)*Rgp and DiameterPrimary = 2*sqrt(5/3)*Rgp
766        //case example: Rgp=22, RadiusP= 28, and DiameterP=56 A
767        //multiply result of this by DiameterP to have real world Rg of aggregate...
768       
769        make/Free/N=(DimSize(MassFractalAggregate, 0)) DistanceWave2
770        DistanceWave2 = MassFractalAggregate[p][0]^2 + MassFractalAggregate[p][1]^2 + MassFractalAggregate[p][2]^2
771        variable tmpVal= sum(DistanceWave2)/numpnts(DistanceWave2)
772       
773        print "Rg in Primary Size units is : "+num2str(sqrt(tmpVal))
774        print "Multiply this by Diameter of primary particle =  2*sqrt(5/3)*Rgp "
775        return sqrt(tmpVal)
776end
777
778
779//**************************************************************************************************************************************************************
780
781//******************************************************************************************************************************************************
782//******************************************************************************************************************************************************
783
784static Function IR3A_Model1DIntensity()
785       
786        variable recalculate3D = 0
787        DFref oldDf= GetDataFolderDFR()
788
789        //follow this: IR3P_Calculate1DDataFile()
790        //decide which data - if table has selected data, use that, else current data in the tool
791        string selection = IR3A_FindSelectedAggData()
792        variable isInWorkFolder = 0
793        if(strlen(selection)<1 || stringmatch(selection, "root:Packages:*"))
794                isInWorkFolder =1
795        endif
796        Wave/Z MassFractalAggregate=$(selection)
797        if(!WaveExists(MassFractalAggregate))
798                        Print "MassFractalAggregate 3D wave does not exist, cannot do anything"
799                        return 0                        //end on user.
800        endif
801        setDataFOlder $(GetWavesDataFolder(MassFractalAggregate, 1 ))
802        string OldNote=note(MassFractalAggregate)
803        variable AggSize = str2num(StringByKey("RxRgPrimaryValue", OldNote, "=", ";"))          //recorded is large dimension Rg in A, R = sqrt(5/3) * Rg
804        variable Level1Rg = str2num(StringByKey("Rprimary", OldNote, "=", ";"))                                 //recorded is small dimension Rg in A
805        variable level1Radius =         sqrt(5/3)*Level1Rg                                                                                                                      //Rg^2 = 3 * R^2 / 5
806        variable PrimarySize = 2* sqrt(5/3)*Level1Rg                                                                                                                    //diameter of primary particle in Angstroms, = also side of the simple cubic lattice these particles are sitting on...
807        //NoteText="Mass Fractal Aggregate created="+date()+", "+time()+";z="+num2str(DegreeOfAggregation)+";StickingProbability="+num2str(StickingProbability)+";R="+num2str(RValue)+";Rprimary="+num2str(RgPrimary)+";p="+num2str(pValue)
808        //NoteText+=";RxRgPrimaryValue="+num2str(RValue*RgPrimary)+";s="+num2str(sValue)+";df="+num2str(IN2G_roundSignificant(dfValue,3))+";dmin="+num2str(IN2G_roundSignificant(dminValue,3))+";c="+num2str(IN2G_roundSignificant(cValue,3))+";True Sticking Probability="+num2str(100*DegreeOfAggregation/AttemptValue)+";"
809       
810        //pick the parameters...
811        variable voxelSize = PrimarySize/10
812        variable IsoValue = 0.1
813        variable Qmin = 0.5 * pi/AggSize
814        variable Qmax = 6 * pi/PrimarySize
815        variable NumQSteps = 200
816        variable PrimarySphereRadius = 6                //this is RADIUS of sphere in voxels, not in angstroms.
817        //Internally, each particle volume is made 10xx10x10, 8 seems to be value when the spheres are exactly touching in xy direction, 9 when in xyz direction.
818        Wave/Z ThreeDVoxelGram = Wave3DwithPrimaryShrunk                //this is voxelgram shrunk to min size 
819        if(!WaveExists(ThreeDVoxelGram) || isInWorkFolder || recalculate3D)
820                //convert to voxelgram
821                IR3T_ConvertToVoxelGram(MassFractalAggregate, PrimarySphereRadius)
822                Wave ThreeDVoxelGram = Wave3DwithPrimaryShrunk                  //this is voxelgram shrunk to min size...
823        endif
824
825        //Calculate pdf intensity
826        //TODO: check these are set correctly: VoxelSize, NumRSteps
827        SetScale /P x, 0, VoxelSize, "A" ,ThreeDVoxelGram
828        SetScale /P y, 0, VoxelSize, "A" ,ThreeDVoxelGram
829        SetScale /P z, 0, VoxelSize, "A" ,ThreeDVoxelGram
830        IR3T_CreatePDFIntensity(ThreeDVoxelGram, IsoValue,  Qmin, Qmax, NumQSteps)
831        //append to graph...
832        Wave PDFQWv
833        Wave PDFIntensityWv
834        IR3A_Append1DInMassFracAgg(PDFIntensityWv,PDFQWv)
835
836        Wave/Z ThreeDVoxelGram = Wave3DwithPrimary              //this is voxelgram with even number of rows/columns/layers.   
837        SetScale /P x, 0, VoxelSize, "A" ,ThreeDVoxelGram
838        SetScale /P y, 0, VoxelSize, "A" ,ThreeDVoxelGram
839        SetScale /P z, 0, VoxelSize, "A" ,ThreeDVoxelGram
840        IR3T_CalcAutoCorelIntensity(ThreeDVoxelGram, IsoValue,  Qmin, Qmax, NumQSteps)
841        //these are autocorrelation calculated intensities...
842        Wave AutoCorIntensityWv
843        Wave AutoCorQWv
844        IR3A_Append1DInMassFracAgg(AutoCorIntensityWv,AutoCorQWv)
845
846        //display the intensity.
847//      DoWIndow IR1_LogLogPlotU
848//      if(V_Flag)
849//              DoWIndow/F IR1_LogLogPlotU
850//              Wave IntWave = root:Packages:Irena_UnifFit:OriginalIntensity
851//              Wave Qwave = root:Packages:Irena_UnifFit:OriginalQvector
852//              if(!WaveExists(IntWave))
853//                      setDataFOlder OldDf
854//                      //this is weird, this should exit in this case...
855//                      return 0
856//              endif
857//              variable InvarModel=areaXY(PDFQWv, PDFIntensityWv )
858//              variable InvarData=areaXY(Qwave, IntWave )
859//              PDFIntensityWv*=InvarData/InvarModel
860//              CheckDisplayed /W=IR1_LogLogPlotU PDFIntensityWv
861//              if(V_flag==0)
862//                      AppendToGraph/W=IR1_LogLogPlotU  PDFIntensityWv vs PDFQWv
863//              endif
864//              ModifyGraph lstyle(PDFIntensityWv)=9,lsize(PDFIntensityWv)=3,rgb(PDFIntensityWv)=(1,16019,65535)
865//              ModifyGraph mode(PDFIntensityWv)=4,marker(PDFIntensityWv)=19
866//              ModifyGraph msize(PDFIntensityWv)=3
867//      else
868//              DoWIndow MassFractalAggDataPlot
869//              if(V_Flag)
870//                      DoWIndow/F MassFractalAggDataPlot
871//              else
872//                      IR3A_Display1DIntensity()
873//              endif
874//              variable InvarModel=areaXY(PDFQWv, PDFIntensityWv )
875//              variable InvarData=areaXY(Qwave, IntWave )
876//              PDFIntensityWv*=InvarData/InvarModel
877//              CheckDisplayed /W=MassFractalAggDataPlot PDFIntensityWv
878//              if(V_flag==0)
879//                      AppendToGraph/W=MassFractalAggDataPlot PDFIntensityWv vs PDFQWv
880//              endif
881//               
882//              ModifyGraph mode(PDFIntensityWv)=0,rgb(PDFIntensityWv)=(0,0,0)         
883//              ModifyGraph log=1
884//              ModifyGraph mirror=1
885//              ShowInfo
886//              String LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Intensity"
887//              Label left LabelStr
888//              LabelStr= "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Q [A\\S-1\\M\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"]"
889//              Label bottom LabelStr
890////            string LegendStr="\\F"+IN2G_LkUpDfltStr("FontType")+"\\Z"+IN2G_LkUpDfltVar("LegendSize")+"\\s(PDFIntensityWv) Mass Fractal Model intensity"
891////            Legend/W=MassFractalAggDataPlot/N=text0/J/F=0/A=MC/X=32.03/Y=38.79 LegendStr
892       
893
894        setDataFOlder OldDf
895       
896end
897
898//******************************************************************************************************************************************************
899//******************************************************************************************************************************************************
900static Function IR3A_Display3DAggregate(AppendToNotebook)
901        variable AppendToNotebook
902
903        string selection = IR3A_FindSelectedAggData()
904        if(strlen(selection)>0)
905                KillWIndow/Z MassFractalAggregateView
906                Wave/Z MassFractalAggregate=$(selection)
907                if(WaveExists(MassFractalAggregate))
908                        IR3A_GizmoViewScatterPlot(MassFractalAggregate)
909                        IR3A_DisplayAggNotebook(MassFractalAggregate, AppendToNotebook)
910                        DoWIndow MassFractalAggDataPlot
911                        if(V_FLag)
912                                AutoPositionWindow /M=1 /R=MassFractalAggDataPlot MassFractalAggregateView
913                                AutoPositionWindow /M=0 /R=MassFractalAggregateView Summary3DAggregate                         
914                        endif
915                else
916                        Print "MassFractalAggregate 3D wave does not exist, cannot do anything"
917                endif
918        endif
919
920end
921//******************************************************************************************************************************************************
922//******************************************************************************************************************************************************
923Function/S IR3A_FindSelectedAggData()
924        //locates which Aggregate data user selected - or returns curren tmodel in Packages.
925        wave/Z/T Stored3DAggregatesPaths= root:Packages:AggregateModeling:Stored3DAggregatesPaths
926        if(!WaveExists(Stored3DAggregatesPaths))
927                abort
928        endif
929        wave/T Stored3DAggregates= root:Packages:AggregateModeling:Stored3DAggregates
930        variable i
931        ControlInfo /W=FractalAggregatePanel StoredAggregates
932        if(V_Value<0)           //no selection...
933                Wave/Z MassFractalAggregate=root:Packages:AggregateModeling:MassFractalAggregate
934                if(WaveExists(MassFractalAggregate))
935                        return "root:Packages:AggregateModeling:MassFractalAggregate"                                   
936                else
937                        Print "MassFractalAggregate 3D wave does not exist, cannot do anything"
938                endif
939        endif
940        string selection
941        if(V_Value<numpnts(Stored3DAggregatesPaths))
942                selection = Stored3DAggregatesPaths     [V_Value]
943        else
944                selection = ""
945        endif
946        if(stringMatch(selection,"Current model"))
947                                Wave/Z MassFractalAggregate=root:Packages:AggregateModeling:MassFractalAggregate
948                                if(WaveExists(MassFractalAggregate))
949                                        return "root:Packages:AggregateModeling:MassFractalAggregate"                                   
950                                else
951                                        Print "MassFractalAggregate 3D wave does not exist, cannot do anything"
952                                endif
953        else
954                                Wave/Z MassFractalAggregate=$(selection+":MassFractalAggregate")
955                                if(WaveExists(MassFractalAggregate))
956                                        return selection+":MassFractalAggregate"                                       
957                                else
958                                        Print "MassFractalAggregate 3D wave does not exist, cannot do anything"
959                                endif
960        endif
961
962end     
963//******************************************************************************************************************************************************
964//******************************************************************************************************************************************************
965//******************************************************************************************************************************************************
966static Function IR3A_DisplayAggNotebook(MassFractalAggregate, AppendToNotebook)
967        wave MassFractalAggregate
968        variable AppendToNotebook
969       
970       
971        String nb = "Summary3DAggregate"
972        DoWIndow Summary3DAggregate
973        if(!V_Flag)
974                NewNotebook/N=$nb/F=1/V=1/K=3/ENCG={1,1}/W=(721,68,1397,644)
975                Notebook $nb defaultTab=36, magnification=125
976                Notebook $nb showRuler=1, rulerUnits=2, updating={1, 1}
977                Notebook $nb newRuler=Normal, justification=0, margins={0,0,468}, spacing={0,0,0}, tabs={}, rulerDefaults={"Helvetica",11,0,(0,0,0)}
978                Notebook $nb ruler=Normal, fSize=16, fStyle=1, text="Summary of Mass Fractal Aggregates\r"
979                Notebook $nb fSize=-1, fStyle=-1, text="\r"
980        else
981                DoWIndow/F Summary3DAggregate
982        endif
983        if(AppendToNotebook)
984                //and append results for aggregate we are looking at
985                NVAR Target_dmin=root:Packages:AggregateModeling:BrFract_dmin
986                NVAR Target_c=root:Packages:AggregateModeling:BrFract_c
987                NVAR Target_df=root:Packages:AggregateModeling:BrFract_df
988                string ResultsLocation = GetWavesDataFolder(MassFractalAggregate, 1 )
989                string OldNote=note(MassFractalAggregate)
990                Notebook $nb selection={endOfFile,endOfFile}
991                Notebook $nb ruler=Normal, text="\r"
992                Notebook $nb ruler=Normal, text="*************************************************************************************************   \r"
993                Notebook $nb ruler=Normal, text="Mass Fractal 3D Aggregate from : "+ResultsLocation+"\r"
994                Notebook $nb ruler=Normal, text="Date & time generated : "+StringByKey("Mass Fractal Aggregate created", OldNote, "=", ";")+"\r"
995                Notebook $nb ruler=Normal, text="****      *********      *************    *******    \r"
996                Notebook $nb ruler=Normal, text="Degree of Aggregation z : \t\t\t"+StringByKey("z", OldNote, "=", ";")+"\r"
997                Notebook $nb ruler=Normal, text="Aggregate Size R : \t\t\t\t\t\t"+StringByKey("R", OldNote, "=", ";")+"\r"
998                Notebook $nb ruler=Normal, text="Short circuit path length p : \t\t"+StringByKey("p", OldNote, "=", ";")+"\r"
999                Notebook $nb ruler=Normal, text="Connective path length s : \t\t\t"+StringByKey("s", OldNote, "=", ";")+"\r"
1000                Notebook $nb ruler=Normal, text="Mass fractal dimension df : \t\t"+StringByKey("df", OldNote, "=", ";")+"\r"
1001                Notebook $nb ruler=Normal, text="Minimum dimension dmin : \t\t"+StringByKey("dmin", OldNote, "=", ";")+"\r"
1002                Notebook $nb ruler=Normal, text="Connectivity dimension c : \t\t\t"+StringByKey("c", OldNote, "=", ";")+"\r"
1003                Notebook $nb ruler=Normal, text="****      *********      *************    *******    \r"
1004                Notebook $nb ruler=Normal, text="Sticking Probability :\t\t\t\t\t\t"+StringByKey("StickingProbability", OldNote, "=", ";")+"\r"
1005                Notebook $nb ruler=Normal, text="Sticking Method : \t\t\t\t\t\t\t"+StringByKey("StickingMethod", OldNote, "=", ";")+"\r"
1006                Notebook $nb ruler=Normal, text="Multi Particle Attraction : \t\t\t\t"+StringByKey("MultiParticleAttraction", OldNote, "=", ";")+"\r"
1007                Notebook $nb ruler=Normal, text="True Sticking Probability : \t\t\t"+StringByKey("True Sticking Probability", OldNote, "=", ";")+"\r"
1008                Notebook $nb ruler=Normal, text="Maximum Path Length : \t\t\t\t"+StringByKey("MaximumPathLength", OldNote, "=", ";")+"\r"
1009                Notebook $nb ruler=Normal, text="Max Number Of Paths Per End : \t\t"+StringByKey("MaxNumberOfPathsPerEnd", OldNote, "=", ";")+"\r"
1010                Notebook $nb ruler=Normal, text="Number Of End particles : \t\t\t"+StringByKey("NumberOfEnds", OldNote, "=", ";")+"\r"
1011                Notebook $nb ruler=Normal, text="******          Target    &    resulting values         ************   \r"
1012                Notebook $nb ruler=Normal, text="dmin : \t\t Target : "+num2str(IN2G_roundSignificant(Target_dmin,3))+"\t\tResult : "+StringByKey("dmin", OldNote, "=", ";")+"\r"
1013                Notebook $nb ruler=Normal, text="c : \t\t Target : "+num2str(IN2G_roundSignificant(Target_c,3))+"\t\tResult : "+StringByKey("c", OldNote, "=", ";")+"\r"
1014                Notebook $nb ruler=Normal, text="df : \t\t Target : "+num2str(IN2G_roundSignificant(Target_df,3))+"\t\tResult : "+StringByKey("df", OldNote, "=", ";")+"\r"
1015                Notebook $nb ruler=Normal, text="*************************************************************************************************   \r"
1016                //      NoteText="Mass Fractal Aggregate created="+date()+", "+time()+";z="+num2str(DegreeOfAggregation)+";StickingProbability="+num2str(StickingProbability)+";StickingMethod="+num2str(AllowedNearDistance)
1017                //      NoteText+=";R="+num2str(RValue)+";Rprimary="+num2str(RgPrimary)+";p="+num2str(pValue)
1018                //      NoteText+=";RxRgPrimaryValue="+num2str(RValue*PrimaryDiameter)+";s="+num2str(sValue)+";df="+num2str(IN2G_roundSignificant(dfValue,4))+";dmin="+num2str(IN2G_roundSignificant(dminValue,4))
1019                //      NoteText+=";c="+num2str(IN2G_roundSignificant(cValue,4))+";True Sticking Probability="+num2str(100*DegreeOfAggregation/AttemptValue)
1020                //      NoteText+=";MultiParticleAttraction="+MultiParticleAttraction+";MaximumPathLength="+num2str(MaxPathLength)+";MaxNumberOfPathsPerEnd="+num2str(MaxNumPaths)+";NumberOfEnds="+num2str(NumStarts)+";"
1021
1022                DoWIndow MassFractalAggregateView
1023                if(V_Flag)
1024                        Notebook $nb scaling={80,80}, picture={MassFractalAggregateView, -5, 1 }
1025                        Notebook $nb ruler=Normal, text="\r"
1026                endif
1027                DoWIndow MassFractalAggDataPlot
1028                if(V_Flag)
1029                        Notebook $nb scaling={75,75}, picture={MassFractalAggDataPlot, 2, 1 }
1030                        Notebook $nb ruler=Normal, text="\r"
1031                endif
1032        endif
1033        //      R - aggregate size
1034        //      df - Mass fractal dimension of the aggregate
1035        //      p - short circuit path length
1036        //      s - connective path length
1037        //      dmin - minimum dimension of the aggregate
1038        //      c - connectivity dimension of the aggregate
1039        //      s - connective path length of the aggregate
1040
1041end
1042
1043//******************************************************************************************************************************************************
1044//******************************************************************************************************************************************************
1045//******************************************************************************************************************************************************
1046
1047
1048
1049static Function IR3A_Grow1MassFractAgreg()
1050
1051        DFref oldDf= GetDataFolderDFR()
1052
1053        IR3A_InitializeMassFractAgg()
1054        SetDataFolder root:Packages:AggregateModeling
1055        KillWIndow/Z MassFractalAggregateView
1056        KillWIndow/Z MassFractalAggDataPlot
1057        KillWindow/Z AggStoredResultsOverview
1058               
1059        NVAR DegreeOfAggregation=root:Packages:AggregateModeling:DegreeOfAggregation
1060        NVAR StickingProbability=root:Packages:AggregateModeling:StickingProbability
1061        NVAR NumberOfTestPaths=root:Packages:AggregateModeling:NumberOfTestPaths
1062        NVAR AllowedNearDistance=root:Packages:AggregateModeling:AllowedNearDistance
1063                // Get the starting position of the aggregate
1064        Make/n=(DegreeOfAggregation,3)/O MassFractalAggregate=0         // It starts at 0,0,0
1065        Make/n=(DegreeOfAggregation,4)/O endpoints                                                      //List of end points
1066        make/N=(DegreeOfAggregation)/O Distances                                                        // Distance between existing particles & new one. Needed by MakeAgg
1067        variable StartTicks=ticks
1068        variable Failed
1069        print time()+"  Started Growing Aggregate and evaluation its structure "
1070        IR3A_MakeAgg(DegreeOfAggregation,MassFractalAggregate,StickingProbability,AllowedNearDistance)          // Agg is made with DegreeOfAggregation particles
1071        //Failed = IR3A_Ends(MassFractalAggregate, breakOnFail)
1072        //if(!Failed)
1073                //IR3A_Reted(endpoints)
1074                //IR3A_Path(NumberOfTestPaths)
1075        IR3A_EvaluateAggregateUsingMT()
1076        Failed = IR3A_CalculateParametersMT()
1077        if(!failed)
1078                IR3A_GizmoViewScatterPlot(MassFractalAggregate)
1079        else
1080                Print "Failed to grow meaningful aggregate)"
1081        endif
1082        print time()+"  Finished, done in "+num2str((ticks-StartTicks)/60)+" seconds"   
1083        //else
1084        //      print time()+"  Failed, Aggregate is too compact "     
1085        //endif
1086        setDataFOlder OldDf
1087        return Failed
1088End
1089//******************************************************************************************************************************************************
1090//******************************************************************************************************************************************************
1091//this is code to evaluate paths and statistics using multithreading and lot more unique paths...
1092
1093//******************************************************************************************************************************************************
1094//******************************************************************************************************************************************************
1095static FUnction IR3A_EvaluateAggregateUsingMT()
1096
1097        DFref oldDf= GetDataFolderDFR()
1098        SetDataFolder root:Packages:AggregateModeling
1099        print "****   Evaluating Aggregate structure    *****"
1100        print "This may take anywhere from seconds to hours, depending on size and complexity of the aggregate."
1101        wave ListOfNeighbors
1102        wave NumberOfNeighbors
1103        NVAR NumberOfTestPaths = root:Packages:AggregateModeling:NumberOfTestPaths
1104        NVAR DegreeOfAggregation = root:Packages:AggregateModeling:DegreeOfAggregation
1105
1106        variable timerRefNum=StartMSTimer
1107                //Make/O/N=(dimsize(ListOfNeighbors,0)) ListOfStarts
1108                //Duplicate/O ListOfNeighbors, ListOfNeighborsForStarts
1109                //      variable i, ij
1110                //      For(i=numpnts(NumberOfNeighbors)-1;i>=0;i-=1)           //this iterates over all points in the aggregate, it goes in order of how points were added.
1111                //              if(NumberOfNeighbors[i]>1.5)                            //this point is not end, delete.
1112                //                      DeletePoints/M=0 i, 1, ListOfStarts, ListOfNeighborsForStarts
1113                //              else
1114                //                      ListOfStarts[i]=i
1115                //              endif
1116                //      endFOR 
1117
1118        //the loop above can done easier by this:
1119        //extract indexes for end points. ListOfStarts = endpoints.
1120        Extract /O /INDX NumberOfNeighbors, ListOfStarts, NumberOfNeighbors<2
1121        Duplicate/O ListOfStarts, ListOfNeighborsForStarts
1122        //and these are neighbors fo each end point.
1123        ListOfNeighborsForStarts= ListOfNeighbors[ListOfStarts[p][0]]
1124
1125        //at this moment we have list of start points ListOfStarts and list of next points for each start point.
1126        //we can now start independent thread for each of pairs of ListOfStarts[i] as prior point and NumberOfNeighborsForStarts[i] as current point.
1127        Make/WAVE/O/N=(dimsize(ListOfStarts,0)) ListOfUniquePathListWvs
1128       
1129        multithread ListOfUniquePathListWvs = IR3A_MT_WalkPathThread(ListOfNeighbors,NumberOfNeighbors,ListOfStarts[p], NumberOfTestPaths, DegreeOfAggregation)
1130
1131        variable microSeconds = StopMSTimer(timerRefNum)
1132
1133        Print "***** Evaluted  aggregate in sec : "+num2str(microSeconds/1e6)   //takes needless time..
1134
1135        setDataFOlder OldDf
1136
1137end
1138//******************************************************************************************************************************************************
1139//******************************************************************************************************************************************************
1140
1141static Function IR3A_CalculateParametersMT()
1142
1143        DFref oldDf= GetDataFolderDFR()
1144        SetDataFolder root:Packages:AggregateModeling
1145       
1146        Wave/WAVE ListOfUniquePathListWvs
1147        WAVE ListOfStarts = root:Packages:AggregateModeling:ListOfStarts
1148        NVAR dfValue
1149        NVAR RValue
1150        NVAR pValue
1151        NVAR sValue
1152        NVAR cValue
1153        NVAR dminValue
1154        NVAR TrueStickingProbability
1155        NVAR StickingProbability=root:Packages:AggregateModeling:StickingProbability
1156        NVAR RgPrimary=root:Packages:AggregateModeling:RgPrimary
1157        NVAR AllowedNearDistance
1158        NVAR AttemptValue
1159        NVAR NumberOfTestPaths = root:Packages:AggregateModeling:NumberOfTestPaths
1160        NVAR AllowedNearDistance=root:Packages:AggregateModeling:AllowedNearDistance
1161        NVAR RxRgPrimaryValue=root:Packages:AggregateModeling:RxRgPrimaryValue
1162        NVAR DegreeOfAggregation=root:Packages:AggregateModeling:DegreeOfAggregation
1163        SVAR MultiParticleAttraction = root:Packages:AggregateModeling:MultiParticleAttraction //"Neutral;Attractive;Repulsive;Not allowed;"
1164
1165        variable failed=0
1166        //here we will evaluate all paths, they arrive using wave of waves : ListOfUniquePathListWvs
1167        //each item  in this wave is wavereference to free wave in memory with one path.
1168        variable i,j,ij, SumPL2, SumPL, NumPathsFOund, AveragePath, MaxPathLength, MaxNumPaths, tmpNumUsefulPaths
1169        MaxPathLength = 0
1170        Wave ListOfStarts
1171        For(ij=0;ij<numpnts(ListOfStarts);ij+=1)
1172                Wave UniquePathList=ListOfUniquePathListWvs[ij]
1173                //calculate sum(path lengths^2)/sum(path Lengths)
1174                //collect also length of longest path and number of paths per starting point.
1175                tmpNumUsefulPaths = 0
1176                For(i=0;i<DimSize(UniquePathList,1);i+=1)
1177                        //get one path from the table
1178                        Duplicate /free/R=[][i]  UniquePathList, TempWv
1179                        //extract non NaN points from this.
1180                        extract/Free TempWv, TempWv, TempWv<65534
1181                        //if it is longer than 2, add to list of waves
1182                        if(numpnts(TempWv)>2)
1183                                SumPL+=numpnts(TempWv)
1184                                SumPL2+=numpnts(TempWv)^2
1185                                NumPathsFound+=1
1186                                MaxPathLength = max(MaxPathLength, numpnts(TempWv))
1187                                tmpNumUsefulPaths+=1
1188                        endif
1189                endfor
1190                //this is max number of paths we have found in the system.
1191                MaxNumPaths = max(MaxNumPaths,tmpNumUsefulPaths)
1192        endfor
1193        //weighted avergae path.
1194        AveragePath = SumPL2/SumPL
1195        //this is called also pValue
1196        pValue = AveragePath
1197        //this is useful to evaluate also...
1198        variable PathperEndPoint=NumPathsFound/numpnts(ListOfStarts)
1199        //this will release from memory those free waves with paths, saves Igor memory..
1200        KillWaves/Z ListOfUniquePathListWvs
1201       
1202        //Next we need to evaluate size of the aggregate.
1203        Wave/Z NumNeighbors=root:Packages:AggregateModeling:NumberOfNeighbors
1204        if(WaveExists(NumNeighbors))
1205                Wave Agg=root:Packages:AggregateModeling:MassFractalAggregate
1206                //next we extract from list of neighbors only the short = 1 items items. These are end points.
1207                Extract /FREE /INDX NumNeighbors, AggregateEndsIndex, NumNeighbors<2
1208                // AggregateEndsIndex indexes where are end points in Aggregate
1209                variable NumEnds=numpnts(AggregateEndsIndex)
1210                variable numcomb=binomial(NumEnds,2)
1211                variable REnd, Rsum, cnt, FInx, SInx
1212                make/Free/N=(numcomb) EndsDistances
1213                For(i=0;i<NumEnds;i+=1)
1214                        For(j=i+1;j<NumEnds;j+=1)
1215                                FInx=AggregateEndsIndex[i]
1216                                SInx=AggregateEndsIndex[j]
1217                                EndsDistances[cnt]=sqrt((Agg[FInx][0]-Agg[SInx][0])^2+(Agg[FInx][1]-Agg[SInx][1])^2+(Agg[FInx][2]-Agg[SInx][2])^2)
1218                                cnt+=1
1219                        endfor
1220                endfor
1221                //For(i=0;i<numcomb;i+=1)
1222                //      REnd+=EndsDistances[i]^2
1223                //      Rsum+=EndsDistances[i]
1224                //endfor
1225                duplicate/Free EndsDistances, EndsDistances2
1226                EndsDistances2=EndsDistances^2
1227                REnd = sum(EndsDistances2)
1228                Rsum=sum(EndsDistances)
1229                REnd/=RSum
1230                RValue=Rend
1231                dfValue = log(DegreeOfAggregation)/log(REnd)
1232                cValue=ln(DegreeOfAggregation)/ln(AveragePath)
1233                dminValue=dfValue/cValue
1234                sValue=(exp(ln(DegreeOfAggregation)/dminValue))
1235        endif   
1236
1237        TrueStickingProbability = 100*DegreeOfAggregation/AttemptValue
1238        variable PrimaryDiameter = 2*sqrt(5/3)*RgPrimary
1239        variable NumStarts=numpnts(ListOfStarts)
1240        RxRgPrimaryValue = RValue*PrimaryDiameter
1241        print "***** Results listing ******"
1242        print "Total number of end points is : "+num2str(NumStarts)
1243        if(NumStarts<4)
1244                print "Warning : This particle is too compact (too few ends) to make any sense"
1245                failed = 1
1246        endif
1247        //print "Possible number of paths is : " +num2str(NumStarts*(NumStarts-1))
1248        print "Total number of evaluated paths is : "+num2str(NumPathsFound)
1249        print "Paths per end point is : "+num2str(PathperEndPoint)
1250        if(MaxNumPaths > NumberOfTestPaths-2)
1251                print "Warning : Some path/end point numbers were limited by max endpoint choice. Results may not be valid, increase No of test paths for these conditions"
1252        endif
1253        print "Maximum length path is: "+num2str(MaxPathLength)
1254        print "Weighted Average Path ="+num2str(pValue)
1255        Print "Aggregate Rg [relative] = "+num2str(REnd)
1256        Print "R primary particles [A] = "+num2str(PrimaryDiameter/2)
1257        Print "Aggregate R [Angstroms] = "+num2str(RValue*PrimaryDiameter)
1258        Print "z = "+num2str(DegreeOfAggregation)
1259        Print "p = "+num2str(pValue)
1260        Print "c = "+num2str(cValue)
1261        Print "s = "+num2str(sValue)
1262        Print "df = "+num2str(dfValue)
1263        Print "dmin = "+num2str(dminValue)
1264        Print "True Sticking Probability = "+num2str(100*DegreeOfAggregation/AttemptValue)+"%"
1265        //Print "R [primary particles]= "+num2str(REnd)
1266
1267        //appned note to MassFractalAggregate
1268        string NoteText
1269        NoteText="Mass Fractal Aggregate created="+date()+", "+time()+";z="+num2str(DegreeOfAggregation)+";StickingProbability="+num2str(StickingProbability)+";StickingMethod="+num2str(AllowedNearDistance)
1270        NoteText+=";R="+num2str(RValue)+";Rprimary="+num2str(RgPrimary)+";p="+num2str(pValue)
1271        NoteText+=";RxRgPrimaryValue="+num2str(RValue*PrimaryDiameter)+";s="+num2str(sValue)+";df="+num2str(IN2G_roundSignificant(dfValue,4))+";dmin="+num2str(IN2G_roundSignificant(dminValue,4))
1272        NoteText+=";c="+num2str(IN2G_roundSignificant(cValue,4))+";True Sticking Probability="+num2str(100*DegreeOfAggregation/AttemptValue)
1273        NoteText+=";MultiParticleAttraction="+MultiParticleAttraction+";MaximumPathLength="+num2str(MaxPathLength)+";MaxNumberOfPathsPerEnd="+num2str(MaxNumPaths)+";NumberOfEnds="+num2str(NumStarts)+";"
1274
1275        Wave MassFractalAggregate
1276        Note MassFractalAggregate, NoteText
1277       
1278        setDataFOlder OldDf
1279        return failed
1280end
1281
1282//******************************************************************************************************************************************************
1283//******************************************************************************************************************************************************
1284//
1285Threadsafe Function/Wave IR3A_MT_WalkPathThread(ListOfNeighbors,NumberOfNeighbors,startingPoint, NumberOfTestPaths, DegreeOfAggregation)
1286        wave ListOfNeighbors, NumberOfNeighbors
1287        variable startingPoint, NumberOfTestPaths, DegreeOfAggregation
1288        //Advanced topics, Wave Reference MultiThread Example
1289        DFREF dfSav= GetDataFolderDFR()
1290        // Create a free data folder and set it as the current data folder
1291        SetDataFolder NewFreeDataFolder()
1292        variable/g UniquePathListInx
1293        UniquePathListInx = 0
1294        make/Free/W/U/N=(DegreeOfAggregation,NumberOfTestPaths) UniquePathList
1295        UniquePathList = NaN
1296        String PriorPathList
1297        variable currentPoint
1298        PriorPathList=num2str(startingPoint)+";"                                //this is the path start as list...
1299        currentPoint = ListOfNeighbors[startingPoint][0]                //this is next point since end point has only one neighbor... 
1300        IR3A_MT_NextPathStep(UniquePathList, ListOfNeighbors,NumberOfNeighbors,currentPoint, startingPoint, PriorPathList)
1301        // Restore the current data folder
1302        SetDataFolder dfSav
1303        return UniquePathList
1304end
1305//******************************************************************************************************************************************************
1306//******************************************************************************************************************************************************
1307
1308Threadsafe Function IR3A_MT_NextPathStep(UniquePathList, ListOfNeighbors,NumberOfNeighbors,currentPoint, priorPoint, PriorPathList)
1309//Function MT_NextPathStep(ListOfNeighbors,NumberOfNeighbors,currentPoint, priorPoint, PriorPathList, OrderNumber)
1310        wave UniquePathList, ListOfNeighbors
1311        wave NumberOfNeighbors
1312        variable priorPoint, currentPoint
1313        string PriorPathList
1314
1315        variable pointsFound, i, NewCurrentPoint, NewpriorPoint
1316        string CurrentPathList
1317       
1318        //here we need to check, if the new point is already in the list.
1319        //if it is, this is loop and we shoudl simply return back and refuse to go here. May have to make more decisions later...
1320        if(StringMatch(PriorPathList, "*;"+num2str(currentPoint)+";*"))
1321                //print "This was circle:"+PriorPathList+", repeting point would be :"+num2str(currentPoint)
1322                return 0
1323        endif
1324        NVAR UniquePathListInx
1325        if(UniquePathListInx>=dimsize(UniquePathList,1))
1326                //we have reached max number of paths allowed for this start point...
1327                return 0
1328        endif
1329       
1330        PriorPathList+=num2str(currentPoint)+";"
1331        CurrentPathList= PriorPathList
1332#if(IgorVersion()>8.99)
1333        MatrixOP/Free ListOfNeighborsRow = zapNaNs(replace(row(ListOfNeighbors,currentPoint),priorPoint,NaN)))
1334        //how many points are left?
1335        pointsFound = numpnts(ListOfNeighborsRow)
1336#else
1337        MatrixOP/Free ListOfNeighborsRow = replace(row(ListOfNeighbors,currentPoint),priorPoint,NaN))
1338        //how many points are not Nans?
1339        Extract /FREE /O ListOfNeighborsRow, ListOfNeighborsRow, numtype(ListOfNeighborsRow)==0
1340        pointsFound = numpnts(ListOfNeighborsRow)
1341#endif
1342        //decisions what to do...
1343        if(pointsFound==0)              //this is end point, I just came from the only neighbor this has
1344                //print "One complete path is : "+CurrentPathList
1345                IR3A_MT_WriteOutPathString(UniquePathList, CurrentPathList)
1346                //what we really need here is write out the path into some kind of final container when we get here.
1347        elseif(pointsFound==1)          //this is connecting segment, we need to go in next segment and see what is there
1348                NewpriorPoint = currentPoint
1349                NewCurrentPoint = ListOfNeighborsRow[0]         //this is new point left.
1350                //go in next step
1351                IR3A_MT_NextPathStep(UniquePathList, ListOfNeighbors,NumberOfNeighbors,NewCurrentPoint, NewpriorPoint, CurrentPathList)
1352        elseif(pointsFound<4)                                                           //this is junction with up to 3 new neighbors, let's take all here...
1353                For(i=0;i<pointsFound;i+=1)
1354                        NewpriorPoint = currentPoint
1355                        NewCurrentPoint = ListOfNeighborsRow[i]         //this is new point left.
1356                        CurrentPathList = PriorPathList
1357                        IR3A_MT_NextPathStep(UniquePathList, ListOfNeighbors,NumberOfNeighbors,NewCurrentPoint, NewpriorPoint, CurrentPathList)
1358                endfor
1359        else                                                                                            //more points than 3, let's take first three from here... 
1360                For(i=0;i<3;i+=1)
1361                        NewpriorPoint = currentPoint
1362                        //NewCurrentPoint = ListOfNeighborsRow[floor(pointsFound*(enoise(0.5)+0.5))]            //this is new point left.
1363                        NewCurrentPoint = ListOfNeighborsRow[i]                                                         //this is new point left.
1364                        CurrentPathList = PriorPathList
1365                        IR3A_MT_NextPathStep(UniquePathList, ListOfNeighbors,NumberOfNeighbors,NewCurrentPoint, NewpriorPoint, CurrentPathList)
1366                endfor
1367        endif   
1368end
1369//******************************************************************************************************************************************************
1370//******************************************************************************************************************************************************
1371
1372Threadsafe Function IR3A_MT_WriteOutPathString(UniquePathList, PathStr)
1373        Wave UniquePathList
1374        string PathStr
1375       
1376        NVAR UniquePathListInx
1377       
1378        if(DimSize(UniquePathList, 0)<ItemsInList(PathStr))// || UniquePathListInx > (DimSize(UniquePathList,1)-1))
1379                //the length here is now set to DegreeOfAggregation which means that we cannot have longer path than number of particles...
1380                //variable WvMaxlength, WvNewlength
1381                //WvMaxlength = DimSize(UniquePathList, 0)
1382                //WvNewlength = max(WvMaxlength, ItemsInList(PathStr)+10)
1383                //redimension/N=(WvNewlength,-1) UniquePathList
1384        endif
1385        //now, we need to write values in...
1386        UniquePathList[][UniquePathListInx] = NaN       
1387        UniquePathList[][UniquePathListInx] = str2num(StringFromList(p,PathStr, ";"))   
1388
1389        UniquePathListInx = UniquePathListInx+1
1390end
1391
1392//******************************************************************************************************************************************************
1393//******************************************************************************************************************************************************
1394
1395
1396//******************************************************************************************************************************************************
1397//******************************************************************************************************************************************************
1398//******************************************************************************************************************************************************
1399
1400static Function IR3A_GrowAggregate()
1401        DFref oldDf= GetDataFolderDFR()
1402
1403        IR3A_InitializeMassFractAgg()
1404        SetDataFolder root:Packages:AggregateModeling
1405
1406        NVAR DegreeOfAggregation=root:Packages:AggregateModeling:DegreeOfAggregation
1407        NVAR StickingProbability=root:Packages:AggregateModeling:StickingProbability
1408        NVAR AllowedNearDistance=root:Packages:AggregateModeling:AllowedNearDistance
1409        variable DegreeOfAggregationL=DegreeOfAggregation
1410        VARIABLE StickingProbabilityL=StickingProbability
1411        Prompt DegreeOfAggregationL, "Enter the size of the aggregate (250)"
1412        Prompt StickingProbabilityL, "Enter StickingProbabilitying probability (1 - 100)"
1413        DoPrompt/Help="Basic parameters" "Input model parameters", DegreeOfAggregationL, StickingProbabilityL
1414        DegreeOfAggregation=DegreeOfAggregationL
1415        StickingProbability=StickingProbabilityL
1416        // Get the starting position of the aggregate
1417        Make/n=(DegreeOfAggregation,3)/O MassFractalAggregate=0         // It starts at 0,0,0
1418        make/N=(DegreeOfAggregation)/O Distances        // Distance between existing particles & new one. Needed by MakeAgg
1419        IR3A_MakeAgg(DegreeOfAggregation,MassFractalAggregate,StickingProbability,AllowedNearDistance)          // Agg is made with DegreeOfAggregation particles
1420        setDataFOlder OldDf
1421End
1422//******************************************************************************************************************************************************
1423//******************************************************************************************************************************************************
1424static Function IR3A_InitializeMassFractAgg()
1425
1426        //IN2G_PrintDebugStatement(IrenaDebugLevel, 5,"")
1427        DFref oldDf= GetDataFolderDFR()
1428
1429        NewDataFolder/O/S root:Packages
1430        NewDataFolder/O/S root:Packages:AggregateModeling
1431        IR3A_MakeNBROffsetList()
1432        string/g ListOfVariables
1433        string/g ListOfStrings
1434        //here define the lists of variables and strings needed, separate names by ;...
1435        ListOfVariables="DegreeOfAggregation;StickingProbability;NumberOfTestPaths;BoxSize;RgPrimary;AllowedNearDistance;NUmberOfTestAggregates;"
1436        ListOfVariables+="pValue;dfValue;RValue;RxRgPrimaryValue;cValue;dminValue;sValue;AttemptValue;TrueStickingProbability;"
1437        ListOfVariables+="SelectedLevel;SelectedQlevel;SelectedBlevel;CurrentResults;StoredResults;"
1438        ListOfVariables+="BrFract_G2;BrFract_Rg2;BrFract_B2;BrFract_P2;BrFract_G1;BrFract_Rg1;BrFract_B1;BrFract_P1;BrFract_dmin;"
1439        ListOfVariables+="BrFract_c;BrFract_z;BrFract_fBr;BrFract_fM;BrFract_df;"
1440        ListOfStrings="SlectedBranchedLevels;Model;BrFract_ErrorMessage;MultiParticleAttraction;"
1441        Make/O/N=1/T Stored3DAggregates, Stored3DAggregatesPaths
1442        Make/O/N=1 Stored3DAggSelections
1443        Wave/T Stored3DAggregates
1444        Stored3DAggregates[0] = "Current model"
1445        variable i
1446        //and here we create them
1447        for(i=0;i<itemsInList(ListOfVariables);i+=1)   
1448                IN2G_CreateItem("variable",StringFromList(i,ListOfVariables))
1449        endfor                                                                                         
1450        for(i=0;i<itemsInList(ListOfStrings);i+=1)     
1451                IN2G_CreateItem("string",StringFromList(i,ListOfStrings))
1452        endfor 
1453        NVAR DegreeOfAggregation
1454        if(DegreeOfAggregation<10)
1455                DegreeOfAggregation = 250
1456        endif
1457        NVAR StickingProbability
1458        if(StickingProbability<1 || StickingProbability>100)
1459                StickingProbability = 75
1460        endif
1461        NVAR NumberOfTestPaths
1462        if(NumberOfTestPaths<1000)
1463                NumberOfTestPaths = 2500
1464        endif
1465        NVAR RgPrimary
1466        if(RgPrimary<5)         //Rg of primary partice (voxel size in the model).
1467                RgPrimary = 10
1468        endif
1469        NVAR AllowedNearDistance
1470        if(AllowedNearDistance<0.9 || AllowedNearDistance>3.1)
1471                AllowedNearDistance = 3                         //nearest neighbor distacne squared, for in line neighbors (x,y,z dir) = 1^2, for in plane neighbors is 2 (1^2+1^2), and for in body neighbors is 3 
1472        endif
1473        NVAR CurrentResults
1474        NVAR StoredResults
1475        if(CurrentResults+StoredResults!=1)
1476                StoredResults =1
1477                CurrentResults = 0
1478        endif
1479        SVAR Model
1480        Model = "Branched mass fractal"
1481        SVAR MultiParticleAttraction
1482        if(strlen(MultiParticleAttraction)<1)
1483                MultiParticleAttraction = "Neutral"
1484        endif
1485
1486
1487        NVAR NUmberOfTestAggregates
1488        if(NUmberOfTestAggregates<5)
1489                NUmberOfTestAggregates = 5
1490        endif
1491//      NVAR gdf, gR
1492//      variable/g gp = mom2,gc=ln(DegreeOfAggregation)/ln(gp),gdmin=gdf/gc,gs=round(exp(ln(DegreeOfAggregation)/gdmin))
1493               
1494        setDataFOlder OldDf
1495end
1496//******************************************************************************************************************************************************
1497//******************************************************************************************************************************************************
1498//******************************************************************************************************************************************************
1499
1500static Function IR3A_MakeNBROffsetList()
1501
1502        Make/n=(26,3)/O nghbrOfsetList
1503        Wave nghbrOfsetList = nghbrOfsetList
1504                nghbrOfsetList[0][0] =+1;nghbrOfsetList[0][1] = 0;nghbrOfsetList[0][2] = 0
1505                nghbrOfsetList[1][0] =-1;nghbrOfsetList[1][1] = 0;nghbrOfsetList[1][2] = 0
1506                nghbrOfsetList[2][0] = 0;nghbrOfsetList[2][1] =+1;nghbrOfsetList[2][2] = 0
1507                nghbrOfsetList[3][0] = 0;nghbrOfsetList[3][1] =-1;nghbrOfsetList[3][2] = 0
1508                nghbrOfsetList[4][0] =+1;nghbrOfsetList[4][1] =+1;nghbrOfsetList[4][2] = 0
1509                nghbrOfsetList[5][0] =-1;nghbrOfsetList[5][1] =-1;nghbrOfsetList[5][2] = 0
1510                nghbrOfsetList[6][0] =-1;nghbrOfsetList[6][1] =+1;nghbrOfsetList[6][2] = 0
1511                nghbrOfsetList[7][0] =+1;nghbrOfsetList[7][1] =-1;nghbrOfsetList[7][2] = 0
1512                nghbrOfsetList[8][0] = 0;nghbrOfsetList[8][1] = 0;nghbrOfsetList[8][2] =+1
1513                nghbrOfsetList[9][0] =+1;nghbrOfsetList[9][1] = 0;nghbrOfsetList[9][2] =+1
1514                nghbrOfsetList[10][0]=-1;nghbrOfsetList[10][1]= 0;nghbrOfsetList[10][2]=+1
1515                nghbrOfsetList[11][0]= 0;nghbrOfsetList[11][1]=+1;nghbrOfsetList[11][2]=+1
1516                nghbrOfsetList[12][0]= 0;nghbrOfsetList[12][1]=-1;nghbrOfsetList[12][2]=+1
1517                nghbrOfsetList[13][0]=+1;nghbrOfsetList[13][1]=+1;nghbrOfsetList[13][2]=+1
1518                nghbrOfsetList[14][0]=-1;nghbrOfsetList[14][1]=-1;nghbrOfsetList[14][2]=+1
1519                nghbrOfsetList[15][0]=-1;nghbrOfsetList[15][1]=+1;nghbrOfsetList[15][2]=+1
1520                nghbrOfsetList[16][0]=+1;nghbrOfsetList[16][1]=-1;nghbrOfsetList[16][2]=+1
1521                nghbrOfsetList[17][0]= 0;nghbrOfsetList[17][1]= 0;nghbrOfsetList[17][2]=-1
1522                nghbrOfsetList[18][0]=+1;nghbrOfsetList[18][1]= 0;nghbrOfsetList[18][2]=-1
1523                nghbrOfsetList[19][0]=-1;nghbrOfsetList[19][1]= 0;nghbrOfsetList[19][2]=-1
1524                nghbrOfsetList[20][0]= 0;nghbrOfsetList[20][1]=+1;nghbrOfsetList[20][2]=-1
1525                nghbrOfsetList[21][0]= 0;nghbrOfsetList[21][1]=-1;nghbrOfsetList[21][2]=-1
1526                nghbrOfsetList[22][0]=+1;nghbrOfsetList[22][1]=+1;nghbrOfsetList[22][2]=-1
1527                nghbrOfsetList[23][0]=-1;nghbrOfsetList[23][1]=-1;nghbrOfsetList[23][2]=-1
1528                nghbrOfsetList[24][0]=-1;nghbrOfsetList[24][1]=+1;nghbrOfsetList[24][2]=-1
1529                nghbrOfsetList[25][0]=+1;nghbrOfsetList[25][1]=-1;nghbrOfsetList[25][2]=-1
1530
1531end
1532//******************************************************************************************************************************************************
1533//Function ProfilingRun()
1534//      DFref oldDf= GetDataFolderDFR()
1535
1536//      SetDataFolder root:Packages:AggregateModeling
1537//      AggMod_Initialize()
1538//
1539//      NVAR gGL
1540//      NVAR DegreeOfAggregation
1541//      NVAR StickingProbability
1542//      NVAR NumberOfTestPaths
1543//     
1544//      //variable GL=gGL, DegreeOfAggregation=DegreeOfAggregation, StickingProbability=StickingProbability, NumberOfTestPaths=NumberOfTestPaths
1545//      //Prompt DegreeOfAggregation, "Enter the size of the aggregate (250)"   // Degree of aggregation, z
1546//      //Prompt StickingProbability, "Enter StickingProbabilitying probability (1 - 100)"      // SP = 100% for DLA; less for RLA
1547//      //Prompt NumberOfTestPaths,"Enter number of paths (10000)."             // More paths = more accuracy
1548//      //DoPrompt/Help="Basic parameters" "Input model parameters", DegreeOfAggregation, StickingProbability, NumberOfTestPaths
1549//
1550//      gGL=500;DegreeOfAggregation=500;StickingProbability=40;NumberOfTestPaths=500;
1551//      // Get the starting position of the aggregate
1552//      Make/n=(DegreeOfAggregation,3)/O Agg=0          // It starts at 0,0,0
1553//      Make/n=(DegreeOfAggregation,4)/O endpoints
1554//      make/N=(DegreeOfAggregation)/O Distances        // Distance between existing particles & new one. Needed by MakeAgg
1555//      variable StartTicks=ticks
1556//      print "Started Run All"
1557//              MakeAgg(DegreeOfAggregation,Agg,StickingProbability)            // Agg is made with DegreeOfAggregation particles
1558//              Ends(agg)
1559//              Reted(endpoints)
1560//              Path(NumberOfTestPaths)
1561//              Execute("GizmoViewAggregate()")
1562//      print "Finished, done in "+num2str((ticks-StartTicks)/60)+" seconds"   
1563//      setDataFOlder OldDf
1564//
1565//end
1566
1567//******************************************************************************************************************************************************
1568//******************************************************************************************************************************************************
1569//******************************************************************************************************************************************************
1570//******************************************************************************************************************************************************
1571//******************************************************************************************************************************************************
1572
1573static Function IR3A_MakeAgg(DegreeOfAggregation,MassFractalAggregate,StickingProbability, AllowedNearDistance)
1574        variable DegreeOfAggregation,StickingProbability, AllowedNearDistance
1575        wave MassFractalAggregate
1576       
1577        DFref oldDf= GetDataFolderDFR()
1578
1579        SetDataFolder root:Packages:AggregateModeling
1580        variable NumParticles=dimsize(MassFractalAggregate,0)
1581        //these are waves which must exist...
1582        Wave Distances = root:Packages:AggregateModeling:Distances
1583        NVAR AttemptValue = root:Packages:AggregateModeling:AttemptValue
1584        AttemptValue=1
1585        NVAR BoxSize = root:Packages:AggregateModeling:BoxSize
1586        SVAR MultiParticleAttraction = root:Packages:AggregateModeling:MultiParticleAttraction //"Neutral;Attractive;Repulsive;Not allowed;"
1587        variable StickingProbability1, StickingProbabilityM1,StickingProbabilityM2, StickingProbabilityLoc
1588        if(StringMatch(MultiParticleAttraction, "Neutral"))
1589                StickingProbability1= StickingProbability
1590                StickingProbabilityM1= StickingProbability
1591                StickingProbabilityM2= StickingProbability
1592        elseif(StringMatch(MultiParticleAttraction, "Attractive"))
1593                StickingProbability1= StickingProbability
1594                StickingProbabilityM1= (StickingProbability+100)/2
1595                StickingProbabilityM2= (StickingProbability+300)/4
1596        elseif(StringMatch(MultiParticleAttraction, "Repulsive"))
1597                StickingProbability1 = StickingProbability
1598                StickingProbabilityM1 = (StickingProbability+10)/2
1599                StickingProbabilityM2 = (StickingProbability+30)/4
1600        elseif(StringMatch(MultiParticleAttraction, "Not allowed"))
1601                StickingProbability1 = StickingProbability
1602                StickingProbabilityM1 = 1
1603                StickingProbabilityM2 = 0
1604        else
1605                StickingProbability1 = StickingProbability
1606                StickingProbabilityM1 = StickingProbability
1607                StickingProbabilityM2 = StickingProbability
1608        endif
1609
1610
1611        make/Free/N=(NumParticles,3) CurSite
1612        make/Free/N=(NumParticles) tmpCol
1613        //new recording of particles...
1614        KillWaves/Z ListOfNeighbors, NumberOfNeighbors, ListOfStarts
1615        make/O/N=(NumParticles,26)/S ListOfNeighbors
1616        make/O/N=(NumParticles)/B/U NumberOfNeighbors
1617        ListOfNeighbors = Nan
1618        NumberOfNeighbors = 0
1619       
1620        variable chcnt=1,px,py,pz,aggct=1,cnt,con,stuck,dim,choice,wall,Rd=16,GL=16,farpoint,index=0, tmpVal
1621        variable i
1622       
1623       
1624       
1625       
1626        PauseUpdate
1627        Do
1628                // resize box based on size of Agg
1629                farpoint = max(MassFractalAggregate[aggct-1][0], MassFractalAggregate[aggct-1][1],MassFractalAggregate[aggct-1][2],farpoint)
1630                //For(i=0;i<3;i+=1)
1631                //      if(MassFractalAggregate[aggct-1][i]>farpoint)
1632                //              farpoint=MassFractalAggregate[aggct-1][i]
1633                //      endif
1634                //endfor
1635                GL=2*abs(farpoint)+10
1636                BoxSize=GL             
1637                // initialize particle on a random wall of the box
1638                wall=0;choice=0
1639                if(aggct<64)    // choose random wall on smaller box for low aggct
1640                        do
1641                                wall=abs(round(enoise(7)))
1642                        while(wall==0 || wall==7)
1643                        if(wall==1)
1644                                px=-1+Rd/2;py=round(enoise(Rd/2));pz=round(enoise(Rd/2))
1645                        endif
1646                        if(wall==2)
1647                                px=1-Rd/2;py=round(enoise(Rd/2));pz=round(enoise(Rd/2))
1648                        endif
1649                        if(wall==3)
1650                                px=round(enoise(Rd/2));py=-1+Rd/2;pz=round(enoise(Rd/2))
1651                        endif
1652                        if(wall==4)
1653                                px=round(enoise(Rd/2));py=1-Rd/2;pz=round(enoise(Rd/2))
1654                        endif
1655                        if(wall==5)
1656                                px=round(enoise(Rd/2));py=round(enoise(Rd/2));pz=-1+Rd/2
1657                        endif
1658                        if(wall==6)
1659                                px=round(enoise(Rd/2));py=round(enoise(Rd/2));pz=1-Rd/2
1660                        endif
1661                else            // choose random wall on normal box
1662                        do
1663                                wall=abs(round(enoise(7)))
1664                        while(wall==0 || wall==7)
1665                        if(wall==1)
1666                                px=-1+GL/2;py=round(enoise(GL/2));pz=round(enoise(GL/2))
1667                        endif
1668                        if(wall==2)
1669                                px=1-GL/2;py=round(enoise(GL/2));pz=round(enoise(GL/2))
1670                        endif
1671                        if(wall==3)
1672                                px=round(enoise(GL/2));py=-1+GL/2;pz=round(enoise(GL/2))
1673                        endif
1674                        if(wall==4)
1675                                px=round(enoise(GL/2));py=1-GL/2;pz=round(enoise(GL/2))
1676                        endif
1677                        if(wall==5)
1678                                px=round(enoise(GL/2));py=round(enoise(GL/2));pz=-1+GL/2
1679                        endif
1680                        if(wall==6)
1681                                px=round(enoise(GL/2));py=round(enoise(GL/2));pz=1-GL/2
1682                        endif
1683                endif
1684                // Move the particle until it a) hits the chain or b) leaves the box, if b) you need to have it reenter the box at a mirror position.
1685                do      //move 1 step in any direction
1686                        //choice=0
1687                        //do
1688                        //choice=abs(round(enoise(7)))
1689                        choice = floor(1 + mod(abs(enoise(100*6)),6))   
1690                        //while(choice==0 || choice==7)
1691                        if(choice==1)
1692                                px+=1
1693                        endif
1694                        if(choice==2)   
1695                                px-=1
1696                        endif
1697                        if(choice==3)   
1698                                py+=1
1699                        endif
1700                        if(choice==4)
1701                                py-=1
1702                        endif
1703                        if(choice==5)   
1704                                pz+=1
1705                        endif
1706                        if(choice==6)   
1707                                pz-=1
1708                        endif
1709                        //if you leave the box (a likely event) then go to the other side and reenter the box (mirror)
1710                        If(aggct<64)    // use smaller box for low number of particles
1711                                if(px>Rd/2)
1712                                        px=px-Rd
1713                                endif
1714                                if(px<-Rd/2)
1715                                        px=px+Rd
1716                                endif
1717                                if(py>Rd/2)
1718                                        py=py-Rd
1719                                endif
1720                                if(py<-Rd/2)
1721                                        py=py+Rd
1722                                endif
1723                                if(pz>Rd/2)
1724                                        pz=pz-Rd
1725                                endif
1726                                if(pz<-Rd/2)
1727                                        pz=pz+Rd
1728                                endif   
1729                        else            // use normal box for higher number of particles
1730                                if(px>GL/2)
1731                                        px=px-GL
1732                                endif
1733                                if(px<-GL/2)
1734                                        px=px+GL
1735                                endif
1736                                if(py>GL/2)
1737                                        py=py-GL
1738                                endif
1739                                if(py<-GL/2)
1740                                        py=py+GL
1741                                endif
1742                                if(pz>GL/2)
1743                                        pz=pz-GL
1744                                endif
1745                                if(pz<-GL/2)
1746                                        pz=pz+GL
1747                                endif
1748                        endif
1749                        cnt=0;con=0
1750                        // check how many neighboring sites are occupied
1751                        //this is by far the longest step in the whole procedure
1752                        //basically, we are looking for how many neighbors px,py,pz position has
1753                        //this is already much better than before, but it would really be nice to find better way of doing this. We are doing this A LOT. With every particle move, so it is done many, many times.
1754                        Multithread Distances[0,aggct] = ((px-MassFractalAggregate[p][0])^2 + (py-MassFractalAggregate[p][1])^2 + (pz-MassFractalAggregate[p][2])^2)           
1755                        //      Multithread helps, in my test case reduces time by ~60%
1756                        //this is slower...
1757                        //CurSite[][0]=px
1758                        //CurSite[][1]=py
1759                        //CurSite[][2]=pz
1760                        //MatrixOp/O/Free/NTHR=0 Distances = sumRows(powR((MassFractalAggregate-CurSite),2))
1761                        variable MaxDistance = 1.05*AllowedNearDistance                                                 //1 - in line nearest neighbor (1 step), 2 is two step nearest neighbor and 3 is nearest neighbor in any direction (3 step nearest neighbor).
1762                        Histogram/B={0.5,MaxDistance,2}/R=[0,aggct]/Dest=DistHist Distances             //histogram - bin[0] is from 0.5 - 3.1, max allowed distance^2 is 3
1763                        con = DistHist[0]                                                                                                               // this is number of nearest neighbors with distance below sqrt(3)
1764                        //another method, suggested by WM, but is slower, much slower ...
1765                        //                      CurSite[0][0]=px
1766                        //                      CurSite[0][1]=py
1767                        //                      CurSite[0][2]=pz
1768                        //                      MatrixOp/Free TestWv = catRows(CurSite,agg)
1769                        //                      FPClustering /DSO TestWv
1770                        //                      Wave M_DistanceMap
1771                        //                      //Edit/K=1 root:M_DistanceMap
1772                        //                      MatrixOp/Free Distances=row(M_DistanceMap,0)
1773                        //                      //Edit/K=1 root:Distances
1774                        //                      Histogram/B={0.1,1.8,2}/Dest=DistHist Distances
1775                        //                      //Edit/K=1 root:DistHist
1776                        //                      con = DistHist[0]
1777                        // end of neighbor counting.
1778                        //choice=0
1779                        if(con>=1)      // particle can StickingProbability if there is at least 1 neighbor
1780                                if(con>=3)
1781                                        StickingProbabilityLoc=StickingProbabilityM2
1782                                elseif(con==2)
1783                                        StickingProbabilityLoc=StickingProbabilityM1
1784                                else
1785                                        StickingProbabilityLoc=StickingProbability1
1786                                endif
1787
1788                                do
1789                                        // apply StickingProbabilitying probability between 1% and 100%
1790                                        choice=floor(1 + mod(abs(enoise(100*99)),100))          //generates random integer from 1 to 99
1791                                        choice-=1
1792                                        AttemptValue+=1
1793                                        if(choice<=StickingProbabilityLoc)
1794                                                stuck=1
1795                                                con=0
1796                                        else
1797                                                con-=1
1798                                                stuck=0
1799                                        endif
1800                                while(con>=1)
1801                        //keep moving if alone or rejected by StickingProbabilitying probability
1802                        else
1803                                stuck=0
1804                        endif
1805                        //if the particle StickingProbabilitys, add it to the aggregate
1806                        variable steps=trunc(DegreeOfAggregation/10)
1807                        steps = max(steps,100)
1808                        if(stuck==1 && IR3A_IsPXYZNOTinList3DWave(MassFractalAggregate,px,py,pz, aggct))        //added here to make sure we do not accidentally add existing particle.
1809                                if(mod(aggct,steps)<1) ///round(DegreeOfAggregation/50))==aggct/round(DegreeOfAggregation/50))
1810                                        Print time()+"  Added "+num2str(aggct)+" particles to the aggregate  "  //takes needless time..
1811                                endif
1812                                MassFractalAggregate[aggct][0]=px
1813                                MassFractalAggregate[aggct][1]=py
1814                                MassFractalAggregate[aggct][2]=pz
1815                                //record here particle and its neighbors.
1816                                //who is neighbor? that who has Distances[numparticle]<MaxDistance
1817                                //ListOfNeighbors[aggct] needs to get ^^ added to it.
1818                                //NumberOfNeighbors
1819                                IR3A_AddToNeighborList(ListOfNeighbors,NumberOfNeighbors, Distances, MaxDistance, aggct)
1820                                aggct+=1
1821                        endif
1822                while(stuck==0)
1823                stuck=0 //reset stuck flag
1824        While(aggct<DegreeOfAggregation)        // stop aggregate growth when there are DegreeOfAggregation particles in Agg
1825        Print time()+"  Created Aggregate with "+num2str(aggct)+" particles in it"      //takes needless time..
1826        setDataFOlder OldDf
1827
1828End
1829//******************************************************************************************************************************************************
1830Function IR3A_AddToNeighborList(ListOfNeighbors,NumberOfNeighbors, Distances, MaxDistance, aggct)
1831        wave ListOfNeighbors, NumberOfNeighbors, Distances
1832        variable MaxDistance, aggct
1833       
1834        variable i
1835        For(i=0;i<aggct;i+=1)
1836                if(Distances[i]<MaxDistance)    //this is neighbor!     
1837                        //i is now old particle as neighbor
1838                        //aggct is new particle
1839                        ListOfNeighbors[aggct][NumberOfNeighbors[aggct]]=i
1840                        NumberOfNeighbors[aggct]+=1
1841                        ListOfNeighbors[i][NumberOfNeighbors[i]]=aggct
1842                        NumberOfNeighbors[i]+=1
1843                endif
1844        endfor
1845       
1846end
1847//******************************************************************************************************************************************************
1848
1849//static
1850// Function IR3A_MakeAgg(DegreeOfAggregation,MassFractalAggregate,StickingProbability, AllowedNearDistance)
1851//      variable DegreeOfAggregation,StickingProbability, AllowedNearDistance
1852//      wave MassFractalAggregate
1853//     
1854//      DFref oldDf= GetDataFolderDFR()
1855//
1856//      SetDataFolder root:Packages:AggregateModeling
1857//      Wave Distances
1858//      make/Free/N=(dimsize(MassFractalAggregate,0),3) CurSite
1859//      make/Free/N=(dimsize(MassFractalAggregate,0)) tmpCol
1860//      variable chcnt=1,px,py,pz,aggct=1,cnt,con,stuck,dim,choice,wall,Rd=16,GL=16,farpoint,index=0, tmpVal
1861//      NVAR AttemptValue
1862//      AttemptValue=1
1863//      NVAR BoxSize
1864//      PauseUpdate
1865//      Do
1866//              // resize box based on size of Agg
1867//              index=0
1868//              do
1869//                      if(MassFractalAggregate[aggct-1][index]>farpoint)
1870//                              farpoint=MassFractalAggregate[aggct-1][index]
1871//                      endif
1872//                      index+=1
1873//              while(index<3)
1874//              GL=2*abs(farpoint)+10
1875//              BoxSize=GL             
1876//              // initialize particle on a random wall of the box
1877//              wall=0;choice=0
1878//              if(aggct<64)    // choose random wall on smaller box for low aggct
1879//                      do
1880//                              wall=abs(round(enoise(7)))
1881//                      while(wall==0 || wall==7)
1882//                      if(wall==1)
1883//                              px=-1+Rd/2;py=round(enoise(Rd/2));pz=round(enoise(Rd/2))
1884//                      endif
1885//                      if(wall==2)
1886//                              px=1-Rd/2;py=round(enoise(Rd/2));pz=round(enoise(Rd/2))
1887//                      endif
1888//                      if(wall==3)
1889//                              px=round(enoise(Rd/2));py=-1+Rd/2;pz=round(enoise(Rd/2))
1890//                      endif
1891//                      if(wall==4)
1892//                              px=round(enoise(Rd/2));py=1-Rd/2;pz=round(enoise(Rd/2))
1893//                      endif
1894//                      if(wall==5)
1895//                              px=round(enoise(Rd/2));py=round(enoise(Rd/2));pz=-1+Rd/2
1896//                      endif
1897//                      if(wall==6)
1898//                              px=round(enoise(Rd/2));py=round(enoise(Rd/2));pz=1-Rd/2
1899//                      endif
1900//              else            // choose random wall on normal box
1901//                      do
1902//                              wall=abs(round(enoise(7)))
1903//                      while(wall==0 || wall==7)
1904//                      if(wall==1)
1905//                              px=-1+GL/2;py=round(enoise(GL/2));pz=round(enoise(GL/2))
1906//                      endif
1907//                      if(wall==2)
1908//                              px=1-GL/2;py=round(enoise(GL/2));pz=round(enoise(GL/2))
1909//                      endif
1910//                      if(wall==3)
1911//                              px=round(enoise(GL/2));py=-1+GL/2;pz=round(enoise(GL/2))
1912//                      endif
1913//                      if(wall==4)
1914//                              px=round(enoise(GL/2));py=1-GL/2;pz=round(enoise(GL/2))
1915//                      endif
1916//                      if(wall==5)
1917//                              px=round(enoise(GL/2));py=round(enoise(GL/2));pz=-1+GL/2
1918//                      endif
1919//                      if(wall==6)
1920//                              px=round(enoise(GL/2));py=round(enoise(GL/2));pz=1-GL/2
1921//                      endif
1922//              endif
1923//              // Move the particle until it a) hits the chain or b) leaves the box, if b) you need to have it reenter the box at a mirror position.
1924//              do      //move 1 step in any direction
1925//                      //choice=0
1926//                      //do
1927//                      //choice=abs(round(enoise(7)))
1928//                      choice = floor(1 + mod(abs(enoise(100*6)),6))   
1929//                      //while(choice==0 || choice==7)
1930//                      if(choice==1)
1931//                              px+=1
1932//                      endif
1933//                      if(choice==2)   
1934//                              px-=1
1935//                      endif
1936//                      if(choice==3)   
1937//                              py+=1
1938//                      endif
1939//                      if(choice==4)
1940//                              py-=1
1941//                      endif
1942//                      if(choice==5)   
1943//                              pz+=1
1944//                      endif
1945//                      if(choice==6)   
1946//                              pz-=1
1947//                      endif
1948//                      //if you leave the box (a likely event) then go to the other side and reenter the box (mirror)
1949//                      If(aggct<64)    // use smaller box for low number of particles
1950//                              if(px>Rd/2)
1951//                                      px=px-Rd
1952//                              endif
1953//                              if(px<-Rd/2)
1954//                                      px=px+Rd
1955//                              endif
1956//                              if(py>Rd/2)
1957//                                      py=py-Rd
1958//                              endif
1959//                              if(py<-Rd/2)
1960//                                      py=py+Rd
1961//                              endif
1962//                              if(pz>Rd/2)
1963//                                      pz=pz-Rd
1964//                              endif
1965//                              if(pz<-Rd/2)
1966//                                      pz=pz+Rd
1967//                              endif   
1968//                      else            // use normal box for higher number of particles
1969//                              if(px>GL/2)
1970//                                      px=px-GL
1971//                              endif
1972//                              if(px<-GL/2)
1973//                                      px=px+GL
1974//                              endif
1975//                              if(py>GL/2)
1976//                                      py=py-GL
1977//                              endif
1978//                              if(py<-GL/2)
1979//                                      py=py+GL
1980//                              endif
1981//                              if(pz>GL/2)
1982//                                      pz=pz-GL
1983//                              endif
1984//                              if(pz<-GL/2)
1985//                                      pz=pz+GL
1986//                              endif
1987//                      endif
1988//                      cnt=0;con=0
1989//                      // check how many neighboring sites are occupied
1990//                      //this is by far the longest step in the whole procedure
1991//                      //basically, we are looking for how many neighbors px,py,pz position has
1992//                      //this is already much better than before, but it would really be nice to find better way of doing this. We are doing this A LOT. With every particle move, so it is done many, many times.
1993//                      Multithread Distances[0,aggct] = ((px-MassFractalAggregate[p][0])^2 + (py-MassFractalAggregate[p][1])^2 + (pz-MassFractalAggregate[p][2])^2)            //      Multithread helps, in my test case reduces time by ~60%
1994//                      //this is slower...
1995//                      //CurSite[][0]=px
1996//                      //CurSite[][1]=py
1997//                      //CurSite[][2]=pz
1998//                      //MatrixOp/O/Free/NTHR=0 Distances = sumRows(powR((MassFractalAggregate-CurSite),2))
1999//                      variable MaxDistance = 1.05*AllowedNearDistance                 //1 - in line nearest neighbor (1 step), 2 is two step nearest neighbor and 3 is nearest neighbor in any direction (3 step nearest neighbor).
2000//                      Histogram/B={0.5,MaxDistance,2}/R=[0,aggct]/Dest=DistHist Distances                     //histogram - bin[0] is from 0.5 - 3.1, max allowed distance^2 is 3
2001//                      con = DistHist[0]                                                                                                                                       // this is number of nearest neighbors with distance below sqrt(3)
2002//                      //another method, suggested by WM, but is slower, much slower ...
2003//                      //                      CurSite[0][0]=px
2004//                      //                      CurSite[0][1]=py
2005//                      //                      CurSite[0][2]=pz
2006//                      //                      MatrixOp/Free TestWv = catRows(CurSite,agg)
2007//                      //                      FPClustering /DSO TestWv
2008//                      //                      Wave M_DistanceMap
2009//                      //                      //Edit/K=1 root:M_DistanceMap
2010//                      //                      MatrixOp/Free Distances=row(M_DistanceMap,0)
2011//                      //                      //Edit/K=1 root:Distances
2012//                      //                      Histogram/B={0.1,1.8,2}/Dest=DistHist Distances
2013//                      //                      //Edit/K=1 root:DistHist
2014//                      //                      con = DistHist[0]
2015//                      // end of neighbor counting.
2016//                      //choice=0
2017//                      if(con>=1)      // particle can StickingProbability if there is at least 1 neighbor
2018//                              do
2019//                                      // apply StickingProbabilitying probability between 1% and 100%
2020//                                      choice=floor(1 + mod(abs(enoise(100*99)),100))          //generates random integer from 1 to 99
2021//                                      choice-=1
2022//                                      AttemptValue+=1
2023//                                      if(choice<=StickingProbability)
2024//                                              stuck=1
2025//                                              con=0
2026//                                      else
2027//                                              con-=1
2028//                                              stuck=0
2029//                                      endif
2030//                              while(con>=1)
2031//                      //keep moving if alone or rejected by StickingProbabilitying probability
2032//                      else
2033//                              stuck=0
2034//                      endif
2035//                      //if the particle StickingProbabilitys, add it to the aggregate
2036//                      variable steps=trunc(DegreeOfAggregation/10)
2037//                      steps = max(steps,100)
2038//                      if(stuck==1 && IR3A_IsPXYZNOTinList3DWave(MassFractalAggregate,px,py,pz, aggct))        //added here to make sure we do nto accidentally add existing particle.
2039//                              if(mod(aggct,steps)<1) ///round(DegreeOfAggregation/50))==aggct/round(DegreeOfAggregation/50))
2040//                                      Print time()+"  Added "+num2str(aggct)+" particles to the aggregate  "  //takes needless time..
2041//                              endif
2042//                              MassFractalAggregate[aggct][0]=px
2043//                              MassFractalAggregate[aggct][1]=py
2044//                              MassFractalAggregate[aggct][2]=pz
2045//                              aggct+=1
2046//                      endif
2047//              while(stuck==0)
2048//              stuck=0 //reset stuck flag
2049//      While(aggct<DegreeOfAggregation)        // stop aggregate growth when there are DegreeOfAggregation particles in Agg
2050//      Print time()+"  Created Aggregate with "+num2str(aggct)+" particles in it"      //takes needless time..
2051//      setDataFOlder OldDf
2052//
2053//End
2054//******************************************************************************************************************************************************
2055//******************************************************************************************************************************************************
2056//static
2057Function IR3A_IsPXYZNOTinList3DWave(A3DWaveList,px,py,pz, MaxPoints)
2058        wave A3DWaveList
2059        variable px,py,pz, MaxPoints
2060        variable i
2061        FOr(i=0;i<MaxPoints;i+=1)
2062                if((abs(A3DWaveList[i][0]-px)+abs(A3DWaveList[i][1]-py)+abs(A3DWaveList[i][2]-pz))<0.2)
2063                        return 0
2064                endif
2065        endfor
2066        return 1
2067       
2068end
2069//******************************************************************************************************************************************************
2070//******************************************************************************************************************************************************
2071//******************************************************************************************************************************************************
2072//******************************************************************************************************************************************************
2073//******************************************************************************************************************************************************
2074
2075//static Function IR3A_Ends(MassFractalAggregate, breakOnFail)
2076//      wave MassFractalAggregate
2077//      variable breakOnFail
2078//
2079//      DFref oldDf= GetDataFolderDFR()
2080//      SetDataFolder root:Packages:AggregateModeling
2081//      Wave/Z ListOfNeighbors = root:Packages:AggregateModeling:ListOfNeighbors
2082//      //this is list of each particle neighbors. Any one which has only one neigbor, is end point.
2083//      //we want to create wave containing row numbers for rows, where Column[2] is =0 
2084//      if(WaveExists(ListOfNeighbors))
2085//              make/O/N=(DimSize(ListOfNeighbors, 0)) endpoints
2086//              variable i
2087//              endpoints = numtype(ListOfNeighbors[p][1])==2 ? p : NaN
2088//              Extract /O endpoints, endpoints, numtype(endpoints)==0
2089//      else
2090//              //this is old method, which for soem reason even does nto work as far as I can say...
2091//              NVAR DegreeOfAggregation
2092//              Wave endpoints
2093//     
2094//              variable cnt=0,ncnt=0,con=0,endcnt=0,dim=0
2095//              Make/n=(26,3)/Free nghbr
2096//              Wave/Z nghbrOfsetList
2097//              if(!WaveExists( nghbrOfsetList))
2098//                      IR3A_MakeNBROffsetList()
2099//              endif
2100//              do
2101//                      // define neighbors for each point in agg
2102//                      nghbr = nghbrOfsetList[p][q] + MassFractalAggregate[cnt][q]
2103//                      ncnt=0;con=0
2104//                      do
2105//                              dim=0
2106//                              do
2107//                                      if(nghbr[dim][0]==MassFractalAggregate[ncnt][0]&&nghbr[dim][1]==MassFractalAggregate[ncnt][1]&&nghbr[dim][2]==MassFractalAggregate[ncnt][2])
2108//                                              con+=1
2109//                                      endif
2110//                                      dim+=1
2111//                              while(dim<26)
2112//                              ncnt+=1
2113//                      while(ncnt<DegreeOfAggregation)
2114//                      // it's an endpoint if there is exactly 1 neighboring point
2115//                      // record position in x, y, z and then record index in agg
2116//                      if(con==1)
2117//                              endpoints[endcnt][0]=MassFractalAggregate[cnt][0];endpoints[endcnt][1]=MassFractalAggregate[cnt][1];endpoints[endcnt][2]=MassFractalAggregate[cnt][2];endpoints[endcnt][3]=cnt
2118//                              endcnt+=1
2119//                      endif
2120//                      cnt+=1
2121//                      //Print cnt
2122//              while(cnt<DegreeOfAggregation)
2123//              //remove extra rows from endpoints wave
2124//              DeletePoints endcnt,DegreeOfAggregation, endpoints
2125//      endif
2126//      variable Failed = 0
2127//      if(dimsize(endpoints,0)<5)
2128//              Failed = 1
2129//              if(breakOnFail)
2130//                      DoALert/T="FailedGrowth" 0, "Not enough endpoints found, too compact particle"
2131//              endif
2132//      endif
2133//      setDataFOlder OldDf
2134//      Print time()+"  Finished running Find Ends"     //takes needless time..
2135//      return Failed
2136//End
2137//******************************************************************************************************************************************************
2138//******************************************************************************************************************************************************
2139//******************************************************************************************************************************************************
2140//******************************************************************************************************************************************************
2141//******************************************************************************************************************************************************
2142// 3-16-2021 verification - yields same numbers as original code.
2143//static Function IR3A_Reted(endpoints)
2144//      wave endpoints
2145//      // calculate longest end-to-end distance for each combination of endpoints and its square for weight-averaged end-to-end distance       
2146//      DFref oldDf= GetDataFolderDFR()
2147//
2148//      SetDataFolder root:Packages:AggregateModeling
2149//      NVAR DegreeOfAggregation=root:Packages:AggregateModeling:DegreeOfAggregation
2150//      NVAR RValue=root:Packages:AggregateModeling:RValue
2151//      NVAR dfValue=root:Packages:AggregateModeling:dfValue
2152//      NVAR RgPrimary=root:Packages:AggregateModeling:RgPrimary
2153//      NVAR RxRgPrimaryValue=root:Packages:AggregateModeling:RxRgPrimaryValue
2154//
2155//      variable endnum=DimSize(endpoints,0), numcomb=binomial(endnum,2), cnt=0,endadd=0,ecnt=1,REnd=0,Rsum=0,Retend=0,RAve=0,rem=endnum
2156//      Make/n=(numcomb,1)/o enddist=0
2157//      Make/n=(endnum-3,7)/o Rlarge=0
2158//      do      // determine end-to-end distance between all endpoints
2159//              do
2160//                      enddist[endadd]=sqrt((endpoints[cnt][0]-endpoints[ecnt][0])^2+(endpoints[cnt][1]-endpoints[ecnt][1])^2+(endpoints[cnt][2]-endpoints[ecnt][2])^2)
2161//                      if(endnum-cnt>3)       
2162//                              if(enddist[endadd]>Rlarge[cnt][6])
2163//                                      Rlarge[cnt][0]=endpoints[cnt][0];Rlarge[cnt][1]=endpoints[cnt][1];Rlarge[cnt][2]=endpoints[cnt][2]
2164//                                      Rlarge[cnt][3]=endpoints[ecnt][0];Rlarge[cnt][4]=endpoints[ecnt][1];Rlarge[cnt][5]=endpoints[ecnt][2]
2165//                                      Rlarge[cnt][6]=enddist[endadd]
2166//                              endif
2167//                      endif
2168//                      ecnt+=1
2169//                      endadd+=1
2170//              while(ecnt<endnum)
2171//              cnt+=1
2172//              ecnt=cnt+1
2173//      while(ecnt<endnum)
2174//      cnt=0
2175//      do      // calculate longest end-to-end distance for each combination of endpoints and its square for weight-averaged end-to-end distance
2176//              REnd+=(Rlarge[cnt][6])*(Rlarge[cnt][6])
2177//              RSum+=(Rlarge[cnt][6])
2178//              cnt+=1
2179//      while(cnt<endnum-3)
2180//      RAve=RSum/cnt
2181//      REnd/=RSum
2182//      cnt=0
2183//      // Print and record R, df
2184//     
2185//      Print "R = "+num2str(REnd)
2186//      Print "df= "+num2str(log(DegreeOfAggregation)/log(REnd))
2187//      RValue=Rend
2188//      dfValue=log(DegreeOfAggregation)/log(REnd)
2189//      RxRgPrimaryValue = REnd*RgPrimary
2190//      setDataFOlder OldDf
2191//End
2192//******************************************************************************************************************************************************
2193//******************************************************************************************************************************************************
2194//******************************************************************************************************************************************************
2195//******************************************************************************************************************************************************
2196//******************************************************************************************************************************************************
2197//
2198//static Function IR3A_Path(NumberOfTestPaths)
2199//      variable NumberOfTestPaths
2200//
2201//      DFref oldDf= GetDataFolderDFR()
2202//
2203//      SetDataFolder root:Packages:AggregateModeling
2204//      print time()+"  Started parameters evaluation, Calculating...   this takes longest time for these functions "
2205//      wave MassFractalAggregate
2206//      Wave endpoints
2207//      NVAR DegreeOfAggregation
2208//      NVAR AttemptValue
2209//      variable minPathLengthAccepted
2210//      minPathLengthAccepted = floor(min(DegreeOfAggregation^(1/2), 20))               //minimum length of path accpeted as evaluation path, for large particles original value is crazily large.
2211//      Make/n=(26,3)/Free nghbr
2212//      Wave/Z nghbrOfsetList
2213//      if(!WaveExists(nghbrOfsetList))
2214//              IR3A_MakeNBROffsetList()
2215//      endif
2216//      Make /n=(DegreeOfAggregation,26)/O NeighborList=nan
2217//      Make/n=(DegreeOfAggregation,1)/Free Pathing=nan
2218//      Make/n=(NumberOfTestPaths,7)/Free Paths=nan
2219//      variable endnum=DimSize(endpoints,0), choice=0, cnt=0,ncnt=0,nnum=0,dim=0,con=0,pcnt=0,nans=0,flag=0,ecnt=0,highp,p1,p2,mom1,mom2,startpoint
2220//      Make/n=(endnum,1)/O LongPath=0
2221//      variable i
2222//      MatrixOp/Free EndPointsDistances=col(endpoints,3)
2223//      variable startTime=ticks
2224//      For(cnt=0;cnt<DegreeOfAggregation;cnt+=1)
2225//              nghbr = nghbrOfsetList[p][q] + MassFractalAggregate[cnt][q]                                     //this generates list of 26 neighbor positions, if any, around the agg[cnt] particle
2226//              nnum=0
2227//              For(ncnt=0;ncnt<DegreeOfAggregation;ncnt+=1)
2228//                      For(dim=0;dim<26;dim+=1)
2229//                              if(nghbr[dim][0]==MassFractalAggregate[ncnt][0]&&nghbr[dim][1]==MassFractalAggregate[ncnt][1]&&nghbr[dim][2]==MassFractalAggregate[ncnt][2])
2230//                                      NeighborList[cnt][nnum]=ncnt                    //cnt is external loop, cnt = 0 ...DegreeOfAggregation
2231//                                      nnum+=1
2232//                              endif
2233//                      endfor
2234//              endfor
2235//      endfor
2236//      cnt=0
2237//      do      // snake through paths
2238//              pcnt=0
2239//              startpoint=floor(mod(abs(enoise(100*(endnum-1))),endnum))               //generates random integer from 0 to endnum
2240//              // record starting position and starting index in agg
2241//              Paths[cnt][0]=endpoints[startpoint][0]
2242//              Paths[cnt][1]=endpoints[startpoint][1]
2243//              Paths[cnt][2]=endpoints[startpoint][2]
2244//              Pathing[0]=endpoints[startpoint][3]
2245//              choice=0;con=0
2246//              do      //pick a random neighbor from list
2247//                      ncnt=0
2248//                      nans=0
2249//                      flag=0
2250//                      //this is faster...     
2251//                      For(i=0;i<8;i+=1)
2252//                              if(numtype(NeighborList[Pathing[pcnt]][i])==2)
2253//                                      nans=8-i
2254//                                      break
2255//                              endif
2256//                      endfor
2257//                      variable target= 8-nans-1
2258//                      choice=floor(mod(abs(enoise(100*target)),target+1))             //generates random number from non NaNs entries...             
2259//                      ncnt=0
2260//                      do      // check to see if value already in path
2261//                              if(NeighborList[Pathing[pcnt]][choice]==Pathing[ncnt])
2262//                                      flag+=1
2263//                              endif
2264//                              ncnt+=1
2265//                      while(ncnt<pcnt)
2266//                      if(flag>=2)     // is there a loop in the aggregate?
2267//                              con=1
2268//                      else            // add normally
2269//                              pcnt+=1                         //this pcnt is sometimes too large, not sure how to avoid it...
2270//                              if(pcnt>=DegreeOfAggregation)
2271//                                      Abort "Wrong value in IR3A_Path, try again."
2272//                              endif
2273//                              Pathing[pcnt]=NeighborList[Pathing[pcnt-1]][choice]
2274//                      endif
2275//                      ncnt=1
2276//                      ecnt=0
2277//                      // this is about 2x faster...  EndPointsDistances was created above just once, it does not change.
2278//                      FindValue /V=(Pathing[pcnt])/T=0.1 EndPointsDistances
2279//                      if(V_value>-0.1)
2280//                              con=1
2281//                      endif
2282//                      // (2) done...                 
2283//              while(con!=1)
2284//              if(pcnt>minPathLengthAccepted)                          // only interested in longer paths that span aggregate; throws away short paths
2285//                      Paths[cnt][3]=MassFractalAggregate[Pathing[pcnt]][0]
2286//                      Paths[cnt][4]=MassFractalAggregate[Pathing[pcnt]][1]
2287//                      Paths[cnt][5]=MassFractalAggregate[Pathing[pcnt]][2]
2288//                      Paths[cnt][6]=pcnt+1
2289//                      if(Paths[cnt][6]>LongPath[startpoint])
2290//                              LongPath[startpoint]=Paths[cnt][6]     
2291//                      endif
2292//                      cnt+=1
2293//                      //print "Evaluated path length of "+num2str(pcnt)
2294//              else
2295//                      //print "Discarded path length of "+num2str(pcnt)
2296//                      continue
2297//              endif
2298//              if(mod(cnt,500)==0)
2299//                      Print time()+"  Working... Evaluated "+num2str(cnt)+" Paths through the system" //takes needless time..
2300//              endif
2301//      while(cnt<NumberOfTestPaths)
2302//      ncnt=0;highp=0;p1=0;p2=0
2303//      do      // determine weight-averaged percolation pathway
2304//              if(LongPath[ncnt]>highp)
2305//                      highp=LongPAth[ncnt]
2306//              endif
2307//              p1+=(LongPath[ncnt])
2308//              p2+=(LongPath[ncnt])*(LongPath[ncnt])
2309//              ncnt+=1
2310//      while(ncnt<endnum)
2311//      mom2=round(p2/p1);mom1=round(p1/endnum)
2312//      // print results
2313//      NVAR dfValue
2314//      NVAR RValue
2315//      NVAR pValue
2316//      NVAR sValue
2317//      NVAR cValue
2318//      NVAR dminValue
2319//      NVAR TrueStickingProbability
2320//      NVAR StickingProbability=root:Packages:AggregateModeling:StickingProbability
2321//      NVAR RgPrimary=root:Packages:AggregateModeling:RgPrimary
2322//      NVAR AllowedNearDistance=root:Packages:AggregateModeling:AllowedNearDistance
2323//      pValue = mom2
2324//      cValue=ln(DegreeOfAggregation)/ln(pValue)
2325//      dminValue=dfValue/cValue
2326//      sValue=round(exp(ln(DegreeOfAggregation)/dminValue))
2327//      TrueStickingProbability = 100*DegreeOfAggregation/AttemptValue
2328//      NVAR RxRgPrimaryValue=root:Packages:AggregateModeling:RxRgPrimaryValue
2329//      variable PrimaryDiameter = 2*sqrt(5/3)*RgPrimary
2330//      RxRgPrimaryValue = RValue*PrimaryDiameter
2331//      Print "R [primary particles]= "+num2str(RValue)
2332//      Print "R [Angstroms] = "+num2str(RValue*PrimaryDiameter)
2333//      Print "z = "+num2str(DegreeOfAggregation)
2334//      Print "p = "+num2str(pValue)
2335//      Print "s = "+num2str(sValue)
2336//      Print "df = "+num2str(dfValue)
2337//      Print "dmin = "+num2str(dminValue)
2338//      Print "c = "+num2str(cValue)
2339//      Print "True Sticking Probability = "+num2str(100*DegreeOfAggregation/AttemptValue)+"%"
2340//      //appned note to MassFractalAggregate
2341//      string NoteText
2342//      NoteText="Mass Fractal Aggregate created="+date()+", "+time()+";z="+num2str(DegreeOfAggregation)+";StickingProbability="+num2str(StickingProbability)+";StickingMethod="+num2str(AllowedNearDistance)+";R="+num2str(RValue)+";Rprimary="+num2str(RgPrimary)+";p="+num2str(pValue)
2343//      NoteText+=";RxRgPrimaryValue="+num2str(RValue*PrimaryDiameter)+";s="+num2str(sValue)+";df="+num2str(IN2G_roundSignificant(dfValue,4))+";dmin="+num2str(IN2G_roundSignificant(dminValue,4))+";c="+num2str(IN2G_roundSignificant(cValue,4))+";True Sticking Probability="+num2str(100*DegreeOfAggregation/AttemptValue)+";"
2344//      Note MassFractalAggregate, NoteText
2345//      setDataFolder OldDf
2346//End
2347
2348//******************************************************************************************************************************************************
2349//******************************************************************************************************************************************************
2350//******************************************************************************************************************************************************
2351//******************************************************************************************************************************************************
2352//******************************************************************************************************************************************************
2353//******************************************************************************************************************************************************
2354//******************************************************************************************************************************************************
2355
2356
2357static Function IR3A_GetResults()
2358        DFref oldDf= GetDataFolderDFR()
2359
2360        SetDataFolder root:Packages:AggregateModeling
2361        NVAR RgPrimary=root:Packages:AggregateModeling:RgPrimary
2362        NVAR RxRgPrimaryValue=root:Packages:AggregateModeling:RxRgPrimaryValue
2363        NVAR DegreeOfAggregation
2364        NVAR RValue
2365        NVAR pValue
2366        NVAR dfValue
2367        NVAR dminValue
2368        NVAR cValue
2369        NVAR sValue
2370        NVAR AttemptValue
2371        NVAR AllowedNearDistance
2372        cValue=ln(DegreeOfAggregation)/ln(pValue)
2373        dminValue=dfValue/cValue
2374        sValue=round(exp(ln(DegreeOfAggregation)/dminValue))
2375        Print "R [A]= "+num2str(RxRgPrimaryValue)
2376        Print "df = "+num2str(dfValue)
2377        Print "z = "+num2str(DegreeOfAggregation)
2378        Print "p= "+num2str(pValue)
2379        Print "s = "+num2str(sValue)
2380        Print "dmin = "+num2str(dminValue)
2381        Print "c = "+num2str(cValue)
2382        Print "True Sticking Probability = "+num2str(100*DegreeOfAggregation/AttemptValue)+"%"
2383        print "Sticking method = "+num2str(AllowedNearDistance)
2384        setDataFOlder OldDf
2385End
2386
2387//******************************************************************************************************************************************************
2388//******************************************************************************************************************************************************
2389//******************************************************************************************************************************************************
2390//******************************************************************************************************************************************************
2391//******************************************************************************************************************************************************
2392//******************************************************************************************************************************************************
2393//******************************************************************************************************************************************************
2394
2395
2396static Function IR3A_GizmoViewScatterPlot(ScatterPlotWave) : GizmoPlot
2397        wave ScatterPlotWave
2398
2399        variable mnd = IR3A_GetAggMaxSize(ScatterPlotWave)
2400        variable boxSize=mnd+2          //make the box larger...
2401        DoWIndow MassFractalAggregateView
2402        if(V_Flag)
2403                DoWIndow/F MassFractalAggregateView
2404                ModifyGizmo setOuterBox={-1*boxSize,boxSize,-1*boxSize,boxSize,-1*boxSize,boxSize}
2405                ModifyGizmo scalingOption=0
2406                ModifyGizmo modifyObject=Particle,objectType=Sphere,property={radius,1/(4+mnd)}
2407        else
2408                PauseUpdate             // building window...
2409                // Building Gizmo 7 window...
2410                NewGizmo/K=1/T="Mass Fractal Aggregate View"/W=(35,45,550,505)
2411                DoWIndow/C MassFractalAggregateView
2412                ModifyGizmo startRecMacro=700
2413                AppendToGizmo Scatter=ScatterPlotWave,name=scatter0
2414                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ scatterColorType,0}
2415                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ markerType,0}
2416                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ sizeType,0}
2417                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ rotationType,0}
2418                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ Shape,2}
2419                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ size,1}
2420                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ color,1,0,0,1}
2421                //sphere as object definition:
2422                AppendToGizmo sphere={1/mnd,25,25},name=Particle
2423                ModifyGizmo modifyObject=Particle,objectType=Sphere,property={colorType,1}
2424                ModifyGizmo modifyObject=Particle,objectType=Sphere,property={color,0.000015,0.600000,0.304250,1.000000}
2425                AppendToGizmo attribute diffuse={0.5,0.5,0.5,1,1032},name=diffuse0
2426                ModifyGizmo attributeType=diffuse,modifyAttribute={diffuse0,0.733333,0.733333,0.733333,1,1032}
2427                ModifyGizmo modifyObject=Particle,objectType=Sphere,property={radius,1/(4+mnd)}
2428                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ Shape,7}
2429                ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ objectName,Particle}
2430                AppendToGizmo Axes=boxAxes,name=axes0
2431                ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
2432                ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
2433                ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,3}
2434                ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,3}
2435                ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,3}
2436                ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
2437                AppendToGizmo light=Directional,name=LightFor3dView
2438                ModifyGizmo modifyObject=LightFor3dView,objectType=light,property={ position,0.000000,0.000000,-1.000000,0.000000}
2439                ModifyGizmo modifyObject=LightFor3dView,objectType=light,property={ direction,0.000000,0.000000,-1.000000}
2440                ModifyGizmo modifyObject=LightFor3dView,objectType=light,property={ ambient,0.400000,0.400000,0.400000,1.000000}
2441                ModifyGizmo modifyObject=LightFor3dView,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
2442                ModifyGizmo modifyObject=LightFor3dView,objectType=light,property={ diffuse,0.933333,0.933333,0.933333,1.000000}
2443                ModifyGizmo modifyObject=LightFor3dView,objectType=light,property={ position,-0.6392,0.7354,0.2250,0.0000}
2444                ModifyGizmo modifyObject=LightFor3dView,objectType=light,property={ direction,-0.6392,0.7354,0.2250}
2445                ModifyGizmo setDisplayList=0, object=LightFor3dView
2446                ModifyGizmo setDisplayList=1, object=scatter0
2447                ModifyGizmo setDisplayList=2, object=axes0
2448                ModifyGizmo currentGroupObject=""
2449                //ModifyGizmo SETQUATERNION={0.535993,-0.191531,-0.283415,0.771818}
2450                //now scale this thing...
2451                ModifyGizmo setOuterBox={-1*boxSize,boxSize,-1*boxSize,boxSize,-1*boxSize,boxSize}
2452                ModifyGizmo scalingOption=0
2453                Modifygizmo enhance
2454                //give user tools to work with
2455                //ModifyGizmo showInfo
2456                //ModifyGizmo infoWindow={350,550,822,320}
2457                ModifyGizmo resumeUpdates
2458                ModifyGizmo endRecMacro
2459                ModifyGizmo SETQUATERNION={-0.041312,-0.884834,-0.102589,0.452588}
2460        endif
2461        DoWindow FractalAggregatePanel
2462        if(V_Flag)
2463                AutoPositionWindow /M=0/R=FractalAggregatePanel MassFractalAggregateView
2464        endif
2465        DoWindow MassFractalAggDataPlot
2466        if(V_Flag)
2467                AutoPositionWindow /M=1/R=MassFractalAggDataPlot MassFractalAggregateView
2468        endif
2469
2470EndMacro
2471//******************************************************************************************************************************************************
2472//******************************************************************************************************************************************************
2473//******************************************************************************************************************************************************
2474//******************************************************************************************************************************************************
2475//******************************************************************************************************************************************************
2476//******************************************************************************************************************************************************
2477//******************************************************************************************************************************************************
2478
2479
2480static Function IR3A_GetAggMaxSize(MassFractalAggregate)
2481        wave MassFractalAggregate
2482       
2483        WaveStats/Q MassFractalAggregate
2484
2485        variable MaxNeeded=max(V_max, abs(V_min) )
2486       
2487        return MaxNeeded
2488end
2489
2490//******************************************************************************************************************************************************
2491//******************************************************************************************************************************************************
2492//******************************************************************************************************************************************************
2493//******************************************************************************************************************************************************
2494//******************************************************************************************************************************************************
2495//******************************************************************************************************************************************************
2496//******************************************************************************************************************************************************
2497
2498
2499Function IR3A_PopMenuProc(pa) : PopupMenuControl
2500        STRUCT WMPopupAction &pa
2501
2502        NVAR SelectedLevel=root:Packages:AggregateModeling:SelectedLevel
2503        SVAR SlectedBranchedLevels=root:Packages:AggregateModeling:SlectedBranchedLevels
2504        NVAR SelectedQlevel=root:Packages:AggregateModeling:SelectedQlevel
2505        NVAR SelectedBlevel=root:Packages:AggregateModeling:SelectedBlevel
2506
2507        switch( pa.eventCode )
2508                case 2: // mouse up
2509                        Variable popNum = pa.popNum
2510                        String popStr = pa.popStr
2511                        string CtrlName=pa.ctrlName
2512                        if(stringMatch(CtrlName,"AvailableLevels"))
2513                                SelectedLevel = str2num(popStr[0,0])
2514                                SlectedBranchedLevels=popStr
2515                                IR3A_CalculateBranchedMassFr()
2516                        endif
2517                        if(stringMatch(CtrlName,"IntensityDataName")||stringMatch(CtrlName,"SelectDataFolder"))
2518                                IR2C_PanelPopupControl(pa)             
2519                                IR3A_SetControlsInPanel()       
2520                                IR3A_FindAvailableLevels()
2521                                IR3A_ClearVariables()
2522                                IR3A_CalculateBranchedMassFr()
2523                        endif   
2524                        if(stringMatch(CtrlName,"AllowedNearDistance"))
2525                                NVAR AllowedNearDistance=root:Packages:AggregateModeling:AllowedNearDistance
2526                                AllowedNearDistance = popNum
2527                        endif
2528                        if(stringMatch(CtrlName,"NUmberOfTestAggregates"))
2529                                NVAR NUmberOfTestAggregates=root:Packages:AggregateModeling:NUmberOfTestAggregates
2530                                NUmberOfTestAggregates = str2num(popStr)
2531                        endif
2532                       
2533                        if(stringMatch(CtrlName,"MParticleAttraction"))
2534                                SVAR MultiParticleAttraction=root:Packages:AggregateModeling:MultiParticleAttraction
2535                                MultiParticleAttraction = popStr
2536                        endif
2537                        break
2538        endswitch
2539
2540        return 0
2541End
2542//******************************************************************************************************************************************************
2543//******************************************************************************************************************************************************
2544//******************************************************************************************************************************************************
2545//******************************************************************************************************************************************************
2546//******************************************************************************************************************************************************
2547//******************************************************************************************************************************************************
2548//******************************************************************************************************************************************************
2549
2550
2551static Function IR3A_FindAvailableLevels()
2552       
2553        NVAR/Z UseCurrentResults=root:Packages:AggregateModeling:CurrentResults
2554        DoWIndow FractalAggregatePanel
2555        if(!NVAR_Exists(UseCurrentresults) || !V_Flag)
2556                return 0
2557        endif
2558       
2559        NVAR UseStoredResults=root:Packages:AggregateModeling:StoredResults
2560        String quote = "\""
2561
2562        SVAR Model=root:Packages:AggregateModeling:Model
2563        variable LNumOfLevels, i
2564       
2565        if(UseCurrentResults)
2566                NVAR/Z NumberOfLevels = root:Packages:Irena_UnifFit:NumberOfLevels
2567                if(!NVAR_Exists(NumberOfLevels))
2568                        return 0                        //no Unified Fit at all...
2569                endif
2570                LNumOfLevels = NumberOfLevels
2571        else
2572                LNumOfLevels =IR3A_ReturnNoteNumValue("NumberOfModelledLevels")
2573        endif
2574        string AvailableLevels=""
2575        if(stringmatch(Model,"Branched mass fractal")) 
2576//              if(LNumOfLevels>=1)
2577//                      AvailableLevels+=num2str(1)+";"
2578//              endif
2579                For(i=2;i<=LNumOfLevels;i+=1)
2580                        AvailableLevels+=num2str(i)+"/"+num2str(i-1)+";"//+num2str(i)+";"
2581                endfor
2582        else
2583                AvailableLevels=""
2584//              For(i=1;i<=LNumOfLevels;i+=1)
2585//                      AvailableLevels+=num2str(i)+";"
2586//              endfor
2587        endif
2588        string OnlyNumLevels=AvailableLevels
2589//      if(stringmatch(Model,"TwoPhase*"))     
2590//              AvailableLevels+="Range;All;"
2591//      endif   
2592//      if(stringmatch(Model,"Invariant*"))     
2593//              AvailableLevels+="Range;"
2594//      endif   
2595        AvailableLevels = quote + AvailableLevels + quote
2596        OnlyNumLevels = quote + OnlyNumLevels + quote
2597        NVAR SelectedQlevel=root:Packages:AggregateModeling:SelectedQlevel
2598        NVAR SelectedBlevel=root:Packages:AggregateModeling:SelectedBlevel
2599        SVAR SlectedBranchedLevels=root:Packages:AggregateModeling:SlectedBranchedLevels
2600        string loQStr="---"
2601        if(SelectedQlevel>0)
2602                loQStr=num2str(SelectedQlevel)
2603        endif
2604        string hiQStr="---"
2605        if(SelectedBlevel>0)
2606                hiQStr=num2str(SelectedBlevel)
2607        endif
2608        string AvLevStr="---"
2609        if(stringmatch(AvailableLevels,"*"+SlectedBranchedLevels+"*"))
2610                AvLevStr=SlectedBranchedLevels
2611        endif
2612        PopupMenu AvailableLevels,win=FractalAggregatePanel,mode=1,popvalue= AvLevStr, value=#AvailableLevels
2613//      PopupMenu SelectedQlevel,win=FractalAggregatePanel,mode=1,popvalue=loQStr, value=#OnlyNumLevels
2614//      PopupMenu SelectedBlevel,win=FractalAggregatePanel,mode=1,popvalue=hiQStr, value=#OnlyNumLevels
2615end
2616//******************************************************************************************************************************************************
2617//******************************************************************************************************************************************************
2618//******************************************************************************************************************************************************
2619//******************************************************************************************************************************************************
2620//******************************************************************************************************************************************************
2621//******************************************************************************************************************************************************
2622//******************************************************************************************************************************************************
2623
2624static  Function IR3A_ReturnNoteNumValue(KeyWord)
2625        string KeyWord
2626       
2627        variable LUKVal
2628        SVAR DataFolderName= root:Packages:AggregateModeling:DataFolderName
2629        SVAR IntensityWaveName = root:Packages:AggregateModeling:IntensityWaveName
2630       
2631        Wave/Z LkpWv=$(DataFolderName+IntensityWaveName)
2632        if(!WaveExists(LkpWv))
2633                return NaN
2634        endif
2635       
2636        LUKVal = NumberByKey(KeyWord, note(LkpWv)  , "=",";")
2637        return LUKVal
2638       
2639 end
2640//******************************************************************************************************************************************************
2641//******************************************************************************************************************************************************
2642//******************************************************************************************************************************************************
2643//******************************************************************************************************************************************************
2644//******************************************************************************************************************************************************
2645//******************************************************************************************************************************************************
2646//******************************************************************************************************************************************************
2647
2648
2649static Function IR3A_CalculateBranchedMassFr()
2650
2651        NVAR SelectedLevel = root:Packages:AggregateModeling:SelectedLevel
2652        SVAR SlectedBranchedLevels=root:Packages:AggregateModeling:SlectedBranchedLevels
2653        NVAR UseCurrentResults=root:Packages:AggregateModeling:CurrentResults
2654        NVAR UseStoredResults=root:Packages:AggregateModeling:StoredResults
2655
2656        NVAR BrFract_G2=root:Packages:AggregateModeling:BrFract_G2
2657        NVAR BrFract_Rg2=root:Packages:AggregateModeling:BrFract_Rg2
2658        NVAR BrFract_B2=root:Packages:AggregateModeling:BrFract_B2
2659        NVAR BrFract_P2=root:Packages:AggregateModeling:BrFract_P2
2660        NVAR BrFract_G1=root:Packages:AggregateModeling:BrFract_G1
2661        NVAR BrFract_Rg1=root:Packages:AggregateModeling:BrFract_Rg1
2662        NVAR BrFract_B1=root:Packages:AggregateModeling:BrFract_B1
2663        NVAR BrFract_P1=root:Packages:AggregateModeling:BrFract_P1
2664        SVAR BrFract_ErrorMessage=root:Packages:AggregateModeling:BrFract_ErrorMessage
2665        NVAR BrFract_dmin=root:Packages:AggregateModeling:BrFract_dmin
2666        NVAR BrFract_c=root:Packages:AggregateModeling:BrFract_c
2667        NVAR BrFract_z=root:Packages:AggregateModeling:BrFract_z
2668        NVAR BrFract_fBr=root:Packages:AggregateModeling:BrFract_fBr
2669        NVAR BrFract_fM=root:Packages:AggregateModeling:BrFract_fM
2670        NVAR BrFract_df=root:Packages:AggregateModeling:BrFract_df
2671       
2672        BrFract_ErrorMessage = ""
2673       
2674        if(stringMatch(SlectedBranchedLevels,"2/1"))
2675                SelectedLevel = 2
2676        elseif(stringMatch(SlectedBranchedLevels,"3/2"))
2677                SelectedLevel = 3
2678        elseif(stringMatch(SlectedBranchedLevels,"4/3"))
2679                SelectedLevel = 4
2680        elseif(stringMatch(SlectedBranchedLevels,"5/4"))
2681                SelectedLevel = 5
2682        endif
2683
2684        if(SelectedLevel>=2)
2685                if(UseCurrentResults)
2686                        NVAR gG2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"G")
2687                        BrFract_G2 = gG2
2688                        NVAR gRg2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"Rg")
2689                        BrFract_Rg2 = gRg2
2690                        NVAR gB2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"B")
2691                        BrFract_B2 = gB2
2692                        NVAR gP2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"P")
2693                        BrFract_P2 = gP2
2694                        NVAR gG1=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel-1)+"G")
2695                        BrFract_G1 = gG1
2696                        NVAR gRg1=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel-1)+"Rg")
2697                        BrFract_Rg1 = gRg1
2698                        NVAR gB1=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel-1)+"B")
2699                        BrFract_B1 = gB1
2700                        NVAR gP1=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel-1)+"P")
2701                        BrFract_P1 = gP1
2702                else
2703                        //look up from wave note...
2704                        BrFract_G2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"G")
2705                        BrFract_Rg2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"Rg")
2706                        BrFract_B2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"B")
2707                        BrFract_P2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"P")
2708                        BrFract_G1 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel-1)+"G")
2709                        BrFract_Rg1 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel-1)+"Rg")
2710                        BrFract_B1 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel-1)+"B")
2711                        BrFract_P1 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel-1)+"P")
2712                endif
2713        elseif(SelectedLevel==1)
2714                if(UseCurrentResults)
2715                        NVAR gG2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"G")
2716                        BrFract_G2 = gG2
2717                        NVAR gRg2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"Rg")
2718                        BrFract_Rg2 = gRg2
2719                        NVAR gB2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"B")
2720                        BrFract_B2 = gB2
2721                        NVAR gP2=$("root:Packages:Irena_UnifFit:Level"+num2str(SelectedLevel)+"P")
2722                        BrFract_P2 = gP2
2723                        BrFract_G1 =0
2724                        BrFract_Rg1 = 0
2725                        BrFract_B1 =0
2726                        BrFract_P1 = 0
2727                else
2728                        //look up from wave note...
2729                        BrFract_G2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"G")
2730                        BrFract_Rg2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"Rg")
2731                        BrFract_B2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"B")
2732                        BrFract_P2 = IR3A_ReturnNoteNumValue("Level"+num2str(SelectedLevel)+"P")
2733                        BrFract_G1 = 0
2734                        BrFract_Rg1 = 0
2735                        BrFract_B1 = 0
2736                        BrFract_P1 = 0
2737                endif
2738        else
2739                        BrFract_G2 = 0
2740                        BrFract_Rg2 = 0
2741                        BrFract_B2 = 0
2742                        BrFract_P2 = 0
2743                        BrFract_G1 = 0
2744                        BrFract_Rg1 = 0
2745                        BrFract_B1 = 0
2746                        BrFract_P1 = 0
2747        endif
2748        if(strlen(SlectedBranchedLevels)>1)
2749                BrFract_dmin  =BrFract_B2*BrFract_Rg2^(BrFract_P2)/(exp(gammln(BrFract_P2/2))*BrFract_G2)
2750                BrFract_c  =BrFract_P2/(BrFract_B2*BrFract_Rg2^(BrFract_P2)/(exp(gammln(BrFract_P2/2))*BrFract_G2))
2751                BrFract_z  =BrFract_G2/BrFract_G1 + 1                   //Greg, 11-24-2018: It should be G2/G1 +1  Karsten figured that out.  If G2 is 0 you still have one primary particle.
2752                BrFract_fBr =(1-(BrFract_G2/BrFract_G1)^(1/(BrFract_P2/(BrFract_B2*BrFract_Rg2^(BrFract_P2)/(exp(gammln(BrFract_P2/2))*BrFract_G2)))-1))
2753                BrFract_fM  = (1-(BrFract_G2/BrFract_G1)^(1/((BrFract_B2*BrFract_Rg2^(BrFract_P2)/(exp(gammln(BrFract_P2/2))*BrFract_G2)))-1))
2754                BrFract_df  = BrFract_P2
2755        else
2756                BrFract_dmin  =BrFract_B2*BrFract_Rg2^(BrFract_P2)/(exp(gammln(BrFract_P2/2))*BrFract_G2)
2757                BrFract_c  =BrFract_P2/(BrFract_B2*BrFract_Rg2^(BrFract_P2)/(exp(gammln(BrFract_P2/2))*BrFract_G2))
2758                BrFract_z  =Nan
2759                BrFract_fBr =NaN
2760                BrFract_fM  = NaN
2761        endif
2762       
2763        if(BrFract_c<0.96)
2764                BrFract_ErrorMessage =  "The mass fractal is too polydisperse to analyse, c < 1"
2765        else
2766                if(BrFract_c>=0.96 && BrFract_c<=1.04)//this should be in the range of 1 say +- .02
2767                        BrFract_ErrorMessage = "THIS IS A LINEAR CHAIN WITH NO BRANCHES!"
2768                endif
2769                if(BrFract_c>=3)
2770                        BrFract_ErrorMessage =  "There is a problem with the fit, c must be less than 3"
2771                endif
2772                if(BrFract_dmin>=3)
2773                        BrFract_ErrorMessage = "There is a problem with the fit since dmin must be less than 3"
2774                endif
2775                if(BrFract_dmin>=0.96 && BrFract_dmin<=1.04 )//this should be in the range of 1 say +- .02
2776                        BrFract_ErrorMessage = "This is a regular object, i.e.  c=1 rod, c=2 disk, c=3 sphere, etc."
2777                endif
2778        endif
2779       
2780        if(strlen(BrFract_ErrorMessage)>0)
2781                SetVariable BrFract_ErrorMessage, win=FractalAggregatePanel,  labelBack=(65535,49151,49151)
2782                beep
2783        else
2784                SetVariable BrFract_ErrorMessage, win=FractalAggregatePanel,  labelBack=0
2785        endif
2786
2787end
2788//******************************************************************************************************************************************************
2789//******************************************************************************************************************************************************
2790//******************************************************************************************************************************************************
2791//******************************************************************************************************************************************************
2792//******************************************************************************************************************************************************
2793//******************************************************************************************************************************************************
2794//******************************************************************************************************************************************************
2795
2796
2797Function IR3A_CheckProc(cba) : CheckBoxControl
2798        STRUCT WMCheckboxAction &cba
2799
2800        switch( cba.eventCode )
2801                case 2: // mouse up
2802                        Variable checked = cba.checked
2803                        NVAR CurrentResults=root:packages:AggregateModeling:CurrentResults
2804                        NVAR StoredResults=root:packages:AggregateModeling:StoredResults
2805                        SVAR Model=root:Packages:AggregateModeling:Model
2806                        if(stringMatch(cba.ctrlName,"CurrentResults"))
2807                                StoredResults=!CurrentResults
2808                        endif
2809                        if(stringMatch(cba.ctrlName,"StoredResults"))
2810                                CurrentResults=!StoredResults
2811                        endif
2812                        IR3A_UpdatePanelValues()
2813                        break
2814        endswitch
2815
2816        return 0
2817End
2818//******************************************************************************************************************************************************
2819//******************************************************************************************************************************************************
2820//******************************************************************************************************************************************************
2821Function IR3A_UpdatePanelValues()
2822
2823                        IR3A_SetControlsInPanel()       
2824                        IR3A_FindAvailableLevels()
2825                        IR3A_ClearVariables()
2826                        IR3A_CalculateBranchedMassFr()
2827                        IR3A_Create3DAggListForListbox()
2828
2829end
2830
2831//******************************************************************************************************************************************************
2832//******************************************************************************************************************************************************
2833//******************************************************************************************************************************************************
2834//******************************************************************************************************************************************************
2835
2836
2837static Function IR3A_SetControlsInPanel()
2838
2839                DOWindow FractalAggregatePanel
2840                if(!V_Flag)
2841                        return 0
2842                endif
2843                SVAR Model=root:Packages:AggregateModeling:Model
2844                NVAR CurrentResults = root:packages:AggregateModeling:CurrentResults
2845               
2846                PopupMenu SelectDataFolder, win=FractalAggregatePanel, disable=CurrentResults
2847                PopupMenu IntensityDataName, win=FractalAggregatePanel, disable=CurrentResults
2848                Setvariable FolderMatchStr, win=FractalAggregatePanel, disable=CurrentResults
2849
2850end
2851
2852//******************************************************************************************************************************************************
2853//******************************************************************************************************************************************************
2854//******************************************************************************************************************************************************
2855//******************************************************************************************************************************************************
2856//******************************************************************************************************************************************************
2857//******************************************************************************************************************************************************
2858//******************************************************************************************************************************************************
2859
2860static Function IR3A_ClearVariables()
2861
2862
2863        NVAR BrFract_G2=root:Packages:AggregateModeling:BrFract_G2
2864        NVAR BrFract_Rg2=root:Packages:AggregateModeling:BrFract_Rg2
2865        NVAR BrFract_B2=root:Packages:AggregateModeling:BrFract_B2
2866        NVAR BrFract_P2=root:Packages:AggregateModeling:BrFract_P2
2867        NVAR BrFract_G1=root:Packages:AggregateModeling:BrFract_G1
2868        NVAR BrFract_Rg1=root:Packages:AggregateModeling:BrFract_Rg1
2869        NVAR BrFract_B1=root:Packages:AggregateModeling:BrFract_B1
2870        NVAR BrFract_P1=root:Packages:AggregateModeling:BrFract_P1
2871        SVAR BrFract_ErrorMessage=root:Packages:AggregateModeling:BrFract_ErrorMessage
2872        NVAR BrFract_dmin=root:Packages:AggregateModeling:BrFract_dmin
2873        NVAR BrFract_c=root:Packages:AggregateModeling:BrFract_c
2874        NVAR BrFract_z=root:Packages:AggregateModeling:BrFract_z
2875        NVAR BrFract_fBr=root:Packages:AggregateModeling:BrFract_fBr
2876        NVAR BrFract_fM=root:Packages:AggregateModeling:BrFract_fM
2877        NVAR BrFract_df=root:Packages:AggregateModeling:BrFract_df
2878       
2879                BrFract_G2 = 0
2880                BrFract_Rg2 = 0
2881                BrFract_B2 = 0
2882                BrFract_P2 = 0
2883                BrFract_G1 = 0
2884                BrFract_Rg1 = 0
2885                BrFract_B1 = 0
2886                BrFract_P1 = 0
2887                BrFract_ErrorMessage=""
2888                BrFract_dmin=0
2889                BrFract_c=0
2890                BrFract_z=0
2891                BrFract_fBr=0
2892                BrFract_fM=0
2893                BrFract_df=0
2894
2895end
2896//******************************************************************************************************************************************************
2897//******************************************************************************************************************************************************
2898//******************************************************************************************************************************************************
2899//******************************************************************************************************************************************************
2900//******************************************************************************************************************************************************
2901//******************************************************************************************************************************************************
2902//******************************************************************************************************************************************************
2903
2904
2905Function IR3A_PanelButtonProc(ba) : ButtonControl
2906        STRUCT WMButtonAction &ba
2907
2908        variable failed = 0
2909        switch( ba.eventCode )
2910                case 2: // mouse up
2911                        // click code here
2912                        if(stringmatch(ba.ctrlName,"Grow1AggAll"))
2913                                failed = IR3A_Grow1MassFractAgreg()
2914                                if(!failed)
2915                                        IR3A_Calculate1DIntensity()
2916                                        IR3A_Display3DAggregate(1)
2917                                        IR3A_Create3DAggListForListbox()
2918                                endif
2919                        endif
2920                        if(stringmatch(ba.ctrlName,"GrowNAggAll"))
2921                                NVAR NUmberOfTestAggregates=root:Packages:AggregateModeling:NUmberOfTestAggregates
2922                                variable i
2923                                For(i=0;i<NUmberOfTestAggregates;i+=1)
2924                                        failed = IR3A_Grow1MassFractAgreg()
2925                                        if(!failed)
2926                                                IR3A_Calculate1DIntensity()
2927                                                IR3A_StoreCurrentMassFractAgreg()
2928                                                IR3A_Display3DAggregate(1)
2929                                                IR3A_Create3DAggListForListbox()
2930                                                print ">>>      Grown "+num2str(i+1)+" aggregate out of "+num2str(NUmberOfTestAggregates)
2931                                        endif
2932                                endfor
2933                                print "   >>>>     Done Growing "+num2str(NUmberOfTestAggregates)+" aggregates <<<<    "
2934                        endif
2935                        if(stringmatch(ba.ctrlName,"Display3DMassFracGizmo"))
2936                                IR3A_Display3DAggregate(0)
2937                                IR3A_Create3DAggListForListbox()
2938                        endif
2939                        if(stringmatch(ba.ctrlName,"SaveAggregateData"))
2940                                IR3A_StoreCurrentMassFractAgreg()
2941                                IR3A_Create3DAggListForListbox()
2942                        endif
2943                        if(stringmatch(ba.ctrlName,"Display3DMFASummary"))
2944                                Wave/Z MassFractalAggregate = root:Packages:AggregateModeling:MassFractalAggregate
2945                                if(WaveExists(MassFractalAggregate))
2946                                        IR3A_DisplayAggNotebook(MassFractalAggregate,0)
2947                                endif
2948                                IR3A_Create3DAggListForListbox()
2949                        endif
2950                        if(stringmatch(ba.ctrlName,"Calculate1DIntensity"))
2951                                IR3A_Calculate1DIntensity()
2952                                IR3A_Create3DAggListForListbox()
2953                        endif
2954                        if(stringmatch(ba.ctrlName,"Model1DIntensity"))
2955                                print "This takes a lot of time as it uses Monte Carlo method to calculate PDF and converts that to intensity. Also, results are noisy as single particle has very poor sampling."
2956                                IR3A_Model1DIntensity()
2957                                IR3A_Create3DAggListForListbox()
2958                        endif
2959                        if(stringmatch(ba.ctrlName,"Display1DData"))
2960                                IR3A_Display1DIntensity()
2961                                IR3A_Create3DAggListForListbox()
2962                        endif
2963                        if(StringMatch(ba.ctrlName, "GetHelp" ))
2964                                //Open www manual with the right page
2965                                IN2G_OpenWebManual("Irena/3DAggregate.html")
2966                        endif
2967                        if(StringMatch(ba.ctrlName, "CompareStoredResults" ))
2968                                IR3A_CompareStoredResults()
2969                        endif
2970                        if(StringMatch(ba.ctrlName, "DeleteStoredResults" ))
2971                                IR3A_DeleteStoredResults()
2972                        endif
2973
2974
2975                        break
2976                case -1: // control being killed
2977                        break
2978        endswitch
2979
2980        return 0
2981End
2982
2983//******************************************************************************************************************************************************
2984//******************************************************************************************************************************************************
2985//******************************************************************************************************************************************************
2986//******************************************************************************************************************************************************
2987//******************************************************************************************************************************************************
2988static Function IR3A_DeleteStoredResults()
2989
2990        DoAlert /T="Did you thnk about this?" 1, "You will delete all stored Mass Fractal AGgregates, really want to do it?"
2991
2992        if(V_FLag==1)
2993                DFref oldDf= GetDataFolderDFR()
2994
2995                KillWIndow/Z MassFractalAggregateView
2996                KillWIndow/Z MassFractalAggDataPlot
2997                KillWindow/Z AggStoredResultsOverview
2998                       
2999                KillDataFolder /Z root:MassFractalAggregates:
3000                IR3A_Create3DAggListForListbox()
3001                setDataFOlder OldDf
3002        endif
3003end
3004
3005//******************************************************************************************************************************************************
3006//******************************************************************************************************************************************************
3007//******************************************************************************************************************************************************
3008//******************************************************************************************************************************************************
3009//******************************************************************************************************************************************************
3010//******************************************************************************************************************************************************
3011//******************************************************************************************************************************************************
3012       
3013static Function IR3A_CompareStoredResults()
3014
3015        DFref oldDf= GetDataFolderDFR()
3016
3017
3018        DOWindow FractalAggregatePanel
3019        if(!V_Flag)
3020                return 0
3021        endif
3022        wave/Z Stored3DAggSelections= root:Packages:AggregateModeling:Stored3DAggSelections
3023        if(!WaveExists(Stored3DAggSelections))
3024                abort
3025        endif
3026        ///wave/T Stored3DAggregates= root:Packages:AggregateModeling:Stored3DAggregates
3027        //wave/T Stored3DAggregatesPaths= root:Packages:AggregateModeling:Stored3DAggregatesPaths
3028        variable NumOfFolders
3029        string CurrentList, tempStr
3030        if(DataFolderExists("root:MassFractalAggregates"))
3031                setDataFolder root:MassFractalAggregates
3032                CurrentList=stringByKey("FOLDERS",DataFolderDir(1),":")+","
3033                NumOfFolders = ItemsInList(CurrentList,",")+1
3034        else
3035                CurrentList=""
3036                NumOfFolders = 1
3037        endif
3038       
3039        if(NumOfFolders>0)
3040                redimension/N=(NumOfFolders) Stored3DAggSelections
3041                Stored3DAggSelections = 0
3042                setDataFolder root:Packages:AggregateModeling:
3043                make/O/N=(NumOfFolders-1) IndexStoredResWave, dMinStoredResWave, cValStoredResWave, dfStoredResWave
3044                make/O/N=(NumOfFolders-1) dMinTarget, cValTarget, dfTarget
3045                Wave IndexStoredResWave
3046                Wave dMinStoredResWave
3047                Wave cValStoredResWave
3048                Wave dfStoredResWave
3049                setDataFolder root:MassFractalAggregates
3050        else
3051                return 0
3052        endif   
3053               
3054        variable i
3055        if(NumOfFolders>1)
3056                For(i=0;i<NumOfFolders-1;i+=1)
3057                        tempStr = StringFromList(i, CurrentList, ",")
3058                        IndexStoredResWave[i] = i
3059                        cValStoredResWave[i] = IR3A_Return3DAggParamVal("root:MassFractalAggregates:"+tempStr,"cval")
3060                        dMinStoredResWave[i] = IR3A_Return3DAggParamVal("root:MassFractalAggregates:"+tempStr,"dMin")
3061                        dfStoredResWave[i] = IR3A_Return3DAggParamVal("root:MassFractalAggregates:"+tempStr,"df")
3062                endfor 
3063        endif
3064        variable MinVal, MaxVal
3065        //create plot of the three against the target values...
3066        NVAR dminTarg = root:Packages:AggregateModeling:BrFract_dmin
3067        NVAR cTarg = root:Packages:AggregateModeling:BrFract_c
3068        NVAR dfTarg = root:Packages:AggregateModeling:BrFract_df
3069        dMinTarget = dminTarg
3070        cValTarget = cTarg
3071        dfTarget = dfTarg
3072        KillWIndow/Z AggStoredResultsOverview
3073       
3074        Display/W=(695,66,1292,720)/K=1/N=AggStoredResultsOverview dMinStoredResWave vs IndexStoredResWave as "Overview of the Stored 3D Aggregates"
3075        AppendToGraph dMinTarget vs IndexStoredResWave
3076        AppendToGraph/L=cAxis cValStoredResWave,cValTarget vs IndexStoredResWave
3077        AppendToGraph/L=dfAxis dfStoredResWave,dfTarget vs IndexStoredResWave
3078        ModifyGraph mode(dMinStoredResWave)=3,mode(cValStoredResWave)=3,mode(dfStoredResWave)=3
3079        ModifyGraph marker(dMinStoredResWave)=19,marker(cValStoredResWave)=17,marker(dfStoredResWave)=26
3080        ModifyGraph lStyle(dMinTarget)=3,lStyle(cValTarget)=3,lStyle(dfTarget)=3
3081        ModifyGraph rgb(cValStoredResWave)=(1,16019,65535),rgb(cValTarget)=(1,16019,65535)
3082        ModifyGraph rgb(dfStoredResWave)=(3,52428,1),rgb(dfTarget)=(3,52428,1)
3083        ModifyGraph lblPosMode(cAxis)=3
3084        ModifyGraph lblPos(left)=63,lblPos(cAxis)=64,lblPos(dfAxis)=66
3085        ModifyGraph lblLatPos(left)=-8,lblLatPos(cAxis)=-6,lblLatPos(dfAxis)=-6
3086        ModifyGraph freePos(cAxis)=0
3087        ModifyGraph freePos(dfAxis)=0
3088        ModifyGraph axisEnab(left)={0,0.31}
3089        ModifyGraph axisEnab(cAxis)={0.35,0.64}
3090        ModifyGraph axisEnab(dfAxis)={0.68,1}
3091        ModifyGraph mirror=1
3092        ModifyGraph tick(bottom)=2
3093        ModifyGraph tick(left)=1
3094        ModifyGraph tick(cAxis)=1
3095        ModifyGraph tick(dfAxis)=1
3096       
3097        Label left "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"d\\Bmin"
3098        Label bottom "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Stored model index"
3099        Label cAxis "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"c value"
3100        Label dfAxis "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"df"
3101        wavestats/Q dMinStoredResWave
3102        MinVal=min(V_min,dminTarg)*0.95
3103        MaxVal=max(V_max,dminTarg)*1.05
3104        if(numtype(dminTarg)==0)
3105                SetAxis left MinVal,MaxVal
3106        else
3107                SetAxis/A left
3108        endif
3109       
3110        //
3111        wavestats/Q cValStoredResWave
3112        MinVal=min(V_min,cTarg)*0.95
3113        MaxVal=max(V_max,cTarg)*1.05
3114        if(numtype(cTarg)==0)
3115                SetAxis cAxis MinVal,MaxVal
3116        else
3117                SetAxis/A cAxis
3118        endif
3119        //
3120        wavestats/Q dfStoredResWave
3121        MinVal=min(V_min,dfTarg)*0.95
3122        MaxVal=max(V_max,dfTarg)*1.05
3123        if(numtype(dfTarg)==0)
3124                SetAxis dfAxis MinVal,MaxVal
3125        else
3126                SetAxis/A dfAxis
3127        endif
3128
3129        DrawLine 0,0.67,1,0.67
3130        DrawLine 0,0.34,1,0.34
3131        string LegendText="\\F"+IN2G_LkUpDfltStr("FontType")+"\\Z"+IN2G_LkUpDfltVar("LegendSize")+"\\s(dMinStoredResWave) d\\Bmin\r\\s(dMinTarget) d\\Bmin\\M  target\r\\s(cValStoredResWave) c\r\\s(cValTarget) c target\r\\s(dfStoredResWave) d\\Bf\r\\s(dfTarget) d\\Bf\\M target"
3132        Legend/C/N=LegendText/J/A=LT LegendText
3133        setDataFOlder OldDf
3134end     
3135//******************************************************************************************************************************************************
3136//******************************************************************************************************************************************************
3137
3138static Function IR3A_Return3DAggParamVal(PathToFolder, ParName)
3139        string PathToFOlder, ParName
3140        DFref oldDf= GetDataFolderDFR()
3141
3142        SetDataFolder PathToFolder
3143        variable RetValue
3144        if(StringMatch(ParName, "cval" ))
3145                NVAR cval= cValue
3146                RetValue = cval
3147        elseif(StringMatch(ParName, "dmin"))
3148                NVAR dMin=DminValue
3149                RetValue = dMin
3150        elseif(StringMatch(ParName, "df"))
3151                NVAR df = dfValue
3152                RetValue = df
3153        else
3154                RetValue = 0
3155        endif
3156        setDataFOlder OldDf
3157        return RetValue
3158end
3159
3160//******************************************************************************************************************************************************
3161//******************************************************************************************************************************************************
3162//******************************************************************************************************************************************************
3163//******************************************************************************************************************************************************
3164//******************************************************************************************************************************************************
3165
3166
3167static Function IR3A_StoreCurrentMassFractAgreg()
3168        DFref oldDf= GetDataFolderDFR()
3169
3170        SetDataFolder root:Packages:AggregateModeling
3171        NVAR OldDegreeOfAggregation=root:Packages:AggregateModeling:DegreeOfAggregation
3172        NVAR OldStickingProbability=root:Packages:AggregateModeling:StickingProbability
3173        NVAR OldRValue=root:Packages:AggregateModeling:RValue
3174        NVAR OldpValue=root:Packages:AggregateModeling:pValue
3175        NVAR OlddfValue=root:Packages:AggregateModeling:dfValue
3176        NVAR OlddminValue=root:Packages:AggregateModeling:dminValue
3177        NVAR OldcValue=root:Packages:AggregateModeling:cValue
3178        NVAR OldsValue=root:Packages:AggregateModeling:sValue
3179        NVAR OldAttemptValue=root:Packages:AggregateModeling:AttemptValue
3180        NVAR AllowedNearDistance=root:Packages:AggregateModeling:AllowedNearDistance
3181        Wave/Z MSF=root:Packages:AggregateModeling:MassFractalAggregate
3182        if(WaveExists(MSF))
3183                string NewFolderName
3184                NewFolderName = "MFA_DOA_"+num2str(OldDegreeOfAggregation)+"_Stick_"+num2str(OldStickingProbability)+"_StickMeth_"+num2str(AllowedNearDistance)+"_"
3185                NewDataFolder/O/S root:MassFractalAggregates
3186                NewFolderName = UniqueName(NewFolderName, 11, 0 )
3187                NewDataFolder/O/S $(NewFolderName)
3188                Duplicate/O MSF, MassFractalAggregate
3189                variable/g DegreeOfAggregation = OldDegreeOfAggregation
3190                variable/g StickingProbability = OldStickingProbability
3191                variable/g RValue = OldRValue
3192                variable/g pValue = OldpValue
3193                variable/g dfValue = OlddfValue
3194                variable/g dminValue = OlddminValue
3195                variable/g cValue = OldcValue
3196                variable/g sValue = OldsValue
3197                variable/g AttemptValue = OldAttemptValue               
3198                variable/g StickingMethod = AllowedNearDistance
3199        endif
3200        setDataFOlder OldDf
3201
3202
3203end
3204//******************************************************************************************************************************************************
3205//******************************************************************************************************************************************************
3206//******************************************************************************************************************************************************
3207//******************************************************************************************************************************************************
3208
3209//******************************************************************************************************************************************************
3210//******************************************************************************************************************************************************
3211//                      POV/PDB import/export/modeling code
3212//******************************************************************************************************************************************************
3213//******************************************************************************************************************************************************
3214//******************************************************************************************************************************************************
3215
3216Function IR3P_POVPDBPanel()
3217        PauseUpdate             // building window...
3218        NewPanel /K=1 /W=(5,20,395,680) as "POV/PDB Panel"
3219        DoWindow/C POVPDBPanel
3220        TitleBox MainTitle title="\Zr200POV/PDB panel",pos={20,0},frame=0,fstyle=3, fixedSize=1,font= "Times New Roman", size={350,24},anchor=MC,fColor=(0,0,52224)
3221        Button GetHelp,pos={305,45},size={80,15},fColor=(65535,32768,32768), proc=IR3P_POVPDBButtonProc,title="Get Help", help={"Open www manual page for this tool"}   //<<< fix button to help!!!
3222        //COPY FROM IR2U_UnifiedEvaPanelFnct()
3223        Checkbox  UseForPOV, pos={50,40}, size={50,15}, variable =root:packages:POVPDBImport:UseForPOV
3224        Checkbox  UseForPOV, title="Use for POV",mode=1,proc=IR3P_POVPDBCheckProc
3225        Checkbox  UseForPOV, help={"Select of you want to import POV files from SAXSMorph"}
3226        Checkbox  UseForPDB, pos={200,40}, size={50,15}, variable =root:packages:POVPDBImport:UseForPDB
3227        Checkbox  UseForPDB, title="Use for PDB",mode=1,proc=IR3P_POVPDBCheckProc
3228        Checkbox  UseForPDB, help={"Select of you want to import PDB files from SAXSMorph"}
3229
3230        TitleBox Info1 title="\Zr120Select where to put the data",pos={60,65},frame=0,fstyle=1, fixedSize=1,size={300,20},fColor=(0,0,52224)
3231        SetVariable NewFolderName,value= root:Packages:POVPDBImport:NewFolderName
3232        SetVariable NewFolderName,pos={15,90},size={200,20},title="root:",noproc, help={"Type in new folder name"}
3233        Button CreateFolder,pos={250,87},size={100,20},proc=IR3P_POVPDBButtonProc,title="Create Folder", help={"Create Folder for data"}
3234
3235        SetVariable CurrentFolderName,value= root:Packages:POVPDBImport:CurrentFolderName
3236        SetVariable CurrentFolderName,pos={15,120},size={280,20},title="Current Folder Name",noproc, help={"Current FOlder name to use"}
3237
3238        TitleBox FakeLine1 title=" ",fixedSize=1,size={330,3},pos={16,145},frame=0,fColor=(0,0,52224), labelBack=(0,0,52224)
3239        TitleBox Info2 title="\Zr120Import data",pos={60,160},frame=0,fstyle=1, fixedSize=1,size={300,20},fColor=(0,0,52224)
3240
3241
3242        SetVariable voxelSize,value= root:packages:POVPDBImport:voxelSize, frame=1
3243        SetVariable voxelSize,pos={15,190},noproc, help={"WHat is voxel Size of the imported 3D structure? "}
3244        SetVariable voxelSize title="Voxel Size [A]                ",size={200,17},limits={1,100,1}
3245
3246
3247        Button Import3DData,pos={30,220},size={150,20},proc=IR3P_POVPDBButtonProc,title="Import 3D Data", help={"Import 3D data in this folder"}
3248        Button Display3DData,pos={30,250},size={150,20},proc=IR3P_POVPDBButtonProc,title="Display 3D Data", help={"Display 3D data in this folder"}
3249        Button ImportIntQData,pos={30,300},size={150,20},proc=IR3P_POVPDBButtonProc,title="Import Int/Q Data", help={"Import 1D data in this folder"}
3250        Button DisplayIntQData,pos={30,330},size={150,20},proc=IR3P_POVPDBButtonProc,title="Display Int/Q Data", help={"Display 1D data in this folder"}
3251        Button CalculateIntQData,pos={30,360},size={150,20},proc=IR3P_POVPDBButtonProc,title="Calculate Int/Q Data", help={"Calculate 1D data and append"}
3252
3253end
3254
3255//******************************************************************************************************************************************************
3256//******************************************************************************************************************************************************
3257//*****************************************************************************************************************
3258//******************************************************************************************************************************************************
3259//******************************************************************************************************************************************************
3260Function IR3P_POVPDBCheckProc(cba) : CheckBoxControl
3261        STRUCT WMCheckboxAction &cba
3262
3263        switch( cba.eventCode )
3264                case 2: // mouse up
3265                        Variable checked = cba.checked
3266                        NVAR UseForPOV=root:packages:POVPDBImport:UseForPOV
3267                        NVAR UseForPDB=root:packages:POVPDBImport:UseForPDB
3268                        if(stringMatch(cba.ctrlName,"UseForPOV"))
3269                                UseForPDB=!UseForPOV
3270                        endif
3271                        if(stringMatch(cba.ctrlName,"UseForPDB"))
3272                                UseForPOV=!UseForPDB
3273                        endif
3274                        break
3275        endswitch
3276
3277        return 0
3278End
3279//******************************************************************************************************************************************************
3280//******************************************************************************************************************************************************
3281//******************************************************************************************************************************************************
3282//******************************************************************************************************************************************************
3283Function IR3P_POVPDBButtonProc(ba) : ButtonControl
3284        STRUCT WMButtonAction &ba
3285
3286        switch( ba.eventCode )
3287                case 2: // mouse up
3288                        // click code here
3289                        DFref oldDf= GetDataFolderDFR()
3290
3291                        setDataFolder root:Packages:POVPDBImport       
3292                        if(StringMatch(ba.ctrlName, "CreateFolder" ))
3293                                        IR3P_CreateFolder()     
3294                        endif
3295                        if(StringMatch(ba.ctrlName, "GetHelp" ))
3296                                        print "Fix IR3P_POVPDBButtonProc to do what it is suppose to do..."
3297                        endif
3298                        if(StringMatch(ba.ctrlName, "Import3DData" ))
3299                                        IR3P_Read3DDataFile()
3300                        endif
3301                        if(StringMatch(ba.ctrlName, "Display3DData" ))
3302                                        IR3P_POV3DDataGizmo()
3303                        endif
3304                        if(StringMatch(ba.ctrlName, "ImportIntQData" ))
3305                                        IR3P_Read1DDataFile()
3306                        endif
3307                        if(StringMatch(ba.ctrlName, "DisplayIntQData" ))
3308                                        IR3P_DIsplay1DDataFile()
3309                        endif
3310                        if(StringMatch(ba.ctrlName, "CalculateIntQData" ))
3311                                        IR3P_Calculate1DDataFile()
3312                        endif
3313                       
3314                        setDataFolder oldDF             
3315                        break
3316                case -1: // control being killed
3317                        break
3318        endswitch
3319
3320        return 0
3321End
3322//******************************************************************************************************************************************************
3323//******************************************************************************************************************************************************
3324//******************************************************************************************************************************************************
3325//******************************************************************************************************************************************************
3326//******************************************************************************************************************************************************
3327Function IR3P_InitializePOVPDB()
3328
3329        //IN2G_PrintDebugStatement(IrenaDebugLevel, 5,"")
3330        DFref oldDf= GetDataFolderDFR()
3331
3332        NewDataFolder/O/S root:Packages
3333        NewDataFolder/O/S root:Packages:POVPDBImport
3334        string/g ListOfVariables
3335        string/g ListOfStrings
3336        //here define the lists of variables and strings needed, separate names by ;...
3337        ListOfVariables="UseForPOV;UseForPDB;VoxelSize;"
3338        ListOfStrings="NewFolderName;CurrentFolderName;"
3339        variable i
3340        //and here we create them
3341        for(i=0;i<itemsInList(ListOfVariables);i+=1)   
3342                IN2G_CreateItem("variable",StringFromList(i,ListOfVariables))
3343        endfor                                                                                         
3344        for(i=0;i<itemsInList(ListOfStrings);i+=1)     
3345                IN2G_CreateItem("string",StringFromList(i,ListOfStrings))
3346        endfor 
3347               
3348        NVAR VoxelSize
3349        if(VoxelSize<1)
3350                VoxelSize = 1           //default to 1A
3351        endif
3352        setDataFOlder OldDf
3353end
3354//******************************************************************************************************************************************************
3355//******************************************************************************************************************************************************
3356////******************************************************************************************************************************************************
3357//
3358//Window POV3DData() : GizmoPlot
3359//      PauseUpdate             // building window...
3360//      // Building Gizmo 8 window...
3361//      NewGizmo/T="POV Imported 3D data"/W=(919,45,1434,505)
3362//      ModifyGizmo startRecMacro=700
3363//      ModifyGizmo scalingOption=63
3364//      AppendToGizmo isoSurface=root:test:centersWave,name=isoSurface0
3365//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ surfaceColorType,1}
3366//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineColorType,0}
3367//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidthType,0}
3368//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ fillMode,2}
3369//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidth,1}
3370//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ isoValue,0.5}
3371//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ frontColor,1,0,0,1}
3372//      ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ backColor,0,0,1,1}
3373//      ModifyGizmo modifyObject=isoSurface0,objectType=Surface,property={calcNormals,1}
3374//      AppendToGizmo light=Directional,name=light0
3375//      ModifyGizmo modifyObject=light0,objectType=light,property={ position,-0.241800,-0.664500,0.707100,0.000000}
3376//      ModifyGizmo modifyObject=light0,objectType=light,property={ direction,-0.241800,-0.664500,0.707100}
3377//      ModifyGizmo modifyObject=light0,objectType=light,property={ ambient,0.133000,0.133000,0.133000,1.000000}
3378//      ModifyGizmo modifyObject=light0,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
3379//      AppendToGizmo Axes=boxAxes,name=axes0
3380//      ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
3381//      ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
3382//      ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
3383//      AppendToGizmo attribute specular={1,1,0,1,1032},name=specular0
3384//      AppendToGizmo attribute shininess={5,20},name=shininess0
3385//      ModifyGizmo setDisplayList=0, object=light0
3386//      ModifyGizmo setDisplayList=1, attribute=shininess0
3387//      ModifyGizmo setDisplayList=2, attribute=specular0
3388//      ModifyGizmo setDisplayList=3, object=isoSurface0
3389//      ModifyGizmo setDisplayList=4, opName=clearColor, operation=clearColor, data={0.8,0.8,0.8,1}
3390//      ModifyGizmo setDisplayList=5, object=axes0
3391//      ModifyGizmo autoscaling=1
3392//      ModifyGizmo currentGroupObject=""
3393//      ModifyGizmo showInfo
3394//      ModifyGizmo infoWindow={1436,23,2253,322}
3395//      ModifyGizmo endRecMacro
3396//      ModifyGizmo SETQUATERNION={-0.134543,0.320522,-0.081181,0.934117}
3397//EndMacro
3398//******************************************************************************************************************************************************
3399//******************************************************************************************************************************************************
3400//******************************************************************************************************************************************************
3401
3402
3403
3404
3405//******************************************************************************************************************************************************
3406//******************************************************************************************************************************************************
3407Function IR3P_CreateFolder()   
3408        SVAR NewFolderName=root:packages:POVPDBImport:NewFolderName
3409        SVAR CurrentFolderName=root:packages:POVPDBImport:CurrentFolderName
3410        if(Strlen(NewFolderName)>2)
3411                setDataFOlder root:
3412                NewDataFOlder/O/S $(PossiblyQuoteName(NewFolderName))
3413                CurrentFolderName = GetDataFolder(1)
3414        else
3415                Abort "No Folder name exists, type in name first"
3416        endif
3417end
3418//******************************************************************************************************************************************************
3419//******************************************************************************************************************************************************
3420Function IR3P_Read1DDataFile()
3421
3422        SVAR CurrentFolderName=root:packages:POVPDBImport:CurrentFolderName
3423        SetDataFOlder $(CurrentFolderName)
3424        IR1I_KillAutoWaves()
3425        KillWaves/Z Qwave, IntWave
3426        LoadWave/Q/A/D/G/D
3427        Wave Wave0
3428        Wave Wave1
3429        Rename Wave0, Qwave
3430        Rename Wave1, IntWave
3431end
3432///*************************************************************************************************************************************
3433///*************************************************************************************************************************************
3434Function IR3P_DIsplay1DDataFile()
3435        DOWIndow/K/Z POV1DGraph
3436        SVAR CurrentFolderName=root:packages:POVPDBImport:CurrentFolderName
3437        SetDataFOlder $(CurrentFolderName)
3438        Wave IntWave
3439        Wave QWave
3440        Display/K=1/W=(100,50,600,550) IntWave vs QWave as "POV 1D Data display"
3441        DoWindow/C POV1DGraph
3442        Label left "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Intensity"
3443        Label bottom "\\Z"+IN2G_LkUpDfltVar("AxisLabelSize")+"Qvector"
3444        ModifyGraph log=1
3445end
3446
3447///*************************************************************************************************************************************
3448///*************************************************************************************************************************************
3449Function IR3P_Calculate1DDataFile()
3450        SVAR CurrentFolderName=root:packages:POVPDBImport:CurrentFolderName
3451        SetDataFOlder $(CurrentFolderName)
3452        Wave ThreeDVoxelGram = POVVoxelWave
3453        NVAR voxelSize = root:packages:POVPDBImport:voxelSize
3454        print "Warning = IR3P_Calculate1DDataFile needs to get VoxelSize fixzed, it is fixerd at 2A"
3455        //variable NumRSteps=300
3456        variable IsoValue = 0.5
3457        variable Qmin = 0.001
3458        variable Qmax = 0.6
3459        variable NumQSteps = 200
3460        setScale/P x, 0, VoxelSize, ThreeDVoxelGram
3461        setScale/P y, 0, VoxelSize, ThreeDVoxelGram
3462        setScale/P z, 0, VoxelSize, ThreeDVoxelGram
3463        IR3T_CreatePDFIntensity(ThreeDVoxelGram, IsoValue, Qmin, Qmax, NumQSteps)
3464        Wave PDFQWv
3465        Wave PDFIntensityWv
3466        Wave Qwave
3467        Wave IntWave
3468        //and this does not work for odd number of rows/columns/...
3469        //IR3T_CalcAutoCorelIntensity(ThreeDVoxelGram, IsoValue, Qmin, Qmax, NumQSteps)
3470        //these are autocorrelation calculated intensities...
3471        //Wave AutoCorIntensityWv
3472        //Wave AutoCorQWv
3473        variable InvarModel=areaXY(PDFQWv, PDFIntensityWv )
3474        variable InvarData=areaXY(Qwave, IntWave )
3475        PDFIntensityWv*=InvarData/InvarModel
3476        //InvarModel=areaXY(AutoCorQWv, AutoCorIntensityWv )
3477        //AutoCorIntensityWv*=InvarData/InvarModel
3478       
3479        DOWIndow POV1DGraph
3480        if(V_Flag)
3481                DoWIndow/F POV1DGraph
3482                CheckDisplayed /W=POV1DGraph PDFIntensityWv
3483                if(V_flag==0)
3484                        AppendToGraph/W=POV1DGraph  PDFIntensityWv vs PDFQWv
3485                endif
3486                //CheckDisplayed /W=POV1DGraph AutoCorIntensityWv
3487                //if(V_flag==0)
3488                //      AppendToGraph/W=POV1DGraph  AutoCorIntensityWv vs AutoCorQWv
3489                //endif
3490                ModifyGraph lstyle(PDFIntensityWv)=9,lsize(PDFIntensityWv)=3,rgb(PDFIntensityWv)=(1,16019,65535)
3491                ModifyGraph mode(PDFIntensityWv)=4,marker(PDFIntensityWv)=19
3492                ModifyGraph msize(PDFIntensityWv)=3
3493                //ModifyGraph lsize(AutoCorIntensityWv)=3,rgb(AutoCorIntensityWv)=(3,52428,1)
3494                Legend/C/N=text0/A=MC
3495        endif
3496end
3497///*************************************************************************************************************************************
3498///*************************************************************************************************************************************
3499Function IR3P_Read3DDataFile()
3500        //variable dimx,dimy,dimz
3501       
3502        NVAR UseForPOV = root:Packages:POVPDBImport:UseForPOV
3503        NVAR UseForPDB = root:Packages:POVPDBImport:UseForPDB
3504        if(UseForPOV+UseForPDB !=1)
3505                UseForPDB=0
3506                UseForPOV=1
3507        endif
3508        if(UseForPOV)
3509                SVAR CurrentFolderName=root:packages:POVPDBImport:CurrentFolderName
3510                //SetDataFOlder $(CurrentFolderName)
3511                Variable refNum,err=0
3512                variable FInalSize
3513                OPEN/R/F=".POV"/M="Find POV file" refNum
3514                if(strlen(S_fileName)>0)
3515                        Make/O/U/B/Free/n=((500),(500),(500)) centersWave
3516                        centersWave=0
3517                        String lineStr
3518                        Variable count=0
3519                        do
3520                                FreadLine refNum,lineStr
3521                                if(strlen(lineStr)<=0)
3522                                        break
3523                                endif
3524                                if(strsearch(lineStr,"sphere",0)>=0)
3525                                        IR3P_POVprocessAtomLine(lineStr,count,centersWave)
3526                                        count+=1
3527                                endif
3528                        while(err==0)
3529                        Close refNum
3530                        FinalSize = count^(1/3)
3531                        //print FinalSize
3532                        Redimension/N=(FInalSize,FInalSize,FInalSize) centersWave
3533                        Duplicate/O centersWave, $(CurrentFolderName+"POVVoxelWave")
3534                endif
3535        elseif(UseForPDB)
3536                print "Finish PDB in IR3P_Read3DDataFile"
3537        endif
3538end
3539///*************************************************************************************************************************************
3540///*************************************************************************************************************************************
3541///*************************************************************************************************************************************
3542///*************************************************************************************************************************************
3543Function IR3P_POVprocessAtomLine(lineStr,count, destWv)
3544        String lineStr
3545        Variable count
3546        Wave destWv
3547       
3548        Variable n1,n2,n3,n4,n5,xx,yy,zz
3549        String s1,s2,s3,s4,s5
3550        //sphere{<50,48,33>, 1.5, 1 }
3551        lineStr = ReplaceString("{<", lineStr, ",")
3552        lineStr = ReplaceString(">,", lineStr, ",")
3553        lineStr = ReplaceString("} \r", lineStr, ",")
3554        lineStr = ReplaceString(" ", lineStr, "")
3555        sscanf lineStr,"sphere,%i,%i,%i,%f,%f,",n1,n2,n3,n4,n5
3556        variable FillVal = 0
3557        if(n5>0)
3558                FillVal = 1
3559        endif
3560        destWv[(n1-1)][(n2-1)][(n3-1)]= FillVal
3561End
3562//******************************************************************************************************************************************************
3563//******************************************************************************************************************************************************
3564//******************************************************************************************************************************************************
3565
3566Function IR3P_POV3DDataGizmo() : GizmoPlot
3567        DoWIndow/K/Z POV3DData
3568        PauseUpdate             // building window...
3569        SVAR CurrentFolderName=root:packages:POVPDBImport:CurrentFolderName
3570        Wave Imported3DPOVWave = $(CurrentFolderName+"POVVoxelWave")
3571        // Building Gizmo 8 window...
3572        NewGizmo/T="POV Imported 3D data"/W=(919,45,1434,505)
3573        DoWindow/C POV3DData
3574        ModifyGizmo startRecMacro=700
3575        ModifyGizmo scalingOption=63
3576        AppendToGizmo isoSurface=Imported3DPOVWave,name=ImportedPOV
3577        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ surfaceColorType,1}
3578        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ lineColorType,0}
3579        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ lineWidthType,0}
3580        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ fillMode,2}
3581        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ lineWidth,1}
3582        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ isoValue,0.5}
3583        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ frontColor,1,0,0,1}
3584        ModifyGizmo ModifyObject=ImportedPOV,objectType=isoSurface,property={ backColor,0,0,1,1}
3585        ModifyGizmo modifyObject=ImportedPOV,objectType=Surface,property={calcNormals,1}
3586        AppendToGizmo light=Directional,name=light0
3587        ModifyGizmo modifyObject=light0,objectType=light,property={ position,-0.241800,-0.664500,0.707100,0.000000}
3588        ModifyGizmo modifyObject=light0,objectType=light,property={ direction,-0.241800,-0.664500,0.707100}
3589        ModifyGizmo modifyObject=light0,objectType=light,property={ ambient,0.133000,0.133000,0.133000,1.000000}
3590        ModifyGizmo modifyObject=light0,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
3591        AppendToGizmo Axes=boxAxes,name=axes0
3592        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
3593        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
3594        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
3595        AppendToGizmo attribute specular={1,1,0,1,1032},name=specular0
3596        AppendToGizmo attribute shininess={5,20},name=shininess0
3597        ModifyGizmo setDisplayList=0, object=light0
3598        ModifyGizmo setDisplayList=1, attribute=shininess0
3599        ModifyGizmo setDisplayList=2, attribute=specular0
3600        ModifyGizmo setDisplayList=3, object=ImportedPOV
3601        ModifyGizmo setDisplayList=4, opName=clearColor, operation=clearColor, data={0.8,0.8,0.8,1}
3602        ModifyGizmo setDisplayList=5, object=axes0
3603        ModifyGizmo autoscaling=1
3604        ModifyGizmo currentGroupObject=""
3605        ModifyGizmo showInfo
3606        ModifyGizmo infoWindow={1436,23,2253,322}
3607        ModifyGizmo endRecMacro
3608        ModifyGizmo SETQUATERNION={-0.134543,0.320522,-0.081181,0.934117}
3609EndMacro
3610///*************************************************************************************************************************************
3611///*************************************************************************************************************************************
3612///*************************************************************************************************************************************
3613//
3614//Function IR3P_Convert3DMatrixToList(My3DVoxelWave, threshVal)
3615//      wave My3DVoxelWave
3616//      variable threshVal
3617//     
3618//      variable maxlen=DimSize(My3DVoxelWave, 0 )*DimSize(My3DVoxelWave,1)*DimSize(My3DVoxelWave, 2)
3619//     
3620//      Make/O/N=(maxlen,3) MyScatterWave
3621//      MyScatterWave = NaN
3622//      variable i, j, k, indx=0
3623//      For(i=0;i<DimSize(My3DVoxelWave, 0 );i+=1)
3624//              For(j=0;j<DimSize(My3DVoxelWave, 1 );j+=1)
3625//                      For(k=0;k<DimSize(My3DVoxelWave, 2 );k+=1)
3626//                              if(My3DVoxelWave[i][j][k]<threshVal)
3627//                                      MyScatterWave[indx][0]=i
3628//                                      MyScatterWave[indx][1]=j
3629//                                      MyScatterWave[indx][2]=k
3630//                                      indx+=1
3631//                              endif
3632//                      endfor
3633//              endfor
3634//      endfor
3635//      DeletePoints/M=0 indx, (maxlen-indx),  MyScatterWave
3636//     
3637//end
3638///*************************************************************************************************************************************
3639///*************************************************************************************************************************************
3640///*************************************************************************************************************************************
3641 
Note: See TracBrowser for help on using the repository browser.