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 [(OverTilde[ν])]= (cm^(-1)))
ω=2πν= 2π c (OverTilde[ν])
(OverTilde[ν])=(ω/(2π c))

rads2wn= (1./(2.*pi*c))*(1./base2centi);

Frequency transformations

Frequency to wavenumber ([ν]=1/s=Hz  to [(OverTilde[ν])]= (cm^(-1)))
ν=  c (OverTilde[ν])
(OverTilde[ν])=(ν/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((OverTilde[ε])(ω)))/(c n(ω)))=((ω Im((OverTilde[ε])(ω)))/(c  Re((Sqrt[(OverTilde[ε])(ω)]))))
 (OverTilde[ε])(ω)=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(ω)]","\!\(\*FractionBox[\(2  ω\\ Im[\*SqrtBox[\(\*SubscriptBox[\(ε\), \(∞\)] + 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]