source: sans/SANSReduction/branches/kline_29MAR07/Put in User Procedures/SANS_Reduction_v5.00/FIT_Ops.ipf @ 82

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

2nd pass of pulling out NCNR calls. Also cleared a lot of deadwood from the code, removing depricated functions that were mostly already commented out.

File size: 36.2 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=ReducedDataFileList(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=ReducedDataFileList(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
735
736// window recreation macro to draw the fit/rpa panel
737//globals and data folders must be present before drawing panel
738//
739Window FitRPAPanel()
740
741        String angst = root:myGlobals:gAngstStr
742        PauseUpdate; Silent 1           // building window...
743        NewPanel /W=(250,266,591,579)/K=1
744        ModifyPanel cbRGB=(32768,54528,65280)
745        ModifyPanel fixedSize=1
746        SetDrawLayer UserBack
747        SetDrawEnv fstyle= 1
748        DrawText 81,19,"Select Experimental Data"
749        SetDrawEnv fstyle= 1
750        DrawText 97,102,"q-range to fit ("+angst+"^-1)"
751        SetDrawEnv fstyle= 1
752        DrawText 87,239,"Select the fit parameters"
753        SetDrawEnv fillpat= 0
754        DrawRect 1,103,338,224
755        SetDrawEnv fillpat= 0
756        DrawRect 1,20,337,83
757        SetDrawEnv fillpat= 0
758        DrawRect 2,241,337,275
759//      Button PathButton,pos={6,26},size={80,20},proc=FitRPAPickPathButton,title="Pick Path"
760//      Button PathButton,help={"Select the local path to the folder containing your SANS data"}
761//      SetVariable setPath,pos={95,29},size={240,17},title="Path:"
762//      SetVariable setPath,help={"The current path to the local folder with SANS data"}
763//      SetVariable setPath,fSize=10
764//      SetVariable setPath,limits={0,0,0},value= root:myGlobals:FITRPA:gPathStr
765        PopupMenu filePopup,pos={8,30},size={96,21},proc=FilePopMenuProc,title="Files"
766        PopupMenu filePopup,help={"Select the data file to load."}
767        PopupMenu filePopup,mode=5,popvalue="none",value= #"root:myGlobals:FITRPA:gDataPopList"
768        SetVariable lambda,pos={111,250},size={120,18},title="Lambda ("+angst+")"
769        SetVariable lambda,help={"This sets the wavelength for the multiple scattering corrections."}
770        SetVariable lambda,limits={0,10,0},value= root:myGlobals:FITRPA:gLambda
771        Button GoFit,pos={60,286},size={80,20},proc=DoFITRPA,title="Do the Fit"
772        Button GoFit,help={"This button will do the specified fit using the selections in this panel"}
773        SetVariable lolim,pos={82,113},size={134,28},title="Lower Limit"
774        SetVariable lolim,help={"Enter the lower q-limit to perform the fit ("+angst+"^-1)"}
775        SetVariable lolim,limits={0,5,0},value= root:myGlobals:FITRPA:gLolim
776        SetVariable uplim,pos={80,140},size={134,28},title="Upper Limit"
777        SetVariable uplim,help={"Enter the upper q-limit to perform the fit ("+angst+"^-1)"}
778        SetVariable uplim,limits={0,5,0},value= root:myGlobals:FITRPA:gUplim
779        CheckBox RPA_check0,pos={64,198},size={190,20},title="Use cursor range from FitWindow"
780        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
781        PopupMenu model,pos={3,249},size={101,21},title="Standard"
782        PopupMenu model,help={"This popup selects which standard should be used to fit this data"}
783        PopupMenu model,mode=1,popvalue="B",value= #"\"B;C;AS\""
784        Button sh_all,pos={82,168},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range"
785        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."}
786        Button loadButton,pos={20,55},size={70,20},proc=FITRPA_Load_Proc,title="Load File"
787        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."}
788        Button helpButton,pos={270,55},size={25,20},proc=showFITHelp,title="?"
789        Button helpButton,help={"Show help file for RPA fitting"}
790        Button DoneButton,pos={200,286},size={50,20},proc=FITRPADoneButton,title="Done"
791        Button DoneButton,help={"This button will close the panel and the associated graph"}
792EndMacro
793
794
795Proc FITRPADoneButton(ctrlName) : ButtonControl
796        String ctrlName
797        DoWindow/K FitWindow
798        DoWindow/K FitRPAPanel
799end
800
801//dispatches the fit to the appropriate model
802//and the appropriate range, based on selections in the panel
803//
804Proc DoFITRPA(ctrlName) : ButtonControl
805        String ctrlName
806        //
807        String cleanLastFileName = "root:"+CleanupName(root:myGlobals:gLastFileName,0)
808        String tempName
809
810        tempName=cleanLastFileName+"_q"
811        Duplicate/O $tempName xAxisWave
812        tempName=cleanLastFileName+"_i"
813        Duplicate/O $tempName yAxisWave
814        tempName=cleanLastFileName+"_s"
815        Duplicate/O $tempName yErrWave
816        Duplicate/O $tempName yWtWave
817        Duplicate/O $tempName residWave
818        yWtWave = 1/yErrWave
819       
820        String xlabel = "q (A^-1)"
821        String ylabel = "Intensity"
822        //Check to see if the FitWindow exists
823        //Plot the data in a FitWindow
824        If(WinType("FitWindow") == 0)
825                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
826                ModifyGraph mode=3,standoff=0,marker=8
827                ErrorBars yAxisWave Y,wave=(yErrWave,yErrWave)
828                DoWindow/C FitWindow
829                ShowInfo
830        else
831                //window already exists, just bring to front for update
832                DoWindow/F FitWindow
833                // remove old text boxes
834                TextBox/K/N=text_1
835                TextBox/K/N=text_2
836                TextBox/K/N=text_3
837        endif
838       
839        //see if the user wants to use the data specified by the cursors - else use numerical values
840        Variable xlow,xhigh,ylow,yhigh,yes_cursors
841        ControlInfo/W=FitRPAPanel RPA_check0            //V_value = 1 if it is checked, meaning yes, use cursors
842        yes_cursors = V_value
843
844        ControlInfo/W=FitRPAPanel lolim
845        xlow = V_value
846        ControlInfo/W=FitRPAPanel uplim
847        xhigh = V_value
848        if(yes_cursors)
849                xlow = xAxisWave[xcsr(A)]
850                xhigh = xAxisWave[xcsr(B)]
851                if(xlow > xhigh)
852                        xhigh = xlow
853                        xlow = xAxisWave[xcsr(B)]
854                endif
855//              Print "xlow,xhigh = ",xlow,xhigh
856                ylow = yAxisWave[xcsr(A)]
857                yhigh = yAxisWave[xcsr(B)]
858                if(ylow > yhigh)
859                        ylow=yhigh
860                        yhigh = yAxisWave[xcsr(A)]
861                endif
862        else
863                FindLevel/P/Q xAxisWave, xlow
864                if(V_flag == 1)                 //level NOT found
865                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
866                        DoWindow/K FitWindow
867                        Abort
868                endif
869                Cursor/P A, yAxisWave,trunc(V_LevelX)+1
870                ylow = yAxisWave[V_LevelX]
871                FindLevel/P/Q xAxisWave, xhigh
872                if(V_flag == 1)
873                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
874                        DoWindow/K FitWindow
875                        Abort
876                endif
877                Cursor/P B, yAxisWave,trunc(V_LevelX)
878                yhigh = yAxisWave[V_LevelX]
879                if(ylow > yhigh)
880                        yhigh=ylow
881                        ylow = yAxisWave[V_levelX]
882                endif
883        endif   //if(V_value)
884        SetAxis bottom,xlow,xhigh
885       
886//      print "ylow,yhigh",ylow,yhigh
887       
888        //Get the rest of the data from the panel
889        //such as which standard, the wavelength
890        ControlInfo/W=FitRPAPanel model
891
892        //find the model name
893        String modelName = S_value
894       
895        Variable first_guess, seglength,iabs,iarb,thick
896        Make/D/O/N=2 fitParams
897       
898        seglength = 6.8
899       
900        first_guess = 1.0
901        fitParams[0] = first_guess
902        fitParams[1] = seglength
903
904        If (cmpstr(modelName,"B")==0)
905                iabs = BStandardFunction(fitParams,xlow)
906                fitParams[0] = yhigh/iabs
907//              Print fitParams[0],fitParams[1]
908                FuncFit BStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
909                iarb = BStandardFunction(fitParams, 0.0)
910                iabs = iarb/fitParams[0]
911                thick = 0.153
912        endif
913        If (cmpstr(modelName,"C")==0)
914                iabs = CStandardFunction(fitParams,xlow)
915                fitParams[0] = yhigh/iabs
916                FuncFit CStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
917                iarb = CStandardFunction(fitParams, 0.0)
918                iabs = iarb/fitParams[0]
919                thick= 0.153
920        endif
921        If (cmpstr(modelName,"AS")==0)
922                iabs = ASStandardFunction(fitParams,xlow)
923                fitParams[0] = yhigh/iabs
924                FuncFit ASStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
925                iarb = ASStandardFunction(fitParams, 0.0)
926                iabs = iarb/fitParams[0]
927                thick = 0.1
928        endif
929        ModifyGraph rgb(fit_yAxisWave)=(0,0,0)
930        Label left ylabel
931        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead       
932
933        ControlInfo/W=FitRPAPanel lambda
934       
935        Variable cor_mult = 1.0 + 2.2e-4*V_Value^2
936       
937        //WAVE W_coef=W_coef
938        //WAVE W_sigma=W_sigma
939        String textstr_1,textstr_2,textstr_3 = ""
940        textstr_1 = "Scaling Parameter: "+num2str(fitParams[0])+" ± "+num2str(W_sigma[0])
941        textstr_1 += "\rSegment Length: "+num2str(fitParams[1])+" ± "+num2str(W_sigma[1])
942        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3))
943       
944        textstr_2 = "Cross section at q=0:  Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
945        textstr_2 += "\rData extrapolated to q=0: Im(0) = "+num2str(iarb)+" Counts/(10\S8\M  Mon cts)"
946        textstr_2 += "\rData corrected for multiple scattering: I(0) = "+num2str(iarb/cor_mult)+" Counts/(10\S8\M  Mon cnts)"
947       
948        textstr_3 = "In the ABS protocol, "
949        textstr_3 += "\rStandard Thickness, d = "+num2str(thick)+"cm"
950        textstr_3 += "\rI(0), Iarb(0) = "+num2str(iarb/cor_mult)+"Counts/(10\S8\M Mon cts)"
951        textstr_3 += "\rStandard Cross Section, Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
952        TextBox/K/N=text_1
953        TextBox/K/N=text_2
954        TextBox/K/N=text_3
955        TextBox/N=text_2/A=RT textstr_2
956        TextBox/N=text_3/A=RC textstr_3
957        TextBox/N=text_1/A=RB textstr_1
958       
959End
960
961//loads the file selected in the popup for fitting with POL
962//standard functions. Reads the wavelength from the header, using
963//6 A as the default
964//plots the data in FitWindow after reading the file
965//updates lambda and full q-range  on the Panel
966//
967Proc FITRPA_Load_Proc(ctrlName): ButtonControl
968        String ctrlName
969        //Load the data
970        String tempName="",partialName=""
971        Variable err
972        ControlInfo $"filePopup"
973        //find the file from the partial filename
974        If( (cmpstr(S_value,"")==0) || (cmpstr(S_value,"none")==0) )
975                //null selection, or "none" from any popup
976                Abort "no file selected in popup menu"
977        else
978                //selection not null
979                partialName = S_value
980                //Print partialName
981        Endif
982        //get a valid file based on this partialName and catPathName
983        tempName = FindValidFilename(partialName)
984
985        Variable lambdaFromFile=GetLambdaFromReducedData(tempName)
986        Variable/G root:myGlobals:FITRPA:gLambda = lambdaFromFile
987        Print "Lambda in file read as:", lambdaFromFile
988       
989        //prepend path to tempName for read routine
990        PathInfo catPathName
991        tempName = S_path + tempName
992       
993        //load in the data (into the root directory)
994        LoadOneDDataWithName(tempName)
995        //Print S_fileName
996        //Print tempName
997       
998        String cleanLastFileName = "root:"+CleanupName(root:myGlobals:gLastFileName,0)
999
1000        tempName=cleanLastFileName+"_q"
1001        Duplicate/o $tempName xAxisWave
1002        tempName=cleanLastFileName+"_i"
1003        Duplicate/o $tempName yAxisWave
1004        tempName=cleanLastFileName+"_s"
1005        Duplicate/o $tempName yErrWave
1006        Variable xmin, xmax
1007        WaveStats/Q xAxisWave
1008        root:myGlobals:FITRPA:gLolim=V_min
1009        root:myGlobals:FITRPA:gUplim=V_max
1010        ControlUpdate/W=FITRPAPanel/A
1011       
1012        //Check to see if the FitWindow exists
1013        //Plot the data in a FitWindow
1014        If(WinType("FitWindow") == 0)
1015                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
1016                ModifyGraph mode=3,standoff=0,marker=8
1017                ErrorBars yAxisWave Y,wave=(yErrWave,yErrWave)
1018                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1019                DoWindow/C FitWindow
1020                ShowInfo
1021        else
1022                //window already exists, just bring to front for update
1023                DoWindow/F FitWindow
1024                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1025        endif
1026        // remove old text boxes
1027        TextBox/K/N=text_1
1028        TextBox/K/N=text_2
1029        TextBox/K/N=text_3
1030        RemoveFromGraph/W=fitWindow /Z fit_yAxisWave
1031        SetAxis/A
1032       
1033        //put cursors on the graph at first and last points
1034        Cursor/P A  yAxisWave  0
1035        Cursor/P B yAxisWave (numpnts(yAxisWave) - 1)
1036End
1037
1038//Fitting function for the POL-B standards
1039//
1040Function BStandardFunction(parameterWave, x)
1041        Wave parameterWave; Variable x
1042       
1043        //Model parameters
1044        Variable KN=4.114E-3,CH=9.613E4, CD=7.558E4, NH=1872, ND=1556
1045        Variable INC=0.32, CHIV=2.2E-6
1046        //Correction based on absolute flux measured 5/93
1047        Variable CORR = 1.1445
1048       
1049        //Local variables
1050        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1051       
1052        //Calculate the function here
1053        ap2=parameterWave[1]^2
1054        qrgh = x*sqrt(nh*ap2/6)
1055        qrgd = x*sqrt(nd*ap2/6)
1056        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1057        iabs = corr*iabs_rpa + inc
1058       
1059        //return the result
1060        return parameterWave[0]*iabs
1061       
1062End
1063
1064//Fitting function for the POL-C standards
1065//
1066Function CStandardFunction(parameterWave, x)
1067        Wave parameterWave; Variable x
1068       
1069        //Model parameters
1070        Variable KN=4.114E-3,CH=2.564E5, CD=1.912E5, NH=4993, ND=3937
1071        Variable INC=0.32, CHIV=2.2E-6
1072        //Correction based on absolute flux measured 5/93
1073        Variable CORR = 1.0944
1074       
1075        //Local variables
1076        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1077       
1078        //Calculate the function here
1079        ap2=parameterWave[1]^2
1080        qrgh = x*sqrt(nh*ap2/6)
1081        qrgd = x*sqrt(nd*ap2/6)
1082        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1083        iabs = corr*iabs_rpa + inc
1084       
1085        //return the result
1086        return parameterWave[0]*iabs
1087       
1088End
1089
1090//fitting function for the POL-AS standards
1091//
1092Function ASStandardFunction(parameterWave, x)
1093        Wave parameterWave; Variable x
1094       
1095        //Model parameters
1096        Variable KN=64.5,CH=1.0, CD=1.0, NH=766, ND=766
1097        Variable INC=0.32, CHIV=0.0
1098       
1099        //Local variables
1100        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1101       
1102        //Calculate the function here
1103        ap2=parameterWave[1]^2
1104        qrgh = x*sqrt(nh*ap2/6)
1105
1106//The following lines were commented out in the fortran function
1107        //qrgd = x*sqrt(nd*ap2/6)
1108        //iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1109
1110        iabs_rpa = kn*FIT_dbf(qrgh)
1111        iabs = iabs_rpa + inc
1112       
1113        //return the result
1114        return parameterWave[0]*iabs
1115       
1116End
1117
1118//Debye Function used for polymer standards
1119Function FIT_dbf(rgq)
1120        Variable rgq
1121        Variable x
1122       
1123        x=rgq*rgq
1124        if (x < 5.0E-3)
1125                return 1.0 - x/3 + x^2/12
1126        else
1127                return 2*(exp(-x) + x - 1)/x^2
1128        endif
1129End
1130
Note: See TracBrowser for help on using the repository browser.