Clear["Global`*"];

(* import data from Desktop *)
path="C:\\Users\\evansr2\\";
SetDirectory[path<>"\\Desktop"];
FileNames[];
Data=Import["Improving the Quantitative Interpretation of diffusion-ordered NMR spectra - original data.xlsx"];

(*setting up some of the parameters of the data set and some constants*)
ndp=Length[data];
na=6.022 10^23;
kb=8.3144/na;
t=298.15;
nsolvents=5;

kB=1.380662 10^-23;
nA=6.022045 10^23;
temp=298.15;

(*some values for the presentation of data*)
imagesize={600,600};
square=Graphics[{EdgeForm[Thick],White,Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
squareB=Graphics[{Black, Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]}];
solidsquare=Graphics[Polygon[{{-.01,-.01},{.01,-.01},{.01,.01},{-.01,.01}}]];
diamond=Graphics[{EdgeForm[Thick],White,Polygon[{{-.01,.0},{.0,.01},{.01,.0},{.0,-.01}}]}];
soliddiamond=Graphics[Polygon[{{-.01,.0},{.0,.01},{.01,.0},{.0,-.01}}]];
triangle=Graphics[{EdgeForm[Thick],White,Polygon[{{0,1},{Sin[2 Pi/3],Cos[2 Pi/3]},{Sin[4 Pi/3],Cos[4 Pi/3]}}]}];
circle=Graphics[{EdgeForm[{Thick,Black}],Thickness[0.1],Circle[{0,0},1]}];disk=Graphics[Disk[{0,0},1]];
plotmarkers={{diamond,.023},{square,.023},{triangle,.03},{solidsquare,.022},{soliddiamond,.03}};
plotmarkersD={squareB,0.02};
drange={{0,41},{0,41}};
drange2={{0,50},{0,50}};

(*obtains data from original data excel sheet*)
(*original Dexp*)
Dexp={Data[[1, 2, 4]]};
For[i=3, i≤110,i++, Dexp=AppendTo[Dexp, Data[[1, i,4]]]];

(*original visc data*)
visc={Data[[1, 2, 8]]};
For[i=3, i≤110,i++, visc=AppendTo[visc, Data[[1, i, 8]]]];
visc;

(*original MW data*)
MW={Data[[1, 2, 1]]};
For[i=3, i≤110,i++, MW=AppendTo[MW, Data[[1, i, 1]]]];
MW;

(*original MWs data*)
MWs={Data[[1, 2, 7]]};
For[i=3, i≤110,i++, MWs=AppendTo[MWs, Data[[1, i, 7]]]];
MWs;

(*original density data*)
dens={Data[[1, 2, 2]]};
For[i=3, i≤110,i++, dens=AppendTo[dens, Data[[1, i, 2]]]];
dens;

(*riginal packing data - whether sample is solid or liquid*)
pack={Data[[1, 2, 3]]};
For[i=3, i≤110,i++, pack=AppendTo[pack, Data[[1, i,3]]]];
pack;

data=Transpose[{MW , dens, pack,visc,MWs, Dexp 10^10 }];
ndp=Length[data]
adexptvsdse=Array[dexptvsdses,nsolvents];

(* Comparison of experimental and calculated Diffusion coefficients *)

(*calculate using SE method- data includes trimesic acid, old viscosities*)
(* dexptvsdses - array of experimental D against Stokes-Einstein D, calculated with simple M -> R relationship  *)

j=0;  prevmw=0;
For[i=1,i≤ndp,i++,rest=(3 data[[i,1]] data[[i,3]]10^-6/(data[[i,2]] 4 Pi na))^(1/3);If[data[[i,5]]!=prevmw,j++;dexptvsdses[j]={{10^10 kb t/(6 Pi data[[i,4]] 10^-3 rest),data[[i,6]]}};,AppendTo[dexptvsdses[j],{10^10 kb t/(6 Pi data[[i,4]] 10^-3 rest),data[[i,6]]}];];prevmw=data[[i,5]];];

(* Plot D comparisons *)
(* Experimental D vs Stokes-Einstein D *)

adexptvsdseplot=ListPlot[Tooltip[adexptvsdse],PlotRange->drange,AspectRatio->0.88,ImageSize->imagesize,PlotMarkers->plotmarkers,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.01],Frame->{{True,True},{True,True}},FrameStyle->{{Thick, Thick},{Thick, Thick}},FrameLabel->{Style["(Subscript[D, est]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False],Style["(Subscript[D, exp]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False]},LabelStyle->{24,Black, FontFamily->"Helvetica"}];
names=Graphics[{{Text[Style["(Subscript[CDCl, 3])",18,FontFamily->Helvetica],{32.25,13}],Text[Style["(Subscript[D, 2])O",18,SingleLetterItalics->False,FontFamily->Helvetica],{31.5,11}],Text[Style["[(Subscript[d, 6])]DMSO",18,SingleLetterItalics->False,FontFamily->Helvetica],{33.80,9}],Text[Style["[(Subscript[d, 8])]toluene",18,SingleLetterItalics->False,FontFamily->Helvetica],{33.95,7}],Text[Style["[(Subscript[d, 4])]methanol",18,SingleLetterItalics->False,FontFamily->Helvetica],{34.65,5}]}}];
legendmarkers=ListPlot[Table[{{28,17-(i+1)*2}},{i,1,5}],PlotRange->drange,AspectRatio->1,ImageSize->imagesize,PlotMarkers->plotmarkers,PlotStyle->{Thickness[0.003]},Axes->None];

slope45=Plot[Tooltip[x,"theoretical line"],{x,0,41},AspectRatio->1,ImageSize->imagesize,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.003],LabelStyle->18];
slope45v2=Plot[Tooltip[x,"theoretical line"],{x,0,50},AspectRatio->1,ImageSize->imagesize,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.003],LabelStyle->18];
Print[Style["Experimental D vs Stokes-Einstein D",18,FontFamily->Helvetica]];
Show[adexptvsdseplot,slope45,legendmarkers,names]

