| >  | 
						
						 with(LinearAlgebra): 
						with(plots): 
						with(Statistics): 
						 | 
					 
				
			 
			   
			The principle is always the same: 
			    1/   Let L a straight line which is either defined by its two ending points (xkvd_hline) or taken as the default [0, 0], [1, 0] line. 
			          For xkvd_hline the given line L is firstly rotate to be aligned with the horizontal axis. 
			 
			    2/   Let P1, ..., PN N points on L. Each Pn writes [xn, yn] 
			 
			    3/   A random perturbation rn is added yo the values y1, ..., yN 
			 
			    4/   A stationnary random process RP, with gaussian correlation function is used to build a smooth curve passing through the points 
			          (x1, y1+r1), ..., (xN, yN+rN) (procedure KG where "KG" stands for "Kriging") 
			 
			    5/   The result is drawn or mapped to some predefined shape : 
			                  xkcd_hist, 
			                  xkcd_polyline, 
			                  xkcd_circle 
			 
			    6/   A procedure xkcd_func is also provided to draw functions defined by an explicit relation. 
			  
			
				
					
						| >  | 
						
						 KG := proc(X, Y, psi, sigma) 
						  local NX, DX, K, mu, k, y: 
						  NX := numelems(X); 
						  DX := < seq(Vector[row](NX, X), k=1..NX) >^+: 
						  K  := (sigma, psi) -> evalf( sigma^2 *~ exp~(-((DX - DX^+) /~ psi)^~2) ): 
						  mu := add(Y) / NX; 
						  k  := (x, sigma, psi) -> evalf( convert(sigma^2 *~ exp~(-((x -~ X ) /~ psi)^~2), Vector[row]) ): 
						  y  := mu + k(x, sigma, psi) . (K(sigma, psi))^(-1) . convert((Y -~ mu), Vector[column]): 
						  return y 
						end proc: 
						 
						 
						xkcd_hline := proc(p1::list, p2::list, a::nonnegative, lc::positive, col) 
						  # p1 : first ending point 
						  # p2 : second ending point 
						  # a  : amplitude of the random perturbations 
						  # lc : correlation length 
						  # col: color 
						  local roll, NX, LX, X, Z: 
						  roll := rand(-1.0 .. 1.0): 
						  NX   := 10: 
						  LX   := p2[1]-p1[1]: 
						  X    := [seq(p1[1]..p2[1], LX/(NX-1))]: 
						  Z    := [p1[2], seq(p1[2]+a*roll(), k=1..NX-1)]: 
						  return plot(KG(X, Z, lc*LX, 1), x=min(X)..max(X), color=col, scaling=constrained): 
						end proc: 
						 
						 
						xkcd_line := proc(L::list, a::nonnegative, lc::positive, col, {lsty::integer:=1}) 
						  # L  : list which contains the two ending point 
						  # a  : amplitude of the random perturbations 
						  # lc : correlation length 
						  # col: color 
						  local T, roll, NX, DX, DY, LX, A, m, M, X, Z, P: 
						  T    := (a, x0, y0, l) ->  
						             plottools:-transform( 
						               (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ] 
						             ): 
						  roll := rand(-1.0 .. 1.0): 
						  NX   := 5: 
						  DX   := L[2][1]-L[1][1]: 
						  DY   := L[2][2]-L[1][2]: 
						  LX := sqrt(DX^2+DY^2): 
						  if DX <> 0 then  
						     A := arcsin(DY/LX): 
						  else 
						     A:= Pi/2: 
						  end if: 
						  X := [seq(0..1, 1/(NX-1))]: 
						  Z := [ seq(a*roll(), k=1..NX)]: 
						  P := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained, linestyle=lsty): 
						  return T(A, op(L[1]), LX)(P) 
						end proc: 
						 
						 
						xkcd_func := proc(f, r::list, NX::posint, a::positive, lc::positive, col) 
						  # f  : function to draw 
						  # r  : plot range 
						  # NX : number of equidistant "nodes" in the range r (boundaries included) 
						  # a  : amplitude of the random perturbations 
						  # lc : correlation length 
						  # col: color 
						  local roll, F, LX, Pf, Xf, Zf: 
						  roll := rand(-1.0 .. 1.0): 
						  F    := unapply(f, indets(f, name)[1]); 
						  LX   := r[2]-r[1]: 
						  Pf   := [seq(r[1]..r[2], LX/(NX-1))]: 
						  Xf   := Pf +~ [seq(a*roll(), k=1..numelems(Pf))]: 
						  Zf   := F~(Pf) +~ [seq(a*roll(), k=1..numelems(Pf))]: 
						  return plot(KG(Xf, Zf, lc*LX, 1), x=min(Xf)..max(Xf), color=col): 
						end proc: 
						 
						 
						 
						 
						xkcd_hist := proc(H, ah, av, ax, ay, lch, lcv, lcx, lcy, colh, colxy) 
						  # H   : Histogram 
						  # ah  : amplitude of the random perturbations on the horizontal boundaries of the bins 
						  # av  : amplitude of the random perturbations on the vertical boundaries of the bins 
						  # ax  : amplitude of the random perturbations on the horizontal axis 
						  # ay  : amplitude of the random perturbations on the vertical axis 
						  # lch : correlation length on the horizontal boundaries of the bins 
						  # lcv : correlation length on the vertical boundaries of the bins 
						  # lcx : correlation length on the horizontal axis 
						  # lcy : correlation length on the vertical axis 
						  # colh: color of the histogram 
						  # col : color of the axes 
						  local data, horiz, verti, horizontal_lines, vertical_lines, po, rpo, p1, p2: 
						  data  := op(1..-2, op(1, H)): 
						  verti := sort( [seq(data[n][3..4][], n=1..numelems([data]))] , key=(x->x[1]) ); 
						  verti := verti[1],  
						           map( 
						                n -> if verti[n][2] > verti[n+1][2] then  
						                        verti[n]  
						                      else  
						                        verti[n+1]  
						                      end if,  
						                [seq(2..numelems(verti)-2,2)] 
						           )[],  
						           verti[-1]; 
						  horiz := seq(data[n][[4, 3]], n=1..numelems([data])): 
						 
						  horizontal_lines := NULL: 
						  for po in horiz do 
						    horizontal_lines := horizontal_lines, xkcd_hline(po[1], po[2], ah, lch, colh): 
						  end do: 
						 
						  vertical_lines := NULL: 
						  for po in [verti] do 
						    rpo := po[[2, 1]]: 
						    vertical_lines := vertical_lines, xkcd_hline([0, rpo[2]], rpo, av, lcv, colh): 
						  end do: 
						 
						  p1 := [2*verti[1][1]-verti[2][1], 0]: 
						  p2 := [2*verti[-1][1]-verti[-2][1], 0]: 
						 
						  return  
						    display( 
						      horizontal_lines,  
						      T~([vertical_lines]),  
						      xkcd_hline(p1, p2, ax, lcx, colxy),  
						      T(xkcd_hline([0, 0], [1.2*max(op~(2, [verti])), 0], ay, lcy, colxy)),  
						      axes=none,  
						      scaling=unconstrained 
						    ); 
						end proc: 
						 
						 
						xkcd_polyline := proc(L::list, a::nonnegative, lc::positive, col) 
						  # xkcd_polyline reduces to xkcd_line whebn L has 2 elements 
						  # L  : List of points 
						  # a  : amplitude of the random perturbations 
						  # lc : correlation length 
						  # col: color 
						  local T, roll, NX, n, DX, DY, LX, A, m, M, X, Z, P: 
						  T    := (a, x0, y0, l) ->  
						             plottools:-transform( 
						               (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ] 
						             ): 
						  roll := rand(-1.0 .. 1.0): 
						  NX   := 5: 
						  for n from 1 to numelems(L)-1 do 
						    DX   := L[n+1][1]-L[n][1]: 
						    DY   := L[n+1][2]-L[n][2]: 
						    LX := sqrt(DX^2+DY^2): 
						    if DX <> 0 then  
						      A := evalf(arcsin(abs(DY)/LX)): 
						      if DX >= 0 and DY <= 0 then A := -A end if: 
						      if DX <= 0 and DY >  0 then A := Pi-A end if: 
						      if DX <= 0 and DY <= 0 then A := Pi+A end if: 
						    else 
						      A:= Pi/2: 
						      if DY < 0 then A := 3*Pi/2 end if: 
						    end if: 
						    if n=1 then 
						      X := [seq(0..1, 1/(NX-1))]: 
						      Z := [seq(a*roll(), k=1..NX)]: 
						    else 
						      X := [0    , seq(1/(NX-1)..1, 1/(NX-1))]: 
						      Z := [Z[NX], seq(a*roll(), k=1..NX-1)]: 
						    end if: 
						    P    := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained): 
						    P||n := T(A, op(L[n]), LX)(P): 
						  end do; 
						  return seq(P||n, n=1..numelems(L)-1) 
						end proc: 
						 
						 
						xkcd_circle := proc(a::nonnegative, lc::positive, r::positive, cent::list, col) 
						  # a   : amplitude of the random perturbations 
						  # lc  : correlation length 
						  # r   : redius of the circle 
						  # cent: center of the circle 
						  # col : color 
						  local roll, NX, LX, X, Z, xkg, A: 
						  roll := rand(-1.0 .. 1.0): 
						  NX   := 10: 
						  X    := [seq(0..1, 1/(NX-1))]: 
						  Z    := [0, seq(a*roll(), k=1..NX-1)]:  
						  xkg  := KG(X, Z, lc, 1): 
						  A    := Pi*roll(): 
						  return plot([cent[1]+r*(1+xkg)*cos(2*Pi*x+A), cent[2]+r*(1+xkg)*sin(2*Pi*x+A), x=0..1], color=col) 
						end proc: 
						 
						T := plottools:-transform((x,y) -> [y, x]): 
						  
						 | 
					 
				
			 
			
				
					
						| >  | 
						
						 # Axes plot 
						 
						x_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black): 
						y_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black): 
						display(  
						  x_axis,  
						  T(y_axis), 
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
			
				
					
						| >  | 
						
						 # A simple function 
						 
						f := 1+10*(x/5-1)^2: 
						F := xkcd_func(f, [0.5, 9.5], 6, 0.3, 0.4, red): 
						 
						display(  
						  x_axis,  
						  T(y_axis),  
						  F,  
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
			
				
					
						| >  | 
						
						 # An histogram 
						 
						S := Sample(Normal(0,1),100): 
						H := Histogram(S, maxbins=6): 
						xkcd_hist(H,   0, 0.02, 0.001, 0.01,   1, 0.1, 0.01, 1,   red, black) 
						 | 
					 
				
			 
			
			
				
					
						| >  | 
						
						 # Axes plus grid with two red straight lines 
						 
						r := rand(-0.1 .. 0.1): 
						 
						x_axis := xkcd_line([[-2, 0], [10, 0]], 0.01, 0.2, black): 
						y_axis := xkcd_line([[0, -2], [0, 10]], 0.01, 0.2, black): 
						d1     := xkcd_line([[-1, 1], [9, 9]] , 0.01, 0.2, red): 
						d2     := xkcd_line([[-1, 9], [9, -1]], 0.01, 0.2, red): 
						display(  
						  x_axis, y_axis, 
						  seq( xkcd_line([[-2+0.3*r(), u+0.3*r()], [10+0.3*r(), u+0.3*r()]], 0.005, 0.5, gray), u in [seq(-1..9, 2)]), 
						  seq( xkcd_line([[u+0.3*r(), -2+0.3*r()], [u+0.3*r(), 10+0.3*r()]], 0.005, 0.5, gray), u in [seq(-1..9, 2)]), 
						  d1, d2, 
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
			
				
					
						| >  | 
						
						 # Axes and a couple of polygonal lines 
						 
						d1 := xkcd_polyline([[0, 0], [1, 3], [3, 5], [7, 1], [9, 7]], 0.01, 1, red): 
						d2 := xkcd_polyline([[0, 9], [2, 8], [5, 2], [8, 3], [10, -1]], 0.01, 1, blue): 
						 
						display(  
						  x_axis, y_axis, 
						  d1, d2, 
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
			
				
					
						| >  | 
						
						 # A few polygonal shapes 
						 
						display(  
						  xkcd_polyline([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]], 0.01, 1, red), 
						  xkcd_polyline([[1/3, 1/3], [2/3, 1/3], [2/3, 4/3], [-1, 4/3], [1/3, 1/3]], 0.01, 1, blue), 
						  xkcd_polyline([[-1/3, -1/3], [4/3, 1/2], [1/2, 1/2], [1/2,-1], [-1/3, -1/3]], 0.01, 1, green), 
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
			
				
					
						| >  | 
						
						 # A few circles 
						 
						cols  := [red, green, blue, gold, black]:                                # colors 
						cents := convert( Statistics:-Sample(Uniform(-1, 3), [5, 2]), listlist): # centers 
						radii := Statistics:-Sample(Uniform(1/2, 2), 5):                         # radii 
						lcs   := Statistics:-Sample(Uniform(0.2, 0.7), 5):                       # correlations lengths 
						 
						display(  
						  seq( 
						    xkcd_circle(0.02, lcs[n], radii[n], cents[n], cols[n]), 
						    n=1..5 
						  ), 
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
				
					
						| 
						   
						 | 
						  | 
					 
				
			 
			
				
					
						| >  | 
						
						 # A 3D drawing 
						 
						x_axis := xkcd_line([[0, 0], [5, 0]], 0.01, 0.2, black): 
						y_axis := xkcd_line([[0, 0], [4, 2]], 0.01, 0.2, black): 
						z_axis := xkcd_line([[0, 0], [0, 5]], 0.01, 0.2, black): 
						 
						f1 := 4*cos(x/6)-1: 
						F1 := xkcd_func(f1, [0.5, 5], 6, 0.001, 0.8, red): 
						F2 := xkcd_line([[0.5, eval(f1, x=0.5)], [3, 4]], 0.01, 0.2, red): 
						f3 := 4*cos((x-2)/6): 
						F3 := xkcd_func(f3, [3, 7], 6, 0.001, 0.8, red): 
						F4 := xkcd_line([[5, eval(f1, x=5)], [7, eval(f3, x=7)]], 0.01, 0.2, red): 
						 
						dx := xkcd_line([[2, 1], [4, 1]], 0.01, 0.2, gray, lsty=3): 
						dy := xkcd_line([[2, 0], [4, 1]], 0.01, 0.2, gray, lsty=3): 
						dz := xkcd_line([[4, 1], [4, 3]], 0.01, 0.2, gray, lsty=3): 
						 
						po := xkcd_circle(0.02, 0.3, 0.1, [4, 3], blue): 
						 
						# Numerical value come from "probe info + copy/paste" 
						 
						nvect   := xkcd_polyline([[4, 3], [4.57, 4.26], [4.35, 4.1], [4.57, 4.26], [4.58, 4.02]], 0.01, 1, blue): 
						tg1vect := xkcd_polyline([[4, 3], [4.78, 2.59], [4.49, 2.87], [4.78, 2.59], [4.46, 2.57]], 0.01, 1, blue): 
						tg2vect := xkcd_polyline([[4, 3], [4.79, 3.35], [4.70, 3.13], [4.79, 3.35], [4.46, 3.35]], 0.01, 1, blue): 
						rec1    := xkcd_polyline([[4.118, 3.286], [4.365, 3.396], [4.257, 3.108]], 0.01, 1, blue): 
						rec2    := xkcd_polyline([[4.257, 3.108], [4.476, 2.985], [4.259, 2.876]], 0.01, 1, blue): 
						 
						 
						 
						display(  
						  x_axis, y_axis, z_axis, 
						  F1, F2, F3, F4, 
						  dx, dy, dz, 
						  po, 
						  nvect, tg1vect, tg2vect, rec1, rec2, 
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
				
					
						| 
						   
						 | 
						  | 
					 
				
			 
			
				
					
						| >  | 
						
						 # Arrow 
						 
						d1 := xkcd_polyline([[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, -0.05]], 0.01, 1, red): 
						 
						 
						T := (a, x0, y0, l) ->  
						             plottools:-transform( 
						               (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ] 
						             ): 
						 
						 
						display(  
						  seq( T(2*Pi*n/10, 0.5, 0, 1/2)( 
						           display( 
						              xkcd_polyline( 
						                  [[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, -0.05]],  
						                  0.01,  
						                  1,  
						                  ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]) 
						               ) 
						           ) 
						        ),  
						       n=1..10 
						  ), 
						  axes=none,  
						  scaling=constrained  
						) 
						 | 
					 
				
			 
			
			
			 |