# Adapted for numpy/ma/cdms2 by convertcdms.py
import MV2
import cdms2
import genutil
import unidata
import vcs
import numpy
from vcs import VCS_validation_functions
thermo_objects = []


def Es(T, method=None):
    """Computes saturated pressure in Pa given T in K, using the method:
    1: Hyland-Wexler formulation, polynomial coeff (absolute norm)
    2: Wexler formulation
    3: Hyland-Wexler formulation, polynomial coeff (relative error norm)
    4: classic Goff Gratch equation
    5: 6.112*numpy.ma.exp(17.67*tempc/(tempc+243.5))

    Default is method 1

    Note: 1 and 2 use method 3 where T is not : 173.15 < T < 473.15
    ref for 1, 2 and 3:
    <NAME>., Journal of Applied Met., Vol 31, Dec 1992
    ( http://ams.allenpress.com/perlserv/?request=get-document&\
doi=10.1175%2F1520-0450(1992)031%3C1507%3APFTSVP%3E2.0.CO%3B2&ct=1 )
    """
    if method is None:
        method = 1

    if method == 1:
        # Put in C
        x = T - 273.15
        # Water vapor
        c0 = 0.611220713E03
        c1 = 0.443944344E02
        c2 = 0.143195336E01
        c3 = 0.263350515E-01
        c4 = 0.310636053E-03
        c5 = 0.185218710E-05
        c6 = 0.103440324E-07
        c7 = -0.468258100E-10
        c8 = 0.466533033E-13
        eswat = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (c5 + x * (c6 +
                                                                            x * (c7 + x * c8)))))))
        # ice
        c0 = .611153246E03
        c1 = .503261230E02
        c2 = .188595709E01
        c3 = .422115970E-01
        c4 = .620376691E-03
        c5 = .616082536E-05
        c6 = .405172828E-07
        c7 = .161492905E-09
        c8 = .297886454E-12
        esice = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (c5 + x * (c6 +
                                                                            x * (c7 + x * c8)))))))
        # Combine
        es = MV2.where(MV2.less(T, 273.15), esice, eswat)
        # Overwrite values outside valid range with method 2
        mn, mx = genutil.minmax(T)
        if mn < 173.16 or mx > 473.15:
            es2 = Es(T, method=2)
            es = MV2.where(MV2.less(T, 173.16), es2, es)
            es = MV2.where(MV2.greater(T, 473.15), es2, es)
    elif method == 2:
        # over water
        g0 = -0.29912729E4
        g1 = -0.60170128E4
        g2 = 0.1887643854E2
        g3 = -0.28354721E-1
        g4 = 0.17838301E-4
        g5 = -0.84150417E-9
        g6 = 0.44412543E-12
        g7 = 0.2858487E1
        # over ice
        k0 = -0.58653696e4
        k1 = 0.2224103300E2
        k2 = 0.13749042E-1
        k3 = -0.34031775E-4
        k4 = 0.26967687E-7
        k5 = 0.6918651
        esice = (k0 + (k1 + k5 * MV2.log(T) + (k2 + (
                       k3 + k4 * T) * T) * T) * T) / T  # over ice
        eswat = (g0 + (g1 + (g2 + g7 * MV2.log(T) + (g3 + (g4 + (
                             g5 + g6 * T) * T) * T) * T) * T) * T) / T ** 2  # over water
        es = MV2.where(MV2.less(T, 273.15), esice, eswat)
        es = MV2.exp(es)
    elif method == 3:
        # Put in C
        x = T - 273.15
        # Water vapor
        c0 = 0.611213476E03
        c1 = 0.444007856E02
        c2 = 0.143064234E01
        c3 = 0.264461437E-01
        c4 = 0.305930558E-03
        c5 = 0.196237241E-05
        c6 = 0.892344772E-08
        c7 = -0.373208410E-10
        c8 = 0.209339997E-13
        eswat = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (
                          c5 + x * (c6 + x * (c7 + x * c8)))))))
        # ice
        c0 = .611123516E03
        c1 = .503109514E02
        c2 = .1888369801E01
        c3 = .420547422E-01
        c4 = .614396778E-03
        c5 = .602780717E-05
        c6 = .387940929E-07
        c7 = .149436277E-09
        c8 = .262655803E-12
        esice = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (
                          c5 + x * (c6 + x * (c7 + x * c8)))))))
        # Combine
        es = MV2.where(MV2.less(T, 273.15), esice, eswat)
        # Overwrite values outside valid range with method 2
        mn, mx = genutil.minmax(T)
        if mn < 173.16 or mx > 473.15:
            es2 = Es(T, method=2)
            es = MV2.where(MV2.less(T, 173.16), es2, es)
            es = MV2.where(MV2.greater(T, 473.15), es2, es)
    elif method == 4:
        est = 101324.6  # Pa
        Ts = 373.16 / T
        a = -7.90298
        b = 5.02808
        c = -1.3816E-7
        d = 11.344
        f = 8.1328E-3
        h = -3.49149
        maxexp = int(numpy.log10(numpy.finfo(numpy.float).max))
        minexp = 1 - a
        es = a * (Ts - 1.)
        es = es + b * numpy.ma.log10(Ts)
        A = d * (1. - Ts)
        A = numpy.ma.masked_greater(A, maxexp)
        A = numpy.ma.masked_less(A, minexp)
        es = es + c * (numpy.ma.power(10, A) - 1.)
        A = h * (Ts - 1.)
        A = numpy.ma.masked_greater(A, maxexp)
        A = numpy.ma.masked_less(A, minexp)
        es = es + f * (numpy.ma.power(10, A) - 1.)
        es = est * numpy.ma.power(10, es)
    elif method == 5:
        tempc = T - 273.15
        es = 611.2 * numpy.ma.exp(17.67 * tempc / (tempc + 243.5))
    return es


def getPvalues(T):
    """computes pressure values from level axis,
    if not designated level then use last axis...
    """
    P = T.getLevel()
    if P is None:
        P = T.getAxis(-1)
    P = P.clone()
    try:  # tries to convert to Pa
        for i in range(len(P)):
            u = unidata.udunits(P[i], P.units)
            P[i] = u.to("Pa").value
    except Exception:
        pass
    return P[:]


def Ws(T, P=None, method=None):
    """ Compute saturated mixing ratio at T,P, cpomputing es using the method:
    T in K
    P in Pa
    """
    es = Es(T, method)
    if P is None:
        P = getPvalues(T)

    return .622 * es / (P - es)


def gammaw(T, rh, P=None, method=None):
    """
    Function to calculate the moist adiabatic lapse rate (K/Pa) based
    on the temperature, pressure, and rh of the environment.
    T in K
    rh in %
    P in Pa
    if P is not passed will be optaines from T level axis (or last axis)
    """
    if P is None:
        P = getPvalues(T)
    ws = Ws(T, P, method=method)
    w = rh * ws / 100
    Tv = T * (1.0 + 0.6 * w)
    latent = 1000 * (2502.2 - 2.43089 * (T - 273.15))
    A = 1.0 + latent * ws / (287 * T)
    B = 1.0 + 0.622 * latent * latent * ws / (1005 * 287 * T * T)
    Density = P / (287 * Tv)  # in SI
    lapse = (A / B) / (1005. * Density)
    return lapse


def LiftWet(To, Po, Pend, deltaP=1000., method=None):
    """
    Lift a parcel moist adiabatically from Po to Pend.
    Initial temperature is To in K.
    To in K
    Po in Pa
    detail # number to increment P slowly
    returns T at all P (every deltaP Pa), and P
    """
    if method is None:
        method = 5  # MUCH faster
    temp = [To, ]
    pres = [Po, ]
    while (pres[-1] - deltaP >= Pend):
        tmp = temp[-1] - deltaP * gammaw(temp[-1], 100, pres[-1] - deltaP / 2, method=method)
        temp.append(tmp)
        pres.append(pres[-1] - deltaP)
    return MV2.array(temp, id='T'), MV2.array(pres, id='P')


def Dewpoint(T, rh, method=None):
    """Computes dewpoint profile from T and rh
    Usage:
      Td = Dewpoint(T,rh,method=None)
    Where:
     T is in K
     rh is in % (relative humidity)
     method used to compute Es:
         1: Hyland-Wexler formulation, polynomial coeff (absolute norm) (default)
         2: Wexler formulation
         3: Hyland-Wexler formulation, polynomial coeff (relative error norm)
         4: classic Goff Gratch equation
         5: 6.112*numpy.ma.exp(17.67*tempc/(tempc+243.5))



    Assumes constant latent temperature (except over ice or water)
    """
    To = 273.15     # K
    Rv = 461.53     # J/(K*kg)
    eo = 611        # Pa
    L = MV2.ones(T.shape) * 2.502E6   # J/kg
    L = MV2.where(MV2.less(T, To), 2.83E6, L)   # J/kg
    es = Es(T, method=method)  # in Pa
    e = rh * es / 100.
    Td = 1. / (1. / To - Rv / L * (e / eo))
    Td.setAxisList(T.getAxisList())
    return Td


class Gth(object):

    def TP2XY(self, T, P):
        x, y = self._TP2XY(T, P)
        try:
            x = float(x)
            y = float(y)
        except:
            pass
        return x, y

    def XY2TP(self, X, Y):
        t, p = self._XY2TP(X, Y)
        try:
            t = float(t)
            p = float(p)
        except:
            pass
        return t, p

    def setdiagramdefs(self):
        # Figures out which kind of diagram we're drawing
        if self.type[:2] == 'sk':
            self._TP2XY = self.TP2XYskewT
            self._XY2TP = self.XY2TPskewT
        elif self.type[:2] == 'em':
            self._TP2XY = self.TP2XYemagram
            self._XY2TP = self.XY2TPemagram
        elif self.type[:2] == 'te':
            self._TP2XY = self.TP2XYtephigram
            self._XY2TP = self.XY2TPtephigram
        elif self.type[:2] == 'st':
            self._TP2XY = self.TP2XYstuve
            self._XY2TP = self.XY2TPstuve
        elif self.type[:2] == 'cu':
            # user defined diagram
            pass
        else:
            raise ValueError('Error type: ' + repr(self.type) + ' is not a valid diagram type')

        # Gets the xs for each corner
        dwx1 = self.datawc_x1 + 273.16
        dwx2 = self.datawc_x2 + 273.16
        dwy1 = self.datawc_y1 * 100.
        dwy2 = self.datawc_y2 * 100.
        self.xmin, self.ymin = self.TP2XY(dwx1, dwy1)
        self.xmax, dum = self.TP2XY(dwx2, dwy1)
        dum, self.ymax = self.TP2XY(dwx2, dwy2)

    def __del__(self):
        self.x.removeobject(self._isotherms)
        self.x.removeobject(self._isothermsfilled)
        self.x.removeobject(self._isobars)
        self.x.removeobject(self._dryadiabats)
        self.x.removeobject(self._pseudoadiabats)
        self.x.removeobject(self._mixingratio)
        self.x.removeobject(self._line)

    def __init__(self, x=None, type='skewT', name='default'):
        self.displays = []
        if x is None:
            self.x = vcs.init()
        else:
            self.x = x

# if name in thermo_objects:
#             raise Exception,"Error thermo object %s already exists" % name
# thermo_objects.append(name)
        self.gname = 'Gth'
        self._name = name
        # Set some constants
        self.k = 0.286
        self.R = 287.04
        self.cp = 1004.
