1. Quick start in MATLAB - SSRL



User Guide for GIWAXS-SIIRkit?data analysis suiteRefer to J. Appl. Cryst. (2020); : scattering intensity, indexing and refraction calculation toolkit for grazing-incidence wide-angle X-ray scattering of organic materialsV. Savikhin, H.-G. Steinrück, R.-Z. Liang, B. A. Collins, S. D. Oosterhout, P. M. Beaujuge and M. F. ToneyContents TOC \o "1-3" \h \z \u 1. Quick start in MATLAB PAGEREF _Toc16362327 \h 42. Complex refractive index calculator PAGEREF _Toc16362328 \h 63. Scattering intensity analysis PAGEREF _Toc16362329 \h 7A) Calculating intensity versus a single variable PAGEREF _Toc16362330 \h 8B) Calculating intensity versus two variables PAGEREF _Toc16362331 \h 9C) Entering a list of samples to extract a normalization factor PAGEREF _Toc16362332 \h 9D) Error bar analysis for a single sample PAGEREF _Toc16362333 \h 104. Scattering peak position analysis PAGEREF _Toc16362334 \h 125. Peak indexing tool PAGEREF _Toc16362335 \h 14A) Extracting peaks PAGEREF _Toc16362336 \h 14B) Uploading peaks from a file PAGEREF _Toc16362337 \h 16C) Indexing peaks: “Step 1” and “Step 2” PAGEREF _Toc16362338 \h 16D) Viewing results PAGEREF _Toc16362339 \h 176. Tips for navigating output data and customizing code PAGEREF _Toc16362340 \h 21A) General plot manipulation in MATLAB PAGEREF _Toc16362341 \h 21i. Change plot appearance PAGEREF _Toc16362342 \h 21ii. Obtain data from a figure PAGEREF _Toc16362343 \h 21iii. Obtain data from table PAGEREF _Toc16362344 \h 23B) Modifying default values of GUI elements PAGEREF _Toc16362345 \h 24C) Useful modifications to scattering intensity analysis PAGEREF _Toc16362346 \h 24i. Changing the beam profile in space PAGEREF _Toc16362347 \h 24ii. Changing the beam divergence distribution PAGEREF _Toc16362348 \h 26ii. Adding or modifying shortcut substrate refractive indices PAGEREF _Toc16362349 \h 26iii. Changing the convergence condition for error bar calculation or error bar definition PAGEREF _Toc16362350 \h 27iv. Isolating scattering intensity from one layer of the film PAGEREF _Toc16362351 \h 28v. Changing the number of points simulated PAGEREF _Toc16362352 \h 29D) Useful modifications to peak indexing tool PAGEREF _Toc16362353 \h 29i. Changing Qxy-Qz image file type PAGEREF _Toc16362354 \h 29ii. Adding a polar angle “sin(χ)” intensity correction to Qxy-Qz image PAGEREF _Toc16362355 \h 30iii. Adding a “weight” parameter to peak indexing position error PAGEREF _Toc16362356 \h 31iv. Changing refinement parameters for “Step 1” and “Step 2” PAGEREF _Toc16362357 \h 321. Quick start in MATLABThis toolbox runs in MATLAB. It does not use any advanced toolboxes and should be fully functional using the basic model of MATLAB, version 2012a or newer. If problems are encountered, please e-mail mftoney@slac.stanford.edu or vsavikhi@.More detail about the calculations performed can be found in the functionality provided in this toolbox can be accessed through one of four files:I-indexing/IMain.m (for crystal structure calculation)SI-scattering intensity/SIMain.m (for scattering intensity normalization and error calculation)R-refraction/RMain.m (for refractive shift calculation)R-refraction/RefractiveIndexCalculator.m (for refractive index calculator)There is also a standalone tool under I-indexing/EquivalentCells.m which can be used calculate a unit cell that is equivalent to a user-specified unit cell but with different unit cell vectors. This tool does not have a GUI.All the other files included in the .zip folder are also necessary for functionality and should NOT be deleted.To launch one of these files, either double-click to open the file in MATLAB or use the File->Open dropdown from the MATLAB interface. This will open the raw code in the “Editor” window shown in REF _Ref9186622 \h Figure 1A. Launch the file by clicking on the green arrow in the top toolbox (circled in REF _Ref9186622 \h Figure 1A), or by hitting F5. If a pop-up window like the one shown in REF _Ref9186622 \h Figure 1B appears, select “Change Folder” to continue.Figure SEQ Figure \* ARABIC 1. A) MATLAB Editor window, showing the IMain.m tool. The tool can be launched by clicking the green arrow (outlined), or by hitting F5. B) Possible pop-up window on first launch; click “Change Folder.”If an error occurs, a message will most likely appear in red font in the MATLAB command window. An error may also result in unexpected output (such as no data in an output plot) which may or may not be accompanied by an orange warning message. This should not happen but if it does, please send an email including the error or warning message and screenshot of your GUI to vsavikhi@.If a calculation is taking an unexpectedly long time, it is usually possible to terminate the calculation (leaving all windows open) by hitting Ctrl+C. More time-consuming calculations typically display a progress bar.For the peak indexing toolset, we have included a few example files in the “Sample Files” folder that can be used to test the code. 2. Complex refractive index calculatorThe refractive index calculator is found in the file R-refraction/RefractiveIndexCalculator.m. When launched, the graphical user interface (GUI) shown in REF _Ref9187310 \h Figure 2A appears. The beam energy and material density should be adjusted as needed. The user must additionally provide an atomic composition OR an average atomic number, mass, and attenuation coefficient. These two modes are toggled with the radiobuttons ?. The atomic composition should be entered using standard periodic table atomic labels followed by numbers, which may include decimal places. In addition, parentheses can be used to group parts of the atomic composition together. Parentheses may be nested. For example, (C2H4(C3H6)4.3S)3 →2C+4H+3C+6H?4.3+S?3=44.7C+89.4H+3S.Additional symbols and spaces can be added to improve legibility. An additional window will pop up to confirm that the structure was read correctly, as shown in REF _Ref9187310 \h Figure 2C.Launching the tool using an atomic composition will automatically update the average atomic number, mass, and attenuation coefficient. The resultant β and δ, where n=1-δ+iβ, are displayed in the main window, as shown in REF _Ref9187310 \h Figure 2B. The corresponding f1 and f2 are also shown.When launching the tool using the second mode, f1 and f2 are not available. β is calculated from the attenuation coefficient and δ is approximated.Figure SEQ Figure \* ARABIC 2. A) Interface of the refractive index calculator on startup. B) Interface after launching, populated with results. C) Additional pop-up seen on launch, which can be used to confirm correct reading of chemical structure. 3. Scattering intensity analysisThe intensity analysis tool is found in SI-scattering intensity/SIMain.m. When launched, the GUI shown in REF _Ref9187859 \h Figure 3A appears. Each of the variables should be updated as needed, with the exception of unused ± boxes which will be disabled. These variables are described in REF _Ref9189619 \h Table 1.incident angle (°)=the angle between the direct beam and the substrate surfacesample pos (mm)=the sample position along the beam direction where “0” means that the center of the sample aligns with the center of the beam.beam energy (keV)=the energy of the incident beamsamp length (mm)=sample size along the beam directionhoriz FWHM (um)=the horizontal width of the beam in micrometers μmsamp width (mm)=sample size perpendicular to the beam directionvert FWHM (um)=the vertical width of the beamdensity (g/cm^3)=average density of the sampledivergence (±°)=the angular divergence of the beam (assumed Gaussian)thickness (nm)=thickness of scattering material (e.g. thickness of a film which is on top of a substrate)beam shape=the shape of the beam along the horizontal and vertical directions: drop-down includes Gaussian, Lorentzian, and rectangularTable SEQ Table \* ARABIC 1. Descriptions of variables specified in scattering intensity interface.In addition, the β and δ (“Beta” and “Delta”) of both film and substrate are specified. A dropdown is included to automatically fill in the substrate refractive index based off calculated values for several common substrate materials: Silicon, SiO2, microscope slide (composed of 73.3% SiO2, 13.5% Na2O, 8.9% CaO, and 4.4% MgO by weight, and a density of 2.5 g/cm3), Si3N4, and ITO (Indium Tin Oxide: composed of 74% In, 18% O2, and 8% Sn by weight, and a density of 7.14 g/cm3). These are given at the standard 11-3 beamline energy of 12.735 keV. Instructions for adding additional shortcut materials are given in Section 6C. The film refractive index should be specified for the “center” density value and will be scaled appropriately if density is selected as a variable. On the other hand, the film refractive index will NOT be scaled with beam energy, as this can have a nonlinear effect on the scattering factors f1 and f2. In practice, beam energy does not vary significantly during an experiment, so this should not affect the accuracy of simulation in most cases. Figure SEQ Figure \* ARABIC 3. A) Interface for intensity calculator tool on startup. B) First output window with X-axis set to “Incident angle” and Y axis set to “none”, with an incident angle range of 0.12±0.10° (0.02 to 0.22°) and other parameters as indicated in A. The dependence of scattering intensity on incident angle is shown with and without beam divergence taken into account. C) Second output window from the simulation in B, showing an image of the average beam footprint (right panel) and the dependence of footprint on incident angle (left panel). D) Final output window from the simulation in B, showing the overall dependence of normalization factor on incident angle (right panel) and the % change of scattering intensity from the nominal incident angle (left panel). The nominal value is shown with a black x.A) Calculating intensity versus a single variableTo display a 1D plot select “None” for the “X axis” dropdown and the variable of interest for the “Y axis” dropdown, or vice versa. In either case the variable will be plotted along the X axis and intensity will be plotted along the Y axis. Intensity will be calculated for values of the variable of interest specified by the ± box. For example, specifying an incident angle of 0.12 ± 0.1 will plot intensity for incident angle between 0.02° and 0.22°. Press “Go!” to launch the simulation.Each simulation will create two output windows, plus one additional figure if the incident angle is selected as a variable. The first window (when incident angle is selected), shown in REF _Ref9187859 \h Figure 3B, is a plot of absorbed intensity versus incident angle, with an ideal non-divergent beam and with divergence if specified. The next window, shown in REF _Ref9187859 \h Figure 3C, contains two panels: the first contains a top-down image of the beam footprint on the sample, and the second shows the fraction of beam which hits the sample versus the selected variable. The final window, shown in REF _Ref9187859 \h Figure 3D, also contains two panels: the first contains a plot of the normalization factor versus the selected variable and the second shows the % change in normalization factor compared to the nominal variable value.B) Calculating intensity versus two variablesThe change in intensity versus two variables can be visualized similarly to the procedure for 1 variable. Both an X axis and Y axis variable should be selected from the dropdown, and their ± boxes filled. This gives similar output figures to the 1-variable procedure but with a color-mapped intensity image. This type of image can help compare the impact of two parameters on scattering intensity. An example of the final output figure is given below.Figure SEQ Figure \* ARABIC 4. Output window with “X axis” set as incident angle for 0.12±0.10°, “Y axis” set as sample length for 20±5 mm, and other parameters as specified in REF _Ref9187859 \h Figure 3A.C) Entering a list of samples to extract a normalization factorA normalization factor calculator can be launched from the main intensity GUI by pushing the “Enter sample list…” button. This will pop up a new GUI figure that is formatted like a table, where each column represents a single sample, as shown in REF _Ref9190430 \h Figure 5. For this function, it is assumed that the beam dimensions and substrate are constant between samples, and these values are taken from the main GUI. The user defines an incident angle, beam energy, sample length/width/thickness, and β and δ for each sample (the top 7 rows). There is no row for changing sample density: this change is accounted for by updating β and δ values. Up to 50 samples at a time can be simulated.For convenience, if a parameter is the same for several samples, it only needs to be entered for the first sample. The parameter will be repeated to fill blank table cells. The first column is automatically filled in from the values in the main GUI.Once all samples are defined, the user should click once outside of the table and once on the “Go!” button to generate normalization values. Results are given in rows 9 and 10. The raw normalization factor is shown as the first result, followed by the normalization factor rescaled by the median value in the sample set. The second value is useful as a normalization factor because it is approximately 1, and using this normalization factor won’t significantly change the overall intensity. However, this value may change if the sample set is changed or expanded, whereas the raw scattering intensity will remain constant.Figure SEQ Figure \* ARABIC 5. Interface for “Enter sample list…” tool launched from the main scattering intensity GUI. Interface is shown after simulation has been run, which has populated the results in rows 9 and 10.D) Error bar analysis for a single sampleSystematic error due to changes in sample properties can be avoided by calculating an accurate normalization factor. However, uncertainty in input parameters will contribute error bars to each sample intensity.Simulating the expected error bar for a single sample uses the same main GUI, but instead of plotting scattering intensity versus variable(s), five variables are allowed to vary at once. These variables are incident angle and sample position, length, density, and thickness. Sample width is not considered as this is usually much larger than the beam width and shouldn’t affect scattering intensity in most cases. To access the error bar calculation mode, simply select “none” for both “X axis” and “Y axis” variables. An example setup is shown in REF _Ref9191376 \h Figure 6A.After launching the simulation by clicking “Go!”, a popup window will prompt the user for a time limit in seconds for the Monte Carlo calculation, as shown in REF _Ref9191376 \h Figure 6B. Since our method requires adding datapoints until convergence is reached, the calculation can become quite time-consuming on some computers—for example, a 10-minute simulation is not unexpected. However even a simulation that does not fully converge can give a ballpark error bar measurement which may be useful for quick comparisons. A progress window, shown in REF _Ref9191376 \h Figure 6C, will show the number of iterations completed and the fraction of the time limit that has passed.Once the time limit is reached or the convergence condition is satisfied, two windows give the result. The first window contains a graph of calculated error versus the number of iterations, and a histogram of intensity values for the final set of iterations, so that the user can confirm convergence. This figure additionally notes whether or not convergence was reached in the time limit and the final number of iterations. An example is shown in REF _Ref9191376 \h Figure 6E. The second window, shown in REF _Ref9191376 \h Figure 6D, is a simple message with the final % error, taken as the average error for the last 1000 points of the error versus # of iterations graph. In other words, if N total iterations are needed for convergence, the % error measured for N-1000 through N iterations is considered—this minimizes the effect of adding an “outlier” iteration. This final error is also shown as a horizontal line in the left panel of REF _Ref9191376 \h Figure 6E.Figure SEQ Figure \* ARABIC 6. A) Sample setup to launch an error bar calculation. B) Prompt to enter a cutoff time limit to avoid excessive computational time. C) Progress window showing the number of iterations completed and the fraction of the time liit that has passed. D) Window showing the final % error result of the calculation. E) Window showing the convergence of the error measurement as iterations are added (left panel) and a histogram of the final set of intensities (right panel).4. Scattering peak position analysisThe peak position analysis tool can be launched from R-refraction/RMain.m. The appearance of the interface on startup is shown in REF _Ref9194396 \h Figure 7A. The GUI is designed in a table format similar to the GUI for sample set scattering intensity analysis in Section 3C. Each peak is represented by a single row in the table. The calculation requires the user to specify a wavelength in Angstroms, an incident angle, and a film complex refractive index δ and β for each sample. For convenience, if a set of samples have the same parameter, that parameter only needs to be specified for the first sample and it will be replicated in all the blank rows below.The tool has two modes:Arrow pointing right ->: starting from a user-defined set of Qxy and Qz peak positions, calculate the actual d-spacings and polar angles.Arrow pointing left <-: starting from a user-defined set of d-spacings and polar angles, calculate the observed/measured Qxy and QZ.The two modes can be toggled by clicking on the arrow button at the top of the GUI. Both modes also calculate the 2θ, Qxy, and Qz values that would be observed without refraction. They also both calculate the position of the second peak due to reflection off of the substrate.After filling in the required parameters, click once outside the table and then click the “Go!” button to launch the simulation. The result is a plot with the peak positions in Qxy-Qz space, as shown in REF _Ref9194396 \h Figure 7C, and the table will be filled out, as shown in REF _Ref9194396 \h Figure 7B.The plot shows the positions of the directly scattered peak (circles) and reflected peak (dots) as observed on a detector during an actual measurement, as well as the Bragg’s Law position of the peak that would be observed with no refraction (x’s). Each peak is assigned a unique color to aid in identification.Figure SEQ Figure \* ARABIC 7. A) Interface for refractive shift calculator on startup. B) Interface after an example simulation. C) The result of B, showing the seven peak positions in Qxy-Qz space.5. Peak indexing toolThe peak indexing toolbox can be launched from I-indexing/IMain.m. The main interface contains buttons leading to each tool, as shown in REF _Ref9264551 \h Figure 8. Once crystal structure calculation is complete, it will also display the unit cell parameters in the space at the bottom of the window.Figure SEQ Figure \* ARABIC 8. Interface of the indexing tool on startup.A) Extracting peaksIn order to launch the peak finder tool, click the first button in the main GUI: “Find peak positions.” The peak finder can also be used as a standalone tool by launching I-indexing/peakfinder_waxtools.m. This will create the new window shown in REF _Ref9267002 \h Figure 9A. The GUI contains controls on the left side, a space for the Qxy-Qz image in the middle, and a table of extracted peak positions on the left side. Note, this tool requires an image which is already corrected to Qxy-Qz axes.The first step is to define the dimensions of the Qxy-Qz image. These are consistent with the dimensions given in WAXtools. The starting point is defined as the upper left corner and the ending point is the bottom right corner. All dimensions are in ?-1. For files that have not been generated by WAXtools, please see section 6C). Once dimensions have been defined, an image can be selected using the “Select data file…” button. This pops up a simple file selection window.The Qxy-Qz image will be displayed in the center panel. To select a peak, zoom in on the peak region using the zoom-in/zoom-out buttons at the top of the figure and clicking + dragging a box around the peak, as shown in REF _Ref9267002 \h Figure 9B. It may take longer to extract a peak from larger areas, so the peak fitting area should be kept as small as possible. Once the peak of interest has been isolated, “Fit peak on current region” will perform the fitting operation. This will create a new window, as shown in REF _Ref9267002 \h Figure 9C, with the rescaled peak overlaid with a fit topography map. The peak center is defined with a white x. The location of the peak center is also noted at the top of the figure.If the user is satisfied with the accuracy of this fit, the peak position can be added to the table on the left by hitting the “Add to table” button in the peak finder GUI. This also records the average FWHM calculated for the peak and the total peak intensity. The figure showing the fit result can be closed.Once all the peaks are fitted in this way, the table of peak positions can either be copied into a different location (select cells by clicking and dragging and use “Ctrl+C” to copy), or saved using the “Save peak positions…” button in the peak finder GUI. This will save the table data in a simple .txt document.Figure SEQ Figure \* ARABIC 9. A) Interface for peak finding tool on startup. B) Interface midway through peak fitting with several peak positions populating the table on the right. The next peak is selected by clicking and dragging with the magnifying tool around the peak area and clicking “Fit peak on current region.” C) Example of a fitted peak result. Original data is shown using a color scaling, fit result is shown with topography lines, and the fitted peak position is shown with a white x. The user can inspect the result before accepting the peak using the “Add to table” button.B) Uploading peaks from a fileWhen the peak finder tool is used, peak positions are automatically transferred to the main GUI. Otherwise, peak positions can be uploaded by selecting the “Upload peak positions” button in the main GUI. This uses a simple file selection window. Input should have at least two columns where the first column specifies Qxy values and the second column specifies Qz values. This can be in the form of a .txt file, as generated by the peak finder tool, or an .xlsx (Excel) file. The uploaded data will be shown in a new window so that the user can ensure accuracy of file reading. This new window can then be closed.C) Indexing peaks: “Step 1” and “Step 2”Indexing is performed in a two-step process. Both steps are opened by clicking the respective buttons in the main GUI, which both create a new window. In the following figures we index the peaks found in the Sample Files/SM1_peakpositions.txt file.The “Step 1” GUI has controls on the left side and a table with peak Qxy positions on the right side, as shown in REF _Ref9267921 \h Figure 10A. The table is automatically populated with the uploaded peak positions. The table also has columns for h and k values predicted for each peak, Qxy positions based on the prediction, error between measured and predicted Qxy position (ΔQxy), and calculated a||*, b||*, and γ||* (these last three are only displayed in the top row). These are populated after the simulation is run.At the top of the GUI, the user selects three peaks on which to base the initial calculation. These are numbered sequentially in the order they appear in the table. The first two peaks may have their (hk) indices restricted using the drop-down menus under “First peak restriction” and “Second peak restriction.” The first peak drop-down menu includes options “01”, “01 or 11”, “01, 11, 02”, and “01, 11, 02, 12.” The second peak drop-down menu includes these as well as a few more. These restrictions include all ± combinations of h and k. As an example, a restriction of “01 or 11” allows (hk) values of (01), (0-1), (10), (-10), (11), (1-1), (-11), (-1-1).Finally, the user can restrict the h and k values to be less than a certain value for the entire set of peaks. This can be used to keep indices as low as possible, which corresponds to keeping the unit cell as small as possible.With all restrictions in place, click “Calculate!” to launch the calculation. The program first finds the combination of (hk) values for the first three peaks which minimizes Qxy error. Then it runs two refinements to try and further minimize error through small modifications to the calculated a||*, b||*, and γ||* values. Note that even if (hk) was restricted for certain peaks for the initial calculation, the peaks may later be indexed with different (hk) if this minimizes overall error after refinement.The final peak errors are displayed in a new window as a function of Qxy position, as shown in REF _Ref9267921 \h Figure 10B, as well as being shown in the table.The “Step 2” GUI is similar to the “Step 1” GUI: input parameters on the left and a table of peak positions on the right, as shown in REF _Ref9267921 \h Figure 10C. In addition to Qxy and Qz values, the output table includes columns for h,k,? values of each peak, predicted Qz positions, Cartesian ΔQ error between measured and predicted peak positions, and a⊥*, b⊥*, and c* values (top row only). These are populated after calculation.Input parameters for “Step 2” include two peak selections for the initial calculation and an initial c* value. They also include restrictions on ? for the two selected peaks. A restriction of ?<=1 means that ? can be -1, 0, or 1. After calculation, a plot of Cartesian peak position error is displayed similarly to Step 1, as shown in REF _Ref9267921 \h Figure 10D.Figure SEQ Figure \* ARABIC 10. A) Interface for “Step 1,” set up to perform the first indexing step for the sample file SM1_peakpositions.txt. B) Result of running A: window showing ΔQxy versus Qxy for all peaks. C) Interface for “Step 2,” set up to perform the second indexing step. D) Result of running C: window showing ΔQ versus Qxy for all peaks.D) Viewing resultsOnce “Step 1” and “Step 2” are completed with satisfactory peak position error, the results can be viewed in several ways from the main GUI. Either the “View peaks” or “View unit cell” buttons will calculate the unit cell based on the indexing parameters and display results in the main GUI, as shown in REF _Ref9269655 \h Figure 11A.The “View peaks” button will create a window illustrating the measured peak positions (blue circles) and the calculated scattering pattern for the crystal structure (red dots), as shown in REF _Ref9269655 \h Figure 11B. It also displays the indices of each measured peak. Zooming in on a plot area can help resolve the overlapping text. In cases where multiple calculated scattering peaks could be responsible for a measured scattering peak, the indices of the closest one are displayed. Figure SEQ Figure \* ARABIC 11. A) Interface of indexing tool after indexing is completed, with unit cell parameters displayed. B) Result of “View peaks” tool, showing experimental and calculated peak positions. C) Result of “View unit cell” tool. The top left panel can be manipulated in 3D by using the “Rotate 3D” button in the top toolbox or under Tools→Rotate 3D. D) Result of “Modify peaks” tool. Checking the box in the top right corner allows changes in this window to propagate toward the IMain window (this overrides the original indexed unit cell parameters). E) Result of “Image overlay” tool, after loading an image file.The “View unit cell” button will pop out a figure illustrating the real-space unit cell, as shown in REF _Ref9269655 \h Figure 11C. At the top left is a 3D model of the unit cell which can be rotated by toggling the button at the top of the figure and clicking+dragging. The top right plot is an illustration of the unit cell as viewed from above the sample. The bottom two plots show the b-c and a-c planes, where the x-axis is along the substrate. These are viewed from a direction normal to the unit cell plane rather than a vertical cross-section of the film: i.e., they are not projections. Thus the y-axis is not necessarily the out-of-plane direction.The “Modify peaks” button launches a tool that can be used to modify unit cell parameters (a, b, c, α, β, γ) and observe how this changes the scattering peak pattern, compared to the measured peak positions, in real time. This window is shown in REF _Ref9269655 \h Figure 11D. Note that either “View peaks” or “View unit cell” should be used first to calculate the unit cell parameters. The unit cell can be modified on the left side with the sliders or by entering values directly into the bottom left boxes. The peak pattern is displayed on the right and updated continuously. Checking the “Override main GUI” box will cause the program to “forget” the unit cell calculated by indexing using “Step 1” and “Step 2” and instead use the new unit cell to calculate additional output figures. The “Modify peaks” tool can also be used independently without first calculating a unit cell or importing peak positions. Furthermore this tool can be launched separately from I-indexing/CellSlider.m (in this case, the “Override main GUI” toggle does nothing).The “Image Overlay” button is a hybrid of other features on the main interface. This GUI allows the user to open the original Qxy-Qz image (similar to “Find peak positions”), load peak positions if not already loaded (similar to “Upload Peak Positions”), overlay the peaks from the crystal structure (similar to “View Peaks”), and make modifications to the unit cell (similar to “Modify peaks”). It also gives a table indicating the closest possible simulated scattering pattern peak for each user-defined peak location (simulated Qxy, Qz, h, k, ?, and ΔQ between experimental and simulated peaks). This file can also be launched independently from I-indexing/IndexingImageOverlay_waxtools.m.Finally, a bonus file, I-indexing/EquivalentCell.m, may be useful for finding unit cell parameters for equivalent unit cells. This is a stand-alone script. For example, if the original unit cell has unit cell vectors (a1, a2, a3), the user can request the unit cell parameters (a, b, c,α, β, γ) for an equivalent cell with vectors (a1, a2+a3, a3). This is done by modifying the code at the top of the EquivalentCell.m file:% user inputs original unit cell values here% a,b,c in Angstrom% alpha,beta,gamma in degreesa=8.89147; alpha=117.579;b=36.5443; beta=65.6386;c=15.345; gamma=155.682; % user inputs three peaks to serve as new unit cell heren1='a1'; n2='a2+a3'; n3='a3';Running the code will create three windows. The first, shown in REF _Ref9270604 \h Figure 12, will show lattice points for 2x2x2 unit cells along with the original a1,a2,a3 vectors and the new n1,n2,n3 vectors. Each lattice point will be labeled with a position in terms of the original vectors. This plot can be rotated using the button. The second window contains a textbox with the new unit cell parameters. The third window shows the planes of the new unit cell equivalently to the “View Unit Cell” button (see REF _Ref9269655 \h Figure 11C).Figure SEQ Figure \* ARABIC 12. Result of running I-indexing/EquivalentCell.m. Blue circles show positions of unit cell vertices for 2x2x2 original unit cells. Blue lines show the original unit cell vectors, and red lines show the new unit cell vectors.6. Tips for navigating output data and customizing codeThis software package is completely open-source, which means that a user with some programming experience can customize parts of the code as needed. We include thorough commenting in our code to aid in this process. Here we list a few modifications that may be of interest.Line numbers given are approximate and subject to change as the code evolves. We suggest finding the code of interest by using the search function (“Edit->Find and Replace…” or Ctrl+F in the Editor window).A) General plot manipulation in MATLABi. Change plot appearanceThe appearance of any output plot can be modified using the MATLAB command window or the toolbars available in the figure window. These toolbars can be turned on from the “View->” drop-down menu.Figure Toolbar (appears above figure): simple zoom and rotation, add data label, toggle colorbars and legends, rearrange figure elements (turn on button), change text in labels by double-clickingCamera Toolbar (appears above figure): some useful functions for 3D figuresPlot Edit Toolbar (appears above figure): add annotations, arrows, and simple shapes; also modify text appearanceFigure Palette: add annotations, arrows, and simple shapesPlot Browser: show or hide datasetsProperty Editor (appears below figure): text font/size/alignment, line width/type/color, marker size/type/color, axis limits, colormapsMany of these options are also available from the dropdown menus.ii. Obtain data from a figureFigures can be copied to the clipboard using the “Edit->Copy Figure” dropdown menu. They can also be saved using “File->Save”, which saves as a MATLAB figure (.fig), or “File->Save As…” to save as a figure such as a bitmap or TIFF. The benefit to saving the file as a MATLAB figure is that its appearance can be modified and the raw data can be extracted later.Unfortunately, MATLAB does not contain a tool to easily extract raw data from a figure, and this must be done programmatically. To extract raw data from a 1D plot, first make the plot active by clicking on it. Then go to the MATLAB command window and type or paste the following commands:objs=get(gca,'Children');xdata=get(objs,'XData');ydata=get(objs,'YData');names= get(objs,'DisplayName');openvar('xdata')openvar('ydata')openvar('names')The “xdata” (x values) and “ydata” (corresponding y values) will appear as tabs in a Variable Editor window, which should resemble an Excel worksheet, as shown in REF _Ref9271078 \h Figure 13A. Similarly to Excel, data in cells can be highlighted (by clicking+dragging) and copied (Ctrl+C) to a new location for further processing. If there is only one curve plotted, the data will appear in a row in the Variable Editor. The data can be flipped from a row into a column if preferred by entering:xdata=xdata’;ydata=ydata’;On the other hand, if there are multiple curves plotted, “xdata” and “ydata” will instead contain text in each cell that says something like <1x101 double>. Each cell represents one curve. The legend label for each curve is given in the corresponding cell under the tab labeled “names.” Access the data for each curve by double clicking on its cell, which will open the data in a new tab labeled (for example) “xdata{1,1}”. The data can be switched from row to column by entering:xdata{1,1}=xdata{1,1}’;ydata{1,1}=ydata{1,1}’;where {1,1} should be replaced with the label on the data of interest.To extract raw data from a 2D image, first make the image active by clicking on it. In the MATLAB command window, type:objs=get(gca,'Children');xdata=get(objs,'XData');ydata=get(objs,'YData');Cdata=get(objs,'CData');openvar('xdata')openvar('ydata')openvar('Cdata')Here, “xdata” will contain a row specifying the x locations of the data points and “ydata” will specify the y locations of the data points (these can be turned into columns as before). “Cdata” will contain a 2D matrix giving the intensity throughout the image. This is shown in REF _Ref9271078 \h Figure 13B.Repeating the procedure for another figure will overwrite the data from the previous figure, unless new variable names are used.Figure SEQ Figure \* ARABIC 13. A) MATLAB main window showing commands and result for extracting data from a plot. B) Result for extracting data from a 2D image.iii. Obtain data from tableUnfortunately, tables in MATLAB GUIs do not always allow a user to copy and paste a portion of data. A workaround is to copy the table data into the command window, similarly to the previous section. To do this, click on the figure containing the table to make it active, and enter the following into the command window:t= findall(gcf,'type','uitable');dats=get(t,'Data');openvar('dats');This should open the table in the variable editor, where copy+paste works as expected.B) Modifying default values of GUI elementsFor each GUI, default input values have been set based on typical values. These can be modified permanently so that the user doesn’t have to re-enter them each time. The default values are usually set within the first ~100 lines of the relevant file. For example, in SI-scattering intensity/SIMain.m, some default parameter values are given near line 42:incang=uicontrol('Style','edit','String','0.12','Units','normalized','Position',[.26 .74 .1 .08]);ene=uicontrol('Style','edit','String','12.735','Units','normalized','Position',[.26 .64 .1 .08]);horiz=uicontrol('Style','edit','String','150','Units','normalized','Position',[.26 .54 .1 .08]);vert=uicontrol('Style','edit','String','50','Units','normalized','Position',[.26 .44 .1 .08]);div=uicontrol('Style','edit','String','0.014','Units','normalized','Position',[.26 .34 .1 .08]);New default values can be set by modifying the highlighted numbers.C) Useful modifications to scattering intensity analysisi. Changing the beam profile in spaceThe default scattering intensity GUI offers a Gaussian, Lorentzian, or rectangular lineshape, which is used to calculate the beam footprint. This may not be representative of an actual beam shape. To add a custom beam shape, first add this option to the dropdown menu. The code to do this is located in SIscattering intensity/SIMain.m near line 60. The original code is:shape=uicontrol('Style','popupmenu','String', ...{'Elliptical (Gaussian edge)','Elliptical (Lorentzian)','Elliptical (hard edge)','Rectangular (hard edge)'}, ...'Units','normalized','Position',[.26 .24 .24 .08]);Change this code by adding your own custom shape as a fourth option (any text can be used), without changing the original three options:shape=uicontrol('Style','popupmenu','String', ...{'Elliptical (Gaussian edge)','Elliptical (Lorentzian)','Elliptical (hard edge)','Rectangular (hard edge)',’Custom’}, ...'Units','normalized','Position',[.26 .24 .24 .08]);Next, the beam shape should be defined in the “footprintcalc” function. This is located near line 1330 of SIMain.m. Find the following code and add the highlighted code:if type==1 % elliptical, gaussian beam=gauss(subsX,0,FWHMhor / (2*sqrt(2*log(2))) ) * gauss(subsY,0,FWHMvert / (2*sqrt(2*log(2))) ); I=trapz(subsY,trapz(subsX,beam,1),2);elseif type==2 % elliptical, lorentzian...elseif type==5 % custom beam shape beam=[your beam footprint definition here]; I=trapz(subsY,trapz(subsX,beam,1),2);endThe beam footprint should be defined as a 2-dimensional array where the x axis is defined by subsX (perpendicular to beam direction) and the y axis is defined by subsY (parallel to beam direction). The variables subsX and subsY are defined earlier in the code such that 0,0 is the beam center and each axis is terminated at the edges of the sample. For example, if the sample is 10 mm wide, 20 mm long, and offset from the beam center by 3 mm, “subsx” spans from -5 mm to +5 mm and “subsy” spans from -7 mm to +13 mm. Furthermore, subsX is a column vector while subsY is a row vector. FWHMhor and FWHMvert record the horizontal and vertical beam widths after projection onto the substrate. Due to the projection geometry, FWHMvert runs along the x axis (along the beam direction) and FWHMhor runs along the y axis. The variable “I” integrates the total intensity that falls on the sample surface. To use a fractional beam measurement, “I” should be divided by the total beam intensity: this is not strictly necessary for calculation, but is useful to maintain consistency between experiments that may have different beam intensities.A custom beam shape can also be defined based on a function, which will closely follow the code for a Gaussian or Lorentzian beam (type 1 or type 2). In this code, a 2D beam is created by multiplying a vertical vector by a horizontal vector (an outer product). For example, a beam with a pyramid structure can be defined based on a new function “triangle”:function I = triangle( x,center,FWHM ) % Generate a triangle shape with a total area of 1. I=zeros(size(x)); for k=1:length(x) if x(k)<center-FWHM || x(k)>center+FWHM I(k)=0; elseif x(k)<center I(k)=1 + 1/FWHM .* (x(k)-center); elseif x(k)==center I(k)=1; else I(k)=1 - 1/FWHM .* (x(k)-center); end end I=I./FWHM;endThen, the new code to be added in SIMain.m will be:elseif type==5 % triangle beam= triangle(subsX,0,FWHMhor) * triangle(subsY,0,FWHMvert); I=trapz(subsY,trapz(subsX,beam,1),2);Otherwise, a beamshape can be defined numerically based on an image taken of the beam. For example, consider an array A that records beam intensity versus x (horizontal) and y (vertical), as measured perpendicular to the beam direction (a cross-section of the beam). The size of each pixel in x (in mm) is given as xpixel and the size of each pixel in y (in mm) is given as ypixel. The center of the beam is assumed to be at the center of the array. “incangle” is the incident angle in radians, which is available within the “footprintcalc” function. Then the following code will transform A into a beam footprint:A=###; % your 2D data herexpixel=#; % mm/pixelypixel=#; % mm/pixel % define x, y vectors of your data[sx,sy]=size(A);xdata=linspace(-xpixel*sx/2,xpixel*sx/2,sx)';ydata=linspace(-ypixel*sy/2,ypixel*sy/2,sy);% project onto substrateydata= ydata ./ sin(incangle); % turn x, y axis vectors into arraysxdata2=xdata*ones(size(ydata));ydata2=ones(size(xdata))*ydata;subsX2=subsX*ones(size(subsY));subsY2=ones(size(subsX))*subsY; % define beam: =0 outside provided databeam=zeros(size(subsX2));% interpolate your data where availablevalidpoints=subsX2>min(xdata) & subsX2<max(xdata) & ... subsY2>min(ydata) & subsY2<max(ydata);beam(validpoints)=interp2(xdata2',ydata2',A',subsX2(validpoints),subsY2(validpoints));I=trapz(subsY,trapz(subsX,beam,1),2);This code should create a beam profile that is compatible with the rest of the program.ii. Changing the beam divergence distributionThe divergence of the beam (or angular distribution) is assumed to have a Gaussian lineshape by default. This can be changed in SIMain.m in the “scattintensitycalc” function. The relevant code is near line 1264:lineshape=gauss(-5*G:stepsize:5*G,0, G ); % normalize by summed intensity to avoid changing the intensity of the% I(angle) curve when performing convolutionlineshape=lineshape/sum(lineshape);Here, “lineshape” represents a plot of beam intensity versus angle which is convolved with the plot of total absorbed intensity versus angle (a broadening effect). Because of this convolution, it is important that the angular step size between points is the same for both plots. This is defined in the variable “stepsize.” The variable G relates to the user-defined beam divergence through FWHM= G?2ln2= beam divergence. In the original code, the user-defined beam divergence is G, and the curve is calculated out to ±5*G for accuracy. The lineshape is also normalized by its sum (NOT integral) to avoid changing the magnitude of the absorbed intensity when the convolution is applied.The angular distribution can be changed by modifying these two lines, as long as the same “stepsize” is used. For example, “gauss” in the above code can be simply replaced with “lorentz” for a Lorentzian distribution, though in this case the variable L should be used such that FWHM= L?2=beam divergence.iii. Adding or modifying shortcut substrate refractive indicesTo add substrate refractive indices, first these options must be added to the dropdown menu. The code to do this is located in SIMain.m near line 96. The original code, and its modification (highlighted), is:dropsub=uicontrol('Style','popupmenu','Units','normalized','Position',[.54 .09 .2 .08],'String',... {'Select substrate...','Silicon','SiO2','Microscope slide','Si3N4','ITO (In2O3/SnO2)',’Custom’}, ... 'Callback',{@changesubs});Next, the refractive indices should be added to the “changesubs” function. This code is near line 191:switch S case 2; set(SBeta,'String',2.8244e-8); set(SDelta,'String',3.0066e-6); % Silicon case 3; set(SBeta,'String',1.7833e-8); set(SDelta,'String',3.4156e-6); % SiO2 case 4; set(SBeta,'String',2.0404e-8); set(SDelta,'String',3.2095e-6); % Microscope slide case 5; set(SBeta,'String',2.4638e-8); set(SDelta,'String',4.0855e-6); % SiN case 6; set(SBeta,'String',2.9854e-7); set(SDelta,'String',8.0313e-6); % ITO case 7; set(SBeta,'String',###); set(SDelta,'String',###); % Custom endThe appropriate substrate β (“SBeta”) and δ (“SDelta”) values should be added for the new substrate at the energy of interest, at the ### location.The default refractive indices for the existing substrates can be changed in this code as well. This is necessary if the energy is changed significantly from the beamline 11-3 standard of 12.735 keV.iv. Changing the convergence condition for error bar calculation or error bar definitionThe Monte Carlo calculation for error bars continues to add points until the user-defined time limit or the convergence condition is reached. Convergence conditions are defined in the SIMain.m, in the calc0Derror function, near line 950:% track the error versus number of samples that have been tested. Error% is measured as the standard deviation divided by the average x100%.errortrack(k+1)=std(Iall)/mean(Iall)*100; % start testing for convergence once at least 1000 points have been% consideredif length(errortrack)>=1000 % fit a line to the last 1000 points p=polyfit(1:1000,errortrack(end-999:end),1); % changeerror1=slope of the line fit to the last 1000 points. % "Convergence" means that the line should have flattened out and % slope is small. changeerror1=abs(p(1)); % changeerror3=standard deviation of the last 1000 points. % "Convergence" means that the error shouldn't be bouncing around % too much as more points are added. changeerror3=std(errortrack(end-999:end)); % errorbar=average value of the last 1000 points. Potentially more % accurate than the last value only, which could be affected by the % addition of an "outlier" sample. error1=mean(errortrack(end-999:end)); % if 1000 points haven't yet been reached, just update error to the% most recent valueelse error1=mean(errortrack);endThe “errortrack” variable contains an array of the measured error bar versus the number of iterations: e.g., with 100 datapoints the error bar is errortrack(100), with 200 datapoints it is errortrack(200), etc. Convergence is reached when adding more datapoints does not significantly change the error bar, which means the errortrack variable should be approximately constant. To test this, we check the slope of a line fit to the last 1000 points of errortrack (changeerror1) as well as the standard deviation of the last 1000 points (changeerror3). We require that the slope be less than 0.0001% (or a total change of 0.1% over the last 1000 points) and the standard deviation is less than 0.01%: this is specified at the top of the “while” loop near line 920 and again near line 998:while ( changeerror1>1e-4 || changeerror3>1e-2) && (toc<timelimit)…if changeerror1>1e-4 || changeerror3>1e-2; title('Warning: convergence not reached in time limit');A user can manually alter the convergence condition in several ways. Redefine the “changeerror” variable near line 950 (for example, from standard deviation to the difference between the last two points errortrack(end)-errortrack(end-1)):Define a new “changeerror4” convergence condition variable near line 950, and add this condition to the “while” loop statement near 920 and the “if” statement near line 998:(changeerror1>1e-4 || changeerror3>1e-2 || changeerror4>###)(### should contain additional convergence condition)Change the convergence condition in the “while” loop statement near 920 and the “if” statement near line 998 without changing the convergence variablesThe output error bar value is defined as the average value for the last 1000 points of errortrack, and errortrack is defined as the standard deviation / average value of all intensity points. This can be changed by redefining “errortrack” or “error1” near line 950 as needed. For example, to measure the variance instead of the standard deviation, change the “errortrack” definition to:errortrack(k+1)=var(Iall)/mean(Iall)*100;Or, to change the output error bar to the error calculated for the full set of iterations, change the “error1” definition to:error1=errortrack(end);v. Isolating scattering intensity from one layer of the filmIn some cases, a thin film has a different surface structure versus a bulk structure. For example, a 100 nm film may have a face-on bulk orientation but an edge-on orientation in the top 5 nm. Then the scattering intensity from the top portion will follow beam intensity integrated for the top of the film rather than the entire film. To examine scattering intensity from a specific region of the film, the user can manually change the integration region near line 1295 of SIMain.m:z=linspace(0,sampthick,1000);...intens=trapz(z,abs(EE).^2,2);This line performs an integration (trapz) of the scattering intensity (EE2) for all z, where z is defined as a vector from 0 to the sample thickness (sampthick) with 1000 points. Here, “sampthick” is given in Angstrom; z=0 is the top surface of the film; z=sampthick is the bottom surface of the film (substrate interface). To change the limits of integration, change the highlighted text accordingly. For example, to integrate over just the top 5 nm (50 Angstrom), this text will be:z=linspace(0,50,1000);To exclude the top 5 nm, the text will be:z=linspace(50,sampthick,1000);vi. Changing the number of points simulatedThe number of points for a scattering intensity versus 1 or 2 variables is defined at the beginning of the functions “calc2Derror” and “calc1Derror” respectively. By default, the program simulates 101x101 points for 2 variables or 301 points for 1 variable, but a user may manually request a larger or smaller parameter space. For 2D output (where both an X axis variable and a Y axis variable are selected), the code near line 203 should be modified:% define the sizes of the variables. Must be odd!szx=101;szy=101;Here, “szx” specifies the number of points simulated for the X axis and “szy” specifies the number of points fo the Y axis. For 1D output (where an X axis variable OR a Y axis variable are selected), the code near line 632 should be modified:% set up size of arraysz=301;“sz” is the number of points simulated. In order to include the nominal value, an odd number should be entered. While any number of points can be specified, adding points can significantly increase computational time, especially in the 2D case.D) Useful modifications to peak indexing tooli. Changing Qxy-Qz image file typeThe default image format in the peak finder tool (I-indexing/peakfinder_waxtools.m) and the image overlay tool (I-indexing/IndexingImageOverlay_waxtools.m) is a .tif image generated by WAXtools (using the “save as tif” button). We also include alternate files (I-indexing/peakfinder_wxdiff.m and IndexingImageOverlay_wxdiff.m) which can handle .bin images generated by the wxdiff software (using “File>Save binary image”). The GUIs for these versions have slightly modified inputs to reflect the parameters readily available in wxdiff, as shown in REF _Ref9360968 \h Figure 14.Figure SEQ Figure \* ARABIC 14. Interface of the wxdiff version of the Peak Finder tool on startup.Other than the modified inputs and a different file import protocol, the two versions are identical. The main GUI can be modified to use the wxdiff version by default by changing the following lines in Iindexing/IMain.m (near line 10):ppeaks=uicontrol('Style','pushbutton','String','Find peak positions','Units','normalized','Position',[.1 .9 .8 .08],'Callback',{@peakfinder_waxtools});…overlay=uicontrol('Style','pushbutton','String','Image overlay','Units','normalized','Position',[.52 .4 .38 .08],'Callback',{@IndexingImageOverlay_waxtools});Change peakfinder_waxtools to peakfinder_wxdiff, and IndexingImageOverlay_waxtools to IndexingImageOverlay_wxdiff.If software other than WAXtools or wxdiff is used to generate Qxy-Qz files, the file import procedure may need to be modified. First, try loading the file using both the WAXtools and the wxdiff version—one or the other may work for a new file type. If neither works, the file import script in the “getf” function of the files will need to be modified. In peakfinder_waxtools.m/IndexingImageOverlay_waxtools.m, the relevant lines (near line 65/99) uses the MATLAB “importdata” function:A=double(importdata([pname fname]))';In peakfinder_wxdiff.m/IndexingImageOverlay_wxdiff.m, the relevant lines (near line 76/110) use the MATLAB “fopen” and “fread” functions:fid=fopen([pname fname]);A=reshape(fread(fid,inf,'uint16'),length(qx),length(qy)); fclose(fid);Some trial and error may be needed to determine the appropriate import procedure for a file type. The writer of the software used to generate the image file may be able to assist.ii. Adding a polar angle “sin(χ)” intensity correction to Qxy-Qz imageThe polar angle correction can be turned on simply by uncommenting (removing %s, or highlighting and pressing Ctrl+T, or highlighting and selecting the “Text->Uncomment” dropdown menu) from the following code from each file (peakfinder_waxtools.m and IndexingImageOverlay_waxtools.m or, if using wxdiff to generate files, their _wxdiff.m equivalents):% chi=zeros(size(A));% for xx=1:length(qx)% for yy=1:length(qy)% chi(xx,yy)=atan(qx(xx)./qy(yy));% end% end% A=A.*abs(sin(chi));While this will make extracted peak intensities consistent with integration over 3D space, this may make it more difficult to distinguish peaks at polar angles near out-of-plane. Alternatively, the peak intensities recorded in the table may be corrected (without changing the way the image looks) by adding the following highlighted code in peakfinder_waxtools.m (or _wxdiff.m) near line 101 (or 110):dats=get(tab,'Data');pkint=pkint * cos( atan(pky/pkx) );dats=[dats; pkx pky pkwid pkint];set(tab,'Data',dats);peakpos=dats;iii. Adding a “weight” parameter to peak indexing position errorWhen indexing a scattering pattern, it may be more important to have an accurate peak position for certain peaks than others. For example, a weak and broad peak at high Q may be less important than a well-defined and high intensity peak at low Q. By default, our code weighs all peaks equally, but a user may want to place more emphasis on some peaks than others. The peaks may be weighted by their widths, their intensity, or some other user-defined parameter.A weight can be added in both I-indexing/Step1.m and I-indexing/Step2.m. The easiest way to do this is to add the weight parameter to the peak position file as a new column, if it does not already exist. This new column will be imported into the “peakpos” variable. Peaks generated and saved by the peak finder GUI have their average FWHM and intensity saved as columns 3 and 4 respectively in the peakpos variable, which is defined to be “global” (can be accessed from within any function in the same file). In Step1.m, error can be weighted in two places: during the initial selection of (hk) for the three user-selected peaks, and during the refinement process. This code is found near line 169 (“calc” function) and line 265 (“refine” function) respectively:err(xx)=err(xx) + min(abs(Qxy(n).^2-qguess));…for n=1:length(a) % for this a,b,gam combination, add up error between each % experimental Qxy and calculated Qxy for its assigned (hk) % err=sum(error^2) err(n)=sum( ( Qxy - sqrt(a(n)^2*h.^2+b(n)^2*k.^2 + 2*a(n)*b(n).*h.*k.*cosd(gam(n))) ).^2 ); % err=sqrt( sum(error^2)/numpoints ) = RMS error err(n)=sqrt(err(n)/length(Qxy)); % update waitbar waitbar(n/length(a));endAt the first location, a weight parameter can be added simply by multiplying by the calculated error for the nth peak:err(xx)=err(xx) + min(abs(Qxy(n).^2-qguess)) * peakpos(n,3);Here, “3” specifies the third column in peakpos, but any other column or some arithmetic combination of columns can be used as desired.At the second location, a weight parameter can be added as follows (here, the “n” index no longer refers to peak number):err(n)=sum( ( Qxy - sqrt(a(n)^2*h.^2+b(n)^2*k.^2 + 2*a(n)*b(n).*h.*k.*cosd(gam(n))) ).^2 .* peakpos(:,3));In Step2.m, a similar procedure is followed. The relevant lines are near 126 and 226:err(xx)=err(xx)+min( sqrt( (Qxy(n)-Qxyguess).^2 + (Qz(n)-Qzguess).^2 ) ) * peakpos(n,3) );…for m=1:length(Qxy) % error=minimized Cartesian distance % err=sum(error^2) err(n)=err(n) + min( sqrt( (Qxy(m)-Qxyguess).^2 + (Qz(m)-Qzguess).^2 ) .* peakpos(m,3) ) ^2 ;endiv. Changing refinement parameters for “Step 1” and “Step 2”Both Step1.m and Step2.m refine variables by testing all combinations within certain boundaries. The refinement function inputs are:For Step1.m: refine(dats,arange,asteps,brange,bsteps,gamrange,gamsteps)“dats” is the array shown in the table, which includes the currently calculated values for a||*, b||*, and γ||*. “arange" and “asteps” indicate that a||* will be sampled in the range of the current value ± arange in “asteps” number of steps. b||* and γ||* follow a similar format with “brange/bsteps” and “gamrange/gamsteps” respectively.For Step2.m: refine(dats,hguess,kguess,lguess,Qxyguess,arange,asteps,brange,bsteps,crange,csteps)“dats” is again the array in the table. hguess, kguess, and lguess are lists of possible values for h, k, and ?, and Qxyguess contains the corresponding Qxy peak positions calculated for the a||*, b||*, and γ||* defined Step1.m (which are assumed to be accurate). Similarly to the function in Step1.m, “arange" and “asteps” define the sampling range and number of steps for a⊥*; and the same for b⊥* and c* (“brange/bsteps” and “crange/csteps” respectively).The default refinements for Step1.m are given near line 213 of that file:% refine a, b, gamma to minimize error. Each refinement updates the% data table, but (hk) values remain constant. % Each refinement pops up its own waitbar.dats=refine(dats,.01,40,.01,40,1,40, 1,2);dats=refine(dats,.001,40,.001,40,.1,40, 2,2);Two refinements are performed. The first refinement varies a||* and b||* within ±0.01 ?-1 in 40 steps and varies γ||* within ±1° in 40 steps. The second refinement uses ±0.001 ?-1 and ±0.1°. The user can modify this refinement process by adding or taking away refinement steps or by changing either the ranges or step sizes for the refinement steps. The computational time can increase significantly with more steps: it scales with asteps·bsteps·gamsteps.The default refinements for Step2.m are given near line 104 of that file:% refine values% refine apar, bpar, c* to minimize error. Each refinement updates the% data table. Each pops up its own waitbar.dats=refine(dats,hguess,kguess,lguess,Qxyguess,.02,20,.02,20,.1,20, 1,3);dats=refine(dats,hguess,kguess,lguess,Qxyguess,.005,20,.005,20,.05,20, 2,3);dats=refine(dats,hguess,kguess,lguess,Qxyguess,.001,20,.001,20,.01,20, 3,3);Three refinements are performed. The first varies a⊥* and b⊥* within ±0.02 ?-1 in 20 steps and c* within ±0.1 ?-1 in 20 steps. The second uses ±0.005 ?-1 for a⊥* and b⊥* and ±0.05 ?-1 for c*. The third uses ±0.001 ?-1 for a⊥* and b⊥* and ±0.01 ?-1 for c*. Again, the user can add or take away refinement steps or change the ranges/number of steps. Again, computational time scales as asteps·bsteps·csteps. ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download