(*Export["RE_DvsSE.png",Show[adexptvsdseplot,slope45,legendmarkers,names]];
Export["RE_DvsSE.eps",Show[adexptvsdseplot,slope45,legendmarkers,names]];*)

expvspredall=Flatten[adexptvsdse,1];
dratios=Table[expvspredall[[i,2]]/expvspredall[[i,1]],{i,1,Length[expvspredall]}];
Print["RMS percentage deviation from Stokes-Einstein model  ",100 (Sum[(dratios[[i]]-1)^2,{i,1,Length[expvspredall]}]/Length[expvspredall])^(1/2)]; 

(*calculate using SEGWE method- data includes trimesic acid, old viscosities, original fitted density = 619.065 kg/m^3 *)

Clear[model]
α=(mws/mw)^(1/3);
model=10^10(kB temp((1.5α)+(1/(1+α))))/(6Pi η (0.003mw/(4 Pi ρeff nA))^(1/3));
temp=298.15;

ρeffopt=0.619065;

atestfitdata=Array[testfitdata,nsolvents];
j=0;  prevmw=0;
For[i=1,i≤ndp,i++,If[data[[i,5]]!=prevmw,j++;testfitdata[j]= {{model/.{mw->data[[i,1]]  10^-3,mws->data[[i,5]]  10^-3,η->  data[[i,4]]  10^-3,ρeff->ρeffopt},data[[i,6]]}};,AppendTo[testfitdata[j],{model/.{mw->data[[i,1]]  10^-3,mws->data[[i,5]]  10^-3,η->  data[[i,4]]  10^-3,ρeff->ρeffopt},data[[i,6]]}];];prevmw=data[[i,5]];];

DvsDSEGWE=ListPlot[Tooltip[atestfitdata], PlotRange->drange,AspectRatio->0.88,ImageSize->imagesize, PlotMarkers->plotmarkers,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.01],Frame->{{True,True},{True,True}},FrameStyle->{{Thick, Thick},{Thick, Thick}},FrameLabel->{Style["(Subscript[D, est]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False],Style["(Subscript[D, exp]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False]},LabelStyle->{24,Black, FontFamily->"Helvetica"}];
namesDvsD=Graphics[{{Text[Style["original data set",24,SingleLetterItalics->False,FontFamily->"Helvetica"],{32,6.13}]}}];
DvDplot=Show[DvsDSEGWE, legendmarkers,names,slope45]

(*Export["RE_DvsSE.png",Show[adexptvsdseplot,slope45,legendmarkers,names]];
Export["RE_DvsSE.eps",Show[adexptvsdseplot,slope45,legendmarkers,names]];*)


