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

Last change on this file since 564 was 564, checked in by srkline, 13 years ago

A variety of changes...

  • fixed model strings for Sum_Model
  • /Z flag for RTHelp, Transmission help, FIT help, MRED help
  • added help and done buttons to SASCALC
  • added comments in MonteCarlo?
  • restored the 2-line buttons in USANS panels, now supported in Igor 6.1 (COR_Graph and Main_USANS)
  • added simple text notebook of UCALC configuration, and a save of PNG of the whole panel
  • added help and done buttons to UCALC
  • added links in SANS_Reduction_Help to the simluation help file (not added to SVN yet)
File size: 37.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.0
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.