# self.L=2.5E3

        # World coordiantes in C/hPa
        self._datawc_x1 = -25.
        self._datawc_x2 = 50.
        self._datawc_y1 = 1000.
        self._datawc_y2 = 100.

        # Skewness (for skewT)
        self._skewness = -35.

        # Sets the diagram type
        self._type = type
        self.setdiagramdefs()

        # Level of detail
        self._detail = 25

        # Draw mixing ratio and pseudoadiabats up to
        self._Pmaxmixingratio = 200

        # Do we draw all the lines on the diagram ? (default, yes)
        self._drawisothermsfilled = 1
        self._drawisotherms = 1
        self._drawisobars = 1
        self._drawdryadiabats = 1
        self._drawpseudoadiabats = 1
        self._drawmixingratio = 1

        # Creates the isolines
        self._isotherms = self.x.createisoline()
        self._isobars = self.x.createisoline()
        self._dryadiabats = self.x.createisoline()
        self._pseudoadiabats = self.x.createisoline()
        self._mixingratio = self.x.createisoline()
        self._isothermsfilled = self.x.createisofill()

        # Isofill stuff for temperatures
        self._isothermsfilled.levels = [[-200, -190], [-180, -170], [-160, -150],
                                        [-140, -130], [-120, -110], [-100, -90],
                                        [-80, -70], [-60, -50], [-40, -30],
                                        [-20, -10], [0, 10], [20, 30], [40, 50],
                                        [60, 70], [80, 90], [100, 110],
                                        [120, 130], [140, 150], [160, 170],
                                        [180, 190], [200, 210]]
#         self.isothermsfilled.levels=([-240.0, -200.0], [-160.0, -120.0], [-80.0, -40.0], [0.0, 40.0], [80.0, 120.0])
        self._isothermsfilled.fillareacolors = ["light green", ] * 50
        self._isothermsfilled.xticlabels1 = {}
        self._isothermsfilled.yticlabels1 = {}

        # Isotherms
        self._isotherms.level = list(range(-200, 200, 10))
        self._isotherms.linecolors = ['brown', ]
        self._isotherms.linetypes = ['solid']
        self._isotherms.label = 'y'

        # Isobars default settings
        self._isobars.linecolors = ['brown']
        self._isobars.linetypes = ['solid']
        self._isobars.label = 'n'

        # Dry adiabatics
        self._dryadiabats.label = 'y'
        self._dryadiabats.linecolors = ['brown', ]
        self._dryadiabats.linetypes = ['solid']
# self._dryadiabats.level=range(-40,170,10)

        # Pseudoadiabatics
        self._pseudoadiabats.linecolors = ["DarkGreen", ]
        self._pseudoadiabats.linetypes = ['solid']
        self._pseudoadiabats.label = 'y'
# self._pseudoadiabats.level=vcs.mkevenlevels(-40,40,40)

        # Mixing ratio
        self._mixingratio.linetypes = ['long-dash']
        self._mixingratio.linecolors = ["DarkGreen", ]
        self._mixingratio.level = [.1, .2, .4, .6, .8, 1., 1.5, 2., 2.5,
                                   3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20,
                                   24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 68]
        self._mixingratio.label = 'y'

        # Windbarb
        self._windbarbsscales = [25., 10., 5.]

        # Lines for plot
        self._line = self.x.createline()

    def list(self):
        print('----------ThermodynamicDiagram (Gth) member (attribute) listings----------')
        print('graphic method = Gth')
        print('name =', self.name)
        print('type =', self.type)
        if self.type[:2] == 'sk':
            print('skewness:', self.skewness)
        print('detail =', self.detail)
        print('datawc_x1 =', self.datawc_x1)
        print('datawc_x2 =', self.datawc_x2)
        print('datawc_y1 =', self.datawc_y1)
        print('datawc_y2 =', self.datawc_y2)
        print('Pmaxmixingratio =', self.Pmaxmixingratio)
        print('isotherms =', self.isotherms.name)
        print('drawisotherms =', self.drawisotherms)
        print('isothermsfilled =', self.isothermsfilled.name)
        print('drawisothermsfilled =', self.drawisothermsfilled)
        print('isobars =', self.isobars.name)
        print('drawisobars =', self.drawisobars)
        print('dryadiabats =', self.dryadiabats.name)
        print('drawdryadiabats =', self.drawdryadiabats)
        print('pseudoadiabats =', self.pseudoadiabats.name)
        print('drawpseudoadiabats =', self.drawpseudoadiabats)
        print('mixingratio =', self.mixingratio.name)
        print('drawmixingratio =', self.drawmixingratio)