expvspredall=Flatten[atestfitdata,1];
dratios=Table[expvspredall[[i,2]]/expvspredall[[i,1]],{i,1,Length[expvspredall]}];
Print["RMS percentage deviation from Stokes-Einstein model  ",100 (Sum[(dratios[[i]]-1)^2,{i,1,Length[expvspredall]}]/Length[expvspredall])^(1/2)]; 

(*calculate using SEGWE method- data includes trimesic acid, arrhenius viscosities, original fitted density = 619.065 kg/m^3 *)

(*arrhenius viscosity data*)
Avisc={Data[[1, 2, 11]]};
For[i=3, i≤110,i++, Avisc=AppendTo[Avisc, Data[[1, i, 11]]]];
Avisc;

Clear[Data, DvsDSEGWE, atestfitdata,testfitdata];
data=Transpose[{MW , dens, pack,Avisc,MWs, Dexp 10^10 }];
ndp=Length[data]
adexptvsdse=Array[dexptvsdses,nsolvents];

Clear[model]
α=(mws/mw)^(1/3);
model=10^10(kB temp((1.5α)+(1/(1+α))))/(6Pi η (0.003mw/(4 Pi ρeff nA))^(1/3));
temp=298.15;

ρeffopt=0.619065;

atestfitdata=Array[testfitdata,nsolvents];
j=0;  prevmw=0;
For[i=1,i≤ndp,i++,If[data[[i,5]]!=prevmw,j++;testfitdata[j]= {{model/.{mw->data[[i,1]]  10^-3,mws->data[[i,5]]  10^-3,η->  data[[i,4]]  ,ρeff->ρeffopt},data[[i,6]]}};,AppendTo[testfitdata[j],{model/.{mw->data[[i,1]]  10^-3,mws->data[[i,5]]  10^-3,η->  data[[i,4]]  ,ρeff->ρeffopt},data[[i,6]]}];];prevmw=data[[i,5]];];


DvsDSEGWE=ListPlot[Tooltip[atestfitdata], PlotRange->drange,AspectRatio->0.88,ImageSize->imagesize, PlotMarkers->plotmarkers,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.01],Frame->{{True,True},{True,True}},FrameStyle->{{Thick, Thick},{Thick, Thick}},FrameLabel->{Style["(Subscript[D, est]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False],Style["(Subscript[D, exp]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False]},LabelStyle->{24,Black, FontFamily->"Helvetica"}];
namesDvsD=Graphics[{{Text[Style["original data set",24,SingleLetterItalics->False,FontFamily->"Helvetica"],{32,6.13}]}}];
DvDplot=Show[DvsDSEGWE, legendmarkers,names,slope45]

(*Export["RE_DvsSE.png",Show[adexptvsdseplot,slope45,legendmarkers,names]];
Export["RE_DvsSE.eps",Show[adexptvsdseplot,slope45,legendmarkers,names]];*)

expvspredall=Flatten[atestfitdata,1];
dratios=Table[expvspredall[[i,2]]/expvspredall[[i,1]],{i,1,Length[expvspredall]}];
Print["RMS percentage deviation from Stokes-Einstein model  ",100 (Sum[(dratios[[i]]-1)^2,{i,1,Length[expvspredall]}]/Length[expvspredall])^(1/2)]; 

(*image for TOC*)

SEplotTOC=ListPlot[Tooltip[adexptvsdse],PlotRange->drange2,AspectRatio->0.88,ImageSize->imagesize,PlotMarkers->plotmarkersD,PlotStyle->{Black,Thickness[0.003]},AxesStyle->Thickness[0.01],Frame->{{True,True},{True,True}},FrameStyle->{{Thick, Thick},{Thick, Thick}},FrameLabel->{Style["(Subscript[D, est]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False],Style["(Subscript[D, exp]) / (10^(-10)) (m^2) (s^(-2))   ⟶",SingleLetterItalics->False]},LabelStyle->{24,Black, FontFamily->"Helvetica"}];
namesDvsD=Graphics[{{Text[Style["original data set",24,SingleLetterItalics->False,FontFamily->"Helvetica"],{32,6.13}]}}];
SETOC=Show[SEplotTOC, slope45v2]

Export["TOC_DvsSE.tiff",SETOC];
Export["TOC_DvsSE.svg",SETOC];