in=[0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85];
r1 = [0, 4, 7, 11, 15, 18, 22, 25, 29, 32, 35, 37, 41, 43, 45, 46, 48, 48];
r2 = [0, 4, 7, 11, 15, 18, 22, 25, 30, 32, 35, 38, 41, 43, 45, 46, 48, 49];
r3 = [0, 4, 7, 11, 15, 18, 22, 25, 28, 32, 35, 38, 41, 43, 45, 46, 47, 48];
dr1 = [0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2];
dr2 = [0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2];
dr3 = [0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2];

average = 0.3333333*(sind(r1)+sind(r2)+sind(r3));

for i = 1:18
  temp = cosd(r1(i))*cosd(r1(i))*dr1(i)*dr1(i);
  temp += cosd(r2(i))*cosd(r2(i))*dr2(i)*dr2(i);
  temp += cosd(r3(i))*cosd(r3(i))*dr3(i)*dr3(i);
  error(i) = sqrt(temp) * 3.1415927/(3.* 180.);
end

sinin = sind(in);

p = polyfit(sinin, average,1);
plot(sinin,polyval(p,sinin), 'r-');
hold on;
errorbar (sinin, average, error)
h = legend ('numerical fit', 'Data', 'location', 'southeast');
set (h, "fontsize", 16);
xlabel('sin \theta_{in}', 'FontSize', 16)
ylabel('sin  \theta_r', 'FontSize', 16)
axis('tight')
set(gca,'FontSize',16)