source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/FIT_Ops.ipf @ 570

Revision 570, 37.3 KB checked in by srkline, 4 years ago (diff)

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5//**************************
6//procedures for creating and initializing the FIT panel
7//global variables (numerical only) are kept in root:myGlobals:FIT folder,
8//created as needed
9//
10//also contains code for the FIT_RPA panel and its fitting procedures
11//
12//**************************
13//
14
15//main procedures to open the panel, initializing the data folder and global variables
16//as necessary. All are kept in a :FIT subfolder to avoid overlap with other variables
17//
18// To use any of the fit functions in the FIT panel, I(Q) data must already
19// be loaded into memory (using "plot" from the 1-D data operations
20// ** this may be useful to change in the future, replacing the 3 popups
21// with a list box - allowing the user to pick/load the data from the fit panel
22// and not offering any choice of q/i/s waves to use. (more consistent with the operation
23// of the FIT/RPA panel)
24//
25Proc OpenFitPanel()
26        If(WinType("FitPanel") == 0)
27                //create the necessary data folder
28                NewDataFolder/O root:myGlobals:FIT
29                //initialize the values
30                Variable/G root:myGlobals:FIT:gLolim = 0.02
31                Variable/G root:myGlobals:FIT:gUplim = 0.04
32                Variable/G root:myGlobals:FIT:gExpA = 1
33                Variable/G root:myGlobals:FIT:gExpB = 1
34                Variable/G root:myGlobals:FIT:gExpC = 1
35                Variable/G root:myGlobals:FIT:gBack = 0
36                String/G root:myGlobals:FIT:gDataPopList = "none"
37                FitPanel()
38        else
39                //window already exists, just bring to front for update
40                DoWindow/F FitPanel
41                CheckBox check0,value=0         //deselect the checkbox to use cursors
42        endif
43        //pop the file menu
44        FIT_FilePopMenuProc("",1,"")
45End
46
47//the actual window recreation macro to draw the fit panel. Globals and data folder must
48// already be initialized
49Window FitPanel()
50        String angst = root:myGlobals:gAngstStr
51        PauseUpdate; Silent 1           // building window...
52        NewPanel /W=(461,46,735,455)/K=1
53        ModifyPanel cbRGB=(32768,54615,65535), fixedSize=1
54        SetDrawLayer UserBack
55        DrawText 56,20,"Select Experimental Data"
56        DrawText 66,138,"q-range to fit ("+angst+"^-1)"
57        DrawText 42,268,"Select the y and x-axis scaling"
58        DrawLine 1,21,271,21
59        DrawLine -1,272,273,272
60        DrawLine -1,140,272,140
61//      PopupMenu xwave,pos={13,31},size={158,19},title="q-values"
62//      PopupMenu xwave,help={"Select the experimental q-values. If no data is available, use Macro LoadQISData or LoadQSIGData."}
63//      PopupMenu xwave,mode=1,value= #"WaveList(\"*q\",\";\",\"\")"
64        PopupMenu ywave,pos={13,40},size={154,19},title="Data File"
65        PopupMenu ywave,help={"Select the experimental intensity values"}
66        PopupMenu ywave,mode=1,value=root:myGlobals:FIT:gDataPopList,proc=FIT_FilePopMenuProc
67        Button loadButton,pos={25,80},size={130,20},proc=FIT_Load_Proc,title="Load and Plot File"
68        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."}
69        Button helpButton,pos={200,80},size={25,20},proc=showFITHelp,title="?"
70        Button helpButton,help={"Show help file for linearized fitting"}
71//      PopupMenu ewave,pos={13,87},size={186,19},title="Std. Deviation"
72//      PopupMenu ewave,help={"Select the standard deviation of the intensity"}
73//      PopupMenu ewave,mode=1,value= #"WaveList(\"*s\",\";\",\"\")"
74        PopupMenu ymodel,pos={20,281},size={76,19},title="y-axis"
75        PopupMenu ymodel,help={"This popup selects how the y-axis will be linearized based on the chosen data"}
76        PopupMenu ymodel,mode=1,value= #"\"I;log(I);ln(I);1/I;I^a;Iq^a;I^a q^b;1/sqrt(I);ln(Iq);ln(Iq^2)\""
77        Button GoFit,pos={60,367},size={70,20},proc=DispatchModel,title="Do the Fit"
78        Button GoFit,help={"This button will do the specified fit using the selections in this panel"}
79        Button DoneButton,pos={180,367},size={50,20},proc=FITDoneButton,title="Done"
80        Button DoneButton,help={"This button will close the panel and the associated graph"}
81        SetVariable lolim,pos={64,147},size={134,17},title="Lower Limit"
82        SetVariable lolim,help={"Enter the lower q-limit to perform the fit ("+angst+"^-1)"}
83        SetVariable lolim,limits={0,5,0},value= root:myGlobals:FIT:gLolim
84        SetVariable uplim,pos={63,169},size={134,17},title="Upper Limit"
85        SetVariable uplim,help={"Enter the upper q-limit to perform the fit ("+angst+"^-1)"}
86        SetVariable uplim,limits={0,5,0},value= root:myGlobals:FIT:gUplim
87        SetVariable expa,pos={13,311},size={80,17},title="pow \"a\""
88        SetVariable expa,help={"This sets the exponent \"a\" for some y-axis formats. The value is ignored if the model does not use an adjustable exponent"}
89        SetVariable expa,limits={-2,10,0},value= root:myGlobals:FIT:gExpA
90        SetVariable expb,pos={98,311},size={80,17},title="pow \"b\""
91        SetVariable expb,help={"This sets the exponent \"b\" for some x-axis formats. The value is ignored if the model does not use an adjustable exponent"}
92        SetVariable expb,limits={0,10,0},value= root:myGlobals:FIT:gExpB
93        PopupMenu xmodel,pos={155,280},size={79,19},title="x-axis"
94        PopupMenu xmodel,help={"This popup selects how the x-axis will be linearized given the chosen data"}
95        PopupMenu xmodel,mode=1,value= #"\"q;log(q);q^2;q^c\""
96        CheckBox check0,pos={18,223},size={240,20},title="Use cursor range from FitWindow"
97        CheckBox check0,help={"Checking this will perform a fit between the cursors on the graph in FitWindow and ignore the numerical limits typed above"},value=0
98        SetVariable back,pos={70,338},size={139,17},title="background"
99        SetVariable back,help={"This constant background value will be subtracted from the experimental intensity before fitting is done"}
100        SetVariable back,limits={-Inf,Inf,0},value= root:myGlobals:FIT:gBack
101        SetVariable expc,pos={182,310},size={80,17},title="pow \"c\""
102        SetVariable expc,help={"This sets the exponent \"c\" for some x-axis formats. The value is ignored if the model does not use \"c\" as an adjustable exponent"}
103        SetVariable expc,limits={-10,10,0},value= root:myGlobals:FIT:gExpC
104        Button sh_all,pos={65,193},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range"
105        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."}
106EndMacro
107
108
109Proc FITDoneButton(ctrlName): ButtonControl
110        String ctrlName
111        DoWindow/K FitWindow
112        DoWindow/K FitPanel
113end
114
115Proc showFITHelp(ctrlName): ButtonControl
116        String ctrlName
117        DisplayHelpTopic/K=1/Z "SANS Data Reduction Tutorial[Fit Lines to Your Data]"
118        if(V_flag !=0)
119                DoAlert 0,"The SANS Data Reduction Tutorial Help file could not be found"
120        endif
121end
122
123//Loads the selected file for fitting
124//graphs the data as needed
125Proc FIT_Load_Proc(ctrlName): ButtonControl
126        String ctrlName
127       
128        //Load the data
129        String tempName="",partialName=""
130        Variable err
131        ControlInfo $"ywave"
132        //find the file from the partial filename
133        If( (cmpstr(S_value,"")==0) || (cmpstr(S_value,"none")==0) )
134                //null selection, or "none" from any popup
135                Abort "no file selected in popup menu"
136        else
137                //selection not null
138                partialName = S_value
139                //Print partialName
140        Endif
141        //get a valid file based on this partialName and catPathName
142        tempName = FindValidFilename(partialName)
143
144        //prepend path to tempName for read routine
145        PathInfo catPathName
146       
147        tempName = S_path + tempName
148       
149        //load in the data (into the its own folder)
150        A_LoadOneDDataWithName(tempName,0)
151        //Print S_fileName
152        //Print tempName
153       
154        String cleanLastFileName = CleanupName(root:Packages:NIST:gLastFileName,0)
155
156        tempName="root:"+cleanLastFileName+":"+cleanLastFileName
157        Duplicate/O $(tempName+"_q") xAxisWave
158        //tempName=cleanLastFileName+"_i"
159        Duplicate/O $(tempName+"_i") yAxisWave
160        //tempName=cleanLastFileName+"_s"
161        Duplicate/O $(tempName+"_s") yErrWave
162
163        //Plot, and adjust the scaling to match the axis scaling set by the popups
164        Rescale_Data()
165       
166End
167
168//gets a valid file list (simply not the files with ".SAn" in the name)
169//
170Function FIT_FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
171        String ctrlName
172        Variable popNum
173        String popStr
174       
175        String tempStr=ReducedDataFileList(ctrlName)
176        if(strlen(tempStr)==0)
177                tempStr = "Pick the data path"
178        Endif
179        String/G root:myGlobals:FIT:gDataPopList =tempStr               //function is in NSORT.ipf
180        ControlUpdate ywave
181       
182End
183
184//direct porting of the fit program from the VAX, with no corrections and only
185//minor modifications of additional  linearizations and the option
186//to subtract a constant (q-independent) background value before doing any
187//of the fits. The original data on disk (and as loaded) is never modified, all
188//manipulation is done from a copy of the data.
189
190//button procedure to show the entire axis range of the data, rather than just
191//the fitted range, which is the default display after a fit is performed
192//
193Function ShowAllButtonProc(ctrlName) : ButtonControl
194        String ctrlName
195
196        //bring the FitWindow to the front and Autoscale the axes
197        DoWindow/F FitWindow
198        SetAxis/A
199End
200
201// function that takes the current dataset (already loaded)
202// and replots it based on the X/Y axis scaling selected in the popups
203// (does not fit the data)
204//
205Function Rescale_Data()
206       
207        //Scaling exponents and background value
208        Variable pow_a,pow_b,pow_c,bkg
209        ControlInfo/W=FitPanel expa
210        pow_a = V_value
211        ControlInfo/W=FitPanel expb
212        pow_b = V_value
213        ControlInfo/W=FitPanel expc
214        pow_c = V_value
215        ControlInfo/W=FitPanel back
216        bkg = V_value
217       
218//check for physical limits on exponent values
219// if bad values found, alert, and reset to good values so the rescaling can continue
220        NVAR gA = root:myGlobals:FIT:gExpA
221        NVAR gB = root:myGlobals:FIT:gExpB
222        NVAR gC = root:myGlobals:FIT:gExpC
223        if((pow_a < -2) || (pow_a > 10))
224                DoAlert 0,"Exponent a must be in the range (-2,10) - the exponent a has been reset to 1"
225                gA = 1
226        endif
227        if((pow_b < 0) || (pow_b > 10))
228                DoAlert 0,"Exponent b must be in the range (0,10) - the exponent b has been reset to 1"
229                gB = 1
230        endif
231        //if q^c is the x-scaling, c must be be within limits and also non-zero
232        ControlInfo/W=FitPanel xModel
233        If (cmpstr("q^c",S_Value) == 0)
234                if(pow_c == 0)
235                        DoAlert 0,"Exponent c must be non-zero, c has been reset to 1"
236                        gC = 1
237                endif
238                if((pow_c < -10) || (pow_c > 10))
239                        DoAlert 0,"Exponent c must be in the range (-10,10), c has been reset to 1"
240                        gC = 1
241                endif
242        endif
243       
244        //do the rescaling of the data
245        // get the current experimental q, I, and std dev. waves (as they would be loaded )
246
247        //ControlInfo/W=FitPanel ywave
248        //get the filename from the global as it's loaded, rather from the popup - as version numbers
249        // do cause problems here. This global is also used later in this function
250        SVAR gLastFileName = root:Packages:NIST:gLastFileName
251        String tempStr = CleanupName(gLastFileName,0)
252        String tempName = "root:"+tempStr+":"+tempStr
253       
254        Wave xw = $(tempName+"_q")
255        Wave yw = $(tempName+"_i")
256        Wave ew = $(tempName+"_s")
257       
258        //variables set for each model to control look of graph
259        Variable xlow,xhigh,ylow,yhigh,yes_cursors
260        String xlabel,ylabel,xstr,ystr
261        //check for proper y-scaling selection, make the necessary waves
262        ControlInfo/W=FitPanel yModel
263        ystr = S_Value
264//      print "ystr = ",ystr
265        do
266                // make the new yaxis waves, including weighting wave
267                Duplicate/O yw yAxisWave,yErrWave,yWtWave,residWave
268                //subtract the background value from yAxisWave before doing any rescaling
269                yAxisWave = yw - bkg
270               
271                If (cmpstr("I",S_Value) == 0)
272                        SetScale d 0,0,"1/cm",yAxisWave
273                        yErrWave = ew
274                        yWtWave = 1/yErrWave
275                        yAxisWave = yAxisWave
276                        ylabel = "I(q)"
277                        break   
278                endif
279                If (cmpstr("ln(I)",S_Value) == 0)
280                        SetScale d 0,0,"",yAxisWave
281                        yErrWave = ew/yAxisWave
282                        yWtWave = 1/yErrWave
283                        yAxisWave = ln(yAxisWave)
284                        ylabel = "ln(I)"
285                        break   
286                endif
287                If (cmpstr("log(I)",S_Value) == 0)
288                        SetScale d 0,0,"",yAxisWave
289                        yErrWave = ew/(2.30*yAxisWave)
290                        yWtWave = 1/yErrWave
291                        yAxisWave = log(yAxisWave)
292                        ylabel = "log(I)"
293                        break   
294                endif
295                If (cmpstr("1/I",S_Value) == 0)
296                        SetScale d 0,0,"",yAxisWave
297                        yErrWave = ew/yAxisWave^2
298                        yWtWave = 1/yErrWave
299                        yAxisWave = 1/yAxisWave
300                        ylabel = "1/I"
301                        break
302                endif
303                If (cmpstr("I^a",S_Value) == 0)
304                        SetScale d 0,0,"",yAxisWave
305                        yErrWave = ew*abs(pow_a*(yAxisWave^(pow_a-1)))
306                        yWtWave = 1/yErrWave
307                        yAxisWave = yAxisWave^pow_a
308                        ylabel = "I^"+num2str(pow_a)
309                        break
310                endif
311                If (cmpstr("Iq^a",S_Value) == 0)
312                        SetScale d 0,0,"",yAxisWave
313                        yErrWave = ew*xw^pow_a
314                        yWtWave = 1/yErrWave
315                        yAxisWave = yAxisWave*xw^pow_a
316                        ylabel = "I*q^"+num2str(pow_a)
317                        break
318                endif
319                If (cmpstr("I^a q^b",S_Value) == 0)
320                        SetScale d 0,0,"",yAxisWave
321                        yErrWave = ew*abs(pow_a*(yAxisWave^(pow_a-1)))*xw^pow_b
322                        yWtWave = 1/yErrWave
323                        yAxisWave = yAxisWave^pow_a*xw^pow_b
324                        ylabel = "I^" + num2str(pow_a) + "q^"+num2str(pow_b)
325                        break
326                endif
327                If (cmpstr("1/sqrt(I)",S_Value) == 0)
328                        SetScale d 0,0,"",yAxisWave
329                        yErrWave = 0.5*ew*yAxisWave^(-1.5)
330                        yWtWave = 1/yErrWave
331                        yAxisWave = 1/sqrt(yAxisWave)
332                        ylabel = "1/sqrt(I)"
333                        break
334                endif
335                If (cmpstr("ln(Iq)",S_Value) == 0)
336                        SetScale d 0,0,"",yAxisWave
337                        yErrWave =ew/yAxisWave
338                        yWtWave = 1/yErrWave
339                        yAxisWave = ln(xw*yAxisWave)
340                        ylabel = "ln(q*I)"
341                        break
342                endif
343                If (cmpstr("ln(Iq^2)",S_Value) == 0)
344                        SetScale d 0,0,"",yAxisWave
345                        yErrWave = ew/yAxisWave
346                        yWtWave = 1/yErrWave
347                        yAxisWave = ln(xw*xw*yAxisWave)
348                        ylabel = "ln(I*q^2)"
349                        break
350                endif
351                //more ifs for each case
352               
353                // if selection not found, abort
354                DoAlert 0,"Y-axis scaling incorrect. Aborting"
355                Abort
356        while(0)        //end of "case" statement for y-axis scaling
357       
358       
359        //check for proper x-scaling selection
360        Variable low,high
361       
362        ControlInfo/W=FitPanel lolim
363        low = V_value
364        ControlInfo/W=FitPanel uplim
365        high = V_value
366        if ((high<low) || (high==low))
367                DoAlert 0,"Unphysical fitting limits - re-enter better values"
368                Abort
369        endif
370       
371        SVAR aStr = root:myGlobals:gAngstStr
372       
373        ControlInfo/W=FitPanel xModel
374        xstr = S_Value
375        do
376                // make the new yaxis wave
377                Duplicate/o xw xAxisWave
378                If (cmpstr("q",S_Value) == 0)   
379                        SetScale d 0,0,aStr+"^-1",xAxisWave
380                        xAxisWave = xw
381                        xlabel = "q"
382                        xlow = low
383                        xhigh = high
384                        break   
385                endif
386                If (cmpstr("q^2",S_Value) == 0)
387                        SetScale d 0,0,aStr+"^-2",xAxisWave
388                        xAxisWave = xw*xw
389                        xlabel = "q^2"
390                        xlow = low^2
391                        xhigh = high^2
392                        break   
393                endif
394                If (cmpstr("log(q)",S_Value) == 0)     
395                        SetScale d 0,0,"",xAxisWave
396                        xAxisWave = log(xw)
397                        xlabel = "log(q)"
398                        xlow = log(low)
399                        xhigh = log(high)
400                        break   
401                endif
402                If (cmpstr("q^c",S_Value) == 0)
403                        SetScale d 0,0,"",xAxisWave
404                        xAxisWave = xw^pow_c
405                        xlabel = "q^"+num2str(pow_c)
406                        xlow = low^pow_c
407                        xhigh = high^pow_c
408                        break
409                endif
410       
411                //more ifs for each case
412               
413                // if selection not found, abort
414                DoAlert 0,"X-axis scaling incorrect. Aborting"
415                Abort
416        while(0)        //end of "case" statement for x-axis scaling
417
418        //plot the data
419       
420        String cleanLastFileName = "root:"+CleanupName(gLastFileName,0)
421        If(WinType("FitWindow") == 0)
422                Display /W=(5,42,480,400)/K=1 yAxisWave vs xAxisWave
423                ModifyGraph mode=3,standoff=0,marker=8
424                ErrorBars/T=0 yAxisWave Y,wave=(yErrWave,yErrWave)
425                ModifyGraph opaque(yAxisWave)=1
426                DoWindow/C FitWindow
427        else
428                //window already exists, just bring to front for update
429                DoWindow/F FitWindow
430                // remove old text boxes
431                TextBox/K/N=text_1
432                TextBox/K/N=text_2
433                TextBox/K/N=text_3
434        endif
435        SetAxis/A
436        ModifyGraph tickUnit=1          //suppress tick units in labels
437        TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
438        //clear the old fit from the window, if it exists
439        RemoveFromGraph/W=FitWindow/Z fit_yAxisWave
440       
441        // add the cursors if desired...       
442        //see if the user wants to use the data specified by the cursors - else use numerical values
443       
444        ControlInfo/W=FitPanel check0           //V_value = 1 if it is checked, meaning yes, use cursors
445        yes_cursors = V_value
446
447        DoWindow/F FitWindow
448        ShowInfo
449        if(yes_cursors)
450                xlow = xAxisWave[xcsr(A)]
451                xhigh = xAxisWave[xcsr(B)]
452                if(xlow > xhigh)
453                        xhigh = xlow
454                        xlow = xAxisWave[xcsr(B)]
455                endif
456//              Print xlow,xhigh
457        else
458                FindLevel/P/Q xAxisWave, xlow
459                if(V_flag == 1)                 //level NOT found
460                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
461                        //DoWindow/K FitWindow
462                        //Abort
463                endif
464                Cursor/P A, yAxisWave,trunc(V_LevelX)+1
465                ylow = V_LevelX
466                FindLevel/P/Q xAxisWave, xhigh
467                if(V_flag == 1)
468                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
469                        //DoWindow/K FitWindow
470                        //Abort
471                endif
472                Cursor/P B, yAxisWave,trunc(V_LevelX)
473                yhigh = V_LevelX
474        endif   //if(V_value)
475        //SetAxis bottom,xlow,xhigh
476        //SetAxis left,ylow,yhigh
477        Label left ylabel
478        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead       
479
480End
481
482
483//button procedure that is activated to "DotheFit"
484//the panel is parsed for proper fitting limits
485// the appropriate linearization is formed (in the Rescale_Data() function)
486// and the fit is done,
487//and the results are plotted
488// function works in root level data folder (where the loaded 1-d data will be)
489Function DispatchModel(GoFit) : ButtonControl
490        String GoFit
491
492        //check for the FitWindow - to make sure that there is data to fit
493        If(WinType("FitWindow") == 0)           //if the window doesn't exist
494                Abort "You must Load and Plot a File before fitting the data"
495        endif
496        // rescale the data, to make sure it's as selected on the panel
497        Rescale_Data()
498       
499        // now go do the fit
500       
501// get the current low and high q values for fitting
502        Variable low,high
503       
504        ControlInfo/W=FitPanel lolim
505        low = V_value
506        ControlInfo/W=FitPanel uplim
507        high = V_value
508        if ((high<low) || (high==low))
509                DoAlert 0,"Unphysical fitting limits - re-enter better values"
510                Abort
511        endif
512
513        //try including residuals on the graph /R=residWave, explicitly place on new axis
514        //if only /R used, residuals are automatically placed on graph
515       
516        CurveFit line yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D 
517        //CurveFit line yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave  /R /D 
518        ModifyGraph rgb(fit_yAxisWave)=(0,0,0)
519        ModifyGraph lsize(fit_yAxisWave)=2
520// annotate graph, filtering out special cases of Guinier fits
521// Text Boxes must be used, since ControlBars on graphs DON'T print out
522       
523        // need access to Global wave, result of fit
524        //ystr and xstr are the axis strings - filter with a do-loop
525        String ystr="",xstr=""
526        SVAR gLastFileName = root:Packages:NIST:gLastFileName
527        SVAR aStr = root:myGlobals:gAngstStr
528        String tmpStr = CleanupName(gLastFileName,0)
529        //ControlInfo/W=FitPanel ywave
530        Wave xw = $("root:"+tmpStr+":"+tmpStr+"_q")
531        ControlInfo/W=FitPanel yModel
532        ystr = S_Value
533        ControlInfo/W=FitPanel xModel
534        xstr = S_Value
535       
536        WAVE W_coef=W_coef
537        WAVE W_sigma=W_sigma
538        String textstr_1,textstr_2,textstr_3 = ""
539        Variable rg,rgerr,minfit,maxfit
540       
541        textstr_1 = "Slope = " + num2str(W_coef[1]) + " ± " + num2str(W_sigma[1])
542        textstr_1 += "\rIntercept = " + num2str(W_coef[0]) + " ± " + num2str(W_sigma[0])
543        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3))
544       
545        minfit = xw[xcsr(A)]
546        maxfit = xw[xcsr(B)]
547        textstr_2 = "Qmin =  " + num2str(minfit)
548        textstr_2 += "\rQmax =  " + num2str(maxfit)
549       
550        //model-specific calculations - I(0), Rg, etc.
551        //put these in textstr_3, at bottom
552        do
553                If (cmpstr("I",ystr) == 0)
554                        textstr_3 = "I(q=0) =  "  + num2str(W_coef[0])
555                        break   
556                endif
557                If (cmpstr("ln(I)",ystr) == 0)
558                        textstr_3 = "I(q=0) =  "  + num2str(exp(W_coef[0]))
559                        if(cmpstr("q^2",xstr) == 0)     //then a Guinier plot for a sphere (3-d)
560                                rg = sqrt(-3*W_coef[1])
561                                rgerr = 3*W_sigma[1]/(2*rg)
562                                textstr_3 += "\rRg ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
563                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
564                                break
565                        endif
566                        break   
567                endif
568                If (cmpstr("log(I)",ystr) == 0)
569                        if(cmpstr("log(q)",xstr) !=0 )  //extrapolation is nonsense
570                                textstr_3 = "I(q=0) =  "  + num2str(10^(W_coef[0]))
571                        endif
572                        break   
573                endif
574                If (cmpstr("1/I",ystr) == 0)
575                        textstr_3 = "I(q=0) =  "  + num2str(1/W_coef[0])+ " ± " + num2str(1/(W_coef[0] - W_sigma[0])-1/W_coef[0])
576                        break
577                endif
578                If (cmpstr("I^a",ystr) == 0)
579                        //nothing
580                        break
581                endif
582                If (cmpstr("Iq^a",ystr) == 0)
583                        //nothing
584                        break
585                endif
586                If (cmpstr("I^a q^b",ystr) == 0)
587                        //nothing
588                        break
589                endif
590                If (cmpstr("1/sqrt(I)",ystr) == 0)
591                        textstr_3 = "I(q=0) =  "  + num2str((W_coef[0])^2)
592                        break
593                endif
594                If (cmpstr("ln(Iq)",ystr) == 0)
595                        //nothing
596                        if(cmpstr("q^2",xstr) == 0)     //then a x-sect Guinier plot for a rod (2-d)
597                                // rg now is NOT the radius of gyration, but the x-sect DIAMETER
598                                rg = 4*sqrt(-W_coef[1])
599                                rgerr = 8*W_sigma[1]/rg
600                                textstr_3 = "Rod diameter ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
601                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
602                                break
603                        endif
604                        break
605                endif
606                If (cmpstr("ln(Iq^2)",ystr) == 0)
607                        //nothing
608                        if(cmpstr("q^2",xstr) == 0)     //then a 1-d Guinier plot for a sheet
609                                // rg now is NOT the radius of gyration, but the thickness
610                                rg = sqrt(-12*W_coef[1])
611                                rgerr = 6*W_sigma[1]/(2*rg)
612                                textstr_3 = "Platelet thickness ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
613                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
614                                break
615                        endif
616                        break
617                endif
618               
619        while(0)
620        //kill the old textboxes, if they exist
621        TextBox/W=FitWindow/K/N=text_1
622        TextBox/W=FitWindow/K/N=text_2
623        TextBox/W=FitWindow/K/N=text_3
624        // write the new text boxes
625        TextBox/W=FitWindow/N=text_1/A=LT textstr_1
626        TextBox/W=FitWindow/N=text_2/A=LC textstr_2
627        If (cmpstr("",textstr_3) != 0)          //only display textstr_3 if it isn't null
628                TextBox/W=FitWindow/N=text_3/A=LB textstr_3
629        endif
630       
631        //adjust the plot range to reflect the actual fitted range
632        //cursors are already on the graph, done by Rescale_Data()
633        AdjustAxisToCursors()
634       
635End
636
637// adjusts both the x-axis scaling  and y-axis scaling to the cursor range
638// **cursors are already on the graph, done by Rescale_Data()
639//
640// will expand the scale to show an extra 5 points in each direction (if available)
641Function AdjustAxisToCursors()
642
643        DoWindow/F FitWindow
644        WAVE xAxisWave = root:xAxisWave
645        WAVE yAxisWave = root:yAxisWave
646        Variable xlow,xhigh,ylow,yhigh,yptlow,ypthigh
647        Variable extraPts = 5, num=numpnts(xAxisWave)
648       
649        String csrA = CsrInfo(A ,"FitWindow")
650        String csrB = CsrInfo(B ,"FitWindow")
651       
652        //x-levels, these are monotonic
653        Variable ptLow,ptHigh,tmp
654        ptLow = NumberByKey("POINT", csrA ,":" ,";")
655        ptHigh = NumberByKey("POINT", csrB ,":" ,";")
656        if(ptLow > ptHigh)
657                tmp= ptLow
658                ptLow=ptHigh
659                ptHigh=tmp
660        endif
661
662        // keep extended point range in bounds
663        ptLow = (ptLow-extraPts) >= 0 ? ptLow-extraPts : 0
664        ptHigh = (ptHigh+extraPts) <= (num-1) ? ptHigh + extraPts : num-1
665       
666        xlow = xAxisWave[ptLow]
667        xhigh = xAxisWave[ptHigh]
668//old way
669//      xlow = xAxisWave[xcsr(A)]
670//      xhigh = xAxisWave[xcsr(B)]
671//      if(xlow > xhigh)
672//              xhigh = xlow
673//              xlow = xAxisWave[xcsr(B)]
674//      endif
675       
676        //y-levels (old way)
677//      FindLevel/P/Q xAxisWave, xlow
678//      if(V_flag == 1)                 //level NOT found
679//              DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
680//      endif
681//      yptlow = V_LevelX
682//      FindLevel/P/Q xAxisWave, xhigh
683//      if(V_flag == 1)
684//              DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
685//      endif
686//      ypthigh = V_LevelX
687
688//      Print xlow,xhigh,yptlow,ypthigh
689//      Print yAxisWave[yptlow],yAxisWave[ypthigh]
690       
691
692        // make sure ylow/high are in the correct order, since the slope could be + or -
693        yhigh = max(yAxisWave[ptlow],yAxisWave[pthigh])
694        ylow = min(yAxisWave[ptlow],yAxisWave[pthigh])
695       
696//      Print ptLow,ptHigh
697//      print xlow,xhigh
698//      print ylow,yhigh
699       
700        SetAxis bottom,xlow,xhigh
701        SetAxis left ylow,yhigh
702       
703End
704
705
706//****************************************
707//procedures for creating and initializing the FITRPA panel
708//global variables (numerical only) are kept in root:myGlobals:FITRPA folder,
709//created as needed
710//
711// very similar in function to the FIT panel
712//
713Proc OpenFitRPAPanel()
714        If(WinType("FitRPAPanel") == 0)
715                //create the necessary data folder
716                NewDataFolder/O root:myGlobals:FITRPA
717                //initialize the values
718                Variable/G root:myGlobals:FITRPA:gLolim = 0.02
719                Variable/G root:myGlobals:FITRPA:gUplim = 0.04
720                Variable/G root:myGlobals:FITRPA:gBack = 0
721                Variable/G root:myGlobals:FITRPA:gLambda = 6.0
722                PathInfo/S catPathName
723                String localpath = S_path
724                if (V_flag == 0)
725                //path does not exist - no folder selected
726                        String/G root:myGlobals:FITRPA:gPathStr = "no folder selected"
727                else
728                        String/G root:myGlobals:FITRPA:gPathStr = localpath
729                endif
730                String/G    root:myGlobals:FITRPA:gDataPopList = "none"
731                FitRPAPanel()
732        else
733                //window already exists, just bring to front for update
734                DoWindow/F FitRPAPanel
735        endif
736        //pop the menu
737        FilePopMenuProc("",1,"")
738End
739
740//used on the fit/rpa panel to select the path for the data
741// - automatically pops the file list after the new  path selection
742Function FITRPAPickPathButton(ctrlName) : ButtonControl
743        String ctrlName
744
745        Variable err = PickPath()               //sets global path value
746        SVAR pathStr = root:myGlobals:gCatPathStr
747       
748        //set the global string for NSORT to the selected pathname
749        String/G root:myGlobals:FITRPA:gPathStr = pathStr
750       
751        //call the popup menu proc's to re-set the menu choices
752        FilePopMenuProc("filePopup",1,"")
753       
754End
755
756//gets a valid file list (simply not the files with ".SAn" in the name)
757//
758Function FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
759        String ctrlName
760        Variable popNum
761        String popStr
762
763        String tempStr=ReducedDataFileList(ctrlName)
764        if(strlen(tempStr)==0)
765                tempStr = "Pick the data path"
766        Endif
767        String/G root:myGlobals:FITRPA:gDataPopList =   tempStr //function is in NSORT.ipf
768        ControlUpdate filePopup
769       
770End
771
772
773
774// window recreation macro to draw the fit/rpa panel
775//globals and data folders must be present before drawing panel
776//
777Window FitRPAPanel()
778
779        String angst = root:myGlobals:gAngstStr
780        PauseUpdate; Silent 1           // building window...
781        NewPanel /W=(250,266,591,579)/K=1
782        ModifyPanel cbRGB=(32768,54528,65280)
783        ModifyPanel fixedSize=1
784        SetDrawLayer UserBack
785        SetDrawEnv fstyle= 1
786        DrawText 81,19,"Select Experimental Data"
787        SetDrawEnv fstyle= 1
788        DrawText 97,102,"q-range to fit ("+angst+"^-1)"
789        SetDrawEnv fstyle= 1
790        DrawText 87,239,"Select the fit parameters"
791        SetDrawEnv fillpat= 0
792        DrawRect 1,103,338,224
793        SetDrawEnv fillpat= 0
794        DrawRect 1,20,337,83
795        SetDrawEnv fillpat= 0
796        DrawRect 2,241,337,275
797//      Button PathButton,pos={6,26},size={80,20},proc=FitRPAPickPathButton,title="Pick Path"
798//      Button PathButton,help={"Select the local path to the folder containing your SANS data"}
799//      SetVariable setPath,pos={95,29},size={240,17},title="Path:"
800//      SetVariable setPath,help={"The current path to the local folder with SANS data"}
801//      SetVariable setPath,fSize=10
802//      SetVariable setPath,limits={0,0,0},value= root:myGlobals:FITRPA:gPathStr
803        PopupMenu filePopup,pos={8,30},size={96,21},proc=FilePopMenuProc,title="Files"
804        PopupMenu filePopup,help={"Select the data file to load."}
805        PopupMenu filePopup,mode=5,popvalue="none",value= #"root:myGlobals:FITRPA:gDataPopList"
806        SetVariable lambda,pos={111,250},size={120,18},title="Lambda ("+angst+")"
807        SetVariable lambda,help={"This sets the wavelength for the multiple scattering corrections."}
808        SetVariable lambda,limits={0,10,0},value= root:myGlobals:FITRPA:gLambda
809        Button GoFit,pos={60,286},size={80,20},proc=DoFITRPA,title="Do the Fit"
810        Button GoFit,help={"This button will do the specified fit using the selections in this panel"}
811        SetVariable lolim,pos={82,113},size={134,28},title="Lower Limit"
812        SetVariable lolim,help={"Enter the lower q-limit to perform the fit ("+angst+"^-1)"}
813        SetVariable lolim,limits={0,5,0},value= root:myGlobals:FITRPA:gLolim
814        SetVariable uplim,pos={80,140},size={134,28},title="Upper Limit"
815        SetVariable uplim,help={"Enter the upper q-limit to perform the fit ("+angst+"^-1)"}
816        SetVariable uplim,limits={0,5,0},value= root:myGlobals:FITRPA:gUplim
817        CheckBox RPA_check0,pos={64,198},size={190,20},title="Use cursor range from FitWindow"
818        CheckBox RPA_check0,help={"Checking this will perform a fit between the cursors on the graph in FitWindow and ignore the numerical limits typed above"},value=0
819        PopupMenu model,pos={3,249},size={101,21},title="Standard"
820        PopupMenu model,help={"This popup selects which standard should be used to fit this data"}
821        PopupMenu model,mode=1,popvalue="B",value= #"\"B;C;AS\""
822        Button sh_all,pos={82,168},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range"
823        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."}
824        Button loadButton,pos={20,55},size={70,20},proc=FITRPA_Load_Proc,title="Load File"
825        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."}
826        Button helpButton,pos={270,55},size={25,20},proc=showFITHelp,title="?"
827        Button helpButton,help={"Show help file for RPA fitting"}
828        Button DoneButton,pos={200,286},size={50,20},proc=FITRPADoneButton,title="Done"
829        Button DoneButton,help={"This button will close the panel and the associated graph"}
830EndMacro
831
832
833Proc FITRPADoneButton(ctrlName) : ButtonControl
834        String ctrlName
835        DoWindow/K FitWindow
836        DoWindow/K FitRPAPanel
837end
838
839//dispatches the fit to the appropriate model
840//and the appropriate range, based on selections in the panel
841//
842Proc DoFITRPA(ctrlName) : ButtonControl
843        String ctrlName
844        //
845        String cleanLastFileName = CleanupName(root:Packages:NIST:gLastFileName,0)
846        String tmpStr = "root:"+cleanLastFileName+":"+cleanLastFileName
847
848        Duplicate/O $(tmpStr+"_q") xAxisWave
849        Duplicate/O $(tmpStr+"_i") yAxisWave
850        Duplicate/O $(tmpStr+"_s") yErrWave,yWtWave,residWave
851
852        yWtWave = 1/yErrWave
853       
854        String xlabel = "q (A^-1)"
855        String ylabel = "Intensity"
856        //Check to see if the FitWindow exists
857        //Plot the data in a FitWindow
858        If(WinType("FitWindow") == 0)
859                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
860                ModifyGraph mode=3,standoff=0,marker=8
861                ErrorBars/T=0 yAxisWave Y,wave=(yErrWave,yErrWave)
862                DoWindow/C FitWindow
863                ShowInfo
864        else
865                //window already exists, just bring to front for update
866                DoWindow/F FitWindow
867                // remove old text boxes
868                TextBox/K/N=text_1
869                TextBox/K/N=text_2
870                TextBox/K/N=text_3
871        endif
872       
873        //see if the user wants to use the data specified by the cursors - else use numerical values
874        Variable xlow,xhigh,ylow,yhigh,yes_cursors
875        ControlInfo/W=FitRPAPanel RPA_check0            //V_value = 1 if it is checked, meaning yes, use cursors
876        yes_cursors = V_value
877
878        ControlInfo/W=FitRPAPanel lolim
879        xlow = V_value
880        ControlInfo/W=FitRPAPanel uplim
881        xhigh = V_value
882        if(yes_cursors)
883                xlow = xAxisWave[xcsr(A)]
884                xhigh = xAxisWave[xcsr(B)]
885                if(xlow > xhigh)
886                        xhigh = xlow
887                        xlow = xAxisWave[xcsr(B)]
888                endif
889//              Print "xlow,xhigh = ",xlow,xhigh
890                ylow = yAxisWave[xcsr(A)]
891                yhigh = yAxisWave[xcsr(B)]
892                if(ylow > yhigh)
893                        ylow=yhigh
894                        yhigh = yAxisWave[xcsr(A)]
895                endif
896        else
897                FindLevel/P/Q xAxisWave, xlow
898                if(V_flag == 1)                 //level NOT found
899                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
900                        DoWindow/K FitWindow
901                        Abort
902                endif
903                Cursor/P A, yAxisWave,trunc(V_LevelX)+1
904                ylow = yAxisWave[V_LevelX]
905                FindLevel/P/Q xAxisWave, xhigh
906                if(V_flag == 1)
907                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
908                        DoWindow/K FitWindow
909                        Abort
910                endif
911                Cursor/P B, yAxisWave,trunc(V_LevelX)
912                yhigh = yAxisWave[V_LevelX]
913                if(ylow > yhigh)
914                        yhigh=ylow
915                        ylow = yAxisWave[V_levelX]
916                endif
917        endif   //if(V_value)
918        SetAxis bottom,xlow,xhigh
919       
920//      print "ylow,yhigh",ylow,yhigh
921       
922        //Get the rest of the data from the panel
923        //such as which standard, the wavelength
924        ControlInfo/W=FitRPAPanel model
925
926        //find the model name
927        String modelName = S_value
928       
929        Variable first_guess, seglength,iabs,iarb,thick
930        Make/D/O/N=2 fitParams
931       
932        seglength = 6.8
933       
934        first_guess = 1.0
935        fitParams[0] = first_guess
936        fitParams[1] = seglength
937
938        If (cmpstr(modelName,"B")==0)
939                iabs = BStandardFunction(fitParams,xlow)
940                fitParams[0] = yhigh/iabs
941//              Print fitParams[0],fitParams[1]
942                FuncFit BStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
943                iarb = BStandardFunction(fitParams, 0.0)
944                iabs = iarb/fitParams[0]
945                thick = 0.153
946        endif
947        If (cmpstr(modelName,"C")==0)
948                iabs = CStandardFunction(fitParams,xlow)
949                fitParams[0] = yhigh/iabs
950                FuncFit CStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
951                iarb = CStandardFunction(fitParams, 0.0)
952                iabs = iarb/fitParams[0]
953                thick= 0.153
954        endif
955        If (cmpstr(modelName,"AS")==0)
956                iabs = ASStandardFunction(fitParams,xlow)
957                fitParams[0] = yhigh/iabs
958                FuncFit ASStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
959                iarb = ASStandardFunction(fitParams, 0.0)
960                iabs = iarb/fitParams[0]
961                thick = 0.1
962        endif
963        ModifyGraph rgb(fit_yAxisWave)=(0,0,0)
964        Label left ylabel
965        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead       
966
967        ControlInfo/W=FitRPAPanel lambda
968       
969        Variable cor_mult = 1.0 + 2.2e-4*V_Value^2
970       
971        //WAVE W_coef=W_coef
972        //WAVE W_sigma=W_sigma
973        String textstr_1,textstr_2,textstr_3 = ""
974        textstr_1 = "Scaling Parameter: "+num2str(fitParams[0])+" ± "+num2str(W_sigma[0])
975        textstr_1 += "\rSegment Length: "+num2str(fitParams[1])+" ± "+num2str(W_sigma[1])
976        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3))
977       
978        textstr_2 = "Cross section at q=0:  Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
979        textstr_2 += "\rData extrapolated to q=0: Im(0) = "+num2str(iarb)+" Counts/(10\S8\M  Mon cts)"
980        textstr_2 += "\rData corrected for multiple scattering: I(0) = "+num2str(iarb/cor_mult)+" Counts/(10\S8\M  Mon cnts)"
981       
982        textstr_3 = "In the ABS protocol, "
983        textstr_3 += "\rStandard Thickness, d = "+num2str(thick)+"cm"
984        textstr_3 += "\rI(0), Iarb(0) = "+num2str(iarb/cor_mult)+"Counts/(10\S8\M Mon cts)"
985        textstr_3 += "\rStandard Cross Section, Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
986        TextBox/K/N=text_1
987        TextBox/K/N=text_2
988        TextBox/K/N=text_3
989        TextBox/N=text_2/A=RT textstr_2
990        TextBox/N=text_3/A=RC textstr_3
991        TextBox/N=text_1/A=RB textstr_1
992       
993End
994
995//loads the file selected in the popup for fitting with POL
996//standard functions. Reads the wavelength from the header, using
997//6 A as the default
998//plots the data in FitWindow after reading the file
999//updates lambda and full q-range  on the Panel
1000//
1001Proc FITRPA_Load_Proc(ctrlName): ButtonControl
1002        String ctrlName
1003        //Load the data
1004        String tempName="",partialName=""
1005        Variable err
1006        ControlInfo $"filePopup"
1007        //find the file from the partial filename
1008        If( (cmpstr(S_value,"")==0) || (cmpstr(S_value,"none")==0) )
1009                //null selection, or "none" from any popup
1010                Abort "no file selected in popup menu"
1011        else
1012                //selection not null
1013                partialName = S_value
1014                //Print partialName
1015        Endif
1016        //get a valid file based on this partialName and catPathName
1017        tempName = FindValidFilename(partialName)
1018
1019        Variable lambdaFromFile=GetLambdaFromReducedData(tempName)
1020        Variable/G root:myGlobals:FITRPA:gLambda = lambdaFromFile
1021        Print "Lambda in file read as:", lambdaFromFile
1022       
1023        //prepend path to tempName for read routine
1024        PathInfo catPathName
1025        tempName = S_path + tempName
1026       
1027        //load in the data (into the root directory)
1028        A_LoadOneDDataWithName(tempName,0)
1029        //Print S_fileName
1030        //Print tempName
1031       
1032        String cleanLastFileName = CleanupName(root:Packages:NIST:gLastFileName,0)
1033        String tmpStr = "root:"+cleanLastFileName+":"+cleanLastFileName
1034
1035        Duplicate/o $(tmpStr+"_q") xAxisWave
1036        Duplicate/o $(tmpStr+"_i") yAxisWave
1037        Duplicate/o $(tmpStr+"_s") yErrWave
1038       
1039        Variable xmin, xmax
1040        WaveStats/Q xAxisWave
1041        root:myGlobals:FITRPA:gLolim=V_min
1042        root:myGlobals:FITRPA:gUplim=V_max
1043        ControlUpdate/W=FITRPAPanel/A
1044       
1045        //Check to see if the FitWindow exists
1046        //Plot the data in a FitWindow
1047        If(WinType("FitWindow") == 0)
1048                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
1049                ModifyGraph mode=3,standoff=0,marker=8
1050                ErrorBars/T=0 yAxisWave Y,wave=(yErrWave,yErrWave)
1051                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1052                DoWindow/C FitWindow
1053                ShowInfo
1054        else
1055                //window already exists, just bring to front for update
1056                DoWindow/F FitWindow
1057                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1058        endif
1059        // remove old text boxes
1060        TextBox/K/N=text_1
1061        TextBox/K/N=text_2
1062        TextBox/K/N=text_3
1063        RemoveFromGraph/W=fitWindow /Z fit_yAxisWave
1064        SetAxis/A
1065       
1066        //put cursors on the graph at first and last points
1067        Cursor/P A  yAxisWave  0
1068        Cursor/P B yAxisWave (numpnts(yAxisWave) - 1)
1069End
1070
1071//Fitting function for the POL-B standards
1072//
1073Function BStandardFunction(parameterWave, x)
1074        Wave parameterWave; Variable x
1075       
1076        //Model parameters
1077        Variable KN=4.114E-3,CH=9.613E4, CD=7.558E4, NH=1872, ND=1556
1078        Variable INC=0.32, CHIV=2.2E-6
1079        //Correction based on absolute flux measured 5/93
1080        Variable CORR = 1.1445
1081       
1082        //Local variables
1083        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1084       
1085        //Calculate the function here
1086        ap2=parameterWave[1]^2
1087        qrgh = x*sqrt(nh*ap2/6)
1088        qrgd = x*sqrt(nd*ap2/6)
1089        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1090        iabs = corr*iabs_rpa + inc
1091       
1092        //return the result
1093        return parameterWave[0]*iabs
1094       
1095End
1096
1097//Fitting function for the POL-C standards
1098//
1099Function CStandardFunction(parameterWave, x)
1100        Wave parameterWave; Variable x
1101       
1102        //Model parameters
1103        Variable KN=4.114E-3,CH=2.564E5, CD=1.912E5, NH=4993, ND=3937
1104        Variable INC=0.32, CHIV=2.2E-6
1105        //Correction based on absolute flux measured 5/93
1106        Variable CORR = 1.0944
1107       
1108        //Local variables
1109        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1110       
1111        //Calculate the function here
1112        ap2=parameterWave[1]^2
1113        qrgh = x*sqrt(nh*ap2/6)
1114        qrgd = x*sqrt(nd*ap2/6)
1115        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1116        iabs = corr*iabs_rpa + inc
1117       
1118        //return the result
1119        return parameterWave[0]*iabs
1120       
1121End
1122
1123//fitting function for the POL-AS standards
1124//
1125Function ASStandardFunction(parameterWave, x)
1126        Wave parameterWave; Variable x
1127       
1128        //Model parameters
1129        Variable KN=64.5,CH=1.0, CD=1.0, NH=766, ND=766
1130        Variable INC=0.32, CHIV=0.0
1131       
1132        //Local variables
1133        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1134       
1135        //Calculate the function here
1136        ap2=parameterWave[1]^2
1137        qrgh = x*sqrt(nh*ap2/6)
1138
1139//The following lines were commented out in the fortran function
1140        //qrgd = x*sqrt(nd*ap2/6)
1141        //iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1142
1143        iabs_rpa = kn*FIT_dbf(qrgh)
1144        iabs = iabs_rpa + inc
1145       
1146        //return the result
1147        return parameterWave[0]*iabs
1148       
1149End
1150
1151//Debye Function used for polymer standards
1152Function FIT_dbf(rgq)
1153        Variable rgq
1154        Variable x
1155       
1156        x=rgq*rgq
1157        if (x < 5.0E-3)
1158                return 1.0 - x/3 + x^2/12
1159        else
1160                return 2*(exp(-x) + x - 1)/x^2
1161        endif
1162End
Note: See TracBrowser for help on using the repository browser.