Wind Visualization
Last Modified: Nov 4, 2018 7:41 AM
Making a visualization of the wind, similar to some that have appeared online for every hurricane recently.
JSL below with the parameters turned down for better interactive performance.
The model the code is using for the wind field is designed to make a pretty visualization and does not represent reality.
The vectors for the video look like this; the visualization changes when you drag the vectors.
NTRAIL = 50; // <<< try 3 to 50 or so -- longer trail of fading dots
NPARTICLES = 1000; // try 100...10000. interacts with NBINS, need less bins with less particles
SLOWER=50; // try 20 to 200 or so, smaller divisor makes bigger steps -- interacts with smoothing that affects cornering ability
SMOOTHING = 9; // try 1..10 or so to smooth the corners; smaller is more accurate, larger is prettier
NBINS=15; // try 2 to 19 or so. determines the granularity of the spawning locations where low densities are found
WINDFIELDSIZE = 40; // granularity of field, interacts with smoothing
PARTICLESIZE = 1; // try 0..10 or so may affect rendering speed
FADEPOWER = 1.5; // power tweaks fade curve
FADEMIN=20/255; // bottom end of fade curve 1..FADEMIN - darkest color that displays
graphxsize=1920/2;
graphysize=1080/2;
//outdir="F:\weather\windfieldpics\";
//==================
// flow visualization -- make a wind field
dx = J( WINDFIELDSIZE, WINDFIELDSIZE, 0 );
dy = J( WINDFIELDSIZE, WINDFIELDSIZE, 0 );
Parallel Assign( {SLOWER=SLOWER}, dx[r( -1, 1 ), c( -1, 1 )] = Sin( r ) / SLOWER );
Parallel Assign( {SLOWER=SLOWER}, dy[r( -1, 1 ), c( -1, 1 )] = Sin( c ) / SLOWER ); //
windlist = {{{0.171875, 0.0458333333333334}, {0.15125, -0.715}}, {{0.171875, -0.9075},
{0.763125, -0.348333333333333}}, {{-0.859375, -0.504166666666667}, {-0.1925, -0.9075
}}, {{-0.385, 0.889166666666667}, {-0.955625, 0.119166666666667}}, {{0.78375, -0.165
}, {0.268125, 0.128333333333333}}, {{-0.240625, 0.0458333333333334}, {-0.78375,
-0.22}}, {{0.86625, -0.0183333333333333}, {0.2475, 0.8525}}, {{-0.144375,
-0.705833333333333}, {-0.2475, -0.0641666666666667}}};
/*
{
{{.10,.20},{.50,.30}},
{{.20,.20},{.50,.40}},
{{.30,.20},{.50,.50}},
{{.40,.20},{.50,.60}},
{{.50,.20},{.50,.70}},
{{.60,.20},{.50,.80}},
{{.70,.20},{.50,.90}}
};
*/
caprow=.;
capcol=.;
ReleaseFunction = Function( {x, y},
{},
caprow = capcol = .;
// put the wind vectors into the dx/dy matrices
Parallel Assign(
{SLOWER = SLOWER, windlist = windlist},
dx[r( -1, 1 ), c( -1, 1 )] = (x = 0 ; totalscale = 0 ; For( iwind = 1, iwind <= N Items( windlist ), iwind += 1,
cx = (windlist[iwind][1][1] + windlist[iwind][2][1]) / 2;
cy = (windlist[iwind][1][2] + windlist[iwind][2][2]) / 2;
dx = windlist[iwind][2][1] - windlist[iwind][1][1]; // positive toward arrow
scale = 1 / (Sqrt( (c - cx) ^ 2 + (r - cy) ^ 2 ) + .001);
totalscale += scale;
x += dx * scale;
) ; x / totalscale / SLOWER)
);
Parallel Assign(
{SLOWER = SLOWER, windlist = windlist},
dy[r( -1, 1 ), c( -1, 1 )] = (y = 0 ; totalscale = 0 ; For( iwind = 1, iwind <= N Items( windlist ), iwind += 1,
cx = (windlist[iwind][1][1] + windlist[iwind][2][1]) / 2;
cy = (windlist[iwind][1][2] + windlist[iwind][2][2]) / 2;
dy = windlist[iwind][2][2] - windlist[iwind][1][2]; // positive toward arrow
scale = 1 / (Sqrt( (c - cx) ^ 2 + (r - cy) ^ 2 ) + .001);
totalscale += scale;
y += dy * scale;
) ; y / totalscale / SLOWER)
);
);
ReleaseFunction(0,0);
New Window( "Example",
Graph Box(
X Scale( -1.1, 1.1 ),
Y Scale( -1.1, 1.1 ),
Pen Size( 4 );
For( iwind = 1, iwind <= N Items( windlist ), iwind += 1,
Arrow( windlist[iwind][1], windlist[iwind][2] )
);
Mousetrap(
Function( {x, y},
{irow, nearrow = ., icol, nearcol = ., dist, neardist = 1e9},
If( Is Missing( caprow ),
For( irow = 1, irow <= N Items( windlist ), irow += 1,
For( icol = 1, icol <= 2, icol += 1,
dist = Sqrt( (windlist[irow][icol][1] - x) ^ 2 + (windlist[irow][icol][2] - y) ^ 2 );
If( dist < neardist,
nearrow = irow;
nearcol = icol;
neardist = dist;
);
)
);
If( !Is Missing( nearrow ),
caprow = nearrow;
capcol = nearcol;
);
);
If( !Is Missing( caprow ) & !Is Missing( capcol ),
windlist[caprow][capcol][1] = x;
windlist[caprow][capcol][2] = y;
);
)//move
, //
ReleaseFunction(x,y)//release
);
)
);
// ===============================
SPREAD=1/NBINS; // mu for random normal spawn location
density = J( NBINS, NBINS, 0 );
//density[1, 1] = 1; // first loc max needs something to work with
_DX_ = 1;
_DY_ = 2;
_COLOR_ = 3;
_TRAILIDX_ = 4;
_TRAIL_ = 5; // NTRAIL begins here
_X_ = _TRAIL_;
_Y_ = _TRAIL_ + 1;
OUTOFPLAY = 99; // well beyond the graph's edge
DFRACTION=SMOOTHING;
NFRACTION = DFRACTION-1;
tran = J( NTRAIL, 1, 0 );
Parallel Assign( {FADEPOWER=FADEPOWER,FADEMIN=FADEMIN^(1/FADEPOWER)}, tran[r( 1, FADEMIN ), c] = r ^ FADEPOWER );
_STRAIL_ = 2; // stride x y pairs
particles = J( NPARTICLES , _TRAIL_ - 1 + _STRAIL_ * NTRAIL, OUTOFPLAY ); // dx,dy,color,xy trail has alpha fade
particles[0, _DX_] = 0;
particles[0, _DY_] = 0;
particles[0, _TRAILIDX_] = NTRAIL; // used to fade in, fade out
ExpectedParticlesPerBin = NPARTICLES/nitems(density);
New Window( "Example",
g = Graph Box(
Frame Size( graphxsize, graphysize ),
X Scale( -1.1, 1.1 ),
Y Scale( -1.1, 1.1 ),
suppressaxes,
<<( PARTICLESIZE ),
<<( RGB Color( 0, .1, .2 ) ), // contrast with "white", below
minDenLocs = Loc( density == Min( density ) );
minDenYLocs = Floor( (minDenLocs - 1) / N Cols( Density ) + 1 );
minDenXLocs = minDenLocs - (minDenYLocs - 1) * N Cols( Density );
yymin = .5;
yymax = N Rows( Density ) + .5;
minDenYLocs = 2 * ((minDenYLocs - yymin) / (yymax - yymin)) - 1;
xxmin = .5;
xxmax = N Cols( Density ) + .5;
minDenXLocs = 2 * ((minDenXLocs - xxmin) / (xxmax - xxmin)) - 1;
density[0, 0] = 0;
For( i = 1, i <= N Rows( particles ), i += 1,
fadeIdx = particles[i, _TRAILIDX_];
If(
fadeIdx > NTRAIL, // point is in the fading sequence
// shift old particles out
particles[i, (_TRAIL_ + (_STRAIL_ * 1)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL) - 1)] = particles[i, (_TRAIL_ + (_STRAIL_ * 0)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL - 1) - 1)];
nn = (particles[i, _TRAILIDX_] -= 1); // counting time while previous point fades away
If( nn == NTRAIL, // restart this completely faded point
rr = Random Integer( 1, N Rows( minDenXLocs ) );
tx = particles[i, _X_ + _STRAIL_ * (NTRAIL - 1)] = Max( -.999, Min( .999, Random Normal( minDenXLocs[rr], SPREAD ) ) ); // next to last
ty = particles[i, _Y_ + _STRAIL_ * (NTRAIL - 1)] = Max( -.999, Min( .999, Random Normal( minDenYLocs[rr], SPREAD ) ) );
// index to field data
idx = 1 + Round( (N Cols( dx ) - 1) * ((tx + 1) / 2) ); // tx+1 is 0...2
idy = 1 + Round( (N Rows( dx ) - 1) * ((ty + 1) / 2) ); // divide by 2 is 0...1 becomes 0..nr-1...round, +1 becomes 1...nr
particles[i, _DX_] = dx[idy, idx]; // dx,dy is the field speed. init the particles DX,DY
particles[i, _DY_] = dy[idy, idx]; // to the correct value
// the values just stored at NTRAIL-2 will display, slightly bright, on this cycle
// and move to NTRAIL-1 (dimmest) on next cycle
);//
, /* elseif */ fadeIdx >= 1, // point is in the start sequence, or started
If( fadeIdx > 1, // still fading in
(particles[i, _TRAILIDX_] -= 1);// for next round
particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx - 2)] = particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx-1)] + particles[i, _DX_::_DY_];
,// else...sticks at 1
// shift old particles out
particles[i, (_TRAIL_ + (_STRAIL_ * 1)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL) - 1)] = particles[i, (_TRAIL_ + (_STRAIL_ * 0)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL - 1) - 1)];
particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx - 1)] = particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx-0)] + particles[i, _DX_::_DY_];
);
If(
!(-1 < particles[i, _X_ + _STRAIL_ * (fadeIdx - 1)] < 1) | !(-1 < particles[i, _Y_ + _STRAIL_ * (fadeIdx - 1)] < 1) |
Sqrt( particles[i, _DX_] ^ 2 + particles[i, _DY_] ^ 2 ) < .0001,
// fix up particles leaving graph
// leaving because a 99,99 was incremented by dx/dy
particles[i, _TRAILIDX_] = 2 * NTRAIL;// enter fade-out sequence
particles[i, _X_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
particles[i, _Y_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
particles[i, _DX_] = 0;
particles[i, _DY_] = 0; //
, // else
// index to field data
idx = 1 + Round( (N Cols( dx ) - 1) * ((particles[i, _X_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
idy = 1 + Round( (N Rows( dx ) - 1) * ((particles[i, _Y_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
// new deltas based on nearest field delta
particles[i, _DX_] = (NFRACTION * particles[i, _DX_] + dx[idy, idx]) / DFRACTION;
particles[i, _DY_] = (NFRACTION * particles[i, _DY_] + dy[idy, idx]) / DFRACTION;
// nearest density bin...increment
yyy = 1 + Round( (N Rows( density ) - 1) * ((particles[i, _Y_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
xxx = 1 + Round( (N Cols( density ) - 1) * ((particles[i, _X_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
q=((density[yyy, xxx] += 1))/ExpectedParticlesPerBin;
if( q*randomuniform()>2,
particles[i, _TRAILIDX_] = 2 * NTRAIL;// enter fade-out sequence
particles[i, _X_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
particles[i, _Y_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
particles[i, _DX_] = 0;
particles[i, _DY_] = 0; //
);
);//
, // else bad _TRAILIDX_
Throw( "???" )
);
);
For( i = NTRAIL, i >= 1, i -= 1,
Transparency( tran[i] );
Marker( Color State( "white" ), particles[0, _X_ + _STRAIL_ * (i - 1)],/*flip -*/ particles[0, _Y_ + _STRAIL_ * (i - 1)] );
);
)
);
wait(.01);
picno = 1e6;
while(1, g <<<<;wait(0);
//pic = newimage(((g[framebox(1)]<<getimage)<<getpixels)[3::1082,3::1922]);
//pic<<saveimage(outdir||char(picno+=1)||".png","png");
);
Comments