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

Last change on this file since 412 was 412, checked in by ajj, 14 years ago

More reorg.

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