Unverified Commit 2aeadff3 authored by tomrink's avatar tomrink Committed by GitHub
Browse files

Merge pull request #146 from tomrink/master

Runge-Kutta (RK4) determination of parcel trajectory path
parents cb9985bf ad472e46
......@@ -189,6 +189,232 @@ public class Trajectory {
}
public void forwardRK4(FlowInfo info, float[][] flow_values, float[][] flow_valuesB, float[][] flow_valuesC, byte[][] color_values, GriddedSet spatial_set, FlatField terrain, int direction, float timeStep)
throws VisADException, RemoteException {
if (offGrid) return;
clrDim = color_values.length;
float[] intrpClr = new float[clrDim];
this.spatial_set = spatial_set;
this.manifoldDimension = spatial_set.getManifoldDimension();
float[][] flowLoc = new float[3][1];
float[][] flowVec = new float[3][1];
float[] intrpFlow = new float[3];
float[] k1 = new float[3];
float[] k2 = new float[3];
float[] k3 = new float[3];
float[] k4 = new float[3];
float[] xyz = new float[3];
if (indices[0] != null) {
java.util.Arrays.fill(intrpFlow, 0f);
java.util.Arrays.fill(intrpClr, 0);
// k1 -------------------------------------------
k1[0] = 0;
k1[1] = 0;
k1[2] = 0;
for (int j=0; j<indices[0].length; j++) {
int idx = indices[0][j];
flowLoc[0][0] = info.spatial_values[0][idx];
flowLoc[1][0] = info.spatial_values[1][idx];
flowLoc[2][0] = info.spatial_values[2][idx];
flowVec[0][0] = flow_values[0][idx];
flowVec[1][0] = flow_values[1][idx];
flowVec[2][0] = flow_values[2][idx];
float[][] del = TrajectoryManager.computeDisplacement(info, flowLoc, flowVec, timeStep);
k1[0] += weights[0][j]*(direction)*del[0][0];
k1[1] += weights[0][j]*(direction)*del[1][0];
k1[2] += weights[0][j]*(direction)*del[2][0];
}
// k2 ----------------------------------
xyz[0] = startPts[0] + k1[0]/2;
xyz[1] = startPts[1] + k1[1]/2;
xyz[2] = startPts[2] + k1[2]/2;
int[][] indicesK = new int[1][];
float[][] weightsK = new float[1][];
if (manifoldDimension == 2) {
float[][] pts2D = new float[2][1];
pts2D[0][0] = xyz[0];
pts2D[1][0] = xyz[1];
spatial_set.valueToInterp(pts2D, indicesK, weightsK, guess2D);
}
else if (manifoldDimension == 3) {
float[][] pts3D = new float[3][1];
pts3D[0][0] = xyz[0];
pts3D[1][0] = xyz[1];
pts3D[2][0] = xyz[2];
spatial_set.valueToInterp(pts3D, indicesK, weightsK, guess3D);
}
if (indicesK[0] != null) {
k2[0] = 0;
k2[1] = 0;
k2[2] = 0;
for (int j=0; j<indicesK[0].length; j++) {
int idx = indicesK[0][j];
flowLoc[0][0] = info.spatial_values[0][idx];
flowLoc[1][0] = info.spatial_values[1][idx];
flowLoc[2][0] = info.spatial_values[2][idx];
flowVec[0][0] = flow_valuesB[0][idx];
flowVec[1][0] = flow_valuesB[1][idx];
flowVec[2][0] = flow_valuesB[2][idx];
float[][] del = TrajectoryManager.computeDisplacement(info, flowLoc, flowVec, timeStep);
k2[0] += weightsK[0][j]*(direction)*del[0][0];
k2[1] += weightsK[0][j]*(direction)*del[1][0];
k2[2] += weightsK[0][j]*(direction)*del[2][0];
}
}
// k3 -----------------------------------
xyz[0] = startPts[0] + k2[0]/2;
xyz[1] = startPts[1] + k2[1]/2;
xyz[2] = startPts[2] + k2[2]/2;
indicesK = new int[1][];
weightsK = new float[1][];
if (manifoldDimension == 2) {
float[][] pts2D = new float[2][1];
pts2D[0][0] = xyz[0];
pts2D[1][0] = xyz[1];
spatial_set.valueToInterp(pts2D, indicesK, weightsK, guess2D);
}
else if (manifoldDimension == 3) {
float[][] pts3D = new float[3][1];
pts3D[0][0] = xyz[0];
pts3D[1][0] = xyz[1];
pts3D[2][0] = xyz[2];
spatial_set.valueToInterp(pts3D, indicesK, weightsK, guess3D);
}
if (indicesK[0] != null) {
k3[0] = 0;
k3[1] = 0;
k3[2] = 0;
for (int j=0; j<indicesK[0].length; j++) {
int idx = indicesK[0][j];
flowLoc[0][0] = info.spatial_values[0][idx];
flowLoc[1][0] = info.spatial_values[1][idx];
flowLoc[2][0] = info.spatial_values[2][idx];
flowVec[0][0] = flow_valuesB[0][idx];
flowVec[1][0] = flow_valuesB[1][idx];
flowVec[2][0] = flow_valuesB[2][idx];
float[][] del = TrajectoryManager.computeDisplacement(info, flowLoc, flowVec, timeStep);
k3[0] += weightsK[0][j]*(direction)*del[0][0];
k3[1] += weightsK[0][j]*(direction)*del[1][0];
k3[2] += weightsK[0][j]*(direction)*del[2][0];
}
}
// k4 ----------------------------------------
xyz[0] = startPts[0] + k3[0];
xyz[1] = startPts[1] + k3[1];
xyz[2] = startPts[2] + k3[2];
indicesK = new int[1][];
weightsK = new float[1][];
if (manifoldDimension == 2) {
float[][] pts2D = new float[2][1];
pts2D[0][0] = xyz[0];
pts2D[1][0] = xyz[1];
spatial_set.valueToInterp(pts2D, indicesK, weightsK, guess2D);
}
else if (manifoldDimension == 3) {
float[][] pts3D = new float[3][1];
pts3D[0][0] = xyz[0];
pts3D[1][0] = xyz[1];
pts3D[2][0] = xyz[2];
spatial_set.valueToInterp(pts3D, indicesK, weightsK, guess3D);
}
if (indicesK[0] != null) {
k4[0] = 0;
k4[1] = 0;
k4[2] = 0;
for (int j=0; j<indicesK[0].length; j++) {
int idx = indicesK[0][j];
flowLoc[0][0] = info.spatial_values[0][idx];
flowLoc[1][0] = info.spatial_values[1][idx];
flowLoc[2][0] = info.spatial_values[2][idx];
flowVec[0][0] = flow_valuesC[0][idx];
flowVec[1][0] = flow_valuesC[1][idx];
flowVec[2][0] = flow_valuesC[2][idx];
float[][] del = TrajectoryManager.computeDisplacement(info, flowLoc, flowVec, timeStep);
k4[0] += weightsK[0][j]*(direction)*del[0][0];
k4[1] += weightsK[0][j]*(direction)*del[1][0];
k4[2] += weightsK[0][j]*(direction)*del[2][0];
}
}
stopPts[0] = startPts[0] + (1f/6)*k1[0] + (1f/3)*k2[0] + (1f/3)*k3[0] + (1f/6)*k4[0];
stopPts[1] = startPts[1] + (1f/6)*k1[1] + (1f/3)*k2[1] + (1f/3)*k3[1] + (1f/6)*k4[1];
stopPts[2] = startPts[2] + (1f/6)*k1[2] + (1f/3)*k2[2] + (1f/3)*k3[2] + (1f/6)*k4[2];
xyz[0] = stopPts[0];
xyz[1] = stopPts[1];
xyz[2] = stopPts[2];
if (manifoldDimension == 2) {
float[][] pts2D = new float[2][1];
pts2D[0][0] = xyz[0];
pts2D[1][0] = xyz[1];
spatial_set.valueToInterp(pts2D, indices, weights, guess2D);
}
else if (manifoldDimension == 3) {
float[][] pts3D = new float[3][1];
pts3D[0][0] = xyz[0];
pts3D[1][0] = xyz[1];
pts3D[2][0] = xyz[2];
spatial_set.valueToInterp(pts3D, indices, weights, guess3D);
if (terrain != null) {
adjustFlowAtTerrain(terrain, color_values);
}
}
if (indices[0] != null) {
for (int j=0; j<indices[0].length; j++) {
int idx = indices[0][j];
intrpClr[0] += weights[0][j]*ShadowType.byteToFloat(color_values[0][idx]);
intrpClr[1] += weights[0][j]*ShadowType.byteToFloat(color_values[1][idx]);
intrpClr[2] += weights[0][j]*ShadowType.byteToFloat(color_values[2][idx]);
if (clrDim == 4) {
intrpClr[3] += weights[0][j]*ShadowType.byteToFloat(color_values[3][idx]);
}
}
stopColor[0] = ShadowType.floatToByte(intrpClr[0]);
stopColor[1] = ShadowType.floatToByte(intrpClr[1]);
stopColor[2] = ShadowType.floatToByte(intrpClr[2]);
if (clrDim == 4) {
stopColor[3] = ShadowType.floatToByte(intrpClr[3]);
}
}
addPair(startPts, stopPts, startColor, stopColor);
uVecPath[0] = stopPts[0] - startPts[0];
uVecPath[1] = stopPts[1] - startPts[1];
uVecPath[2] = stopPts[2] - startPts[2];
float mag = (float) Math.sqrt(uVecPath[0]*uVecPath[0] + uVecPath[1]*uVecPath[1] + uVecPath[2]*uVecPath[2]);
uVecPath[0] /= mag;
uVecPath[1] /= mag;
uVecPath[2] /= mag;
update();
}
}
private void addPair(float[] startPt, float[] stopPt, byte[] startColor, byte[] stopColor) {
indexes[npairs] = trajMan.getCoordinateCount();
......
......@@ -87,13 +87,15 @@ public class TrajectoryManager {
boolean manualIntrpPts;
boolean trajDoIntrp = true;
boolean trajCachingEnabled = false;
boolean doHysplit = true;
//boolean doHysplit = false;
//boolean doRK4 = true;
float trcrSize = 1f;
boolean trcrEnabled;
boolean terrainFollowEnabled;
int numIntrpPts;
int trajSkip;
TrajectoryParams.SmoothParams smoothParams;
TrajectoryParams.Method method;
int direction;
int trajForm = LINE; // Default
float cylWidth = 0.01f;
......@@ -107,6 +109,9 @@ public class TrajectoryManager {
float[] intrpU_1;
float[] intrpV_1;
float[] intrpW_1;
float[] intrpU_2;
float[] intrpV_2;
float[] intrpW_2;
CubicInterpolator uInterp;
CubicInterpolator vInterp;
CubicInterpolator wInterp;
......@@ -148,6 +153,7 @@ public class TrajectoryManager {
numIntrpPts = trajParams.getNumIntrpPts();
trajSkip = trajParams.getStartSkip();
smoothParams = trajParams.getSmoothParams();
method = trajParams.getMethod();
direction = trajParams.getDirection();
startPts = trajParams.getStartPoints();
trajDoIntrp = trajParams.getDoIntrp();
......@@ -251,11 +257,20 @@ public class TrajectoryManager {
intrpV = new float[numSpatialPts];
intrpW = new float[numSpatialPts];
if (doHysplit) {
if (method == TrajectoryParams.Method.HySplit) {
intrpU_1 = new float[numSpatialPts];
intrpV_1 = new float[numSpatialPts];
intrpW_1 = new float[numSpatialPts];
}
else if (method == TrajectoryParams.Method.RK4) {
intrpU_1 = new float[numSpatialPts];
intrpV_1 = new float[numSpatialPts];
intrpW_1 = new float[numSpatialPts];
intrpU_2 = new float[numSpatialPts];
intrpV_2 = new float[numSpatialPts];
intrpW_2 = new float[numSpatialPts];
}
uInterp = new CubicInterpolator(trajDoIntrp, numSpatialPts);
vInterp = new CubicInterpolator(trajDoIntrp, numSpatialPts);
......@@ -362,6 +377,9 @@ public class TrajectoryManager {
}
float timeStep = (float) timeSteps[i]/numIntrpPts;
if (!trajDoIntrp && (method == TrajectoryParams.Method.RK4)) {
timeStep *= 2;
}
if ((timeAccum >= trajRefreshInterval)) { // for non steady state trajectories (refresh frequency)
trajectories = new ArrayList<Trajectory>();
......@@ -424,7 +442,7 @@ public class TrajectoryManager {
int numTrajectories = trajectories.size();
reset();
for (int ti=0; ti<numIntrpPts; ti++) { // additional points per domain time step
double dst = (x1 - x0)/numIntrpPts;
double xt = x0 + dst*ti;
......@@ -434,30 +452,50 @@ public class TrajectoryManager {
uInterp.interpolate(xt, intrpU);
vInterp.interpolate(xt, intrpV);
wInterp.interpolate(xt, intrpW);
if (doHysplit) { // NOAA HySplit
if (ti == numIntrpPts-1) {
System.arraycopy(values1[0], 0, intrpU_1, 0, intrpU_1.length);
System.arraycopy(values1[1], 0, intrpV_1, 0, intrpV_1.length);
System.arraycopy(values1[2], 0, intrpW_1, 0, intrpW_1.length);
}
else {
if (method == TrajectoryParams.Method.RK4) { // Runge-Kutta
double step = dst;
if (!trajDoIntrp) {
step = 2*dst;
}
if (k == dataDomainLength-2 && ti > numIntrpPts-2) {
continue;
}
uInterp.interpolate(xt+step/2, intrpU_1);
vInterp.interpolate(xt+step/2, intrpV_1);
wInterp.interpolate(xt+step/2, intrpW_1);
uInterp.interpolate(xt+step, intrpU_2);
vInterp.interpolate(xt+step, intrpV_2);
wInterp.interpolate(xt+step, intrpW_2);
for (int t=0; t<numTrajectories; t++) {
Trajectory traj = trajectories.get(t);
traj.currentTimeIndex = direction*i;
traj.currentTime = direction*times[i];
traj.forwardRK4(info, new float[][] {intrpU, intrpV, intrpW},
new float[][] {intrpU_1, intrpV_1, intrpW_1},
new float[][] {intrpU_2, intrpV_2, intrpW_2},
color_values, spatialSetTraj, terrain, direction, timeStep);
}
}
else {
if (method == TrajectoryParams.Method.HySplit) { // NOAA HySplit
uInterp.interpolate(xt+dst, intrpU_1);
vInterp.interpolate(xt+dst, intrpV_1);
wInterp.interpolate(xt+dst, intrpW_1);
}
intrpU = mean(intrpU, intrpU_1);
intrpV = mean(intrpV, intrpV_1);
intrpW = mean(intrpW, intrpW_1);
}
intrpU = mean(intrpU, intrpU_1);
intrpV = mean(intrpV, intrpV_1);
intrpW = mean(intrpW, intrpW_1);
}
for (int t=0; t<numTrajectories; t++) {
Trajectory traj = trajectories.get(t);
traj.currentTimeIndex = direction*i;
traj.currentTime = direction*times[i];
traj.forward(info, new float[][] {intrpU, intrpV, intrpW}, color_values, spatialSetTraj, terrain, direction, timeStep);
for (int t=0; t<numTrajectories; t++) {
Trajectory traj = trajectories.get(t);
traj.currentTimeIndex = direction*i;
traj.currentTime = direction*times[i];
traj.forward(info, new float[][] {intrpU, intrpV, intrpW}, color_values, spatialSetTraj, terrain, direction, timeStep);
}
}
} // inner time loop (time interpolation)
......@@ -2328,43 +2366,49 @@ public class TrajectoryManager {
String propStr = null;
propStr = prop.getProperty("CylinderWidthFactor");
if (propStr != null) {
float fac = Float.valueOf(propStr);
float fac = Float.valueOf(propStr.trim());
trajParams.setCylinderWidth(trajParams.getCylinderWidth()*fac);
}
propStr = prop.getProperty("RibbonWidthFactor");
if (propStr != null) {
float fac = Float.valueOf(propStr);
float fac = Float.valueOf(propStr.trim());
trajParams.setRibbonWidthFactor(fac);
}
propStr = prop.getProperty("ManualIntrpPts");
if (propStr != null) {
trajParams.setManualIntrpPts(Boolean.valueOf(propStr));
trajParams.setManualIntrpPts(Boolean.valueOf(propStr.trim()));
}
propStr = prop.getProperty("TrajDoIntrp");
if (propStr != null) {
trajParams.setDoIntrp(Boolean.valueOf(propStr));
trajParams.setDoIntrp(Boolean.valueOf(propStr.trim()));
}
propStr = prop.getProperty("TerrainFollow");
if (propStr != null) {
trajParams.setTerrainFollowing(Boolean.valueOf(propStr));
trajParams.setTerrainFollowing(Boolean.valueOf(propStr.trim()));
}
propStr = prop.getProperty("NumIntrpPts");
if (propStr != null) {
trajParams.setNumIntrpPts(Integer.valueOf(propStr));
trajParams.setNumIntrpPts(Integer.valueOf(propStr.trim()));
}
propStr = prop.getProperty("TrajRefreshInterval");
if (propStr != null) {
trajParams.setTrajRefreshInterval(Double.valueOf(propStr));
trajParams.setTrajRefreshInterval(Double.valueOf(propStr.trim()));
}
propStr = prop.getProperty("StartSkip");
if (propStr != null) {
trajParams.setStartSkip(Integer.valueOf(propStr.trim()));
}
propStr = prop.getProperty("TrajForm");
if (propStr != null) {
propStr = propStr.trim();
if (propStr.equals("LINE")) {
trajParams.setTrajectoryForm(TrajectoryParams.LINE);
}
......@@ -2381,6 +2425,20 @@ public class TrajectoryManager {
trajParams.setTrajectoryForm(TrajectoryParams.POINT);
}
}
propStr = prop.getProperty("Method");
if (propStr != null) {
propStr = propStr.trim();
if (propStr.equals("HYSPLIT")) {
trajParams.setMethod(TrajectoryParams.Method.HySplit);
}
else if (propStr.equals("RK4")) {
trajParams.setMethod(TrajectoryParams.Method.RK4);
}
else if (propStr.equals("EULER")) {
trajParams.setMethod(TrajectoryParams.Method.Euler);
}
}
is.close();
}
......
......@@ -48,6 +48,12 @@ public class TrajectoryParams {
}
}
public static enum Method {
HySplit,
RK4,
Euler;
}
public static final int LINE = 0;
public static final int RIBBON = 1;
public static final int CYLINDER = 2;
......@@ -81,6 +87,8 @@ public class TrajectoryParams {
// terrain (lower boundary) Implicit: meters above MSL
FlatField terrain = null;
Method method = Method.HySplit; //Default
public TrajectoryParams() {
}
......@@ -108,6 +116,7 @@ public class TrajectoryParams {
this.zStartSkip = params.getZStartSkip();
this.terrain = params.getTerrain();
this.terrainFollowEnabled = params.getTerrainFollowing();
this.method = params.getMethod();
}
public TrajectoryParams(double trajVisibilityTimeWindow, double trajRefreshInterval, int numIntrpPts, int startSkip, SmoothParams smoothParams) {
......@@ -211,6 +220,10 @@ public class TrajectoryParams {
this.markerEnabled = yesno;
}
public void setMethod(Method method) {
this.method = method;
}
public void setCachingEnabled(boolean yesno) {
this.cachingEnabled = yesno;
}
......@@ -291,6 +304,10 @@ public class TrajectoryParams {
return terrainFollowEnabled;
}
public Method getMethod() {
return method;
}
public void setStartPoints(float[][] startPts) {
this.startPoints = startPts;
this.startType = Display.DisplaySpatialCartesianTuple;
......@@ -380,6 +397,9 @@ public class TrajectoryParams {
else if (this.terrainFollowEnabled != trajParams.terrainFollowEnabled) {
return false;
}
else if (this.method != trajParams.method) {
return false;
}
}
return true;
}
......
......@@ -1427,7 +1427,15 @@ System.out.println("Texture.BASE_LEVEL_LINEAR = " + Texture.BASE_LEVEL_LINEAR);
VisADGeometryArray[] auxArray = new VisADGeometryArray[2];
ArrayList<float[]> achrArrays = null;
for (int k=0; k<dataDomainLength; k++) {
// We don't really have a time interval defined beyond the next to last point
int computeLength = dataDomainLength-1;
if (!trajParams.getDoIntrp()) {
if (trajParams.getMethod() == TrajectoryParams.Method.RK4) { // need 3 time steps
computeLength = dataDomainLength-2;
}
}
for (int k=0; k<computeLength; k++) {
int i = (direction < 0) ? ((dataDomainLength-1) - k) : k;
FlowInfo info = flowInfoList.get(i);
......
......@@ -6,28 +6,28 @@ import java.util.Arrays;
public class CubicInterpolator {
LUDecomposition solver;
private LUDecomposition solver;
double[][] solution = null;
private double[][] solution = null;
double x0 = 0;
double x1 = 0;
double x2 = 0;
double x0_last = 0;
double x0_save;
private double x0 = 0;
private double x1 = 0;
private double x2 = 0;
private double x0_last = 0;
private double x0_save;
float[] values0 = null;
float[] values1 = null;
float[] values2 = null;
float[] values0_last = null;
float[] values0_save = null;
private float[] values0 = null;
private float[] values1 = null;
private float[] values2 = null;
private float[] values0_last = null;
private float[] values0_save = null;
int numSpatialPts = 1;
private int numSpatialPts = 1;
boolean doIntrp = true;
private boolean doIntrp = true;
boolean[] needed = null;
boolean[] computed = null;
private boolean[] needed = null;
private boolean[] computed = null;
public CubicInterpolator(boolean doIntrp, int numSpatialPts) {
this.doIntrp = doIntrp;
......@@ -63,6 +63,9 @@ public class CubicInterpolator {
else if (xt == x1) {
System.arraycopy(values1, 0, interpValues, 0, numSpatialPts);
}
else if (xt == x2) {