/* traj.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"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "globals.h"
#include "grid.h"
#include "proj.h"
#define TRUNC(X) ( (float) (int) (X) )
/*
* Initialize the trajectory tracing module.
* Return: 1 = ok, 0 = error
*/
int init_traj( Context ctx )
{
int i, j;
float lata, lona, latb, lonb, latc, lonc;
float d, us, vs;
float lat0, lon0, lat1, lon1, lat2, lon2;
float midrow, midcol;
int var = 0;
/* the index of any wind variable */
if (ctx->dpy_ctx->TrajU>=0) var = ctx->dpy_ctx->TrajU;
else if (ctx->dpy_ctx->TrajV>=0) var= ctx->dpy_ctx->TrajV;
else if (ctx->dpy_ctx->TrajW>=0) var= ctx->dpy_ctx->TrajW;
/* Compute initial trajectory step and length values */
switch (ctx->Projection) {
case PROJ_GENERIC:
ctx->TrajStep = 1.0;
ctx->TrajLength = 1.0;
break;
default:
/* TODO: verify this is ok: (seems to work) */
/* This is tricky: compute distance, in meters, from left to */
/* right edge of domain. Do it in two steps to prevent "wrap- */
/* around" when difference in longitudes > 180 degrees. */
midrow = (float) ctx->Nr / 2.0;
midcol = (float) ctx->Nc / 2.0;
rowcol_to_latlon( ctx, -1, -1, midrow, 0.0, &lat0, &lon0 );
rowcol_to_latlon( ctx, -1, -1, midrow, midcol, &lat1, &lon1 );
rowcol_to_latlon( ctx, -1, -1, midrow, (float) (ctx->Nc-1), &lat2, &lon2 );
d = earth_distance( lat0, lon0, lat1, lon1 )
+ earth_distance( lat1, lon1, lat2, lon2 );
ctx->TrajStep = TRUNC( (d / 100.0) / 25.0 );
/* the above was: (float) ctx->Elapsed[1] / 2.0;*/
ctx->TrajLength = 5.0 * ctx->TrajStep;
/* the above was: (float) ctx->Elapsed[1]; */
break;
}
/* These values are set by the user through the trajectory widget */
/* They are just multipliers of the internal values TrajStep and */
/* TrajLength: */
ctx->dpy_ctx->UserTrajStep = ctx->dpy_ctx->UserTrajLength = 1.0;
/*
* Compute m/s to boxes/s scaling factors for U and V components
*/
switch (ctx->Projection) {
case PROJ_GENERIC:
/* for a generic projection, we assume U, and V velocities are in */
/* X per second, where X is the same units used for NorthBound, */
/* WestBound, ctx->RowInc, and ctx->ColInc */
us = 1.0 / ctx->ColInc;
vs = -1.0 / ctx->RowInc;
for (i=0;i<ctx->Nr;i++) {
for (j=0;j<ctx->Nc;j++) {
ctx->Uscale[i][j] = us;
ctx->Vscale[i][j] = vs;
}
}
break;
case PROJ_LINEAR:
case PROJ_LAMBERT:
case PROJ_STEREO:
case PROJ_ROTATED:
case PROJ_CYLINDRICAL:
case PROJ_SPHERICAL:
case PROJ_MERCATOR:
for (i=0;i<ctx->Nr;i++) {
for (j=0;j<ctx->Nc;j++) {
float ii = (float) i;
float jj = (float) j;
/* Compute U scale */
if (j==0) {
rowcol_to_latlon( ctx, 0, var, ii, jj, &lata, &lona );
rowcol_to_latlon( ctx, 0, var, ii, jj+1.0, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else if (j==ctx->Nc-1) {
rowcol_to_latlon( ctx, 0, var, ii, jj-1.0, &lata, &lona );
rowcol_to_latlon( ctx, 0, var, ii, jj, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else {
rowcol_to_latlon( ctx, 0, var, ii, jj-1.0, &lata, &lona );
rowcol_to_latlon( ctx, 0, var, ii, jj, &latb, &lonb );
rowcol_to_latlon( ctx, 0, var, ii, jj+1.0, &latc, &lonc );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
if (latc > 89.9) latc = 89.9;
if (latc < -89.9) latc = -89.9;
d = (earth_distance( lata,lona, latb,lonb)
+ earth_distance( latb,lonb, latc,lonc)) / 2.0;
}
ctx->Uscale[i][j] = 1.0 / d;
/* Compute V scale */
if (i==0) {
rowcol_to_latlon( ctx, 0, ctx->dpy_ctx->TrajV, ii, jj, &lata, &lona );
rowcol_to_latlon( ctx, 0, ctx->dpy_ctx->TrajV, ii+1.0, jj, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else if (j==ctx->Nc-1) {
rowcol_to_latlon( ctx, 0, ctx->dpy_ctx->TrajV, ii-1.0, jj, &lata, &lona );
rowcol_to_latlon( ctx, 0, ctx->dpy_ctx->TrajV, ii, jj, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else {
rowcol_to_latlon( ctx, 0, ctx->dpy_ctx->TrajV, ii-1.0, jj, &lata, &lona );
rowcol_to_latlon( ctx, 0, ctx->dpy_ctx->TrajV, ii, jj, &latb, &lonb );
rowcol_to_latlon( ctx, 0, ctx->dpy_ctx->TrajV, ii+1.0, jj, &latc, &lonc );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
if (latc > 89.9) latc = 89.9;
if (latc < -89.9) latc = -89.9;
d = (earth_distance( lata,lona, latb,lonb)
+ earth_distance( latb,lonb, latc,lonc)) / 2.0;
}
ctx->Vscale[i][j] = -1.0 / d; /* Note negative!!! */
}
}
break;
default:
printf("Error in init_traj: Projection=%d\n", ctx->Projection);
return 0;
}
/*
* Compute m/s to boxes/s scaling factors for W component
*/
switch (ctx->VerticalSystem) {
case VERT_GENERIC:
for (i=0;i<ctx->MaxNl;i++) {
ctx->Wscale[i] = 1.0 / ctx->LevInc;
}
break;
case VERT_EQUAL_KM:
for (i=0;i<ctx->MaxNl;i++) {
ctx->Wscale[i] = 1.0 / (ctx->LevInc * 1000.0);
}
break;
case VERT_NONEQUAL_MB:
case VERT_NONEQUAL_KM:
for (i=0;i<ctx->MaxNl;i++) {
if (i==0) {
float hgt1 = gridlevel_to_height( ctx, 1.0 );
float hgt0 = gridlevel_to_height( ctx, 0.0 );
float diff = hgt1-hgt0;
if (fabs(diff) < 0.000001) diff = 0.000001;
ctx->Wscale[i] = 1.0 / (diff * 1000.0);
/* ctx->Wscale[i] = 1.0 / ((hgt1-hgt0) * 1000.0); */
}
else if (i==ctx->MaxNl-1) {
float hgt1 = gridlevel_to_height( ctx, (float)(ctx->MaxNl-1) );
float hgt0 = gridlevel_to_height( ctx, (float)(ctx->MaxNl-2) );
float diff = hgt1-hgt0;
if (fabs(diff) < 0.000001) diff = 0.000001;
ctx->Wscale[i] = 1.0 / (diff * 1000.0);
/* ctx->Wscale[i] = 1.0 / ((hgt1-hgt0) * 1000.0); */
}
else {
float a, b;
float hgt2 = gridlevel_to_height( ctx, (float)(i+1) );
float hgt1 = gridlevel_to_height( ctx, (float) i );
float hgt0 = gridlevel_to_height( ctx, (float)(i-1) );
float diffa = hgt1-hgt0;
float diffb = hgt2-hgt1;
if (fabs(diffa) < 0.000001) diffa = 0.000001;
if (fabs(diffb) < 0.000001) diffb = 0.000001;
a = 1.0 / (diffa * 1000.0);
b = 1.0 / (diffb * 1000.0);
/*
float a = 1.0 / ((hgt1-hgt0) * 1000.0);
float b = 1.0 / ((hgt2-hgt1) * 1000.0);
*/
ctx->Wscale[i] = (a+b) / 2.0;
}
}
break;
default:
printf("Error in init_traj: ctx->VerticalSystem=%d\n", ctx->VerticalSystem);
return 0;
}
return 1;
}
/*
* Return U, V, and W for a discrete position in time and space.
* Return: 1 = ok, 0 = missing data.
*/
static int get_discrete_uvw( Context ctx, int time, int r, int c, int l,
float *uptr, float *vptr, float *wptr )
{
float u = get_grid_value( ctx, time, ctx->dpy_ctx->TrajU, r, c, l );
float v = get_grid_value( ctx, time, ctx->dpy_ctx->TrajV, r, c, l );
float w = get_grid_value( ctx, time, ctx->dpy_ctx->TrajW, r, c, l );
if (IS_MISSING(u) || IS_MISSING(v) || IS_MISSING(w)) {
/* missing */
return 0;
}
else {
/* return u,v,w in boxes/sec */
*uptr = u * ctx->Uscale[r][c];
*vptr = v * ctx->Vscale[r][c];
*wptr = w * ctx->Wscale[l];
return 1;
}
}
static int get_discrete_uv( Context ctx, int time, int r, int c, int l,
float *uptr, float *vptr, float *wptr )
{
float u = get_grid_value( ctx, time, ctx->dpy_ctx->TrajU, r, c, l );
float v = get_grid_value( ctx, time, ctx->dpy_ctx->TrajV, r, c, l );
if (IS_MISSING(u) || IS_MISSING(v)) {
/* missing */
return 0;
}
else {
/* return u,v,w in boxes/sec */
*uptr = u * ctx->Uscale[r][c];
*vptr = v * ctx->Vscale[r][c];
*wptr = 0;
return 1;
}
}
int init_trajPRIME( Display_Context dtx )
{
int i, j;
float lata, lona, latb, lonb, latc, lonc;
float d, us, vs;
float lat0, lon0, lat1, lon1, lat2, lon2;
float midrow, midcol;
int var = 0;
dtx->UserTrajStep = dtx->UserTrajLength = 1.0;
switch (dtx->Projection) {
case PROJ_GENERIC:
dtx->TrajStep = 1.0;
dtx->TrajLength = 1.0;
break;
default:
/* TODO: verify this is ok: (seems to work) */
/* This is tricky: compute distance, in meters, from left to */
/* right edge of domain. Do it in two steps to prevent "wrap- */
/* around" when difference in longitudes > 180 degrees. */
midrow = (float) dtx->Nr / 2.0;
midcol = (float) dtx->Nc / 2.0;
rowcolPRIME_to_latlon( dtx, -1, -1, midrow, 0.0, &lat0, &lon0 );
rowcolPRIME_to_latlon( dtx, -1, -1, midrow, midcol, &lat1, &lon1 );
rowcolPRIME_to_latlon( dtx, -1, -1, midrow, (float) (dtx->Nc-1), &lat2, &lon2 );
d = earth_distance( lat0, lon0, lat1, lon1 )
+ earth_distance( lat1, lon1, lat2, lon2 );
dtx->TrajStep = TRUNC( (d / 100.0) / 25.0 );
/* the above was: (float) dtx->Elapsed[1] / 2.0;*/
dtx->TrajLength = 5.0 * dtx->TrajStep;
/* the above was: (float) dtx->Elapsed[1]; */
break;
}
/*
* Compute m/s to boxes/s scaling factors for U and V components
*/
switch (dtx->Projection) {
case PROJ_GENERIC:
/* for a generic projection, we assume U, and V velocities are in */
/* X per second, where X is the same units used for NorthBound, */
/* WestBound, dtx->RowInc, and dtx->ColInc */
us = 1.0 / dtx->ColInc;
vs = -1.0 / dtx->RowInc;
for (i=0;i<dtx->Nr;i++) {
for (j=0;j<dtx->Nc;j++) {
dtx->Uscale[i][j] = us;
dtx->Vscale[i][j] = vs;
}
}
break;
case PROJ_LINEAR:
case PROJ_LAMBERT:
case PROJ_STEREO:
case PROJ_ROTATED:
case PROJ_CYLINDRICAL:
case PROJ_SPHERICAL:
case PROJ_MERCATOR:
for (i=0;i<dtx->Nr;i++) {
for (j=0;j<dtx->Nc;j++) {
float ii = (float) i;
float jj = (float) j;
/* Compute U scale */
if (j==0) {
rowcolPRIME_to_latlon( dtx, 0, var, ii, jj, &lata, &lona );
rowcolPRIME_to_latlon( dtx, 0, var, ii, jj+1.0, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else if (j==dtx->Nc-1) {
rowcolPRIME_to_latlon( dtx, 0, var, ii, jj-1.0, &lata, &lona );
rowcolPRIME_to_latlon( dtx, 0, var, ii, jj, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else {
rowcolPRIME_to_latlon( dtx, 0, var, ii, jj-1.0, &lata, &lona );
rowcolPRIME_to_latlon( dtx, 0, var, ii, jj, &latb, &lonb );
rowcolPRIME_to_latlon( dtx, 0, var, ii, jj+1.0, &latc, &lonc );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
if (latc > 89.9) latc = 89.9;
if (latc < -89.9) latc = -89.9;
d = (earth_distance( lata,lona, latb,lonb)
+ earth_distance( latb,lonb, latc,lonc)) / 2.0;
}
dtx->Uscale[i][j] = 1.0 / d;
/* Compute V scale */
if (i==0) {
rowcolPRIME_to_latlon( dtx, 0, dtx->TrajV, ii, jj, &lata, &lona );
rowcolPRIME_to_latlon( dtx, 0, dtx->TrajV, ii+1.0, jj, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else if (j==dtx->Nc-1) {
rowcolPRIME_to_latlon( dtx, 0, dtx->TrajV, ii-1.0, jj, &lata, &lona );
rowcolPRIME_to_latlon( dtx, 0, dtx->TrajV, ii, jj, &latb, &lonb );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
d = earth_distance( lata, lona, latb, lonb );
}
else {
rowcolPRIME_to_latlon( dtx, 0, dtx->TrajV, ii-1.0, jj, &lata, &lona );
rowcolPRIME_to_latlon( dtx, 0, dtx->TrajV, ii, jj, &latb, &lonb );
rowcolPRIME_to_latlon( dtx, 0, dtx->TrajV, ii+1.0, jj, &latc, &lonc );
/* WLH 3-21-97 */
if (lata > 89.9) lata = 89.9;
if (lata < -89.9) lata = -89.9;
if (latb > 89.9) latb = 89.9;
if (latb < -89.9) latb = -89.9;
if (latc > 89.9) latc = 89.9;
if (latc < -89.9) latc = -89.9;
d = (earth_distance( lata,lona, latb,lonb)
+ earth_distance( latb,lonb, latc,lonc)) / 2.0;
}
dtx->Vscale[i][j] = -1.0 / d; /* Note negative!!! */
}
}
break;
default:
printf("Error in init_trajPRIME: Projection=%d\n", dtx->Projection);
return 0;
}
/*
* Compute m/s to boxes/s scaling factors for W component
*/
switch (dtx->VerticalSystem) {
case VERT_GENERIC:
for (i=0;i<dtx->MaxNl;i++) {
dtx->Wscale[i] = 1.0 / dtx->LevInc;
}
break;
case VERT_EQUAL_KM:
for (i=0;i<dtx->MaxNl;i++) {
dtx->Wscale[i] = 1.0 / (dtx->LevInc * 1000.0);
}
break;
case VERT_NONEQUAL_MB:
case VERT_NONEQUAL_KM:
for (i=0;i<dtx->MaxNl;i++) {
if (i==0) {
float hgt1 = gridlevelPRIME_to_height( dtx, 1.0 );
float hgt0 = gridlevelPRIME_to_height( dtx, 0.0 );
float diff = hgt1-hgt0;
if (fabs(diff) < 0.000001) diff = 0.000001;
dtx->Wscale[i] = 1.0 / (diff * 1000.0);
/* dtx->Wscale[i] = 1.0 / ((hgt1-hgt0) * 1000.0); */
}
else if (i==dtx->MaxNl-1) {
float hgt1 = gridlevelPRIME_to_height( dtx, (float)(dtx->MaxNl-1) );
float hgt0 = gridlevelPRIME_to_height( dtx, (float)(dtx->MaxNl-2) );
float diff = hgt1-hgt0;
if (fabs(diff) < 0.000001) diff = 0.000001;
dtx->Wscale[i] = 1.0 / (diff * 1000.0);
/* dtx->Wscale[i] = 1.0 / ((hgt1-hgt0) * 1000.0); */
}
else {
float a, b;
float hgt2 = gridlevelPRIME_to_height( dtx, (float)(i+1) );
float hgt1 = gridlevelPRIME_to_height( dtx, (float) i );
float hgt0 = gridlevelPRIME_to_height( dtx, (float)(i-1) );
float diffa = hgt1-hgt0;
float diffb = hgt2-hgt1;
if (fabs(diffa) < 0.000001) diffa = 0.000001;
if (fabs(diffb) < 0.000001) diffb = 0.000001;
a = 1.0 / (diffa * 1000.0);
b = 1.0 / (diffb * 1000.0);
/*
float a = 1.0 / ((hgt1-hgt0) * 1000.0);
float b = 1.0 / ((hgt2-hgt1) * 1000.0);
*/
dtx->Wscale[i] = (a+b) / 2.0;
}
}
break;
default:
printf("Error in init_trajPRIME: dtx->VerticalSystem=%d\n", dtx->VerticalSystem);
return 0;
}
return 1;
}
/*
* Given two timesteps and a (row,col,lev) position perform a 4-D
* interpolation of neighboring values to get the U,V, and W wind
* vector components in units of grid boxes/second.
* Input: t0, t1 - the two timesteps
* at, bt - timestep weights (at+bt=1.0)
* row, col, lev - the grid position
* Output: u, v, w - the interpolated and scaled wind vector values
* Return: 0 = missing data was found, 1 = valid values
*/
static int get_uvw( Context ctx, int t0, int t1,
float at, float bt,
float row, float col, float lev,
float *u, float *v, float *w, float sl )
{
int i, j, k;
float bra, ara, bca, aca, bla, ala, aa, ba, ab, bb;
float aaaa,abaa,baaa,bbaa,aaba,abba,baba,bbba,
aaab,abab,baab,bbab,aabb,abbb,babb,bbbb;
float ullll, ullhl, ulhll, ulhhl, uhlll, uhlhl, uhhll, uhhhl,
ulllh, ullhh, ulhlh, ulhhh, uhllh, uhlhh, uhhlh, uhhhh;
float vllll, vllhl, vlhll, vlhhl, vhlll, vhlhl, vhhll, vhhhl,
vlllh, vllhh, vlhlh, vlhhh, vhllh, vhlhh, vhhlh, vhhhh;
float wllll, wllhl, wlhll, wlhhl, whlll, whlhl, whhll, whhhl,
wlllh, wllhh, wlhlh, wlhhh, whllh, whlhh, whhlh, whhhh;
/* calculate weights for interpolating in box */
i = (int) row;
j = (int) col;
k = (int) lev;
bra = row-i;
ara = 1.0-bra;
bca = col-j;
aca = 1.0-bca;
bla = lev-k;
ala = 1.0-bla;
aa = ala*at;
ba = bla*at;
ab = ala*bt;
bb = bla*bt;
/* get 16 weights for 4-d box */
/* note 4 a/b's are in order of j,i,k,time */
aaaa = aca*ara*aa;
abaa = aca*bra*aa;
baaa = bca*ara*aa;
bbaa = bca*bra*aa;
aaba = aca*ara*ba;
abba = aca*bra*ba;
baba = bca*ara*ba;
bbba = bca*bra*ba;
aaab = aca*ara*ab;
abab = aca*bra*ab;
baab = bca*ara*ab;
bbab = bca*bra*ab;
aabb = aca*ara*bb;
abbb = aca*bra*bb;
babb = bca*ara*bb;
bbbb = bca*bra*bb;
/* get u,v,w grid values at 16 corners of 4-d box */
/* note 4 l/h's are in order of j,i,k,time */
/* return 0 if any missing data is found */
if (!sl){
if (!get_discrete_uvw( ctx, t0, i,j,k, &ullll, &vllll, &wllll )) return 0;
if (!get_discrete_uvw( ctx, t0, i,j,k+1, &ullhl, &vllhl, &wllhl )) return 0;
if (!get_discrete_uvw( ctx, t0, i+1,j,k, &ulhll, &vlhll, &wlhll )) return 0;
if (!get_discrete_uvw( ctx, t0, i+1,j,k+1, &ulhhl, &vlhhl, &wlhhl )) return 0;
if (!get_discrete_uvw( ctx, t0, i,j+1,k, &uhlll, &vhlll, &whlll )) return 0;
if (!get_discrete_uvw( ctx, t0, i,j+1,k+1, &uhlhl, &vhlhl, &whlhl )) return 0;
if (!get_discrete_uvw( ctx, t0, i+1,j+1,k, &uhhll, &vhhll, &whhll )) return 0;
if (!get_discrete_uvw( ctx, t0, i+1,j+1,k+1, &uhhhl, &vhhhl, &whhhl )) return 0;
if (!get_discrete_uvw( ctx, t1, i,j,k, &ulllh, &vlllh, &wlllh )) return 0;
if (!get_discrete_uvw( ctx, t1, i,j,k+1, &ullhh, &vllhh, &wllhh )) return 0;
if (!get_discrete_uvw( ctx, t1, i+1,j,k, &ulhlh, &vlhlh, &wlhlh )) return 0;
if (!get_discrete_uvw( ctx, t1, i+1,j,k+1, &ulhhh, &vlhhh, &wlhhh )) return 0;
if (!get_discrete_uvw( ctx, t1, i,j+1,k, &uhllh, &vhllh, &whllh )) return 0;
if (!get_discrete_uvw( ctx, t1, i,j+1,k+1, &uhlhh, &vhlhh, &whlhh )) return 0;
if (!get_discrete_uvw( ctx, t1, i+1,j+1,k, &uhhlh, &vhhlh, &whhlh )) return 0;
if (!get_discrete_uvw( ctx, t1, i+1,j+1,k+1, &uhhhh, &vhhhh, &whhhh )) return 0;
/* interpolate to get wind vector <U,V,W> */
*u = aaaa*ullll+aaab*ulllh+aaba*ullhl+aabb*ullhh+
abaa*ulhll+abab*ulhlh+abba*ulhhl+abbb*ulhhh+
baaa*uhlll+baab*uhllh+baba*uhlhl+babb*uhlhh+
bbaa*uhhll+bbab*uhhlh+bbba*uhhhl+bbbb*uhhhh;
*v = aaaa*vllll+aaab*vlllh+aaba*vllhl+aabb*vllhh+
abaa*vlhll+abab*vlhlh+abba*vlhhl+abbb*vlhhh+
baaa*vhlll+baab*vhllh+baba*vhlhl+babb*vhlhh+
bbaa*vhhll+bbab*vhhlh+bbba*vhhhl+bbbb*vhhhh;
*w = aaaa*wllll+aaab*wlllh+aaba*wllhl+aabb*wllhh+
abaa*wlhll+abab*wlhlh+abba*wlhhl+abbb*wlhhh+
baaa*whlll+baab*whllh+baba*whlhl+babb*whlhh+
bbaa*whhll+bbab*whhlh+bbba*whhhl+bbbb*whhhh;
}
else{
/* only one level */
if (!get_discrete_uv( ctx, t0, i,j,k, &ullll, &vllll, &wllll )) return 0;
if (!get_discrete_uv( ctx, t0, i,j,k, &ullhl, &vllhl, &wllhl )) return 0;
if (!get_discrete_uv( ctx, t0, i+1,j,k, &ulhll, &vlhll, &wlhll )) return 0;
if (!get_discrete_uv( ctx, t0, i+1,j,k, &ulhhl, &vlhhl, &wlhhl )) return 0;
if (!get_discrete_uv( ctx, t0, i,j+1,k, &uhlll, &vhlll, &whlll )) return 0;
if (!get_discrete_uv( ctx, t0, i,j+1,k, &uhlhl, &vhlhl, &whlhl )) return 0;
if (!get_discrete_uv( ctx, t0, i+1,j+1,k, &uhhll, &vhhll, &whhll )) return 0;
if (!get_discrete_uv( ctx, t0, i+1,j+1,k, &uhhhl, &vhhhl, &whhhl )) return 0;
if (!get_discrete_uv( ctx, t1, i,j,k, &ulllh, &vlllh, &wlllh )) return 0;
if (!get_discrete_uv( ctx, t1, i,j,k, &ullhh, &vllhh, &wllhh )) return 0;
if (!get_discrete_uv( ctx, t1, i+1,j,k, &ulhlh, &vlhlh, &wlhlh )) return 0;
if (!get_discrete_uv( ctx, t1, i+1,j,k, &ulhhh, &vlhhh, &wlhhh )) return 0;
if (!get_discrete_uv( ctx, t1, i,j+1,k, &uhllh, &vhllh, &whllh )) return 0;
if (!get_discrete_uv( ctx, t1, i,j+1,k, &uhlhh, &vhlhh, &whlhh )) return 0;
if (!get_discrete_uv( ctx, t1, i+1,j+1,k, &uhhlh, &vhhlh, &whhlh )) return 0;
if (!get_discrete_uv( ctx, t1, i+1,j+1,k, &uhhhh, &vhhhh, &whhhh )) return 0;
/* interpolate to get wind vector <U,V,W> */
*u = aaaa*ullll+aaab*ulllh+aaba*ullhl+aabb*ullhh+
abaa*ulhll+abab*ulhlh+abba*ulhhl+abbb*ulhhh+
baaa*uhlll+baab*uhllh+baba*uhlhl+babb*uhlhh+
bbaa*uhhll+bbab*uhhlh+bbba*uhhhl+bbbb*uhhhh;
*v = aaaa*vllll+aaab*vlllh+aaba*vllhl+aabb*vllhh+
abaa*vlhll+abab*vlhlh+abba*vlhhl+abbb*vlhhh+
baaa*vhlll+baab*vhllh+baba*vhlhl+babb*vhlhh+
bbaa*vhhll+bbab*vhhlh+bbba*vhhhl+bbbb*vhhhh;
*w = 0;
}
return 1;
}
/*
* Trace a wind trajectory.
* Input: row0, col0, lev0 - initial location
* time0 - initial timestep
* step - integration step size in seconds
* max - size of vr,vc,vl,vt arrays
* Output: vr, vc, vl - array of trajectory vertices
* vt - array of corresponding trajectory vertex times in seconds
* elapsed since first time step.
* Return: number of elements in vr,vc,vl,vt
*/
int trace( Context ctx, float row0, float col0, float lev0,
int time0, int step, int max,
float vr[], float vc[], float vl[], int vt[] )
{
float row, col, lev;
int time; /* in seconds since first time step */
int itime; /* time step corresponding to 'time' */
int i, vcount, singlelevel;
float maxr, maxc, maxl, minl;
/* grid bounds */
maxr = (float) (ctx->Nr-1);
maxc = (float) (ctx->Nc-1);
maxl = (float) (ctx->Nl[ctx->dpy_ctx->TrajU]-1);
minl = (float) ctx->Variable[ctx->dpy_ctx->TrajU]->LowLev;
/* allow trajs to work if u&v in one level only */
if (maxl == 0 && lev0 == minl){
singlelevel = 1;
}
else{
singlelevel = 0;
}
vcount = max;
/* TRACE TRAJECTORY BACKWARD (put vertices into array in backward order) */
row = row0;
col = col0;
lev = lev0;
time = ctx->Elapsed[time0];
itime = time0;
while (1) {
float at, bt, u, v, w;
/* test if current position is out of bounds */
if (row<0.0 || row>maxr || col<0.0 || col>maxc || lev<0.0 || lev>maxl || lev<minl)
break;
/* record point */
vcount--;
vr[vcount] = row;
vc[vcount] = col;
vl[vcount] = lev;
vt[vcount] = time;
if (vcount==0)
break;
/* test if time is out of bounds */
if (time<0)
break;
/* get u,v,w vector at current location and time */
if (itime+1==ctx->NumTimes) {
if (!get_uvw( ctx, itime, itime, 1.0, 0.0, row, col, lev, &u, &v, &w, singlelevel ))
break;
}
else {
if (ctx->Elapsed[itime]==ctx->Elapsed[itime+1]) {
/* this should never happen unless the file is screwy */
break;
}
at = (float) (ctx->Elapsed[itime+1] - time)
/ (float) (ctx->Elapsed[itime+1]-ctx->Elapsed[itime]);
bt = 1.0 - at;
if (!get_uvw( ctx, itime, itime+1, at, bt, row, col, lev, &u, &v, &w, singlelevel ))
break;
}
/* update position and time */
col -= u * (float) step; /* pos' = pos + velocity * time */
row -= v * (float) step;
lev -= w * (float) step;
time -= step;
if (time<ctx->Elapsed[itime])
itime--;
}
/* MOVE VERTICES FROM END TO BEGINNING OF ARRAY */
i = 0;
while (vcount<max) {
vr[i] = vr[vcount];
vc[i] = vc[vcount];
vl[i] = vl[vcount];
vt[i] = vt[vcount];
vcount++;
i++;
}
vcount = i;
/* TRACE TRAJECTORY FORWARD */
row = row0;
col = col0;
lev = lev0;
time = ctx->Elapsed[time0];
itime = time0;
while (1) {
float at, bt, u, v, w;
/* test if position is out of bounds */
if (row<0.0 || row>maxr || col<0.0 || col>maxc || lev<0.0 || lev>maxl || lev<minl)
break;
/* record point */
vr[vcount] = row;
vc[vcount] = col;
vl[vcount] = lev;
vt[vcount] = time;
vcount++;
if (vcount>=max)
break;
/* test if time is out of bounds */
if (time>=ctx->Elapsed[ctx->NumTimes-1])
break;
/* get u,v,w vector at current location and time */
if (ctx->Elapsed[itime]==ctx->Elapsed[itime+1]) {
/* this should never happen unless the file is screwy */
break;
}
at = (float) (ctx->Elapsed[itime+1] - time)
/ (float) (ctx->Elapsed[itime+1]-ctx->Elapsed[itime]);
bt = 1.0 - at;
if (!get_uvw( ctx, itime, itime+1, at, bt, row, col, lev, &u, &v, &w, singlelevel ))
break;
/* update position and time */
col += u * (float) step; /* pos' = pos + velocity * time */
row += v * (float) step;
lev += w * (float) step;
time += step;
if (time>ctx->Elapsed[itime+1])
itime++;
}
/* print the vertex list */
/*for (i=0;i<vcount;i++) {
printf("%3d: %6.3f %6.3f %6.3f %d\n", i, vr[i], vc[i], vl[i], vt[i] );
}*/
if (vcount > max){
vcount = max;
}
return vcount;
}
#define MAX 5000
#define EPS 0.0000000001
#define WSCALE 150.0
#define RADIUS 10.0
/*
* Convert a trajectory path to a ribbon.
*/
int to_ribbon( int num, float xt[], float yt[], float zt[],
int itint[], float xn[], float yn[], float zn[] )
{
int i, it, ip, il, ih;
int iuu, imm, idd, iu, id;
float ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez;
float cc, dd, ee;
float scale;
float xr[MAX], yr[MAX], zr[MAX];
int itinr[MAX];
scale = 2.0/WSCALE;
for (i=0;i<num;i++) {
xr[i] = xt[i];
yr[i] = yt[i];
zr[i] = zt[i];
itinr[i] = itint[i];
}
/* DO 300 ITRAJ=ITRLOW,ITRHI */
/* IL = LSTART(ITRAJ) */
/* IH = LEND(ITRAJ) */
/* CALCULATE UNFILTERED RIBBON ORIENTATIONS */
il = 0;
ih = num;
for (i=il;i<ih;i++) {
/* DO 10 I=IL,IH */
it = 2*i - 0;
ip = it+1;
/* SHORT TRAJECTORY CASE */
if (ih-il <= 2) {
xt[it] = xr[i];
yt[it] = yr[i];
zt[it] = zr[i];
xn[it] = 1.0;
yn[it] = 0.0;
zn[it] = 0.0;
itint[it] = 0;
xt[ip] = xr[i];
yt[ip] = yr[i];
zt[ip] = zr[i];
xn[ip] = 1.0;
yn[ip] = 0.0;
zn[ip] = 0.0;
itint[ip] = itinr[i];
/* go to 10 */
}
else {
/* LONG TRAJECTORY CASE */
/* RAW CROSS RIBBON VECTORS */
iuu = i+1;
imm = i;
idd = i-1;
if (i == il) {
iuu = i+2;
imm = i+1;
idd = i;
}
if (i == ih-1) {
iuu = i;
imm = i-1;
idd = i-2;
}
ax = xr[iuu]-xr[imm];
ay = yr[iuu]-yr[imm];
az = zr[iuu]-zr[imm];
bx = xr[imm]-xr[idd];
by = yr[imm]-yr[idd];
bz = zr[imm]-zr[idd];
cx = ay*bz-az*by;
cy = az*bx-ax*bz;
cz = ax*by-ay*bx;
cc = sqrt(cx*cx+cy*cy+cz*cz);
if (cc < EPS)
cc = EPS;
cc = 1.0/cc;
xn[ip] = cx*cc;
yn[ip] = cy*cc;
zn[ip] = cz*cc;
/* ALONG RIBBON VECTOR */
iu = i+1;
id = i-1;
if (i == il) id = i;
if (i == ih-1) iu = i;
dx = xr[iu]-xr[id];
dy = yr[iu]-yr[id];
dz = zr[iu]-zr[id];
dd = sqrt(dx*dx+dy*dy+dz*dz);
if (dd < EPS)
dd = EPS;
dd = 1.0/dd;
xn[it] = dx*dd;
yn[it] = dy*dd;
zn[it] = dz*dd;
}
}
/* 10 CONTINUE */
if (ih-il < 2) /* <= 2 ??? */
return 0;
/* PROPOGATE CROSS RIBBON VECTOR */
cx = 0.0;
cy = 0.0;
cz = 0.0;
/* DO 20 I=IL,IH */
for (i=il;i<ih;i++) {
it = 2*i-0;
ip = it+1;
/* ADD (OR SUBTRACT) RAW CROSS RIBBON FROM RUNNING SUM */
ax=xn[ip];
ay=yn[ip];
az=zn[ip];
if (ax*cx+ay*cy+az*cz > 0.0) {
cx=cx+ax;
cy=cy+ay;
cz=cz+az;
}
else {
cx=cx-ax;
cy=cy-ay;
cz=cz-az;
}
/* GET COMPONENT OF (CROSS RIBBON) NORMAL TO ALONG RIBBON */
dx=xn[it];
dy=yn[it];
dz=zn[it];
dd=cx*dx+cy*dy+cz*dz;
cx=cx-dd*dx;
cy=cy-dd*dy;
cz=cz-dd*dz;
/* SCALE CROSS RIBBON TO RADIUS LENGTH */
cc = sqrt(cx*cx+cy*cy+cz*cz);
if (cc < EPS)
cc = EPS;
cc = RADIUS/cc;
cx=cx*cc;
cy=cy*cc;
cz=cz*cc;
xt[it]=cx;
yt[it]=cy;
zt[it]=cz;
}
/* 20 CONTINUE; */
/* NOW DO RIBBONS */
/* DO 40 I=IL,IH */
for (i=il;i<ih;i++) {
it = 2*i - 0;
ip = it+1;
cx=xt[it];
cy=yt[it];
cz=zt[it];
dx=xn[it];
dy=yn[it];
dz=zn[it];
cc = scale/RADIUS;
cx=cx*cc;
cy=cy*cc;
cz=cz*cc;
/* CALCULATE NORMAL TO RIBBON */
ex = cy*dz - cz*dy;
ey = cz*dx - cx*dz;
ez = cx*dy - cy*dx;
ee = sqrt(ex*ex+ey*ey+ez*ez);
if (ee < EPS) ee = EPS;
ee = 1.0/ee;
ex=ex*ee;
ey=ey*ee;
ez=ez*ee;
/* MAKE RIBBON POINTS */
xt[it]=xr[i] - cx;
yt[it]=yr[i] - cy;
zt[it]=zr[i] - cz;
xn[it]=ex;
yn[it]=ey;
zn[it]=ez;
itint[it]=0;
xt[ip]=xr[i] + cx;
yt[ip]=yr[i] + cy;
zt[ip]=zr[i] + cz;
xn[ip]=ex;
yn[ip]=ey;
zn[ip]=ez;
itint[ip]=itinr[i];
}
/* 40 CONTINUE */
return num*2;
}
syntax highlighted by Code2HTML, v. 0.9.1