diff --git a/blizzard.cpp b/blizzard.cpp
new file mode 100644
index 00000000..12f677b7
--- /dev/null
+++ b/blizzard.cpp
@@ -0,0 +1,250 @@
+double randd() { return (rand() % 1000000) / 1000000. + .0000005; }
+
+hyperpoint randomPointIn(int t) {
+  while(true) {
+    hyperpoint h = spin(2*M_PI*(randd()-.5)/t) * tC0(xpush(asinh(randd())));
+    double d =
+      purehepta ? tessf : t == 6 ? hexhexdist : crossf;
+    if(hdist0(h) < hdist0(xpush(-d) * h))
+      return spin(2*M_PI/t * (rand() % t)) * h;
+    }
+  }
+
+struct snowball {
+  transmatrix T;
+  transmatrix global;
+  snowball *prev;
+  snowball *next;
+  double phase;
+  snowball(int t) { T = rgpushxto0(randomPointIn(t)); phase = randd(); }
+  };                                                                
+
+struct blizzardcell {
+  cell *c;
+  int frame;
+  int tmp;
+  transmatrix *gm;
+  char wmap;
+  int inward, outward, ward;
+  int qty[7];
+  vector<snowball*> inorder, outorder;
+  int inid, outid;
+  ~blizzardcell() { for(auto i: inorder) delete i; }
+  };
+
+map<cell*, blizzardcell> blizzardcells;
+
+vector<blizzardcell*> bcells;
+
+int N;
+
+blizzardcell* getbcell(cell *c) {
+  int i = c->aitmp;
+  if(i<0 || i >= N) return NULL;
+  if(bcells[i]->c != c) return NULL;
+  return bcells[i];
+  }
+
+void drawBlizzards() {
+  auto it = blizzardcells.begin();
+  bcells.clear();
+  while(it != blizzardcells.end()) 
+    if(it->second.frame != frameid || !gmatrix.count(it->first))
+      it = blizzardcells.erase(it);
+      else {
+        it->second.c = it->first;
+        bcells.push_back(&it->second);
+        it++;
+        }
+  N = size(bcells);
+  for(int i=0; i<N; i++) {
+    auto& bc = *bcells[i];
+    bc.tmp = bc.c->aitmp,
+    bc.c->aitmp = i;
+    bc.gm = &gmatrix[bc.c];
+    bc.wmap = windmap::at(bc.c);
+    }
+
+  for(int i=0; i<N; i++) {
+    auto& bc = *bcells[i];
+    cell *c = bc.c;
+    bc.inward = bc.outward = 0;
+    for(int i=0; i<c->type; i++) {
+      int& qty = bc.qty[i];
+      qty = 0;
+      cell *c2 = c->mov[i];
+      if(!c2) continue;
+      auto bc2 = getbcell(c2);
+      if(!bc2) continue;
+      int z = (bc2->wmap - bc.wmap) & 255;
+      if(z >= windmap::NOWINDBELOW && z < windmap::NOWINDFROM)
+        bc.outward += qty = z / 8;
+      z = (-z) & 255;
+      if(z >= windmap::NOWINDBELOW && z < windmap::NOWINDFROM)
+        bc.inward += z / 8, qty = -z/8;
+      }
+    bc.ward = max(bc.inward, bc.outward);
+    while(size(bc.inorder) < bc.ward) {
+      auto sb = new snowball(c->type);
+      bc.inorder.push_back(sb);
+      bc.outorder.push_back(sb);
+      }
+    for(auto& sb: bc.inorder) sb->prev = sb->next = NULL;
+    bc.inid = 0;
+    }
+  
+  double at = (ticks % 250) / 250.0;
+
+  for(int i=0; i<N; i++) {
+    auto& bc = *bcells[i];
+    cell *c = bc.c;
+    
+    for(auto sb: bc.inorder)
+      sb->global = (*bc.gm) * sb->T;
+    
+    bc.outid = 0;
+    
+    for(int d=0; d<c->type; d++) for(int k=0; k<bc.qty[d]; k++) {
+      auto& bc2 = *getbcell(c->mov[d]);
+      auto& sball = *bc.outorder[bc.outid++];
+      auto& sball2 = *bc2.inorder[bc2.inid++];
+      sball.next = &sball2;
+      sball2.prev = &sball;
+      
+      hyperpoint t = inverse(sball.global) * tC0(sball2.global);
+      double at0 = at + sball.phase;
+      if(at0>1) at0 -= 1;
+      transmatrix tpartial = sball.global * rspintox(t) * xpush(hdist0(t) * at0);
+      
+      queuepoly(tpartial, shSnowball, 0xFFFFFF80);
+      }
+    }
+
+  for(int ii=0; ii<N; ii++) {
+    auto& bc = *bcells[ii];
+    int B = size(bc.outorder);
+    if(B<2) continue;
+    int i = rand() % B;
+    int j = rand() % (B-1);
+    if(i==j) j++;
+    
+    if(1) {
+      auto& sb1 = *bc.outorder[i];
+      auto& sb2 = *bc.outorder[j];
+      
+      double swapcost = 0;
+      if(sb1.next) swapcost -= hdist(tC0(sb1.global), tC0(sb1.next->global));
+      if(sb2.next) swapcost -= hdist(tC0(sb2.global), tC0(sb2.next->global));
+      if(sb1.next) swapcost += hdist(tC0(sb2.global), tC0(sb1.next->global));
+      if(sb2.next) swapcost += hdist(tC0(sb1.global), tC0(sb2.next->global));
+      if(swapcost < 0) {
+        swap(bc.outorder[i], bc.outorder[j]);
+        swap(sb1.next, sb2.next);
+        if(sb1.next) sb1.next->prev = &sb1;
+        if(sb2.next) sb2.next->prev = &sb2;
+        }
+      }
+    
+    if(1) {
+      auto& sb1 = *bc.inorder[i];
+      auto& sb2 = *bc.inorder[j];
+      
+      double swapcost = 0;
+      if(sb1.prev) swapcost -= hdist(tC0(sb1.global), tC0(sb1.prev->global));
+      if(sb2.prev) swapcost -= hdist(tC0(sb2.global), tC0(sb2.prev->global));
+      if(sb1.prev) swapcost += hdist(tC0(sb2.global), tC0(sb1.prev->global));
+      if(sb2.prev) swapcost += hdist(tC0(sb1.global), tC0(sb2.prev->global));
+      if(swapcost < 0) {
+        swap(bc.inorder[i], bc.inorder[j]);
+        swap(sb1.prev, sb2.prev);
+        if(sb1.prev) sb1.prev->next = &sb1;
+        if(sb2.prev) sb2.prev->next = &sb2;
+        }
+      }
+    
+    auto& sbp = *bc.inorder[i];
+    if(sbp.next && sbp.prev) {
+      double p1 = sbp.next->phase;
+      double p2 = sbp.prev->phase;
+      double d = p2-p1;
+      if(d<=.5) d+=1;
+      if(d>=.5) d-=1;
+      sbp.phase = p1 + d/2;
+      if(sbp.phase >= 1) sbp.phase -= 1;
+      if(sbp.phase < 0) sbp.phase += 1;
+      }
+    }
+
+  for(auto bc: bcells)
+    bc->c->aitmp = bc->tmp;
+  }
+       
+#if 0
+    {
+      int v0 = windmap::at(c);
+      const bool method2 = true;
+      
+      auto a = [] (cell *c) {
+        ld dx=0, dy=0;
+        int v0 = windmap::at(c);
+        forCellIdEx(c2, i, c) {
+          int v1 = windmap::at(c->mov[i]);
+          if(!c2) continue;
+          int xv1 = (v1 - v0) & 255;
+          if(xv1 >= 128) xv1 -= 256;
+          int hdir = displaydir(c, i);
+          dx += sin(hdir * M_PI / 42) * xv1;
+          dy += cos(hdir * M_PI / 42) * xv1;
+          }
+        return atan2(dx, dy);
+        };
+      
+      ld alpha = a(c);
+      
+      double xx = M_PI / (1 + sqrt(5));
+      
+      forCellIdEx(c2, i, c) {
+        int v1 = windmap::at(c2);
+        int xv1 = (v1 - v0) & 255;
+        if(xv1 >= windmap::NOWINDBELOW && xv1 < windmap::NOWINDFROM) {
+        
+          // method2 ? (ticks/16) - v0 : ticks/4;
+          
+          double xcph = (ticks/16.) - v0;
+          xcph -= floor(xcph / 256) * 256;
+          
+          ld beta = a(c2);
+          for(int u=0; u<(method2 ? 16 : 1); u++) {
+            xcph += 16; if(xcph >= 256) xcph -= 256;
+            if(method2 ? xcph < xv1 : true) {
+              int hdir = displaydir(c, i);
+              int aircol = 0xFFFFFF00; // 0xD0D0FFFF;
+              aircol += int(255. * windmap::NOWINDBELOW / xv1); //  (1-exp(-xv1 / 50.)));
+              double at = method2 ? xcph*1./xv1 : xcph / 256.;
+  
+              double ldist = 
+                purehepta ? tessf :
+                c->type == 6 && c2->type == 6 ? hexhexdist :
+                crossf;
+              
+/*              transmatrix alphamatrix = (*Vdp) * spin(M_PI/2) * xpush(sin(u*xx) * ldist/4) * spin(-M_PI/2) * spin(alpha) * xpush(ldist*at);
+              transmatrix betamatrix = (*Vdp) * spin(hdir*M_PI/42) * xpush(ldist) * 
+                spin(-displaydir(c2, c->spn(i))*M_PI/S42) * 
+                spin(M_PI/2) * xpush(-sin(u*xx)*ldist/4) * spin(-M_PI/2) * 
+                spin(beta) * xpush(ldist * (1-at)); 
+              
+              hyperpoint t = inverse(alphamatrix) * tC0(betamatrix);
+              transmatrix tpartial = alphamatrix * rspintox(t) * xpush(hdist0(t) * at * at * (3-2*at)); */
+                
+              poly_outline = OUTLINE_TRANS;
+              transmatrix V0 = spin((hdir) * M_PI / S42);
+              // queuepoly(tpartial, shSnowball, aircol);
+              queuepoly((*Vdp) * V0 * xpush(at*ldist), shSnowball, aircol);
+              poly_outline = OUTLINE_DEFAULT;
+              }          
+            }
+
+          }         
+        }
+      }
+#endif