BrianWeinstein

is_outlier <- function(x) {
    # Computes a bootstrapped confidence interval
    # INPUT: 
    #    x: a numeric vector
    # OUTPUT:
    #    a boolean vector, indicating if each value in x is an outlier
    
    return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

BrianWeinstein

bootstrap_ci <- function(data, sample_size_pct = 0.50, samples = 100, conf_level = 0.95){
    # Computes a bootstrapped confidence interval
    # INPUT: 
    #    data: a numeric vector
    #    sample_size_pct: the percentage of the input data to be used in each bootsrapped sample
    #    samples: the number of samples
    #    conf_level: the desired confidence level
    # OUTPUT:
    #    a bootstrapped conf_level confidence interval
    

BrianWeinstein

# tm::removeSparseTerms attempts to remove sparse terms via slicing a sparse matrix.
# The slicing operation tries to convert the sparse matrix to a dense matrix, but this
# fails if the dense matrix has more than ((2^31) - 1) entries [i.e., if (nrow * ncol) > ((2^31) - 1)]
# 
# The error message is 
# In nr * nc : NAs produced by integer overflow
# 
# Instead of using tm::removeSparseTerms, the following function subsets the sparse matrix directly
# and avoids converting the sparse matrix to a dense one.

BrianWeinstein

ic[x_, y_] := 1 E^(-350 ((x - 1/5)^2 + ( y - 1/3)^2))

solnDir =
 NDSolve[
  {D[u[x, y, t], {t, 2}] == D[u[x, y, t], {x, 2}] + D[u[x, y, t], {y, 2}],
   u[x, y, 0] == ic[x, y],
   (D[u[x, y, t], t] /. t -> 0) == 0,
   u[0, y, t] == ic[0, y],
   u[1, y, t] == ic[1, y],
   u[x, 0, t] == ic[x, 0],

BrianWeinstein

circ = 1; rad = circ/(2 \[Pi]); nRunners = 5;
rList[t_] := {1 t, 2 t, 4 t, 8 t, 9.6 t, 21 t, 31 t, 33 t}[[1 ;; nRunners]]

dist[d\[Theta]_, circ_] := 
  N[circ/2 (TriangleWave[(d\[Theta] - \[Pi]/2)/(2 \[Pi])] + 1)/2]
minDist[runnerList_, circ_] :=
 Table[
  runner = runnerList[[i]];
  other = DeleteCases[runnerList, runner];
  Min[dist[Abs[runner - other], circ]],

BrianWeinstein

x[A1_, A2_, f1_, f2_, p1_, p2_, d1_, d2_, t_] := A1 Sin[t f1 + p1] E^(-d1 t) + A2 Sin[t f2 + p2] E^(-d2 t)
y[A3_, A4_, f3_, f4_, p3_, p4_, d3_, d4_, t_] := A3 Sin[t f3 + p3] E^(-d3 t) + A4 Sin[t f4 + p4] E^(-d4 t)
 
 
Manipulate[
 ParametricPlot[
  {x[A1, A2, f1, f2, p1, p2, d1, d2, t],
   y[A3, A4, f3, f4, p3, p4, d3, d4, t]},
  {t, 0, tmax},
  PlotPoints -> 200, Axes -> False, PlotStyle -> {Thick, Opacity[0.5]}, PlotRange -> All

BrianWeinstein

x[n_, \[Theta]_] := 2*Cos[Pi/(2*n)]*Cos[(1/2)*(\[Theta] + (Pi/n)*(2*Floor[(n*\[Theta])/(2*Pi)] + 1))] - Cos[(Pi/n)*(2*Floor[(n*\[Theta])/(2*Pi)] + 1)]

y[n_, \[Theta]_] := 2*Cos[Pi/(2*n)]*Sin[(1/2)*(\[Theta] + (Pi/n)*(2*Floor[(n*\[Theta])/(2*Pi)] + 1))] - Sin[(Pi/n)*(2*Floor[(n*\[Theta])/(2*Pi)] + 1)]

reuRotate[n_, \[Phi]_] := 
  {pts[n, \[Phi]] = Table[RotationMatrix[\[Phi]] . {x[n, \[Theta]], y[n, \[Theta]]}, {\[Theta], 0, 2*Pi, (2*Pi)/100}]; 
   xmin = Min[pts[n, \[Phi]][[All,1]]]; 
   ymin = Min[pts[n, \[Phi]][[All,2]]];
   xmax = Max[pts[n, \[Phi]][[All,1]]]; 
   ymax = Max[pts[n, \[Phi]][[All,2]]]; 

BrianWeinstein

(* Rutherford model *)

es[t_] := {{Cos[t + Pi], 2.5*Sin[t + Pi], 2*Sin[t + Pi]}, 
           {Cos[t + (4*Pi)/5], 2*Sin[t + (4*Pi)/5], -1.5*Sin[t + (4*Pi)/5]}, 
           {2*Sin[t + (3*Pi)/5], Cos[t + (3*Pi)/5], 2*Sin[t + (3*Pi)/5]}, 
           {2.5*Sin[t + (2*Pi)/5], Cos[t + (2*Pi)/5], -1.5*Sin[t + (2*Pi)/5]}}
Manipulate[Show[
    ParametricPlot3D[Evaluate[es[2*u]], {u, 0, 2*Pi}, PlotStyle -> Directive[Thick, Dotted]], 
    Graphics3D[
           {Specularity[White, 200], 

BrianWeinstein

nc = 15; nr = 3; 
cx = Table[ToExpression[StringJoin["cx", ToString[i]]], {i, 1, nc}]; 
cy = Table[ToExpression[StringJoin["cy", ToString[i]]], {i, 1, nc}]; 
rx = Table[ToExpression[StringJoin["rx", ToString[i]]], {i, 1, nr}]; 
ry = Table[ToExpression[StringJoin["ry", ToString[i]]], {i, 1, nr}]; 
coordList = Flatten[{Transpose[{cx, cy}], Transpose[{rx, ry}]}]; 
cspeed = 1; 
rspeed = 1.1; 

eqns = Flatten[

BrianWeinstein

xmin = -3.6; xmax = 3.6; 

p[x_] := 6 - Sqrt[6^2 - x^2] (* circle *)
p[x_] := x^2/7.5 (* parabola *)

img[t_, rays_] := 
  Show[
    Graphics[
      {Thick, RGBColor[0.243, 0.62, 0.612], 
       Table[Line[{{xi, -t}, {xi, 20}}], {xi, xmin + 0.25, xmax - 0.25, (xmax - xmin - 0.5)/rays}],