hyperrogue/blizzard.cpp

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