source: sans/SANSReduction/trunk/Put in User Procedures/SANS_Reduction_v5.00/FIT_Ops.ipf @ 47

Last change on this file since 47 was 47, checked in by srkline, 16 years ago

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