kepler_tex/tikz/tikz.ml

   1 (* ========================================================================== *)
   2 (* FLYSPECK - BOOK PREPARATION                                                *)
   3 (*                                                                            *)
   4 (* Chapter: Graphics                                                          *)
   5 (* Author: Thomas C. Hales                                                    *)
   6 (* Date: 2011-11-26                                                           *)
   7 (* ========================================================================== *)
   8
   9 (*
  10 Some procedures to facilitate the generation of tikz graphics.
  11 Tikz.execute produces output in /tmp/x.txt
  12 Read in by fig.tex to produce graphics.
  13
  14 This is mostly independent of the flyspeck .hl files, but it might make
  15 light use of HOL's lib.ml functions, such as sort.
  16
  17 Lessons learned:
  18
  19 Tikz is almost totally unusable for 3D graphics.
  20 Generate 3D coordinates with OCAML, then project to 2D at the very end.
  21 Use --plot[smooth] coordinates { ... } rather than try to use tikz elliptical
  22 arc routines.
  23
  24    At the beginning, I used Mathematica to generate tikz files, but
  25    eventually almost everything was done with Ocaml rather than Math'ca.
  26
  27 *)
  28
  29 module Tikz = struct
  30
  31 let map = List.map;;
  32
  33 let filter= List.filter;;
  34
  35 let hd = List.hd;;
  36
  37 (* from HOL LIGHT lib.ml . *)
  38
  39 (* let (o) = fun f g x -> f(g x);; *)
  40
  41 let rec (--) = fun m n -> if m > n then [] else m::((m + 1) -- n);;
  42
  43 let rec funpow n f x =
  44   if n < 1 then x else funpow (n-1) f (f x);;
  45
  46 let rec forall p l =
  47   match l with
  48     [] -> true
  49   | h::t -> p(h) & forall p t;;
  50
  51 let rec mem x lis =
  52   match lis with
  53     [] -> false
  54   | (h::t) -> Pervasives.compare x h = 0 or mem x t;;
  55
  56 let subtract l1 l2 = filter (fun x -> not (mem x l2)) l1;;
  57
  58 let intersect l1 l2 = filter (fun x -> mem x l2) l1;;
  59
  60 let subset l1 l2 = forall (fun t -> mem t l2) l1;;
  61
  62 let rec partition p l =
  63   match l with
  64     [] -> [],l
  65   | h::t -> let yes,no = partition p t in
  66             if p(h) then (if yes == t then l,[] else h::yes,no)
  67             else (if no == t then [],l else yes,h::no);;
  68
  69 let rec sort cmp lis =
  70   match lis with
  71     [] -> []
  72   | piv::rest ->
  73       let r,l = partition (cmp piv) rest in
  74       (sort cmp l) @ (piv::(sort cmp r));;
  75
  76 let rec zip l1 l2 =
  77   match (l1,l2) with
  78         ([],[]) -> []
  79       | (h1::t1,h2::t2) -> (h1,h2)::(zip t1 t2)
  80       | _ -> failwith "zip";;
  81
  82 let rec end_itlist f l =
  83   match l with
  84         []     -> failwith "end_itlist"
  85       | [x]    -> x
  86       | (h::t) -> f h (end_itlist f t);;
  87
  88 let rec assoc a l =
  89   match l with
  90     (x,y)::t -> if Pervasives.compare x a = 0 then y else assoc a t
  91   | [] -> failwith "find";;
  92
  93
  94 (* I/O *)
  95
  96 let output_filestring tmpfile a =
  97   let outs = open_out tmpfile in
  98   let _ = try (Printf.fprintf outs "%s" a)
  99   with _ as t -> (close_out outs; raise t) in
 100    close_out outs ;;
 101
 102 let unsplit d f = function
 103   | (x::xs) ->  List.fold_left (fun s t -> s^d^(f t)) (f x) xs
 104   | [] -> "";;
 105
 106 let join_comma  = unsplit "," (fun x-> x);;
 107
 108 let join_lines  = unsplit "\n" (fun x-> x);;
 109
 110 let join_space  = unsplit " " (fun x-> x);;
 111
 112
 113
 114 (* math *)
 115
 116 let cos = Pervasives.cos;;
 117 let sin = Pervasives.sin;;
 118 let cot x = cos x /. sin x;;
 119 let sqrt = Pervasives.sqrt;;
 120 let pi = 4.0 *. atan(1.0);;
 121 let nth = List.nth;;
 122
 123 (* arg between 0 and 2pi *)
 124
 125 let arg x y = if (y<0.0) then atan2 y x +. 2.0 *. pi else atan2 y x;;
 126
 127 let degree x = 180.0 *. x /. pi;;
 128
 129 let radian x = pi *. x /. 180.0;;
 130
 131 let eta x y z =
 132   let s = (x +. y +. z)/. 2.0 in
 133    x *. y *. z /. ( 4. *. sqrt(s *. (s -. x) *. ( s -. y) *. (s -. z)));;
 134
 135 let orig3 = (0.0,0.0,0.0);;
 136
 137 let orig2 = (0.0,0.0);;
 138
 139 (* vector sum, difference, scalar product, dot product *)
 140
 141 let map3 f (x,y,z) = (f x,f y,f z);;
 142
 143 let map2 f (x,y) = (f x , f y);;
 144
 145 let (+..) (x1,x2) (y1,y2) = (x1+. y1,x2+. y2);;
 146
 147 let (-..)  (x1,x2) (y1,y2) = (x1-. y1,x2-. y2);;
 148
 149 let uminus3 (x1,x2,x3) = (-. x1,-.x2,-.x3);;
 150
 151 let uminus2 (x1,x2) = (-. x1,-.x2);;
 152
 153 let ( %.. ) s (x1,x2) = (s *. x1, s *. x2);;
 154
 155 let ( *.. ) (x1,x2) (y1,y2) = (x1 *. y1 +. x2 *. y2);;
 156
 157 let (+...) (x1,x2,x3) (y1,y2,y3) = (x1 +. y1, x2 +. y2, x3+. y3);;
 158
 159 let (-...) (x1,x2,x3) (y1,y2,y3) = (x1 -. y1, x2 -. y2, x3-. y3);;
 160
 161 let ( %... ) s (x1,x2,x3) = (s *. x1, s *. x2, s*. x3);;
 162
 163 let ( *... ) (x1,x2,x3) (y1,y2,y3) = (x1 *. y1 +. x2 *. y2 +. x3 *. y3);;
 164
 165 let cross (x1,x2,x3) (y1,y2,y3) =
 166   (x2 *. y3 -. x3 *. y2, x3 *. y1 -. x1 *. y3, x1 *. y2 -. x2 *. y1);;
 167
 168 let det3 x y z = x *... (cross y z);;
 169
 170 let det2 (x1,y1) (x2,y2) = (x1 *. y2 -. y1 *. x2);;
 171
 172 let conj (x,y) = (x,-. y);;
 173
 174 let cmul (x1,y1) (x2,y2) = (x1 *. x2 -. y1 *. y2, x1 *. y2 +. x2 *. y1);;
 175
 176 let cinv v = (1.0/. (v *.. v)) %.. (conj v);;
 177
 178 let cdiv u v = cmul u (cinv v);;
 179
 180
 181
 182 let delta1 = (1.0,0.0,0.0);;
 183
 184 let delta2 = (0.0,1.0,0.0);;
 185
 186 let delta3 = (0.0,0.0,1.0);;
 187
 188 let proj e1 e2 x = (x *... e1) , (x *... e2);;
 189
 190 let perp p x  =  x -... (((x *... p) /. (p *... p)) %... p) ;; (* ortho to p *)
 191
 192 let transpose ((a11,a12,a13),(a21,a22,a23),(a31,a32,a33)) =
 193   ((a11,a21,a31),(a12,a22,a32),(a13,a23,a33));;
 194
 195 let transpose2 ((x1,y1),(x2,y2)) = ((x1,x2),(y1,y2));;
 196
 197 let mul3 (e1,e2,e3) x =
 198      (e1 *... x,  e2 *... x, e3 *... x);;
 199
 200 let tuple3 [v1;v2;v3] = (v1,v2,v3);;
 201
 202 let list3 (v1,v2,v3) = [v1;v2;v3];;
 203
 204 let tuple2 [v1;v2] = (v1,v2);;
 205
 206 let list2 (v1,v2) = [v1;v2];;
 207
 208 let norm2 x = sqrt(x *.. x);;
 209
 210 let norm3 x = sqrt(x *... x);;
 211
 212 let normalize3 x = (1.0 /. sqrt(x *... x)) %... x;;
 213
 214 let normalize2 x = (1.0 /. sqrt(x *.. x)) %.. x;;
 215
 216 let dist3 x y =
 217   let z = x -... y in sqrt (z *... z);;
 218
 219 let dist2 x y =
 220   let z = x -.. y in sqrt (z *.. z);;
 221
 222 let rec outer x y =
 223   match x with
 224     | [] -> []
 225     | a::r -> (map (fun i -> (a,i)) y) @ (outer r y);;
 226
 227 let solve33 (m1,m2,m3) c =    (* solve m.x ==c for x by Cramer *)
 228   let d = det3 m1 m2 m3 in
 229   let (t1,t2,t3) = transpose (m1,m2,m3) in
 230    map3 (fun t -> t/. d) (det3 c t2 t3, det3 t1 c t3, det3 t1 t2 c);;
 231
 232 let solve22 (m1,m2) c =
 233   let d = det2 m1 m2 in
 234   let (t1,t2) = transpose2 (m1,m2) in
 235    map2 (fun t -> t/. d) (det2 c t2, det2 t1 c);;
 236
 237 let extreme_point m' =
 238   solve33 m' (map3 (fun m -> 0.5 *. (m *... m)) m');;
 239
 240 let lex3 (i,j,k) (i',j',k') = (i<i') or ((i=i') &&(j<j')) or ((i=i')&&(j=j')&&(k<k'));;
 241
 242 let etaV u v w =
 243   let x = dist2 (u -.. v) orig2 in
 244   let y = dist2 (v -.. w) orig2 in
 245   let z = dist2 (u -.. w) orig2 in
 246    eta x y z;;
 247
 248 (* BUGGED SOMEHOW.
 249 let circum3 (a,b,c) =
 250   let [a';b';c']= map norm3 [a;b;c] in
 251   let e = eta a' b' c' in
 252   let u = a -... c in
 253   let v = b -... c in
 254   let w = (cross u v) in
 255   let n = normalize3 (cross w u) in
 256   let t = sqrt (e *. e -. a' *. a' /. 4.0) in
 257     c +... (0.5 %... u) +... (t %... n);;
 258 *)
 259
 260 let circumcenter (a,b,c) =
 261   let a' = a -.. c in
 262   let b' = b -.. c in
 263   c +.. (0.5 %.. (solve22 (a',b') (a' *.. a', b' *.. b')));;
 264
 265 let frame_of v1 v2 =
 266   let e1 = normalize3 (v1) in
 267   let e2 = normalize3 (perp e1 v2) in
 268   let e3 = cross e1 e2 in
 269     (e1,e2,e3);;
 270
 271 let random_SO3 seed =
 272   let _ = Random.init seed in
 273   let v3 () = tuple3 (map (fun _ -> -1.0 +. Random.float 2.0) [0;0;0]) in
 274     frame_of (v3()) (v3());;
 275
 276
 277 (* TIKZ OUTPUT *)
 278
 279 let ppair (x,y) = Printf.sprintf "(%f,%f)" x y;;
 280
 281 let pcoord s (x,y) =
 282   Printf.sprintf "\\coordinate (%s) at (%f,%f) " s x y;;
 283
 284 let plet s y =
 285   Printf.sprintf "\\pgfmathsetmacro\\%s{%s}" s y;;
 286
 287
 288
 289 (* specific cases *)
 290
 291 (* CLOSE PACKING CHAPTER FIGURES *)
 292
 293 (* SEYIMIE *)
 294
 295 let fcc_fun_domain =
 296   let v1 = (1.0,0.0,0.0) in
 297   let v2 = (0.5,sqrt(3.0)/.2.0,0.0) in
 298   let v3 = (0.5,1.0 /. sqrt(12.0),sqrt(2.0/.3.0)) in
 299   let v12 = v1 +... v2 in
 300   let v23 = v2 +... v3 in
 301   let v13 = v1 +... v3 in
 302   let v123 = v1 +... v2 +... v3 in
 303   let f = frame_of (1.0,0.1,0.1) (0.3,0.5,1.0) in (* 0.3 0.5 1.0 *)
 304   let p v = proj delta1 delta2 (mul3 f v) in
 305   let [w1;w2;w3;w12;w23;w13;w123] = map p [v1;v2;v3;v12;v23;v13;v123] in
 306   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 307     join_lines (map coord [("w1",w1);("w2",w2);("w3",w3);("w12",w12);("w23",w23);("w13",w13);("w123",w123)]);;
 308
 309 (*  figAZGXQWC *)
 310
 311 let tet_oct_ratio =
 312   let v0 = (0.0,0.0,0.0) in
 313   let v1 = (1.0,0.0,0.0) in
 314   let v2 = (0.5,sqrt(3.0)/.2.0,0.0) in
 315   let v3 = (0.5,1.0 /. sqrt(12.0),sqrt(2.0/.3.0)) in
 316   let v12 = 0.5 %... (v1 +... v2) in
 317   let v02 = 0.5 %...( v0 +... v2) in
 318   let v01 = 0.5 %... v1 in
 319   let v03 = 0.5 %... v3 in
 320   let v13 = 0.5 %... (v1 +... v3) in
 321   let v23 = 0.5 %... (v2 +... v3) in
 322   let f = frame_of (1.0,0.1,0.1) (0.3,0.5,1.0) in (* 0.3 0.5 1.0 *)
 323   let p v = proj delta1 delta2 (mul3 f v) in
 324   let  [w0;w1;w2;w3;w12;w02;w01;w03;w13;w23] = map p [v0;v1;v2;v3;v12;v02;v01;v03;v13;v23] in
 325   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 326     join_lines (map coord [("w0",w0);("w1",w1);("w2",w2);("w3",w3);("w12",w12);("w02",w02);("w01",w01);("w03",w03);("w13",w13);("w23",w23)]);;
 327
 328 (* SGIWBEN *)
 329
 330 let fcc_hcp_pattern =
 331   let f = frame_of
 332 (*     (0.3,0.4,0.1) (-0.1,0.1,0.4) in *)
 333      (0.4,0.3,0.1) (-0.2,0.1,0.4) in
 334   let g = mul3 f in
 335   let u = g delta3 in
 336   let v0 = (0.0,0.0,0.0) in
 337   let v1 = g(1.0,0.0,0.0) in
 338   let v2 = g(0.5,sqrt(3.0)/.2.0,0.0) in
 339   let v3 = g(0.5,1.0 /. sqrt(12.0),sqrt(2.0/.3.0)) in
 340   let v4 = v2 -... v1 in
 341   let v5 = v0 -... v1 in
 342   let v6 = v0 -... v2 in
 343   let v7 = v0 -... v4 in
 344   let v8 = v3 +... v5 in
 345   let v9 = v3 +... v6 in
 346   let v10 = v0 -... v3 in
 347   let v11 = v0 -... v8 in
 348   let v12 = v0 -...  v9 in
 349   let n v = v -... (2.0 *. (u *... v) ) %... u in
 350   let v13 = v0 -... n v3 in
 351   let v14 = v0 -... n v8 in
 352   let v15 = v0 -... n v9 in
 353   let p v = proj delta1 delta2 (v) in
 354   let  [w0;w1;w2;w3;w4;w5;w6;w7;w8;w9;w10;w11;w12;w13;w14;w15] = map p
 355     [v0;v1;v2;v3;v4;v5;v6;v7;v8;v9;v10;v11;v12;v13;v14;v15] in
 356   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 357     join_lines (map coord [("w0",w0);("w1",w1);("w2",w2);("w3",w3);
 358      ("w4",w4);("w5",w5);("w6",w6);("w7",w7);("w8",w8);("w9",w9);("w10",w10);("w11",w11);("w12",w12);("w13",w13);("w14",w14);("w15",w15)]);;
 359
 360
 361 let fcc_packing =
 362   let f = frame_of
 363 (*     (0.4,0.3,0.1) (-0.2,0.1,0.4) in  *)
 364      (0.5,0.4,0.) (-. 0.0,0.1,0.4) in
 365   let g = mul3 f in
 366   let u = g delta3 in
 367   let v0 = (0.0,0.0,0.0) in
 368   let v1 = g(1.0,0.0,0.0) in
 369   let v2 = g(0.5,sqrt(3.0)/.2.0,0.0) in
 370   let e1 = g delta1 in
 371   let e2 = g delta2 in
 372   let e3 = u in
 373   let v3 = g(0.5,1.0 /. sqrt(12.0),sqrt(2.0/.3.0)) in
 374   let e12 = e1 +... e2 in
 375   let e13 = e1 +... e3 in
 376   let e23 = e2 +... e3 in
 377   let e123 = e1 +... e2 +... e3 in
 378   let p v = proj delta1 delta2 (v) in
 379   let  [w0;w1;w2;w3;e1;e2;e3;e12;e13;e23;e123] = map p [v0;v1;v2;v3;e1;e2;u;e12;e13;e23;e123] in
 380   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 381     join_lines (map coord [("w0",w0);("w1",w1);("w2",w2);("w3",w3);("e1",e1);("e2",e2);("e3",e3);("e12",e12);("e13",e13);("e23",e23);("e123",e123)]);;
 382
 383 let pascal_packing =
 384   let f = frame_of
 385 (*     (0.5,0.4,0.0) (-0.0,0.1,0.4) in  *)
 386      (0.5,0.4,0.) (-0.2,0.1,0.4) in
 387   let g = mul3 f in
 388   let v0 = (0.0,0.0,0.0) in
 389   let v1 = g(1.0,0.0,0.0) in
 390   let v2 = g(0.5,sqrt(3.0)/.2.0,0.0) in
 391   let v3 = g(0.5,1.0 /. sqrt(12.0),sqrt(2.0/.3.0)) in
 392   let p v = proj delta1 delta2 (v) in
 393   let  [w0;w1;w2;w3] = map p [v0;v1;v2;v3] in
 394   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 395     join_lines (map coord [("w0",w0);("w1",w1);("w2",w2);("w3",w3)]);;
 396
 397
 398
 399
 400 (* TCFVGTS % fig:face-centered-cubic *)
 401
 402 let circle_point r u v =
 403     u +... r %... normalize3 (v -... u);;
 404
 405 let circle_interpolate r s u v1 v2 =
 406   let v = (s %... v2) +... ((1.0 -. s) %... v1) in
 407     circle_point r u v;;
 408
 409 let pcircle r n v u1 u2 label =
 410   let p1 = map (fun s -> circle_interpolate r (float_of_int s /. float_of_int n) v u1 u2) (0--n) in
 411   let q1 = map (proj delta1 delta2) p1 in
 412   let w1 = join_space (map (fun (x,y)-> Printf.sprintf "(%f,%f) " x y) q1) in
 413     Printf.sprintf "\\def\\%s{%s}" label w1 ;;
 414
 415 let cubic_layers =
 416   let f = frame_of (1.0,0.1,0.4) (-0.5,1.0,0.0) in
 417   let g = mul3 f in
 418   let r = sqrt(8.0) in
 419   let v0 =  (0.0,0.0,0.0) in
 420   let v1 = g(r %... delta1) in
 421   let v2 = g(r %... delta2) in
 422   let v3 = g(r %... delta3) in
 423   let v12 = (v1 +... v2) in
 424   let v13 = ( v1 +... v3) in
 425   let v23 = (v2+... v3) in
 426   let v123 = (v1 +... v2 +... v3) in
 427   let vfront = (0.5 %... (v12)) in
 428   let vtop = (0.5 %... (v2 +... v123)) in
 429   let vright = (0.5 %... (v1 +... v123)) in
 430   let p v = proj delta1 delta2 v in
 431   let   [w0;w1;w2;w3;w12;w13;w23;w123;wfront;wtop;wright]  = map p
 432     [v0;v1;v2;v3;v12;v13;v23;v123;vfront;vtop;vright] in
 433   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 434   let cc =   (map coord
 435                           [("w0",w0);("w1",w1);("w2",w2);("w3",w3);("w12",w12);("w13",w13);
 436                            ("w23",w23);("w123",w123);("wfront",wfront);("wtop",wtop);
 437                            ("wright",wright)]) in
 438   let b s = "tcf"^s in
 439   let paths = [(v0,v1,v2,b "a");(v1,v12,v0,b "b");(v12,v2,v1,b "c");
 440               (v2,v0,v12,b "d");
 441               (vfront,v12,v2,b"e");(vfront,v2,v0,b"f");(vfront,v0,v1,b"g");
 442               (vfront,v1,v12,b"h");
 443               (v1,v13,v12,b"i");(v13,v123,v1,b"j");(v123,v12,v13,b"k");
 444               (v12,v1,v123,b"l");
 445               (vright,v123,v12,b"m");(vright,v12,v1,b"n");(vright,v1,v13,b"o");
 446               (vright,v13,v123,b"p");
 447               (v2,v12,v23,b"q");(v12,v123,v2,b"r");(v123,v12,v23,b"s");
 448               (v23,v123,v2,b"t");
 449               (vtop,v12,v123,b"u");(vtop,v123,v23,b"v");(vtop,v23,v2,b"w");
 450               (vtop,v2,v12,b"x");] in
 451   let pc = map (fun (u,v1,v2,s)-> pcircle 1.0 5 u v1 v2 s) paths in
 452     join_lines (cc @ pc);;
 453
 454
 455 (* NTNKMGO *)
 456
 457 let square_layers =
 458   let f = frame_of (1.0,0.1,0.4) (-0.5,1.0,0.0) in
 459   let g = mul3 f in
 460   let v0 =  (0.0,0.0,0.0) in
 461   let v1 = g( delta1) in
 462   let v2 = g( delta2) in
 463   let v3 = g( delta3) in
 464   let p v = proj delta1 delta2 v in
 465   let   [w0;w1;w2;w3]  = map p   [v0;v1;v2;v3] in
 466   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 467   let cc = map coord  [("w0",w0);("w1",w1);("w2",w2);("w3",w3)] in
 468     join_lines (cc);;
 469
 470
 471
 472 (* PQJIJGE *)
 473
 474 let rhombic_dodec =
 475   let f = frame_of (1.0,0.1,0.4) (-0.5,1.0,0.0) in
 476   let g = mul3 f in
 477   let r = 2.0 in
 478   let v0 =  (0.0,0.0,0.0) in
 479   let v1 = g(r %... delta1) in
 480   let v2 = g(r %... delta2) in
 481   let v3 = g(r %... delta3) in
 482   let v12 = (v1 +... v2) in
 483   let v13 = ( v1 +... v3) in
 484   let v23 = (v2+... v3) in
 485   let v123 = (v1 +... v2 +... v3) in
 486   let center = 0.5 %... v123 in
 487   let vfront = (0.5 %... (v12 -... v3)) in
 488   let vtop = (0.5 %... (v2 +... v123)) +... 0.5 %... v2 in
 489   let vright = (0.5 %... (v1 +... v123)) +... 0.5 %... v1 in
 490   let vback = vfront +... 2.0 %... v3 in
 491   let vleft = vright -... 2.0 %... v1 in
 492   let vbottom = vtop -... 2.0 %... v2 in
 493   let p v = proj delta1 delta2 v in
 494   let   [w0;w1;w2;w3;w12;w13;w23;w123;wfront;wtop;wright;wback;wleft;wbottom;wcenter]  = map p
 495     [v0;v1;v2;v3;v12;v13;v23;v123;vfront;vtop;vright;vback;vleft;vbottom;center] in
 496   let coord (s,u) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s (fst u) (snd u) in
 497   let cc =   (map coord
 498                           [("w0",w0);("w1",w1);("w2",w2);("w3",w3);("w12",w12);("w13",w13);
 499                            ("w23",w23);("w123",w123);("wfront",wfront);("wtop",wtop);
 500                            ("wright",wright);
 501                             ("wback",wback);("wleft",wleft);("wbottom",wbottom);
 502                           ("wcenter",wcenter)]) in
 503     join_lines cc;;
 504
 505
 506
 507
 508
 509 (* PACKING CHAPTER FIGURES *)
 510
 511 (* figDEQCVQL *)
 512
 513 let rec mindist2 r w = function
 514   | [] -> r
 515   | l::ls -> if (dist2 l w < r) then mindist2 (dist2 l w) w ls else mindist2 r w ls;;
 516
 517 let randompacking radius seed xdim ydim = (* seed=5 works well *)
 518   let _ = Random.init seed in
 519 (*  let radius = 0.15 in *)
 520   let v () = (Random.float xdim,Random.float ydim) in
 521   let add len ls =
 522     let w = v() in
 523       if (mindist2 100.0 w ls < 2.0 *. radius) or List.length ls > len then ls else w::ls in
 524   let unsat = funpow 40 (add 15) [] in
 525   let sat = funpow 100000 (add 20000) unsat in
 526    (unsat,sat);;
 527
 528 let print_satunsat seed =
 529   let radius = 0.15 in
 530   let (unsat,sat) = randompacking radius seed 2.0 2.0 in
 531   let line d (x,y)  = Printf.sprintf "\\draw[gray,fill=black!%d] (%f,%f) circle (0.15);" d x y in
 532   let punsat = map (line 30) unsat in
 533   let psat = map (line 10) sat in
 534     join_lines (["\\begin{scope}[shift={(0,0)}]"] @ punsat @
 535       ["\\end{scope}\\begin{scope}[shift={(3.0,0)}]"] @ psat @ punsat @ ["\\end{scope}"]);;
 536
 537 (* \figXOHAZWO Voronoi cells of a random saturated packing. Start with Delaunay triangles. *)
 538
 539
 540
 541 let center2 s (i,j,k) =  circumcenter (s i , s j , s k);;
 542 (*
 543   let si = s i -.. s k in
 544   let sj = s j -.. s k in
 545     s k +.. (0.5 %.. (solve22 (si,sj) (si *.. si, sj *.. sj)));;
 546 *)
 547
 548 let sat_triples sat =
 549   let radius = 0.15 in
 550   let r = List.length sat in
 551   let s = nth sat in
 552   let rr = 0--(r-1) in
 553   let allt = outer rr (outer rr rr) in
 554   let triple =  (map (fun (i,(j,k))->(i,j,k)) (filter(fun (i,(j,k))->(i<j && j<k))  allt)) in
 555   let allpair = filter (fun (i,j) -> i< j) (outer rr rr) in
 556   let shortpair = filter (fun (i,j) -> dist2 (s i) (s j) < 4.0 *. radius +. 1.0e-5) allpair in
 557   let ftriple = filter (fun (i,j,k) -> mem (i,j) shortpair && mem(i,k) shortpair && mem(j,k) shortpair) triple in
 558   let fit (i,j,k) =
 559     let c = center2 s (i,j,k) in
 560     let rad = dist2 c (s k) in
 561     let rest = subtract sat [s i;s j;s k] in
 562     let vals = filter (fun v -> dist2 v c < rad) rest in
 563       List.length vals = 0 in
 564     filter fit ftriple;;
 565
 566 let print_satst seed=
 567   let radius = 0.15 in
 568   let (_,sat) =  randompacking radius seed 5.0 2.0 in
 569   let satst = sat_triples sat in
 570   let s = nth sat in
 571   let prs = filter (fun ((i,j,k),(i',j',k')) ->
 572                       List.length (intersect [i;j;k] [i';j';k'])=2 && lex3 (i,j,k) (i',j',k'))
 573     (outer satst satst) in
 574   let pp = map (fun (t, t') ->
 575                   let (x,y) = center2 s t in
 576                   let (x',y') = center2 s t' in
 577                     Printf.sprintf "\\draw (%f,%f) -- (%f,%f) ;" x y x' y') prs in
 578   let smalldot (x,y)  = Printf.sprintf "%c\\draw[gray!10,fill=gray!30] (%f,%f) circle (0.15);\n\\smalldot{%f,%f};" '%' x y x y in
 579   let psmalldot = map (smalldot) sat in
 580     join_lines (psmalldot @ pp);;
 581
 582 (* autoBUGZBTW *)
 583
 584 let print_rogers seed=
 585   let radius = 0.15 in
 586   let (_,sat) = randompacking radius seed 3.0 1.4 in (* 3,1.4 *)
 587   let satst = sat_triples sat in
 588   let s = nth sat in
 589   let coord_triple t = center2 (nth sat) t in
 590   let prs = filter (fun ((i,j,k),(i',j',k')) ->
 591                       List.length (intersect [i;j;k] [i';j';k'])=2 && lex3 (i,j,k) (i',j',k'))
 592     (outer satst satst) in
 593   let pp = map (fun (t, t') ->
 594                   let (x,y) = coord_triple t in
 595                   let (x',y') = coord_triple t' in
 596                     Printf.sprintf "\\draw[very thick] (%f,%f) -- (%f,%f) ;" x y x' y') prs in
 597   let smalldot (x,y)  = Printf.sprintf "%c\\draw[gray!10,fill=gray!30] (%f,%f) circle (0.15);\n\\smalldot{%f,%f};" '%' x y x y in
 598   let psmalldot = map (smalldot) sat in
 599   let draw ((ux,uy),(vx,vy)) = Printf.sprintf "\\draw (%f,%f)--(%f,%f);" ux uy vx vy in
 600   let drawc (ax,ay) (bx,by) (cx,cy) (dx,dy)=
 601       Printf.sprintf "\\draw[fill=gray] (%f,%f)--(%f,%f)--(%f,%f)--(%f,%f)--cycle;" ax ay bx by cx cy dx dy in
 602   let draw_radial (i,j,k) =
 603     let c = coord_triple (i,j,k) in
 604     let (u1,u2,u3)=(s i,s j,s k) in
 605     let vv = map draw  [(c,u1);(c,u2);(c,u3);] in
 606       join_lines vv in
 607   let radial = map draw_radial satst in
 608   let draw_dedge (t,t') =
 609     let c = coord_triple t in
 610     let c'= coord_triple t' in
 611     let (i,j) = tuple2 (intersect (list3 t) (list3 t')) in
 612     let u1 = s i in
 613     let u2 = s j in
 614     let w = u2 -.. u1 in
 615     let d = c -.. u1 in
 616     let d' = c' -.. u1 in
 617       if (det2 d w  *. det2 w d' > 0.0) then draw (u1, u2) else drawc u1 c u2 c' in
 618   let dedge = map draw_dedge prs in
 619     join_lines (psmalldot @ radial @ dedge @ pp);;
 620
 621 (* EVIAIQP
 622 *)
 623 let print_voronoi seed=
 624   let radius = 0.15 in
 625   let (_,sat) = randompacking radius seed 3.0 1.4 in (* 3,1.4 *)
 626   let satst = sat_triples sat in
 627   let coord_triple t = center2 (nth sat) t in
 628   let prs = filter (fun ((i,j,k),(i',j',k')) ->
 629                       List.length (intersect [i;j;k] [i';j';k'])=2 && lex3 (i,j,k) (i',j',k'))
 630     (outer satst satst) in
 631   let pp = map (fun (t, t') ->
 632                   let (x,y) = coord_triple t in
 633                   let (x',y') = coord_triple t' in
 634                     Printf.sprintf "\\draw[very thick] (%f,%f) -- (%f,%f) ;" x y x' y') prs in
 635   let dot (x,y)  = Printf.sprintf "\\smalldot{%f,%f};" x y  in
 636   let psmalldot = map (dot) sat in
 637     join_lines (psmalldot @  pp);;
 638
 639 (* ANNTKZP *)
 640
 641 let print_delaunay seed=
 642   let radius = 0.15 in
 643   let (_,sat) = randompacking radius seed 3.0 1.4 in (* 3,1.4 *)
 644   let satst = sat_triples sat in
 645   let delaunay_edge = List.flatten (map (fun (i,j,k) -> [(i,j);(j,k);(k,i)]) satst) in
 646
 647   let s = nth sat in
 648   let smalldot (x,y)  = Printf.sprintf "\\smalldot{%f,%f};"  x y  in
 649   let psmalldot = map (smalldot) sat in
 650   let draw ((ux,uy),(vx,vy)) = Printf.sprintf "\\draw (%f,%f)--(%f,%f);" ux uy vx vy in
 651   let pdraw = map (fun (i,j) -> draw (s i,s j)) delaunay_edge in
 652     join_lines ( psmalldot @ (* radial @ dedge @  pp @ *) pdraw);;
 653
 654
 655
 656 (* figYAJOTSL *)
 657
 658 let mk_tetrahedron seed =
 659   let rot v = mul3 (random_SO3 seed) ( sqrt(3.0 /. 2.0) %... v) in
 660   let v1 = rot (0.0,0.0,1.0) in
 661   let v2 = rot (sqrt(8.0)/. 3.0,0.0,-. 1.0/. 3.0) in
 662   let v3 = rot (-. sqrt(8.0)/. 6.0, sqrt(2.0/. 3.0), -. 1.0/. 3.0) in
 663   let v4 = rot (-. sqrt(8.0)/. 6.0,-. sqrt(2.0 /. 3.0),-. 1.0/. 3.0) in
 664      (v1,v2,v3,v4);;
 665
 666 let print_tetra  =
 667   let (v1,v2,v3,v4) = (mk_tetrahedron 50) in
 668   let [w1;w2;w3;w4]= map (proj delta1 delta2) [v1;v2;v3;v4] in
 669   let draw ((ux,uy),(vx,vy)) = Printf.sprintf "\\draw (%f,%f)--(%f,%f);" ux uy vx vy in
 670   let drawv ((ux,uy),(vx,vy)) = Printf.sprintf "\\draw[very thick] (%f,%f)--(%f,%f);" ux uy vx vy in
 671   let face (s, (ux,uy),(vx,vy),(wx,wy)) =
 672     Printf.sprintf "\\draw[fill=black!%s] (%f,%f)--(%f,%f)--(%f,%f)--cycle;" s ux uy vx vy wx wy in
 673   let vertex (x,y) = Printf.sprintf "\\smalldot{%f,%f};" x y in
 674   let vertices = map vertex [w1;w2;w3;w4] in
 675   let vv = map draw [(w1,w2);(w1,w3);(w1,w4);(w2,w3);(w3,w4);(w2,w4)] in
 676   let shade = join_lines (map face [("45",w1,w2,w4); ("30",w1,w2,w3);("10",w2,w3,w4)]) in
 677   let edges=  join_lines vv in
 678   let triangulate (u,v,w) =
 679     let c = proj delta1 delta2 (0.3333 %... (u +... v +... w)) in
 680     let [muv;mvw;muw] = map (fun (u,v)-> proj delta1 delta2 (0.5 %... (u +... v)))
 681       [(u,v);(v,w);(u,w)] in
 682     let [pu;pv;pw] = map (proj delta1 delta2) [u;v;w] in
 683     let vv = map draw [(c,pu);(c,pv);(c,pw)] in
 684     let ww = map drawv [(c,muv);(c,mvw);(c,muw)] in
 685       join_lines (vv @ ww) in
 686     join_lines (edges :: shade :: (map triangulate [(v1,v2,v4);(v1,v2,v3);(v2,v3,v4)]) @ vertices);;
 687
 688 (* autoBWEYURN *)
 689
 690 let print_marchal seed=
 691   let radius = 0.15 in
 692   let radius_sqrt2 = 0.212 in
 693   let (_,sat) = randompacking radius seed 3.0 1.4 in (* 3,1.4 *)
 694   let satst = sat_triples sat in
 695   let rr = 0-- (List.length sat - 1) in
 696   let allpair = (outer rr rr) in
 697   let s = nth sat in
 698   let shortpair =filter(fun (i,j)-> (i<j)&& dist2 (s i) (s j) < 2.0 *. radius_sqrt2 -. 1.0e-4)
 699     allpair in
 700   let coord_triple t = center2 (nth sat) t in
 701   let prs = filter (fun ((i,j,k),(i',j',k')) ->
 702                       List.length (intersect [i;j;k] [i';j';k'])=2 && lex3 (i,j,k) (i',j',k'))
 703     (outer satst satst) in
 704   let pp = map (fun (t, t') ->
 705                   let (x,y) = coord_triple t in
 706                   let (x',y') = coord_triple t' in
 707                     Printf.sprintf "\\draw[very thick] (%f,%f) -- (%f,%f) ;" x y x' y') prs in
 708   let line (x,y)  = Printf.sprintf "\\draw[black,fill=black!20] (%f,%f) circle (%f);"
 709     x y radius_sqrt2  in
 710   let dot_line (x,y) = Printf.sprintf "\\smalldot{%f,%f};" x y  in
 711   let psmalldot = map (line) sat in
 712   let dot = map (dot_line) sat in
 713   let draw ((ux,uy),(vx,vy)) = Printf.sprintf "\\draw (%f,%f)--(%f,%f);" ux uy vx vy in
 714   let draw_radial (i,j,k) =
 715     let c = coord_triple (i,j,k) in
 716     let (u1,u2,u3)=(s i,s j,s k) in
 717     let vv = map draw  [(c,u1);(c,u2);(c,u3);] in
 718       join_lines vv in
 719   let radial = map draw_radial satst in
 720   let draw_dedge (t,t') =
 721     let c = coord_triple t in
 722     let c'= coord_triple t' in
 723     let (i,j) = tuple2 (intersect (list3 t) (list3 t')) in
 724     let u1 = s i in
 725     let u2 = s j in
 726     let w = u2 -.. u1 in
 727     let d = c -.. u1 in
 728     let d' = c' -.. u1 in
 729       if (det2 d w  *. det2 w d' > 0.0) then draw (u1, u2) else "%" in
 730   let dedge = map draw_dedge prs in
 731   let fill_tri (s,(ux,uy),(vx,vy),(wx,wy)) =
 732     Printf.sprintf "\\draw[black,fill=black!%s] (%f,%f)--(%f,%f)--(%f,%f)--cycle;"
 733       s ux uy vx vy wx wy in
 734   let rotate u v x =
 735     v +.. cmul (normalize2 (u -.. v)) x in
 736   let drawc (ax,ay) (bx,by) (cx,cy) (dx,dy) =
 737       Printf.sprintf "\\draw[fill=black!35] (%f,%f)--(%f,%f)--(%f,%f)--(%f,%f)--cycle;\n\\draw (%f,%f)--(%f,%f);"
 738         ax ay bx by cx cy dx dy ax ay cx cy in
 739   let draw_cell2 (i,j) =
 740     let u1 = s i in
 741     let u2 = s j in
 742     let r = dist2 u1 u2 in
 743     let h2 = radius_sqrt2 *. radius_sqrt2 -. r *. r /. 4.0 in
 744     let _ = (h2 >= 0.0) or failwith (Printf.sprintf "expected pos %d %d %f" i j h2) in
 745     let h = sqrt(h2) in
 746     let c = rotate u1 u2 (r /. 2.0,h) in
 747     let c' = rotate u1 u2 (r /. 2.0, -. h) in
 748       drawc u1 c u2 c' in
 749   let cell2 = map draw_cell2 shortpair in
 750   let draw_cell3 (i,j,k) =
 751     let (u,v,w) = (s i,s j,s k) in fill_tri ("60",u,v,w) in
 752   let cell3filter (i,j,k) =
 753     etaV (s i) (s j) (s k) < radius *. sqrt(2.0) in
 754   let cell3 = map draw_cell3 (filter cell3filter satst) in
 755     join_lines (psmalldot @ pp @ radial @ dedge @ cell2 @ cell3 @ dot);;
 756
 757 (* FIFJALK *)
 758
 759 let print_ferguson_hales seed=
 760   let radius = 0.15 in
 761   let radius_h = 1.255 *. radius  in
 762   let radius_h2 = 2.51 *. radius  in
 763   let radius_s = sqrt(8.0) *. radius in
 764   let (_,sat) = randompacking radius seed 3.0 1.4 in (* 3,1.4 *)
 765   let s = nth sat in
 766   let satst = sat_triples sat in
 767   let qrtet = filter (fun (i,j,k) -> dist2 (s i) (s j) <  radius_h2 &&
 768                      dist2 (s j) (s k) < radius_h2 &&
 769                      dist2 (s i) (s k) < radius_h2) satst in
 770   let quarter = filter (fun (i,j,k) ->
 771                           let [a;b;c] =
 772                             sort (<) [dist2 (s i) (s j);dist2 (s j) (s k);dist2 (s i) (s k)] in
 773                           (a < radius_h2 && b < radius_h2 && c < radius_s)) satst in
 774   let satst = sat_triples sat in
 775   let rr = 0-- (List.length sat - 1) in
 776   let allpair = (outer rr rr) in
 777   let coord_triple t = center2 (nth sat) t in
 778   let prs = filter (fun ((i,j,k),(i',j',k')) ->
 779                       List.length (intersect [i;j;k] [i';j';k'])=2 && lex3 (i,j,k) (i',j',k'))
 780     (outer satst satst) in
 781   let pp = map (fun (t, t') ->
 782                   let (x,y) = coord_triple t in
 783                   let (x',y') = coord_triple t' in
 784                     Printf.sprintf "\\draw[very thick] (%f,%f) -- (%f,%f) ;" x y x' y') prs in
 785   let shortpair =filter(fun (i,j)-> (i<j)&& dist2 (s i) (s j) < radius_h2 -. 1.0e-4)
 786     allpair in
 787   let circle (x,y)  = Printf.sprintf "\\draw[black,fill=black!20] (%f,%f) circle (%f);"
 788     x y radius_h  in
 789   let dot_line (x,y) = Printf.sprintf "\\smalldot{%f,%f};" x y  in
 790   let pcircle = map (circle) sat in
 791   let dot = map (dot_line) sat in
 792   let fill_tri (s,(ux,uy),(vx,vy),(wx,wy)) =
 793     Printf.sprintf "\\draw[black,fill=black!%s] (%f,%f)--(%f,%f)--(%f,%f)--cycle;"
 794       s ux uy vx vy wx wy in
 795   let rotate u v x =
 796     v +.. cmul (normalize2 (u -.. v)) x in
 797   let drawc (ax,ay) (bx,by) (cx,cy) (dx,dy) =
 798       Printf.sprintf "\\draw[fill=black!35] (%f,%f)--(%f,%f)--(%f,%f)--(%f,%f)--cycle;\n\\draw (%f,%f)--(%f,%f) (%f,%f)--(%f,%f);"
 799         ax ay bx by cx cy dx dy ax ay cx cy bx by dx dy in
 800   let draw_cell2 (i,j) =
 801     let u1 = s i in
 802     let u2 = s j in
 803     let r = dist2 u1 u2 in
 804     let h2 = radius_h *. radius_h -. r *. r /. 4.0 in
 805     let _ = (h2 >= 0.0) or failwith (Printf.sprintf "expected pos %d %d %f" i j h2) in
 806     let h = sqrt(h2) in
 807     let c = rotate u1 u2 (r /. 2.0,h) in
 808     let c' = rotate u1 u2 (r /. 2.0, -. h) in
 809       drawc u1 c u2 c' in
 810   let cell2 = map draw_cell2 shortpair in
 811   let draw_cell3 (i,j,k) =
 812     let (u,v,w) = (s i,s j,s k) in fill_tri ("60",u,v,w) in
 813   let cell3 = map draw_cell3 (qrtet @ quarter) in
 814     join_lines (pp @ pcircle  @ cell2 @ cell3 @ dot);;
 815
 816
 817
 818 (*
 819 SENQMWT
 820 *)
 821
 822 let print_thue seed=
 823   let radius = 0.15 in
 824   let radius_2_div_sqrt3 = radius *. 2.0 /. sqrt(3.0) in
 825   let (_,sat) = randompacking radius seed 3.0 1.4 in (* 3,1.4 *)
 826   let rr = 0-- (List.length sat - 1) in
 827   let allpair = (outer rr rr) in
 828   let s = nth sat in
 829   let shortpair =filter(fun (i,j)-> (i<j)&& dist2 (s i) (s j) < 2.0 *. radius_2_div_sqrt3 -. 1.0e-4)
 830     allpair in
 831   let line (x,y)  = Printf.sprintf "\\draw[black,fill=black!20] (%f,%f) circle (%f);"
 832     x y radius_2_div_sqrt3  in
 833   let dot_line (x,y) = Printf.sprintf "\\smalldot{%f,%f};" x y  in
 834   let psmalldot = map (line) sat in
 835   let dot = map (dot_line) sat in
 836   let rotate u v x =
 837     v +.. cmul (normalize2 (u -.. v)) x in
 838   let drawc (ax,ay) (bx,by) (cx,cy) (dx,dy) =
 839       Printf.sprintf "\\draw[fill=black!35] (%f,%f)--(%f,%f)--(%f,%f)--(%f,%f)--cycle;\n\\draw (%f,%f)--(%f,%f);"
 840         ax ay bx by cx cy dx dy bx by dx dy in
 841   let draw_cell2 (i,j) =
 842     let u1 = s i in
 843     let u2 = s j in
 844     let r = dist2 u1 u2 in
 845     let h2 = radius_2_div_sqrt3 *. radius_2_div_sqrt3 -. r *. r /. 4.0 in
 846     let _ = (h2 >= 0.0) or failwith (Printf.sprintf "expected pos %d %d %f" i j h2) in
 847     let h = sqrt(h2) in
 848     let c = rotate u1 u2 (r /. 2.0,h) in
 849     let c' = rotate u1 u2 (r /. 2.0, -. h) in
 850       drawc u1 c u2 c' in
 851   let cell2 = map draw_cell2 shortpair in
 852     join_lines (psmalldot  (* @  dedge *) @ cell2  @ dot);;
 853
 854
 855
 856 (* figKVIVUOT *)
 857
 858 let kv_inter r2 u v =
 859   let nu = norm3 u in
 860   let nvu = norm3 (v-...u) in
 861   let t = sqrt ( (r2 -. nu *. nu) /. (nvu *. nvu)) in
 862     u +... t %... (v-... u);;
 863
 864 (*
 865 let kv_interp s u v1 v2 =
 866   let r2 = 2.0 in
 867   let v = (s %... v2) +... ((1.0 -. s) %... v1) in
 868   let nu = norm3 u in
 869   let nvu = norm3 (v-...u) in
 870   let t = sqrt ( (r2 -. nu *. nu) /. (nvu *. nvu)) in
 871     u +... t %... (v-... u)
 872   ;;
 873 *)
 874
 875 let kv_interp r2 s u v1 v2 =
 876   let v = (s %... v2) +... ((1.0 -. s) %... v1) in
 877     kv_inter r2 u v;;
 878
 879
 880 let kv_proj rho (x,y,z) =  x %.. (1.0,0.0) +.. y %.. (0.0,1.0) +.. z %.. rho;;
 881
 882 let rx u1 u2 label =
 883   let r2 = 2.0 in
 884   let rho = (0.33,0.66) in
 885   let null3 = (0.0,0.0,0.0) in
 886   let p1 = map (fun s -> kv_interp r2 (float_of_int s /. 5.0) null3 u1 u2) (0--5) in
 887   let q1 = map (kv_proj rho) p1 in
 888   let w1 = join_space (map (fun (x,y)-> Printf.sprintf "(%f,%f) " x y) q1) in
 889     Printf.sprintf "\\def\\kv%s{%s}" label w1 ;;
 890
 891 let col1 =
 892   let a = 1.5 in
 893   let b = 2.0 in
 894   let c = b +. 0.2 in
 895   let bb = sqrt(b*.b -. a*.a) in
 896   let cc = sqrt(c*.c -. b*.b) in
 897   let r2 = 2.0 in
 898   let r = sqrt(r2) in
 899   let omega1 = a %... delta2 in
 900   let omega2 = omega1 +... bb %... delta1 in
 901   let omega3 = omega2 +... cc %... delta3 in
 902   let v1 = r %... normalize3 omega1 in
 903   let v2 = r %... normalize3 omega2 in
 904   let v3 = r %... normalize3 omega3 in
 905   let w12 = rx v1 v2 "oneab" in
 906   let w13 = rx v1 v3 "oneac" in
 907   let w32 = rx v3 v2 "onecb" in
 908     join_lines [w12;w13;w32];;
 909
 910 let col2 =
 911   let a = 1.0 in
 912   let b = 2.0 in
 913   let c = b +. 0.2 in
 914   let bb = sqrt(b*.b -. a*.a) in
 915   let cc = sqrt(c*.c -. b*.b) in
 916   let r2 = 2.0 in
 917   let r = sqrt(r2) in
 918   let omega1 = a %... delta2 in
 919   let omega2 = omega1 +... bb %... delta1 in
 920   let omega3 = omega2 +... cc %... delta3 in
 921   let v13 = kv_inter r2 omega1 omega3 in
 922   let v12 = kv_inter r2 omega1 omega2 in
 923   let v3 = r %... normalize3 omega3 in
 924   let v2 = r %... normalize3 omega2 in
 925   let wab = rx v13 v3 "Bab" in
 926   let wbc = rx v3 v2 "Bbc" in
 927   let wcd = rx v2 v12 "Bcd" in
 928   let wda = rx v12 v13 "Bda" in
 929     join_lines [wab;wbc;wcd;wda];;
 930
 931
 932 let col3 =
 933   let a = 1.0 in
 934   let b = 1.35 in
 935   let c = 2.0 in
 936   let bb = sqrt(b*.b -. a*.a) in
 937   let cc = sqrt(c*.c -. b*.b) in
 938   let r2 = 2.0 in
 939   let r = sqrt(r2) in
 940   let omega1 = a %... delta2 in
 941   let omega2 = omega1 +... bb %... delta1 in
 942   let omega3 = omega2 +... cc %... delta3 in
 943   let v13 = kv_inter r2 omega1 omega3 in
 944   let v23 = kv_inter r2 omega2 omega3 in
 945   let v3 = r %... normalize3 omega3 in
 946   let wab = rx v13 v3 "Cab" in
 947   let wbc = rx v3 v23 "Cbc" in
 948   let wca = rx v23 v13 "Cca" in
 949     join_lines [wab;wbc;wca];;
 950
 951
 952
 953 (*
 954 figDEJKNQK
 955
 956 v0 = {0, 0, 1};
 957 null = {0, 0, 0};
 958 v1 = {0, 0.85, Sqrt[1 - 0.85^2]};
 959 v2 = FarFrom[{1, 0, 0}, Vertex[{v0, 1}, {v1, 1}, {null, 1}]]
 960 v3 = FarFrom[{0, 0, -1}, Vertex[{null, 1}, {v1, 1}, {v2, 1.26}]]
 961 v4 = FarFrom[v1, Vertex[{null, 1}, {v2, 1}, {v3, 1.26}]]
 962 v5 = FarFrom[v2, Vertex[{null, 1}, {v4, 1}, {v3, 1.26}]]
 963 w4 = FarFrom[v1, Vertex[{null, 1}, {v2, 1}, {v3, 1}]]
 964 w5 = FarFrom[v2, Vertex[{null, 1}, {v4, 1}, {v3, 1}]]
 965 w6 = FarFrom[v4, Vertex[{null, 1}, {v3, 1}, {v5, 1}]]
 966 *)
 967
 968
 969 let vws =
 970   let v1=( 0.0, 0.85, 0.526783) in
 971   let v2=( -0.820069, 0.278363, 0.5) in
 972   let v3=( 0.32578, 0.00230249, 0.945443) in
 973   let  v4=( -0.586928, -0.690949, 0.422025) in
 974   let v5=( 0.329363, -0.938133, 0.106892) in
 975   let w4=( -0.409791, -0.617314, 0.671562) in
 976   let w5=(0.40333, -0.826891, 0.391887) in
 977   let  w6=(0.965333, -0.171655, 0.196637) in
 978                        (v1,v2,v3,v4,v5,w4,w5,w6);;
 979
 980 let rxdej(u1,u2,label) =
 981   let r2 = 1.0 in
 982   let rho = (0.0,0.0) in
 983   let null3 = (0.0,0.0,0.0) in
 984   let p1 = map (fun s -> kv_interp r2 (float_of_int s /. 5.0) null3 u1 u2) (0--5) in
 985   let q1 = map (kv_proj rho) p1 in
 986   let w1 = join_space (map (fun (x,y)-> Printf.sprintf "(%f,%f) " x y) q1) in
 987     Printf.sprintf "\\def\\dejk%s{%s}" label w1 ;;
 988
 989 let mkdejA =
 990   let (v1,v2,v3,v4,v5,w4,w5,w6) = vws in
 991   let vv = join_lines (map rxdej
 992       [(v1,v2,"a");(v2,v4,"b");(v4,v5,"c");(v5,w6,"d");(w6,v3,"e");(v3,v1,"f");
 993        (v2,v3,"g");(v4,v3,"h");(v5,v3,"i");(v4,w5,"j");(w5,v3,"k");(v2,w4,"l");
 994       (w4,v3,"m")]) in
 995     vv;;
 996
 997 let mkdejB =
 998   let rho = (0.0,0.0) in
 999   let (v1,v2,v3,v4,v5,w4,w5,w6) = vws in
1000   let a ((x,y),s) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s x y in
1001   let ww = join_lines (map (fun (v,s) -> a ((kv_proj rho v),s))
1002      [(v1,"v1");(v2,"v2");(v3,"v3");(v4,"v4");(v5,"v5");(w4,"w4");(w5,"w5");(w6,"w6")]) in
1003     ww;;
1004
1005
1006 (* QTICQYN
1007
1008 Mathematica:
1009 v0 = {0, 0, 1};
1010 null = {0, 0, 0};
1011 v1 = {0, 0.85, Sqrt[1 - 0.85^2]};
1012 v2 = FarFrom[{1, 0, 0}, Vertex[{v0, 1}, {v1, 1}, {null, 1}]];
1013 v3 = FarFrom[{0, 0, -1}, Vertex[{null, 1}, {v1, 1}, {v2, 1.5}]];
1014 v4 = FarFrom[v1, Vertex[{null, 1}, {v2, 1}, {v3, 1.5}]];
1015 v5 = FarFrom[v2, Vertex[{null, 1}, {v4, 1}, {v3, 1}]];
1016 *)
1017
1018 let qtvv =
1019   let v1=(0.0, 0.85, 0.526783) in
1020   let v2=(-0.820069, 0.278363,   0.5) in
1021   let v3=(0.651104, 0.124197,       0.748758) in
1022   let v4=(-0.572254, -0.688878, 0.444941) in
1023   let v5=(0.421862, -0.796553, 0.433055) in
1024     (v1,v2,v3,v4,v5);;
1025
1026
1027 let rxqt(u1,u2,label) =
1028   let r2 = 1.0 in
1029   let rho = (0.0,0.0) in
1030   let null3 = (0.0,0.0,0.0) in
1031   let p1 = map (fun s -> kv_interp r2 (float_of_int s /. 5.0) null3 u1 u2) (0--5) in
1032   let q1 = map (kv_proj rho) p1 in
1033   let w1 = join_space (map (fun (x,y)-> Printf.sprintf "(%f,%f) " x y) q1) in
1034     Printf.sprintf "\\def\\qt%s{%s}" label w1 ;;
1035
1036 let mkqtA =
1037   let (v1,v2,v3,v4,v5) = qtvv in
1038   let vv = join_lines (map rxqt
1039       [(v1,v2,"a");(v2,v4,"b");(v4,v5,"c");(v5,v3,"d");(v3,v1,"e");(v2,v3,"f");
1040        (v4,v3,"g")]) in
1041     vv;;
1042
1043 let mkqtB =
1044   let rho = (0.0,0.0) in
1045   let (v1,v2,v3,v4,v5) = qtvv in
1046   let a ((x,y),s) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s x y in
1047   let ww = join_lines (map (fun (v,s) -> a ((kv_proj rho v),s))
1048      [(v1,"v1");(v2,"v2");(v3,"v3");(v4,"v4");(v5,"v5")]) in
1049     ww;;
1050
1051 (*
1052 HEABLRG
1053 *)
1054
1055 let vvhe =
1056   let (a,b,c) = (0.45,0.6,0.4) in
1057     map normalize3 [(a,b,c);(-. a,b,c);(-.a,-.b,c);(a,-.b,c)];;
1058
1059 let dualhe =
1060   let [v1;v2;v3;v4] = vvhe in
1061   let nc (u,v) = normalize3 (cross u v) in
1062    map nc [(v1,v2);(v2,v3);(v3,v4);(v4,v1)];;
1063
1064 let mkhe =
1065   let rho = (0.0,0.0) in
1066   let [v1;v2;v3;v4] = vvhe in
1067   let [w1;w2;w3;w4] = dualhe in
1068   let vv = join_lines (map rxqt
1069       [(v1,v2,"a");(v2,v3,"b");(v3,v4,"c");(v4,v1,"d");(w1,w2,"e");(w2,w3,"f");
1070        (w3,w4,"g");(w4,w1,"h")]) in
1071   let a ((x,y),s) = Printf.sprintf "\\coordinate (%s) at (%f,%f);" s x y in
1072   let ww = join_lines (map (fun (v,s) -> a ((kv_proj rho v),s))
1073      [(v1,"v1");(v2,"v2");(v3,"v3");(v4,"v4");(w1,"w1");(w2,"w2");(w3,"w3");
1074       (w4,"w4")]) in
1075     join_lines [ww;vv];;
1076
1077
1078
1079 (* figYAHDBVO *)
1080
1081 let vv i = (72.0*.i +. (-.40.0) +. Random.float 80.0,Random.float 1.5 +. 1.0);;
1082
1083 (* map (vv o float_of_int) [0;1;2;3;4];; *)
1084
1085 let vout =[ (-1.40573615605411817, 2.43152527011496122);
1086    (62.2310998421659392, 1.50101500229540341);
1087    (166.419445012932584, 1.80579527399678152);
1088    (206.27916875919712, 1.73501080990908485);
1089    (293.766402309343221, 1.59228179599956721)];;
1090
1091 let midangles =
1092   let mm = map fst vout in
1093   let suc i = ((i+1) mod 5) in
1094   let mm1 i = nth mm (suc i) +. if (nth mm i < nth mm (suc i)) then 0.0 else 360. in
1095   let mm' i = 0.5 *. (nth  mm i  +. mm1 i ) in
1096   let f =   map mm' (0--4) in
1097   let mm' i = 1.0/. cos( (pi /. 180.0) *.(0.5 *. (mm1 i -. nth mm i))) in
1098   let g = map mm' (0--4) in
1099    zip f g;;
1100
1101 let poly2_extreme =
1102   let m = map (fun (theta,r)-> (r *. cos (radian theta), r*. sin (radian theta))) vout in
1103   let v i = nth m i in
1104   let suc i = v ((i+1) mod 5) in
1105   let inter i = solve22 (v i, suc i) (v i *.. v i, suc i *.. suc i) in
1106     map inter (0--4);;
1107
1108 (* figZXEVDCA *)
1109
1110
1111 let fix_SO3 =  (*  random_SO3 () ;;  *)
1112   ((-0.623709278332764572, -0.768294961660169751, -0.143908262477248777),
1113    (-0.592350038826666592, 0.34445174255424732, 0.728336754910384077),
1114    (-0.510008007411323683, 0.539514456654259789, -0.669937298138706616));;
1115
1116
1117 (* vertices of an icosahedron, vector lengths sqrt3. *)
1118
1119 let icos_vertex =
1120   let sqrt3 = sqrt(3.0) in
1121   let v = sqrt3 %... (1.0,0.0,0.0) in
1122 (*  let d0 = 2.10292 in *)  (* 20 Solid[2,2,2,d0,d0,d0] = 4 Pi *)
1123   let theta = 1.10715 in (* arc[2,2,d0] = theta *)
1124   let ct = cos theta in
1125   let st = sin theta in
1126   let p5 = pi/. 5.0 in
1127   let vv =    map (mul3 fix_SO3)
1128     ( v :: map
1129        (fun i->  (sqrt3 %... (ct, st *. cos (i *. p5), st *. sin (i *. p5))))
1130        [2.0;4.0;6.0;8.0;10.0]) in
1131    (vv @ (map ( uminus3 ) vv));;
1132
1133 let iv  = nth icos_vertex;;
1134
1135 let icos_edge =
1136   let v i = nth icos_vertex i in
1137   let dall = filter (fun (i,j) -> (i<j)) (outer (0--11) (0--11)) in
1138    filter (fun (i,j) -> dist3 (v i) (v j) < 2.2) dall;;  (* note 2.2 cutoff *)
1139
1140 let icos_face =
1141   let micos (i,j) = mem ( i, j ) icos_edge in
1142   let balance (i,(j,k)) = (i,j,k) in
1143     map balance (filter (fun (i,(j,k)) -> micos (i,j) && micos (i,k) && micos (j,k))
1144      (outer (0--11) (outer (0--11) (0--11))));;
1145
1146 let dodec_vertex =
1147   map (fun (i,j,k) -> extreme_point (iv i,iv j,iv k)) icos_face;;  (* voronoi cell vertices *)
1148
1149 let next_icos_face (a,b,u3)=  (* input flag: [a] subset u2 subset u3 *)
1150   let v3 = list3 u3  in
1151   let _ = mem a v3 or failwith "next_dodec_face a" in
1152   let _ = mem b v3 or failwith "next_dodec_face b" in
1153   let ifc = map (tuple3) (filter (subset [a;b]) (map list3 icos_face)) in
1154   let ifc' = subtract ifc [u3] in
1155   let _ = List.length ifc' = 1 or failwith "next_dodec_face c" in
1156   let w3 = nth ifc' 0 in
1157   let cx = subtract (list3 w3) [a;b] in
1158   let _ = List.length cx = 1 or failwith "next_dodec_face d" in
1159   let c = hd cx in
1160   let w2x = filter (fun (i,j)->(i<j)) [(a,c);(c,a)] in
1161   let _ = List.length w2x = 1 or failwith "next_dodec_face e" in
1162     (a,c,w3);;
1163
1164 let icos_vertex_cycle a =
1165   let [i;j;k] = hd (filter (mem a) (map list3 icos_face)) in
1166   let startp = (a,(if (a=i) then j else i),(i,j,k)) in
1167   let t = map (fun i -> funpow i next_icos_face startp) [0;1;2;3;4] in
1168     map (fun (_,_,u) -> u) t;;
1169
1170 let icos_cycles = map icos_vertex_cycle (0--11);;
1171
1172 let ht xxs =
1173   let (_,_,z) = end_itlist ( +... ) xxs in
1174     z;;
1175
1176 let sort_dodec_face =
1177   let t = icos_cycles in
1178   let lookup = zip icos_face dodec_vertex in
1179   let htc cycle = ht (map (fun i -> assoc i lookup) cycle) in
1180   let hts = map htc t in
1181   let z = zip hts t in
1182   let z' = filter (fun (h,_) -> (h>0.0)) z in
1183   let t' = map snd (sort (fun (a,_) (b,_) -> a<b) z') in
1184     t';;
1185
1186 let center_face cycle =
1187   let lookup = zip icos_face dodec_vertex in
1188   let coords = map (fun i -> assoc i lookup) cycle in
1189     (1.0 /. float_of_int (List.length cycle)) %... (end_itlist (+...) coords);;
1190
1191 let pname (i,j,k) = Printf.sprintf "V%d-%d-%d" i j k;;
1192
1193 let print_cycles =
1194   let lookup = zip icos_face (map (proj delta1 delta2) dodec_vertex) in
1195     map (fun (r,(x,y)) -> Printf.sprintf "\\coordinate (%s) at (%f,%f);" (pname r) x y) lookup;;
1196
1197 let print_dodec_face =
1198   let opt = "fill=white" in
1199   let pdraw = Printf.sprintf "\\draw[%s] " opt in
1200   let cycle m = join_space (map (fun s -> Printf.sprintf "(%s)--" ( pname s)) m) in
1201   let s m = pdraw ^ (cycle m) ^ "cycle;" in
1202     map s sort_dodec_face;;
1203
1204 (* let print_dodec = join_lines (print_cycles @ print_dodec_face);; *)
1205
1206
1207
1208 (* printing spherical caps.
1209    Parameters:
1210      R = radius of sphere at the origin
1211      v =(_,_,_) norm 1 vector pointing to center of cap.
1212      theta = arclength on unit sphere of cap.
1213
1214 *)
1215
1216
1217
1218 let frame_cap v =
1219   let v = normalize3 v in
1220   let (x,y,z) = v in
1221   let w = normalize3 (-. y,x,0.0) in
1222   frame_of v w;;
1223
1224 let ellipse_param v rad theta =
1225   let (v,w,u) = frame_cap v in
1226   let q = (rad *. cos theta) %... v in
1227   let qbar = proj delta1 delta2 q in
1228   let p = ((rad *. cos theta) %... v) +... ((rad *. sin theta) %... u) in
1229   let pbar = proj delta1 delta2 p in
1230   let h = dist2 qbar pbar in
1231   let k = rad *. sin theta in
1232     (qbar,h,k);;
1233
1234 let calc_psi theta v =
1235   let (x,y,_) = normalize3 v in
1236   let cospsi = cos theta /. sqrt( x *. x +. y *. y) in
1237     if (abs_float cospsi <= 1.0) then acos cospsi else 0.0;;
1238
1239 let calc_alpha rad psi qbar =
1240   let nqbar = normalize2 qbar in
1241   let rtrue = cmul (rad *. cos psi, rad *. sin psi) nqbar in
1242   let a = normalize2 (rtrue -.. qbar) in
1243     acos (a *.. nqbar);;
1244
1245 let adjust_alpha h k alpha =  (* compensate for TIKZ BUG in arc specs *)
1246   let ca = cos alpha in
1247   let sa = sin alpha in
1248   let t = 1.0 /. sqrt(ca *. ca /. (h *. h) +. sa *. sa /. (k *. k)) in
1249     acos (t *. ca /. h);;
1250
1251 let print_ellipse rad qbar h k psi =
1252   let nqbar = normalize2 qbar in
1253   let r = (rad *. cos psi, rad *. sin psi) in
1254   let (qbx,qby) = qbar in
1255   let (px,py) = cmul r nqbar in
1256   let (px',py') = cmul (conj r) nqbar in
1257   let alpha = adjust_alpha h k (calc_alpha rad psi qbar) in
1258   let endangle = 2.0 *. pi -. alpha in
1259   let rotateAngle = degree (arg qbx qby) in
1260   let cstart = degree (arg px' py') in
1261   let delta = 2.0 *. psi in
1262   let s = "\\draw[ball color=gray!10,shading=ball] (0,0) circle (1);\n\\end{scope}" in
1263     if (psi<0.01) then
1264        Printf.sprintf
1265         "\\begin{scope} \\clip (%f,%f) circle[x radius=%f,y radius=%f,rotate=%f];\n%s"
1266         qbx qby h k rotateAngle s
1267     else
1268       Printf.sprintf
1269         "\\begin{scope}\\clip (%f,%f) arc[x radius=%f,y radius = %f,start angle=%f,end angle=%f,rotate=%f] arc[radius=1.0,start angle=%f,delta angle=%f];\\pgfpathclose;\n%s"
1270         px py h k  (degree alpha) (degree endangle) rotateAngle
1271         cstart (degree delta)
1272         s;;
1273
1274 let print_dodec_ellipse =
1275   let vs = map center_face sort_dodec_face in
1276   let theta = pi /. 6.0 in
1277   let rad = 1.0 in
1278   let one_ellipse v =
1279     let (q,h,k) = ellipse_param v rad theta in
1280     let psi = calc_psi theta v in
1281      print_ellipse rad q h k psi in
1282     map one_ellipse vs;;
1283
1284
1285
1286
1287
1288 (* output *)
1289
1290 let outfilestring = "/tmp/x.txt";;
1291
1292 let execute() =
1293   let outs = open_out outfilestring in
1294   let write_out s = try (Printf.fprintf outs "%s" s)
1295   with _ as t -> (close_out outs; raise t) in
1296   let _ = write_out  "% AUTO GENERATED BY tikz.ml. Do not edit.\n" in
1297   let wrap s s' = Printf.sprintf "\\def\\%s{%s}\n\n\n" s s' in
1298   let add name s = write_out (wrap name s) in
1299   let _ =  add "autoZXEVDCA"
1300     (join_lines (print_cycles @ print_dodec_face @ print_dodec_ellipse)) in
1301   let voronoi_seed = 12345678 in (* 50, 300 good, *)
1302   let _ = add "autoDEQCVQL" (print_satunsat 5) in
1303 (* ok to here. In ML mode next line causes stack overflow. *)
1304   let _ = add "autoXOHAZWO" (print_satst voronoi_seed) in
1305   let _ = add "autoBUGZBTW" (print_rogers voronoi_seed) in
1306   let _ = add "autoORQISJR" (print_rogers 45) in
1307   let _ =  add "autoSENQMWT" (print_thue 45) in
1308   let _ =  add "autoEVIAIQP" (print_voronoi 45) in
1309   let _ =  add "autoANNTKZP" (print_delaunay 45) in
1310   let _ =  add "autoFIFJALK" (print_ferguson_hales 45) in
1311   let _ =  add "autoBWEYURN" (print_marchal 45) in
1312   let _ =  add "autoYAJOTSL" (print_tetra) in
1313   let _ =  add "autoKVIVUOT" (join_lines[col1;col2;col3]) in
1314   let _ =  add "autoDEJKNQK" (mkdejA) in
1315   let _ =  add "autocDEJKNQK" (mkdejB) in
1316   let _ =  add "autoQTICQYN" (join_lines [mkqtA;mkqtB]) in
1317   let _ =  add "autoHEABLRG" (mkhe) in
1318   let _ =  add "autoSEYIMIE" (fcc_fun_domain) in
1319   let _ =  add "autoAZGXQWC" (tet_oct_ratio) in
1320   let _ =  add "autoTCFVGTS" (cubic_layers) in
1321   let _ =  add "autoPQJIJGE" (rhombic_dodec) in
1322   let _ =  add "autoSGIWBEN" (fcc_hcp_pattern) in
1323   let _ =  add "autoDHQRILO" (fcc_packing) in
1324   let _ =  add "autoBDCABIA" (pascal_packing) in
1325   let _ =  add "autoNTNKMGO" (square_layers) in
1326   let _ = close_out outs in
1327   "to save results move /tmp/x.txt to kepler_tex/x.txt"
1328     ;;
1329
1330 (*
1331 let gen2_out =
1332   let wrap s s' = Printf.sprintf "\\def\\%s{%s}\n\n\n" s s' in
1333   let outstring = ref "" in
1334   let add name s = outstring:= !outstring ^ wrap name s in
1335   let _ =  add "autoYAJOTSL" (print_tetra) in
1336   let _ =  add "autoKVIVUOT" (join_lines[col1;col2;col3]) in
1337   let _ =  add "autoDEJKNQK" (mkdejA) in
1338   let _ =  add "autocDEJKNQK" (mkdejB) in
1339   let _ =  add "autoQTICQYN" (join_lines [mkqtA;mkqtB]) in
1340   let _ =  add "autoHEABLRG" (mkhe) in
1341   let _ =  add "autoSEYIMIE" (fcc_fun_domain) in
1342   let _ =  add "autoAZGXQWC" (tet_oct_ratio) in
1343   let _ =  add "autoTCFVGTS" (cubic_layers) in
1344   let _ =  add "autoPQJIJGE" (rhombic_dodec) in
1345   let _ =  add "autoSGIWBEN" (fcc_hcp_pattern) in
1346   let _ =  add "autoDHQRILO" (fcc_packing) in
1347   let _ =  add "autoBDCABIA" (pascal_packing) in
1348   let _ =  add "autoNTNKMGO" (square_layers) in
1349     output_filestring "/tmp/y.txt" (!outstring);;
1350
1351 let genz_out  =
1352   let wrap s s' = Printf.sprintf "\\def\\%s{%s}\n\n\n" s s' in
1353   let outstring = ref "" in
1354   let add name s = outstring:= !outstring ^ wrap name s in
1355   let _ =  add "autoNTNKMGO" (square_layers) in
1356     output_filestring "/tmp/z.txt" (!outstring);;
1357 *)
1358
1359 end;;