Mosspec: Process Mössbauer Spectrum Files

Package Initialization

This section contains the initialization code and the definitions of the functions and variables in this Mosspec package.

Initialization: Check the Mathematica version and start the application

If[!TrueQ[$VersionNumber≥ 6.0],
(* v5 or lower: a wrong version dialog, then close this notebook *)
NotebookPut[Notebook[{Cell["You need Mathematica v. 6 or higher to run this application.","Subsubtitle" ],ButtonBox["\tOK\t",ButtonFrame->"DialogBox",Active->True,ButtonFunction:>(NotebookClose[ButtonNotebook[]]&),ButtonMargins-> 5]
},
ShowCellBracket->False,Saveable->False,
WindowTitle->"Mosspec",
WindowElements-> {},
WindowFloating-> True,
WindowFrame-> "ModalDialog",
WindowSize->{250,100}
]
];NotebookClose[EvaluationNotebook[]],

(* v6 or higher: start CurveFit (if necessary) and start Mosspec *)
Needs["CurveFit`"];
NotebookLocate["package"];
SelectionMove[EvaluationNotebook[],All,CellGroup];
SelectionEvaluateCreateCell[EvaluationNotebook[]];
SelectionMove[EvaluationNotebook[],After,Notebook]
];

Package definition section

The public symbols defined by the package

If[!TrueQ[$VersionNumber≥ 7.0],
BeginPackage["Mosspec`",{"NonlinearRegression`","CurveFit`"}];
$ContextPath=Join[{"NonlinearRegression`","RegressionCommon`"},$ContextPath],
BeginPackage["Mosspec`",{"CurveFit`"}]
];

(* clear everything first, so the definitions may be read multiple times *)
Unprotect["`*"]; 
ClearAll["`*"];

Package usage statements

mosspecPalette::usage="Starts the Mosspec Functions palette.";

LoadMosspecFile::usage="LoadMosspecFile[filename] loads the named file containing Experiment 29 Mössbauer spectrum data saved using the MOSSPEC program.\n\n"<>
"LoadMosspecFile[ ] (no file name given) uses the file name of the last file loaded or saved, as defined by the value of CurveFit's DataFileName.\n\n"<>
"LoadMosspecFile needs CurveFit in order to function.";

FeLineV::usage=
"FeLineV[] returns a list of the standard velocities of the six hyperfine lines in a natural iron Mössbauer spectrum (in mm/sec). \n\n"<>
"FeLineV[(Subscript[v, iso]), (Subscript[v, B]), (Subscript[v, Q]), (Subscript[g, ratio])] returns a list of velocities of the six hyperfine lines in an iron Mössbauer spectrum (in mm/sec) calculated from the parameters given: \n\n"<>
"(Subscript[v, iso]) :\t Isomeric velocity shift (in mm/sec) \n"<>"(Subscript[v, B]) :\t Hyperfine velocity shift sensitivity (in mm/sec) \n"<>"(Subscript[v, Q]) :\t Quadrupole velocity shift sensitivity (in mm/sec) \n"<>"(Subscript[g, ratio]) :\t (SuperStar[g])/g (ratio of gyromagnetic ratios)\n\n"<>" The order of the lines in the list is the same as the order for FeLineIz[].";

FeLineFunction::usage="FeLineFunction[(Subscript[v, iso]), (Subscript[v, B]), (Subscript[v, Q]), (Subscript[g, ratio])][(Subscript[I, z]), (SuperStar[(Subscript[I, z])])] =\n"<>
"\!\(\*TagBox[FormBox[RowBox[{SubscriptBox["v", "iso"], "+", RowBox[{SubscriptBox["v", "B"], " ", RowBox[{"(", RowBox[{SubscriptBox["I", "z"], "-", RowBox[{SubscriptBox["g", "ratio"], " ", SuperscriptBox[SubscriptBox["I", "z"], "*"]}]}], ")"}]}], "+", RowBox[{SubscriptBox["v", "Q"], " ", RowBox[{"(", RowBox[{SuperscriptBox[RowBox[{"(", SuperscriptBox[SubscriptBox["I", "z"], "*"], ")"}], "2"], "-", FractionBox["5", "4"]}], ")"}]}]}],TraditionalForm],TraditionalForm,Editable->True]\)";

FeLineIz::usage=
"FeLineIz[] returns a list of pairs {(Subscript[I, z]), (SuperStar[(Subscript[I, z])])} for the six hyperfine lines in an iron Mössbauer spectrum. The order of the lines in the list is from lowest velocity to highest.";

FeLineQ::usage=
"FeLineQ[] returns a list of quadrupole component values (±1) for the six hyperfine lines in an iron Mössbauer spectrum. The order of the lines in the list is the same as the order for FeLineIz[].";

LineSpectrum::usage="LineSpectrum[pars][x] returns the amplitude at position x of an emission line spectrum determined by the argument list pars. The spectrum has a baseline of 0 and consists of one or more positive-going Lorentzian peaks.\n\n"<>
"The structure of the argument list determines the number of peaks in the spectrum as follows:\n\n"<>
"LineSpectrum[{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]),(Subscript[x, 5]),(Subscript[x, 6])},(Subscript[y, 1]),(Subscript[y, 2]),γ]:  6 lines at positions (Subscript[x, 1]) .. (Subscript[x, 6]), all with width γ. Lines (Subscript[x, 3]) and (Subscript[x, 4]) have height 1, (Subscript[x, 2]) and (Subscript[x, 5]) have height (Subscript[y, 2]), and (Subscript[x, 1]) and (Subscript[x, 6]) have height (Subscript[y, 1]).\n\n"<>
"LineSpectrum[{(Subscript[x, 1]),(Subscript[x, 2])},(Subscript[y, 2]),γ]:  2 lines at positions (Subscript[x, 1]) and (Subscript[x, 2]), both with width γ. Line (Subscript[x, 1]) has height 1, (Subscript[x, 2]) has height (Subscript[y, 2]).\n\n"<>
"LineSpectrum[{(Subscript[x, 1])},γ] or LineSpectrum[(Subscript[x, 1]),γ]:  a single line at position (Subscript[x, 1]) with width γ and height 1. See MossSpectrum.";

