xan/awips2
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge "Issue #2289 Speed up cape calculations. Change-Id: I56f2373779b2973a145b62c4f47b278458fd2d47" into development

Browse source 
Former-commit-id: 7dc787bc00 [formerly 397b5123cb] [formerly 7dc787bc00 [formerly 397b5123cb] [formerly 91a5f2b669 [formerly 53faa91d53be5f71f6d0bee2174ec50fc41405a9]]]
Former-commit-id: 91a5f2b669
Former-commit-id: 7a71c50926 [formerly ded4514e9d]
Former-commit-id: 1a3e6b7501
This commit is contained in:
Richard Peter 2013-08-23 11:23:25 -05:00 • committed by Gerrit Code Review
commit 6189645b1d
4 changed files with 227 additions and 87 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

63
cave/com.raytheon.uf.viz.derivparam/localization/derivedParameters/definitions/cCape.xml
View file

@ -72,4 +72,67 @@
<Method name="Max" levels="MaxCape">
<Field abbreviation="cCape" level="3D" />
</Method>
<!-- For Composite Agl levels it is better to explicitly list the levels
to Average to avoid pulling in extra levels, specifically levels defined at
kft heights. -->
<Method name="Average" levels="0-1kmAgl">
<Field abbreviation="cCape" level="0kmAgl" />
<Field abbreviation="cCape" level="0.5kmAgl" />
<Field abbreviation="cCape" level="1kmAgl" />
</Method>
<Method name="Average" levels="0-2kmAgl">
<Field abbreviation="cCape" level="0kmAgl" />
<Field abbreviation="cCape" level="0.5kmAgl" />
<Field abbreviation="cCape" level="1kmAgl" />
<Field abbreviation="cCape" level="1.5kmAgl" />
<Field abbreviation="cCape" level="2kmAgl" />
</Method>
<Method name="Average" levels="0-3kmAgl">
<Field abbreviation="cCape" level="0kmAgl" />
<Field abbreviation="cCape" level="0.5kmAgl" />
<Field abbreviation="cCape" level="1kmAgl" />
<Field abbreviation="cCape" level="1.5kmAgl" />
<Field abbreviation="cCape" level="2kmAgl" />
<Field abbreviation="cCape" level="2.5kmAgl" />
<Field abbreviation="cCape" level="3kmAgl" />
</Method>
<Method name="Average" levels="0-4kmAgl">
<Field abbreviation="cCape" level="0kmAgl" />
<Field abbreviation="cCape" level="0.5kmAgl" />
<Field abbreviation="cCape" level="1kmAgl" />
<Field abbreviation="cCape" level="1.5kmAgl" />
<Field abbreviation="cCape" level="2kmAgl" />
<Field abbreviation="cCape" level="2.5kmAgl" />
<Field abbreviation="cCape" level="3kmAgl" />
<Field abbreviation="cCape" level="3.5kmAgl" />
<Field abbreviation="cCape" level="4kmAgl" />
</Method>
<Method name="Average" levels="0-5kmAgl">
<Field abbreviation="cCape" level="0kmAgl" />
<Field abbreviation="cCape" level="0.5kmAgl" />
<Field abbreviation="cCape" level="1kmAgl" />
<Field abbreviation="cCape" level="1.5kmAgl" />
<Field abbreviation="cCape" level="2kmAgl" />
<Field abbreviation="cCape" level="2.5kmAgl" />
<Field abbreviation="cCape" level="3kmAgl" />
<Field abbreviation="cCape" level="3.5kmAgl" />
<Field abbreviation="cCape" level="4kmAgl" />
<Field abbreviation="cCape" level="4.5kmAgl" />
<Field abbreviation="cCape" level="5kmAgl" />
</Method>
<Method name="Average" levels="0-6kmAgl">
<Field abbreviation="cCape" level="0kmAgl" />
<Field abbreviation="cCape" level="0.5kmAgl" />
<Field abbreviation="cCape" level="1kmAgl" />
<Field abbreviation="cCape" level="1.5kmAgl" />
<Field abbreviation="cCape" level="2kmAgl" />
<Field abbreviation="cCape" level="2.5kmAgl" />
<Field abbreviation="cCape" level="3kmAgl" />
<Field abbreviation="cCape" level="3.5kmAgl" />
<Field abbreviation="cCape" level="4kmAgl" />
<Field abbreviation="cCape" level="4.5kmAgl" />
<Field abbreviation="cCape" level="5kmAgl" />
<Field abbreviation="cCape" level="5.5kmAgl" />
<Field abbreviation="cCape" level="6kmAgl" />
</Method>
</DerivedParameter>

