(*LASr - April 11 2002 - Dauphas et al. 2022 - dauphas@uchicago.edu*)
mrb85=84.91178974;
mrb87=86.90918053;
msr84=83.913419;
msr86=85.90926073;
msr87=86.90887750;
msr88=87.90561226;
rb8785=0.2783/0.7217; (*87Rb/85Rb ratio*)
sr8886=8.375209; (*88Sr/86Sr ratio of SRM987 for internal normalization *)

(*uncomment to provide here the path to the starting directory*)
(*SetDirectory["/Users/daxuphasu/Desktop/VisitBremen_2022/DATA_FILES_NOTES_CODE"]*)

(*Select the method to handle reduction of transient signal produced by laser ablation*)
mtd=ChoiceDialog["Select a method for handling LA data ",{"Average counts"->"m1","Peak Integration"->"m2","Extrapolation to zero time"->"m3"}]

(*user navigates the directories to find the data file and import it*)
file=SystemDialogInput["FileOpen",Directory[]];
mf=Import[file];
dir=DirectoryName[file]; (*recover the path of the imported file*)
SetDirectory[dir]; (*set the working directory to the path of the imported file*)
fl=FileBaseName[file]; (*recover the file name without the extension*)
 efil=FileNameJoin[{dir,fl<>"_export"<>".xlsx"}]; (*defines the filename for exported data*)
mf=mf[[1]]; (*table containing the data*)

(*this dialog asks the user to define the entries where the time, label, and instensities are using dropdown menus*)
t1=Table[i,{i,1,Dimensions[mf[[{1,2,3,4,5},All]]][[2]]}]; (*only the first 5 entries are displayed for practical reasons*)
t1c=Join[{0},t1];
p0=1;p1=2;p2=3;p3=0;p4=0;p5=0;p6=0;p7=0;p8=0;p9=0;p10=0;
resb=DialogInput[{TextCell["Select the Rb-Sr data (0=not available; leave 85Rb=0 if this is a pure Sr analysis)."],Manipulate[Grid[Transpose[Join[{t1},mf[[1;;5,All]]]],Frame->All,Background->{1->Pink,{p0->Gray,p1->Gray,p2->Gray,p3->Gray,p4->Gray,p5->Gray,p6->Gray,p7->Gray,p8->Gray,p9->Gray,p10->Gray}}],{{p0,1,"Index"},t1c,ControlType->PopupMenu},{{p1,1,"Timestamp"},t1c,ControlType->PopupMenu},{{p2,1,"Label"},t1c,ControlType->PopupMenu},{{p3,1,"84SrF"},t1c,ControlType->PopupMenu},{{p4,1,"85RbF"},t1c,ControlType->PopupMenu},{{p5,1,"86SrF"},t1c,ControlType->PopupMenu},
{{p6,1,"87SrF+RbF"},t1c,ControlType->PopupMenu},
{{p7,1,"88SrF"},t1c,ControlType->PopupMenu},
{{p8,1,"85Rb"},t1c,ControlType->PopupMenu},
{{p9,1,"87Rb+Sr"},t1c,ControlType->PopupMenu},
{{p10,1,"88Sr"},t1c,ControlType->PopupMenu},LocalizeVariables->False,SynchronousUpdating->True],Button["Proceed",DialogReturn[{p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10}]]}];

thod=If[resb[[9]]==0,"Sronly","RbSr"]; (*If the user left 85Rb at 0, then it is a Sr only measurement*)

(*this dialog asks the user to define the header*)
q1=2;
resc=DialogInput[{TextCell["Select the last column of header."],Manipulate[Grid[Join[{{1,2,3,4,5}},Transpose[mf[[1;;5,All]]]],Frame->All,Background->{None,1->Pink,{{{1,Max[t1c]+1},{1,q1}}->Gray}}],{{q1,1,"Header columns"},{1,2,3,4,5},ControlType->PopupMenu},LocalizeVariables->False,SynchronousUpdating->True],Button["Proceed",DialogReturn[{q1}]]}];