MossSpectrum::usage="MossSpectrum[LineSpectrum[pars],depth][x] returns the amplitude at position x of a simulated Mössbauer absorption spectrum corresponding to the spectrum LineSpectrum[pars]. It has a baseline of 1 and has lines reaching to depth below the baseline.";

Calibrate::usage="Calibrate[] performs a velocity calibration using the current CurveFit data set. The data must be of a Natural Iron hyperfine Mössbauer spectrum, and it must be windowed to include only the lower half or upper half of the double Mössbauer spectrum. The two halves of the spectrum may be calibrated independently; each calibration will be saved.\n\n"<>
"The Natural Iron spectrum is fit using the line velocities returned by the call: FeLineV[]. The fit determines the channel-to-velocity conversion used by the function ChanToVelocity[].";

CalLowerQ::usage="CalLowerQ[] returns True if a velocity calibration has been performed on the lower half of the MOSSPEC spectrum. See: Calibrate.";

CalUpperQ::usage="CalUpperQ[] returns True if a velocity calibration has been performed on the upper half of the MOSSPEC spectrum. See: Calibrate.";

ChanToVelocity::usage = "ChanToVelocity[] converts the CurveFit X channel data to velocity (in mm/sec) using the results of Calibrate[], as saved in the list returned by a call to VelocityParameters[]. It uses either the lower half or upper half velocity calibration as appropriate. No adjustments to X sigmas are made, since the velocity calibration transformation introduces a systematic uncertainty.\n\n"<>
"ChanToVelocity[chan] converts the channel number chan to velocity. It uses the lower calibration if chan < 511, otherwise it uses the upper calibration.";

VelocityParameters::usage="VelocityParameters[\"Lower\"] returns the velocity calibration parameters for the lower half of a MOSSPEC spectrum calculated by a call to Calibrate[]. VelocityParameters[\"Upper\"] returns the velocity calibration parameters for the upper half.\n\n"<>
"The returned value is a list of two numbers: {(Subscript[channel, 0]), slope}, where (Subscript[channel, 0]) is the zero-velocity channel number and slope is velocity per channel in mm/sec. If a velocity calibration has not been performed for the selected half of the spectrum, then the returned list is empty: {}. The conversion from channel number to velocity is given by: velocity = slope\[Cross](channel - (Subscript[channel, 0]))\n\n"<>
"VelocityParameters[\"SigmaLower\"] returns the uncertainties in the velocity calibration parameters returned by VelocityParameters[\"Lower\"], and similarly for VelocityParameters[\"SigmaUpper\"]";

VelocitySigma::usage="VelocitySigma[vel,\"Lower\"] and VelocitySigma[vel,\"Upper\"] return the uncertainty in velocity vel resulting from the uncertainties in the velocity calibration parameters of the lower and upper half velocity calibration performed by Calibrate[]. All velocities are in mm/sec.";

FitLines::usage="FitLines[6] attempts to fit a 6-line Mössbauer spectrum to the CurveFit data (Fe, (Subscript[Fe, 2])(Subscript[O, 3])).\n\n"<>
"FitLines[6, {list of 6 peak positions} ] attempts to fit a 6-line Mössbauer spectrum to the CurveFit data, using the list of initial peak positions provided as the 2nd argument.\n\n"<>
"FitLines[2] attempts to fit a 2-line Mössbauer spectrum to the CurveFit data ((Subscript[FeSO, 4])).\n\n"<>
"FitLines[2, {list of 2 peak positions} ] attempts to fit a 2-line Mössbauer spectrum to the CurveFit data, using the list of initial peak positions provided as the 2nd argument.\n\n"<>
"FitLines[1] attempts to fit a single-line Mössbauer spectrum to the CurveFit data (Stainless steel).\n\n"<>
"The X scale must be transformed to velocity using ChanToVelocity[] before calling this function.";

PeakList::usage="PeakList[] returns a list of the positions of the peaks found by a call to FitLines[].";

OtherParamsList::usage="OtherParamsList = other spectrum parameters from a 6-line fit in this order: {g,d,y1,y2,c,h}.";

PeakSigmaList::usage="PeakSigmaList[] returns a list of the sigmas of the positions of the peaks in PeakList[].";

NonlinFit::usage="NonlinFit[Data, Model, ParameterList, Variables, options] provides a Mathematica version-independent nonlinear regression routine. The argument list is the same as is used for NonlinearRegress[] (Mathematica v6) and NonlinearModelFit[] (Mathematica v7).";

ModelPeaks::usage="ModelPeaks[PeakList] fits a list of 6 peak positions (in mm/sec) to the theoretical model given by: \nFeLineFunction[(Subscript[v, iso]), (Subscript[v, B]), (Subscript[v, Q]), (Subscript[g, ratio])][(Subscript[I, z]), (SuperStar[(Subscript[I, z])])]\n\n"<>"ModelPeaks[PeakList,SigmaList] fits the data with a list of 6 peak uncertainties in SigmaList. The functions PeakList[] and PeakSigmaList[] may be used as arguments to ModelPeaks after you have fit a spectrum of 6 lines.";

Constants from Mathematica's Units and PhysicalConstants packages

E0::usage="E0 = The iron Mössbauer γ energy in eV. Mathematica code to get this value:\n\n"<>
"(10^3)(IsotopeData[\"Fe57\",\"ExcitedStateEnergies\"]//First)";
E0=14412.9;

cLight::usage="cLight = The speed of light in meters/sec. Mathematica code to get this value:\n\n"<>
"Needs[\"Units`\"]\n"<>
"Needs[\"PhysicalConstants`\"]\n"<>
"SpeedOfLight/(Meter/Second)";
cLight=299792458;

g::usage="g = The iron ground state gyromagnetic ratio. Mathematica code to get this value:\n\n"<>
"2 IsotopeData[\"Fe57\",\"MagneticMoment\"]";
g=0.1809;

μn::usage="μn = The Nuclear Magneton in eV/Tesla. Mathematica code to get this value:\n\n"<>
"Needs[\"Units`\"]\n"<>
"Needs[\"PhysicalConstants`\"]\n"<>
"Convert[NuclearMagneton,ElectronVolt/Tesla]";
μn=3.15245*(10^(-8));