63
cave/com.raytheon.uf.viz.derivparam/localization/derivedParameters/definitions/cCin.xml
View file

@ -54,4 +54,67 @@
<Field abbreviation="SHx"/>
<ConstantField value="0.0"/>
</Method>
<!-- For Composite Agl levels it is better to explicitly list the levels
to Average to avoid pulling in extra levels, specifically levels defined at
kft heights. -->
<Method name="Average" levels="0-1kmAgl">
<Field abbreviation="cCin" level="0kmAgl" />
<Field abbreviation="cCin" level="0.5kmAgl" />
<Field abbreviation="cCin" level="1kmAgl" />
</Method>
<Method name="Average" levels="0-2kmAgl">
<Field abbreviation="cCin" level="0kmAgl" />
<Field abbreviation="cCin" level="0.5kmAgl" />
<Field abbreviation="cCin" level="1kmAgl" />
<Field abbreviation="cCin" level="1.5kmAgl" />
<Field abbreviation="cCin" level="2kmAgl" />
</Method>
<Method name="Average" levels="0-3kmAgl">
<Field abbreviation="cCin" level="0kmAgl" />
<Field abbreviation="cCin" level="0.5kmAgl" />
<Field abbreviation="cCin" level="1kmAgl" />
<Field abbreviation="cCin" level="1.5kmAgl" />
<Field abbreviation="cCin" level="2kmAgl" />
<Field abbreviation="cCin" level="2.5kmAgl" />
<Field abbreviation="cCin" level="3kmAgl" />
</Method>
<Method name="Average" levels="0-4kmAgl">
<Field abbreviation="cCin" level="0kmAgl" />
<Field abbreviation="cCin" level="0.5kmAgl" />
<Field abbreviation="cCin" level="1kmAgl" />
<Field abbreviation="cCin" level="1.5kmAgl" />
<Field abbreviation="cCin" level="2kmAgl" />
<Field abbreviation="cCin" level="2.5kmAgl" />
<Field abbreviation="cCin" level="3kmAgl" />
<Field abbreviation="cCin" level="3.5kmAgl" />
<Field abbreviation="cCin" level="4kmAgl" />
</Method>
<Method name="Average" levels="0-5kmAgl">
<Field abbreviation="cCin" level="0kmAgl" />
<Field abbreviation="cCin" level="0.5kmAgl" />
<Field abbreviation="cCin" level="1kmAgl" />
<Field abbreviation="cCin" level="1.5kmAgl" />
<Field abbreviation="cCin" level="2kmAgl" />
<Field abbreviation="cCin" level="2.5kmAgl" />
<Field abbreviation="cCin" level="3kmAgl" />
<Field abbreviation="cCin" level="3.5kmAgl" />
<Field abbreviation="cCin" level="4kmAgl" />
<Field abbreviation="cCin" level="4.5kmAgl" />
<Field abbreviation="cCin" level="5kmAgl" />
</Method>
<Method name="Average" levels="0-6kmAgl">
<Field abbreviation="cCin" level="0kmAgl" />
<Field abbreviation="cCin" level="0.5kmAgl" />
<Field abbreviation="cCin" level="1kmAgl" />
<Field abbreviation="cCin" level="1.5kmAgl" />
<Field abbreviation="cCin" level="2kmAgl" />
<Field abbreviation="cCin" level="2.5kmAgl" />
<Field abbreviation="cCin" level="3kmAgl" />
<Field abbreviation="cCin" level="3.5kmAgl" />
<Field abbreviation="cCin" level="4kmAgl" />
<Field abbreviation="cCin" level="4.5kmAgl" />
<Field abbreviation="cCin" level="5kmAgl" />
<Field abbreviation="cCin" level="5.5kmAgl" />
<Field abbreviation="cCin" level="6kmAgl" />
</Method>
</DerivedParameter>

21
edexOsgi/com.raytheon.uf.common.wxmath/src/com/raytheon/uf/common/wxmath/CapeFunc.java
View file