# print 'windbarbs =',self.windbarbs.name
        print('windbarbsscales =', self.windbarbsscales)
        print('linecolor', self.line.color[0])
        print('linewidth', self.line.width[0])
        print('linetype', self.line.type[0])

    __slots__ = [
        '_TP2XY',
        '_XY2TP',
        # 'TP2XYskewT',
        # 'TP2XYemagram',
        # 'TP2XYtephigram',
        # 'TP2XYstuve',
        'gname',
        '_name',
        '_displays',
        '_detail',
        'R',
        'cp',
        'k',
        '_datawc_x1',
        '_datawc_x2',
        '_datawc_y1',
        '_datawc_y2',
        '_type',
        'xmin',
        'xmax',
        'ymin',
        'ymax',
        'x',
        '_Pmaxmixingratio',
        '_isotherms',
        '_drawisotherms',
        '_isothermsfilled',
        '_drawisothermsfilled',
        '_isobars',
        '_drawisobars',
        '_dryadiabats',
        '_drawdryadiabats',
        '_pseudoadiabats',
        '_drawpseudoadiabats',
        '_mixingratio',
        '_drawmixingratio',
        '_windbarbs',
        '_windbarbsscales',
        '_skewness',
        '_line',
        '_linecolor',
        '_linewidth',
        '_linetype',
    ]

    def _getskewness(self):
        return self._skewness

    def _setskewness(self, value):
        VCS_validation_functions.checkName(self, 'skewness', value)
        if self.type[:2] != 'sk':
            raise ValueError('skewness can only be set for skewT type diagrams')
        value = VCS_validation_functions.checkIntFloat(self, 'skewness', value)
        self._skewness = value
        self.setdiagramdefs()
    skewness = property(_getskewness, _setskewness, None, 'Skewness of skewT diagram')

    def _getwindbarbsscales(self):
        return self._windbarbsscales

    def _setwindbarbsscales(self, value):
        value = VCS_validation_functions.checkListOfNumbers(
            self, 'windbarbsscales', value, minelements=3, maxelements=3)
        self._windbarbsscales = value
    windbarbsscales = property(_getwindbarbsscales, _setwindbarbsscales, None,
                               'WindBarbs scaless: Triangles, Full Barbs, Half Barb\ndefault is [25,10,5]')

    def _getdrawpseudoadiabats(self):
        return self._drawpseudoadiabats

    def _setdrawpseudoadiabats(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawpseudoadiabats', value)
        self._drawpseudoadiabats = value
    drawpseudoadiabats = property(_getdrawpseudoadiabats, _setdrawpseudoadiabats, None,
                                  'Draw the pseudo-adiabats ? (1/0) or (y/n)\nAlso called moist-adiabats')

    def _getdrawmixingratio(self):
        return self._drawmixingratio

    def _setdrawmixingratio(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawmixingratio', value)
        self._drawmixingratio = value
    drawmixingratio = property(_getdrawmixingratio, _setdrawmixingratio, None,
                               'Draw the pseudo-adiabats ? (1/0) or (y/n)\nAlso called moist-adiabats')

    def _getdrawdryadiabats(self):
        return self._drawdryadiabats

    def _setdrawdryadiabats(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawdryadiabats', value)
        self._drawdryadiabats = value
    drawdryadiabats = property(_getdrawdryadiabats, _setdrawdryadiabats, None, 'Draw the dryadiabats ? (1/0) or (y/n)')

    def _getdrawisobars(self):
        return self._drawisobars

    def _setdrawisobars(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawisobars', value)
        self._drawisobars = value
    drawisobars = property(_getdrawisobars, _setdrawisobars, None, 'Draw the isobars ? (1/0) or (y/n)')

    def _getdrawisothermsfilled(self):
        return self._drawisothermsfilled

    def _setdrawisothermsfilled(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawisothermsfilled', value)
        self._drawisothermsfilled = value
    drawisothermsfilled = property(_getdrawisothermsfilled, _setdrawisothermsfilled,
                                   None, 'Draw the isothermsfilled ? (1/0) or (y/n)')

    def _getdrawisotherms(self):
        return self._drawisotherms

    def _setdrawisotherms(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawisotherms', value)
        self._drawisotherms = value
    drawisotherms = property(_getdrawisotherms, _setdrawisotherms, None, 'Draw the isotherms ? (1/0) or (y/n)')

    def _getline(self):
        return self._line

    def _setline(self, value):
        value = VCS_validation_functions.checkLine(self, 'line', value)
        if isinstance(value, str):
            self._line = self.x.getline(value)
        else:
            self._line = value
    line = property(
        _getline,
        _setline,
        None,
        'Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getlinecolor(self):
        return self._line.color

    def _setlinecolor(self, value):
        value = VCS_validation_functions.checkColor(self, 'linecolor', value)
        self._line.color = [value]
    linecolor = property(
        _getlinecolor,
        _setlinecolor,
        None,
        'Color Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getlinewidth(self):
        return self._line.width

    def _setlinewidth(self, value):
        if not isinstance(value, (int, float)):
            raise 'Error width must be a number'
        if not 1 < value < 300:
            raise 'Error width must be between 1 and 300'
        self._line.width = [value]
    linewidth = property(
        _getlinewidth,
        _setlinewidth,
        None,
        'Width Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getlinetype(self):
        return self._line.type

    def _setlinetype(self, value):
        self._line.type = [value]
    linetype = property(
        _getlinetype,
        _setlinetype,
        None,
        'Type Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getwindbarbs(self):
        return self._windbarbs

    def _setwindbarbs(self, value):
        value = VCS_validation_functions.checkLine(self, 'windbarbs', value)
        self._windbarbs = value
    windbarbs = property(
        _getwindbarbs,
        _setwindbarbs,
        None,
        'Line graphic method to use to draw the windbarbs\nUse windbarbs.list() to see attributes')

    def _getdryadiabats(self):
        return self._dryadiabats

    def _setdryadiabats(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'dryadiabats', value)
        self._dryadiabats = value
    dryadiabats = property(
        _getdryadiabats,
        _setdryadiabats,
        None,
        'Isoline graphic method to use to draw the dryadiabats\nUse dryadiabats.list() to see attributes')

    def _getpseudoadiabats(self):
        return self._pseudoadiabats

    def _setpseudoadiabats(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'pseudoadiabats', value)
        self._pseudoadiabats = value
    pseudoadiabats = property(
        _getpseudoadiabats,
        _setpseudoadiabats,
        None,
        'Isoline graphic method to use to draw the pseudoadiabats\nUse pseudoadiabats.list() to see attributes')

    def _getmixingratio(self):
        return self._mixingratio

    def _setmixingratio(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'mixingratio', value)
        self._mixingratio = value
    mixingratio = property(
        _getmixingratio,
        _setmixingratio,
        None,
        'Isoline graphic method to use to draw the mixingratio\nUse mixingratio.list() to see attributes')

    def _getisobars(self):
        return self._isobars

    def _setisobars(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'isobars', value)
        self._isobars = value
    isobars = property(
        _getisobars,
        _setisobars,
        None,
        'Isoline graphic method to use to draw the isobars\nUse isobars.list() to see attributes')

    def _getisothermsfilled(self):
        return self._isothermsfilled

    def _setisothermsfilled(self, value):
        value = VCS_validation_functions.checkIsofill(self, 'isothermsfilled', value)
        self._isothermsfilled = value
    isothermsfilled = property(
        _getisothermsfilled,
        _setisothermsfilled,
        None,
        'Isofill graphic method to use to draw the isothermsfilled\nUse isothermsfilled.list() to see attributes')

    def _getisotherms(self):
        return self._isotherms

    def _setisotherms(self, value):
        value = VCS_validation_functions.checkIsoline(value)
        self._isotherms = value
    isotherms = property(
        _getisotherms,
        _setisotherms,
        None,
        'Isoline graphic method to use to draw the isotherms\nUse isotherms.list() to see attributes')

    def _getPmaxmixingratio(self):
        return self._Pmaxmixingratio

    def _setPmaxmixingratio(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'Pmaxmixingratio', value)
        self._Pmaxmixingratio = value
    Pmaxmixingratio = property(_getPmaxmixingratio, _setPmaxmixingratio, None,
                               'Pressure after which mixing ratio and pseudoadiats curves are not drawn anymore (hPa)')

    def _getdetail(self):
        return self._detail

    def _setdetail(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'detail', value)
        self._detail = value
    detail = property(
        _getdetail,
        _setdetail,
        None,
        'Smoothness of contours/graphics\nhigher value means smoother (and slower) plots')

    def _gettype(self):
        return self._type

    def _settype(self, value):
        VCS_validation_functions.checkName(self, 'type', value)
        if not isinstance(value, str):
            raise ValueError('Type should be a string')
        value = value.lower()
        if not value[:2] in ['sk', 'te', 'st', 'cu', 'em']:
            raise ValueError('Error Thermodynamic Diagram type must skewT, Tephigram, Emmagram, Stuve or Custom')
        self._type = value
        self.setdiagramdefs()
    type = property(
        _gettype,
        _settype,
        None,
        'Thermodynamic Diagram type, one of: skewT, Tephigram, Emmagram, Stuve or Custom')

    def _getname(self):
        return self._name

    def _setname(self, value):
        value = VCS_validation_functions.checkname(self, 'name', value)
        self._name = value
    name = property(_getname, _setname, None, 'Name of the Graphic method')

    def _getdatawc_x1(self):
        return self._datawc_x1

    def _setdatawc_x1(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_x1', value)
        self._datawc_x1 = value
        self.setdiagramdefs()
    datawc_x1 = property(_getdatawc_x1, _setdatawc_x1, None, 'Temperature of lower left corner (in C)')

    def _getdatawc_x2(self):
        return self._datawc_x2

    def _setdatawc_x2(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_x2', value)
        self._datawc_x2 = value
        self.setdiagramdefs()
    datawc_x2 = property(_getdatawc_x2, _setdatawc_x2, None, 'Temperature of lower right corner (in C)')

    def _getdatawc_y1(self):
        return self._datawc_y1

    def _setdatawc_y1(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_y1', value)
        self._datawc_y1 = value
        self.setdiagramdefs()
    datawc_y1 = property(_getdatawc_y1, _setdatawc_y1, None, 'Minimum Pressure (in hPa)')

    def _getdatawc_y2(self):
        return self._datawc_y2

    def _setdatawc_y2(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_y2', value)
        self._datawc_y2 = value
        self.setdiagramdefs()
    datawc_y2 = property(_getdatawc_y2, _setdatawc_y2, None, 'Maximum Pressure (in hPa)')

    def _getdisplays(self):
        return self._displays

    def _setdisplays(self, value):
        if not isinstance(value, list):
            raise Exception("displays must be a list of displays")
        for v in value:
            if not isinstance(v, vcs.displayplot.Dp):
                raise Exception("displays must be a list of displays")
        self._displays = value
    displays = property(_getdisplays, _setdisplays, None, "List of displays")

    def plot_pseudo(self, template, bg):
        """
        ## Ok this is pathetic but i can't get it to work any other way!!!
        ## We'll have to do the integration for each levels
        """
        levs = self.pseudoadiabats.level
        levels = []
        for l in levs:
            levels.append(l[0])
        # Ok we now have the required levels
        xs = []
        ys = []
        xstxt = []
        ystxt = []
        txt = []
        for il in range(len(levels)):
            l = levels[il]
            tmpxs = []
            tmpys = []
            t, p = LiftWet(l + 273.15, self.datawc_y1 * 100., self.Pmaxmixingratio * 100., 25000. / self.detail)
            nt = len(t)
            if il % 3 == 0:
                crit1 = .05
                crit2 = .95
            elif il % 3 == 1:
                crit1 = .2
                crit2 = .9
            elif il % 3 == 2:
                crit1 = .32
                crit2 = .85
            for i in range(nt):
                tx, ty = self.TP2XY(t[i], p[i])
                tmpxs.append(tx)
                tmpys.append(ty)
                if ((i - 1.) / nt < crit1 < i / float(nt)) or ((i - 1.) / nt < crit2 < i / float(nt)):
                    xstxt.append(tx)
                    ystxt.append(ty)
                    txt.append(str(l))
            xs.append(tmpxs)
            ys.append(tmpys)
        l = self.x.createline()
        l.color = [self.pseudoadiabats.linecolors[0],]
        l.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
        l.worldcoordinate = [self.xmin, self.xmax, self.ymin, self.ymax]
        l.x = xs
        l.y = ys
        self.displays.append(self.x.plot(l, bg=bg))
        t = self.x.createtext()
        t.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
        t.worldcoordinate = [self.xmin, self.xmax, self.ymin, self.ymax]
        t.color = l.color[0]
        t.x = xstxt
        t.y = ystxt
        t.string = txt
        t.height = 10
        self.displays.append(self.x.plot(t, bg=bg))

    def clear(self):
        self.displays = []
        self.x.clear()

    def plot_TP(self, T, P=None, template=None, bg=0):
        if template is None:
            template = 'thdiags_template'
        elif not isinstance(template, str):
            template = template.name
        try:
            isotemplate = self.x.createtemplate(template)
        except:
            isotemplate = self.x.gettemplate(template)
        self.plotIsos(template, bg=bg)
        if P is None:
            P = getPvalues(T)

        while T.rank() > 1:
            T = T[0]
        xs = []
        ys = []
        for i in range(len(T)):
            t = T[i]
            p = P[i]
            tx, ty = self.TP2XY(t, p)
            try:
                tx = float(tx)
                ty = float(ty)
                if numpy.ma.is_masked(tx) or numpy.ma.is_masked(ty) or numpy.isnan(
                        tx) or numpy.isnan(ty) or numpy.isinf(tx) or numpy.isinf(ty):
                    raise ValueError
                xs.append(tx)
                ys.append(ty)
            except:
                pass
# if not(MV2.isMA(tx) or MV2.isMA(ty)):
# xs.append(tx)
# ys.append(ty)
# else:
# print tx,ty,t,p
# xs.append(tx)
# ys.append(ty)
        if xs == [] and ys == []:
            for d in self.displays:
                d.off = 1
                self.x.remove_display_name(d.name)
            for i in range(len(self.displays)):
                self.displays.pop(-1)

            raise Exception("Error: data are all out of range!")
        iso = self.x.createline()
        iso.type = self.line.type
        iso.width = self.line.width
        iso.color = self.line.color

        iso.worldcoordinate = [self.xmin, self.xmax, self.ymin, self.ymax]
        iso.viewport = [isotemplate.data.x1, isotemplate.data.x2,
                        isotemplate.data.y1, isotemplate.data.y2]
        iso.x = xs
        iso.y = ys
        iso.priority = 3
        self.displays.append(self.x.plot(iso, bg=bg))
        return

    def setIso(self, iso, values=None):
        iso.datawc_x1 = self.xmin
        iso.datawc_x2 = self.xmax
        iso.datawc_y1 = self.ymin
        iso.datawc_y2 = self.ymax
        if values is not None:
            if iso.g_name == 'Gi':
                iso.level = values
                col = [iso.linecolors[0], ] * len(values)
                iso.linecolors = col
                iso.linetypes = [iso.linetypes[0], ] * len(values)
                if iso.text is None:
                    iso.text = ["default", ] * len(values)
                else:
                    iso.text = [iso.text[0], ] * len(values)
                iso.textcolors = col
            else:
                col = [iso.fillareacolors[0], ] * len(values)
                iso.fillareacolors = col
                iso.levels = values
        return

    def plotIsos(self, template=None, bg=0):
        if template is None:
            template = 'thdiags_template'
        elif not isinstance(template, type('')):
            template = template.name
        try:
            isotemplate = self.x.createtemplate(template)
        except:
            isotemplate = self.x.gettemplate(template)

        att = dir(isotemplate)
        for a in att:
            try:
                b = getattr(isotemplate, a)
                setattr(b, 'priority', 0)
            except:
                pass

        isotemplate.data.priority = 1
        isotemplate.box1.priority = 1
        isotemplate.box1.x1 = isotemplate.data.x1
        isotemplate.box1.x2 = isotemplate.data.x2
        isotemplate.box1.y1 = isotemplate.data.y1
        isotemplate.box1.y2 = isotemplate.data.y2

        dX = self.xmax - self.xmin
        dY = self.ymax - self.ymin
        X = numpy.ma.arange(self.detail, dtype=MV2.float)
        X = X * dX / (self.detail - 1) + self.xmin
        Xaxis = cdms2.createAxis(X)
        X = numpy.ma.resize(X, (int(self.detail), int(self.detail)))
        Y = numpy.ma.arange(self.detail, dtype=MV2.float)
        Y = Y * dY / (self.detail - 1) + self.ymin
        Yaxis = cdms2.createAxis(Y)
        Y = numpy.ma.resize(Y, (int(self.detail), int(self.detail)))
        Y = numpy.ma.transpose(Y)

        # Computes T,P on this grid
        T, P = self.XY2TP(X, Y)
        T = MV2.array(T)

        # Potential Temperature
        Tp = T / numpy.ma.power(P / 100000., self.k)

        ws = Ws(T, P)
        # Seems like we need not ws after 600mb
        ws = numpy.ma.masked_where(numpy.ma.less(P, self.Pmaxmixingratio * 100.), ws)

        T = T - 273.16
        Tmin, Tmax = vcs.minmax(T)

        Tvalues = self.isotherms.level
        if Tvalues == [[0.0, 1.0000000200408773e+20]]:
            Tvalues = vcs.mkscale(Tmin, Tmax)
            # Now sets the isothermsfilled
            Tvalues2 = []
            for i in range(len(Tvalues) / 2 - 1):
                t1, t2 = Tvalues[2 * i], Tvalues[2 * i + 1]
                if isinstance(t1, (list, tuple)):
                    t1, t2 = t1[0], t2[0]
                Tvalues2.append([t1, t2])
        else:
            Tvalues2 = self.isothermsfilled.levels
        self.setIso(self.isotherms, Tvalues)
# self.setIso(self.isothermsfilled,Tvalues2)

        P = P / 100.
        Pvalues = vcs.mkscale(self.datawc_y2, self.datawc_y1)
        self.setIso(self.isobars, Pvalues)

        Tp = Tp - 273.16
        Tpvalues = self.dryadiabats.level
        if Tpvalues == [[0.0, 1.0000000200408773e+20]]:
            min, max = vcs.minmax(Tp)
            Tpvalues = vcs.mkscale(min, max)
        self.setIso(self.dryadiabats, Tpvalues)

        # Pseudoadiabats
        Thevalues = self.pseudoadiabats.level
        if Thevalues == [[0.0, 1.0000000200408773e+20]]:
            Thevalues = vcs.mkevenlevels(-40, 40, 40)
        self.setIso(self.pseudoadiabats, Thevalues)

        # Mixing Ratio
        ws = ws * 1000.
        wsvalues = self.mixingratio.level
        if wsvalues == [[0.0, 1.0000000200408773e+20]]:
            min, max = vcs.minmax(ws)
            wsvalues = vcs.mkscale(min, max)
        self.setIso(self.mixingratio, wsvalues)

        # Figures out where to plot the P labels
        dicP = {}
        dicPl = vcs.mklabels(Pvalues)
        for p in Pvalues:
            X, Y = self.TP2XY(self.datawc_x1 + 273.15, p * 100)
            if not numpy.ma.isMA(X):
                dicP[Y] = dicPl[p]
        try:
            ttp = self.x.createtexttable('Plabels')
        except:
            ttp = self.x.gettexttable('Plabels')
        ttp.color = self.isobars.linecolors[0]
        isotemplate.ylabel1.texttable = ttp
        dicT = {}
        Tvalues = list(numpy.ma.array(Tvalues).filled().ravel())
        dicTl = vcs.mklabels(Tvalues)
        for t in Tvalues:
            X, Y = self.TP2XY(t + 273.15, self.datawc_y1 * 100)
            dicT[X] = dicTl[t]
        try:
            ttp = self.x.createtexttable('Tlabels')
        except:
            ttp = self.x.gettexttable('Tlabels')
        ttp.color = self.isotherms.linecolors[0]
        isotemplate.xlabel1.texttable = ttp
        isotemplate.ylabel1.priority = 1
        isotemplate.xlabel1.priority = 1
        isotemplate.data.x1 = isotemplate.data.x1
        isotemplate.data.x2 = isotemplate.data.x2
        isotemplate.data.y1 = isotemplate.data.y1
        isotemplate.data.y2 = isotemplate.data.y2
        self.isotherms.yticlabels1 = dicP
        self.isotherms.xticlabels1 = dicT
        self.isobars.yticlabels1 = {}
        self.isobars.xticlabels1 = {}
        self.dryadiabats.yticlabels1 = {}
        self.dryadiabats.xticlabels1 = {}
        self.pseudoadiabats.yticlabels1 = {}
        self.pseudoadiabats.xticlabels1 = {}
        self.mixingratio.yticlabels1 = {}
        self.mixingratio.xticlabels1 = {}
        self.isothermsfilled.datawc_x1 = self.isothermsfilled.datawc_x1
        self.isothermsfilled.datawc_x2 = self.isothermsfilled.datawc_x2
        self.isothermsfilled.datawc_y1 = self.isothermsfilled.datawc_y1
        self.isothermsfilled.datawc_x2 = self.isothermsfilled.datawc_y2
        # Puts the dims on it
        T.id = 'T'
        T.setAxis(0, Yaxis)
        T.setAxis(1, Xaxis)
        P = MV2.array(P)
        P.id = 'P'
        P.setAxis(0, Yaxis)
        P.setAxis(1, Xaxis)
        Tp = MV2.array(Tp)
        Tp.setAxis(0, Yaxis)
        Tp.setAxis(1, Xaxis)
        ws = MV2.array(ws)
        ws.setAxis(0, Yaxis)
        ws.setAxis(1, Xaxis)

        # plot if self.drawisothermsfilled:
        if self.drawisothermsfilled:
            self.displays.append(self.x.plot(T, isotemplate, self.isothermsfilled, bg=bg))
        if self.drawisotherms:
            self.displays.append(self.x.plot(T, isotemplate, self.isotherms, bg=bg))
        if self.drawisobars:
            self.displays.append(self.x.plot(P, isotemplate, self.isobars, bg=bg))
        if self.drawdryadiabats:
            self.displays.append(self.x.plot(Tp, isotemplate, self.dryadiabats, bg=bg))
        if self.drawpseudoadiabats:
            self.plot_pseudo(isotemplate, bg=bg)
        if self.drawmixingratio:
            self.displays.append(self.x.plot(ws, isotemplate, self.mixingratio, bg=bg))
        return

    # All take args with T in K and P in Pa
    # SkewT-LogP diagram defs
    def TP2XYskewT(self, T, P):
        X = T + self.skewness * numpy.ma.log(P)
        Y = -self.R * numpy.ma.log(P)
        return X, Y

    def XY2TPskewT(self, X, Y):
        P = numpy.ma.exp(-Y / self.R)
        T = X - self.skewness * numpy.ma.log(P)
        return T, P

    # Emagram defs
    def TP2XYemagram(self, T, P):
        X = T
        Y = -self.R * numpy.ma.log(P)
        return X, Y

    def XY2TPemagram(self, X, Y):
        T = X
        P = numpy.ma.exp(-Y / self.R)
        return T, P

    # Tephigram defs
    def TP2XYtephigram(self, T, P):
        X = T
        Tp = T / numpy.ma.power(P / 100000., self.k)
        Y = self.cp * numpy.ma.log(Tp)
        return X, Y

    def XY2TPtephigram(self, X, Y):
        T = X
        Tp = numpy.ma.exp(Y / self.cp)
        P = 100000. * numpy.ma.power(T / Tp, 1. / self.k)
        return T, P

    # Stuve diagram defs
    def TP2XYstuve(self, T, P):
        X = T
        Y = numpy.ma.power(P, self.k)
        return X, Y

    def XY2TPstuve(self, X, Y):
        T = X
        P = numpy.ma.power(Y, 1. / self.k)
        return T, P

    def plot_windbarb(self, u, v, P=None, template=None, type='uv', bg=0):
        if type == 'uv':
            n = numpy.ma.sqrt(u * u + v * v)
            d = numpy.ma.arccos(u / n)
            d = numpy.ma.where(numpy.ma.less(v, 0.), -d, d)

        if P is None:
            P = u.getLevel()
            if P is None:
                P = u.getAxis(-1)
            P = P.clone()
            try:  # tries to convert to Pa
                for i in range(len(P)):
                    uni = unidata.udunits(P[i], P.units)
                    P[i] = uni.to("Pa").value
            except Exception:
                pass
            P = P[:]

        n1 = n / self.windbarbsscales[0]
        n1 = n1.astype('i')
        n2 = (n - n1 * self.windbarbsscales[0]) / self.windbarbsscales[1]
        n2 = n2.astype('i')
        n3 = (n - n1 * self.windbarbsscales[0] - n2 * self.windbarbsscales[1]) / self.windbarbsscales[2]
        n3 = n3.astype('i')

# print 'n1:',n1
# print 'n2:',n2
# print 'n3:',n3

        nitems = len(P)
        for i in range(nitems):
            if not(MV2.isMA(n[i])):
                p = P[i]
                # Figure out the altitude to plot it (y coord) !
                dum, Y = self.TP2XY(273, p)
                lin = self.x.createline()
                lin.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
                lin.worldcoordinate = [-1, 1, -1, 1]
                lin.x = [0, 0]
                lin.y = [-1, 1]
                self.displays.append(self.x.plot(lin, bg=bg))
                lin = self.x.createline()
                lin.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
                lin.worldcoordinate = [-1, 1, self.ymin, self.ymax]
                # Set everything correctly y wise
                rw = 2 / (self.ymax - self.ymin)
                if self.x.islandscape():
                    r = 1.375
                else:
                    r = .72727272
                rv = (template.data.x2 - template.data.x1) / (template.data.y2 - template.data.y1)
                rw = 1. / rw
                x, y = self.make_barb(n[i], d[i], n1[i], n2[i], n3[i], rw * r * rv, Y)
                lin.x = x
                lin.y = y
                lin.type = ['solid']
                self.displays.append(self.x.plot(lin, bg=bg))

    def make_barb(self, n, d, n1, n2, n3, ydeformation, yoff):
        """ Draws a wind barb, the head is at location 0,0, the tail length is 1


        ---------------------------------------------------------------
        \  |  /   /    /
         \ | /   /    /
          \|/   /
           V

           total length is 1
           barbs angle is 60 degrees
           barb depth is .1

        """
        # Intitalize the Drawing, head is at zero
        xs = [0]
        ys = [0]
        xs.append(numpy.cos(d))
        ys.append(numpy.sin(d))

        # Ok now does the BIG barbs
        bangle = 25. / 180. * numpy.pi  # 30 degrees in radians
        bdepth = .3
        bwidth = numpy.tan(bangle) * bdepth

        hinit = 1.
        # Drawing of big barbs
        for i in range(n1):
            # summit
            l = hinit - bwidth  # legth to the middle of big barb
            h = numpy.sqrt(bdepth ** 2 + l ** 2)
            mu = numpy.arctan(bdepth / l)
            ys.append(h * numpy.sin(d + mu))
            xs.append(h * numpy.cos(d + mu))

            # End of barb
            ys.append((l - bwidth) * numpy.sin(d))
            xs.append((l - bwidth) * numpy.cos(d))
            hinit = hinit - 2 * bwidth
        if n1 != 0:
            hinit = hinit + bwidth / 2.

        # Ok now the long sticks
        for i in range(n2):
            l = hinit - bwidth
            # Move to begining of next one
            ys.append(l * numpy.sin(d))
            xs.append(l * numpy.cos(d))

            # Go to the top one
            l = hinit
            h = numpy.sqrt(bdepth ** 2 + l ** 2)
            mu = numpy.arctan(bdepth / l)
            ys.append(h * numpy.sin(d + mu))
            xs.append(h * numpy.cos(d + mu))

            # Back to begining
            ys.append(ys[-2])
            xs.append(xs[-2])

            # Adjust hinit
            hinit = hinit - bwidth / 2.

        # Now the half barbs
        bwidth = bwidth / 2.
        bdepth = bdepth / 2.
        for i in range(n3):
            l = hinit - bwidth
            # Move to begining of next one
            ys.append(l * numpy.sin(d))
            xs.append(l * numpy.cos(d))

            # Go to the top one
            l = hinit
            h = numpy.sqrt(bdepth ** 2 + l ** 2)
            mu = numpy.arctan(bdepth / l)
            ys.append(h * numpy.sin(d + mu))
            xs.append(h * numpy.cos(d + mu))
            # Back to begining
            ys.append(ys[-2])
            xs.append(xs[-2])

            # Adjust hinit
            hinit = hinit - bwidth / 2.

        # Scales and offset
        for i in range(len(ys)):
            # print 'i,x,y:',i,xs[i],ys[i]
            ys[i] *= ydeformation
            ys[i] += yoff

        return xs, ys

# Adapted for numpy/ma/cdms2 by convertcdms.py
import MV2
import cdms2
import genutil
import unidata
import vcs
import numpy
from vcs import VCS_validation_functions
thermo_objects = []


def Es(T, method=None):
    """Computes saturated pressure in Pa given T in K, using the method:
    1: Hyland-Wexler formulation, polynomial coeff (absolute norm)
    2: Wexler formulation
    3: Hyland-Wexler formulation, polynomial coeff (relative error norm)
    4: classic Goff Gratch equation
    5: 6.112*numpy.ma.exp(17.67*tempc/(tempc+243.5))

    Default is method 1

    Note: 1 and 2 use method 3 where T is not : 173.15 < T < 473.15
    ref for 1, 2 and 3:
    <NAME>., Journal of Applied Met., Vol 31, Dec 1992
    ( http://ams.allenpress.com/perlserv/?request=get-document&\
doi=10.1175%2F1520-0450(1992)031%3C1507%3APFTSVP%3E2.0.CO%3B2&ct=1 )
    """
    if method is None:
        method = 1

    if method == 1:
        # Put in C
        x = T - 273.15
        # Water vapor
        c0 = 0.611220713E03
        c1 = 0.443944344E02
        c2 = 0.143195336E01
        c3 = 0.263350515E-01
        c4 = 0.310636053E-03
        c5 = 0.185218710E-05
        c6 = 0.103440324E-07
        c7 = -0.468258100E-10
        c8 = 0.466533033E-13
        eswat = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (c5 + x * (c6 +
                                                                            x * (c7 + x * c8)))))))
        # ice
        c0 = .611153246E03
        c1 = .503261230E02
        c2 = .188595709E01
        c3 = .422115970E-01
        c4 = .620376691E-03
        c5 = .616082536E-05
        c6 = .405172828E-07
        c7 = .161492905E-09
        c8 = .297886454E-12
        esice = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (c5 + x * (c6 +
                                                                            x * (c7 + x * c8)))))))
        # Combine
        es = MV2.where(MV2.less(T, 273.15), esice, eswat)
        # Overwrite values outside valid range with method 2
        mn, mx = genutil.minmax(T)
        if mn < 173.16 or mx > 473.15:
            es2 = Es(T, method=2)
            es = MV2.where(MV2.less(T, 173.16), es2, es)
            es = MV2.where(MV2.greater(T, 473.15), es2, es)
    elif method == 2:
        # over water
        g0 = -0.29912729E4
        g1 = -0.60170128E4
        g2 = 0.1887643854E2
        g3 = -0.28354721E-1
        g4 = 0.17838301E-4
        g5 = -0.84150417E-9
        g6 = 0.44412543E-12
        g7 = 0.2858487E1
        # over ice
        k0 = -0.58653696e4
        k1 = 0.2224103300E2
        k2 = 0.13749042E-1
        k3 = -0.34031775E-4
        k4 = 0.26967687E-7
        k5 = 0.6918651
        esice = (k0 + (k1 + k5 * MV2.log(T) + (k2 + (
                       k3 + k4 * T) * T) * T) * T) / T  # over ice
        eswat = (g0 + (g1 + (g2 + g7 * MV2.log(T) + (g3 + (g4 + (
                             g5 + g6 * T) * T) * T) * T) * T) * T) / T ** 2  # over water
        es = MV2.where(MV2.less(T, 273.15), esice, eswat)
        es = MV2.exp(es)
    elif method == 3:
        # Put in C
        x = T - 273.15
        # Water vapor
        c0 = 0.611213476E03
        c1 = 0.444007856E02
        c2 = 0.143064234E01
        c3 = 0.264461437E-01
        c4 = 0.305930558E-03
        c5 = 0.196237241E-05
        c6 = 0.892344772E-08
        c7 = -0.373208410E-10
        c8 = 0.209339997E-13
        eswat = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (
                          c5 + x * (c6 + x * (c7 + x * c8)))))))
        # ice
        c0 = .611123516E03
        c1 = .503109514E02
        c2 = .1888369801E01
        c3 = .420547422E-01
        c4 = .614396778E-03
        c5 = .602780717E-05
        c6 = .387940929E-07
        c7 = .149436277E-09
        c8 = .262655803E-12
        esice = c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (
                          c5 + x * (c6 + x * (c7 + x * c8)))))))
        # Combine
        es = MV2.where(MV2.less(T, 273.15), esice, eswat)
        # Overwrite values outside valid range with method 2
        mn, mx = genutil.minmax(T)
        if mn < 173.16 or mx > 473.15:
            es2 = Es(T, method=2)
            es = MV2.where(MV2.less(T, 173.16), es2, es)
            es = MV2.where(MV2.greater(T, 473.15), es2, es)
    elif method == 4:
        est = 101324.6  # Pa
        Ts = 373.16 / T
        a = -7.90298
        b = 5.02808
        c = -1.3816E-7
        d = 11.344
        f = 8.1328E-3
        h = -3.49149
        maxexp = int(numpy.log10(numpy.finfo(numpy.float).max))
        minexp = 1 - a
        es = a * (Ts - 1.)
        es = es + b * numpy.ma.log10(Ts)
        A = d * (1. - Ts)
        A = numpy.ma.masked_greater(A, maxexp)
        A = numpy.ma.masked_less(A, minexp)
        es = es + c * (numpy.ma.power(10, A) - 1.)
        A = h * (Ts - 1.)
        A = numpy.ma.masked_greater(A, maxexp)
        A = numpy.ma.masked_less(A, minexp)
        es = es + f * (numpy.ma.power(10, A) - 1.)
        es = est * numpy.ma.power(10, es)
    elif method == 5:
        tempc = T - 273.15
        es = 611.2 * numpy.ma.exp(17.67 * tempc / (tempc + 243.5))
    return es