Begin the Private package definitions and code

Begin["`Private`"];

(* clear everything first, so the definitions may be read multiple times *)
Unprotect["`*"]; 
ClearAll["`*"];

(* set plot lines to be thicker for visibility *)
SetOptions[Plot,{PlotStyle->{Thick,Blue}}];

Managing dialog windows

mosspecPaletteOpen=False;
If[TrueQ[Head[mosspecPaletteNotebook]===NotebookObject],NotebookClose[mosspecPaletteNotebook]];
If[TrueQ[Head[mosspecAskDialog]===NotebookObject],NotebookClose[mosspecAskDialog]];
mosspecHelpNotebook = Notebook[{}];
mosspecAskDialog=Notebook[{}];

Load an old MOSSPEC Data File

LoadMosspecFile[name_String]:=
Block[{dList,data,i,digits},
(* get a list of strings, one for each data line in the file *)
dList=Import[name,"List","IgnoreEmptyLines"->True];
digits=Table[ToString[i],{i,0,9}];
dList=Select[dList,!StringFreeQ[#,digits]&];
(* insert commas to delimit the numeric fields, wrap in { } *)
For[i=61, i>0, i-=8,dList=(StringInsert[#,",",i]&)/@dList];
dList=("{"<>#<>"}"&)/@dList;

(* now interpret each line as a list of values to get data array *)
data=1.0(ToExpression[#]&)/@dList;

(* drop the first column (X indices) and flatten to get just Y's *)
data=(data//Transpose//Rest//Transpose//Flatten);

(* The first element is the acquisition time. Build a header *)
CurveFit`SaveForUndo[];
CurveFit`DataFileName = name;
CurveFit`DataFileHeader = 
"S Mosspec Data File for Experiment 29 \n"<>
"S Acquisition Time "<>ToString[First[data]/1000.]<>" sec \n";

(* Now build the CurveFit data structure *)
CurveFit`yy=Rest[data];
CurveFit`n=Length[CurveFit`yy];
CurveFit`xx=Range[1,CurveFit`n];
CurveFit`sx=CurveFit`sy=ConstantArray[0,CurveFit`n];

Print["Imported "<>ToString[CurveFit`n]<>" values from file:"];
Print[CurveFit`DataFileName];
Print[];
Print["Comment Lines (generated automatically):"];
Print[CurveFit`DataFileHeader];
]

LoadMosspecFile[]:=LoadMosspecFile[CurveFit`DataFileName]

6-line Spectrum Line Velocity Calculations

FeLineIz[]={{-1/2,-3/2},{-1/2,-1/2},{-1/2,+1/2},{+1/2,-1/2},{+1/2,+1/2},{+1/2,+3/2}};

FeLineQ[]={1,-1,-1,-1,-1,1};

(* The natural iron line velocities used for calibration *)
FeLineV[]={-10.657/2,-6.167/2,-1.667/2,1.667/2,6.167/2,10.657/2};

(* General line velocity calculators *)
FeLineFunction[vIso_,vB_,vQ_,gRatio_][Iz_,IzStar_]:=
vIso+vB(Iz-gRatio IzStar)+vQ((IzStar^2)-5/4)
FeLineFunction[p__][{Iz_,IzStar_}]:=FeLineFunction[p][Iz,IzStar]

FeLineV[vIso_,vB_,vQ_,gRatio_]:=FeLineFunction[vIso,vB,vQ,gRatio]/@FeLineIz[]

Line Spectra

(* 6-line spectrum (in emission) *)
LineSpectrum[v:{_,_,_,_,_,_},y1_,y2_,gam_][x_]:=
((CurveFit`fLorentzianC[#,gam,1,0][x]&)/@v).{y1,y2,1,1,y2,y1}

(* 2-line spectrum (in emission) *)
LineSpectrum[{v1_,v2_},y2_,gam_][x_]:=
CurveFit`fLorentzianC[v1,gam,1,0][x]+CurveFit`fLorentzianC[v2,gam,y2,0][x]

(* single line spectrum (in emission) *)
LineSpectrum[{v_},gam_][x_]:=
CurveFit`fLorentzianC[v,gam,1,0][x]

LineSpectrum[v_,gam_][x_]:=
CurveFit`fLorentzianC[v,gam,1,0][x]

(* Mossbauer spectrum (in absorption) *)
MossSpectrum[LineSpectrum[p__],d_][x_]:=1-d LineSpectrum[p][x]

Baseline Correction for 1/(r^2) effect

BaselineCorrection[w_,h_][x_]:=h(1-(x^2)/(w^2))

Convert X scale to velocity

ChanToVelocity[x_]:=Which[
x<511&&CalLowerQ[],Last[VelocityParameters["Lower"]](x-First[VelocityParameters["Lower"]]),
x>=511&&CalUpperQ[],Last[VelocityParameters["Upper"]](x-First[VelocityParameters["Upper"]]),
True,x
]

ChanToVelocity[]:=(
If[!CurveFit`CheckLength[], Abort[]];
Which[
Mean[CurveFit`xx]<511&&CalLowerQ[],
CurveFit`SaveForUndo[];CurveFit`xx=Last[VelocityParameters["Lower"]](CurveFit`xx-First[VelocityParameters["Lower"]]);
CurveFit`sx*=Abs[Last[VelocityParameters["Lower"]]],

Mean[CurveFit`xx]>=511&&CalUpperQ[],
CurveFit`SaveForUndo[];CurveFit`xx=Last[VelocityParameters["Upper"]](CurveFit`xx-First[VelocityParameters["Upper"]]);
CurveFit`sx*=Abs[Last[VelocityParameters["Upper"]]];
{xx,sx,yy,sy}=({xx,sx,yy,sy}//Transpose//Sort//Transpose)
];
)

VelocitySigma[v_,"Lower"]:=If[CalLowerQ[],
Sqrt[(v^2)(((Last[VelocityParameters["SigmaLower"]]/Last[VelocityParameters["Lower"]]))^2)+(((Last[VelocityParameters["Lower"]]First[VelocityParameters["SigmaLower"]]))^2)]
]
VelocitySigma[v_,"Upper"]:=If[CalUpperQ[],
Sqrt[(v^2)(((Last[VelocityParameters["SigmaLower"]]/Last[VelocityParameters["Upper"]]))^2)+(((Last[VelocityParameters["Upper"]]First[VelocityParameters["SigmaUpper"]]))^2)]
]

Mössbauer spectra fitting functions

Nonlinear Regression function (depends on Mathematica version)

If[TrueQ[$VersionNumber≥ 7.0],
(* Version 7 nonlinear fits *)
NonlinFit[p__]:=Block[{model},
model=NonlinearModelFit[p];
{model["EstimatedVariance"],
model["BestFitParameters"],
model["ParameterErrors"]/Sqrt@model["EstimatedVariance"],model["CorrelationMatrix"]}
],
(* Version 6 nonlinear fits *)
NonlinFit[p__]:=Block[{model},
model=NonlinearRegress[p,RegressionReport->{BestFitParameters,AsymptoticCovarianceMatrix,EstimatedVariance,AsymptoticCorrelationMatrix}];
{(EstimatedVariance/.model),(BestFitParameters/.model),Sqrt@(Diagonal@First@(AsymptoticCovarianceMatrix/.model)/(EstimatedVariance/.model)),(AsymptoticCorrelationMatrix/.model)}
]
];

Velocity calibration using Natural Fe spectrum

CalLowerQ[]=False;
CalUpperQ[]=False;
VelocityParameters["Lower"]={};
VelocityParameters["Upper"]={};
VelocityParameters["SigmaLower"]={};
VelocityParameters["SigmaUpper"]={};

Calibrate::uncal="Calibrate only uncalibrated (raw) channel numbers.";
Calibrate::range="Limit the data to either the lower half or upper half of the full spectrum, then Calibrate.";

Calibrate[]:=
Block[{lowerQ,data,wts,results,f,init,x,c,m,s,d,y1,y2,g,h,sc,sm,ss,sd,sy1,sy2,sg,sh,sigs},
(* data checks *)
If[!CurveFit`CheckLength[], Abort[]];
If[Min[CurveFit`xx]<0||Max[CurveFit`xx]<50,
Message[Calibrate::uncal];Abort[]];
If[(Min[CurveFit`xx]<511-50)&&(Max[CurveFit`xx]>511+50),
Message[Calibrate::range];Abort[]];

(* fitting function *)
lowerQ=(Mean[CurveFit`xx]<511);
f:=c(MossSpectrum[LineSpectrum[FeLineV[],y1,y2,g],d][s(x-m)]+BaselineCorrection[9.0,h][s(x-m)]);
data=Transpose[{CurveFit`xx,CurveFit`yy}];
If[Min[CurveFit`sy]>0.0,wts=1/(CurveFit`sy^2)];
init={
{c,1.0Last[CurveFit`yy]},
{m,If[lowerQ,256,768]},
{s,If[lowerQ,9.0/256,-9.0/256]},
{d,.1},
{y1,2.0},
{y2,1.5},
{g,.3},
{h,If[lowerQ,-.01,.01]}
};

(* Print beginning info, then fit *)
CurveFit`Private`funct="Fit to "<>If[lowerQ,"lower","upper"]<>" half of a Natural Fe Mössbauer spectrum to determine channel velocity calibration of the MS-900 spectrometer.";
Print["n = ",CurveFit`n];
Print[CurveFit`Private`funct];
results=
If[Head[wts]===Symbol,
Print["Fit of (x,y)  (unweighted)"];
NonlinFit[data,f,init,x],
Print["Fit of (x,y±\!\(σ\_y\))  "];
NonlinFit[data,f,init,x,Weights->wts]
];
results=ReplacePart[results,3->Thread[Rule[{sc,sm,ss,sd,sy1,sy2,sg,sh},results[[3]](If[Head[wts]===Symbol,Sqrt[results[[1]]],1]) ]]];

(* Assign values to CurveFit functions resulting from the fit *)
Clear[CurveFit`fY,CurveFit`Private`yff,CurveFit`Private`ssy];
CurveFit`fY=Function[{x},Evaluate[f/.results[[2]]]];
CurveFit`Private`yff[x_]:=fY[x];
CurveFit`Private`ssy=CurveFit`sy;

(* Assign values to the appropriate velocity calibration results *)
If[lowerQ,
CalLowerQ[]=True;
VelocityParameters["Lower"]={m,s}/.results[[2]];
VelocityParameters["SigmaLower"]={sm,ss}/.results[[3]],

CalUpperQ[]=True;
VelocityParameters["Upper"]={m,s}/.results[[2]];
VelocityParameters["SigmaUpper"]={sm,ss}/.results[[3]]
];

CurveFit`Private`results=
Grid[{
If[Head[wts]===Symbol,{"Std. Deviation =",Sqrt[results[[1]]]},{"(χ^2) =",results[[1]]}],
{},
{Style["Velocity Calibration:",Bold]},
{"0-velocity channel =",m/.results[[2]],"±",sm/.results[[3]]},
{"slope (mm/s/channel) =",s/.results[[2]],"±",ss/.results[[3]]},
{},
{Style["Fe Spectrum Parameters:",Bold]},
{"Line width (mm/s) =",g/.results[[2]],"±",sg/.results[[3]]},
{"Inner line depth =",d/.results[[2]],"±",sd/.results[[3]]},
{"Depth outer/inner =",y1/.results[[2]],"±",sy1/.results[[3]]},
{"Depth middle/inner =",y2/.results[[2]],"±",sy2/.results[[3]]},
{"Baseline counts =",c/.results[[2]],"±",sc/.results[[3]]},
{"1/(r^2) correction factor =",h/.results[[2]],"±",sh/.results[[3]]}
},Alignment->{{Right,Right,Center,Left},Baseline}];
Print[CurveFit`Private`results];
]

Fit a Fe spectrum

FitLines::uncal="Perform a velocity calibration and then transform X from channels to velocity using ChanToVelocity[].";

FitLines[6,pInit_List:{}]:=
Block[{data,wts,results,f,init,thr,peaks,x,c,p1,p2,p3,p4,p5,p6,y1,y2,g,d,h,sc,sp1,sp2,sp3,sp4,sp5,sp6,sy1,sy2,sg,sd,sh,sigs},
(* data checks *)
If[!CurveFit`CheckLength[], Abort[]];
If[Max[CurveFit`xx]>50,
Message[FitLines::uncal];Abort[]];

(* identify peaks *)
If[Head[pInit]===List&&Length[pInit]==6,
peaks=pInit,
thr=.5(Max[CurveFit`yy]+Min[CurveFit`yy]);
peaks=Mean/@
FindClusters[
Select[Transpose[{CurveFit`xx,CurveFit`yy}],Last[#]<thr&][[All,1]],
6,Method->"Agglomerate"]
];

(* fitting function *)
f:=c(MossSpectrum[LineSpectrum[{p1,p2,p3,p4,p5,p6},y1,y2,g],d][x]+BaselineCorrection[9.0,h][x]);
data=Transpose[{CurveFit`xx,CurveFit`yy}];
If[Min[CurveFit`sy]>0.0,wts=1/(CurveFit`sy^2)];
init=Join[
{{c,1.0Last[CurveFit`yy]}},
Transpose[{{p1,p2,p3,p4,p5,p6},peaks}],
{{d,.1},
{y1,2.0},
{y2,1.5},
{g,.3},
{h,0.}}
];

(* Print beginning info, then fit *)
CurveFit`Private`funct="Fit 6-line Fe Mössbauer spectrum.";
Print["n = ",CurveFit`n];
Print[CurveFit`Private`funct];
results=
If[Head[wts]===Symbol,
Print["Fit of (x,y)  (unweighted)"];
NonlinFit[data,f,init,x],
Print["Fit of (x,y±\!\(σ\_y\))  "];
NonlinFit[data,f,init,x,Weights->wts]
];
results=ReplacePart[results,3->Thread[Rule[{sc,sp1,sp2,sp3,sp4,sp5,sp6,sd,sy1,sy2,sg,sh},results[[3]](If[Head[wts]===Symbol,Sqrt[results[[1]]],1]) ]]];
OtherParamsList={g,d,y1,y2,c,h}/.results[[2]];

(* Assign values to PeakList and PeakSigmaList *)
PeakList[]={p1,p2,p3,p4,p5,p6}/.results[[2]];
PeakSigmaList[]={sp1,sp2,sp3,sp4,sp5,sp6}/.results[[3]];

(* Assign values to CurveFit functions resulting from the fit *)
Clear[CurveFit`fY,CurveFit`Private`yff,CurveFit`Private`ssy];
CurveFit`fY=Function[{x},Evaluate[f/.results[[2]]]];
CurveFit`Private`yff[x_]:=fY[x];
CurveFit`Private`ssy=CurveFit`sy;

CurveFit`Private`results=
Column[{
Row[If[Head[wts]===Symbol,{"Std. Deviation = ",Sqrt[results[[1]]]},{"(χ^2) = ",results[[1]]}]],
,
Style["Peak Positions:",Bold],
Grid[{
{p1/.results[[2]],"±",sp1/.results[[3]],"  ",p2/.results[[2]],"±",sp2/.results[[3]]},
{p3/.results[[2]],"±",sp3/.results[[3]],"  ",p4/.results[[2]],"±",sp4/.results[[3]]},
{p5/.results[[2]],"±",sp5/.results[[3]],"  ",p6/.results[[2]],"±",sp6/.results[[3]]}
},Alignment->{{Right,Center,Left,Center,Right,Center,Left},Baseline}],
,
Style["Other Parameters:",Bold],
Grid[{
{"Line width (mm/s) =",g/.results[[2]],"±",sg/.results[[3]]},
{"Inner line depth =",d/.results[[2]],"±",sd/.results[[3]]},
{"Depth outer/inner =",y1/.results[[2]],"±",sy1/.results[[3]]},
{"Depth middle/inner =",y2/.results[[2]],"±",sy2/.results[[3]]},
{"Baseline counts =",c/.results[[2]],"±",sc/.results[[3]]},
{"1/(r^2) correction factor =",h/.results[[2]],"±",sh/.results[[3]]}
},Alignment->{{Right,Right,Center,Left},Baseline}],
,
"Use PeakList[] to return a list of peak positions.",
"Use PeakSigmaList[] to get their uncertainties."
}];
Print[CurveFit`Private`results];
]

FitLines[2,pInit_List:{}]:=
Block[{data,wts,results,f,init,thr,peaks,x,c,p1,p2,y2,g,d,h,sc,sp1,sp2,sy2,sg,sd,sh,sigs},
(* data checks *)
If[!CurveFit`CheckLength[], Abort[]];
If[Max[CurveFit`xx]>50,
Message[FitLines::uncal];Abort[]];

(* identify peaks *)
If[Head[pInit]===List&&Length[pInit]==2,
peaks=pInit,
thr=.5(Max[CurveFit`yy]+Min[CurveFit`yy]);
peaks=Mean/@
FindClusters[
Select[Transpose[{CurveFit`xx,CurveFit`yy}],Last[#]<thr&][[All,1]],
2,Method->"Agglomerate"]
];

(* fitting function *)
f:=c(MossSpectrum[LineSpectrum[{p1,p2},y2,g],d][x]+BaselineCorrection[9.0,h][x]);
data=Transpose[{CurveFit`xx,CurveFit`yy}];
If[Min[CurveFit`sy]>0.0,wts=1/(CurveFit`sy^2)];
init=Join[
{{c,1.0Last[CurveFit`yy]}},
Transpose[{{p1,p2},peaks}],
{{d,.1},
{y2,1.5},
{g,.3},
{h,0.}}
];

(* Print beginning info, then fit *)
CurveFit`Private`funct="Fit 2-line Fe Mössbauer spectrum.";
Print["n = ",CurveFit`n];
Print[CurveFit`Private`funct];
results=
If[Head[wts]===Symbol,
Print["Fit of (x,y)  (unweighted)"];
NonlinFit[data,f,init,x],
Print["Fit of (x,y±\!\(σ\_y\))  "];
NonlinFit[data,f,init,x,Weights->wts]
];
results=ReplacePart[results,3->Thread[Rule[{sc,sp1,sp2,sd,sy2,sg,sh},results[[3]](If[Head[wts]===Symbol,Sqrt[results[[1]]],1]) ]]];
Clear[OtherParamsList];

(* Assign values to PeakList and PeakSigmaList *)
PeakList[]={p1,p2}/.results[[2]];
PeakSigmaList[]={sp1,sp2}/.results[[3]];

(* Assign values to CurveFit functions resulting from the fit *)
Clear[CurveFit`fY,CurveFit`Private`yff,CurveFit`Private`ssy];
CurveFit`fY=Function[{x},Evaluate[f/.results[[2]]]];
CurveFit`Private`yff[x_]:=fY[x];
CurveFit`Private`ssy=CurveFit`sy;

CurveFit`Private`results=
Column[{
Row[If[Head[wts]===Symbol,{"Std. Deviation = ",Sqrt[results[[1]]]},{"(χ^2) = ",results[[1]]}]],
,
Style["Peak Positions:",Bold],
Grid[{
{p1/.results[[2]],"±",sp1/.results[[3]],"  ",p2/.results[[2]],"±",sp2/.results[[3]]}
},Alignment->{{Right,Center,Left,Center,Right,Center,Left},Baseline}],
,
Style["Other Parameters:",Bold],
Grid[{
{"Line width (mm/s) =",g/.results[[2]],"±",sg/.results[[3]]},
{"Lower line depth =",d/.results[[2]],"±",sd/.results[[3]]},
{"Depth upper/lower =",y2/.results[[2]],"±",sy2/.results[[3]]},
{"Baseline counts =",c/.results[[2]],"±",sc/.results[[3]]},
{"1/(r^2) correction factor =",h/.results[[2]],"±",sh/.results[[3]]}
},Alignment->{{Right,Right,Center,Left},Baseline}],
,
"Use PeakList[] to return a list of peak positions.",
"Use PeakSigmaList[] to get their uncertainties."
}];
Print[CurveFit`Private`results];
]

FitLines[1]:=
Block[{data,wts,results,f,init,peak,x,c,p1,g,d,h,sc,sp1,sg,sd,sh,sigs},
(* data checks *)
If[!CurveFit`CheckLength[], Abort[]];
If[Max[CurveFit`xx]>50,
Message[FitLines::uncal];Abort[]];

(* identify peaks *)
peak=
Mean[Select[
Transpose[{CurveFit`xx,CurveFit`yy}],
Last[#]≤Min[CurveFit`yy]&
][[All,1]] ];

(* fitting function *)
f:=c(MossSpectrum[LineSpectrum[{p1},g],d][x]+BaselineCorrection[9.0,h][x]);
data=Transpose[{CurveFit`xx,CurveFit`yy}];
If[Min[CurveFit`sy]>0.0,wts=1/(CurveFit`sy^2)];
init=Join[
{{c,1.0Last[CurveFit`yy]}},
{{p1,peak}},
{{d,.1},
{g,.3},
{h,0.}}
];

(* Print beginning info, then fit *)
CurveFit`Private`funct="Fit single line Fe Mössbauer spectrum.";
Print["n = ",CurveFit`n];
Print[CurveFit`Private`funct];
results=
If[Head[wts]===Symbol,
Print["Fit of (x,y)  (unweighted)"];
NonlinFit[data,f,init,x],
Print["Fit of (x,y±\!\(σ\_y\))  "];
NonlinFit[data,f,init,x,Weights->wts]
];
results=ReplacePart[results,3->Thread[Rule[{sc,sp1,sd,sg,sh},results[[3]](If[Head[wts]===Symbol,Sqrt[results[[1]]],1]) ]]];
Clear[OtherParamsList];

(* Assign values to PeakList and PeakSigmaList *)
PeakList[]={p1}/.results[[2]];
PeakSigmaList[]={sp1}/.results[[3]];

(* Assign values to CurveFit functions resulting from the fit *)
Clear[CurveFit`fY,CurveFit`Private`yff,CurveFit`Private`ssy];
CurveFit`fY=Function[{x},Evaluate[f/.results[[2]]]];
CurveFit`Private`yff[x_]:=fY[x];
CurveFit`Private`ssy=CurveFit`sy;

CurveFit`Private`results=
Column[{
Row[If[Head[wts]===Symbol,{"Std. Deviation = ",Sqrt[results[[1]]]},{"(χ^2) = ",results[[1]]}]],
,
Style["Peak Position:",Bold],
Grid[{
{p1/.results[[2]],"±",sp1/.results[[3]]}
},Alignment->{{Right,Center,Left},Baseline}],
,
Style["Other Parameters:",Bold],
Grid[{
{"Line width (mm/s) =",g/.results[[2]],"±",sg/.results[[3]]},
{"Line depth =",d/.results[[2]],"±",sd/.results[[3]]},
{"Baseline counts =",c/.results[[2]],"±",sc/.results[[3]]},
{"1/(r^2) correction factor =",h/.results[[2]],"±",sh/.results[[3]]}
},Alignment->{{Right,Right,Center,Left},Baseline}],
,
"Use PeakList[] to return a list of peak positions.",
"Use PeakSigmaList[] to get their uncertainties."
}];
Print[CurveFit`Private`results];
]

Fit peaks to theory

ModelPeaks[peaks:{_,_,_,_,_,_}]:=Block[{data,vIso,vB,vQ,gRatio,svIso,svB,svQ,sgRatio,I1,I2,results},
data=Flatten/@Transpose[{FeLineIz[],peaks}];
results=NonlinFit[data,
FeLineFunction[vIso,vB,vQ,gRatio][I1,I2],
{{vIso,0.0},{vB,4.0},{vQ,0.0},{gRatio,-1/1.75}},
{I1,I2}
];
results=ReplacePart[results,3->Thread[Rule[{svIso,svB,svQ,sgRatio},results[[3]]Sqrt[results[[1]]] ]]];

(* Assign values to CurveFit functions resulting from the fit *)
If[Head[OtherParamsList]===List && Length[OtherParamsList]==6,
Clear[CurveFit`fY,CurveFit`Private`yff,CurveFit`Private`ssy];
CurveFit`fY=Function[{x},
Evaluate[{g,d,y1,y2,c,h}=OtherParamsList;
c(MossSpectrum[LineSpectrum[(FeLineV[vIso,vB,vQ,gRatio]/.results[[2]]),y1,y2,g],d][x]+BaselineCorrection[9.0,h][x])]];
CurveFit`Private`yff[x_]:=fY[x];
CurveFit`Private`ssy=CurveFit`sy;
];

CurveFit`Private`results=
Column[{
Row[{"Std. Deviation = ",Sqrt[results[[1]]]}],
,
Style["Model Parameters:",Bold],
Grid[{
{"Iso shift ((Subscript[v, iso]),mm/s) =",Chop[vIso/.results[[2]]],"±",svIso/.results[[3]]},
{"Hyperfine ((Subscript[v, B]),mm/s) =",vB/.results[[2]],"±",svB/.results[[3]]},
{"(SuperStar[g])/g =",gRatio/.results[[2]],"±",sgRatio/.results[[3]]},
{"-1/((SuperStar[g])/g) =",-1/gRatio/.results[[2]],"±",sgRatio/gRatio^2/.Join[results[[2]],results[[3]]]},
{"Quadrupole ((Subscript[v, Q]),mm/s) =",Chop[vQ/.results[[2]]],"±",svQ/.results[[3]]}
},Alignment->{{Right,Right,Center,Left},Baseline},Frame->True],
,
Style["Peak Positions:",Bold],
Grid[{
Prepend[peaks,"Data:"],
Prepend[
(FeLineV[vIso,vB,vQ,gRatio]/.results[[2]]),"Model:"]
},Alignment->Left]
}];
Print[CurveFit`Private`results];
If[Head[OtherParamsList]===List && Length[OtherParamsList]==6,
Print["Plot Fit Results to see modeled spectrum v. counts data."]
];
]

ModelPeaks[peaks:{_,_,_,_,_,_},sigs:{_,_,_,_,_,_}]:=Block[{data,vIso,vB,vQ,gRatio,svIso,svB,svQ,sgRatio,I1,I2,results,g,d,y1,y2,c,h},
data=Flatten/@Transpose[{FeLineIz[],peaks}];
results=NonlinFit[data,
FeLineFunction[vIso,vB,vQ,gRatio][I1,I2],
{{vIso,0.0},{vB,4.0},{vQ,0.0},{gRatio,-1/1.75}},
{I1,I2},Weights->1/sigs^2
];
results=ReplacePart[results,3->Thread[Rule[{svIso,svB,svQ,sgRatio},results[[3]] ]]];

(* Assign values to CurveFit functions resulting from the fit *)
If[Head[OtherParamsList]===List && Length[OtherParamsList]==6,
Clear[CurveFit`fY,CurveFit`Private`yff,CurveFit`Private`ssy];
CurveFit`fY=Function[{x},
Evaluate[{g,d,y1,y2,c,h}=OtherParamsList;
c(MossSpectrum[LineSpectrum[(FeLineV[vIso,vB,vQ,gRatio]/.results[[2]]),y1,y2,g],d][x]+BaselineCorrection[9.0,h][x])]];
CurveFit`Private`yff[x_]:=fY[x];
CurveFit`Private`ssy=CurveFit`sy;
];

CurveFit`Private`results=
Column[{
Row[{"(χ^2) = ",results[[1]]}],
,
Style["Model Parameters:",Bold],
Grid[{
{"Iso shift ((Subscript[v, iso]),mm/s) =",Chop[vIso/.results[[2]]],"±",svIso/.results[[3]]},
{"Hyperfine ((Subscript[v, B]),mm/s) =",vB/.results[[2]],"±",svB/.results[[3]]},
{"(SuperStar[g])/g =",gRatio/.results[[2]],"±",sgRatio/.results[[3]]},
{"-1/((SuperStar[g])/g) =",-1/gRatio/.results[[2]],"±",sgRatio/gRatio^2/.Join[results[[2]],results[[3]]]},
{"Quadrupole ((Subscript[v, Q]),mm/s) =",Chop[vQ/.results[[2]]],"±",svQ/.results[[3]]}
},Alignment->{{Right,Right,Center,Left},Baseline},Frame->True],
,
Style["Peak Positions:",Bold],
Grid[{
Prepend[peaks,"Data:"],
Prepend[
(FeLineV[vIso,vB,vQ,gRatio]/.results[[2]]),"Model:"]
},Alignment->Left]
}];
Print[CurveFit`Private`results];
If[Head[OtherParamsList]===List && Length[OtherParamsList]==6,
Print["Plot Fit Results to see modeled spectrum v. counts data."]
];
]

Choice Manager

mosspecPalette:=If[!mosspecPaletteOpen,
mosspecPaletteOpen=True;
NotebookClose[mosspecAskDialog];
NotebookClose[mosspecHelpNotebook];
mosspecAskDialog=mosspecAsk;
mosspecPaletteNotebook=CreatePalette[
Grid[
{
{Style["Data Management",Bold],SpanFromLeft},
{
Button[Tooltip["Load Mosspec Data File","Load a MOSSPEC program data file"],
fileloader],SpanFromLeft
},
{
Button[Tooltip["Set to Poisson count data","Set σy = \!\(\*SqrtBox[\"Y\"]\)"],
eval["MakePoisson[Xerrors -> False]"]],SpanFromLeft
},
{
Button[Tooltip["Keep an X-range","Opens a dialog to remove all but a subset of the data"],
subsetX],SpanFromLeft
},
{
Button[Tooltip["Rebin Data","Merge data in adjacent bins to reduce noise."],
rebinner],SpanFromLeft
},
{
Button[Tooltip["Undo","Undo last data operation"],
eval["Undo[]"]],
Button[Tooltip["Redo","Redo last Undo"],
eval["Redo[]"]]
},

{Style["Velocity Calibration",Bold],SpanFromLeft},
{
Button[Tooltip["Calibrate","Calibrate velocities using a Natural Iron spectrum"],
eval["Calibrate[]"]],SpanFromLeft
},
{
Button[Tooltip["Convert to Velocity","Use calibration results to convert X scale to velocity"],
eval["ChanToVelocity[]"]],SpanFromLeft
},

{Style["Fit Lines",Bold],SpanFromLeft},
{
Button[Tooltip["Fit 6 Lines","Fit a 6-line spectrum. Calibrate velocity first."],
eval["FitLines[6]"]],SpanFromLeft
},
{
Button[Tooltip["Fit 2 Lines","Fit a 2-line spectrum. Calibrate velocity first."],
eval["FitLines[2]"]],SpanFromLeft
},
{
Button[Tooltip["Fit 1 Line","Fit a single line spectrum. Calibrate velocity first."],
eval["FitLines[1]"]],SpanFromLeft
},
{
Button[Tooltip["Plot Fit Results","Linear plot of fit results."],
eval["LinearDifferencePlot[]"]],SpanFromLeft
},


{Style["Fit Model to Lines",Bold],SpanFromLeft},
{
Button[Tooltip["Model 6 Lines","Once you've fit a 6-line spectrum, calculate model parameters from line positions."],
eval["ModelPeaks[PeakList[],PeakSigmaList[]]"]],SpanFromLeft
},

{Style["Legacy MOSSPEC data",Bold],SpanFromLeft},
{
Button[Tooltip["Load Old Mosspec File","Load an old-format MOSSPEC program data file"],
oldfileloader],SpanFromLeft
},

{Style["Application Control",Bold],SpanFromLeft},
{
Button[
Tooltip["Help","Opens a notebook with information"],Quiet[NotebookClose[mosspecHelpNotebook]];
mosspecHelpNotebook=
CreateDocument[{Cell[BoxData[{RowBox[{RowBox[{"Print","[",RowBox[{"Style","[",RowBox[{"\"\<Click on a symbol name for its description\>\"",",","\"\<Subsubsection\>\""}],"]"}],"]"}],";"}],"
",RowBox[{"Information","[",RowBox[{"\"\<Mosspec`*\>\"",",",RowBox[{"LongForm","->","False"}]}],"]"}]}],"Input"]

},
WindowTitle->"Mosspec Help",Saveable->False
];
SetSelectedNotebook[mosspecHelpNotebook];
SelectionMove[mosspecHelpNotebook,All,CellContents];
SelectionEvaluateCreateCell[mosspecHelpNotebook],
ImageSize->70
],
Button[Tooltip["Quit","Quit the Mössbauer spectrum application"],
mosspecPaletteOpen=False;
SetOptions[mosspecAskDialog,Visible->True];NotebookClose[mosspecHelpNotebook];
NotebookClose[mosspecPaletteNotebook],
ImageSize->70
]
}
},
Dividers->{False,{False,False,False,False,False,False,True,False,False,True,False,False,False,False,True,False,True,False,True}}
],
WindowTitle->"Mosspec",
Saveable->False,WindowFrameElements->{},
WindowMargins->{{50,Automatic},{Automatic,10}}
]]

