Needs["ErrorBarPlots`"]
<<PlotLegends`;

(* OK, here I am pretending to have taken data to give me something to plot.  I invent a list of times for the x-axis ... *) 
itimes = {10,11,12.55,16.00, 17.00,24+10.15,24+14,24+15.3,24+15.3,24+15.3,24+16,48+10.8,48+10.8,48+10.8,48+12,48+13.01,48+14.5,48+15,48+16};
times = Join[itimes,itimes,itimes];
xtimes = Table[i,{i,0,70}];
Length[times]

(* Here's the rest of the fake data.   I use the actual solar-elevation equation to get elevations ... *)
(* and I use the RandomReal[] function to add fake noise to the y-axes.*)
psi1 = 23.5*π/180;
psi[t_] := psi1*Sin[2*π*(t/24)/365.25];
phi = 35*π/180;
theta[t_] := 2*π*t/24;
eta[t_] := ArcCos[Cos[psi[t]]*Cos[phi]*Cos[theta[t]] + Sin[psi[t]*Sin[phi]]] + RandomReal[NormalDistribution[0,2*π/180]]

(* OK: given times and angles and knowledge of the error bars, we'll organize the data in a way suitable for plotting. *) 
data = Table[{{times[[i]],eta[times[[i]]]},ErrorBar[10/60,2*π/180]},{i,1,Length[times]}];
datanox =  Table[{times[[i]],eta[times[[i]]]},{i,1,Length[times]}];

(* Fit the data to a sine wave, how does that look?   Fiddle with the PlotRange to make it look better.*)
nlf = NonlinearModelFit[datanox,a+b*Sin[2*π*w*(t + c)],{{w,1/24.},{c,-5},{a,1.56},{b,0.9147}},t]
nlf["ParameterTable"]
tmp = Show[
Plot[nlf[t],{t,4,80},AxesLabel->{"time, h","elevation, rad"},LabelStyle->Larger,PlotRange->{1.2,2.6}],
ErrorListPlot[data]]

(* let's grab the best-fit values and error bars. *)
omega = w/. nlf["BestFitParameters"]
omegaerror =  nlf["ParameterErrors"][[1]]
(* and convert them to times .. *)
daylength = 1/omega 
daylengtherror = daylength * (omegaerror/omega)

(* Looks good!  With that, I'll start the error analysis.   What do the "residuals" look like? *)

resid = Table[{{datanox[[i,1]],datanox[[i,2]] - nlf[datanox[[i,1]]]},ErrorBar[10/60,2*π/180]},{i,1,Length[times]}];
ErrorListPlot[resid]

(* the residuals plot is not bad looking, but I'm interested in seeing histogram of the amplitudes of those residuals. *) 
Histogram[resid[[All,1,2]],{-0.12,0.12,0.02}]


(* Unfortunately the "Histogram[]" command is not as flexible as I want.   I really want to add the (random-counting-error) error bars to the heights of those histogram bins.   To do that, I need to obtain the histogram as a list of numbers, not as a plot.   Use "BinCounts":*)
bc = BinCounts[resid[[All,1,2]],{-0.12,0.12,0.02}];
(* the error on each bin's value N is approximately Sqrt[N], so let's build a new ErrorBarPlot dataset which
draws the "bincounts" data, but with Sqrt[bincounts] error bars. *)
residdata = Table[{{-0.12 + 0.01+0.02*i,bc[[i]]},ErrorBar[Max[1,Sqrt[bc[[i]]]]]},{i,1,12}];
(* Clean up this plot and overlay a Gaussian distribution on it.  That goes in the paper ... *) 
myinset = Show[ErrorListPlot[residdata,PlotLabel->"angle (theory-data), rad",PlotStyle->Red],Plot[1.2*PDF[NormalDistribution[0.0,2*π/180],x],{x,-0.12,0.12},PlotStyle->Red]]

(* OK, now that I have a nice data plot, and a nice error-checking plot, combine then using Show, Epilog, and Inset.  This took >1 hr of help-file-reading to get working. *)
tmp = Show[
Plot[nlf[t],{t,4,80},AxesLabel->{"time, h","elevation, rad"},ImageSize->Large,LabelStyle->Larger,PlotRange->{1.2,2.6},Epilog->Inset[Graphics[myinset],Scaled[{.85,.30}],Center,25]],
ErrorListPlot[data]]
Export["example_fig.pdf",tmp]