def getPvalues(T):
    """computes pressure values from level axis,
    if not designated level then use last axis...
    """
    P = T.getLevel()
    if P is None:
        P = T.getAxis(-1)
    P = P.clone()
    try:  # tries to convert to Pa
        for i in range(len(P)):
            u = unidata.udunits(P[i], P.units)
            P[i] = u.to("Pa").value
    except Exception:
        pass
    return P[:]


def Ws(T, P=None, method=None):
    """ Compute saturated mixing ratio at T,P, cpomputing es using the method:
    T in K
    P in Pa
    """
    es = Es(T, method)
    if P is None:
        P = getPvalues(T)

    return .622 * es / (P - es)


def gammaw(T, rh, P=None, method=None):
    """
    Function to calculate the moist adiabatic lapse rate (K/Pa) based
    on the temperature, pressure, and rh of the environment.
    T in K
    rh in %
    P in Pa
    if P is not passed will be optaines from T level axis (or last axis)
    """
    if P is None:
        P = getPvalues(T)
    ws = Ws(T, P, method=method)
    w = rh * ws / 100
    Tv = T * (1.0 + 0.6 * w)
    latent = 1000 * (2502.2 - 2.43089 * (T - 273.15))
    A = 1.0 + latent * ws / (287 * T)
    B = 1.0 + 0.622 * latent * latent * ws / (1005 * 287 * T * T)
    Density = P / (287 * Tv)  # in SI
    lapse = (A / B) / (1005. * Density)
    return lapse