mosspecAsk:=CreateDialog[Grid[{{TextCell["Are you sure you want to end the\nMössbauer spectrum application?"],SpanFromLeft},{DefaultButton[ImageSize->80],CancelButton[mosspecPalette;SetOptions[mosspecAskDialog,Visible->False],ImageSize->80]}}],WindowFloating->True,Visible->False,WindowTitle->"Quit Mosspec"]

fileloader:=eval[
		"With["<>
		"  {name = SystemDialogInput[\"FileOpen\","<>
			"{DataFileName,{\"data files\"->"<>
			"{\"*.dat\"},\"all files\"->{\"*\"}}}]"<>
			"},\n"<>
		"  If[ name =!= $Canceled, \n"<>
		"    LoadFile[name]]\n"<>
		"]"
		]

oldfileloader:=eval[
"With["<>"  {name = SystemDialogInput[\"FileOpen\",CurveFit`DataFileName]},\n"<>"  If[ name =!= $Canceled, \n"<>
"    LoadMosspecFile[name]]\n"<>
"]"
]

subsetX:=CurveFit`Palette`Private`XReval["XRangeKeep",SymbolName[CurveFit`Palette`Private`plotfcn],SymbolName[CurveFit`plotLogs[CurveFit`Palette`Private`plotfcn][[1]]],"\"Set the X values for the range you wish to keep.\""]

