/* stream.c */ /* * Vis5D system for visualizing five dimensional gridded data sets. * Copyright (C) 1990 - 2000 Bill Hibbard, Johan Kellum, Brian Paul, * Dave Santek, and Andre Battaiola. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * As a special exception to the terms of the GNU General Public * License, you are permitted to link Vis5D with (and distribute the * resulting source and executables) the LUI library (copyright by * Stellar Computer Inc. and licensed for distribution with Vis5D), * the McIDAS library, and/or the NetCDF library, where those * libraries are governed by the terms of their own licenses. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * */ #include "../config.h" /* 2-D streamline making function */ #include #include #include #include "memory.h" #include "globals.h" #include "proj.h" #define GU( R, C ) ( ugrid[ (R) * nc + (C) ] ) #define GV( R, C ) ( vgrid[ (R) * nc + (C) ] ) #define MARKARROW( R, C ) ( markarrow[ (C) * nrstart + (R) ] ) #define MARKSTART( R, C ) ( markstart[ (C) * nrstart + (R) ] ) #define MARKEND( R, C ) ( markend[ (C) * nrend + (R) ] ) #define START2ROW( IRLOW ) ( ((float) nr-1.0) * ((float) (IRLOW)+0.5) / (float) nrstart ) #define START2COL( ICLOW ) ( ((float) nc-1.0) * ((float) (ICLOW)+0.5) / (float) ncstart ) #define ROW2ARROW( ROW ) ( (int) (nrarrow * (ROW) / ((float) nr-1.0) ) ) #define COL2ARROW( COL ) ( (int) (ncarrow * (COL) / ((float) nc-1.0) ) ) #define ROW2START( ROW ) ( (int) (nrstart * (ROW) / ((float) nr-1.0) ) ) #define COL2START( COL ) ( (int) (ncstart * (COL) / ((float) nc-1.0) ) ) #define ROW2END( ROW ) ( (int) (nrend * (ROW) / ((float) nr-1.0) ) ) #define COL2END( COL ) ( (int) (ncend * (COL) / ((float) nc-1.0) ) ) #define LENGTH 0.02 /* * trace one stream line for 2-D u and v arrays. * Note that the input arrays ugrid & vgrid should be in * column-major (FORTRAN) order. * * Input: ugrid, vgrid - the 2-D wind arrays. * nr, nc - size of 2-D array in rows and columns. * vr, vc - arrays to put streamline vertices * dir - direction 1.0=forward, -1.0=backward * maxv - size of vx, vy arrays * numv - pointer to int to return number of vertices in vx, vy * markarrow - mark array for arrows * markstart - mark array for starting streamlines * markend - mark array for ending streamlines * nrarrow, ncarrow - size of markarrow * nrstart, ncstart - size of markstart * nrend, ncend - size of markend * row, col - start location for streamline * step - step size for streamline * rowlength, collength - scale constants for arrows * irend, icend - location of most recent end box * Return: 1 = ok * 0 = no more memory for streamlines */ int stream_trace( Context ctx, float ugrid[], float vgrid[], int nr, int nc, float dir, float vr[], float vc[], int maxv, int *numv, char *markarrow, char *markstart, char *markend, int nrarrow, int ncarrow, int nrstart, int ncstart, int nrend, int ncend, float row, float col, float step, float rowlength, float collength, int irend, int icend) { int ir, ic, ire, ice; int ira, ica, irs, ics; int nend, num; float prevrow, prevcol; float a, c, ac, ad, bc, bd; float u, v; num = *numv; nend = 0; while (1) { float ubd, ubc, uad, uac, vbd, vbc, vad, vac; /* interpolate velocity at row, col */ ir = row; ic = col; a = row - ir; c = col - ic; ac = a*c; ad = a*(1.0-c); bc = (1.0-a)*c; bd = (1.0-a)*(1.0-c); ubd = GU(ir, ic); ubc = GU(ir, ic+1); uad = GU(ir+1, ic); uac = GU(ir+1, ic+1); vbd = GV(ir, ic); vbc = GV(ir, ic+1); vad = GV(ir+1, ic); vac = GV(ir+1, ic+1); /* terminate stream if missing data */ if (IS_MISSING(ubd) || IS_MISSING(ubc) || IS_MISSING(uad) || IS_MISSING(uac) || IS_MISSING(vbd) || IS_MISSING(vbc) || IS_MISSING(vad) || IS_MISSING(vac)) break; u = bd * ubd + bc * ubc + ad * uad + ac * uac; v = bd * vbd + bc * vbc + ad * vad + ac * vac; /* WLH 5 July 96 u = bd * GU(ir, ic) + bc * GU(ir, ic+1) + ad * GU(ir+1, ic) + ac * GU(ir+1, ic+1); v = bd * GV(ir, ic) + bc * GV(ir, ic+1) + ad * GV(ir+1, ic) + ac * GV(ir+1, ic+1); */ /* scale velocity */ /* WLH 7-3-96 u = step * ctx->dpy_ctx->Uscale[ir][ic] * u; v = step * ctx->dpy_ctx->Vscale[ir][ic] * v; */ u = step * u; v = step * v; /* test for too many line segments */ if (num > maxv-2) { deallocate( ctx, markarrow, nrstart * ncstart * sizeof(char) ); deallocate( ctx, markstart, nrstart * ncstart * sizeof(char) ); deallocate( ctx, markend, nrend * ncend * sizeof(char) ); *numv = num; return 0; } /* propogate streamline */ prevrow = row; prevcol = col; row += dir * v; col += dir * u; /* terminate stream if out of grid */ if (row < 0 || col < 0 || row >= nr-1 || col >= nc-1) { break; } ire = ROW2END(row); ice = COL2END(col); /* terminate stream if enters marked end box */ if (ire != irend || ice != icend) { irend = ire; icend = ice; if (irend < 0 || irend >= nrend || icend < 0 || icend >= ncend) { printf("bad 2: irend = %d icend = %d\n", irend, icend); } if (MARKEND(irend, icend) == 1) { break; } MARKEND(irend, icend) = 1; nend = 0; } /* terminate stream if too many steps in one end box */ nend++; if (nend > 100) { break; } /* make line segment */ vr[num] = prevrow; vc[num++] = prevcol; vr[num] = row; vc[num++] = col; /* mark start box */ irs = ROW2START(row); ics = COL2START(col); if (irs < 0 || irs >= nrstart || ics < 0 || ics >= ncstart) { printf("bad 3: irs = %d ics = %d\n", irs, ics); } if (MARKSTART(irs, ics) == 0) { MARKSTART(irs, ics) = 1; } /* check for need to draw arrow head */ ira = ROW2ARROW(row); ica = COL2ARROW(col); if (MARKARROW(ira, ica) == 0) { double rv, cv, v; /* test for too many line segments */ if (num > maxv-4) { deallocate( ctx, markarrow, nrstart * ncstart * sizeof(char) ); deallocate( ctx, markstart, nrstart * ncstart * sizeof(char) ); deallocate( ctx, markend, nrend * ncend * sizeof(char) ); *numv = num; return 0; } MARKARROW(ira, ica) = 1; rv = dir * (row - prevrow); cv = dir * (col - prevcol); v = sqrt(rv*rv + cv*cv); if (v > 0.000000001) { rv = rv / v; cv = cv / v; } vr[num] = row; vc[num++] = col; vr[num] = row - (rv + cv) * rowlength; vc[num++] = col + (rv - cv) * collength; vr[num] = row; vc[num++] = col; vr[num] = row + (cv - rv) * rowlength; vc[num++] = col - (cv + rv) * collength; } } /* end while (forward) */ *numv = num; return 1; } /* * Compute stream lines for 2-D u and v arrays. * Note that the input arrays ugrid & vgrid should be in * row-major order. * * Input: ugrid, vgrid - the 2-D wind arrays. * nr, nc - size of 2-D array in rows and columns. * density - relative density of streamlines. * vr, vc - arrays to put streamline vertices * maxv - size of vx, vy arrays * numv - pointer to int to return number of vertices in vx, vy * Return: 1 = ok * 0 = error (out of memory) */ int stream( Context ctx, float ugrid[], float vgrid[], int nr, int nc, float density, float vr[], float vc[], int maxv, int *numv) { int irstart, icstart, irend, icend, ir, ic; int nrarrow, ncarrow, nrstart, ncstart, nrend, ncend; int num; char *markarrow, *markstart, *markend; float row, col; float step, rowlength, collength; float dir; /* initialize vertex counts */ num = 0; /* density calculations */ if (density < 0.5) density = 0.5; if (density > 2.0) density = 2.0; nrarrow = 15.0001 * density; ncarrow = 15.0001 * density; nrstart = 15.0001 * density; ncstart = 15.0001 * density; nrend = 4 * nrstart; ncend = 4 * ncstart; rowlength = LENGTH * nr / density; collength = LENGTH * nc / density; /* WLH - use shorter step for higher density */ step = ctx->TrajStep / density; /* allocate mark arrays */ markarrow = (char *) allocate( ctx, nrstart * ncstart * sizeof(char) ); if (!markarrow) return 0; markstart = (char *) allocate( ctx, nrstart * ncstart * sizeof(char) ); if (!markstart) return 0; markend = (char *) allocate( ctx, nrend * ncend * sizeof(char) ); if (!markend) return 0; /* initialize mark array to zeros */ memset( markstart, 0, nrstart * ncstart * sizeof(char) ); memset( markend, 0, nrend * ncend * sizeof(char) ); /* only draw arrows in every ninth box */ memset( markarrow, 1, nrstart * ncstart * sizeof(char) ); for (ir = 1; ir= nrend || icend < 0 || icend >= ncend) { printf("bad 1: irend = %d icend = %d\n", irend, icend); } MARKEND(irend, icend) = 1; dir = 1.0; if (stream_trace( ctx, ugrid, vgrid, nr, nc, dir, vr, vc, maxv, &num, markarrow, markstart, markend, nrarrow, ncarrow, nrstart, ncstart, nrend, ncend, row, col, step, rowlength, collength, irend, icend) == 0) { *numv = num; return 1; } /* now trace streamline backward */ row = START2ROW(irstart); col = START2COL(icstart); irend = ROW2END(row); icend = COL2END(col); if (irend < 0 || irend >= nrend || icend < 0 || icend >= ncend) { printf("bad 3: irend = %d icend = %d\n", irend, icend); } MARKEND(irend, icend) = 1; dir = -1.0; if (stream_trace( ctx, ugrid, vgrid, nr, nc, dir, vr, vc, maxv, &num, markarrow, markstart, markend, nrarrow, ncarrow, nrstart, ncstart, nrend, ncend, row, col, step, rowlength, collength, irend, icend) == 0) { *numv = num; return 1; } } /* end if */ } /* end for */ } /* end for */ /* deallocate mark arrays */ deallocate( ctx, markarrow, nrstart * ncstart * sizeof(char) ); deallocate( ctx, markstart, nrstart * ncstart * sizeof(char) ); deallocate( ctx, markend, nrend * ncend * sizeof(char) ); *numv = num; return 1; }