@ -56,6 +56,7 @@ import static java.lang.Math.sqrt;
* ------------ ---------- ----------- --------------------------
* Jun 03, 2013 2043 bsteffen Ported from meteolib C
* Aug 13, 2013 2262 njensen Moved from deriv params
* Aug 22, 2013 2289 bsteffen Performance improvements.
*
* </pre>
*
@ -128,11 +129,17 @@ public class CapeFunc {
if (Double.isNaN(pc) || Double.isNaN(pp) || Double.isNaN(tve)) {
tvp[k] = Double.NaN;
} else {
double t0 = tvc * pow(pp / pc, kapa);
/*
* The following line was originally pow(pp /pc, kapa) but
* this is much faster and accurate enough(15 decimal
* places).
*/
double pr = exp(kapa * log(pp / pc));
double t0 = tvc * pr;
if (pp > pc) {
tvp[k] = t0;
} else {
double td = tec * pow(pp / pc, kapa);
double td = tec * pr;
tvp[k] = td = temp_of_te(td, pp);
if (usetv > 0) {
tvp[k] *= pp
@ -396,11 +403,17 @@ public class CapeFunc {
tvp[k] = Double.NaN;
} else {
pp1 = pp;
double t0 = tvc * pow(pp / pc, kapa);
/*
* The following line was originally pow(pp /pc, kapa) but
* this is much faster and accurate enough(15 decimal
* places).
*/
double pr = exp(kapa * log(pp / pc));
double t0 = tvc * pr;
if (pp > pc) {
tvp[k] = t0;
} else {
double td = tec * pow(pp / pc, kapa);
double td = tec * pr;
tvp[k] = td = temp_of_te(td, pp);
if (usetv > 0) {
tvp[k] *= pp

167
edexOsgi/com.raytheon.uf.common.wxmath/src/com/raytheon/uf/common/wxmath/TempOfTe.java
View file

@ -22,6 +22,8 @@ package com.raytheon.uf.common.wxmath;
import static com.raytheon.uf.common.wxmath.AdiabeticTemperature.adiabatic_te;
import static java.lang.Math.sqrt;
import java.util.Arrays;
/**
* This routine calculates the saturation tempurature of an equivalent
* temperature at given pressure using the adiabatic definition
@ -34,6 +36,10 @@ import static java.lang.Math.sqrt;
* ------------ ---------- ----------- --------------------------
* Jun 03, 2013 2043 bsteffen Ported from meteolib C
* Aug 13, 2013 2262 njensen Moved from deriv params
* Aug 21, 2013 2289 bsteffen Add more pressure levels to TeTable.
* Remove redundant adiabatic_te calls.
* Use binary search in Arrays class.
* Return table values when possible.
*
* </pre>
*
@ -47,116 +53,111 @@ public class TempOfTe {
private static final int tmax = 333;
private static final int nval = 1 + tmax - tmin;
private static final int nt = 1 + tmax - tmin;
private static final double[] Te1000 = new double[nval];
private static final int pmin = 100;
private static final double[] Te850 = new double[nval];
private static final int pmax = 1000;
private static final double[] Te700 = new double[nval];
private static final int pspace = 25;
private static final double[] Te600 = new double[nval];
private static final int np = 1 + (pmax - pmin) / pspace;
private static final double[] Te500 = new double[nval];
private static final float[][] TeTable = new float[np][nt];
private static final double[] Te350 = new double[nval];
private static final double[] Te200 = new double[nval];
static {
int i = 0;
for (int t = tmin; t <= tmax; t += 1) {
Te1000[i] = adiabatic_te(t, 1000);
Te850[i] = adiabatic_te(t, 850);
Te700[i] = adiabatic_te(t, 700);
Te600[i] = adiabatic_te(t, 600);
Te500[i] = adiabatic_te(t, 500);
Te350[i] = adiabatic_te(t, 350);
Te200[i] = adiabatic_te(t, 200);
for (int p = pmin; p <= pmax; p += pspace) {
int j = 0;
for (int t = tmin; t <= tmax; t += 1) {
TeTable[i][j] = (float) adiabatic_te(t, p);
j += 1;
}
i += 1;
}
}
public static double temp_of_te(double te, double press) {
double[] TeLookup = null;
double base;
/* find correct table, check for beyond bounds of table */
if (press <= 250) {
TeLookup = Te200;
base = 200;
} else if (press <= 400) {
TeLookup = Te350;
base = 350;
} else if (press <= 550) {
TeLookup = Te500;
base = 500;
} else if (press <= 650) {
TeLookup = Te600;
base = 600;
} else if (press <= 750) {
TeLookup = Te700;
base = 700;
} else if (press <= 900) {
TeLookup = Te850;
base = 850;
} else {
TeLookup = Te1000;
base = 1000;
/* find the closest table row */
int pIndex = (int) ((press - pmin) / pspace);
if (pIndex < 0) {
pIndex = 0;
}
if (pIndex >= np) {
pIndex = np - 1;
}
float[] TeLookup = TeTable[pIndex];
if (te < TeLookup[1]) {
return te;
}
if (te >= TeLookup[nval - 1]) {
if (te >= TeLookup[nt - 1]) {
return Double.NaN;
}
double base = pmin + pIndex * pspace;
/* use table to get first guesses for value of temp */
double t1 = tmin;
double t2 = (int) te;
if (t2 > tmax) {
t2 = tmax;
}
double t;
while (t2 - t1 >= 3) {
t = (int) ((t1 + t2) / 2);
if (TeLookup[(int) t - tmin] > te) {
t2 = t;
} else if (TeLookup[(int) t - tmin] < te) {
t1 = t;
} else {
if (t1 < t - 1) {
t1 = t - 1;
}
if (t2 > t + 1) {
t2 = t + 1;
}
break;
}
int indexHigh = (int) te - tmin + 1;
if (indexHigh > nt) {
indexHigh = nt - 1;
}
double w = sqrt(base / press);
t1 = (1 - w) * TeLookup[(int) t1 - tmin] + w * t1;
t2 = (1 - w) * TeLookup[(int) t2 - tmin] + w * t2;
int indexLow = Arrays.binarySearch(TeLookup, 0, indexHigh, (float) te);
if (indexLow >= 0) {
if (base == press) {
return tmin + indexLow;
} else {
indexHigh = indexLow + 1;
indexLow = indexLow - 1;
}
} else {
indexHigh = -(indexLow) - 1;
indexLow = -(indexLow) - 2;
}
double tempLow = indexLow + tmin;
double tempHigh = indexHigh + tmin;
double diffLow = te - TeLookup[indexLow];
double diffHigh = TeLookup[indexHigh] - te;
if (base != press) {
/* use weight to take into account pressure difference */
double weight = sqrt(base / press);
tempLow = (1 - weight) * TeLookup[indexLow] + weight * tempLow;
diffLow = te - adiabatic_te(tempLow, press);
tempHigh = (1 - weight) * TeLookup[indexHigh] + weight * tempHigh;
diffHigh = adiabatic_te(tempHigh, press) - te;
}
if (diffLow < 0.01 && diffLow > -0.01) {
return tempLow;
}
if (diffHigh < 0.01 && diffHigh > -0.01) {
return tempHigh;
}
/* Iterate to find the exact solution */
double d1 = te - adiabatic_te(t1, press);
double d2 = adiabatic_te(t2, press) - te;
w = d2 / (d1 + d2);
t = w * t1 + (1 - w) * t2;
double d = adiabatic_te(t, press) - te;
int i = 0;
while (i++ < 10) {
if (d > 0.01) {
d2 = d;
t2 = t;
} else if (d < -0.01) {
d1 = -d;
t1 = t;
double weight = diffHigh / (diffLow + diffHigh);
double temp = weight * tempLow + (1 - weight) * tempHigh;
double diff = adiabatic_te(temp, press) - te;
for (int i = 0; i <= 10; i += 1) {
if (diff > 0.01) {
diffHigh = diff;
tempHigh = temp;
} else if (diff < -0.01) {
diffLow = -diff;
tempLow = temp;
} else {
break;
}
w = d2 / (d1 + d2);
t = w * t1 + (1 - w) * t2;
d = adiabatic_te(t, press) - te;
weight = diffHigh / (diffLow + diffHigh);
temp = weight * tempLow + (1 - weight) * tempHigh;
diff = adiabatic_te(temp, press) - te;
}
return t;
return temp;
}
}