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

Last change on this file since 1171 was 1171, checked in by ilavsky, 9 months ago

Add Search for optimal solution to 3D aggregates

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