/* Data from a photo of the splatter pattern made by a falling bottle of nail polish. A 1 cm grid of lines was drawn over the photo and the the values represent coordinates (in cm) from an arbitrary origin. Data and analysis from Wicklin (2015), "The spiral of splatter", The DO Loop, http://blogs.sas.com/content/iml/2015/06/11/spiral-of-splatter.html published June 11, 2015. */ data Spiral; input Segment x y; if _N_ < 28 | _N_ > 58 then Surface = "Wall"; else Surface = "Trim"; /* reverse order to parameterize by polar angle, theta */ Seg = 85 - _N_ + 1; datalines; 1 0 0 2 0.25 0.5 3 0.8 0.6 4 1.2 0.6 5 1.7 0.5 6 2.2 0.2 7 2.5 0 8 2.8 -0.6 9 3 -1 10 3.1 -1.5 11 3.2 -1.9 12 3.1 -2.5 13 3 -2.9 14 2.8 -3.3 15 2.6 -3.8 16 2.3 -4.2 17 1.8 -4.6 18 1.4 -4.8 19 0.9 -5.2 20 0 -5.4 21 -0.5 -5.5 22 -1 -5.6 23 -1 -5.6 24 -1.5 -5.6 25 -2 -5.5 26 -2.5 -5.3 27 -3 -5.2 28 -3.5 -5.1 29 -3.9 -4.9 30 -4.3 -4.6 31 -4.7 -4.3 32 -5.1 -3.9 33 -5.5 -3.5 34 -5.7 -3.1 35 -6 -2.7 36 -6.3 -2.2 37 -6.5 -1.7 38 -6.7 -1.2 39 -6.8 -0.7 40 -6.9 -0.2 41 -7 0.4 42 -7 0.9 43 -7 1.4 44 -7 1.9 45 -7 2.4 46 -6.9 2.8 47 -6.8 3.4 48 -6.7 3.9 49 -6.5 4.5 50 -6.4 5 51 -6.2 5.3 52 -6 5.7 53 -5.7 6.1 54 -5.5 6.5 55 . . 56 -4.7 7.2 57 -4.3 7.7 58 -3.9 8 59 -3.8 9.7 60 -3.3 9.9 61 -2.9 10.2 62 -2.5 10.4 63 -2.1 10.7 64 -1.7 10.9 65 -1.2 11.1 66 -0.7 11.3 67 -0.2 11.5 68 0.2 11.7 69 0.7 11.8 70 1.2 11.9 71 1.8 12.1 72 2.2 12.1 73 2.8 12.2 74 3.2 12.3 74 . . 75 . . 76 . . 77 . . 78 5.7 12.5 79 . . 80 6.6 12.5 81 . . 82 7.8 12.5 83 . . 84 . . 85 9.4 12.4 ; proc sort data=Spiral; by seg; run; title "Splatter Pattern"; title2 "Nail Polish on Wall"; ods graphics / width=400px height=400px; proc sgplot data=Spiral; scatter x=x y=y / group=Surface nomissinggroup; xaxis grid values=(-8 to 12 by 2); yaxis grid values=(-6 to 14 by 2); keylegend / location=inside position=bottomright opaque across=1; run; title; proc iml; pi = constant("pi"); theta = do(0,6*pi, pi/24)`; /* simulate data */ a = 1; b = -0.3063; x = a*cos(theta)#exp(b*theta); y = a*sin(theta)#exp(b*theta); call scatter(x,y); /* Define module that returns arctan on [-pi, infinity), assuming that a curve is parametrized by the polar angle. See "Computing polar angles from coordinate data" http://blogs.sas.com/content/iml/2015/06/10/polar-angle-curve.html ?*/ start IncrAngle(x, y); theta = atan2(y,x); /* angle from center to pt */ dif = dif(theta); idx = loc(dif^=. & dif<0); n = nrow(dif); pi = constant("pi"); do while( ncol(idx)>0 ); s = idx[1]; theta[s:n] = 2*pi + theta[s:n]; dif = dif(theta); idx = loc(dif^=. & dif<0); end; return(theta); finish; /* test it */ phi = IncrAngle(x,y); /* angle from center to pt */ print (max(abs(theta-phi))); /***** MAIN COMPUTATION *****/ /* model logarithmic spiral */ use Spiral; read all var {x y}; close; /* read observed data */ start ObjF(p) global(x, y); x0=p[1]; y0=p[2]; a=p[3]; b=p[4]; /* parameters to optimize */ xc = x-x0; yc = y-y0; /* center the data */ theta = IncrAngle(xc, yc); /* angle from center to pt */ u = a*cos(theta)#exp(b*theta); /* predicted values */ v = a*sin(theta)#exp(b*theta); d2 = (xc-u)##2 + (yc-v)##2; /* sqr radial distance to prediction */ return( sum(d2) ); finish; /* guess */ p = {1 -1 26 -0.5}; /* initial guess for (x0, y0) */ opt = {0 1}; /* options=(min, print level) */ call nlptr(rc,result,"ObjF", p, opt); /* find optimum (x0, y0) */ print rc,result[c={"x0" "y0" "a" "b"} label="Parameter Estimates"]; /* save parametric curve to data set "Fit" */ x0=result[1];y0= result[2]; a= result[3];b= result[4]; xc = x-x0; yc = y-y0; /* center the data */ theta = IncrAngle(xc, yc); /* angle from center to pt */ u = x0 + a*cos(theta)#exp(b*theta); /* predicted values */ v = y0 + a*sin(theta)#exp(b*theta); create Fit var {"theta" "u" "v"}; append; close; quit; /* merge observed data and fit */ data All; set Spiral Fit; run; title "Fit Logarithmic Spiral to Splatter Pattern"; proc sgplot data=All; scatter x=x y=y / group=Surface nomissinggroup; series x=u y=v / legendlabel="Fit"; xaxis grid values=(-8 to 12 by 2); yaxis grid values=(-6 to 14 by 2); keylegend / location=inside position=bottomright opaque across=1; run; /* Attempt physics-based model. Assume that a bottle of half-height r a) Is falling straight down from height h0 with no initial velocity b) Is rotating with angular velocity w about its center so that the linear velocity of the mouth/nozzle is r*w. The height of the bottle is given by h(t) = h0 - gt^2 / 2 [Eqn 1] where g = 9.8 m/s^2 is the acceleration due to gravity. The bottle hits the ground when t = sqrt(2*h0/g). Assume that a drop that emerges from the the nozzle falls at the same rate as the bottle. If the drop emerges at time t, it travels in the tangential direction for Td = sqrt(2*h0/g)-t [Eqn 2] seconds before it hits the ground. (? Not sure about Eqn 2. How long does it take a drop to hit the ground? Do I need to account for the downward speed of the bottle?) In that time it has traveled d(t) = r*w*Td meters [Eqn 3]. The (x,y) coordinates of the mouth of the bottle is located at (r*cos(wt), r*sin(wt)) [Eqn 4] so the unit vector in the tangent direction is (-sin(wt), cos(wt)) [Eqn 5] */ proc iml; pi = constant("pi"); r = 0.04; /* bottle about 8 cm high, so r=0.04 m */ w = 3*2*pi; /* assume 3 rotations per second */ g = 9.8; /* accel due to gravity */ h0 = 1; /* initial height is 1 m */ t = do(0, sqrt(2*h0/g), 0.01); /* bottle falls in this time */ nozzle = r*cos(w*t) // r*sin(w*t); /* position of nozzle */ tangent = -sin(w*t) // cos(w*t); /* unit vector in tangent direction */ Td = sqrt(2*h0/g) - t; /* ? drop in air this long ??? */ x = nozzle + r*w*Td#tangent; ods graphics / width=320px height=400px; title "Path of Liquid from Rotating Bottle"; title2; call series(x[1,], x[2,]) grid={x y} yvalues=do(-0.2, 0.4, 0.1) other="refline 0 / axis=x; refline 0 / axis=y;";