def LiftWet(To, Po, Pend, deltaP=1000., method=None):
    """
    Lift a parcel moist adiabatically from Po to Pend.
    Initial temperature is To in K.
    To in K
    Po in Pa
    detail # number to increment P slowly
    returns T at all P (every deltaP Pa), and P
    """
    if method is None:
        method = 5  # MUCH faster
    temp = [To, ]
    pres = [Po, ]
    while (pres[-1] - deltaP >= Pend):
        tmp = temp[-1] - deltaP * gammaw(temp[-1], 100, pres[-1] - deltaP / 2, method=method)
        temp.append(tmp)
        pres.append(pres[-1] - deltaP)
    return MV2.array(temp, id='T'), MV2.array(pres, id='P')


def Dewpoint(T, rh, method=None):
    """Computes dewpoint profile from T and rh
    Usage:
      Td = Dewpoint(T,rh,method=None)
    Where:
     T is in K
     rh is in % (relative humidity)
     method used to compute Es:
         1: Hyland-Wexler formulation, polynomial coeff (absolute norm) (default)
         2: Wexler formulation
         3: Hyland-Wexler formulation, polynomial coeff (relative error norm)
         4: classic Goff Gratch equation
         5: 6.112*numpy.ma.exp(17.67*tempc/(tempc+243.5))



    Assumes constant latent temperature (except over ice or water)
    """
    To = 273.15     # K
    Rv = 461.53     # J/(K*kg)
    eo = 611        # Pa
    L = MV2.ones(T.shape) * 2.502E6   # J/kg
    L = MV2.where(MV2.less(T, To), 2.83E6, L)   # J/kg
    es = Es(T, method=method)  # in Pa
    e = rh * es / 100.
    Td = 1. / (1. / To - Rv / L * (e / eo))
    Td.setAxisList(T.getAxisList())
    return Td


class Gth(object):

    def TP2XY(self, T, P):
        x, y = self._TP2XY(T, P)
        try:
            x = float(x)
            y = float(y)
        except:
            pass
        return x, y

    def XY2TP(self, X, Y):
        t, p = self._XY2TP(X, Y)
        try:
            t = float(t)
            p = float(p)
        except:
            pass
        return t, p

    def setdiagramdefs(self):
        # Figures out which kind of diagram we're drawing
        if self.type[:2] == 'sk':
            self._TP2XY = self.TP2XYskewT
            self._XY2TP = self.XY2TPskewT
        elif self.type[:2] == 'em':
            self._TP2XY = self.TP2XYemagram
            self._XY2TP = self.XY2TPemagram
        elif self.type[:2] == 'te':
            self._TP2XY = self.TP2XYtephigram
            self._XY2TP = self.XY2TPtephigram
        elif self.type[:2] == 'st':
            self._TP2XY = self.TP2XYstuve
            self._XY2TP = self.XY2TPstuve
        elif self.type[:2] == 'cu':
            # user defined diagram
            pass
        else:
            raise ValueError('Error type: ' + repr(self.type) + ' is not a valid diagram type')

        # Gets the xs for each corner
        dwx1 = self.datawc_x1 + 273.16
        dwx2 = self.datawc_x2 + 273.16
        dwy1 = self.datawc_y1 * 100.
        dwy2 = self.datawc_y2 * 100.
        self.xmin, self.ymin = self.TP2XY(dwx1, dwy1)
        self.xmax, dum = self.TP2XY(dwx2, dwy1)
        dum, self.ymax = self.TP2XY(dwx2, dwy2)

    def __del__(self):
        self.x.removeobject(self._isotherms)
        self.x.removeobject(self._isothermsfilled)
        self.x.removeobject(self._isobars)
        self.x.removeobject(self._dryadiabats)
        self.x.removeobject(self._pseudoadiabats)
        self.x.removeobject(self._mixingratio)
        self.x.removeobject(self._line)

    def __init__(self, x=None, type='skewT', name='default'):
        self.displays = []
        if x is None:
            self.x = vcs.init()
        else:
            self.x = x

# if name in thermo_objects:
#             raise Exception,"Error thermo object %s already exists" % name
# thermo_objects.append(name)
        self.gname = 'Gth'
        self._name = name
        # Set some constants
        self.k = 0.286
        self.R = 287.04
        self.cp = 1004.
# self.L=2.5E3

        # World coordiantes in C/hPa
        self._datawc_x1 = -25.
        self._datawc_x2 = 50.
        self._datawc_y1 = 1000.
        self._datawc_y2 = 100.

        # Skewness (for skewT)
        self._skewness = -35.

        # Sets the diagram type
        self._type = type
        self.setdiagramdefs()

        # Level of detail
        self._detail = 25

        # Draw mixing ratio and pseudoadiabats up to
        self._Pmaxmixingratio = 200

        # Do we draw all the lines on the diagram ? (default, yes)
        self._drawisothermsfilled = 1
        self._drawisotherms = 1
        self._drawisobars = 1
        self._drawdryadiabats = 1
        self._drawpseudoadiabats = 1
        self._drawmixingratio = 1

        # Creates the isolines
        self._isotherms = self.x.createisoline()
        self._isobars = self.x.createisoline()
        self._dryadiabats = self.x.createisoline()
        self._pseudoadiabats = self.x.createisoline()
        self._mixingratio = self.x.createisoline()
        self._isothermsfilled = self.x.createisofill()

        # Isofill stuff for temperatures
        self._isothermsfilled.levels = [[-200, -190], [-180, -170], [-160, -150],
                                        [-140, -130], [-120, -110], [-100, -90],
                                        [-80, -70], [-60, -50], [-40, -30],
                                        [-20, -10], [0, 10], [20, 30], [40, 50],
                                        [60, 70], [80, 90], [100, 110],
                                        [120, 130], [140, 150], [160, 170],
                                        [180, 190], [200, 210]]
#         self.isothermsfilled.levels=([-240.0, -200.0], [-160.0, -120.0], [-80.0, -40.0], [0.0, 40.0], [80.0, 120.0])
        self._isothermsfilled.fillareacolors = ["light green", ] * 50
        self._isothermsfilled.xticlabels1 = {}
        self._isothermsfilled.yticlabels1 = {}

        # Isotherms
        self._isotherms.level = list(range(-200, 200, 10))
        self._isotherms.linecolors = ['brown', ]
        self._isotherms.linetypes = ['solid']
        self._isotherms.label = 'y'

        # Isobars default settings
        self._isobars.linecolors = ['brown']
        self._isobars.linetypes = ['solid']
        self._isobars.label = 'n'

        # Dry adiabatics
        self._dryadiabats.label = 'y'
        self._dryadiabats.linecolors = ['brown', ]
        self._dryadiabats.linetypes = ['solid']
# self._dryadiabats.level=range(-40,170,10)

        # Pseudoadiabatics
        self._pseudoadiabats.linecolors = ["DarkGreen", ]
        self._pseudoadiabats.linetypes = ['solid']
        self._pseudoadiabats.label = 'y'