gd=Select[resb,#≠0&]; (*if 0s are left, it means that it is a Sr-only measurement and the Rb entries are eliminated*)

(*1=index,2=time,3=label,4=84SrF,5=85RbF,6=86SrF,7=87SrF,8=88SrF,9=85Rb,10=87Sr,11=88Sr -columns 9 to 11 are not in the table if it is a Sr-only measurement*)
mfb=mf[[All,gd]]; (*select the columns that we are interested in *)

mfb=mfb[[resc[[1]]+1;;,All]]; (*remove the header*)
tic=AbsoluteTime[#]&/@mfb[[All,2]];(*vector containing the time of each analysis. The date time stamp is converted to seconds*)
tic=tic-tic[[1]];
mfb[[All,2]]=tic;


(*{1"Index",2"Timestamp",3"Label",4"Static 2:85Rb",5"Static 2:87Sr",6"Static 2:88Sr"}*)
If[thod=="RbSr",rb=mfb[[All,{1,2,3,9,10,11}]]]; (*extract the columns containing Rb*)
(*{1"Index",2"Timestamp",3"Label",4"Static 1:84Sr.F",5"Static 1:85Rb.F",6"Static 1:86Sr.F",7"Static 1:87Sr.F",8"Static 1:88Sr.F"}*)
sr=mfb[[All,1;;8]]; (*extract the columns containing Sr*)

If[thod=="RbSr",rb=Select[rb,#[[4]]!=""&]]; (*remove all empty lines for Rb*)
sr=Select[sr,#[[4]]!=""&]; (*remove all empty lines for Sr*)

nm=DeleteDuplicates[sr[[All,1]]]; (*give the names of all distinct spot labels*)

nc=Dimensions[nm][[1]]; (*number of data points*)

nsr=nm; (*this is a vector that will contain in each cell the matrix of data corresponding to Sr for a particular sample*)
If[thod=="RbSr",nrb=nm]; (*this is a vector that will contain in each cell the matrix of data corresponding to Rb for a particular sample*)
For[i=1,i≤nc,i++,
nsr[[i]]=Select[sr,#[[1]]==nm[[i]]&]; 
If[thod=="RbSr",nrb[[i]]=Select[rb,#[[1]]==nm[[i]]&]]; 
]

lg[x_,v_]:=Module[{out},out=Nearest[v,x][[1]];out]; (*function that finds the nearest value to x in vector v*)

(*datared is a function that does the data reduction for a sample*)
datared[msr_(*sr data for a particular sample*),mrb_(*rb data for a particular sample*)]:=
Module[{timv=Sort[Join[msr[[All,2]],
mrb[[All,2]],1]],brb,bsr,drb,dsr},
idx=msr[[1,1]];
nam=msr[[1,3]];
zero=Min[timv];
rg=Max[timv];
(*use default values here that correspond to the boundaries in each sample*)
l1=zero+26.6;l2=zero+31.6;l3=zero+77; 
(* mrb={1"Timestamp",2"Label","",3"Static 2:85Rb",4"Static 2:87Sr",5"Static 2:88Sr"}*)
(*msr={1"Timestamp",2"Label","",3"Static 1:84Sr.F",4"Static 1:85Rb.F",5"Static 1:86Sr.F",6"Static 1:87Sr.F",7"Static 1:88Sr.F"}*)
(*the user uses interactive sliders to move the boundaries of the acquisition*)
res=DialogInput[{TextCell["Set the limits for OPZ, start, end. Then click Proceed."],Manipulate[ListPlot[{Transpose[{msr[[All,2]],Log[msr[[All,8]]]}],If[thod=="RbSr",Transpose[{mrb[[All,2]],Log[mrb[[All,4]]]}],{,}]},PlotRange->Full,AxesLabel->{"time (sec)","cps; Red=Sr, Blue=Rb"},PlotLabel->StringForm["Index=`1`, Name=`2`, `3`% completed.",idx,nam,NumberForm[1.0*o/nc*100,{3,1}]],PlotStyle->{Red,Blue},GridLines->{{{lg[l1,timv],Red},{lg[l2,timv],Blue},{lg[l3,timv],Blue}},None}],{l1,zero,rg,Appearance->"Open"},{l2,zero,rg,Appearance->"Open"},{l3,zero,rg,Appearance->"Open"},LocalizeVariables->False,SynchronousUpdating->True],Button["Proceed",DialogReturn[{lg[l1,timv],lg[l2,timv],lg[l3,timv]}]]}];

If[thod=="RbSr",brb=Select[mrb,#[[2]]≤res[[1]]&]];(*table that contains the blank of rb*)
bsr=Select[msr,#[[2]]≤res[[1]]&];(*table that contains the blank of sr*)
If[thod=="RbSr",drb=Select[mrb,res[[2]]<=#[[2]]≤res[[3]]&]];(*table that contains the data of rb*)
dsr=Select[msr,res[[2]]<=#[[2]]≤res[[3]]&];(*table that contains the data of sr*)

If[thod=="RbSr",brbt={Length[brb]*idx,0,Length[brb]*nam,Total[brb[[All,4]]],Total[brb[[All,5]]],Total[brb[[All,6]]]}/Length[brb]]; (*the cps on the blank are averaged for Rb*)
bsrt={Length[bsr]*idx,0,Length[bsr]*nam,Total[bsr[[All,4]]],Total[bsr[[All,5]]],Total[bsr[[All,6]]],Total[bsr[[All,7]]],Total[bsr[[All,8]]]}/Length[bsr]; (*the cps on the blank are averaged for Sr*)

If[thod=="RbSr",crb=Length[drb]];
csr=Length[dsr];

If[thod=="RbSr",drb=drb-ConstantArray[brbt,Length[drb]]]; (*remove blank cps from Rb; a column of 0 is added to the blank so the time is not affected by removal of the date to another date)*)
If[thod=="RbSr",For[i=1,i<=Length[drb],drb[[i,1]]=idx;drb[[i,3]]=nam;i++]]; (*puts back the name that was deleted*)
dsr=dsr-ConstantArray[bsrt,Length[dsr]]; (*remove blank cps from Sr*)
For[i=1,i<=Length[dsr],dsr[[i,3]]=nam;dsr[[i,1]]=idx;i++]; (*puts back the name that was deleted*)

(*If the method uses mean cps*)
If[mtd=="m1",
(*format rbg {1"index",2"Label",3"Static 2:85Rb",4"Static 2:87Sr",5"Static 2:88Sr"}*)
If[thod=="RbSr",rbg={Length[drb]*idx,Length[drb]*nam,Total[drb[[All,4]]],Total[drb[[All,5]]],Total[drb[[All,6]]]}/Length[drb]];
(*format srg {1"Index",2"Label",3"Static 1:84Sr.F",4"Static 1:85Rb.F",5"Static 1:86Sr.F",6"Static 1:87Sr.F",7"Static 1:88Sr.F"}*)
srg={Length[drb]*idx,Length[dsr]*nam,Total[dsr[[All,4]]],Total[dsr[[All,5]]],Total[dsr[[All,6]]],Total[dsr[[All,7]]],Total[dsr[[All,8]]]}/Length[dsr];
If[thod=="RbSr",rb85b=(rbg[[4]]-rbg[[5]]*7/82.58)/rb8785];
If[thod=="RbSr",rb85=rbg[[3]]];
(*krf83=srg[[2]];*)
srf84=srg[[3]];
rbf85=srg[[4]];
srf86=srg[[5]];
srf87=srg[[6]];
srf88=srg[[7]]
];

(*integration of the peak over the same time window*)
(*interpolation order of 1 means that the data are extrapolated using linear segments*)
If[mtd=="m2",
(*format drb {1"index",2"Timestamp",3"Label","",4"Static 2:85Rb",5"Static 2:87Sr",6"Static 2:88Sr"}*)
If[thod=="RbSr",frb85=Interpolation[drb[[All,{2,4}]],Method->"Spline",InterpolationOrder->1]]; (*2 column matrix time,85Rb*)
If[thod=="RbSr",fs87b=Interpolation[drb[[All,{2,5}]],Method->"Spline",InterpolationOrder->1]];
If[thod=="RbSr",fs88b=Interpolation[drb[[All,{2,6}]],Method->"Spline",InterpolationOrder->1]];
(*format dsr {1"Index",2"Timestamp",3"Label",4"Static 1:84Sr.F",5"Static 1:85Rb.F",6"Static 1:86Sr.F",7"Static 1:87Sr.F",8"Static 1:88Sr.F"}*)
fsrf84=Interpolation[dsr[[All,{2,4}]],Method->"Spline",InterpolationOrder->1];
frbf85=Interpolation[dsr[[All,{2,5}]],Method->"Spline",InterpolationOrder->1];
fsrf86=Interpolation[dsr[[All,{2,6}]],Method->"Spline",InterpolationOrder->1];
fsrf87=Interpolation[dsr[[All,{2,7}]],Method->"Spline",InterpolationOrder->1];
fsrf88=Interpolation[dsr[[All,{2,8}]],Method->"Spline",InterpolationOrder->1];

If[thod=="RbSr",s87b=NIntegrate[fs87b[x],{x,res[[2]],res[[3]]}]];
If[thod=="RbSr",s88b=NIntegrate[fs88b[x],{x,res[[2]],res[[3]]}]];
If[thod=="RbSr",rb85b=(s87b-s88b*7/82.58)/rb8785];
If[thod=="RbSr",rb85=NIntegrate[frb85[x],{x,res[[2]],res[[3]]}]];
srf84=NIntegrate[fsrf84[x],{x,res[[2]],res[[3]]}];
rbf85=NIntegrate[frbf85[x],{x,res[[2]],res[[3]]}];
srf86=NIntegrate[fsrf86[x],{x,res[[2]],res[[3]]}];
srf87=NIntegrate[fsrf87[x],{x,res[[2]],res[[3]]}];
srf88=NIntegrate[fsrf88[x],{x,res[[2]],res[[3]]}];
];

If[mtd=="m3", (*extrapolate to time zero using parabola*)
(*format drb {1"index",2"Timestamp",3"Label","",4"Static 2:85Rb",5"Static 2:87Sr",6"Static 2:88Sr"}*)
If[thod=="RbSr",frb85=Fit[drb[[All,{2,4}]],{1,x,x^2},x]]; (*2 column matrix time,85Rb*)
If[thod=="RbSr",fs87b=Fit[drb[[All,{2,5}]],{1,x,x^2},x]];
If[thod=="RbSr",fs88b=Fit[drb[[All,{2,6}]],{1,x,x^2},x]];
(*format dsr {1"Index",2"Timestamp",3"Label",4"Static 1:84Sr.F",5"Static 1:85Rb.F",6"Static 1:86Sr.F",7"Static 1:87Sr.F",8"Static 1:88Sr.F"}*)
fsrf84=Fit[dsr[[All,{2,4}]],{1,x,x^2},x];
frbf85=Fit[dsr[[All,{2,5}]],{1,x,x^2},x];
fsrf86=Fit[dsr[[All,{2,6}]],{1,x,x^2},x];
fsrf87=Fit[dsr[[All,{2,7}]],{1,x,x^2},x];
fsrf88=Fit[dsr[[All,{2,8}]],{1,x,x^2},x];

If[thod=="RbSr",s87b=fs87b/.x->res[[2]]];
If[thod=="RbSr",s88b=fs88b/.x->res[[2]]];
If[thod=="RbSr",rb85b=(s87b-s88b*7/82.58)/rb8785];
If[thod=="RbSr",rb85=frb85/.x->res[[2]]];
srf84=fsrf84/.x->res[[2]];
rbf85=frbf85/.x->res[[2]];
srf86=fsrf86/.x->res[[2]];
srf87=fsrf87/.x->res[[2]];
srf88=fsrf88/.x->res[[2]];
];


(*{1"Index",2"Timestamp",3"Label",4"Static 1:84Sr.F",5"Static 1:85Rb.F",6"Static 1:86Sr.F",7"Static 1:87Sr.F",8"Static 1:88Sr.F"}*)


(*calculation of errors from standard deviation of cycles; the c suffix stands for cycle*)
(*rb85b is rb85 calculated from 87 after subtraction of 87Sr based on 88Sr*)
If[thod=="RbSr",frb85bc=Interpolation[(drb[[All,{2,5}]]-drb[[All,{2,6}]]*7/82.58)/rb8785,Method->"Spline",InterpolationOrder->1]]; (*interpolationorder of 1 means that they are interpolates using straight lines*)
If[thod=="RbSr",i85rbbc=frb85c[#]&/@dsr[[All,2]]]; (*the 85Rb is interpolated to the time when Sr is measured for calculation of errors*)
If[thod=="RbSr",frb85c=Interpolation[drb[[All,{2,4}]],Method->"Spline",InterpolationOrder->1]];
If[thod=="RbSr",i85rbc=frb85c[#]&/@dsr[[All,2]]]; (*the 85Rb is interpolated to the time when Sr is measured for calculation of errors*)
bc=Log[(dsr[[All,8]]/dsr[[All,6]])/sr8886]/Log[msr88/msr86];
i84src=dsr[[All,4]];
i86src=dsr[[All,6]];
i87src=dsr[[All,7]]-dsr[[All,5]]*rb8785*(mrb87/mrb85)^bc;
i88src=dsr[[All,8]];
r8486nc=(i84src/i86src)/(msr84/msr86)^bc;
r8786nc=(i87src/i86src)/(msr87/msr86)^bc;
If[thod=="RbSr",rb87bsr86c=i85rbbc*rb8785/i86src]; (*based on 85Rb calculated from the signal at 87*)
If[thod=="RbSr",rb87sr86c=i85rbc*rb8785/i86src];
e8786N=StandardDeviation[r8786nc]/Length[r8786nc]^0.5; (*stdev of mean; std error*)
e8486N=StandardDeviation[r8486nc]/Length[r8486nc]^0.5;
If[thod=="RbSr",erb87sr86=StandardDeviation[rb87sr86c]/Length[rb87sr86c]^0.5];
If[thod=="RbSr",erb87bsr86=StandardDeviation[rb87bsr86c]/Length[rb87bsr86c]^0.5];
If[thod=="RbSr",rho=Correlation[rb87sr86c,r8786nc]];

b=Log[(srf88/srf86)/sr8886]/Log[msr88/msr86];
delta88=((srf88/srf86)/sr8886-1)*1000;
i84sr=srf84;
i86sr=srf86;
i87sr=srf87-rbf85*rb8785*(mrb87/mrb85)^b;
i88sr=srf88;
If[thod=="RbSr",i85rb=rb85];
If[thod=="RbSr",i85rbb=rb85b];

r8486n=(i84sr/i86sr)/(msr84/msr86)^b;
r8786n=(i87sr/i86sr)/(msr87/msr86)^b;
delta84n=(r8486n/0.05655-1)*1000;
delta87n=(r8786n/0.709699-1)*1000;
If[thod=="RbSr",rb87sr86=i85rb*rb8785/i86sr];
If[thod=="RbSr",rb87bsr86=i85rbb*rb8785/i86sr];

{idx,nam,If[thod=="RbSr",crb,"NA"],csr,If[thod=="RbSr",i85rb,"NA"],i88sr,r8486n,(*InverseCDF[StudentTDistribution[csr-1],0.975]*)2*e8486N,If[thod=="RbSr",rb87sr86,"NA"],If[thod=="RbSr",2*erb87sr86,"NA"],If[thod=="RbSr",rho,"NA"],If[thod=="RbSr",rb87bsr86,"NA"],If[thod=="RbSr",2*erb87bsr86,"NA"],r8786n,2*e8786N,delta87n,delta88,delta84n}
]

istab=Table[0,{j,1,nc},{i,1,19}]; (*table where the results will be exported*)
For[o=0,o<nc,o++;
isres=datared[nsr[[o]],If[thod=="RbSr",nrb[[o]],{}]];
istab[[o]]=isres;
Print[o,isres];
]

(*index,label,delta88,r8486n,delta84n,r8786n,delta87n,rb87sr86}*)
istab=Prepend[istab,{"1.Index","2.Label","3.# pt Rb","4.# pt Sr","5.Tot cts 85Rb","6.Tot cts 88Sr","7.84Sr/86Sr norm","8.2sigma mean 84Sr/86Sr norm","9.87Rb/86Sr (from 85Rb signal)","10.2sigma mean 87Rb/86Sr","11.Rho (correlation coefficient between 87Sr/86Sr norm and 87Rb/86Sr)","12.rb87Bsr86 (from 87Rb=87i-88i*87Sr/88Sr)","13.95%CI 87Rb/86Sr (from 87i and 88i)","14.87Sr/86Sr norm","15.2sigma mean 87Sr/86Sr norm","16.d87Sr norm (permil, relative to SRM610 0.709699)","17.d88Sr (permil, relative to SRM987=8.37861)","18.d84Sr norm (permil relative to 0.05655)"}];
ListPlot[Transpose[{istab[[All,11]],istab[[All,13]]}]]
Export[efil,istab];