rebinner := insert[
"(* change the arguments of the following function \n"<>"   to perform the data rebinning you want, \n"<>"   then evaluate this cell.      *)\n\n"<>
"Rebin[ \n"<>"  3,   (* change to the number of channels you want to combine *)\n"<>"  Averaging -> True,   (* change to False to make Y values the total of the binned channel Y's *)\n"<>
"  Xerrors -> False \n"<>"]\n\n"<>
"(* Use Undo[] if you don't like the results. *)"
]

insert[command_String, nb_NotebookObject]:=(
SelectionMove[nb,After,Notebook];
NotebookWrite[nb,command,None];
SelectionMove[nb, All, Cell];
SetOptions[NotebookSelection[nb], 
	GeneratedCell -> False];
)
insert[command_String]:= insert[command, InputNotebook[]]

eval[command_String, nb_NotebookObject]:=(
SelectionMove[nb,After,Notebook];
NotebookWrite[nb,command,None];
SelectionMove[nb, All, Cell];
SetOptions[NotebookSelection[nb], 
	GeneratedCell -> False];
SelectionEvaluateCreateCell[nb];
)
eval[command_String]:= eval[command, InputNotebook[]]

TextHint[s_String, nb_NotebookObject]:=(
SelectionMove[nb,After,Notebook];
NotebookWrite[nb,Cell[s,"Text"],None];
SelectionMove[nb, All, Cell];
SetOptions[NotebookSelection[nb], 
	GeneratedCell -> False];
)
TextHint[s_String]:=TextHint[s,InputNotebook[]]

End[];
EndPackage[];

Start the appication

mosspecPalette;

Some Useful Constant Definitions

User Area