# self._pseudoadiabats.level=vcs.mkevenlevels(-40,40,40)

        # Mixing ratio
        self._mixingratio.linetypes = ['long-dash']
        self._mixingratio.linecolors = ["DarkGreen", ]
        self._mixingratio.level = [.1, .2, .4, .6, .8, 1., 1.5, 2., 2.5,
                                   3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20,
                                   24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 68]
        self._mixingratio.label = 'y'

        # Windbarb
        self._windbarbsscales = [25., 10., 5.]

        # Lines for plot
        self._line = self.x.createline()

    def list(self):
        print('----------ThermodynamicDiagram (Gth) member (attribute) listings----------')
        print('graphic method = Gth')
        print('name =', self.name)
        print('type =', self.type)
        if self.type[:2] == 'sk':
            print('skewness:', self.skewness)
        print('detail =', self.detail)
        print('datawc_x1 =', self.datawc_x1)
        print('datawc_x2 =', self.datawc_x2)
        print('datawc_y1 =', self.datawc_y1)
        print('datawc_y2 =', self.datawc_y2)
        print('Pmaxmixingratio =', self.Pmaxmixingratio)
        print('isotherms =', self.isotherms.name)
        print('drawisotherms =', self.drawisotherms)
        print('isothermsfilled =', self.isothermsfilled.name)
        print('drawisothermsfilled =', self.drawisothermsfilled)
        print('isobars =', self.isobars.name)
        print('drawisobars =', self.drawisobars)
        print('dryadiabats =', self.dryadiabats.name)
        print('drawdryadiabats =', self.drawdryadiabats)
        print('pseudoadiabats =', self.pseudoadiabats.name)
        print('drawpseudoadiabats =', self.drawpseudoadiabats)
        print('mixingratio =', self.mixingratio.name)
        print('drawmixingratio =', self.drawmixingratio)
