diff --git a/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Cape.py b/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Cape.py
index cc3c7bc057..6222903a10 100644
--- a/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Cape.py
+++ b/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Cape.py
@@ -18,11 +18,7 @@
# further licensing information.
##
-# determine the location of the python ported meteolib
-# TODO put the meteolib location in the java MasterDerivScript instead of having
-# to manually determine its location like this
-import sys
-from numpy import zeros
+from com.raytheon.uf.viz.derivparam.python.function import CapeFunc
def execute(*args):
return __execute(*args)[0]
@@ -93,104 +89,9 @@ def __execute(*args):
pressureValue = pressure
pressure = zeros(temperatureValues.shape, temperatureValues.dtype)
pressure[:] = pressureValue
-
- return capeFunc(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, upperTerminationPressure)
-
-import numpy as np
-import ctypes as ct
-from com.raytheon.edex.meteoLib import MeteoLibUtil
-
-def capeFunc(usetv, p_dat, tve_dat, p0, th0, sh0, ptop = None):
- """ Use the native capeFunc function
- """
-
- # define the c_float_ptr type
- c_float_ptr = ct.POINTER(ct.c_float)
-
- # determine the input dimensions
- dataShape = tve_dat.shape
-
- nx = dataShape[2] if len(dataShape) > 2 else None
- ny = dataShape[1] if len(dataShape) > 1 else None
- nz = dataShape[0]
-
- gridArea = nx * ny
-
- # flatten all input arrays
- p_dat = np.copy(p_dat)
- tve_dat = np.copy(tve_dat)
- p0 = np.copy(p0)
- th0 = np.copy(th0)
- sh0 = np.copy(sh0)
-
- p_dat[np.isnan(p_dat)] = 1e37
- tve_dat[np.isnan(tve_dat)] = 1e37
- p0[np.isnan(p0)] = 1e37
- th0[np.isnan(th0)] = 1e37
- sh0[np.isnan(sh0)] = 1e37
-
- p_dat.resize((nz, nx * ny,))
- tve_dat.resize((nz, nx * ny,))
- p0.resize((p0.size,))
- th0.resize((th0.size,))
- sh0.resize((sh0.size,))
-
- if ptop != None:
- ptop.resize((ptop.size,))
-
-
- # load the library
- meteoLibPath = MeteoLibUtil.getSoPath()
- meteoLib = np.ctypeslib.load_library(meteoLibPath,"")
- capeFunc = meteoLib.capeFunc if ptop == None else meteoLib.capeFuncTop
-
- # "define" the capeFunc signature
- capeFunc.restype = ct.c_int # return type
- capeFunc.argtypes = [ct.c_float,
- ct.POINTER(c_float_ptr),
- ct.POINTER(c_float_ptr),
- c_float_ptr,
- c_float_ptr,
- c_float_ptr,
- ct.c_int,
- ct.c_int,
- ct.c_int,
- ct.c_int,
- c_float_ptr,
- c_float_ptr]
-
- if ptop != None:
- capeFunc.argtypes.append(c_float_ptr)
-
- # result arrays
- cape_dat = np.zeros(gridArea,p_dat.dtype)
- cin_dat = np.zeros(gridArea,p_dat.dtype)
-
- capeFuncArgs = [ct.c_float(usetv),
-
- # get c_style pointers to the 2D input arrays
- (c_float_ptr*len(p_dat))(*[row.ctypes.data_as(c_float_ptr) for row in p_dat]),
- (c_float_ptr*len(tve_dat))(*[row.ctypes.data_as(c_float_ptr) for row in tve_dat]),
-
- p0.ctypes.data_as(c_float_ptr),
- th0.ctypes.data_as(c_float_ptr),
- sh0.ctypes.data_as(c_float_ptr),
- ct.c_int(nx),
- ct.c_int(nx),
- ct.c_int(ny),
- ct.c_int(nz),
- cape_dat.ctypes.data_as(c_float_ptr),
- cin_dat.ctypes.data_as(c_float_ptr)]
-
- if ptop != None:
- capeFuncArgs.append(ptop.ctypes.data_as(c_float_ptr))
-
- val = capeFunc(*capeFuncArgs)
- if val == 1 :
- raise MemoryError('Cape derived parameter ran out of memory')
- exit
- # resize the cape and cin data to the appropriate grid size
- cape_dat.resize((ny,nx))
- cin_dat.resize((ny,nx))
-
- return cape_dat, cin_dat
+ if upperTerminationPressure is None:
+ threeDshape = pressureValues.shape
+ return CapeFunc.capeFunc(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__
+ else:
+ threeDshape = pressureValues.shape
+ return CapeFunc.capeFuncTop(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, upperTerminationPressure, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__
diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/AdiabeticTemperature.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/AdiabeticTemperature.java
new file mode 100644
index 0000000000..27e88a0457
--- /dev/null
+++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/AdiabeticTemperature.java
@@ -0,0 +1,51 @@
+/**
+ * This software was developed and / or modified by Raytheon Company,
+ * pursuant to Contract DG133W-05-CQ-1067 with the US Government.
+ *
+ * U.S. EXPORT CONTROLLED TECHNICAL DATA
+ * This software product contains export-restricted data whose
+ * export/transfer/disclosure is restricted by U.S. law. Dissemination
+ * to non-U.S. persons whether in the United States or abroad requires
+ * an export license or other authorization.
+ *
+ * Contractor Name: Raytheon Company
+ * Contractor Address: 6825 Pine Street, Suite 340
+ * Mail Stop B8
+ * Omaha, NE 68106
+ * 402.291.0100
+ *
+ * See the AWIPS II Master Rights File ("Master Rights File.pdf") for
+ * further licensing information.
+ **/
+package com.raytheon.uf.viz.derivparam.python.function;
+
+import static java.lang.Math.exp;
+
+/**
+ * This routine calculates the equivalent tempurature of a temperature and
+ * pressure using the adiabatic definition, assuming saturation put a fudge
+ * factor into L/cp to get agreement of moist adiabats with a published
+ * thermodynamic diagram *
+ *
+ *
+ *
+ * SOFTWARE HISTORY
+ *
+ * Date Ticket# Engineer Description
+ * ------------ ---------- ----------- --------------------------
+ * Jun 3, 2013 2043 bsteffen Ported from meteolib C
+ *
+ *
+ *
+ * @author bsteffen
+ * @version 1.0
+ */
+
+public class AdiabeticTemperature {
+
+ public static double adiabatic_te(double temp, double press) {
+ double e = exp(26.660820 - 0.0091379024 * temp - 6106.396 / temp);
+ e = 0.622 * e / (press - e);
+ return temp * exp(2740.0 * e / temp);
+ }
+}
diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CapeFunc.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CapeFunc.java
new file mode 100644
index 0000000000..dfe0b72abc
--- /dev/null
+++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CapeFunc.java
@@ -0,0 +1,627 @@
+/**
+ * This software was developed and / or modified by Raytheon Company,
+ * pursuant to Contract DG133W-05-CQ-1067 with the US Government.
+ *
+ * U.S. EXPORT CONTROLLED TECHNICAL DATA
+ * This software product contains export-restricted data whose
+ * export/transfer/disclosure is restricted by U.S. law. Dissemination
+ * to non-U.S. persons whether in the United States or abroad requires
+ * an export license or other authorization.
+ *
+ * Contractor Name: Raytheon Company
+ * Contractor Address: 6825 Pine Street, Suite 340
+ * Mail Stop B8
+ * Omaha, NE 68106
+ * 402.291.0100
+ *
+ * See the AWIPS II Master Rights File ("Master Rights File.pdf") for
+ * further licensing information.
+ **/
+package com.raytheon.uf.viz.derivparam.python.function;
+
+import static com.raytheon.uf.viz.derivparam.python.function.AdiabeticTemperature.adiabatic_te;
+import static com.raytheon.uf.viz.derivparam.python.function.TempOfTe.temp_of_te;
+import static java.lang.Math.exp;
+import static java.lang.Math.log;
+import static java.lang.Math.pow;
+import static java.lang.Math.sqrt;
+import jep.INumpyable;
+
+/**
+ * We input theta and specific humidity for initial parcel because these are
+ * things that can be arithemitically averaged for a mixed layer. If usetv=1,
+ * buoyancy is done with virtual temp, if usetv=0 with temp. If usetv=0, must
+ * supply temps in tve_dat.
+ *
+ * There are a few small changes from the original C version. First all math is
+ * done using doubles instead of floats since that is what the java.lang.Math
+ * library uses and it gets more accurate results. Second all uses of 1e37 have
+ * been replaces with Double.NaN to be more consistant with the rest of derived
+ * parameters. Finally the loops have been reordered so that each grid cell is
+ * iterated in the outermost loop which allows all intermediary calculations to
+ * be stored as doubles instead of float[nx*ny] which dramatically reduces the
+ * memory footprint.
+ *
+ *
+ *
+ * SOFTWARE HISTORY
+ *
+ * Date Ticket# Engineer Description
+ * ------------ ---------- ----------- --------------------------
+ * Jun 3, 2013 2043 bsteffen Ported from meteolib C
+ *
+ *
+ *
+ * @author bsteffen
+ * @version 1.0
+ */
+
+public class CapeFunc {
+
+ private static final float c0 = 26.66082f;
+
+ private static final float c1 = 0.0091379024f;
+
+ private static final float c2 = 6106.396f;
+
+ private static final float c_1 = 223.1986f;
+
+ private static final float c_2 = 0.0182758048f;
+
+ private static final float kapa = 0.286f;
+
+ private static final float kapa_1 = 3.498257f;
+
+ public static CapeCinPair capeFunc(float usetv, float[] p_dat,
+ float[] tve_dat, float[] p0, float[] th0, float[] sh0, int nx,
+ int ny, int nz) {
+ int n2 = nx * ny;
+
+ float[] cin = new float[n2];
+ float[] cap = new float[n2];
+
+ // Calculate the parcel equivalent temp, virtual temp, and press at LCL.
+ // Make working copy of sfc press, use as press below current 3d
+ // pressure.
+ for (int i = 0; i < n2; i += 1) {
+ double tec = 0;
+ double tvc = 0;
+ double pc = 0;
+ double pp0 = p0[i];
+ double pp1 = pp0;
+ if (Double.isNaN(pp0) || Double.isNaN(th0[i])
+ || Double.isNaN(sh0[i])
+ || sh0[i] < 0.0005) {
+ tec = tvc = pc = Double.NaN;
+ } else {
+ double t0 = th0[i] * pow(pp0 / 1000, kapa);
+ double b = c0 - log(pp0 / (622.0 / sh0[i] + 0.378));
+ double td = (b - sqrt(b * b - c_1)) / c_2;
+ double tdc = td - (t0 - td)
+ * (-0.37329638 + 41.178204 / t0 + 0.0015945203 * td);
+ pc = pp0 * pow(tdc / t0, kapa_1);
+ tec = adiabatic_te(tdc, pc);
+ tvc = td * (1 + usetv * 0.000608 * sh0[i]);
+ }
+
+ // Initialize md and pmd, which will be pressure of and max Te
+ // delta.
+ double md = 0;
+ double pmd = 0;
+
+ // Now calculate the virtual temperature of the parcel at the
+ // pressures in the input data. Then difference it from the
+ // environmental temp, which has been tweaked to not be cooler than
+ // dry adiabatic from the parcel start. Record the level of max
+ // parcel difference.
+ double[] tvp = new double[nz];
+
+ for (int k = 0; k < nz; k += 1) {
+ float pp = p_dat[k * n2 + i];
+ float tve = tve_dat[k * n2 + i];
+ if (Double.isNaN(pc) || Double.isNaN(pp) || Double.isNaN(tve)) {
+ tvp[k] = Double.NaN;
+ } else {
+ double t0 = tvc * pow(pp / pc, kapa);
+ if (pp > pc) {
+ tvp[k] = t0;
+ } else {
+ double td = tec * pow(pp / pc, kapa);
+ tvp[k] = td = temp_of_te(td, pp);
+ if (usetv > 0) {
+ tvp[k] *= pp
+ / (pp - exp(25.687958917 - c1 * td - c2
+ / td));
+ }
+ }
+ if (tve < t0) {
+ tvp[k] -= t0;
+ } else {
+ tvp[k] -= tve;
+ }
+ if (pp > pc || tvp[k] < md) {
+ continue;
+ }
+ md = tvp[k];
+ pmd = pp;
+ }
+ }
+
+ // This loop performs the actual cape and cin calculation. Here we
+ // will reuse storage for virt temp, equiv temp, and max delta for
+ // prev parcel temp, neg and pos. neg and pos are pending negative
+ // and positive contributions we have not yet added into the cape
+ // and cin yet.
+ double neg = 0;
+ double pos = 0;
+ cin[i] = cap[i] = Float.NaN;
+
+ double tvp1 = tvp[0];
+ pp1 = p_dat[0];
+ for (int k = 1; k < nz; k += 1) {
+ float pp = p_dat[k * n2 + i];
+ if (Double.isNaN(pp0)) {
+ continue;
+ } else if (Double.isNaN(pp1) || Double.isNaN(tvp1)) {
+ ;
+ } else if (pp >= pp1 || Double.isNaN(tvp[k])) {
+ continue;
+ } else if (pp >= pp0) {
+ ;
+ } else {
+ // Now we finally have the data we need for calculating
+ // the cape/cin contribution for this layer.
+ if (Double.isNaN(cap[i])) {
+ cap[i] = cin[i] = 0;
+ }
+ if (pmd == 0) {
+ continue; // No parcel delta>0, we're done.
+ }
+
+ // First deal with possibility of bottom lvl being below the
+ // initial parcel.
+ double dlnp;
+ double dn;
+ if (pp1 > pp0) {
+ dlnp = log(pp0 / pp);
+ dn = 0;
+ } else {
+ dlnp = log(pp1 / pp);
+ dn = dlnp * 287 * tvp1;
+ }
+
+ // Now deal with the fact that not allowing superadiabatic
+ // layers means no cape below condensation pressure.
+ double up;
+ if (pp1 >= pc) {
+ if (dn > 0) {
+ dn = 0;
+ }
+ if (tvp[k] <= 0) {
+ up = dlnp * 287 * tvp[k];
+ } else if (pp >= pc) {
+ up = 0;
+ } else {
+ up = log(pc / pp) * 287 * tvp[k];
+ }
+ } else {
+ up = dlnp * 287 * tvp[k];
+ }
+
+ // Deal with where the break point is.
+ double b = up * dn >= 0 ? 0.5 : up / (up - dn);
+ up *= b;
+ dn *= (1 - b);
+
+ // Now consider this layer's contribution, taking into
+ // account transitions between positive and negative
+ // acceleration.
+ if (up == 0 && dn == 0) {
+ ;
+ // Continuing deceleration.
+ } else if (up <= 0
+ && (dn < 0 || dn == 0 && (pp < pmd || pos == 0))) {
+ neg -= up + dn;
+
+ // Continuing upward acceleration.
+ } else if (up >= 0
+ && (dn > 0 || dn == 0 && (pp < pmd || neg == 0))) {
+ pos += up + dn;
+ if (pp > pmd && cap[i] + pos <= cin[i] + neg) {
+ ; // no net cape and below max delta
+ } else if (pp > pmd || cap[i] == 0) {
+ // below max delta or cape uninitialized
+ cap[i] += pos;
+ cin[i] += neg;
+ neg = pos = 0;
+ } else if (pos >= neg) {
+ // cape initialized and net positive contribution
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ } else if (up > 0 && dn <= 0) {
+ // Transition to upward acceleration.
+ neg += -dn;
+ if (pp1 <= pmd) {
+ // above max delta, only use net pos contribution
+ pos += up;
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ } else if (pp <= pmd) {
+ // straddle max delta, force cape initialization
+ if (cap[i] == 0) {
+ cin[i] += neg;
+ cap[i] += pos;
+ } else if (neg > pos) {
+ cin[i] += neg - pos;
+ } else {
+ cap[i] += pos - neg;
+ }
+ cap[i] += up;
+ neg = pos = 0;
+ } else if (cap[i] + pos + up <= cin[i] + neg) {
+ // no net cape to this point
+ if (cap[i] + pos > 0) {
+ // reinitialize if there was cape before
+ cin[i] -= cap[i] + pos;
+ pos = cap[i] = 0;
+ }
+ cin[i] += neg;
+ pos += up;
+ neg = 0;
+ } else if (cap[i] == 0) { // initialize cape
+ cap[i] += pos + up;
+ cin[i] += neg;
+ neg = pos = 0;
+ } else { // what remains, only use net pos contribution
+ pos += up;
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ }
+ } else {
+ // Transition to decceleration.
+ pos += dn;
+ if (pp1 <= pmd) {
+ // above max delta, only use net pos contribution
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ neg += -up;
+ } else if (cap[i] + pos <= cin[i] + neg - up) {
+ // no net cape to this point
+ if (cap[i] > 0) {
+ // reinitialize if there was cape before
+ cin[i] -= cap[i] + pos;
+ pos = cap[i] = 0;
+ }
+ cin[i] += neg - up;
+ pos = neg = 0;
+ } else if (cap[i] == 0) { // initialize cape
+ cap[i] += pos;
+ cin[i] += neg - up;
+ neg = pos = 0;
+ } else { // what remains, only use net pos contribution
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ neg += -up;
+ }
+ }
+ }
+ // Make current layer top next layer bottom.
+ tvp1 = tvp[k];
+ pp1 = pp;
+ }
+ }
+ return new CapeCinPair(ny, nx, cap, cin);
+ }
+
+ // In this version we stop the computation at some arbitrary upper level,
+ // ptop.
+ public static CapeCinPair capeFuncTop(float usetv, float[] p_dat,
+ float[] tve_dat, float[] p0, float[] th0, float[] sh0,
+ float[] ptop, int nx, int ny, int nz) {
+ int n2 = nx * ny;
+
+ float[] cin = new float[n2];
+ float[] cap = new float[n2];
+
+ // Calculate the parcel equivalent temp, virtual temp, and press at LCL.
+ // Make working copy of sfc press, use as press below current 3d
+ // pressure.
+ for (int i = 0; i < n2; i += 1) {
+ double tec = 0;
+ double tvc = 0;
+ double pc = 0;
+ double pp0 = p0[i];
+ double pp1 = pp0;
+ double pfin = ptop[i];
+ if (Double.isNaN(pp0) || Double.isNaN(th0[i])
+ || Double.isNaN(sh0[i]) || sh0[i] < 0.0005
+ || pp0 < pfin) {
+ tec = tvc = pc = Double.NaN;
+ } else {
+ double t0 = th0[i] * pow(pp0 / 1000, kapa);
+ double b = c0 - log(pp0 / (622.0 / sh0[i] + 0.378));
+ double td = (b - sqrt(b * b - c_1)) / c_2;
+ double tdc = td - (t0 - td)
+ * (-0.37329638 + 41.178204 / t0 + 0.0015945203 * td);
+ pc = pp0 * pow(tdc / t0, kapa_1);
+ tec = adiabatic_te(tdc, pc);
+ tvc = td * (1 + usetv * 0.000608 * sh0[i]);
+ }
+
+ // Initialize md and pmd, which will be pressure of and max Te
+ // delta.
+ double md = 0;
+ double pmd = 0;
+
+ // Now calculate the virtual temperature of the parcel at the
+ // pressures in the input data. Then difference it from the
+ // environmental temp, which has been tweaked to not be cooler than
+ // dry adiabatic from the parcel start. Record the level of max
+ // parcel difference.
+ double[] tvp = new double[nz];
+
+ for (int k = 0; k < nz; k += 1) {
+ float pp = p_dat[k * n2 + i];
+ float tve = tve_dat[k * n2 + i];
+ if (Double.isNaN(pc) || Double.isNaN(pp) || Double.isNaN(tve)
+ || pp1 <= pfin) {
+ tvp[k] = Double.NaN;
+ } else {
+ pp1 = pp;
+ double t0 = tvc * pow(pp / pc, kapa);
+ if (pp > pc) {
+ tvp[k] = t0;
+ } else {
+ double td = tec * pow(pp / pc, kapa);
+ tvp[k] = td = temp_of_te(td, pp);
+ if (usetv > 0) {
+ tvp[k] *= pp
+ / (pp - exp(25.687958917 - c1 * td - c2
+ / td));
+ }
+ }
+ if (tve < t0) {
+ tvp[k] -= t0;
+ } else {
+ tvp[k] -= tve;
+ }
+ if (pp > pc || pp < pfin || tvp[k] < md) {
+ continue;
+ }
+ md = tvp[k];
+ pmd = pp;
+ }
+ }
+
+ // This loop performs the actual cape and cin calculation. Here we
+ // will reuse storage for virt temp, equiv temp, and max delta for
+ // prev parcel temp, neg and pos. neg and pos are pending negative
+ // and positive contributions we have not yet added into the cape
+ // and cin yet.
+ double neg = 0;
+ double pos = 0;
+ cin[i] = cap[i] = Float.NaN;
+
+ double tvp1 = tvp[0];
+ pp1 = p_dat[i];
+ for (int k = 1; k < nz; k += 1) {
+ float pp = p_dat[k * n2 + i];
+ if (Double.isNaN(pp0)) {
+ continue;
+ } else if (Double.isNaN(pp1) || Double.isNaN(tvp1)) {
+ ;
+ } else if (pp >= pp1 || Double.isNaN(tvp[k])) {
+ continue;
+ } else if (pp >= pp0) {
+ ;
+ } else {
+ // Now we finally have the data we need for calculating
+ // the cape/cin contribution for this layer.
+ if (Double.isNaN(cap[i])) {
+ cap[i] = cin[i] = 0;
+ }
+ if (pmd == 0) {
+ continue; // No parcel delta>0, we're done.
+ }
+
+ // First deal with possibility of bottom lvl being below the
+ // initial parcel and/or hitting the top of the computation.
+ double dlnp = 0;
+ double dn;
+ if (pp < pfin) {
+ double b = log(pp1 / pp);
+ dlnp = log(pp1 / pfin);
+ tvp[k] = tvp1 + (dlnp / b) * (tvp[k] - tvp1);
+ }
+ if (pp1 > pp0) {
+ if (pp < pfin) {
+ dlnp = log(pp0 / pfin);
+ } else {
+ dlnp = log(pp0 / pp);
+ }
+ dn = 0;
+ } else {
+ if (pp >= pfin) {
+ dlnp = log(pp1 / pp);
+ }
+ dn = dlnp * 287 * tvp1;
+ }
+
+ // Now deal with the fact that not allowing superadiabatic
+ // layers means no cape below condensation pressure.
+ double up;
+ if (pp1 >= pc) {
+ if (dn > 0) {
+ dn = 0;
+ }
+ if (tvp[k] <= 0) {
+ up = dlnp * 287 * tvp[k];
+ } else if (pp >= pc) {
+ up = 0;
+ } else if (pp < pfin) {
+ up = log(pc / pfin) * 287 * tvp[k];
+ } else {
+ up = log(pc / pp) * 287 * tvp[k];
+ }
+ } else {
+ up = dlnp * 287 * tvp[k];
+ }
+
+ // Deal with where the break point is.
+ double b = up * dn >= 0 ? 0.5 : up / (up - dn);
+ up *= b;
+ dn *= (1 - b);
+
+ // Now consider this layer's contribution, taking into
+ // account transitions between positive and negative
+ // acceleration.
+ if (up == 0 && dn == 0) {
+ ;
+ // Continuing deceleration.
+ } else if (up <= 0
+ && (dn < 0 || dn == 0 && (pp < pmd || pos == 0))) {
+ neg -= up + dn;
+
+ // Continuing upward acceleration.
+ } else if (up >= 0
+ && (dn > 0 || dn == 0 && (pp < pmd || neg == 0))) {
+ pos += up + dn;
+ if (pp > pmd && cap[i] + pos <= cin[i] + neg) {
+ ; // no net cape and below max delta
+ } else if (pp > pmd || cap[i] == 0) {
+ // below max delta or cape uninitialized
+ cap[i] += pos;
+ cin[i] += neg;
+ neg = pos = 0;
+ } else if (pos >= neg) {
+ // cape initialized and net positive contribution
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ } else if (up > 0 && dn <= 0) {
+ // Transition to upward acceleration.
+ neg += -dn;
+ if (pp1 <= pmd) {
+ // above max delta, only use net pos contribution
+ pos += up;
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ } else if (pp <= pmd) {
+ // straddle max delta, force cape initialization
+ if (cap[i] == 0) {
+ cin[i] += neg;
+ cap[i] += pos;
+ } else if (neg > pos) {
+ cin[i] += neg - pos;
+ } else {
+ cap[i] += pos - neg;
+ }
+ cap[i] += up;
+ neg = pos = 0;
+ } else if (cap[i] + pos + up <= cin[i] + neg) {
+ // no net cape to this point
+ if (cap[i] + pos > 0) {
+ // reinitialize if there was cape before
+ cin[i] -= cap[i] + pos;
+ pos = cap[i] = 0;
+ }
+ cin[i] += neg;
+ pos += up;
+ neg = 0;
+ } else if (cap[i] == 0) { // initialize cape
+ cap[i] += pos + up;
+ cin[i] += neg;
+ neg = pos = 0;
+ } else { // what remains, only use net pos contribution
+ pos += up;
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ }
+ } else {
+ // Transition to decceleration.
+ pos += dn;
+ if (pp1 <= pmd) {
+ // above max delta, only use net pos contribution
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ neg += -up;
+ } else if (cap[i] + pos <= cin[i] + neg - up) {
+ // no net cape to this point
+ if (cap[i] > 0) {
+ // reinitialize if there was cape before
+ cin[i] -= cap[i] + pos;
+ pos = cap[i] = 0;
+ }
+ cin[i] += neg - up;
+ pos = neg = 0;
+ } else if (cap[i] == 0) { // initialize cape
+ cap[i] += pos;
+ cin[i] += neg - up;
+ neg = pos = 0;
+ } else { // what remains, only use net pos contribution
+ if (pos >= neg) {
+ cap[i] += (pos - neg);
+ neg = pos = 0;
+ }
+ neg += -up;
+ }
+ }
+ }
+ // Make current layer top next layer bottom.
+ tvp1 = tvp[k];
+ pp1 = pp;
+ }
+ }
+ return new CapeCinPair(ny, nx, cap, cin);
+ }
+
+ public static class CapeCinPair implements INumpyable {
+
+ private final int nx;
+
+ private final int ny;
+
+ private final float[] cape;
+
+ private final float[] cin;
+
+ public CapeCinPair(int nx, int ny, float[] cape, float[] cin) {
+ this.nx = nx;
+ this.ny = ny;
+ this.cape = cape;
+ this.cin = cin;
+ }
+
+ @Override
+ public Object[] getNumPy() {
+ return new Object[] { cape, cin };
+ }
+
+ @Override
+ public int getNumpyX() {
+ return nx;
+ }
+
+ @Override
+ public int getNumpyY() {
+ return ny;
+ }
+
+ }
+
+}
diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/TempOfTe.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/TempOfTe.java
new file mode 100644
index 0000000000..1b59a72929
--- /dev/null
+++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/TempOfTe.java
@@ -0,0 +1,161 @@
+/**
+ * This software was developed and / or modified by Raytheon Company,
+ * pursuant to Contract DG133W-05-CQ-1067 with the US Government.
+ *
+ * U.S. EXPORT CONTROLLED TECHNICAL DATA
+ * This software product contains export-restricted data whose
+ * export/transfer/disclosure is restricted by U.S. law. Dissemination
+ * to non-U.S. persons whether in the United States or abroad requires
+ * an export license or other authorization.
+ *
+ * Contractor Name: Raytheon Company
+ * Contractor Address: 6825 Pine Street, Suite 340
+ * Mail Stop B8
+ * Omaha, NE 68106
+ * 402.291.0100
+ *
+ * See the AWIPS II Master Rights File ("Master Rights File.pdf") for
+ * further licensing information.
+ **/
+package com.raytheon.uf.viz.derivparam.python.function;
+
+import static com.raytheon.uf.viz.derivparam.python.function.AdiabeticTemperature.adiabatic_te;
+import static java.lang.Math.sqrt;
+
+/**
+ * This routine calculates the saturation tempurature of an equivalent
+ * temperature at given pressure using the adiabatic definition
+ *
+ *
+ *
+ * SOFTWARE HISTORY
+ *
+ * Date Ticket# Engineer Description
+ * ------------ ---------- ----------- --------------------------
+ * Jun 3, 2013 2043 bsteffen Ported from meteolib C
+ *
+ *
+ *
+ * @author bsteffen
+ * @version 1.0
+ */
+
+public class TempOfTe {
+
+ private static final int tmin = 193;
+
+ private static final int tmax = 333;
+
+ private static final int nval = 1 + tmax - tmin;
+
+ private static final double[] Te1000 = new double[nval];
+
+ private static final double[] Te850 = new double[nval];
+
+ private static final double[] Te700 = new double[nval];
+
+ private static final double[] Te600 = new double[nval];
+
+ private static final double[] Te500 = new double[nval];
+
+ 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);
+ 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;
+ }
+ if (te < TeLookup[1]) {
+ return te;
+ }
+ if (te >= TeLookup[nval - 1]) {
+ return Double.NaN;
+ }
+ /* 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;
+ }
+ }
+
+ double w = sqrt(base / press);
+ t1 = (1 - w) * TeLookup[(int) t1 - tmin] + w * t1;
+ t2 = (1 - w) * TeLookup[(int) t2 - tmin] + w * t2;
+
+ /* 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;
+ } else {
+ break;
+ }
+ w = d2 / (d1 + d2);
+ t = w * t1 + (1 - w) * t2;
+ d = adiabatic_te(t, press) - te;
+ }
+ return t;
+ }
+}