Complex Fit Function Influences
Options & Clear
Clear["Global`*"]
SetOptions[$FrontEnd,DynamicEvaluationTimeout->20]
<<CompiledFunctionTools`
SetSystemOptions["CompileOptions"->"CompileReportExternal"->True]
On[Compile::noinfo]
dirname=""
dir=SetDirectory[dirname]
Factors
Constants
N[I]
N[π]
N[Pi]
N[E]
[c] = m/s
c=299792458.;
pi=3.14159265; (*N[Pi,9]*)
invpi= 0.318309886; (*1/N[Pi,9]*)
Prefix Transformation
centi
base2centi=10.^(2.);
centi2base=10.^(-2.);
milli
base2milli=10.^(3.);
milli2base=10.^(-3.);
micro
base2micro=10.^(6.);
micro2base=10.^(-6.);
nano
base2nano=10.^(9.);
nano2base=10.^(-9.);
pico
base2pico=10.^(12.);
pico2base=10.^(-12.);
femto
base2femto=10.^(15.);
femto2base=10.^(-15.);
Kilo
base2kilo=10.^(-3.);
kilo2base=10.^(3.);
Mega
base2mega=10.^(-6.);
mega2base=10.^(6.);
Giga
base2giga=10.^(-9.);
giga2base=10.^(9.);
Tera
base2tera=10.^(-12.);
tera2base=10.^(12.);
Conversion factors
Angular frequency transformations
Angular frequency to frequency ([ω]= rad/s  to [ν ] = 1/s = Hz)
ω = 2πν
ν = ω/(2π)
rads2Hz= 1./(2.*pi);
Angular frequency to wavenumbers ([ω]=rad/s to [Overscript[ν, ~]]= cm^-1)
ω=2πν= 2π c Overscript[ν, ~]
Overscript[ν, ~]=ω/(2π c)
rads2wn= (1./(2.*pi*c))*(1./base2centi);
Frequency transformations
Frequency to wavenumber ([ν]=1/s=Hz  to [Overscript[ν, ~]]= cm^-1)
ν=  c Overscript[ν, ~]
Overscript[ν, ~]=ν/c
Hz2wn=(1./c)*(1./base2centi);
Frequency to angular frequency (Hz to rad/s)
ω = 2πν
Hz2rads= 2.pi;
Wavenumber transformations
Wavenumber to angular frequency (cm^-1 to rad/s)
wn2rads=2.*pi*c * (1./centi2base)
Wavenumber to frequency (cm^-1 to 1/s=Hz)
wn2Hz= c* (1./centi2base);
Functions
If each function is transformed by itself and then added up the low frequency part gets too high. 
 α(ω)=(ω Im(Overscript[ε, ~](ω)))/(c n(ω))=(ω Im(Overscript[ε, ~](ω)))/(c  Re(Sqrt[Overscript[ε, ~](ω)]))
 Overscript[ε, ~](ω)=const+∑Subscript[f, i](ω)
α(ω)=(ω Im(Subscript[ε, ∞]+∑Subscript[f, i](ω)))/(c  Re(Sqrt[Subscript[ε, ∞]+∑Subscript[f, i](ω)]))=(ω Im(∑Subscript[f, i](ω)))/(c  Re(Sqrt[Subscript[ε, ∞]+∑Subscript[f, i](ω)]))=(ω/c*∑Im[Subscript[f, i](ω)])*(Re[Sqrt[Subscript[ε, ∞]+∑Subscript[f, i](ω)]])^-1
As such the fit to the frequency-dependent absorption coefficient cannot be simply made by a superposition of the fit function, as it is done for the response function in OKE. Instead, the superposition itself has to be transformed for the proper results.
The Imaginary part of the superposition can then be compared.
imaginary error function
The error function for complex arguments, erf(z), can be approximated (see Stegun or old notebooks for details) and holds well for arguments around -10 to 10. This approximation can be simplified erf(z)=erf(x+Iy) with x=0.
erf(0+I y)=(I y)/π+(2I)/π Σ (E^(-n^2/4)* sinh[n*y])/n
The imaginary error function, as used in the gaussian definition, is given by erfi(z)= -I erf(I z). Given that the argument passed to the imaginary error function in a complete antisymmetrized gaussian will always be real (cause we don't do imaginary frequencies), the argument passed onto the error function will always have a real part of 0 and the aforementioned approximation can be used.
As such:
erfi(Subscript[x, i]+I*0)= -I erf(I*(Subscript[x, i]+I*0))=-I erf(I Subscript[x, i]-0)=-I erf(I Subscript[x, i])= -I *((I Subscript[x, i])/π+(2I)/π Σ (E^(-n^2/4)* sinh[n*Subscript[x, i]])/n)=Subscript[x, i]/π+2/π Σ (E^(-n^2/4)* sinh[n*Subscript[x, i]])/n
erfiI0CD=Compile[{{xI0CD,_Real},{nI0CD,_Real},{inversepi,_Real}},((xI0CD*inversepi)+(2.*inversepi*Sum[(Exp[-0.25*(n^2)]*Sinh[n*xI0CD]/n),{n,1.,nI0CD,1.}]))];
complete Brownian oscillator
brownian[afBO_,o0fBO_,gammafBO_,omegafBO_]:=(afBO*(o0fBO^2.))/((o0fBO^2.)-omegafBO(omegafBO+(I 2.gammafBO)));
complete antisymmetric Gaussian
Cutoffs are inserted to avoid errors from numbers becoming too big/small. They are justified once the values involved become negligible, especially since an approximation is involved anyway and the data fitted has some assumptions and a not too high signal to noise.
The cutoff depends on the value becoming less than a thousandth and as such negligible, due to only minimal changes in the tail vs computational efficiency
For the exponential, if \!\(TraditionalForm\`x >= \ 2.63\), then the term goes to zero. Use the absolute because it is squared anyway
For the product of the exponential and the imaginary error function the cutoff is \!\(TraditionalForm\`x >= \ 11.71\), after which it can go to zero. Use the absolute, since the imaginary error function is uneven. The values will only switch sign, which doesn' t matter
(ω Im[f])/(c Re[Sqrt[f]])
We make the transformation after adding a real constant, so while the imaginary part can go to zero the real part shouldn' t just because of the cutoff.
In theory though Im[f] should go to zero quicker then Re[f] because the exponential is additionally multiplied with the error function, which has a positive exponent in the cases we are talking about.Also, we take the square root before we take the real part, so if we have a negative exponent then it' s size should be reduced again.Not to mention that we also multiply Re[f] by 10^8.
So to sum it up the denominator should end up bigger than the numerator in the cases where we set it to zero.
Even if we try to plot the way specific functions contribute we would still use the total for the real part, so no problems should arise there either.
gaussian=Compile[{{aASG,_Real},{o0ASG,_Real},{gammaASG,_Real},{omegaASG}},With[{x1=((omegaASG - o0ASG)/gammaASG),x2=((omegaASG + o0ASG)/gammaASG)},(Piecewise[{{I*aASG * Exp[-(x1^2.)],Abs[x1]<2.63},{0., Abs[x1]>= 2.63}}]- Piecewise[{{aASG * Exp[-(x1^2.)]* erfiI0CD[x1,20.,0.318309886],Abs[x1]<11.71},{0., Abs[x1]>= 11.71}}] - Piecewise[{{I*aASG * Exp[-(x2^2.)],Abs[x2]<2.63},{0., Abs[x2]>= 2.63}}]+Piecewise[{{aASG * Exp[-(x2^2.)]* erfiI0CD[x2,20.,0.318309886],Abs[x2]<11.71},{0., Abs[x2]>= 11.71}}])],CompilationOptions->{"InlineExternalDefinitions"->True},RuntimeAttributes->{Listable}];
Complete Debye

debye[aD_,tD_,omegaD_]:=aD/(1.-I omegaD tD);
Complete Havriliak - Negami
Subscript[H, α,β](ω)=1/(1+(-Iωτ)^α)^β 
hn[aHN_,tHN_,alphaHN_,betaHN_,omegaHN_]:=aHN/((1.+((-I omegaHN tHN)^(alphaHN)))^(betaHN));
The high frequency tail is unphysical because relaxational processes will follow each other as the time scale gets longer. As such a translational relaxation, for example, should not have a component that is at a higher frequency than a libration. This can be included in the frequency domain through a rise function as follows.

For a time domain exponential rise function

As such the HN becomes (for Im[HN])
Subscript[H, α,β](ω)=Subscript[A, HN]/(1+(-Iωτ)^α)^β-Subscript[A, HN]/(1+(-I(ω+ιγ)τ)^α)^β
	     =Subscript[A, HN]/(1+(-Iωτ)^α)^β-Subscript[A, HN]/(1+(-Iωτ+γτ)^α)^β
This is the correct form!!!!!!!!!!!!!!!!!!!!!!!
hnr[aHNR_,tHNR_,gHNR_,alphaHNR_,betaHNR_,omegaHNR_]:=(aHNR/((1.+((-I omegaHNR tHNR)^(alphaHNR)))^(betaHNR)))-(aHNR/((1.+((-I omegaHNR tHNR+gHNR tHNR)^(alphaHNR)))^(betaHNR)));
If α, the cole - cole broadening parameter, is very small the function, even with a rise correction, has a very long tail (due to the broadness of the cole - cole which approaches linearity if the scale is not big enough).
Combined fit function
If[dattype =="OKE",
fitfunc[x_,total_]:=Im[total],
If[dattype =="IR",
fitfunc[x_,total_]:=(((rads2Hz*(x*wn2rads))*Im[total])/(c*Re[Sqrt[total]]))
]
];
fittotal[con_,mat1_,mat2_,mat3_]:=Total[Join[{con},mat1,mat2,mat3]];
Parameters
eINF=1.5^2;
Gaussian/Brownian
aP= 1.;
o0P= 100.;
gP= 30.;
Havriliak - Negami
aHNP=1.;
tauHNP=2.;
gHNP=200.;
alphaHNP=1.;
betaHNP=1.;
Plotting
Gaussian/Brownian
boImReGraph=Plot[{Im[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]],Re[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]},{f,0,400},
PlotLegends->Placed[LineLegend[{"Im[B(ω)]","Re[B(ω)]"},LegendMarkerSize->20,LabelStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold}],{Right,Top}],
PlotPoints->500,Frame->True,
FrameLabel->{"Wavenumber [cm^-1]","ε'/ε''"},
RotateLabel->True,BaseStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold},
FrameStyle->Directive[Black,Thick],FrameTicks->None,
LabelStyle->Black,
PlotStyle->{{RGBColor[188/256,4/256,4/256],Thickness[.01]},{RGBColor[32/256,57/256,95/256],Thickness[.01]},{RGBColor[77/256,135/256,60/256],Thickness[.01]},{RGBColor[223/256,145/256,47/256],Thickness[.01]},{RGBColor[85/256,150/256,205/256],Thickness[.01]},{RGBColor[151/256,184/256,51/256],Thickness[.01]},{RGBColor[242/256,207/256,61/256],Thickness[.01]}},
ImageSize->{QuantityMagnitude[UnitConvert[Quantity[120/1000,"Meters"],"PrinterPoints" ]],QuantityMagnitude[UnitConvert[Quantity[90/1000,"Meters"],"PrinterPoints" ]]},Background->White,ImageMargins-> {{5,10},{0,0}},
PlotRangeClipping-> False,
Epilog->Inset[Style["(a)",14],Scaled[{-0.05,.98}]]]
Export[StringJoin[dirname,"Fig SI2a BO_Im_Re.png"],boImReGraph,ImageSize->1476,ImageResolution->500]
Plot[{Sqrt[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]],Im[Sqrt[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]],Re[Sqrt[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]]},{f,0,400},PlotStyle->{{Blue,Thick},{Dashed},{Dashed}}];
Plot[(((rads2Hz*(f*wn2rads))*Im[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]])/c),{f,0,400}];
tOKEtransB=Table[Im[eINF+brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]],{f,0,400}];
tIRtransB=Table[((2*(f*wn2rads)*Im[Sqrt[eINF+brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]])/c),{f,0,400}];
scaleBOKE=1/Max[tOKEtransB]
scaleBIR=1/Max[tIRtransB]
boIRShiftGraph=Plot[{(Im[eINF+brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]*scaleBOKE),(((2*(f*wn2rads)*Im[Sqrt[eINF+brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]])/c)*scaleBIR)},{f,0,400},
PlotLegends->Placed[LineLegend[{"Im[Subscript[ε, ∞]+B(ω)]","(ω*Im[B(ω)])/(c*Re Sqrt[Subscript[ε, ∞]+B(ω)])"},LegendMarkerSize->20,LabelStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold}],{Right,Top}],
PlotPoints->500,Frame->True,
FrameLabel->{"Wavenumber [cm^-1]","ε'/ε''"},
RotateLabel->True,BaseStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold},
FrameStyle->Directive[Black,Thick],FrameTicks->None,
LabelStyle->Black,
PlotStyle->{{RGBColor[188/256,4/256,4/256],Thickness[.01]},{RGBColor[32/256,57/256,95/256],Thickness[.01]},{RGBColor[77/256,135/256,60/256],Thickness[.01]},{RGBColor[223/256,145/256,47/256],Thickness[.01]},{RGBColor[85/256,150/256,205/256],Thickness[.01]},{RGBColor[151/256,184/256,51/256],Thickness[.01]},{RGBColor[242/256,207/256,61/256],Thickness[.01]}},
ImageSize->{QuantityMagnitude[UnitConvert[Quantity[120/1000,"Meters"],"PrinterPoints" ]],QuantityMagnitude[UnitConvert[Quantity[90/1000,"Meters"],"PrinterPoints" ]]},Background->White,ImageMargins-> {{5,10},{0,0}},PlotRange->Full,
PlotRangeClipping-> False,
Epilog->Inset[Style["(b)",14],Scaled[{-0.05,.98}]]]
Export[StringJoin[dirname,"Fig SI2b BO_IRShift.png"],boIRShiftGraph,ImageSize->1476,ImageResolution->500]
Table[{f,((((rads2Hz*(f*wn2rads))*Im[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]])/(c*Re[Sqrt[brownian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]]))*scaleB)},{f,0,400}];
asgImReGraph=Plot[{Im[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]],Re[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]},{f,0,400},
PlotLegends->Placed[LineLegend[{"Im[G(ω)]","Re[G(ω)]"},LegendMarkerSize->20,LabelStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold}],{Right,Top}],
PlotPoints->500,Frame->True,
FrameLabel->{"Wavenumber [cm^-1]","ε'/ε''"},
RotateLabel->True,BaseStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold},
FrameStyle->Directive[Black,Thick],FrameTicks->None,
LabelStyle->Black,
PlotStyle->{{RGBColor[188/256,4/256,4/256],Thickness[.01]},{RGBColor[32/256,57/256,95/256],Thickness[.01]},{RGBColor[77/256,135/256,60/256],Thickness[.01]},{RGBColor[223/256,145/256,47/256],Thickness[.01]},{RGBColor[85/256,150/256,205/256],Thickness[.01]},{RGBColor[151/256,184/256,51/256],Thickness[.01]},{RGBColor[242/256,207/256,61/256],Thickness[.01]}},
ImageSize->{QuantityMagnitude[UnitConvert[Quantity[120/1000,"Meters"],"PrinterPoints" ]],QuantityMagnitude[UnitConvert[Quantity[90/1000,"Meters"],"PrinterPoints" ]]},Background->White,ImageMargins-> {{5,10},{0,0}},PlotRange->Full,
PlotRangeClipping-> False,
Epilog->Inset[Style["(a)",14],Scaled[{-0.05,.98}]]]
Export[StringJoin[dirname,"Fig SI1a ASG_Im_Re.png"],asgImReGraph,ImageSize->1476,ImageResolution->500]
Plot[{Sqrt[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]],Im[Sqrt[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]],Re[Sqrt[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]]},{f,0,400},PlotStyle->{{Blue,Thick},{Dashed},{Dashed}}];
Plot[(((rads2Hz*(f*wn2rads))*Im[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]])/c),{f,0,400}, PlotRange->All];
scaleG=1./Max[Table[((2*(f*wn2rads)*Im[Sqrt[eINF+gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]])/c),{f,0,400}]]
asgIRShiftGraph=Plot[{Im[eINF+gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]],(((2*(f*wn2rads)*Im[Sqrt[eINF+gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]])/c)*scaleG)},{f,0,400},
PlotLegends->Placed[LineLegend[{"Im[Subscript[ε, ∞]+G(ω)]","(2ω Im[Sqrt[Subscript[ε, ∞]+G(ω)]])/c"},LegendMarkerSize->20,LabelStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold}],{Right,Top}],
PlotPoints->500,Frame->True,
FrameLabel->{"Wavenumber [cm^-1]","ε'/ε''"},
RotateLabel->True,BaseStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold},
FrameStyle->Directive[Black,Thick],FrameTicks->None,
LabelStyle->Black,
PlotStyle->{{RGBColor[188/256,4/256,4/256],Thickness[.01]},{RGBColor[32/256,57/256,95/256],Thickness[.01]},{RGBColor[77/256,135/256,60/256],Thickness[.01]},{RGBColor[223/256,145/256,47/256],Thickness[.01]},{RGBColor[85/256,150/256,205/256],Thickness[.01]},{RGBColor[151/256,184/256,51/256],Thickness[.01]},{RGBColor[242/256,207/256,61/256],Thickness[.01]}},
ImageSize->{QuantityMagnitude[UnitConvert[Quantity[120/1000,"Meters"],"PrinterPoints" ]],QuantityMagnitude[UnitConvert[Quantity[90/1000,"Meters"],"PrinterPoints" ]]},Background->White,ImageMargins-> {{5,10},{0,0}},PlotRange->Full,
PlotRangeClipping-> False,
Epilog->Inset[Style["(b)",14],Scaled[{-0.05,.98}]]]
Export[StringJoin[dirname,"Fig SI1b ASG_IRShift.png"],asgIRShiftGraph,ImageSize->1476,ImageResolution->500]
Table[{f,((((rads2Hz*(f*wn2rads))*Im[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]])/(c*Re[Sqrt[gaussian[aP,(o0P*wn2rads),(gP*wn2rads),(f*wn2rads)]]]))*scaleG)},{f,0,400}];
HN
hnrLINGraph=Plot[{Im[hnr[aHNP,(tauHNP*(pico2base/Hz2rads)),(gHNP*wn2rads),1.,1.,(f*wn2rads)]],Im[hnr[aHNP,(tauHNP*(pico2base/Hz2rads)),(gHNP*wn2rads),.5,1.,(f*wn2rads)]],Im[hnr[aHNP,(tauHNP*(pico2base/Hz2rads)),(gHNP*wn2rads),1.,.5,(f*wn2rads)]]},{f,0,400},
PlotLegends->Placed[LineLegend[{"Subscript[H, 1,1](ω)","Subscript[H, 0.5,1](ω)","Subscript[H, 1,0.5](ω)"},LegendMarkerSize->20,LabelStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold}],{Right,Top}],
PlotPoints->500,Frame->True,
FrameLabel->{"Wavenumber [cm^-1]","ε''"},
RotateLabel->True,BaseStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold},
FrameStyle->Directive[Black,Thick],FrameTicks->None,
LabelStyle->Black,
PlotStyle->{{RGBColor[188/256,4/256,4/256],Thickness[.01]},{RGBColor[32/256,57/256,95/256],Thickness[.01]},{RGBColor[77/256,135/256,60/256],Thickness[.01]},{RGBColor[223/256,145/256,47/256],Thickness[.01]},{RGBColor[85/256,150/256,205/256],Thickness[.01]},{RGBColor[151/256,184/256,51/256],Thickness[.01]},{RGBColor[242/256,207/256,61/256],Thickness[.01]}},
ImageSize->{QuantityMagnitude[UnitConvert[Quantity[120/1000,"Meters"],"PrinterPoints" ]],QuantityMagnitude[UnitConvert[Quantity[90/1000,"Meters"],"PrinterPoints" ]]},Background->White,ImageMargins-> {{5,10},{0,0}},PlotRange->Full,
PlotRangeClipping-> False,
Epilog->Inset[Style["(a)",14],Scaled[{-0.05,.98}]]]
Export[StringJoin[dirname,"Fig SI3a HNR_LinPlot.png"],hnrLINGraph,ImageSize->1476,ImageResolution->500]
hnrLOGGraph=LogLinearPlot[{Im[hnr[aHNP,(tauHNP*(pico2base/Hz2rads)),(gHNP*wn2rads),1.,1.,(f*wn2rads)]],Im[hnr[aHNP,(tauHNP*(pico2base/Hz2rads)),(gHNP*wn2rads),.5,1.,(f*wn2rads)]],Im[hnr[aHNP,(tauHNP*(pico2base/Hz2rads)),(gHNP*wn2rads),1.,.5,(f*wn2rads)]]},{f,1,400},
PlotLegends->Placed[LineLegend[{"Subscript[H, 1,1](ω)","Subscript[H, 0.5,1](ω)","Subscript[H, 1,0.5](ω)"},LegendMarkerSize->20,LabelStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold}],{Right,Top}],
PlotPoints->500,Frame->True,
FrameLabel->{"Wavenumber [cm^-1]","ε''"},
RotateLabel->True,BaseStyle->{FontFamily->"Myriad Pro",FontSize->14,FontWeight->Bold},
FrameStyle->Directive[Black,Thick],FrameTicks->None,
LabelStyle->Black,
PlotStyle->{{RGBColor[188/256,4/256,4/256],Thickness[.01]},{RGBColor[32/256,57/256,95/256],Thickness[.01]},{RGBColor[77/256,135/256,60/256],Thickness[.01]},{RGBColor[223/256,145/256,47/256],Thickness[.01]},{RGBColor[85/256,150/256,205/256],Thickness[.01]},{RGBColor[151/256,184/256,51/256],Thickness[.01]},{RGBColor[242/256,207/256,61/256],Thickness[.01]}},
ImageSize->{QuantityMagnitude[UnitConvert[Quantity[120/1000,"Meters"],"PrinterPoints" ]],QuantityMagnitude[UnitConvert[Quantity[90/1000,"Meters"],"PrinterPoints" ]]},Background->White,ImageMargins-> {{5,10},{0,0}},PlotRange->Full,
PlotRangeClipping-> False,
Epilog->Inset[Style["(b)",14],Scaled[{-0.05,.98}]]]
Export[StringJoin[dirname,"Fig SI3b HNR_LogPlot.png"],hnrLOGGraph,ImageSize->1476,ImageResolution->500]