# print 'windbarbs =',self.windbarbs.name
        print('windbarbsscales =', self.windbarbsscales)
        print('linecolor', self.line.color[0])
        print('linewidth', self.line.width[0])
        print('linetype', self.line.type[0])

    __slots__ = [
        '_TP2XY',
        '_XY2TP',
        # 'TP2XYskewT',
        # 'TP2XYemagram',
        # 'TP2XYtephigram',
        # 'TP2XYstuve',
        'gname',
        '_name',
        '_displays',
        '_detail',
        'R',
        'cp',
        'k',
        '_datawc_x1',
        '_datawc_x2',
        '_datawc_y1',
        '_datawc_y2',
        '_type',
        'xmin',
        'xmax',
        'ymin',
        'ymax',
        'x',
        '_Pmaxmixingratio',
        '_isotherms',
        '_drawisotherms',
        '_isothermsfilled',
        '_drawisothermsfilled',
        '_isobars',
        '_drawisobars',
        '_dryadiabats',
        '_drawdryadiabats',
        '_pseudoadiabats',
        '_drawpseudoadiabats',
        '_mixingratio',
        '_drawmixingratio',
        '_windbarbs',
        '_windbarbsscales',
        '_skewness',
        '_line',
        '_linecolor',
        '_linewidth',
        '_linetype',
    ]

    def _getskewness(self):
        return self._skewness

    def _setskewness(self, value):
        VCS_validation_functions.checkName(self, 'skewness', value)
        if self.type[:2] != 'sk':
            raise ValueError('skewness can only be set for skewT type diagrams')
        value = VCS_validation_functions.checkIntFloat(self, 'skewness', value)
        self._skewness = value
        self.setdiagramdefs()
    skewness = property(_getskewness, _setskewness, None, 'Skewness of skewT diagram')

    def _getwindbarbsscales(self):
        return self._windbarbsscales

    def _setwindbarbsscales(self, value):
        value = VCS_validation_functions.checkListOfNumbers(
            self, 'windbarbsscales', value, minelements=3, maxelements=3)
        self._windbarbsscales = value
    windbarbsscales = property(_getwindbarbsscales, _setwindbarbsscales, None,
                               'WindBarbs scaless: Triangles, Full Barbs, Half Barb\ndefault is [25,10,5]')

    def _getdrawpseudoadiabats(self):
        return self._drawpseudoadiabats

    def _setdrawpseudoadiabats(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawpseudoadiabats', value)
        self._drawpseudoadiabats = value
    drawpseudoadiabats = property(_getdrawpseudoadiabats, _setdrawpseudoadiabats, None,
                                  'Draw the pseudo-adiabats ? (1/0) or (y/n)\nAlso called moist-adiabats')

    def _getdrawmixingratio(self):
        return self._drawmixingratio

    def _setdrawmixingratio(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawmixingratio', value)
        self._drawmixingratio = value
    drawmixingratio = property(_getdrawmixingratio, _setdrawmixingratio, None,
                               'Draw the pseudo-adiabats ? (1/0) or (y/n)\nAlso called moist-adiabats')

    def _getdrawdryadiabats(self):
        return self._drawdryadiabats

    def _setdrawdryadiabats(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawdryadiabats', value)
        self._drawdryadiabats = value
    drawdryadiabats = property(_getdrawdryadiabats, _setdrawdryadiabats, None, 'Draw the dryadiabats ? (1/0) or (y/n)')

    def _getdrawisobars(self):
        return self._drawisobars

    def _setdrawisobars(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawisobars', value)
        self._drawisobars = value
    drawisobars = property(_getdrawisobars, _setdrawisobars, None, 'Draw the isobars ? (1/0) or (y/n)')

    def _getdrawisothermsfilled(self):
        return self._drawisothermsfilled

    def _setdrawisothermsfilled(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawisothermsfilled', value)
        self._drawisothermsfilled = value
    drawisothermsfilled = property(_getdrawisothermsfilled, _setdrawisothermsfilled,
                                   None, 'Draw the isothermsfilled ? (1/0) or (y/n)')

    def _getdrawisotherms(self):
        return self._drawisotherms

    def _setdrawisotherms(self, value):
        value = VCS_validation_functions.checkOnOff(self, 'drawisotherms', value)
        self._drawisotherms = value
    drawisotherms = property(_getdrawisotherms, _setdrawisotherms, None, 'Draw the isotherms ? (1/0) or (y/n)')

    def _getline(self):
        return self._line

    def _setline(self, value):
        value = VCS_validation_functions.checkLine(self, 'line', value)
        if isinstance(value, str):
            self._line = self.x.getline(value)
        else:
            self._line = value
    line = property(
        _getline,
        _setline,
        None,
        'Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getlinecolor(self):
        return self._line.color

    def _setlinecolor(self, value):
        value = VCS_validation_functions.checkColor(self, 'linecolor', value)
        self._line.color = [value]
    linecolor = property(
        _getlinecolor,
        _setlinecolor,
        None,
        'Color Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getlinewidth(self):
        return self._line.width

    def _setlinewidth(self, value):
        if not isinstance(value, (int, float)):
            raise 'Error width must be a number'
        if not 1 < value < 300:
            raise 'Error width must be between 1 and 300'
        self._line.width = [value]
    linewidth = property(
        _getlinewidth,
        _setlinewidth,
        None,
        'Width Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getlinetype(self):
        return self._line.type

    def _setlinetype(self, value):
        self._line.type = [value]
    linetype = property(
        _getlinetype,
        _setlinetype,
        None,
        'Type Line graphic method to use to draw the profile\nUse line.list() to see attributes')

    def _getwindbarbs(self):
        return self._windbarbs

    def _setwindbarbs(self, value):
        value = VCS_validation_functions.checkLine(self, 'windbarbs', value)
        self._windbarbs = value
    windbarbs = property(
        _getwindbarbs,
        _setwindbarbs,
        None,
        'Line graphic method to use to draw the windbarbs\nUse windbarbs.list() to see attributes')

    def _getdryadiabats(self):
        return self._dryadiabats

    def _setdryadiabats(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'dryadiabats', value)
        self._dryadiabats = value
    dryadiabats = property(
        _getdryadiabats,
        _setdryadiabats,
        None,
        'Isoline graphic method to use to draw the dryadiabats\nUse dryadiabats.list() to see attributes')

    def _getpseudoadiabats(self):
        return self._pseudoadiabats

    def _setpseudoadiabats(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'pseudoadiabats', value)
        self._pseudoadiabats = value
    pseudoadiabats = property(
        _getpseudoadiabats,
        _setpseudoadiabats,
        None,
        'Isoline graphic method to use to draw the pseudoadiabats\nUse pseudoadiabats.list() to see attributes')

    def _getmixingratio(self):
        return self._mixingratio

    def _setmixingratio(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'mixingratio', value)
        self._mixingratio = value
    mixingratio = property(
        _getmixingratio,
        _setmixingratio,
        None,
        'Isoline graphic method to use to draw the mixingratio\nUse mixingratio.list() to see attributes')

    def _getisobars(self):
        return self._isobars

    def _setisobars(self, value):
        value = VCS_validation_functions.checkIsoline(self, 'isobars', value)
        self._isobars = value
    isobars = property(
        _getisobars,
        _setisobars,
        None,
        'Isoline graphic method to use to draw the isobars\nUse isobars.list() to see attributes')

    def _getisothermsfilled(self):
        return self._isothermsfilled

    def _setisothermsfilled(self, value):
        value = VCS_validation_functions.checkIsofill(self, 'isothermsfilled', value)
        self._isothermsfilled = value
    isothermsfilled = property(
        _getisothermsfilled,
        _setisothermsfilled,
        None,
        'Isofill graphic method to use to draw the isothermsfilled\nUse isothermsfilled.list() to see attributes')

    def _getisotherms(self):
        return self._isotherms

    def _setisotherms(self, value):
        value = VCS_validation_functions.checkIsoline(value)
        self._isotherms = value
    isotherms = property(
        _getisotherms,
        _setisotherms,
        None,
        'Isoline graphic method to use to draw the isotherms\nUse isotherms.list() to see attributes')

    def _getPmaxmixingratio(self):
        return self._Pmaxmixingratio

    def _setPmaxmixingratio(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'Pmaxmixingratio', value)
        self._Pmaxmixingratio = value
    Pmaxmixingratio = property(_getPmaxmixingratio, _setPmaxmixingratio, None,
                               'Pressure after which mixing ratio and pseudoadiats curves are not drawn anymore (hPa)')

    def _getdetail(self):
        return self._detail

    def _setdetail(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'detail', value)
        self._detail = value
    detail = property(
        _getdetail,
        _setdetail,
        None,
        'Smoothness of contours/graphics\nhigher value means smoother (and slower) plots')

    def _gettype(self):
        return self._type

    def _settype(self, value):
        VCS_validation_functions.checkName(self, 'type', value)
        if not isinstance(value, str):
            raise ValueError('Type should be a string')
        value = value.lower()
        if not value[:2] in ['sk', 'te', 'st', 'cu', 'em']:
            raise ValueError('Error Thermodynamic Diagram type must skewT, Tephigram, Emmagram, Stuve or Custom')
        self._type = value
        self.setdiagramdefs()
    type = property(
        _gettype,
        _settype,
        None,
        'Thermodynamic Diagram type, one of: skewT, Tephigram, Emmagram, Stuve or Custom')

    def _getname(self):
        return self._name

    def _setname(self, value):
        value = VCS_validation_functions.checkname(self, 'name', value)
        self._name = value
    name = property(_getname, _setname, None, 'Name of the Graphic method')

    def _getdatawc_x1(self):
        return self._datawc_x1

    def _setdatawc_x1(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_x1', value)
        self._datawc_x1 = value
        self.setdiagramdefs()
    datawc_x1 = property(_getdatawc_x1, _setdatawc_x1, None, 'Temperature of lower left corner (in C)')

    def _getdatawc_x2(self):
        return self._datawc_x2

    def _setdatawc_x2(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_x2', value)
        self._datawc_x2 = value
        self.setdiagramdefs()
    datawc_x2 = property(_getdatawc_x2, _setdatawc_x2, None, 'Temperature of lower right corner (in C)')

    def _getdatawc_y1(self):
        return self._datawc_y1

    def _setdatawc_y1(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_y1', value)
        self._datawc_y1 = value
        self.setdiagramdefs()
    datawc_y1 = property(_getdatawc_y1, _setdatawc_y1, None, 'Minimum Pressure (in hPa)')

    def _getdatawc_y2(self):
        return self._datawc_y2

    def _setdatawc_y2(self, value):
        value = VCS_validation_functions.checkIntFloat(self, 'datawc_y2', value)
        self._datawc_y2 = value
        self.setdiagramdefs()
    datawc_y2 = property(_getdatawc_y2, _setdatawc_y2, None, 'Maximum Pressure (in hPa)')

    def _getdisplays(self):
        return self._displays

    def _setdisplays(self, value):
        if not isinstance(value, list):
            raise Exception("displays must be a list of displays")
        for v in value:
            if not isinstance(v, vcs.displayplot.Dp):
                raise Exception("displays must be a list of displays")
        self._displays = value
    displays = property(_getdisplays, _setdisplays, None, "List of displays")

    def plot_pseudo(self, template, bg):
        """
        ## Ok this is pathetic but i can't get it to work any other way!!!
        ## We'll have to do the integration for each levels
        """
        levs = self.pseudoadiabats.level
        levels = []
        for l in levs:
            levels.append(l[0])
        # Ok we now have the required levels
        xs = []
        ys = []
        xstxt = []
        ystxt = []
        txt = []
        for il in range(len(levels)):
            l = levels[il]
            tmpxs = []
            tmpys = []
            t, p = LiftWet(l + 273.15, self.datawc_y1 * 100., self.Pmaxmixingratio * 100., 25000. / self.detail)
            nt = len(t)
            if il % 3 == 0:
                crit1 = .05
                crit2 = .95
            elif il % 3 == 1:
                crit1 = .2
                crit2 = .9
            elif il % 3 == 2:
                crit1 = .32
                crit2 = .85
            for i in range(nt):
                tx, ty = self.TP2XY(t[i], p[i])
                tmpxs.append(tx)
                tmpys.append(ty)
                if ((i - 1.) / nt < crit1 < i / float(nt)) or ((i - 1.) / nt < crit2 < i / float(nt)):
                    xstxt.append(tx)
                    ystxt.append(ty)
                    txt.append(str(l))
            xs.append(tmpxs)
            ys.append(tmpys)
        l = self.x.createline()
        l.color = [self.pseudoadiabats.linecolors[0],]
        l.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
        l.worldcoordinate = [self.xmin, self.xmax, self.ymin, self.ymax]
        l.x = xs
        l.y = ys
        self.displays.append(self.x.plot(l, bg=bg))
        t = self.x.createtext()
        t.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
        t.worldcoordinate = [self.xmin, self.xmax, self.ymin, self.ymax]
        t.color = l.color[0]
        t.x = xstxt
        t.y = ystxt
        t.string = txt
        t.height = 10
        self.displays.append(self.x.plot(t, bg=bg))

    def clear(self):
        self.displays = []
        self.x.clear()

    def plot_TP(self, T, P=None, template=None, bg=0):
        if template is None:
            template = 'thdiags_template'
        elif not isinstance(template, str):
            template = template.name
        try:
            isotemplate = self.x.createtemplate(template)
        except:
            isotemplate = self.x.gettemplate(template)
        self.plotIsos(template, bg=bg)
        if P is None:
            P = getPvalues(T)

        while T.rank() > 1:
            T = T[0]
        xs = []
        ys = []
        for i in range(len(T)):
            t = T[i]
            p = P[i]
            tx, ty = self.TP2XY(t, p)
            try:
                tx = float(tx)
                ty = float(ty)
                if numpy.ma.is_masked(tx) or numpy.ma.is_masked(ty) or numpy.isnan(
                        tx) or numpy.isnan(ty) or numpy.isinf(tx) or numpy.isinf(ty):
                    raise ValueError
                xs.append(tx)
                ys.append(ty)
            except:
                pass
# if not(MV2.isMA(tx) or MV2.isMA(ty)):
# xs.append(tx)
# ys.append(ty)
# else:
# print tx,ty,t,p
# xs.append(tx)
# ys.append(ty)
        if xs == [] and ys == []:
            for d in self.displays:
                d.off = 1
                self.x.remove_display_name(d.name)
            for i in range(len(self.displays)):
                self.displays.pop(-1)

            raise Exception("Error: data are all out of range!")
        iso = self.x.createline()
        iso.type = self.line.type
        iso.width = self.line.width
        iso.color = self.line.color

        iso.worldcoordinate = [self.xmin, self.xmax, self.ymin, self.ymax]
        iso.viewport = [isotemplate.data.x1, isotemplate.data.x2,
                        isotemplate.data.y1, isotemplate.data.y2]
        iso.x = xs
        iso.y = ys
        iso.priority = 3
        self.displays.append(self.x.plot(iso, bg=bg))
        return

    def setIso(self, iso, values=None):
        iso.datawc_x1 = self.xmin
        iso.datawc_x2 = self.xmax
        iso.datawc_y1 = self.ymin
        iso.datawc_y2 = self.ymax
        if values is not None:
            if iso.g_name == 'Gi':
                iso.level = values
                col = [iso.linecolors[0], ] * len(values)
                iso.linecolors = col
                iso.linetypes = [iso.linetypes[0], ] * len(values)
                if iso.text is None:
                    iso.text = ["default", ] * len(values)
                else:
                    iso.text = [iso.text[0], ] * len(values)
                iso.textcolors = col
            else:
                col = [iso.fillareacolors[0], ] * len(values)
                iso.fillareacolors = col
                iso.levels = values
        return

    def plotIsos(self, template=None, bg=0):
        if template is None:
            template = 'thdiags_template'
        elif not isinstance(template, type('')):
            template = template.name
        try:
            isotemplate = self.x.createtemplate(template)
        except:
            isotemplate = self.x.gettemplate(template)

        att = dir(isotemplate)
        for a in att:
            try:
                b = getattr(isotemplate, a)
                setattr(b, 'priority', 0)
            except:
                pass

        isotemplate.data.priority = 1
        isotemplate.box1.priority = 1
        isotemplate.box1.x1 = isotemplate.data.x1
        isotemplate.box1.x2 = isotemplate.data.x2
        isotemplate.box1.y1 = isotemplate.data.y1
        isotemplate.box1.y2 = isotemplate.data.y2

        dX = self.xmax - self.xmin
        dY = self.ymax - self.ymin
        X = numpy.ma.arange(self.detail, dtype=MV2.float)
        X = X * dX / (self.detail - 1) + self.xmin
        Xaxis = cdms2.createAxis(X)
        X = numpy.ma.resize(X, (int(self.detail), int(self.detail)))
        Y = numpy.ma.arange(self.detail, dtype=MV2.float)
        Y = Y * dY / (self.detail - 1) + self.ymin
        Yaxis = cdms2.createAxis(Y)
        Y = numpy.ma.resize(Y, (int(self.detail), int(self.detail)))
        Y = numpy.ma.transpose(Y)

        # Computes T,P on this grid
        T, P = self.XY2TP(X, Y)
        T = MV2.array(T)

        # Potential Temperature
        Tp = T / numpy.ma.power(P / 100000., self.k)

        ws = Ws(T, P)
        # Seems like we need not ws after 600mb
        ws = numpy.ma.masked_where(numpy.ma.less(P, self.Pmaxmixingratio * 100.), ws)

        T = T - 273.16
        Tmin, Tmax = vcs.minmax(T)

        Tvalues = self.isotherms.level
        if Tvalues == [[0.0, 1.0000000200408773e+20]]:
            Tvalues = vcs.mkscale(Tmin, Tmax)
            # Now sets the isothermsfilled
            Tvalues2 = []
            for i in range(len(Tvalues) / 2 - 1):
                t1, t2 = Tvalues[2 * i], Tvalues[2 * i + 1]
                if isinstance(t1, (list, tuple)):
                    t1, t2 = t1[0], t2[0]
                Tvalues2.append([t1, t2])
        else:
            Tvalues2 = self.isothermsfilled.levels
        self.setIso(self.isotherms, Tvalues)
# self.setIso(self.isothermsfilled,Tvalues2)

        P = P / 100.
        Pvalues = vcs.mkscale(self.datawc_y2, self.datawc_y1)
        self.setIso(self.isobars, Pvalues)

        Tp = Tp - 273.16
        Tpvalues = self.dryadiabats.level
        if Tpvalues == [[0.0, 1.0000000200408773e+20]]:
            min, max = vcs.minmax(Tp)
            Tpvalues = vcs.mkscale(min, max)
        self.setIso(self.dryadiabats, Tpvalues)

        # Pseudoadiabats
        Thevalues = self.pseudoadiabats.level
        if Thevalues == [[0.0, 1.0000000200408773e+20]]:
            Thevalues = vcs.mkevenlevels(-40, 40, 40)
        self.setIso(self.pseudoadiabats, Thevalues)

        # Mixing Ratio
        ws = ws * 1000.
        wsvalues = self.mixingratio.level
        if wsvalues == [[0.0, 1.0000000200408773e+20]]:
            min, max = vcs.minmax(ws)
            wsvalues = vcs.mkscale(min, max)
        self.setIso(self.mixingratio, wsvalues)

        # Figures out where to plot the P labels
        dicP = {}
        dicPl = vcs.mklabels(Pvalues)
        for p in Pvalues:
            X, Y = self.TP2XY(self.datawc_x1 + 273.15, p * 100)
            if not numpy.ma.isMA(X):
                dicP[Y] = dicPl[p]
        try:
            ttp = self.x.createtexttable('Plabels')
        except:
            ttp = self.x.gettexttable('Plabels')
        ttp.color = self.isobars.linecolors[0]
        isotemplate.ylabel1.texttable = ttp
        dicT = {}
        Tvalues = list(numpy.ma.array(Tvalues).filled().ravel())
        dicTl = vcs.mklabels(Tvalues)
        for t in Tvalues:
            X, Y = self.TP2XY(t + 273.15, self.datawc_y1 * 100)
            dicT[X] = dicTl[t]
        try:
            ttp = self.x.createtexttable('Tlabels')
        except:
            ttp = self.x.gettexttable('Tlabels')
        ttp.color = self.isotherms.linecolors[0]
        isotemplate.xlabel1.texttable = ttp
        isotemplate.ylabel1.priority = 1
        isotemplate.xlabel1.priority = 1
        isotemplate.data.x1 = isotemplate.data.x1
        isotemplate.data.x2 = isotemplate.data.x2
        isotemplate.data.y1 = isotemplate.data.y1
        isotemplate.data.y2 = isotemplate.data.y2
        self.isotherms.yticlabels1 = dicP
        self.isotherms.xticlabels1 = dicT
        self.isobars.yticlabels1 = {}
        self.isobars.xticlabels1 = {}
        self.dryadiabats.yticlabels1 = {}
        self.dryadiabats.xticlabels1 = {}
        self.pseudoadiabats.yticlabels1 = {}
        self.pseudoadiabats.xticlabels1 = {}
        self.mixingratio.yticlabels1 = {}
        self.mixingratio.xticlabels1 = {}
        self.isothermsfilled.datawc_x1 = self.isothermsfilled.datawc_x1
        self.isothermsfilled.datawc_x2 = self.isothermsfilled.datawc_x2
        self.isothermsfilled.datawc_y1 = self.isothermsfilled.datawc_y1
        self.isothermsfilled.datawc_x2 = self.isothermsfilled.datawc_y2
        # Puts the dims on it
        T.id = 'T'
        T.setAxis(0, Yaxis)
        T.setAxis(1, Xaxis)
        P = MV2.array(P)
        P.id = 'P'
        P.setAxis(0, Yaxis)
        P.setAxis(1, Xaxis)
        Tp = MV2.array(Tp)
        Tp.setAxis(0, Yaxis)
        Tp.setAxis(1, Xaxis)
        ws = MV2.array(ws)
        ws.setAxis(0, Yaxis)
        ws.setAxis(1, Xaxis)

        # plot if self.drawisothermsfilled:
        if self.drawisothermsfilled:
            self.displays.append(self.x.plot(T, isotemplate, self.isothermsfilled, bg=bg))
        if self.drawisotherms:
            self.displays.append(self.x.plot(T, isotemplate, self.isotherms, bg=bg))
        if self.drawisobars:
            self.displays.append(self.x.plot(P, isotemplate, self.isobars, bg=bg))
        if self.drawdryadiabats:
            self.displays.append(self.x.plot(Tp, isotemplate, self.dryadiabats, bg=bg))
        if self.drawpseudoadiabats:
            self.plot_pseudo(isotemplate, bg=bg)
        if self.drawmixingratio:
            self.displays.append(self.x.plot(ws, isotemplate, self.mixingratio, bg=bg))
        return

    # All take args with T in K and P in Pa
    # SkewT-LogP diagram defs
    def TP2XYskewT(self, T, P):
        X = T + self.skewness * numpy.ma.log(P)
        Y = -self.R * numpy.ma.log(P)
        return X, Y

    def XY2TPskewT(self, X, Y):
        P = numpy.ma.exp(-Y / self.R)
        T = X - self.skewness * numpy.ma.log(P)
        return T, P

    # Emagram defs
    def TP2XYemagram(self, T, P):
        X = T
        Y = -self.R * numpy.ma.log(P)
        return X, Y

    def XY2TPemagram(self, X, Y):
        T = X
        P = numpy.ma.exp(-Y / self.R)
        return T, P

    # Tephigram defs
    def TP2XYtephigram(self, T, P):
        X = T
        Tp = T / numpy.ma.power(P / 100000., self.k)
        Y = self.cp * numpy.ma.log(Tp)
        return X, Y

    def XY2TPtephigram(self, X, Y):
        T = X
        Tp = numpy.ma.exp(Y / self.cp)
        P = 100000. * numpy.ma.power(T / Tp, 1. / self.k)
        return T, P

    # Stuve diagram defs
    def TP2XYstuve(self, T, P):
        X = T
        Y = numpy.ma.power(P, self.k)
        return X, Y

    def XY2TPstuve(self, X, Y):
        T = X
        P = numpy.ma.power(Y, 1. / self.k)
        return T, P

    def plot_windbarb(self, u, v, P=None, template=None, type='uv', bg=0):
        if type == 'uv':
            n = numpy.ma.sqrt(u * u + v * v)
            d = numpy.ma.arccos(u / n)
            d = numpy.ma.where(numpy.ma.less(v, 0.), -d, d)

        if P is None:
            P = u.getLevel()
            if P is None:
                P = u.getAxis(-1)
            P = P.clone()
            try:  # tries to convert to Pa
                for i in range(len(P)):
                    uni = unidata.udunits(P[i], P.units)
                    P[i] = uni.to("Pa").value
            except Exception:
                pass
            P = P[:]

        n1 = n / self.windbarbsscales[0]
        n1 = n1.astype('i')
        n2 = (n - n1 * self.windbarbsscales[0]) / self.windbarbsscales[1]
        n2 = n2.astype('i')
        n3 = (n - n1 * self.windbarbsscales[0] - n2 * self.windbarbsscales[1]) / self.windbarbsscales[2]
        n3 = n3.astype('i')

# print 'n1:',n1
# print 'n2:',n2
# print 'n3:',n3

        nitems = len(P)
        for i in range(nitems):
            if not(MV2.isMA(n[i])):
                p = P[i]
                # Figure out the altitude to plot it (y coord) !
                dum, Y = self.TP2XY(273, p)
                lin = self.x.createline()
                lin.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
                lin.worldcoordinate = [-1, 1, -1, 1]
                lin.x = [0, 0]
                lin.y = [-1, 1]
                self.displays.append(self.x.plot(lin, bg=bg))
                lin = self.x.createline()
                lin.viewport = [template.data.x1, template.data.x2, template.data.y1, template.data.y2]
                lin.worldcoordinate = [-1, 1, self.ymin, self.ymax]
                # Set everything correctly y wise
                rw = 2 / (self.ymax - self.ymin)
                if self.x.islandscape():
                    r = 1.375
                else:
                    r = .72727272
                rv = (template.data.x2 - template.data.x1) / (template.data.y2 - template.data.y1)
                rw = 1. / rw
                x, y = self.make_barb(n[i], d[i], n1[i], n2[i], n3[i], rw * r * rv, Y)
                lin.x = x
                lin.y = y
                lin.type = ['solid']
                self.displays.append(self.x.plot(lin, bg=bg))

    def make_barb(self, n, d, n1, n2, n3, ydeformation, yoff):
        """ Draws a wind barb, the head is at location 0,0, the tail length is 1


        ---------------------------------------------------------------
        \  |  /   /    /
         \ | /   /    /
          \|/   /
           V

           total length is 1
           barbs angle is 60 degrees
           barb depth is .1

        """
        # Intitalize the Drawing, head is at zero
        xs = [0]
        ys = [0]
        xs.append(numpy.cos(d))
        ys.append(numpy.sin(d))

        # Ok now does the BIG barbs
        bangle = 25. / 180. * numpy.pi  # 30 degrees in radians
        bdepth = .3
        bwidth = numpy.tan(bangle) * bdepth

        hinit = 1.
        # Drawing of big barbs
        for i in range(n1):
            # summit
            l = hinit - bwidth  # legth to the middle of big barb
            h = numpy.sqrt(bdepth ** 2 + l ** 2)
            mu = numpy.arctan(bdepth / l)
            ys.append(h * numpy.sin(d + mu))
            xs.append(h * numpy.cos(d + mu))

            # End of barb
            ys.append((l - bwidth) * numpy.sin(d))
            xs.append((l - bwidth) * numpy.cos(d))
            hinit = hinit - 2 * bwidth
        if n1 != 0:
            hinit = hinit + bwidth / 2.

        # Ok now the long sticks
        for i in range(n2):
            l = hinit - bwidth
            # Move to begining of next one
            ys.append(l * numpy.sin(d))
            xs.append(l * numpy.cos(d))

            # Go to the top one
            l = hinit
            h = numpy.sqrt(bdepth ** 2 + l ** 2)
            mu = numpy.arctan(bdepth / l)
            ys.append(h * numpy.sin(d + mu))
            xs.append(h * numpy.cos(d + mu))

            # Back to begining
            ys.append(ys[-2])
            xs.append(xs[-2])

            # Adjust hinit
            hinit = hinit - bwidth / 2.

        # Now the half barbs
        bwidth = bwidth / 2.
        bdepth = bdepth / 2.
        for i in range(n3):
            l = hinit - bwidth
            # Move to begining of next one
            ys.append(l * numpy.sin(d))
            xs.append(l * numpy.cos(d))

            # Go to the top one
            l = hinit
            h = numpy.sqrt(bdepth ** 2 + l ** 2)
            mu = numpy.arctan(bdepth / l)
            ys.append(h * numpy.sin(d + mu))
            xs.append(h * numpy.cos(d + mu))
            # Back to begining
            ys.append(ys[-2])
            xs.append(xs[-2])

            # Adjust hinit
            hinit = hinit - bwidth / 2.

        # Scales and offset
        for i in range(len(ys)):
            # print 'i,x,y:',i,xs[i],ys[i]
            ys[i] *= ydeformation
            ys[i] += yoff

        return xs, ys

# Adapted for numpy/ma/cdms2 by convertcdms.py Computes saturated pressure in Pa given T in K, using the method:     1: Hyland-Wexler formulation, polynomial coeff (absolute norm)     2: Wexler formulation     3: Hyland-Wexler formulation, polynomial coeff (relative error norm)     4: classic Goff Gratch equation     5: 6.112*numpy.ma.exp(17.67*tempc/(tempc+243.5))      Default is method 1      Note: 1 and 2 use method 3 where T is not : 173.15 < T < 473.15     ref for 1, 2 and 3:     <NAME>., Journal of Applied Met., Vol 31, Dec 1992     ( http://ams.allenpress.com/perlserv/?request=get-document&\ doi=10.1175%2F1520-0450(1992)031%3C1507%3APFTSVP%3E2.0.CO%3B2&ct=1 ) # Put in C # Water vapor # ice # Combine # Overwrite values outside valid range with method 2 # over water # over ice # over ice # over water # Put in C # Water vapor # ice # Combine # Overwrite values outside valid range with method 2 # Pa computes pressure values from level axis,     if not designated level then use last axis... # tries to convert to Pa Compute saturated mixing ratio at T,P, cpomputing es using the method:     T in K     P in Pa Function to calculate the moist adiabatic lapse rate (K/Pa) based     on the temperature, pressure, and rh of the environment.     T in K     rh in %     P in Pa     if P is not passed will be optaines from T level axis (or last axis) # in SI Lift a parcel moist adiabatically from Po to Pend.     Initial temperature is To in K.     To in K     Po in Pa     detail # number to increment P slowly     returns T at all P (every deltaP Pa), and P # MUCH faster Computes dewpoint profile from T and rh     Usage:       Td = Dewpoint(T,rh,method=None)     Where:      T is in K      rh is in % (relative humidity)      method used to compute Es:          1: Hyland-Wexler formulation, polynomial coeff (absolute norm) (default)          2: Wexler formulation          3: Hyland-Wexler formulation, polynomial coeff (relative error norm)          4: classic Goff Gratch equation          5: 6.112*numpy.ma.exp(17.67*tempc/(tempc+243.5))        Assumes constant latent temperature (except over ice or water) # K # J/(K*kg) # Pa # J/kg # J/kg # in Pa # Figures out which kind of diagram we're drawing # user defined diagram # Gets the xs for each corner # if name in thermo_objects: #             raise Exception,"Error thermo object %s already exists" % name # thermo_objects.append(name) # Set some constants # self.L=2.5E3 # World coordiantes in C/hPa # Skewness (for skewT) # Sets the diagram type # Level of detail # Draw mixing ratio and pseudoadiabats up to # Do we draw all the lines on the diagram ? (default, yes) # Creates the isolines # Isofill stuff for temperatures #         self.isothermsfilled.levels=([-240.0, -200.0], [-160.0, -120.0], [-80.0, -40.0], [0.0, 40.0], [80.0, 120.0]) # Isotherms # Isobars default settings # Dry adiabatics # self._dryadiabats.level=range(-40,170,10) # Pseudoadiabatics # self._pseudoadiabats.level=vcs.mkevenlevels(-40,40,40) # Mixing ratio # Windbarb # Lines for plot # print 'windbarbs =',self.windbarbs.name # 'TP2XYskewT', # 'TP2XYemagram', # 'TP2XYtephigram', # 'TP2XYstuve', ## Ok this is pathetic but i can't get it to work any other way!!!         ## We'll have to do the integration for each levels # Ok we now have the required levels # if not(MV2.isMA(tx) or MV2.isMA(ty)): # xs.append(tx) # ys.append(ty) # else: # print tx,ty,t,p # xs.append(tx) # ys.append(ty) # Computes T,P on this grid # Potential Temperature # Seems like we need not ws after 600mb # Now sets the isothermsfilled # self.setIso(self.isothermsfilled,Tvalues2) # Pseudoadiabats # Mixing Ratio # Figures out where to plot the P labels # Puts the dims on it # plot if self.drawisothermsfilled: # All take args with T in K and P in Pa # SkewT-LogP diagram defs # Emagram defs # Tephigram defs # Stuve diagram defs # tries to convert to Pa # print 'n1:',n1 # print 'n2:',n2 # print 'n3:',n3 # Figure out the altitude to plot it (y coord) ! # Set everything correctly y wise Draws a wind barb, the head is at location 0,0, the tail length is 1           ---------------------------------------------------------------         \  |  /   /    /          \ | /   /    /           \|/   /            V             total length is 1            barbs angle is 60 degrees            barb depth is .1 # Intitalize the Drawing, head is at zero # Ok now does the BIG barbs # 30 degrees in radians # Drawing of big barbs # summit # legth to the middle of big barb # End of barb # Ok now the long sticks # Move to begining of next one # Go to the top one # Back to begining # Adjust hinit # Now the half barbs # Move to begining of next one # Go to the top one # Back to begining # Adjust hinit # Scales and offset # print 'i,x,y:',i,xs[i],ys[i]