src/FDM/YASim/Turbulence.cpp

   1 #include "Turbulence.hpp"
   2 #include "Math.hpp"
   3
   4 namespace yasim {
   5 // Typical velocity spectrum: MIN 0.017 MAX 0.72 AVG 0.30 RMS 0.33
   6
   7 // Maximum conceivable turbulence flow, in m/s.  In practice, most
   8 // generated turbulence fields top out at about 70% of this number.
   9 const float MAX_TURBULENCE = 20;
  10
  11 // Rate, in "meters" per second, of the default time axis motion.
  12 // This will be multiplied by the rate-hz property to get the actual
  13 // time axis offset.  A setting of 2.0 causes the maximum frequency
  14 // component to arrive at 1 Hz.
  15 const float BASE_RATE = 2.0;
  16
  17 // Power to which the input magnitude (always in the range [0:1]) is
  18 // raised to get a coefficient for the turbulence velocity.  Setting
  19 // this to 1.0 makes the scale linear.  Increasing it makes it
  20 // curvier, with a sharp increase at the high end of the scale.
  21 const double MAGNITUDE_EXP = 2.0;
  22
  23 // How many generations are "meaningful" (i.e., not part of the normal
  24 // wind computation).  Decreasing this number will reallocate
  25 // bandwidth to the higher frequency components.  A turbulence field
  26 // will swing between maximal values over a distance of approximately
  27 // 2^(MEANINGFUL_GENS-1).
  28 const int MEANINGFUL_GENS = 7;
  29
  30 static const float FT2M = 0.3048;
  31
  32 // 8 x 32 s-box used by hashrand.  Read out of /dev/random on my Linux
  33 // box; not analyzed for linearity or coverage issues.
  34 static unsigned int SBOX[] = {
  35     0x0a881716, 0x20daa8ee, 0x61eb7d78, 0x46164e74, 0x39ab9d9d, 0x633a33f6,
  36     0x437c821d, 0x60a66f29, 0xc4ae45ab, 0x9a5cb3ce, 0x4a43606a, 0x56802c3c,
  37     0xe40d5e25, 0xa0297f41, 0x0457671e, 0xf167ab77, 0x571276db, 0x8b644e02,
  38     0xd5cfc592, 0x2331bfa2, 0xf9dfe7c1, 0xce9e7583, 0xfb133c29, 0x951c31c9,
  39     0x8e61b24e, 0xddf37570, 0x938c3b72, 0xaf907468, 0x98b77ac7, 0xe6edd515,
  40     0xa01f3600, 0xeafea5ad, 0x83fcce08, 0xe2e9fa9d, 0xd87727bb, 0x1945ea4c,
  41     0x831d295f, 0xa796ed85, 0xaa907b24, 0x69b25f12, 0xd4b27868, 0xdcde40f5,
  42     0x0e9e6def, 0x348a4702, 0x298389c8, 0xce405b63, 0x2e36d5a3, 0xf0569882,
  43     0x3beb3219, 0xf2366b9a, 0x69576cca, 0xd2725b8b, 0x6016d6f3, 0x728142ca,
  44     0x448b9f47, 0xe600cd4e, 0xac45d524, 0x0e32531b, 0x425d7b55, 0xc65cd9dc,
  45     0x58d7f9f1, 0x19f49822, 0x6786f2d3, 0x57844748, 0x523de4a3, 0x01079655,
  46     0x6dccea89, 0xb59278f2, 0x13a27e83, 0x19bcfc69, 0x4cff4bf5, 0xb18a3441,
  47     0x1e235c5e, 0xa1b47a42, 0x3bee8a5a, 0xa0962594, 0xa9b1ce4c, 0xb00399c8,
  48     0x83749847, 0x42c666e7, 0x08b81e57, 0xf7eee24b, 0x66720817, 0x3983f5f8,
  49     0x4999a817, 0x94fabd7a, 0x7aa775ef, 0xf6c1adcb, 0x5f32a695, 0x813ecf7e,
  50     0x66615fd5, 0xc0012e15, 0x051dd97e, 0xe6ee2803, 0x2449663c, 0x4024d59c,
  51     0xcb70a774, 0xacd3db94, 0x1845484e, 0xc803ef3c, 0x0662876f, 0x8794fe30,
  52     0xf0f0d16a, 0x41c065b8, 0xff9d5fc7, 0xa4237394, 0x8656614d, 0x26be5da9,
  53     0xb32bc625, 0xf215cc58, 0xc1e21848, 0xb97fe9fc, 0xbb28ef04, 0xde88eb23,
  54     0xe0623033, 0xa3df9e9c, 0xe9b95887, 0x3a4ab03b, 0x1cba812e, 0x174b4b37,
  55     0x0074d24b, 0xe5668d09, 0xf11a070a, 0x2884252b, 0x911149ea, 0x20dab459,
  56     0x89573d33, 0x68c2711d, 0x2b8e9781, 0xf007567b, 0x9761c8fa, 0x574d3a4e,
  57     0xa2ac28dd, 0x924f2211, 0xb0a91028, 0x83a90487, 0xf22cf6f8, 0x17a5dcfe,
  58     0x10497534, 0x27dd1316, 0x94a34815, 0x276e11ee, 0xead1d779, 0x0bfd4f20,
  59     0x45f2228f, 0x35d21bf8, 0x121336c0, 0x43a6538b, 0x55e950dc, 0x88a80871,
  60     0xfda9f61e, 0x5c76d120, 0x2eb8338f, 0x5193bb8e, 0x30a6995c, 0x500505a8,
  61     0x7b214f6a, 0x6a74558d, 0x040d0716, 0x4452846b, 0xd0a0e838, 0xead282e0,
  62     0x6bc6465c, 0xcb4ab107, 0xab990ed7, 0x72a1fe7b, 0x06901fdf, 0x18f90739,
  63     0x8cd2b861, 0xaea9d40c, 0x2dcf7c18, 0x45979e8a, 0x10393f0d, 0x3209d7c9,
  64     0x2c71378f, 0x908a692a, 0xc0e63b24, 0x05de3118, 0xfc974436, 0x1be44823,
  65     0x03de2f3d, 0x66cfb6e4, 0x52727bfc, 0xa7b93651, 0xd7b9035f, 0xfac28d33,
  66     0x59bb4457, 0xeede4004, 0x175ad747, 0x7808d123, 0xc9c97de8, 0x0c26ca26,
  67     0x75e62e96, 0xc8376e97, 0xf2ee6baa, 0x6a885f88, 0x352f92ab, 0x4143f4a4,
  68     0xb1cca58c, 0xe8fbea94, 0x5c306621, 0xfbe64c32, 0xa1ed285d, 0xca7395cf,
  69     0x4eed31a5, 0x31e39fee, 0x7951c585, 0x23434811, 0xfc103036, 0xef001b3c,
  70     0x499f5f34, 0x5f7f38f4, 0x0206d8c5, 0xcc3ee4f1, 0xbc0b485c, 0x4e4f5829,
  71     0x05ee6e6d, 0xc82d5353, 0x44892bec, 0x22984b53, 0x8a6374d1, 0x0850c3f9,
  72     0x0c06ae88, 0x2dfdc126, 0xd1edacdc, 0x1d8dbd39, 0xdeff2db8, 0xd557278d,
  73     0x7e9e3740, 0x49a1ecb5, 0x43f7b391, 0x50b6b9ef, 0x46b9b8f8, 0xd3f5f6d2,
  74     0x8d453b88, 0xc0ba5333, 0x5ab92e37, 0x6e7620a4, 0x8eb9795a, 0x30355a84,
  75     0xf5e4ad33, 0x7d0b4df2, 0xe0f3e2a1, 0xa466f0e6, 0x39a19c9a, 0x1b284524,
  76     0x854f8b3b, 0x02d10b6c, 0x44fb5d9d, 0x60c29fec, 0xda35244a, 0xb5ce6653,
  77     0xfd8356ad, 0xff88d46b, 0x23fd1d16, 0xdc0be23c };
  78
  79 // Random hash function on 32 bit integers.  Works by XORing the input
  80 // word with s-box values looked up from each input byte.  This is
  81 // pretty much the simplest "good" hash function of this type.  The
  82 // instruction count is very low; depending on cache behavior with the
  83 // 1024 byte s-box table, it may or may not be the fastest.
  84 inline unsigned int Turbulence::hashrand(unsigned int i)
  85 {
  86     i ^= SBOX[(i>> 0) & 0xff];
  87     i ^= SBOX[(i>> 8) & 0xff];
  88     i ^= SBOX[(i>>16) & 0xff];
  89     i ^= SBOX[(i>>24) & 0xff];
  90     return i;
  91 }
  92
  93 // 32 bit integer to [0:1] (safe with 64 bit ints)
  94 static inline float i2fu(unsigned int i) { return (1.0/0xffffffffu) * i; }
  95
  96 // 32 bit integer to [-1:1] (safe with 64 bit ints)
  97 static inline float i2fs(unsigned int i) { return 2 * i2fu(i) - 1; }
  98
  99 // Similar conversions, for 8 bit unsigned bytes
 100 static inline float c2fu(unsigned char c) { return (c+0.5)*(1.0/256); }
 101 static inline unsigned char f2cu(float f) {
 102     int c = (int)(f * 256);
 103     return c == 256 ? 255 : c;
 104 }
 105
 106 inline void Turbulence::turblut(int x, int y, float* out)
 107 {
 108     x = x >= _sz ? x - _sz : x;
 109     y = y >= _sz ? y - _sz : y;
 110
 111     unsigned char* turb = _data + 3*(y*_sz+x);
 112     out[0] = c2fu(turb[0]) * (_x1 - _x0) + _x0;
 113     out[1] = c2fu(turb[1]) * (_y1 - _y0) + _y0;
 114     out[2] = c2fu(turb[2]) * (_z1 - _z0) + _z0;
 115 }
 116
 117 void Turbulence::setMagnitude(double mag)
 118 {
 119     _mag = Math::pow(mag, MAGNITUDE_EXP);
 120 }
 121
 122 void Turbulence::update(double dt, double rate)
 123 {
 124     _timeOff += BASE_RATE * dt * rate;
 125 }
 126
 127 void Turbulence::offset(float* offset)
 128 {
 129     for(int i=0; i<3; i++)
 130         _off[i] += offset[i];
 131 }
 132
 133 void Turbulence::getTurbulence(double* loc, float alt, float* up,
 134                                float* turbOut)
 135 {
 136     // Convert to integer 2D coordinates; wrap to [0:_sz].
 137     double a = (loc[0] + _off[0]) + (loc[2] + _off[2]);
 138     double b = (loc[1] + _off[1]) + _timeOff;
 139     a -= _sz * Math::floor(a * (1.0/_sz));
 140     b -= _sz * Math::floor(b * (1.0/_sz));
 141     int x = ((int)Math::floor(a))&(_sz-1);
 142     int y = ((int)Math::floor(b))&(_sz-1);
 143
 144     // Convert to fractional interpolation factors
 145     a -= x;
 146     b -= y;
 147
 148     // Do the lookups
 149     float turb00[3], turb10[3], turb01[3], turb11[3];
 150     turblut(x,     y, turb00);
 151     turblut(x+1,   y, turb10);
 152     turblut(x,   y+1, turb01);
 153     turblut(x+1, y+1, turb11);
 154
 155     // Interpolate, add in units
 156     float mag = _mag * MAX_TURBULENCE;
 157     for(int i=0; i<3; i++) {
 158         float avg0 = (1-a)*turb00[i] + a*turb01[i];
 159         float avg1 = (1-a)*turb10[i] + a*turb11[i];
 160         turbOut[i] = mag * ((1-b)*avg0 + b*avg1);
 161     }
 162
 163     // Adjust for altitude effects
 164     if(alt < 300) {
 165         float altmul = 0.5 + (1-0.5) * (alt*(1.0/300));
 166         if(alt < 100) {
 167             float vmul = alt * (1.0/100);
 168             vmul = vmul / altmul; // pre-correct for the pending altmul
 169             float dot = Math::dot3(turbOut, up);
 170             float off[3];
 171             Math::mul3(dot * (vmul-1), up, off);
 172             Math::add3(turbOut, off, turbOut);
 173         }
 174         Math::mul3(altmul, turbOut, turbOut);
 175     }
 176 }
 177
 178 // Associates a random number in the range [-1:1] with a given lattice
 179 // point.
 180 float Turbulence::lattice(unsigned int x, unsigned int y)
 181 {
 182     return i2fs(hashrand((((_seed << _gens) | x) << _gens) | y));
 183 }
 184
 185 // Returns a scale for a vector that normalizes it into a sphere (as
 186 // opposed to cube) space.  This prevents the overscaling of the
 187 // "corner" vectors you get from choosing three random turbulence
 188 // components.
 189 float Turbulence::cubenorm(float x, float y, float z)
 190 {
 191     x = x < 0 ? -x : x;
 192     y = y < 0 ? -y : y;
 193     z = z < 0 ? -z : z;
 194     float max = ((x > y) && (x > z)) ? x : ((y > z) ? y : z);
 195     return max/Math::sqrt(x*x + y*y + z*z);
 196 }
 197
 198 Turbulence::~Turbulence()
 199 {
 200     delete[] _data;
 201 }
 202
 203 Turbulence::Turbulence(int gens, int seed)
 204 {
 205     _gens = gens;
 206     _sz = 1 << (_gens - 1);
 207     _seed = seed;
 208     _mag = 1;
 209     _x0 = _x1 = _y0 = _y1 = _z0 = _z1 = 0;
 210     _timeOff = 0;
 211     _off[0] = _off[1] = _off[2] = 0;
 212
 213     float* xbuf = new float[_sz*_sz];
 214     float* ybuf = new float[_sz*_sz];
 215     float* zbuf = new float[_sz*_sz];
 216
 217     mkimg(xbuf, _sz);
 218     _seed++;
 219     mkimg(ybuf, _sz);
 220     _seed++;
 221     mkimg(zbuf, _sz);
 222
 223     // "Normalize" them to proper spherical magnitudes, and calculate
 224     // range information for the packing.
 225     for(int i=0; i<_sz*_sz; i++) {
 226         float n = cubenorm(xbuf[i], ybuf[i], zbuf[i]);
 227         xbuf[i] *= n;
 228         ybuf[i] *= n;
 229         zbuf[i] *= n;
 230
 231         _x0 = xbuf[i] < _x0 ? xbuf[i] : _x0;
 232         _x1 = xbuf[i] > _x1 ? xbuf[i] : _x1;
 233         _y0 = ybuf[i] < _y0 ? ybuf[i] : _y0;
 234         _y1 = ybuf[i] > _y1 ? ybuf[i] : _y1;
 235         _z0 = zbuf[i] < _z0 ? zbuf[i] : _z0;
 236         _z1 = zbuf[i] > _z1 ? zbuf[i] : _z1;
 237     }
 238
 239     // Pack into 3 byte tuples for storage.
 240     _data = new unsigned char[3*_sz*_sz];
 241     for(int i=0; i<_sz*_sz; i++) {
 242         float x = xbuf[i], y = ybuf[i], z = zbuf[i];
 243         unsigned char* tuple = _data + 3*i;
 244         tuple[0] = f2cu((x - _x0) / (_x1 - _x0));
 245         tuple[1] = f2cu((y - _y0) / (_y1 - _y0));
 246         tuple[2] = f2cu((z - _z0) / (_z1 - _z0));
 247     }
 248
 249     delete[] xbuf;
 250     delete[] ybuf;
 251     delete[] zbuf;
 252 }
 253
 254 // "Integer" turbulence function.  Takes coordinates in the range
 255 // [0:1] expressed as a fraction of 2^32 (works with 64 bit ints too;
 256 // it just doesn't use the whole range). The output range is
 257 // guaranteed to be within [-1:1], with a typical output range of +/-
 258 // 0.6 or so.
 259 float Turbulence::iturb(unsigned int x, unsigned int y)
 260 {
 261     float amplitude = 0.5; // start here, so it all sums to ~1.0
 262     float total = 0;
 263     int wrapmax = 2;
 264     int startgen = _gens - MEANINGFUL_GENS;
 265     for(int g=startgen; g<_gens; g++) {
 266         int xl = x >> (32 - g); // lattice coordinates
 267         int yl = y >> (32 - g);
 268         float xfrac = i2fu(x << g); // interpolation fractions
 269         float yfrac = i2fu(y << g);
 270         xfrac = xfrac*xfrac*(3 - 2*xfrac); // ... as cubics
 271         yfrac = yfrac*yfrac*(3 - 2*yfrac);
 272
 273         float p00 = lattice(xl,   yl); // lattice values
 274         float p01 = lattice(xl,   yl+1);
 275         float p10 = lattice(xl+1, yl);
 276         float p11 = lattice(xl+1, yl+1);
 277
 278         float p0 = p00 * (1-yfrac) + p01 * yfrac;
 279         float p1 = p10 * (1-yfrac) + p11 * yfrac;
 280         float p = p0 * (1-xfrac) + p1 * xfrac;
 281         total += p * amplitude;
 282         amplitude *= 0.5;
 283         wrapmax *= 2;
 284     }
 285     return total;
 286 }
 287
 288 // Converts "real" turbulence coordinates expressed in the range
 289 // [0:_sz] (modulo) to integers and runs them through iturb().
 290 float Turbulence::fturb(double a, double b)
 291 {
 292     a *= 1.0 / _sz;
 293     b *= 1.0 / _sz;
 294     a -= Math::floor(a);
 295     b -= Math::floor(b);
 296     return iturb((unsigned int)(a * 4294967296.0),
 297                  (unsigned int)(b * 4294967296.0));
 298 }
 299
 300 void Turbulence::mkimg(float* buf, int sz)
 301 {
 302     for(int y=0; y<sz; y++)
 303         for(int x=0; x<sz; x++)
 304             buf[y*sz+x] = fturb(x + 0.5, y + 0.5);
 305 }
 306
 307 }; // namespace yasim