Polar coordinate data adjustment

Addition: 95% confidence bands for the aggregate average.

(* Convert from degrees to radians *)
datar = Import("https://pastebin.com/raw/CM7Rj6jC", "Table");
datar((All, 1)) = 2 π datar((All, 1))/360;

(* Fit curves *)
nlm1 = NonlinearModelFit(datar((All, {1, 2})), a Cos(t + ϕ)^2, {a, ϕ}, t);
nlm2 = NonlinearModelFit(datar((All, {1, 3})), a Cos(t + ϕ)^2, {a, ϕ}, t);

(* Plot results *)
mpb1 = Table(Flatten({t, nlm1("MeanPredictionBands")}), {t, 0, 2 π, π/50});
mpb2 = Table(Flatten({t, nlm2("MeanPredictionBands")}), {t, 0, 2 π, π/50});
Show(ListPolarPlot({datar((All, {1, 2})), datar((All, {1, 3}))}, PlotStyle -> PointSize(0.02)),
  PolarPlot({nlm1(t), nlm2(t)}, {t, 0, 2 π}),
  ListPolarPlot({mpb1((All, {1, 2})), mpb1((All, {1, 3})), 
    mpb2((All, {1, 2})), mpb2((All, {1, 3}))},
    PlotStyle -> {{Blue, Dotted}, {Blue, Dotted}, {Orange, Dotted}, {Orange, Dotted}}, Joined -> True))

Data and adjustment with 95% confidence bands for the average