/beam-3dveglab-vlab/src/main/resources/auxdata/VLabImpl.py

#
# Copyright (C) 2010-2014 Netcetera Switzerland (info@netcetera.com)
# 
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 3 of the License, or (at your option)
# any later version.
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
# more details.
# 
# You should have received a copy of the GNU General Public License along
# with this program; if not, see http://www.gnu.org/licenses/
# 
# @(#) $Id: $
#
# Authors: Cyrill Schenkel, Daniel Kueckenbrink, Joshy Cyriac, Marcel Kessler, Jason Brazile
#

##############################################################################
# Three ways to run it:
# 
# 1. embedded within the BEAM processor (normal case)
#
# 2. standalone tests - headless (either jython or python)
#
#    jython -Dpython.path=jcommon-1.0.16.jar:jfreechart-1.0.13.jar VLabImpl.py
#    python VLabImpl.py
#
# 3. standalone with a 'fake' swing-based gui
# 
#    java -jar ${HOME}/beam-4.11/lib/jython-2.5.2.jar -Dvlab.fakebeam=1 -Dpython.path=${HOME}/beam-4.11/lib/jcommon-1.0.16.jar:${HOME}/beam-4.11/lib/jfreechart-1.0.13.jar VLabImpl.py
#
# Those BEAM-supplied jars can also be obtained like this:
#    wget -U "Mozilla/5.0" http://repo1.maven.org/maven2/jfree/jfreechart/1.0.13/jfreechart-1.0.13.jar
#    wget -U "Mozilla/5.0" http://repo1.maven.org/maven2/jfree/jcommon/1.0.16/jcommon-1.0.16.jar
#
#

import sys, math, operator, array, time, struct
from array import array

class VLAB:
  """VLAB contains conf. constants, static utility methods, and test methods"""

  COPYRIGHT_INFO     = 'Copyright (C) 2010-2014 Netcetera Switzerland (info@netcetera.com)'
  PROCESSOR_NAME     = 'BEAM VLab Processor'
  PROCESSOR_SNAME    = 'beam-vlab'
  REQUEST_TYPE       = 'VLAB'
  UI_TITLE           = 'VLab - Processor'
  VERSION_STRING     = '1.0 (20 July 2014)'
  DEFAULT_LOG_PREFIX = 'vlab'
  LOGGER_NAME        = 'beam.processor.vlab'

  D_PRODNAME         = 'vlab_out.dim'
  P_CONDITION        = '.condition'
  P_EXPRESSION       = '.expression'
  P_OUTPUT           = '.output'

  # NOTE: Once released, random number generation should NOT be reproducible
  CONF_RND_REPRODUCE = True

  JCB                = 'JComboBox'
  JTF                = 'JTextField'
  JTD                = (JTF, False)

  K_LIBRAT           = 'librat'
  K_DART             = 'DART'
  K_DUMMY            = 'dummy'

  K_LAEGERN          = 'Laegern'
  K_THARANDT         = 'Tharandt'
  K_RAMI             = 'RAMI'
  K_USER_DEFINED     = 'User defined (UserDefined.obj)'

  K_YES              = 'Yes'
  K_NO               = 'No'

  K_SENTINEL2        = 'MSI (Sentinel 2)'
  K_SENTINEL3        = 'OLCI (Sentinel 3)'
  K_MODIS            = 'MODIS'
  K_MERIS            = 'MERIS'
  K_LANDSAT_OLI      = 'Landsat (OLI)'
  K_LANDSAT_ETM      = 'Landsat (ETM)'

  K_RURAL            = 'Rural'
  K_MARITIME         = 'Maritime'
  K_URBAN            = 'Urban'
  K_TROPOSPHERIC     = 'Tropospheric'

  DBG                = 'debug' 
  INF                = 'info'
  ERR                = 'error'

  plst               = []

  if sys.platform.startswith('java'):
    from java.util.logging import Logger
    from java.util.logging import FileHandler
    from java.util.logging import Level
    logger = Logger.getLogger(LOGGER_NAME)
    # comment these out
    # logfh = FileHandler('%s.log' % LOGGER_NAME)
    # logfh.setLevel(Level.FINEST)
    # logger.addHandler(logfh)
  else:
    import logging
    logger = logging.getLogger(LOGGER_NAME)
    logger.setLevel(logging.DEBUG)
    logch = logging.StreamHandler()
    logch.setLevel(logging.DEBUG)
    logger.addHandler(logch)
    # comment these out
    # logfh = logging.FileHandler('%s.log' % LOGGER_NAME)
    # logfh.setLevel(logging.DEBUG)
    # logger.addHander(logfh)

  model = (
{'Forward Modeling': (
{'Model Selection': (
 ('3D Scene',          '3dScene',             JCB, (K_RAMI, K_LAEGERN, K_THARANDT, K_USER_DEFINED)),
 ('RT Processor',      'RTProcessor',         JCB, (K_DUMMY, K_LIBRAT, K_DART)))},
{'Spectral Characteristics': (
 ('Sensor',            'Sensor',              JCB, (K_SENTINEL2, K_SENTINEL3, K_MODIS, K_MERIS, K_LANDSAT_OLI, K_LANDSAT_ETM)),
 ('Bands',             'Bands',               JTF, '1, 2, 3, 4, 5, 6, 7, 8, 9, 10'))},
{'Viewing Characteristics': (
 ('Zenith',            'ViewingZenith',       JTF, '20.0'),
 ('Azimuth',           'ViewingAzimuth',      JTF, '0.0'))},
{'Illumination Characteristics':(
 ('Zenith',            'IlluminationZenith',  JTF, '20.0'),
 ('Azimuth',           'IlluminationAzimuth', JTF, '0.0'))},
{'Scene Parameters': (
 ('Pixel Size',        'ScenePixel',          JTD, 'not available yet'),
 ('(Alt A) Filename',  'SceneLocFile',        JTD, 'not available yet'),
 ('(Alt B) XC',        'SceneXC',             JTD, 'not available yet'),
 ('(Alt B) XW',        'SceneXW',             JTD, 'not available yet'),
 ('(Alt B) YC',        'SceneYC',             JTD, 'not available yet'),
 ('(Alt B) YW',        'SceneYW',             JTD, 'not available yet'))},
{'Atmospheric Parameters': (
 ('CO2 Mixing Ratio',  'AtmosphereCO2',       JTF, '1.6'),
 ('Aerosol Profile',   'AtmosphereAerosol',   JCB, (K_RURAL, K_MARITIME, K_URBAN, K_TROPOSPHERIC)),
 ('Water Vapor',       'AtmosphereWater',     JTF, '0.0'),
 ('Ozone Column',      'AtmosphereOzone',     JTF, '300'))},
{'Output Parameters': (
 ('Result file prefix','OutputPrefix',        JTF, 'HET01_DIS_UNI_NIR_20.obj'),
 ('Result Directory',  'OutputDirectory',     JTF, 'dart.rpv.rami.2'),
 ('Image file',        'ImageFile',           JCB, (K_YES, K_NO)),
 ('Ascii file',        'AsciiFile',           JCB, (K_YES, K_NO)))}
)},
{'DHP Simulation': (
{'Model Selection': (
 ('3D Scene',          'DHP_3dScene',         JCB, (K_RAMI, K_LAEGERN, K_THARANDT)),
 (),
 ('Resolution',        'DHP_Resolution',      JTF, '4000000'),
 ())},
{'DHP Location': (
 ('X',                 'DHP_X',               JTD, 'not available yet'),
 ('Y',                 'DHP_Y',               JTD, 'not available yet'))},
{'DHP Properties': (
 ('Zenith',            'DHP_Zenith',          JTF, '0.0'),
 ('Azimuth',           'DHP_Azimuth',         JTF, '0.0'))},
{'DHP Imaging Plane': (
 ('Orientation',       'DHP_Orientation',     JTD, 'not available yet'),
 ('Height(z)',         'DHP_Height',          JTD, 'not available yet'))},
{'Output Parameters': (
 ('Result file prefix','DHP_OutputPrefix',    JTF, 'RAMI00_'),
 ('Result Directory',  'DHP_OutputDirectory', JTF, ''),
 ('Image file',        'DHP_ImageFile',       JCB, (K_YES, K_NO)),
 ('Ascii file',        'DHP_AsciiFile',       JCB, (K_YES, K_NO)))}
)},
)
  # set parameter names
  for tabGroups in model:
    for tabName in tabGroups:
      for group in tabGroups[tabName]:
        for groupName in group:
          for groupTuple in group[groupName]:
            if len(groupTuple) == 4:
              (lbl, nm, typ, vals) = groupTuple
              exec('P_' + nm + ' = nm')
              exec('L_' + nm + ' = lbl')
              plst.append(nm)

  def __init__(self):
    self.cmap = {}
    self.vmap = {}

  def me():
    """Returns name of currently executing method"""
    nm = ''
    try:
      raise ZeroDivisionError
    except ZeroDivisionError:
      nm = sys.exc_info()[2].tb_frame.f_back.f_code.co_name
    return nm+'()'
  me = staticmethod(me)
  def lineSeparator():
    """Return the OS line separator"""
    if sys.platform.startswith('java'):
      from java.lang import System
      if sys.platform.startswith('java1.6'):
        return System.getProperty('line.separator')
      elif sys.platform.startswith('java1.7'):
        return System.lineSeparator()
    else:
      import os
      os.linesep
  lineSeparator = staticmethod(lineSeparator)
  def listdir(path):
    """list files in the directory given by path"""
    if sys.platform.startswith('java'):
      from java.io import File
      array = File(path).list()
      listFile = []
      if array != None:
        for i in xrange(len(array)):
          listFile.append(array[i])
      return listFile
    else:
      import os
      return os.listdir(path)
  listdir = staticmethod(listdir)
  def getenv(key, default=None):
    if sys.platform.startswith('java'):
      from java.lang.System import getenv
      return getenv(key)
    else:
      import os
      return os.getenv(key)
  getenv = staticmethod(getenv)
  def checkFile(fname):
    """open a file if it exists, otherwise die"""
    try:
      fp = open(fname, 'r+')
      return fp
    except IOError, e:
      raise RuntimeError(e)
  checkFile = staticmethod(checkFile)
  def fileExists(fname):
    """check if fname exists as a file"""
    if sys.platform.startswith('java'):
      from java.io import File
      return File(fname).exists()
    else:
      import os
      return os.path.exists(fname)
  fileExists = staticmethod(fileExists)
  def getFullPath(fname):
    """return canonical/absolute path of given fname"""
    if sys.platform.startswith('java'):
      from java.io import File
      return File(fname).getCanonicalPath()
    else:
      import os
      return os.path.abspath(fname)
  getFullPath = staticmethod(getFullPath)
  def copyDir(dname, target):
    """Recursively copy 'dname' directory to 'target'.
       /!\If 'target' already exists it will be removed or overwrited.
    """
    if sys.platform.startswith('java'):
      # import java module
      from java.io import File
      from java.io import FileInputStream
      from java.io import FileOutputStream
      from java.util import Scanner
      from java.lang import String
      dnameFile = File(dname)
      targetFile = File(target)
      # recursive copy
      if dnameFile.isDirectory():
        # Create folder if not exists
        if not targetFile.exists():
          targetFile.mkdir()
        # Copy all content recursively
        for fname in dnameFile.list().tolist():
          VLAB.copyDir(dname + File.separator + fname, target + File.separator + fname)
      else:
        # Read dname file
        istream = FileInputStream(dname)
        scan = Scanner(istream).useDelimiter("\\Z")
        # Test if file is empty
        if scan.hasNextLine():
          content = String(scan.next())
        else:
          content = String("")
        scan.close()
        # Create and write target
        if not targetFile.exists():
          targetFile.createNewFile()
        ostream = FileOutputStream(target)
        ostream.write(content.getBytes())
        ostream.flush()
        ostream.close()
    else:
      import shutil, os
      # remove exisiting target
      if os.path.isdir(target) or os.path.isfile(target):
        shutil.rmtree(target)
      # recursive copy of dnma to target
      shutil.copytree(dname, target)
  copyDir = staticmethod(copyDir)
  def XMLEditNode(fname, nodeName, attributName, value, multiple=False):
    """ Edit a given node (nodeName) in a given XML file (fname)
    attributName and value could be either a list of string or a string
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from javax.xml.transform import TransformerFactory
      from javax.xml.transform import OutputKeys
      from javax.xml.transform.dom import DOMSource
      from javax.xml.transform.stream import StreamResult
      from java.io import File
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(nodeName)
      # Check if we get only one node (as expected)
      if node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      elif node.getLength() == 1:
        nodes = [node.item(0)]
      else:
        if not multiple:
          raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
        else:
          nodes = [ node.item(i) for i in xrange(node.getLength()) ]
      for node in nodes:
        # Modify the node attribute
        if isinstance(attributName, list) and isinstance(value, list):
          for att, val in zip(attributName, value):
            node.setAttribute(att, val)
        elif isinstance(attributName, str) and isinstance(value, str):
          node.setAttribute(attributName, value)
        else:
          raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      transformer = TransformerFactory.newInstance().newTransformer()
      transformer.setOutputProperty(OutputKeys.INDENT, "yes")
      source = DOMSource(tree)
      result = StreamResult(File(fname))
      transformer.transform(source, result)
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      nodes = tree.findall(".//*../%s" % nodeName)
      # Check if we get only one node (as expected)
      if len(nodes) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      elif len(nodes) > 1 and not multiple:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
      for node in nodes:
        # Modify the node attribute
        if isinstance(attributName, list) and isinstance(value, list):
          for att, val in zip(attributName, value):
            node.set(att, val)
        elif isinstance(attributName, str) and isinstance(value, str):
          node.set(attributName, value)
        else:
          raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      tree.write(fname)
  XMLEditNode = staticmethod(XMLEditNode)
  def XMLReplaceNodeContent(fname, parent, subnode, attributName, value, spectralBands=False):
    """ Edit an XML file (fname) and replace the content of a node with subnode(s) (subnode) within attribute(s) and value(s).
    attributName and value could be either a list of string or a string
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from javax.xml.transform import TransformerFactory
      from javax.xml.transform import OutputKeys
      from javax.xml.transform.dom import DOMSource
      from javax.xml.transform.stream import StreamResult
      from java.io import File
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(parent)
      # Check if we get only one node (as expected)
      if node.getLength() > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node.item(0)
      # Remove content node
      while node.hasChildNodes():
        node.removeChild(node.getFirstChild())
      # Modify the node attribute
      elem = tree.createElement(subnode)
      if spectralBands:
        elem.setAttribute("bandNumber", "0")
        elem.setAttribute("spectralDartMode", "0")
      if isinstance(attributName, list) and isinstance(value, list):
        if isinstance(value[0], list):
          for bandNumber, val in enumerate(value):
            if spectralBands:
              elem.setAttribute("bandNumber", str(bandNumber))
            for atr, v in zip(attributName, val):
              elem.setAttribute(atr, v)
            node.appendChild(elem.cloneNode(True))
        else:
          elem = tree.createElement(subnode)
          for atr, v in zip(attributName, value):
            elem.setAttribute(atr, v)
          node.appendChild(elem)
      elif isinstance(attributName, str) and isinstance(value, str):
        elem = tree.createElement(subnode)
        elem.setAttribute(attributName, value)
        node.appendChild(elem)
      else:
        raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      transformer = TransformerFactory.newInstance().newTransformer()
      transformer.setOutputProperty(OutputKeys.INDENT, "yes")
      source = DOMSource(tree)
      result = StreamResult(File(fname))
      transformer.transform(source, result)
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      node = tree.findall(".//*../%s" % parent)
      # Check if we get only one node (as expected)
      if len(node) > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif len(node) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node[0]
      # Remove content of node
      node.clear()
      # Modify the node attribute
      if spectralBands:
        attrib = {"bandNumber":"0", "spectralDartMode":"0"}
      else:
        attrib = {}
      if isinstance(attributName, list) and isinstance(value, list):
        if isinstance(value[0], list):
          for bandNumber, val in enumerate(value):
            if spectralBands:
              attrib["bandNumber"] = str(bandNumber)
            for atr, v in zip(attributName, val):
              attrib[atr] = v
            node.append(ET.Element(subnode, attrib=attrib))
        else:
          for atr, v in zip(attributName, value):
            attrib[atr] = v
          node.append(ET.Element(subnode, attrib=attrib))
      elif isinstance(attributName, str) and isinstance(value, str):
        attrib[attributName] = value
        node.append(ET.Element(subnode, attrib=attrib))
      else:
        raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      tree.write(fname)
  XMLReplaceNodeContent = staticmethod(XMLReplaceNodeContent)
  def XMLAddNode(fname, parent, treeNodes, nodesSetup):
    """ Add a node (subnode) to a given parent in the provided fname file.
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from javax.xml.transform import TransformerFactory
      from javax.xml.transform import OutputKeys
      from javax.xml.transform.dom import DOMSource
      from javax.xml.transform.stream import StreamResult
      from java.io import File
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(parent)
      # Check if we get only one node (as expected)
      if node.getLength() > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node.item(0)
      # Create all the new nodes
      nodes = {parent:node}
      for name in treeNodes:
        nodes[name] = (tree.createElement(name))
        for atr, val in zip(*[ nodesSetup[name][key] for key in nodesSetup[name].iterkeys() if key in ["attribute", "value"] ]):
          nodes[name].setAttribute(atr, val)
      # Add all created node to its parent (previous node in the list)
      for name, elem in nodes.iteritems():
        if name != parent:
          nodes[nodesSetup[name]["parent"]].appendChild(elem)
      # Write new XML tree in fname
      transformer = TransformerFactory.newInstance().newTransformer()
      transformer.setOutputProperty(OutputKeys.INDENT, "yes")
      source = DOMSource(tree)
      result = StreamResult(File(fname))
      transformer.transform(source, result)
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      node = tree.findall(".//*../%s" % parent)
      # Check if we get only one node (as expected)
      if len(node) > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif len(node) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node[0]
      # Create all the new nodes
      nodes = {parent:node}
      for name in treeNodes:
        attrib = {}
        for atr, val in zip(*[ nodesSetup[name][key] for key in nodesSetup[name].iterkeys() if key in ["attribute", "value"] ]):
          attrib[atr] = val
        nodes[name] = ET.Element(name, attrib)
      # Add all created node to its parent (previous node in the list)
      for name, elem in nodes.iteritems():
        if name != parent:
          nodes[nodesSetup[name]["parent"]].append(elem)
      # Write new XML tree in fname
      tree.write(fname)
  XMLAddNode = staticmethod(XMLAddNode)
  def XMLGetNodeAttributeValue(fname, nodeName, attributName):
    """ Return the value of the given attributName for the given nodeName
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from org.w3c.dom import Element
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(nodeName)
      # Check if we get only one node (as expected)
      if node.getLength() > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
      elif node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      else:
        node = node.item(0)
      # Check if that node has attributName
      if node.hasAttribute(attributName):
        return node.getAttributes().getNamedItem(attributName).getNodeValue()
      else:
        raise IOError("Attribute name '%s' not found for node '%s' in file '%s'" % (attributName, nodeName, fname))
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      node = tree.findall(".//*../%s[@%s]" % (nodeName, attributName))
      # Check if we get only one node (as expected)
      if len(node) > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
      elif len(node) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      else:
        node = node[0]
      # Check if that node has an attributName
      if attributName in node.keys():
        return node.get(attributName)
      else:
        raise IOError("Attribute name '%s' not found for node '%s' in file '%s'" % (attributName, nodeName, fname))
  XMLGetNodeAttributeValue = staticmethod(XMLGetNodeAttributeValue)
  def getBandsFromGUI(bands):
    """ Return a DART spectral bands list: ["deltaLambda", "meanLambda"] in micro meter
    In case of several bands the result should be a list of list:
    [["deltaLambda0", "meanLambda0"], ["deltaLambda1", "meanLambda1"], ["deltaLambda2", "meanLambda2"], ...]
    e.g.: [["0.02", "0.56"], ["0.02", "0.58"], ["0.02", "0.60"], ["0.02", "0.62"]]
    """
    # TODO: You should write the spectral band converter here!
    return [["0.02", "0.56"], ["0.02", "0.58"], ["0.02", "0.60"], ["0.02", "0.62"]]
  getBandsFromGUI = staticmethod(getBandsFromGUI)
  class path:
    def exists(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).exists()
      else:
        import os
        return os.path.exists(path)
    exists = staticmethod(exists)
    def isabs(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).isAbsolute()
      else:
        import os
        return os.path.isabs(path)
    isabs = staticmethod(isabs)
    def isdir(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).isDirectory()
      else:
        import os
        return os.path.isdir(path)
    isdir = staticmethod(isdir)
    def join(path, *args):
      if sys.platform.startswith('java'):
        from java.io import File
        f = File(path)
        for a in args:
          g = File(a)
          if g.isAbsolute() or len(f.getPath()) == 0:
            f = g
          else:
            f = File(f, a)
        return f.getPath()
      else:
        import os
        return os.path.join(path, *args)
    join = staticmethod(join)
    def isfile(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).isfile()
      else:
        import os
        return os.path.isfile(path)
    isfile = staticmethod(isfile)
    def normcase(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).getPath()
      else:
        import os
        return os.path.normcase(path)
    normcase = staticmethod(normcase)
    def normpath(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).getCanonicalPath()
      else:
        import os
        return os.path.normpath(path)
    normpath = staticmethod(normpath)
    def split(path):
      if sys.platform.startswith('java'):
        from java.io import File
        f=File(path)
        d=f.getParent()
        if not d:
          if f.isAbsolute():
            d = path
          else:
            d = ""
        return (d, f.getName())
      else:
        import os
        return os.path.split(path)
    split = staticmethod(split)
  def frange(end, start=0, inc=0):
    """a jython-compatible range function with float increments"""
    import math
    if not start:
      start = end + 0.0
      end = 0.0
    else: end += 0.0
    if not inc:
      inc = 1.0
    count = int(math.ceil((start - end) / inc))
    L = [None] * count
    L[0] = end
    for i in (xrange(1,count)):
      L[i] = L[i-1] + inc
    return L
  frange = staticmethod(frange)
  def openFileIfNotExists(filename):
    """open file exclusively"""
    if sys.platform.startswith('java'):
      from java.io import File
      File(filename).getParentFile().mkdirs()
      if File(filename).createNewFile():
        return open(filename, 'w')
      else:
        return None
    else:
      import os
      try:
        # os.O_EXCL => open exclusive => acquire a lock on the file
        fd = os.open(filename, os.O_CREAT|os.O_EXCL|os.O_WRONLY|os.O_TRUNC)
      except:
        return None
      fobj = os.fdopen(fd,'w')
      return fobj
  openFileIfNotExists = staticmethod(openFileIfNotExists)
  def rndInit(seed = None):
    """initialize the randon number generator"""
    if sys.platform.startswith('java'):
      from java.util import Random
      if seed == None:
        return Random()
      else:
        return Random(seed)
    else:
      import random
      random.seed(seed)
      return None
  rndInit = staticmethod(rndInit)
  def rndNextFloat(randState):
    """return the next pseudo-random floating point number in the sequence"""
    if sys.platform.startswith('java'):
      from java.util import Random
      return randState.nextFloat()
    else:
      import random
      return random.random()
  rndNextFloat = staticmethod(rndNextFloat)
  def r2d(v):
    """jython-compatible conversion of radians to degrees"""
    if sys.platform.startswith('java'):
      from java.lang import Math
      return Math.toDegrees(v)
    else:
      return math.degrees(v)
  r2d = staticmethod(r2d)
  def d2r(v):
    """jython-compatible conversion of degrees to radians"""
    if sys.platform.startswith('java'):
      from java.lang import Math
      return Math.toRadians(float(v))
    else:
      return math.radians(v)
  d2r = staticmethod(d2r)
  def mkDirPath(path):
    """create directory (including non-existing parents) for the given path"""
    if sys.platform.startswith('java'):
      from java.io import File
      if not File(path).isDirectory():
        if not File(path).mkdirs():
          raise RuntimeError('failed to mkdir %s' % path)
    else:
      import os
      try:
        os.stat(path)
      except:
        os.makedirs(path)
  mkDirPath = staticmethod(mkDirPath)
  def fPath(d,n):
    """get file path of the file defined by directory d and file name n"""
    if sys.platform.startswith('java'):
      from java.io import File
      return File(d, n).getPath()
    else:
      import os
      return os.path.join(d, n)
  fPath = staticmethod(fPath)
  def savetxt(a, b, fmt=False):
    """save text b in a text file named a"""
    fh = open(a, 'w')
    if not fmt:
      fmt = '%s'
    for row in b:
      for element in row:
        fh.write(fmt % element + ' ')
      fh.write('\n')
    fh.close()
  savetxt = staticmethod(savetxt)
  def osName():
    """return which OS we are currently running on"""
    if sys.platform.startswith('java'):
      from java.lang import System
      oname = System.getProperty('os.name')
    else:
      import os
      oname = os.name
      if not oname.endswith('x'): oname = 'Windows'
    return oname
  osName = staticmethod(osName)
  def expandEnv(instr):
    """potentially expand environment variables in the given string"""
    outstr = instr
    m = {'$HOME':'HOME','%HOMEDRIVE%':'HOMEDRIVE','%HOMEPATH%':'HOMEPATH'}
    for e in m:
      if outstr.find(e) != -1:
        if sys.platform.startswith('java'):
          from java.lang import System
          repl = System.getenv(m[e])
        else:
          import os
          repl = os.getenv(m[e])
        if repl != None:
          outstr = outstr.replace(e, repl)
    return outstr
  expandEnv = staticmethod(expandEnv)

  if sys.platform.startswith('java'):
    from java.lang import Runnable
    class Helper(Runnable):
      def __init__(self, nm, strm, fName):
        self.nm=nm; self.strm=strm; self.fp=None
        if fName != None:
          if VLAB.fileExists(fName):
            self.fp = open(fName, 'w')
          else:
            self.fp = VLAB.openFileIfNotExists(fName)
      def run(self):
        """helper class for slurping up child streams"""
        from java.io import BufferedReader
        from java.io import InputStreamReader
        line = None; br = BufferedReader(InputStreamReader(self.strm))
        line = br.readLine()
        while (line != None):
          if self.fp != None:
            self.fp.write(line + VLAB.lineSeparator())
            self.fp.flush()
          VLAB.logger.info('%s %s' %(self.nm, line.rstrip()))
          line = br.readLine()
        br.close()
        if self.fp != None:
          self.fp.close()
  else:
    def helper(nm, strm):
      """helper method for slurping up child streams"""
      for line in strm: VLAB.logger.info('%s %s' %(nm, line.rstrip()))
      if not strm.closed: strm.close()
    helper = staticmethod(helper)
  def doExec(cmdrec):
    """run the specified external program under windows or unix"""
    cmdLine = []
    osName = VLAB.osName()
    if osName.startswith('Windows'):
      cmd=cmdrec['windows']
      cmdLine = ['cmd', '/c']
    else:
      cmd=cmdrec['linux']
    exe = VLAB.expandEnv(cmd['exe'])
    if not (VLAB.fileExists(exe) or osName.startswith('Windows')):
      raise RuntimeError('Cannot find exe "%s"' % exe)
    cmdLine.append(exe)
    for i in cmd['cmdline']:
      cmdLine.append(VLAB.expandEnv(i))

    VLAB.logger.info('cmdLine is [%s]' % cmdLine)
    if sys.platform.startswith('java'):
      from java.lang import ProcessBuilder
      from java.lang import Thread
      from java.io import BufferedWriter
      from java.io import OutputStreamWriter
      from java.io import File
      pb = ProcessBuilder(cmdLine)
      if 'cwd' in cmd and cmd['cwd'] != None:
        pb.directory(File(VLAB.expandEnv(cmd['cwd'])))
      if 'env' in cmd and cmd['env'] != None:
        env = pb.environment()
        cmdenv = cmd['env']
        for e in cmdenv:
          env[e] = VLAB.expandEnv(cmdenv[e])
      proc = pb.start()
      stdoutfName = None
      if 'stdout' in cmd and cmd['stdout'] != None:
        stdoutfName = VLAB.expandEnv(cmd['stdout'])
      stderrfName = None
      if 'stderr' in cmd and cmd['stderr'] != None:
        stderrfName = VLAB.expandEnv(cmd['stderr'])
      outs = Thread(VLAB.Helper('out', proc.getInputStream(), stdoutfName))
      errs = Thread(VLAB.Helper('err', proc.getErrorStream(), stderrfName))
      outs.start(); errs.start()
      bw = BufferedWriter(OutputStreamWriter(proc.getOutputStream()))
      if 'stdin' in cmd and cmd['stdin'] != None:
        inFile = VLAB.expandEnv(cmd['stdin'])
        VLAB.logger.info('stdin is [%s]' % inFile)
        if 'cwd' in cmd and cmd['cwd'] != None:
          if not VLAB.fileExists(inFile):
            # try pre-pending the cwd
            inFile = VLAB.fPath(VLAB.expandEnv(cmd['cwd']), inFile)
            if not VLAB.fileExists(inFile):
              raise RuntimeError('Cannot find stdin "%s"' % cmd['cwd'])
        fp = open(inFile, 'r')
        for line in fp:
          bw.write(line)
        bw.close()
        fp.close()
      exitCode = proc.waitFor()
      outs.join(); errs.join()
    else:
      import threading, subprocess, os
      if 'cwd' in cmd and cmd['cwd'] != None:
        os.chdir(VLAB.expandEnv(cmd['cwd']))
      if 'env' in cmd and cmd['env'] != None:
        cmdenv = cmd['env']
        for e in cmdenv:
          os.putenv(e, VLAB.expandEnv(cmdenv[e]))
      if 'stdin' in cmd and cmd['stdin'] != None:
        proc = subprocess.Popen(cmdLine, stdin=open(VLAB.expandEnv(cmd['stdin']),'r'),stdout=subprocess.PIPE, stderr=subprocess.PIPE)
      else:
        proc = subprocess.Popen(cmdLine, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
      t1 = threading.Thread(target=VLAB.helper, args=('out', proc.stdout))
      t2 = threading.Thread(target=VLAB.helper, args=('err', proc.stderr))
      t1.start(); t2.start()
      exitCode = proc.wait()
      t1.join(); t2.join()
    VLAB.logger.info('exitCode=%d' % exitCode)
  doExec = staticmethod(doExec)

  def valuesfromfile(path, transpose=False):
    """Returns a 2D array with the values of the csv file at 'path'.

    Keyword arguments:
    transpose -- transpose the matrix with the values

    """
    fp = open(path)
    values = [line.strip().split() for line in fp
              if not line.startswith('#')]
    fp.close()
    values = [[float(value) for value in row] for row in values]
    if transpose:
      values = zip(*values)
    return values
  valuesfromfile = staticmethod(valuesfromfile)
  def fabsa(x):
    """Return the element-wise abs() of the given array."""
    return map(lambda x : math.fabs(x), x)
  fabsa = staticmethod(fabsa)
  def cosa(arcs):
    """Return the element-wise cos() of 'arcs' array given in degrees."""
    return map(lambda x : math.cos(VLAB.d2r(x)), arcs)
  cosa = staticmethod(cosa)
  def replacerectinarray(array, replacement, xul, yul, xlr, ylr):
    """Replace the array with the specified rectangle substituted with
    replacement.
    (array[xul:xlr,yul:ylr])
     +---------------+
     |(xul, yul)     |
     |               |
     |     (xlr, ylr)|
     +---------------+
    """
    ri = 0
    for x in xrange(xul, xlr):
      array[x][yul:ylr] = replacement[ri]
      ri += 1
    return array
  replacerectinarray = staticmethod(replacerectinarray)
  def replaceinarray(haystack, predicate, replacement):
    """Return 'haystack' with 'predicate' matches replaced by 'replacement'"""
    return map(lambda item : {True: replacement, False: item}[predicate(item)],
               haystack)
  replaceinarray = staticmethod(replaceinarray)
  def sqrta(values):
    """Return the element-wise sqrt() of the given array."""
    return map(math.sqrt, values)
  sqrta = staticmethod(sqrta)
  def suba(lista, listb):
    """Subtract the values of a list from the values of another list."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.sub, lista, listb)
  suba = staticmethod(suba)
  def adda(lista, listb):
    """Add the values of a list to the values of another list."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.add, lista, listb)
  adda = staticmethod(adda)
  def diva(lista, listb):
    """Return the element-wise division of 'lista' by 'listb'."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.div, lista, listb)
  diva = staticmethod(diva)
  def mula(lista, listb):
    """Return the element-wise multiplication of 'lista' with 'listb'."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.mul, lista, listb)
  mula = staticmethod(mula)
  def powa(values, exponent):
    """Returns the element-wise exp('exponent') of the given array."""
    if isinstance(exponent, (list, tuple)):
      return map(lambda x, y : x ** y, values, exponent)
    return map(lambda x : x ** exponent, values)
  powa = staticmethod(powa)
  def treemap(fn, tree):
    """Applies `fn' to every value in `tree' which isn't a list and
    returns a list with the same shape as tree and the value of `fn'
    applied to the values in their place.

    """
    result = []
    for node in tree:
      if isinstance(node, (list, tuple)):
        result += [VLAB.treemap(fn, node)]
      else:
        result += [fn(node)]
    return result
  treemap = staticmethod(treemap)
  def makemaskeq(array, value):
    return VLAB.treemap(lambda x : x == value, array)
  makemaskeq = staticmethod(makemaskeq)
  def awhere(mask):
    """Returns the coordinates of the cells which evaluate true."""
    result = []
    if isinstance(mask, (list, tuple)):
      for i, cell in enumerate(mask):
        result += [[i] + sub for sub in VLAB.awhere(cell)
                   if isinstance(sub, (list, tuple))]
      return result
    else:
      if mask:
        return [[]]
      else:
        return []
  awhere = staticmethod(awhere)
  def aclone(tree):
    """Make a deep copy of `tree'."""
    if isinstance(tree, (list, tuple)):
      return list(map(VLAB.aclone, tree))
    return tree
  aclone = staticmethod(aclone)
  def make_chart(title, x_label, y_label, dataset):
    """Helper for creating Charts"""
    if sys.platform.startswith('java'):
      from org.jfree.chart import ChartFactory, ChartFrame, ChartUtilities
      from org.jfree.chart.axis import NumberTickUnit
      from org.jfree.chart.plot import PlotOrientation
      from org.jfree.data.xy import XYSeries, XYSeriesCollection
      chart = ChartFactory.createScatterPlot(title, x_label, y_label, dataset,
                                           PlotOrientation.VERTICAL, True, 
                                           True, False)
      plot = chart.getPlot()
      domain_axis = plot.getDomainAxis()
      domain_axis.setRange(-70, 70)
      range_axis = plot.getRangeAxis()
      range_axis.setRange(0.0, 0.4)
      range_axis.setTickUnit(NumberTickUnit(0.05))
      return chart
    else:
      # TODO: re-merge original python implementation
      raise Exception("original make_chart()")
      return None
  make_chart = staticmethod(make_chart)
  def make_dataset():
    """Dataset helper for supporting chart creation"""
    if sys.platform.startswith('java'):
      from org.jfree.data.xy import XYSeriesCollection
      return XYSeriesCollection()
    else:
      # TODO: re-merge original python implementation
      raise Exception("original make_dataset()")
      return None
  make_dataset = staticmethod(make_dataset)
  def plot(dataset, x, y, label):
    """plot dataset with x and y values and the given label"""
    if sys.platform.startswith('java'):
      from org.jfree.data.xy import XYSeries
      series = XYSeries(label)
      for next_x, next_y in zip(x, y):
        series.add(float(next_x), float(next_y))
      dataset.addSeries(series)
    else:
      # TODO: re-merge original python implementation
      raise Exception("original plot()")
      return None
  plot = staticmethod(plot)
  def save_chart(chart, filename):
    """save the generated chart in the given filename"""
    if sys.platform.startswith('java'):
      from java.io import File
      from org.jfree.chart import ChartUtilities
      ChartUtilities.saveChartAsPNG(File(filename), chart, 800, 600)
    else:
      # TODO: re-merge original python implementation
      raise Exception("save_chart")
      pass
  save_chart = staticmethod(save_chart)
  def maskand(array, mask):
    return map(lambda a, b : a & b, array, mask)
  maskand = staticmethod(maskand)
  def unique(array):
    """return unique values in the given input array"""
    sortedarray = list(array)
    sortedarray.sort()
    result = []
    def addifunique(x, y):
      if x != y:
        result.append(x)
      return y
    result.append(reduce(addifunique, sortedarray))
    return result
  unique = staticmethod(unique)
  def ravel(a):
    def add(a, i):
      if not(isinstance(a, (list, tuple))):
        a = [a]
      if isinstance(i, (list, tuple)):
        return a + ravel(i)
      else:
        return a + [i]
    return reduce(add, a)
  ravel = staticmethod(ravel)
  def min_l_bfgs_b(f, initial_guess, args=(), bounds=None, epsilon=1e-8):
    """The `approx_grad' is not yet implemented, because it's not yet
used."""
    if sys.platform.startswith('java'):
      from java.lang import ClassLoader
      from java.io import File
      for path in sys.path:
        if path.endswith('.jar'):
          libdir = File(path).getParent()
          break
      VLAB.logger.info(libdir)
      # hack taken from: http://blog.cedarsoft.com/2010/11/setting-java-library-path-programmatically
      System.setProperty("java.library.path", libdir)
      syspathfield = ClassLoader.getDeclaredField("sys_paths")
      syspathfield.setAccessible(True);
      syspathfield.set(None, None)
      from lbfgsb import DifferentiableFunction, FunctionValues, Bound, Minimizer
      class function(DifferentiableFunction):
        def __init__(self, function, args=()):
          self.function = function
          self.args = args
        def getValues(self, x):
          f = self.function(x, *self.args)
          g = VLAB.approx_fprime(x, self.function, epsilon, *self.args)
          return FunctionValues(f, g)
      min = Minimizer()
      min.setBounds([Bound(bound[0], bound[1]) for bound in bounds])
      result = min.run(function(f, args), VLAB.ravel(initial_guess))
      return result.point
  min_l_bfgs_b = staticmethod(min_l_bfgs_b)
  def approx_fprime(xk, f, epsilon, *args):
    f0 = f(xk, *args)
    grad = [0. for i in xrange(len(xk))]
    ei = [0. for i in xrange(len(xk))]
    for k in xrange(len(xk)):
      ei[k] = 1.0
      d = map(lambda i : i * epsilon, ei)
      grad[k] = (f(VLAB.adda(xk, d), *args) - f0) / d[k]
      ei[k] = 0.0
    return grad
  approx_fprime = staticmethod(approx_fprime)
  def doLibradtran(args):

    q = {
      'rpv_file':  args['rpv_file'],
      'outfile' :  args['outfile'],
      'infile'  :  args['infile'],
    }

    for k in args:
      q[k] = args[k]

    if args['scene'] == VLAB.K_LAEGERN:
      q['latitude']  = '47.481667'
      q['longitude'] = '-8.394722,17'
    elif args['scene'] == VLAB.K_THARANDT:
      q['latitude']  = '50.9676498'
      q['longitude'] = '-13.520354'

    if args['aerosol'] == VLAB.K_RURAL:
      q['aerosol'] = '1'
    elif args['aerosol'] == VLAB.K_MARITIME:
      q['aerosol'] = '2'
    elif args['aerosol'] == VLAB.K_URBAN:
      q['aerosol'] = '5'
    elif args['aerosol'] == VLAB.K_TROPOSPHERIC:
      q['aerosol'] = '6'

    # TODO: verify these
    if args['sensor'] ==   VLAB.K_SENTINEL2:
      q['wavelength'] = '443 2190'
    elif args['sensor'] == VLAB.K_SENTINEL3:
      q['wavelength'] = '400 2400'
    elif args['sensor'] == VLAB.K_MODIS:
      q['wavelength'] = '405 2155'
    elif args['sensor'] == VLAB.K_MERIS:
      q['wavelength'] = '412 900'
    elif args['sensor'] == VLAB.K_LANDSAT_OLI:
      q['wavelength'] = '433 2300'
    elif args['sensor'] == VLAB.K_LANDSAT_ETM:
      q['wavelength'] = '450 2300'

    if VLAB.osName().startswith('Windows'):
      q['solar_file'] = VLAB.expandEnv('%HOMEDRIVE%%HOMEPATH%\\.beam\\beam-vlab\\auxdata\\libRadtran_win32\\data\\solar_flux\\NewGuey2003.day')
    else:
      q['solar_file'] = VLAB.expandEnv('$HOME/.beam/beam-vlab/auxdata/libRadtran_lin64/data/solar_flux/NewGuey2003.dat')

    sdata = """
solar_file         %s
correlated_k       LOWTRAN
rte_solver         cdisort
rpv_file           %s
deltam             on
nstr               6
zout               TOA
output_user        lambda uu
quiet
""" % (q['solar_file'],
      q['rpv_file'])

    if 'aerosol' in q:
      sdata += "aerosol_default\n"
      sdata += "aerosol_haze       %s\n" % (q['aerosol']) 
    if 'O3' in q:
      sdata += "dens_column        O3 %s\n" % (q['O3']) 
    if 'CO2' in q:
      sdata += "co2_mixing_ratio   %s\n" % (q['CO2']) 
    if 'H2O' in q:
      sdata += "h2o_mixing_ratio   %s\n" % (q['H2O']) 
    if 'umu' in q:
      sdata += "umu                %s\n" % (q['umu']) 
    if 'phi' in q:
      sdata += "phi                %s\n" % (q['phi']) 
    if 'latitude' in q:
      sdata += "latitude           %s\n" % (q['latitude']) 
    if 'longitude' in q:
      sdata += "longitude          %s\n" % (q['longitude']) 
    if 'time' in q:
      sdata += "time               %s\n" % (q['time'])
    if 'sza' in q:
      sdata += "sza                %s\n" % (q['sza'])
    if 'phi0' in q:
      sdata += "phi0               %s\n" % (q['phi0'])
    if 'wavelength' in q:
      sdata += "wavelength         %s\n" % (q['wavelength'])

    fp = open(q['infile'], 'w')
    fp.write(sdata)
    fp.close()

    cmd = {
     'linux' : {
       'cwd'     : '$HOME/.beam/beam-vlab/auxdata/libRadtran_lin64/examples',
       'exe'     : '$HOME/.beam/beam-vlab/auxdata/libRadtran_lin64/bin/uvspec',
       'cmdline' : [],
       'stdin'   : q['infile'],
       'stdout'  : q['outfile'],
       'stderr'  : None,
       'env'     : None
     },
     'windows' : {
       'cwd'     : '%HOMEDRIVE%%HOMEPATH%\\.beam\\beam-vlab\\auxdata\\libRadtran_win32\\examples',
       'exe'     : '%HOMEDRIVE%%HOMEPATH%\\.beam\\beam-vlab\\auxdata\\libRadtran_win32\\uvspec.exe',
       'cmdline' : [],
       'stdin'   : q['infile'],
       'stdout'  : q['outfile'],
       'stderr'  : None,
       'env'     : None
    }
    }
    VLAB.logger.info('%s: spawning libradtran...' % VLAB.me())
    VLAB.doExec(cmd)
  doLibradtran = staticmethod(doLibradtran)


  ###########################################################################
  # 
  # Minimize_* functions...
  #
  # Nelder-Mead simplex minimization of a nonlinear (multivariate) function.
  # 
  # The programming interface is via the minimize() function; see below.
  # 
  # This code has been adapted from the C-coded nelmin.c which was
  # adapted from the Fortran-coded nelmin.f which was, in turn, adapted
  # from the papers
  # 
  #   J.A. Nelder and R. Mead (1965)
  #   A simplex method for function minimization.
  #   Computer Journal, Volume 7, pp 308-313.
  # 
  #   R. O'Neill (1971)
  #   Algorithm AS47. Function minimization using a simplex algorithm.
  #   Applied Statistics, Volume 20, pp 338-345.
  # 
  # and some examples are in
  # 
  #   D.M. Olsson and L.S. Nelson (1975)
  #   The Nelder-Mead Simplex procedure for function minimization.
  #   Technometrics, Volume 17 No. 1, pp 45-51.
  #    
  # For a fairly recent and popular incarnation of this minimizer,
  # see the amoeba function in the famous "Numerical Recipes" text.
  # 
  # P. Jacobs
  # School of Engineering, The University of Queensland
  # 07-Jan-04
  # 
  # Modifications by C. Schenkel
  # Netcetera
  # 31-Oct-13
  # 
  ###########################################################################
  def Minimize_create_new_point(c1, p1, c2, p2):
    """
    Create a new N-dimensional point as a weighting of points p1 and p2.
    """
    p_new = []
    for j in range(len(p1)):
      p_new.append(c1 * p1[j] + c2 * p2[j])
    return p_new
  Minimize_create_new_point = staticmethod(Minimize_create_new_point)

  def Minimize_take_a_step(smplx, Kreflect, Kextend, Kcontract):
    """
    Try to move away from the worst point in the simplex.
  
    The new point will be inserted into the simplex (in place).
    """
    i_low = smplx.lowest()
    i_high = smplx.highest()
    x_high = smplx.vertex_list[i_high]
    f_high = smplx.f_list[i_high]
    # Centroid of simplex excluding worst point.
    x_mid = smplx.centroid(i_high)
    f_mid = smplx.f(x_mid)
    smplx.nfe += 1
  
    # First, try moving away from worst point by
    # reflection through centroid
    x_refl = VLAB.Minimize_create_new_point(1.0+Kreflect, x_mid, -Kreflect, x_high)
    f_refl = smplx.f(x_refl)
    smplx.nfe += 1
    if f_refl < f_mid:
      # The reflection through the centroid is good,
      # try to extend in the same direction.
      x_ext = VLAB.Minimize_create_new_point(Kextend, x_refl, 1.0-Kextend, x_mid)
      f_ext = smplx.f(x_ext)
      smplx.nfe += 1
      if f_ext < f_refl:
        # Keep the extension because it's best.
        smplx.replace_vertex(i_high, x_ext, f_ext)
      else:
        # Settle for the original reflection.
        smplx.replace_vertex(i_high, x_refl, f_refl)
    else:
      # The reflection is not going in the right direction, it seems.
      # See how many vertices are better than the reflected point.
      count = 0
      for i in range(smplx.N+1):
        if smplx.f_list[i] > f_refl: count += 1
      if count <= 1:
        # Not too many points are higher than the original reflection.
        # Try a contraction on the reflection-side of the centroid.
        x_con = VLAB.Minimize_create_new_point(1.0-Kcontract, x_mid, Kcontract, x_high)
        f_con = smplx.f(x_con)
        smplx.nfe += 1
        if f_con < f_high:
          # At least we haven't gone uphill; accept.
          smplx.replace_vertex(i_high, x_con, f_con)
        else:
          # We have not been successful in taking a single step.
          # Contract the simplex about the current lowest point.
          smplx.contract_about_one_point(i_low)
      else:
        # Retain the original reflection because there are many
        # vertices with higher values of the objective function.
        smplx.replace_vertex(i_high, x_refl, f_refl)
    return
  Minimize_take_a_step = staticmethod(Minimize_take_a_step)
  
  def Minimize_minimize(f, x, dx=None, tol=1.0e-6,
         maxfe=300, n_check=20, delta=0.001,
         Kreflect=1.0, Kextend=2.0, Kcontract=0.5, args=()):
    """
    Locate a minimum of the objective function, f.
  
    Input:
    f   : user-specified function f(x)
    x   : list of N coordinates
    args  : Extra arguments passed to f, i.e. ``f(x, *args)''.
    dx    : list of N increments to apply to x when forming
          the initial simplex.  Their magnitudes determine the size
          and shape of the initial simplex.
    tol   : the terminating limit for the standard-deviation
          of the simplex function values.
    maxfe : maximum number of function evaluations that we will allow
    n_check : number of steps between convergence checks
    delta : magnitude of the perturbations for checking a local minimum
          and for the scale reduction when restarting
    Kreflect, Kextend, Kcontract: coefficients for locating the new vertex
  
    Output:
    Returns a tuple consisting of
    [0] a list of coordinates for the best x location,
      corresponding to min(f(x)),
    [1] the function value at that point,
    [2] a flag to indicate if convergence was achieved
    [3] the number of function evaluations and
    [4] the number of restarts (with scale reduction)
    """
    converged = 0
    N = len(x)
    if dx == None:
      dx = [0.1] * N
    smplx = Minimize_NMSimplex(x, dx, f, args)
  
    while (not converged) and (smplx.nfe < maxfe):
      # Take some steps and then check for convergence.
      for i in range(n_check):
        VLAB.Minimize_take_a_step(smplx, Kreflect, Kextend, Kcontract)
      # Pick out the current best vertex.
      i_best = smplx.lowest()
      x_best = list(smplx.get_vertex(i_best))
      f_best = smplx.f_list[i_best]
      # Check the scatter of vertex values to see if we are
      # close enough to call it quits.
      mean, stddev = smplx.f_statistics()
      if stddev < tol:
        # All of the points are close together but we need to
        # test more carefully.
        converged = smplx.test_for_minimum(i_best, delta)
        if not converged:
          # The function evaluations are all very close together
          # but we are not at a true minimum; rescale the simplex.
          smplx.rescale(delta)
    
    return x_best, f_best, converged, smplx.nfe, smplx.nrestarts
  Minimize_minimize = staticmethod(Minimize_minimize)
  
  
  #
  # from here to "VLAB end" is used only for testing
  #

  def fakebye(self, event):
    """fake beam callback for existing the program"""
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: exiting' % me)
    sys.exit()

  def fakeDoProcessing(self, params):
    """fake beam helper for driving processing in test mode"""
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: starting' % me)

    VLAB.logger.info('params: %s' % params)
    if params[VLAB.P_RTProcessor] == VLAB.K_DART:
      rtProcessor = DART()
    elif params[VLAB.P_RTProcessor] == VLAB.K_LIBRAT:
      rtProcessor = LIBRAT()
    elif params[VLAB.P_RTProcessor] == VLAB.K_DUMMY:
      rtProcessor = DUMMY()
    else:
      raise RuntimeError('unknown processor: <%s>' % params[VLAB.P_RTProcessor])
    rtProcessor.doProcessing(None, params)
    VLAB.logger.info('%s : %s' % (me, 'finished'))

  def fakerun(self, event):
    """fake beam callback for running a processor"""
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: starting' % me)
    params = {}
    for i in self.plst:
      if type(self.vmap[i]) == str:
         params[i] = self.cmap[i].text
      elif type(self.vmap[i]) == tuple:
         lst = self.vmap[i]
         params[i] = lst[self.cmap[i].selectedIndex]
    self.fakeDoProcessing(params)

  def fakebeam(self):
    """fake beam Swing GUI"""
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: starting' % me)

    from javax import swing
    from java  import awt

    v = 5; h = 10; self.window = swing.JFrame(VLAB.PROCESSOR_NAME)
    self.window.windowClosing = self.fakebye
    self.window.contentPane.layout = awt.BorderLayout()
    tabbedPane = swing.JTabbedPane()

    self.window.contentPane.add("Center", tabbedPane)
    for tabGroups in VLAB.model:
      for tabName in tabGroups:
        tab = swing.JPanel()
        tab.layout = swing.BoxLayout(tab, swing.BoxLayout.PAGE_AXIS)
        tabbedPane.addTab(tabName, tab)
        for group in tabGroups[tabName]:
          tab.add(swing.JLabel(''))
          p = swing.JPanel()
          p.layout = awt.GridLayout(0, 4);
          p.layout.vgap = v;p.layout.hgap = h
          for groupName in group:
            p.border = swing.BorderFactory.createTitledBorder(groupName)
            for groupTuple in group[groupName]:
              if len(groupTuple) == 4:
                (lbl, nm, typ, vals) = groupTuple

                p.add(swing.JLabel(lbl+':', swing.SwingConstants.RIGHT))
                self.vmap[nm] = vals
                if type(vals) == tuple:
                  self.cmap[nm] = swing.JComboBox(self.vmap[nm])
                else:
                  if type(typ) == tuple:
                    exec('self.cmap[nm] = swing.'+typ[0]+'(self.vmap[nm])')
                    self.cmap[nm].setEditable(typ[1])
                  else:
                    exec('self.cmap[nm] = swing.'+typ+'(self.vmap[nm])')
                self.plst.append(nm)
                p.add(self.cmap[nm])
              else:
                p.add(swing.JLabel(""))
                p.add(swing.JLabel(""))
          tab.add(p)
        # hack
        for i in range(50):
          tab.add(swing.Box.createVerticalGlue())
        tab.add(swing.JLabel(''))
        p = swing.JPanel()
        p.layout = awt.GridLayout(0, 8)
        p.layout.hgap = 4
        for i in range(1,5):
          p.add(swing.JLabel(""))
        p.add(swing.JButton("Run",   actionPerformed=self.fakerun))
        p.add(swing.JButton("Close", actionPerformed=self.fakebye))
        p.add(swing.JButton("Help"))
        tab.add(p)
    self.window.pack(); self.window.show()

  # TESTS
  def selftests(self):
    """run a pre-defined set of tests"""
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: starting' % me)

    # scenario 1
    params = {
      VLAB.P_3dScene             : 'RAMI', 
      VLAB.P_AtmosphereCO2       : '1.6', 
      VLAB.P_ViewingAzimuth      : '0.0', 
      VLAB.P_AtmosphereWater     : '0.0', 
      VLAB.P_OutputPrefix        : 'RAMI_', 
      VLAB.P_ViewingZenith       : '20.0', 
      VLAB.P_AtmosphereAerosol   : 'Rural', 
      VLAB.P_DHP_Zenith          : '20.0', 
      VLAB.P_SceneYW             : '100', 
      VLAB.P_DHP_3dScene         : 'RAMI', 
      VLAB.P_Bands               : '1, 2, 3, 4, 5, 6, 7, 8, 9, 10', 
      VLAB.P_OutputDirectory     : '', 
      VLAB.P_AtmosphereOzone     : '300', 
      VLAB.P_SceneLocFile        : '', 
      VLAB.P_SceneYC             : '100', 
      VLAB.P_DHP_OutputDirectory : '', 
      VLAB.P_DHP_OutputPrefix    : 'RAMI00_', 
      VLAB.P_RTProcessor         : VLAB.K_DUMMY,
      VLAB.P_SceneXW             : '50', 
      VLAB.P_IlluminationAzimuth : '0.0', 
      VLAB.P_Sensor              : 'MSI (Sentinel 2)', 
      VLAB.P_AsciiFile           : 'Yes', 
      VLAB.P_ImageFile           : 'Yes', 
      VLAB.P_SceneXC             : '-50', 
      VLAB.P_DHP_Y               : '0', 
      VLAB.P_DHP_X               : '0', 
      VLAB.P_DHP_ImageFile       : 'Yes', 
      VLAB.P_DHP_AsciiFile       : 'Yes', 
      VLAB.P_DHP_Resolution      : '100x100', 
      VLAB.P_DHP_Azimuth         : '0.0', 
      VLAB.P_IlluminationZenith  : '20.0', 
      VLAB.P_DHP_RTProcessor     : VLAB.K_DUMMY,
      VLAB.P_DHP_Height          : '0', 
      VLAB.P_DHP_Orientation     : '0', 
      VLAB.P_ScenePixel          : '8'
    }
    self.fakeDoProcessing(params)

    # scenario 2
    params[VLAB.P_RTProcessor] = VLAB.K_LIBRAT
    self.fakeDoProcessing(params)

    # scenario 3
    params[VLAB.P_RTProcessor] = VLAB.K_DART
    self.fakeDoProcessing(params)

#### VLAB end ################################################################

class Minimize_NMSimplex:
  """
  Stores the (nonlinear) simplex as a list of lists.

  In an N-dimensional problem, each vertex is a list of N coordinates
  and the simplex consists of N+1 vertices.
  """
  def __init__(self, x, dx, f, args):
    """
    Initialize the simplex.

    Set up the vertices about the user-specified vertex, x,
    and the set of step-sizes dx.
    f is a user-specified objective function f(x).
    """
    self.N = len(x)
    self.vertex_list = []
    self.f_list = []
    self.dx = list(dx)
    self.f = lambda x : f(x, *args)
    self.nfe = 0
    self.nrestarts = 0
    for i in range(self.N + 1):
      p = list(x)
      if i >= 1: p[i-1] += dx[i-1]
      self.vertex_list.append(p)
      self.f_list.append(f(p, *args))
      self.nfe += 1

  def rescale(self, ratio):
    """
    Pick out the current minimum and rebuild the simplex about that point.
    """
    i_min = self.lowest()
    for i in range(self.N):
      self.dx[i] *= ratio
    x = self.get_vertex(i_min)
    self.vertex_list = []
    self.f_list = []
    for i in range(self.N + 1):
      p = list(x)
      if i >= 1: p[i-1] += self.dx[i-1]
      self.vertex_list.append(p)
      self.f_list.append(self.f(p))
      self.nfe += 1
    self.nrestarts += 1
    return
  
  def get_vertex(self, i):
    return list(self.vertex_list[i])

  def replace_vertex(self, i, x, fvalue):
    self.vertex_list[i] = list(x)
    self.f_list[i] = fvalue
    return

  def lowest(self, exclude=-1):
    """
    Returns the index of the lowest vertex, excluding the one specified.
    """
    if exclude == 0:
      indx = 1
    else:
      indx = 0
    lowest_f_value = self.f_list[indx]
    for i in range(self.N + 1):
      if i == exclude: continue
      if self.f_list[i] < lowest_f_value:
        lowest_f_value = self.f_list[i]
        indx = i
    return indx

  def highest(self, exclude=-1):
    """
    Returns the index of the highest vertex, excluding the one specified.
    """
    if exclude == 0:
      indx = 1
    else:
      indx = 0
    highest_f_value = self.f_list[indx]
    for i in range(self.N + 1):
      if i == exclude: continue
      if self.f_list[i] > highest_f_value:
        highest_f_value = self.f_list[i]
        indx = i
    return indx

  def f_statistics(self):
    """
    Returns mean and standard deviation of the vertex fn values.
    """
    sum = 0.0
    for i in range(self.N + 1):
      sum += self.f_list[i]
    mean = sum / (self.N + 1)
    sum = 0.0
    for i in range(self.N +1):
      diff = self.f_list[i] - mean
      sum += diff * diff
    std_dev = math.sqrt(sum / self.N)
    return mean, std_dev

  def centroid(self, exclude=-1):
    """
    Returns the centroid of all vertices excluding the one specified.
    """
    xmid = [0.0]*self.N
    for i in range(self.N + 1):
      if i == exclude: continue
      for j in range(self.N):
        xmid[j] += self.vertex_list[i][j]
    for j in range(self.N):
      xmid[j] /= self.N
    return xmid
  
  def contract_about_one_point(self, i_con):
    """
    Contract the simplex about the vertex i_con.
    """
    p_con = self.vertex_list[i_con]
    for i in range(self.N + 1):
      if i == i_con: continue
      p = self.vertex_list[i]
      for j in range(self.N):
        p[j] = 0.5 * (p[j] + p_con[j])
      self.f_list[i] = self.f(p)
      self.nfe += 1
    return

  def test_for_minimum(self, i_min, delta):
    """
    Perturb the minimum vertex and check that it is a local minimum.
    """
    is_minimum = 1  # Assume it is true and test for failure.
    f_min = self.f_list[i_min]
    for j in range(self.N):
      # Check either side of the minimum, perturbing one
      # coordinate at a time.
      p = self.get_vertex(i_min)
      p[j] += self.dx[j] * delta
      f_p = self.f(p)
      self.nfe += 1
      if f_p < f_min:
        is_minimum = 0
        break
      p[j] -= self.dx[j] * delta * 2
      f_p = self.f(p)
      self.nfe += 1
      if f_p < f_min:
        is_minimum = 0
        break
    return is_minimum

#
# end Minimize_NMSimplex
#

#### DUMMY start #############################################################
class DUMMY:
  """A dummy processor for testing the VLAB plugin """
  def __init__(self):
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: constructor completed...' % me)
  def doProcessing(self, pm, args):
    """do processing for DUMMY processor"""
    me=self.__class__.__name__ +'::'+VLAB.me()

    VLAB.logger.info('%s: doExec() on %s' % (me, args))
    cmd = {
    'linux' : {
      'cwd'     : '$HOME/.beam/beam-vlab/auxdata/dummy_lin64/',
      'exe'     : '$HOME/.beam/beam-vlab/auxdata/dummy_lin64/dummy',
      'cmdline' : [ '-e', '1', '-r', '5' ],
      'stdin'   : None,
      'stdout'  : None,
      'stderr'  : None,
      'env'     : None
    },
    'windows' : {
      'cwd'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dummy_win32',
      'exe'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dummy_win32//dummy.exe',
      'cmdline' : [ '-e', '1', '-r', '5' ],
      'stdin'   : None,
      'stdout'  : None,
      'stderr'  : None,
      'env'     : None
    }
    }
    VLAB.doExec(cmd)

    if (pm != None):
      pm.beginTask("Computing BRF...", 10)
    # ensure at least 1 second to ensure progress popup feedback
    time.sleep(1)
    VLAB.logger.info('%s: finished running...' % me)

#### DUMMY end ###############################################################

#### DART start ##############################################################

##############################################################################
#
# Dart_* integration routines
#
# Authors: Nicolas Lauret, Fabian Schneider
#
# converted to "BEAM-limited jython" by Jason Brazile
#
#==============================================================================
#title            :Dart_DARTDao.py
#description      :This script can be used to recover DART images and spectral bands
#coding       :utf-8
#author           :Nicolas Lauret, (Fabian Schneider)
#date             :20130328
#version          :1.0
#usage            :see Dart_DartImages.py
#==============================================================================

def Dart_enum(**enums):
  return type('Enum', (), enums)

Dart_ProjectionPlane = Dart_enum(SENSOR_PLANE=0, ORTHOPROJECTION=1, BOA_PLANE=2)
Dart_DataUnit = Dart_enum(BRF_TAPP=0, RADIANCE=1)
Dart_DataLevel = Dart_enum(BOA=0, SENSOR=1, TOA=2, ATMOSPHERE_ONLY=3)
Dart_SpectralBandMode = Dart_enum(VISIBLE=0, VISIBLE_AND_THERMAL=1, THERMAL=2)

class Dart_DARTEnv : 
  simulationsDirectory = 'simulations'
  
  inputDirectory = 'input'
  outputDirectory = 'output'
  sequenceDirectory = 'sequence'
  
  libPhaseDirectory = 'lib_phase'
  bandRootDirectory = 'BAND'
  
  brfDirectory = 'BRF'
  tappDirectory = 'Tapp'
  radianceDirectory = 'Radiance'
  
  toaDirectory = 'TOA'
  radiativeBudgetDirectory = 'RADIATIVE_BUDGET'
  
  iterRootDirectory = 'ITER'
  couplDirectory = 'COUPL'
  sensorDirectory = 'SENSOR'
  
  #imageDirectoy = VLAB.path.join('IMAGES_DART', 'Non_Projected_Image')
  imageDirectoy = 'IMAGES_DART'
  orthoProjectedImageDirectory = 'IMAGE_PROJETEE'
  nonProjectedImageDirectory = 'Non_Projected_Image'
  
  brfFile = 'brf'
  radianceFile = 'radiance'
  radianceAtmoFile = 'AtmosphereRadiance'
  tappFile = 'tapp'
  radiativeBudgetFile = 'RadiativeBudget'
  
  propertiesFile = 'simulation.properties.txt'

# Class SpectralBand
class Dart_SpectralBand:
  def __init__(self) :
    self.index = 0
    self.lambdaMin = 0.
    self.lambdaMax = 0.
    self.mode = Dart_SpectralBandMode.VISIBLE
  
  def __str__(self):
    return 'Dart_SpectralBand(' + str(self.index) + ': [' + str(self.lambdaMin) + '; ' + str(self.lambdaMax) + ']; ' + str(self.mode) + ')'
  
  def getCentralWaveLength(self):
    return (self.lambdaMax + self.lambdaMin) / 2.
  
  def getBandWidth(self):
    return (self.lambdaMax - self.lambdaMin)

# Class Direction
class Dart_Direction :
  def __init__(self, theta, phi, solidAngle = 1, angularSector = 1, index = -1) :
    self.theta = theta
    self.phi = phi
    self.solidAngle = solidAngle
    self.angularSector = angularSector
    self.index = index
  
  def __str__(self):
    return 'Dart_Direction(' + str(self.theta) + ', ' + str(self.phi) + ')'
  
  def equal(self, dir2) :
    return abs(self.theta - dir2.theta) < 0.001 and abs(self.phi - dir2.phi) < 0.001

# Class BRF, Tapp, Radiance, per direction
class Dart_BRF :
  def __init__(self) :
    self.valeursParDirections = []
  
  def __str__(self):
    string = '<BDAverage'
    for (direction, moyenne) in self.valeursParDirections :
      string += '\n' + str(direction) + ' : ' + str(moyenne)
    string += '>'
    return string
  
  def addValeurDansDirection (self, direction, valeur) :
    self.valeursParDirections.append((direction, valeur))
  
  def getValueForDirection(self, directionRecherchee) :
    for (directionCourante, valeurCourante) in self.valeursParDirections :
      if (directionCourante.equal(directionRecherchee)) :
        return valeurCourante
  
  def getValue(self, indiceDirection) :
    return self.valeursParDirections[indiceDirection]

# Classe Images par direction
# FIXME Is this class still needed?
class Dart_Images :
  def __init__(self) :
    self.matrice2DparDirections = []
  
  def addMatriceDansDirection (self, direction, matrice) :
    self.matrice2DparDirections.append((direction, matrice))
  
  def getMatrice2D(self, directionRecherchee) :
    for (directionCourante, matriceCourante) in self.matrice2DparDirections :
      if (directionCourante.equal(directionRecherchee)) :
        return matriceCourante

# Classe Bilan Radiatif par type
# FIXME Is this class still needed?
class Dart_BilanRadiatif :
  def __init__(self) :
    self.matriceParType = []
  
  def getTypeBilan(self) :
    typeBilans = []
    for (typeBilanCourant, matriceCourante) in self.matriceParType :
      typeBilans.append(typeBilanCourant)
    return typeBilans
  
  def contientTypeBilan(self, typeBilan) :
    for (typeBilanCourant, matriceCourante) in self.matriceParType :
      if typeBilanCourant == typeBilan :
        return True
    return False

def Dart_getModeFromString(modeString):
  if modeString == 'T':
    return Dart_SpectralBandMode.THERMAL
  elif modeString == 'R+T':
    return Dart_SpectralBandMode.VISIBLE_AND_THERMAL
  else :
    return Dart_SpectralBandMode.VISIBLE

def Dart_getNumMaxIteration(path) :
  dirList = [f for f in VLAB.listdir(path) if (f[0] != '.') and (VLAB.path.isdir(VLAB.path.join(path, f))) and f.startswith(Dart_DARTEnv.iterRootDirectory) and f != Dart_DARTEnv.iterRootDirectory + 'X']
  
  maxIter = -1
  for rep in dirList :
    numIter = int(rep.split(Dart_DARTEnv.iterRootDirectory)[1])
  if numIter > maxIter :
    maxIter = numIter
  return maxIter

def Dart_readProperties(propertyPath):
  # propFile= file(propertyPath, "rU" )
  propFile= file(propertyPath, "r" )
  propDict= dict()
  for propLine in propFile:
    propDef= propLine.strip()
    if len(propDef) == 0:
      continue
    if propDef[0] in ( '!', '#' ):
      continue
    punctuation= [ propDef.find(c) for c in ':= ' ] + [ len(propDef) ]
    found= min( [ pos for pos in punctuation if pos != -1 ] )
    name= propDef[:found].rstrip()
    value= propDef[found:].lstrip(":= ").rstrip()
    propDict[name]= value
  propFile.close()
  return propDict;

def Dart_readPropertiesFromFile(stringFile):
  properties = StringIO.StringIO(stringFile)
  propertiesLines = properties.readlines()
  propertiesDico = {}
  for line in propertiesLines:
    prop,valeur = line.rstrip("\n\r\t").split(':')[0],line.rstrip("\n\r\t").split(':')[1]
    propertiesDico[prop] = valeur
  return propertiesDico

class Dart_DARTSimulation :
  def __init__(self, name, rootSimulations) :
    self.name = name
    self.rootSimulationsDirectory = rootSimulations
  
  def __str__(self):
    return self.name
  
  def getAbsolutePath(self) :
    return VLAB.path.join(self.rootSimulationsDirectory.getAbsolutePath(), self.name)
  
  def isValid(self) :
    if self.name != '' and VLAB.path.exists(self.getAbsolutePath()) :
      dirList=[f for f in VLAB.listdir(self.getAbsolutePath()) if (f[0] != '.') and (VLAB.path.isdir(VLAB.path.join(self.getAbsolutePath(), f)))]
      for directory in dirList:
        # On test si le dossier contient des dossiers output, input et/ou sequence
        # On relance la methode recursivement sur les eventuels autres dossiers (hors input, output et sequence)
        if (directory == Dart_DARTEnv.inputDirectory) or (directory == Dart_DARTEnv.outputDirectory) or (directory == Dart_DARTEnv.sequenceDirectory) :
          return True
    return False
  
  def getProperties(self):
    return Dart_readProperties(VLAB.path.join(self.getAbsolutePath(), Dart_DARTEnv.outputDirectory, Dart_DARTEnv.propertiesFile))
  
  def getSpectralBands(self):   
    properties = self.getProperties()
    nbSpectralBand = int(properties['phase.nbSpectralBands'])
    bandes = []
    for numBand in range(nbSpectralBand):
      bande = Dart_SpectralBand()
      bande.index = numBand
      bande.lambdaMin = float(properties['dart.band' + str(numBand) + '.lambdaMin'])
      bande.lambdaMax = float(properties['dart.band' + str(numBand) + '.lambdaMax'])
      bande.mode = Dart_getModeFromString(properties['dart.band' + str(numBand) + '.mode'])
      bandes.append(bande)
    return bandes
  
  def getSceneGeolocation(self) :
    properties = self.getProperties()
    return (float(properties['dart.maket.latitude']), float(properties['dart.maket.longitude']))
  
  def getDiscretizedDirection(self):    
    properties = self.getProperties()
    nbDD = int(properties['direction.numberOfDirections'])
    discretizedDirections = []
    for i in range(nbDD) :
      discretizedDirections.append(Dart_Direction(float(properties['direction.direction' + str(i) +'.thetaCenter']), float(properties['direction.direction' + str(i) +'.phiCenter']), solidAngle = float(properties['direction.direction' + str(i) +'.omega']), angularSector = int(properties['direction.direction' + str(i) +'.angularSector']), index = i))
    return discretizedDirections
  
  def getUserDirections(self) :
    properties = self.getProperties()
    nbUD = int(properties['direction.nbOfAdditionalSpots'])
    userDirections = []
    for i in range(nbUD) :
      additionalSpotIndex = int(properties['direction.additionalSpot' + str(i) + '.directionIndex'])
      userDirections.append(Dart_Direction(float(properties['direction.direction' + str(additionalSpotIndex) +'.thetaCenter']), float(properties['direction.direction' + str(additionalSpotIndex) +'.phiCenter']), index = additionalSpotIndex))
    return userDirections
  
  def listNomsSequences(self) :
    sequencePath = VLAB.path.join(self.getAbsolutePath(), Dart_DARTEnv.sequenceDirectory)
    listNom = []
    if (VLAB.path.exists(sequencePath)) :
      dirList=[f for f in VLAB.listdir(sequencePath) if (f[0] != '.') and (VLAB.path.isdir(VLAB.path.join(sequencePath, f)))]
      for sequence in dirList :
        nomSequence = sequence.rpartition('_')[0]
        if (not nomSequence in listNom) :
          listNom.append(nomSequence)
    listNom.sort()
    return listNom
  
  def listSequences(self, nom) :
    sequencePath = VLAB.path.join(self.getAbsolutePath(), Dart_DARTEnv.sequenceDirectory)
    listSequence = [f for f in VLAB.listdir(sequencePath) if (f[0] != '.') and (VLAB.path.isdir(VLAB.path.join(sequencePath, f))) and f.startswith(nom + '_')]
    
    listSequence.sort(lambda a, b: cmp(int(a[len(nom)+1:]),int(b[len(nom)+1:])))
    listSimuSequence = [Dart_DARTSimulation(VLAB.path.join(self.name, Dart_DARTEnv.sequenceDirectory, s), self.rootSimulationsDirectory) for s in listSequence]
    return listSimuSequence
  
  def getSequence(self, nomSequence, numeroSequence) :
    return Dart_DARTSimulation(VLAB.path.join(self.name, Dart_DARTEnv.sequenceDirectory, nomSequence + '_' + str(numeroSequence)), self.rootSimulationsDirectory);
  
  def getSubdirectory(self, spectralBand = 0, dataType = Dart_DataUnit.BRF_TAPP, level = Dart_DataLevel.BOA, iteration = 'last') :
    simulationPath = self.getAbsolutePath()
    bandPath = VLAB.path.join(simulationPath, Dart_DARTEnv.outputDirectory, Dart_DARTEnv.bandRootDirectory + str(spectralBand))
    brfPath = ''
    if dataType == Dart_DataUnit.RADIANCE :
      brfPath = VLAB.path.join(bandPath, Dart_DARTEnv.radianceDirectory)
      if (level == Dart_DataLevel.TOA or level == Dart_DataLevel.ATMOSPHERE_ONLY) :
        brfPath = VLAB.path.join(brfPath, Dart_DARTEnv.toaDirectory) # Radiance/TOA
      elif (level == Dart_DataLevel.SENSOR) :
        brfPath = VLAB.path.join(brfPath, Dart_DARTEnv.sensorDirectory) # Radiance/SENSOR
      elif (not iteration == 'last') :
        brfPath = VLAB.path.join(brfPath, iteration) # Radiance/ITER...
      else :
        # Recherche du dossier COUPL
        tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.couplDirectory) # Radiance/COUPL
        if (not VLAB.path.exists(tmpPath)) :
        # Recherche du dossier ITERX
          tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.iterRootDirectory + 'X') # Radiance/ITERX
          if (not VLAB.path.exists(tmpPath)) :
            # Recherche du dossier ITER de plus fort nombre
            tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.iterRootDirectory + str(Dart_getNumMaxIteration(brfPath)))
        brfPath = tmpPath;
    else :
      if (level == Dart_DataLevel.BOA) :
        # Test de la presence d'un dossier brf, sinon tapp
        brfPath = VLAB.path.join(bandPath, Dart_DARTEnv.brfDirectory)
        if (not VLAB.path.exists(brfPath)) :
          # on essaie le dossier tapp
          brfPath = VLAB.path.join(bandPath, Dart_DARTEnv.tappDirectory)
        if (not iteration == 'last') :
          brfPath = VLAB.path.join(brfPath, iteration) # BRF_TAPP/ITER...
        else :
          tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.couplDirectory) # BRF_TAPP/COUPL
          if (not VLAB.path.exists(tmpPath)) :
            # Recherche du dossier ITERX
            tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.iterRootDirectory + 'X') # BRF_TAPP/ITERX
            if (not VLAB.path.exists(tmpPath)) :
              # Recherche du dossier ITER de plus fort nombre
              tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.iterRootDirectory + str(Dart_getNumMaxIteration(brfPath)))
          brfPath = tmpPath;
      elif (level == Dart_DataLevel.SENSOR) :
        brfPath = VLAB.path.join(bandPath, Dart_DARTEnv.sensorDirectory) # SENSOR
      else :
        brfPath = VLAB.path.join(bandPath, Dart_DARTEnv.toaDirectory) # TOA
    
    return brfPath
  
  def getAveragePerDirection(self, spectralBand = 0, dataType = Dart_DataUnit.BRF_TAPP, level = Dart_DataLevel.BOA, iteration = 'last') :
    brfPath = self.getSubdirectory(spectralBand, dataType, level, iteration)
    brf = Dart_BRF()
    if (VLAB.path.exists(brfPath)) :
      # recherche du nom de fichier contenant le brf/tapp/radiance
      if (dataType == Dart_DataUnit.RADIANCE) :
        if (level == Dart_DataLevel.ATMOSPHERE_ONLY):
          brfPath = VLAB.path.join(brfPath, Dart_DARTEnv.radianceAtmoFile)
        else:
          brfPath = VLAB.path.join(brfPath, Dart_DARTEnv.radianceFile)
      else :
        tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.brfFile)
        if (not VLAB.path.exists(tmpPath)) :
          tmpPath = VLAB.path.join(brfPath, Dart_DARTEnv.tappFile)
        brfPath = tmpPath
      
      if (VLAB.path.exists(brfPath) and VLAB.path.isfile(brfPath)) :
        # lecture du fichier et remplissage du conteneur
        fichierBRF = open(brfPath, "r")
        lines = fichierBRF.readlines()
        for line in lines :
          split = line.split()
          brf.addValeurDansDirection(Dart_Direction(float(split[0]), float(split[1])), float(split[2]))
        fichierBRF.close()
    return brf
  
  def getImageNameInDirection(self, direction):
    properties = self.getProperties()
    imageNumber = properties['dart.products.images.number']
    imageIndex = 0
    found = False
    imageName = None
    while (imageIndex < imageNumber and not found):
      directionIndex = int(properties['dart.products.image' + str(imageIndex) + '.directionIndex'])
      if (directionIndex == direction.index):
        imageName = properties['dart.products.image' + str(imageIndex) + '.name']
        found = True
      else:
        imageIndex = imageIndex + 1
    return imageName
  
  def getImageMinMaxInDirection(self, direction, spectralBand = 0, dataType = Dart_DataUnit.BRF_TAPP, level = Dart_DataLevel.BOA, iteration = 'last', projectionPlane=Dart_ProjectionPlane.SENSOR_PLANE) :
    brfPath = self.getSubdirectory(spectralBand, dataType, level, iteration)
    minImage = None
    maxImage = None
    if (VLAB.path.exists(brfPath)) :
      imageDirPath = VLAB.path.join(brfPath, Dart_DARTEnv.imageDirectoy)
      if (projectionPlane == Dart_ProjectionPlane.ORTHOPROJECTION):
        imageDirPath = VLAB.path.join(imageDirPath, Dart_DARTEnv.orthoProjectedImageDirectory)
      elif (projectionPlane == Dart_ProjectionPlane.BOA_PLANE):
        imageDirPath = VLAB.path.join(imageDirPath, Dart_DARTEnv.nonProjectedImageDirectory)
      imageName = self.getImageNameInDirection(direction)
      mprPath = VLAB.path.join(imageDirPath, imageName + '.mpr')
      if (VLAB.path.exists(mprPath)):
        mprFile = open(mprPath, "r")
        lines = mprFile.readlines()
        for line in lines :
          split = line.split('=')
          if (split[0].startswith('MinMax')):
            splitSplit = split[1].split(':')
            minImage = float(splitSplit[0])
            maxImage = float(splitSplit[1])
            mprFile.close()
            return (minImage, maxImage)
        mprFile.close()
    return (minImage, maxImage)

  def getImageInDirection(self, direction, spectralBand = 0, dataType = Dart_DataUnit.BRF_TAPP, level = Dart_DataLevel.BOA, iteration = 'last', projectionPlane=Dart_ProjectionPlane.SENSOR_PLANE):
    brfPath = self.getSubdirectory(spectralBand, dataType, level, iteration)
    VLAB.logger.info('brfPath is %s' % brfPath)
    data = []
    if (VLAB.path.exists(brfPath)):
      imageDirPath = VLAB.path.join(brfPath, Dart_DARTEnv.imageDirectoy)
      if (projectionPlane == Dart_ProjectionPlane.ORTHOPROJECTION):
        imageDirPath = VLAB.path.join(imageDirPath, Dart_DARTEnv.orthoProjectedImageDirectory)
      elif (projectionPlane == Dart_ProjectionPlane.BOA_PLANE):
        imageDirPath = VLAB.path.join(imageDirPath, Dart_DARTEnv.nonProjectedImageDirectory)
      imageName = self.getImageNameInDirection(direction)
      # Recovery of the size of the image, in pixels
      mprPath = VLAB.path.join(imageDirPath, imageName + '.mpr')
      VLAB.logger.info('mprPath is %s' % mprPath)
      if (VLAB.path.exists(mprPath)):
        sizeX = 0
        sizeY = 0
        mprFile = open(mprPath, "r")
        lines = mprFile.readlines()
        for line in lines:
          split = line.split('=')
          if (split[0].startswith('Size')):
            splitSplit = split[1].split(' ')
            sizeX = int(splitSplit[0])
            sizeY = int(splitSplit[1])
        mprFile.close()
        mpSharpPath = VLAB.path.join(imageDirPath, imageName + '.mp#')
        if (VLAB.path.exists(mpSharpPath)):
          mpSharpFile = open(mpSharpPath, "rb")
          for x in range(sizeX):
            line = []
            for y in range(sizeY):
              line.append(struct.unpack('d', mpSharpFile.read(8))[0])
            data.append(line)
          mpSharpFile.close()
    return data

  def getIterMaxPath(self, spectralBand, iterX=True):
    """ Return the iterYY folder where YY is X is exist and iterX is True or the iterMax folder
    """
    # Get spectralBand folder
    iterFolder = VLAB.path.join(self.rootSimulationsDirectory, self.name, "output", spectralBand, "BRF")
    # Get list of iter folders
    iters = [ folder.lstrip('ITER') for folder in VLAB.listdir(iterFolder)
              if folder.startswith('ITER')]
    # Get iterX is exists otherwise select iterMax
    if iterX and 'X' in iters:
      return VLAB.path.join(iterFolder, "ITERX")
    else:
      if 'X' in iters:
        iters.remove('X')
      iters = [ int(i) for i in iters ]
      return VLAB.path.join(iterFolder, "ITER" + str(max(iters)))
        

class Dart_DARTRootSimulationDirectory :
  def __init__(self) :
    pass
  
  def samefile(self, path1, path2) :
    return VLAB.path.normcase(VLAB.path.normpath(path1)) == VLAB.path.normcase(VLAB.path.normpath(path2))
  
  def getAbsolutePath(self) :
    return VLAB.path.join(DART.SDIR)#, Dart_DARTEnv.simulationsDirectory)
  
  def getSimulationsList(self) :
    rootSimulationPath = self.getAbsolutePath()
    # Liste de tous les fichier/dossier non cache du repertoire
    dirList=[f for f in VLAB.listdir(rootSimulationPath) if (f[0] != '.') and (VLAB.path.isdir(VLAB.path.join(rootSimulationPath, f)))]
    dirList.sort()
    listSimulations = []
    for directory in dirList:
      self.listSimulationsInPath(listSimulations, directory)
    return listSimulations
  
  def listSimulationsInPath(self, listSimulations, path) :
    rootPath = path
    if not VLAB.path.isabs(path):
      rootPath = VLAB.path.join(self.getAbsolutePath(), path)
    dirList=[f for f in VLAB.listdir(rootPath) if (f[0] != '.') and (VLAB.path.isdir(VLAB.path.join(rootPath, f)))]
    dirList.sort()
    
    isSimulation = False    
    for directory in dirList:
      # On test si le dossier contient des dossiers output, input et/ou sequence
      # On relance la methode recursivement sur les eventuels autres dossiers (hors input, output et sequence)
      if (directory != Dart_DARTEnv.inputDirectory) and (directory != Dart_DARTEnv.outputDirectory) and (directory != Dart_DARTEnv.sequenceDirectory) :
        self.listSimulationsInPath(listSimulations, VLAB.path.join(path, directory))
      else :
        isSimulation = True
    if isSimulation :
      self.addSimulation(listSimulations, path)
  
  def addSimulation(self, listSimulations, path):
    simulation = self.getSimulation(path)
    if simulation.isValid() :
      listSimulations.append(simulation)
  
  def getSimulation(self, name) :
    rootSimulationPath = self.getAbsolutePath()
    rootPath = name
    if not VLAB.path.isabs(name):
      rootPath = VLAB.path.join(rootSimulationPath, name)
    simulationName = ''
    if not VLAB.path.exists(rootPath) :
      simulationName = name
    elif not self.samefile(rootPath, rootSimulationPath) :
      if VLAB.path.isdir(rootPath):
        (head, tail) = VLAB.path.split(rootPath)
        while not self.samefile(head, rootSimulationPath) and head !='' :
          simulationName = VLAB.path.join(tail, simulationName)
          (head, tail) = VLAB.path.split(head)
        else :
          if head != '' :
            simulationName = VLAB.path.join(tail, simulationName)
          else :
            simulationName = ''
    return Dart_DARTSimulation(simulationName, self)

class Dart_DARTDao :
  def __init__(self) :
    pass
  
  def getRootSimulationDirectory(self) :
    return Dart_DARTRootSimulationDirectory()
    
  def getSimulationsList(self) :
    return self.getRootSimulationDirectory().getSimulationsList()
  
  def getSimulation(self, name) :
    return self.getRootSimulationDirectory().getSimulation(name)

# FIXME Is this class still needed?
def Dart_Bandes(simulationProperties):
  propertiesDico = Dart_readPropertiesFromFile(simulationProperties)
  i = int(propertiesDico['phase.nbSpectralBands'])
  bandes = []
  for j in range(i):
    bande = {}
    bande['spectralBandwidth'] = round(float(propertiesDico['dart.band' + str(j) + '.lambdaMax']) - float(propertiesDico['dart.band' + str(j) + '.lambdaMin']),4)
    a = float(propertiesDico['dart.band' + str(j) + '.lambdaMin']) + bande['spectralBandwidth'] / 2
    bande['centralWavelength'] = round(float(a),4)
    bande['mode'] = propertiesDico['dart.band' + str(j) + '.mode']
    bandes.append(bande)
  return bandes

class Dart_dolibradtran:

  def main(q):
    # generate arguments for doLibradtran
    r = {
      'CO2'      : q[VLAB.P_AtmosphereCO2],
      'H2O'      : q[VLAB.P_AtmosphereWater],
      'O3'       : q[VLAB.P_AtmosphereOzone],
      'scene'    : q['obj'],
      'aerosol'  : q[VLAB.P_AtmosphereAerosol], 
      'sensor'   : q['wb'],
      'rpv_file' : '%s/rami.TOA/rpv.rami.libradtran.dat.all' % LIBRAT.SDIR,
      'infile'   : '%s/%s/input/%s' % (DART.SDIR, q['simulationName'], 'UVSPEC-in.txt'),
      'outfile'  : '%s/%s/output/%s' % (DART.SDIR, q['simulationName'], q['ipfile']),
    }
    # 
    # TODO: loop over something to produce umu, phi ,phi0, sza, etc.
    #
    VLAB.doLibradtran(r)

#==============================================================================
#title            :Dart_DartImages.py
#description      :This script can be used to write a *.bsq binary file and a corresponding ENVI header *.hdr from DART output images
#required scripts :Dart_DARTDao.py
#author           :Fabian Schneider, Nicolas Lauret
#date             :20130328
#version          :1.0
#usage            :see writeDataCube.py
#==============================================================================


class Dart_DartImages :
  def __init__(self) :
    pass

  # http://rightfootin.blogspot.ch/2006/09/more-on-python-flatten.html, based on Mike C. Fletcher's BasicTypes library
  def flatten( self, l, ltypes=(list, tuple)):
    ltype = type(l)
    l = list(l)
    i = 0
    while i < len(l):
      while isinstance(l[i], ltypes):
        if not l[i]:
          l.pop(i)
          i -= 1
          break
        else:
          l[i:i + 1] = l[i]
      i += 1
    return ltype(l)

  def getImagesBandsAsList( self, args ) :

    simulation = Dart_DARTDao().getSimulation( args['di_simName'] )

    imagesList = []
    spectralBandsList = []
    
    if args['di_isSeq']:
      sequencesList = simulation.listSequences( args['di_sequence'] )
      
      for sequence in sequencesList:
        # Recover the number of Spectral Band
        spectralBands = sequence.getSpectralBands()
    
        # Recover the Directions
        if args['ii_isUsrDir']:
          directions = sequence.getUserDirections()
        else:
          directions = simulation.getDiscretizedDirection()
  
        # Recover the images in the defined direction for each spectral band
        for band in spectralBands:
          VLAB.logger.info('spectralBand: %s' % (band))
          imagesList.append(sequence.getImageInDirection( directions[ args['ii_dirNum'] ], band.index, args['ii_dType'], args['ii_iLevel'], args['ii_iter'],  args['ii_projPlane'] ))
          spectralBandsList.append( band )
    else:
      # Recover the number of Spectral Band
      spectralBands = simulation.getSpectralBands()

      # Recover the Directions
      if args['ii_isUsrDir']:
        directions = sequence.getUserDirections()
      else:
        directions = simulation.getDiscretizedDirection()

      # Recover the images in the defined direction for each spectral band
      for band in spectralBands:
        imagesList.append(simulation.getImageInDirection( directions[ args['ii_dirNum'] ], band.index, args['ii_dType'], args['ii_iLevel'], args['ii_iter'],  args['ii_projPlane'] ))
        spectralBandsList.append( band )

    VLAB.logger.info('imagesList has %d entries' % len(imagesList))
    VLAB.logger.info('spectralBandsList has %d entries' % len(spectralBandsList))
    return imagesList, spectralBandsList
  
  def writeDataCube( self, args ) :
    
    imagesList, spectralBandsList = self.getImagesBandsAsList( args )
    
    # write BSQ binary file
    if (len(args['OutputDirectory']) == 0):
      args['OutputDirectory'] = '%s/%s/' % (DART.SDIR, args['di_simName'])
    fname = '%s/%s%s.bsq' % (args['OutputDirectory'], args['OutputPrefix'], args['di_outfname'])
    VLAB.logger.info('writing output bsq %s' % (fname))
    fout = open( fname, 'wb')
    flatarray = array('f', self.flatten( imagesList ))
    flatarray.tofile(fout)
    fout.close()

    # ENVI header information, more infos on: http://geol.hu/data/online_help/ENVI_Header_Format.html
    samples = 0
    lines = 0
    if (len(imagesList) > 0):
      lines = len( imagesList[0] )
      if (len(imagesList[0]) > 0):
        samples = len( imagesList[0][0] )
    bands = len( imagesList )
    headerOffset = 0
    fileType = 'ENVI Standard'
    dataType = 4
    interleave = 'BSQ'
    byteOrder = 0
    xStart = 1
    yStart = 1
    defaultBands = '{0}'
    wavelengthUnits = 'Nanometers'
    wavelength = ''
    fwhm = ''
    for band in spectralBandsList:
      wavelength += str( band.getCentralWaveLength()*1000 ) + ', '
      fwhm += str( band.getBandWidth()*1000 ) + ', '

    # write ENVI header file
    if (len(args['OutputDirectory']) == 0):
      args['OutputDirectory'] = '%s/%s/' % (DART.SDIR, args['di_simName'])
    fname = '%s/%s%s.hdr' % (args['OutputDirectory'], args['OutputPrefix'], args['di_outfname'])
    VLAB.logger.info('writing output hdr %s' % (fname))
    fout = open( fname, 'w')
    fout.writelines( 'ENVI ' + '\n' )
    fout.writelines( 'description = { ' + args['hi_desc'] + ' }' + '\n' )
    fout.writelines( 'samples = ' + str( samples ) + '\n' )
    fout.writelines( 'lines = ' + str( lines ) + '\n' )
    fout.writelines( 'bands = ' + str( bands ) + '\n' )
    fout.writelines( 'header offset = ' + str( headerOffset ) + '\n' )
    fout.writelines( 'file type = ' + fileType + '\n' )
    fout.writelines( 'data type = ' + str( dataType ) + '\n' )
    fout.writelines( 'interleave = ' + interleave + '\n' )
    fout.writelines( 'sensor type = ' + args['hi_sensor'] + '\n' )
    fout.writelines( 'byte order = ' + str( byteOrder ) + '\n' )
    fout.writelines( 'x start = ' + str( xStart ) + '\n' )
    fout.writelines( 'y start = ' + str( yStart ) + '\n' )
    fout.writelines( 'map info = {' + args['hi_projName'] + ', ' + str( args['hi_xRefPixl'] ) + ', ' + str( args['hi_yRefPixl'] ) + ', ' + str( args['hi_xRefCoord'] ) + ', ' + str( args['hi_yRefCoord'] ) + ', ' + str( args['hi_xPixSize'] ) + ', ' + str( args['hi_yPixSize'] ) + ', ' + args['hi_pixUnits'] + ' }' + '\n' )
    fout.writelines( 'default bands = ' + defaultBands + '\n' )
    fout.writelines( 'wavelength units = ' + wavelengthUnits + '\n' )
    fout.writelines( 'wavelength = { ' + wavelength[:-2] + ' }' + '\n' )
    fout.writelines( 'fwhm = { ' + fwhm[:-2] + ' }' + '\n' )
    fout.close()

#==============================================================================
#title            :writeDataCube.py
#description      :This script writes a *.bsq binary file and a corresponding ENVI header *.hdr from DART output images
#required scripts :Dart_DartImages.py + Dart_DARTDao.py
#author           :Fabian Schneider
#date             :20130328
#version          :1.0
#usage            :python writeDataCube.py
#==============================================================================

#
# End of Dart_* integration routines
#

class DART:
  if VLAB.osName().startswith('Windows'):
    SDIR=VLAB.expandEnv('%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dart_local/simulations')
  else:
    SDIR=VLAB.expandEnv('$HOME/.beam/beam-vlab/auxdata/dart_local/simulations')

  """Integration glue for calling external DART programs"""
  def __init__(self):
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: constructor completed...' % me)

  def dartSZA2libradtranSZA(self, sz, sa):
    """ DART uses theta value from 0 to 180 for zenith angle (half a circle) 
        and 0 to 360 for azimuth angle (full circle).
        In the upper sphere the zenith range is from 0 to 90.
        It is necessary to convert to libradtran zenith and azimuth range
        (-90 to 90 for zenith and 0-180 for azimuth).
    """
    if sa > 180.0:
      sz = -sz
      sa = sa - 180.0

    return sz, sa

  def _writeSeqFile(self, args):
    """ Write DART sequence file
    """
    sstr = """<?xml version="1.0" encoding="UTF-8"?>
<DartFile version="1.0">
    <DartSequencerDescriptor sequenceName="sequence;;%s">
        <DartSequencerDescriptorEntries>
            <DartSequencerDescriptorGroup groupName="group_%s">%s
            </DartSequencerDescriptorGroup>
        </DartSequencerDescriptorEntries>
        <DartSequencerPreferences dartLaunched="true"
            demGeneratorLaunched="false" directionLaunched="true"
            displayEnabled="true" hapkeLaunched="false"
            maketLaunched="true" numberParallelThreads="1"
            phaseLaunched="false" prospectLaunched="false"
            triangleFileProcessorLaunched="false"
            vegetationLaunched="false" zippedResults="false"/>
        <DartLutPreferences addedDirection="false" coupl="true"
            generateLUT="false" iterx="true" luminance="true"
            ordre="true" otherIter="true" phiMax="" phiMin=""
            reflectance="true" sensor="true" storeIndirect="false"
            thetaMax="" thetaMin="" toa="true"/>
    </DartSequencerDescriptor>
</DartFile>
"""
    dstr = """
                <DartSequencerDescriptorEntry args="%s"
                    propertyName="%s" type="enumerate"/>"""
    groupstr = ""
    for ent in args['entries']:
      (dargs, prop) = ent
      groupstr += (dstr % (dargs, prop))
    paramfilename = '%s/%s/%s' % (DART.SDIR, args['simulation'], args['fileName'])
    VLAB.logger.info('writing sequence paramfile "%s"' % paramfilename)
    fp = open(paramfilename, 'w')
    fp.write(sstr % (args['fileName'].rstrip(".xml"), args['fileName'].rstrip('.xml'), groupstr))
    fp.close()

  def doProcessing(self, pm, args):
    """do processing for DART processor"""
    me=self.__class__.__name__ +'::'+VLAB.me()

    if (pm != None):
      pm.beginTask("Computing BRF...", 10)
    # ensure at least 1 second to ensure progress popup feedback
    time.sleep(1)

    q = {
      'rpvfile'    : 'angles.rpv.2.dat',
      'simulation' : 'Unknown'
    }

    # overwrite defaults
    for a in args:
      if a == VLAB.P_3dScene:
        if args[a] == VLAB.K_RAMI:
          q['simulation'] = 'HET01_DIS_UNI_NIR_20'
        elif args[a] == VLAB.K_LAEGERN:
          q['simulation'] = 'Laegern'
        elif args[a] == VLAB.K_THARANDT:
          q['simulation'] = "Tharandt"
      elif a == VLAB.P_ViewingAzimuth:
        q['va'] = args[a]
      elif a == VLAB.P_ViewingZenith:
        q['vz'] = args[a]
      elif a == VLAB.P_IlluminationAzimuth:
        q['sa'] = args[a]
      elif a == VLAB.P_IlluminationZenith:
        q['sz'] = args[a]
      elif a == VLAB.P_Bands:
        q['bands'] = args[a]
      elif a == VLAB.P_ScenePixel:
        q['pixelSize'] = args[a]
      else:
        q[a] = args[a]

    # 1. Create the DART scene
    # 1. a. Copy DART original input file to a new folder
    q['simulationName'] = q['simulation'] + "_run"
    VLAB.copyDir(VLAB.path.join(DART.SDIR, q['simulation']), VLAB.path.join(DART.SDIR, q['simulationName']))

    # 1. b. Update the DART input files with parameters from GUI
    # In maket change the pixel size
    maket = VLAB.path.join(DART.SDIR, q['simulationName'], "input", "maket.xml")
    VLAB.XMLEditNode(maket, "CellDimensions", ["x", "z"], [q['pixelSize']]*2)
    # In direction change the Sun viewving angle
    direction = VLAB.path.join(DART.SDIR, q['simulationName'], "input", "directions.xml")
    VLAB.XMLEditNode(direction, "SunViewingAngles",
                     ["sunViewingAzimuthAngle", "sunViewingZenithAngle"],
                     [q['sa'], q['sz']])
    # In direction set ifCosWeighted and numberOfPropagationDirections attributes
    VLAB.XMLEditNode(direction, "Directions",
                     ["ifCosWeighted", "numberOfPropagationDirections"],
                     ["1", "100"])
    # In direction set viewing direction
    VLAB.XMLAddNode(direction, "Directions",
                    ["AddedDirections", "ZenithAzimuth", "Square", "DefineOmega"],
                    {"AddedDirections": {"attribute":
                                         ["directionType", "ifSquareShape", "imageDirection"],
                                         "value": ["0", "1", "1"],
                                         "parent":"Directions"},
                     "ZenithAzimuth": {"attribute":
                                       ["directionAzimuthalAngle", "directionZenithalAngle"],
                                       "value": [q['va'], q['vz']],
                                       "parent":"AddedDirections"},
                     "Square": {"attribute": ["widthDefinition"], "value": ["0"],
                                "parent":"AddedDirections"},
                     "DefineOmega": {"attribute": ["omega"], "value": ["0.001"], "parent":"Square"}
                    }
                   )
    # In phase force the radiance to be stored (could be not selected)
    phase = VLAB.path.join(DART.SDIR, q['simulationName'], "input", "phase.xml")
    VLAB.XMLEditNode(phase, "BrfProductsProperties", "luminanceProducts", "1")
    # In phase change the spectral bands
    if not q['simulationName'].startswith("HET01_DIS_UNI_NIR_20"):
      # TODO: You should addapt the getBandsFromGUI function to be consistent whith what you want
      bands = VLAB.getBandsFromGUI(q['bands'])
      VLAB.XMLReplaceNodeContent(phase, "SpectralIntervals", "SpectralIntervalsProperties",
                                 ["deltaLambda", "meanLambda"],
                                 bands, spectralBands=True)
      # In coeff_diff change multiplicativeFactorForLUT to allways be equal to "0"
      coeffDiff = VLAB.path.join(DART.SDIR, q['simulationName'], "input", "coeff_diff.xml")
      VLAB.XMLEditNode(coeffDiff, "LambertianMulti", "useMultiplicativeFactorForLUT", "0", multiple=False)
    else:
      bands = [None]
      VLAB.logger.info("*********************************************************************************")
      VLAB.logger.info("* WARNING: RAMI scene cannot be multiband according to the RAMI scene definition. *")
      VLAB.logger.info("*********************************************************************************")

    # 2. a. Run DART direction module
    cmd = {
      'linux' : {
        'cwd'     : '$HOME/.beam/beam-vlab/auxdata/dart_lin64/tools/linux',
        'exe'     : '/bin/bash',
        'cmdline' : ['$HOME/.beam/beam-vlab/auxdata/dart_lin64/tools/linux/dart-directions.sh', q['simulationName']],
        'stdin'   : None,
        'stdout'  : None,
        'stderr'  : None,
        'env'     : None,
        },
      'windows'   : {
        'cwd'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dart_win32/tools/windows',
        'exe'     : 'cmd.exe',
        'cmdline' : ['/c', '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dart_win32/tools/windows/dart-directions.bat', q['simulationName']],
        'stdin'   : None,
        'stdout'  : None,
        'stderr'  : None,
        'env'     : None,
       }
    }
    VLAB.logger.info('command: %s' % cmd)
    VLAB.doExec(cmd)

    # 2. b. Get the first two cols (sz, sa) where sz < 70 deg - use as sequence
    zlist = []; alist = []
    for row in VLAB.valuesfromfile('%s/%s/output/directions.txt' % (DART.SDIR, q['simulationName'])):
      (sz, sa, _, _)  = row
      if (sz < 70.0):
        zlist.append(sz); alist.append(sa)
    zparam = ";".join('%.2f' % x for x in zlist)
    aparam = ";".join('%.2f' % x for x in alist)

    # 2. c. Edit simulation directions.xml file to set the number of direcions to 200
    VLAB.XMLEditNode(direction, "Directions", "numberOfPropagationDirections", "200")
 
    # 2. d. Write sequence file
    seqparams = {'fileName'   : 'SunDirections.xml',
                 'simulation' : q['simulationName'],
                 'entries'    : [[aparam, 'Directions.SunViewingAngles.sunViewingAzimuthAngle'],
                                 [zparam, 'Directions.SunViewingAngles.sunViewingZenithAngle']
                                ]
                }
    self._writeSeqFile(seqparams)


    # 3. a. Run DART to pregenerate phase, maket, (3D obj..) before the sequence
    cmd = {
      'linux' : {
        'cwd'     : '$HOME/.beam/beam-vlab/auxdata/dart_lin64/tools/linux',
        'exe'     : '/bin/bash',
        'cmdline' : ['$HOME/.beam/beam-vlab/auxdata/dart_lin64/tools/linux/dart-full.sh', q['simulationName']],
        'stdin'   : None,
        'stdout'  : None,
        'stderr'  : None,
        'env'     : None,
        },
      'windows'   : {
        'cwd'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dart_win32/tools/windows',
        'exe'     : 'cmd.exe',
        'cmdline' : ['/c', '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dart_win32/tools/windows/dart-full.bat', q['simulationName']],
        'stdin'   : None,
        'stdout'  : None,
        'stderr'  : None,
        'env'     : None,
       }
    }
    VLAB.logger.info('command: %s' % cmd)
    VLAB.doExec(cmd)

    # 3. b. Run the sun directions sequence with number of directions=200
    cmd = {
      'linux' : {
        'cwd'     : '$HOME/.beam/beam-vlab/auxdata/dart_lin64/tools/linux',
        'exe'     : '/bin/bash',
        'cmdline' : ['$HOME/.beam/beam-vlab/auxdata/dart_lin64/tools/linux/dart-sequence.sh', q['simulationName'], seqparams['fileName']],
        'stdin'   : None,
        'stdout'  : None,
        'stderr'  : None,
        'env'     : None,
        },
      'windows'   : {
        'cwd'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dart_win32/tools/windows',
        'exe'     : 'cmd.exe',
        'cmdline' : ['/c', '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/dart_win32/tools/windows/dart-sequence.bat', q['simulationName'], seqparams['fileName']],
        'stdin'   : None,
        'stdout'  : None,
        'stderr'  : None,
        'env'     : None,
       }
    }
    VLAB.logger.info('command: %s' % cmd)
    VLAB.doExec(cmd)

    # 4.a. Select the values of sz < 70 to get BANDX.angles.rpv.2.dat
    # for each spectral band in the DART simulation
    rpvlist = None
    rpvFiles = []
    seqDir = VLAB.path.join(DART.SDIR, q['simulationName'], "sequence")
    for bandNumber in xrange(len(bands)):
      rpvlist = []
      for seqSubDir in VLAB.listdir(seqDir):
        dartBRF = VLAB.path.join(Dart_DARTSimulation(seqSubDir, seqDir).getIterMaxPath('BAND' + str(bandNumber)), "brf")
        dirFile = VLAB.path.join(seqDir, seqSubDir, "input", "directions.xml")
        sa = float(VLAB.XMLGetNodeAttributeValue(dirFile, "SunViewingAngles", "sunViewingAzimuthAngle"))
        sz = float(VLAB.XMLGetNodeAttributeValue(dirFile, "SunViewingAngles", "sunViewingZenithAngle"))
        for row in VLAB.valuesfromfile(dartBRF):
          (vz, va, brf)  = row
          if (vz < 70.0):
            # Convert DART output to librattran input
            vz, va = self.dartSZA2libradtranSZA(vz, va)
            rpvlist.append('%.2f\t%.2f\t%2f\t%2f\t%2f' % (sz, sa, vz, va, brf))
      # write the libradtran input file
      fname = VLAB.path.join(DART.SDIR,
                             q['simulationName'],
                             "BAND" + str(bandNumber) + "." + q['rpvfile'])
      rpvFiles.append(fname)
      fp = open(fname, "w")
      fp.write("\n".join(rpvlist))
      fp.close()
    q['rpvFiles'] = rpvFiles

    # 4.b. Run libradtran
    # FIXME: this needs to eventually use the same routine that librat uses
    dolibradtran = Dart_dolibradtran()
    dolibradtran.main(q)

    #
    # collect result into a consolidated data cube
    #
    dargs = {
      'di_simName'   : q['simulationName'],
      'di_isSeq'     : False,
      'di_sequence'  : 'sequence_apex',
      'di_outfname'  : 'DartOutput',
      'ii_iLevel'    : Dart_DataLevel.BOA,
      'ii_isUsrDir'  : False,
      'ii_dirNum'    : 0,
      'ii_dType'     : Dart_DataUnit.RADIANCE,
      'ii_iter'      : 'last',
      'ii_projPlane' : Dart_ProjectionPlane.SENSOR_PLANE,
      'hi_desc'      : 'some text',
      'hi_sensor'    : 'APEX',
      'hi_projName'  : 'Arbitrary',
      'hi_xRefPixl'  : 1,
      'hi_yRefPixl'  : 1,
      'hi_xRefCoord' : '2669660.0000',
      'hi_yRefCoord' : '1259210.0000',
      'hi_xPixSize'  : 1,
      'hi_yPixSize'  : 1,
      'hi_pixUnits'  : 'units=Meters'
    }
    # merge in data cube arguments
    for d in dargs:
      args[d] = dargs[d]

    VLAB.logger.info('args is %s' % args)
    Dart_DartImages().writeDataCube( args )
    VLAB.logger.info('%s: done...' % me)

#### DART end ################################################################

#### LIBRAT start ############################################################
#
# Librat_ integration routines contributed by Mat Disney and Philip Lewis 
# adapted for BEAM-embedded jython by Cyrill Schenkel and Jason Brazile
#
################
class Librat_dobrdf:
  def _writeCamFile(self, camFile, args):

    # defaults
    q = {
      'cam_camera'      : 'simple camera',
      'perspective'     : False,
      'result_image'    : 'result.image',
      'integral_mode'   : 'scattering order',
      'integral'        : 'result',
      'vz'              : 0,
      'va'              : 0,
      'twist'           : 0,
      'look'            : (0., 0., 0.),
      'ideal'           : (100., 100.),
      'boom'            : 1000.,
      'npixels'         : 100000,
      'rpp'             : 1,
      'fov'             : 180
    }

    # overwrite defaults
    for a in args:
      q[a] = args[a]

    cdata = ' camera { \n' \
+ ' %s = "%s";\n' %('camera.name', q['cam_camera']) \
+ ' %s = %s;\n'   %('geometry.zenith', q['vz']) \
+ ' %s = %s;\n'   %('geometry.azimuth', q['va']) \
+ ' %s = "%s";\n' %('result.image', q['result_image']) \
+ ' %s = "%s";\n' %('result.integral.mode', q['integral_mode']) \
+ ' %s = "%s";\n' %('result.integral', q['integral']) \
+ ' %s = %s;\n'   %('samplingCharacteristics.nPixels', q['npixels']) \
+ ' %s = %s;\n'   %('samplingCharacteristics.rpp', q['rpp']) \
+ ' %s = %s;\n'   %('geometry.idealArea', ', '.join(map(str,map('%.1f'.__mod__, q['ideal'])))) \
+ ' %s = %s;\n'   %('geometry.lookat', ', '.join(map(str,map('%.1f'.__mod__,q['look_xyz'])))) \
+ ' %s = %s;\n'   %('geometry.boomlength', q['boom'])

    if q['perspective']:
      cdata += ' %s = %s;\n' %('geometry.perspective', q['perspective'])
    if q['twist']:
      cdata += ' %s = %s;\n' %('geometry.twist', q['twist'])
    if q['fov']:
      cdata += ' %s = %s;\n' %('geometry.fieldOfView', q['fov'])
    if 'lidar' in q:
      if q['lidar']:
        cdata += ' %s = %s;\n' %('lidar.binStep', q['binStep']) \
               + ' %s = %s;\n' %('lidar.binStart', q['binStart']) \
               + ' %s = %s;\n' %('lidar.nBins', q['nBins'])
    cdata += '}'

    if q['fov'] and q['ideal']:
      raise Exception("camera: setting both 'ideal area' and 'fov' is not allowed")

    fp = open(camFile, 'w')
    try:
      fp.write(cdata)
    finally:
      fp.close()

  def _writeLightFile(self, lightFile, args):

    # defaults
    q = {
      'light_camera' : 'simple illumination',
       'sz'          : 0.,
       'sa'          : 0.,
       'twist'       : 0.,
    }

    # overwrite detaults
    for a in args:
      q[a] = args[a]

    ldata = ' camera { \n' \
+ ' %s = "%s";\n'   %('camera.name', q['light_camera']) \
+ ' %s = "%.1f";\n' %('geometry.zenith', float(q['sz'])) \
+ ' %s = "%.1f";\n' %('geometry.azimuth', float(q['sa'])) \
+ ' %s = "%.1f";\n' %('geometry.twist', float(q['twist']))

    key = "sideal"
    if key in q: ldata += '%s = %s\n' %('geometry.idealArea', ', '.join(map(str, map('%.1f'.__mod__, q[key]))))
    key = "slook_xyz"
    if key in q: ldata += '%s = %s\n' %('geometry.lookat', ', '.join(map(str, map('%.1f'.__mod__, q[key]))))
    key = "sboom"
    if key in q: ldata += '%s = %s\n' %('geometry.boom', q[key])
    key = "sperspective"
    if key in q: ldata += '%s = %s\n' %('geometry.perspecitive', q[key])
    key = "sfov"
    if key in q: ldata += '%s = %s\n' %('geometry.fov', q[key])
    ldata += '}'
    fp = open(lightFile, 'w')
    try:
      fp.write(ldata)
    finally:
      fp.close()


  def _writeInputFile(self, inpFile, lightFile, camFile):
    idata = '14' \
   + ' ' \
   + VLAB.getFullPath(camFile) \
   + ' ' \
   + VLAB.getFullPath(lightFile)
    fp = open(inpFile, 'w')
    try:
      fp.write(idata)
    finally:
      fp.close()

  def _writeGrabFile(self, grabFile, args):
    gFilePath = VLAB.getFullPath(grabFile)
    if VLAB.osName().startswith('Windows'): 
      gFilePath = gFilePath.replace("\\", "\\\\")
    # 'cwd'     : '$HOME/.beam/beam-vlab/auxdata/librat_scenes',
    gdata = """cmd = {
  'linux' : {
    'cwd'     : '$HOME/.beam/beam-vlab/auxdata/librat_scenes',
    'exe'     : '$HOME/.beam/beam-vlab/auxdata/librat_lin64/src/start/start',
    'cmdline' : ['-RATv', '-m', '%s', '-RATsensor_wavebands', '$HOME/.beam/beam-vlab/auxdata/librat_scenes/%s', '$HOME/.beam/beam-vlab/auxdata/librat_scenes/%s' ],
    'stdin'   : '%s',
    'stdout'  : '%s',
    'stderr'  : '%s',
    'env'     : {
      'BPMS'  : '$HOME/.beam/beam-vlab/auxdata/librat_lin64/',
  'LD_LIBRARY_PATH' :  '$HOME/.beam/beam-vlab/auxdata/librat_lin64/src/lib',
    }},
  'windows'   : {
    'cwd'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_scenes',
    'exe'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_win32/src/start/ratstart.exe',
    'cmdline' : ['-RATv', '-m', '%s', '-RATsensor_wavebands', '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_scenes/%s', '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_scenes/%s' ],
    'stdin'   : '%s',
    'stdout'  : '%s',
    'stderr'  : '%s',
    'env'     : {
      'BPMS'  : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_win32'
   }}
}
"""
    # hack to allow replacing only %s
    escaped = gdata.replace("%%","\x81\x81").replace("%H","\x81H").replace("%/","\x81/")
    replaced = escaped % \
   (args['sorder'], args['wbfile'], args['objfile'], gFilePath+'.inp', gFilePath+'.out.log', gFilePath+'.err.log', \
    args['sorder'], args['wbfile'], args['objfile'], gFilePath+'.inp', gFilePath+'.out.log', gFilePath+'.err.log')
    gdata = replaced.replace("\x81", "%")
    gdata += 'VLAB.doExec(cmd)\n'
    try:
      fp = open(gFilePath, 'w')
      fp.write(gdata)
    finally:
      fp.close()

  def main(self, args):
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('======> %s' % me)
    for a in args:
      VLAB.logger.info("%s -> %s" % (a, args[a]))

    # set up defaults
    q = {
      'INFINITY'        : 1000000,
      'lightfile'       : 'light.dat',
      'camfile'         : 'camera.dat',
      'wbfile'          : 'wb.test.dat',
      'anglefile'       : 'angles.dat',
      'objfile'         : 'veglab_test.obj',

      'blacksky'        : '',
      'npixels'         : 1000000,
      'image'           : 1,
      'rpp'             : 1,
      'fov'             : False,
      'sorder'          : 100,
      'nice'            : None,
      'boom'            : 100000,
      'ideal'           : False,
      'samplingPattern' : False,
      'mode'            : 'scattering order',
      'opdir'           : 'brdf',
      'result_root'     : 'result',
      'camera_root'     : 'camera',
      'light_root'      : 'light',
      'grabme_root'     : 'grabme',

      'look_xyz'        : (150, 150, 35),
      'location'        : False,
      'vz'              : -1,
      'va'              : -1,
      'sz'              : -1,
      'sa'              : -1,
    }
    for a in args:
      if a == 'look':
        q['look_xyz'] = args[a]
      elif a == 'result':
        q['result_root'] = args[a]
      elif a == 'camera':
        q['camera_root'] = args[a]
      elif a == 'light':
        q['light_root'] = args[a]
      elif a == 'grabme':
        q['grabme_root'] = args[a]
      elif a == 'wb':
        q['wbfile'] = args[a]
      elif a == 'angles':
        q['anglefile'] = args[a]
      elif a == 'obj':
        q['objfile'] = args[a]
      else:
        q[a] = args[a]

    VLAB.mkDirPath('%s/%s' %(LIBRAT.SDIR, q['opdir']))

    angfp = VLAB.checkFile('%s/%s' % (LIBRAT.SDIR, q['anglefile']))
    wbfp  = VLAB.checkFile('%s/%s' % (LIBRAT.SDIR, q['wbfile']))
    objfp = VLAB.checkFile('%s/%s' % (LIBRAT.SDIR, q['objfile']))

    # vz va sz sa
    ang = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, q['anglefile']))

    if 'lookFile' in q:
      lookfp = VLAB.checkFile('%s/%s' % (LIBRAT.SDIR, q['lookFile']))
      q['look_xyz'] = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, q['lookFile']))

    if len(q['look_xyz']) == 3:
      q['look_xyz'] = ((q['look_xyz']),)

    if len(ang) < 4 or (len(ang) > 4 and len(ang[0]) != 4):
      sys.stderr.write("%s: wrong number of fields (%i) in %s - should be 4\n"%(me, len(ang[1]), q['anglefile']))
      raise Exception()

    if len(ang) == 4:
      ang = ((ang),)

    for n, a in enumerate(ang):
      q['vz'] = a[0]
      q['va'] = a[1]
      q['sz'] = a[2]
      q['sa'] = a[3]

      for ll, look in enumerate(q['look_xyz']):
        lightfile = VLAB.fPath('%s/%s' % (LIBRAT.SDIR, q['opdir']), q['light_root'] + '_sz_' + str(q['sz']) + '_sa_' + str(q['sa']) + '_dat')
        ligfp = VLAB.openFileIfNotExists(lightfile)
        if ligfp != None:
          nq = {
            'sz' : q['sz'],
            'sa' : q['sa'],
          }
          self._writeLightFile(lightfile, nq)

        rooty = q['objfile'] + '_vz_' + str(q['vz']) + '_va_' + str(q['va']) + '_sz_' + str(q['sz']) + '_sa_' + str(q['sa']) + '_xyz_' + str(look[0]) + '_' + str(look[1]) + '_' + str(look[2]) + '_wb_' + q['wbfile']
        grabme = VLAB.fPath('%s/%s' % (LIBRAT.SDIR, q['opdir']), q['grabme_root'] + '.' + rooty)
        if 'dhp' in q:
          location = list(look)
          look[2] = q['INFINITY']
          sampling = 'circular'

        if not VLAB.fileExists(grabme):
          grabfp = VLAB.openFileIfNotExists(grabme)
          q['grabme_log'] = grabme + '.log'
          q['camfile'] = VLAB.fPath('%s/%s' % (LIBRAT.SDIR, q['opdir']), q['camera_root'] + '.' + rooty)
          q['oproot'] = VLAB.fPath('%s/%s' % (LIBRAT.SDIR, q['opdir']), q['result_root'] + '.' + rooty)
          nq = {
            'name'             : 'simple camera',
            'vz'               : q['vz'],
            'va'               : q['va'],
            'integral'         : q['oproot'],
            'integral_mode'    : q['mode'],
            'npixels'          : q['npixels'],
            'rpp'              : q['rpp'],
            'boom'             : q['boom'],
            'ideal'            : q['ideal'],
            'look_xyz'         : look,
            'location'         : q['location'],
            'fov'              : q['fov'],
            'samplingPattern'  : q['samplingPattern'],
            'file'             : q['camfile'],
            'grabme_log'       : q['grabme_log'],
            'camfile'          : q['camfile']
          }
          if 'image' in q:
            if 'hips' in q:
              nq['imfile'] = q['oproot'] + '.hips'
            else:
              nq['imfile'] = q['oproot'] + '.bim'
          else:
            nq['imfile'] = 'image'
          self._writeCamFile(nq['camfile'], nq)
          self._writeInputFile(grabme + '.inp', lightfile, nq['camfile'])
          nq = {
            'blacksky' : q['blacksky'],
            'sorder'   : q['sorder'],
            'objfile'  : q['objfile'],
            'wbfile'   : q['wbfile']
          }
          self._writeGrabFile(grabme, nq)
          execfile(grabme)
    VLAB.logger.info('done')

#############################################################################

class Librat_dolibradtran:

  def main(self, q):
    # initial arguments for doLibradtran
    r = {
      'CO2'      : q[VLAB.P_AtmosphereCO2],
      'H2O'      : q[VLAB.P_AtmosphereWater],
      'O3'       : q[VLAB.P_AtmosphereOzone],
      'scene'    : q['obj'],
      'aerosol'  : q[VLAB.P_AtmosphereAerosol], 
      'sensor'   : q['sensor'],
      'rpv_file' : VLAB.path.join(LIBRAT.SDIR, q['rpv'])
    }

    VLAB.mkDirPath('%s/%s' % (LIBRAT.SDIR, q['opdir']))
    angfp = VLAB.checkFile('%s/%s' % (LIBRAT.SDIR, q['angfile']))
    wbfp = VLAB.checkFile('%s/%s' % (LIBRAT.SDIR, q['wbfile']))
    rpvfp = VLAB.checkFile(r['rpv_file'])

    rpv = VLAB.valuesfromfile(r['rpv_file'])
    if len(rpv[0]) != 4: # length of index 1 because index 0 is heading
      sys.stderr.write("%s: rpv file %s wrong no. of cols (should be 4: lambda (nm), rho0, k, theta\n" % (sys.argv[0], rpv_file))
      raise Exception()

    angt = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, q['angfile']))
    wb = [i[1] for i in VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, q['wbfile']))]

    nbands = len(wb)
    wbstep = 1
    if q['v']:
      sys.stderr.write('%s: wbmin = %i, wbmax = %i, wbstep = %i\n'%(sys.argv[0],min(wb),max(wb),wbstep))

    if q['v']:
      # only do all angles if time not specified, if time specified use that to get sza and phi0
      if q['lat'] and ['lon'] and ['time']:
        sys.stderr.write("%s: doing lat lon time, not using sun angles in file %s\n" % (sys.argv[0], q['anglefile']))

    # check for op file if required
    if not VLAB.fileExists('%s/%s' % (LIBRAT.SDIR, q['plotfile'])):
      VLAB.logger.info('%s/%s' % (LIBRAT.SDIR, q['plotfile']))
      plotfilefp = VLAB.openFileIfNotExists('%s/%s' % (LIBRAT.SDIR, q['plotfile']))
    else:
      sys.stderr.write('%s: plotfile %s already exists - move/delete and re-run\n'%(sys.argv[0],q['plotfile']))
      raise Exception('%s: plotfile %s already exists - move/delete and re-run\n'%(sys.argv[0],q['plotfile']))

    # loop over angles and do all wb at once for each angle
    for a, aa in enumerate(angt):
      # need special case for 1 angle only, in which case angt = 4
      if len(angt) == 1:
        vzz = angt[0][0]
        vaa = angt[0][1]
        szz = angt[0][2]
        saa = angt[0][3]
      else:
        # sz, saz generated by user-specified dates and times, lat, lon if that option given
        # AND: need to check angles that we are absolute i.e. no -ves
        vzz = angt[a][0]
        vaa = angt[a][1]
        szz = angt[a][2]
        saa = angt[a][3]

      VLAB.logger.info(str(vzz))
      umu = math.cos(vzz)
      angstr = str(vzz) + '_' + str(vaa) + '_' + str(szz) + '_' + str(saa)
      
      libradtran_ip = VLAB.fPath(LIBRAT.SDIR, 'ip.' + q['root'] + '.' + q['wbfile'] + '_' + angstr)
      libradtran_op = VLAB.fPath(LIBRAT.SDIR, 'op.' + q['root'] + '.' + q['wbfile'] + '_' + angstr)

      if q['v']:
        sys.stderr.write('%s: doing ip file %s\n'%(sys.argv[0],libradtran_ip))
        #sort out zen/azimuth angles i.e. if vz is -ve
      if vzz < 0:
        # doesn't matter as we use umu from above anyway but ...
        vzz *= -1.
      if vaa < 0:
        vaa = 180. - vaa
      if szz < 0:
        szz *= -1.
      if saa < 0:
        saa = 180. - saa

      # add view angles
      def valueorlisttostring(val):
        try:
          return ' '.join(map(str, val))
        except TypeError:
          return str(val)
        
      r['umu'] = valueorlisttostring(umu)
      r['phi'] = valueorlisttostring(vaa)

      # add sz angles
      if q['lat'] and q['lon'] and q['time']:
        r['latitude'] = q['lat']
        r['longitude'] = q['lon']
        r['time'] = q['time']
      else:
        # print out sun angles
        r['sza'] = valueorlisttostring(szz)
        r['phi0'] = valueorlisttostring(saa)

      # write out wavelengths i.e. min and max. Step is determined by step in solar file i.e. 1nm default
      r['wavelength'] = str(int(min(wb))) + ' ' + str(int(max(wb)))
      # run the libradtran command
      r['infile']  = libradtran_ip
      r['outfile'] = libradtran_op
      for k,v in r.items():
        VLAB.logger.info("  [%s] => %s" % (k, v))
      VLAB.doLibradtran(r)

      # now check if a single angle (i.e. angt.size == 4 ) and if so, break out of loop
      if 4 == len(angt) * len(angt[0]) or (q['lat'] and q['lon'] and q['time']): break

      # finish angle loop & srt out collating results into LUT file at end if required
      # write header line i.e. vz va sz sa wb_min -> wb_max
      plotfilefp.write('# vz va sz sa %e %e %e\n' % (min(wb), max(wb), wbstep))
      # plotfilefp.flush()
	
      for f in [f for f in VLAB.listdir('%s/%s' % (LIBRAT.SDIR, q['opdir'])) if f.startswith('op.')]:
        vzz = f.split('_')[-4]
        vaa = f.split('_')[-3]
        szz = f.split('_')[-2]
        saa = f.split('_')[-1]
        d = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, f),transpose=True)
        plotfilefp.write('%s %s %s %s ' % (vzz,vaa,szz,saa))
        plotfilefp.write(' '.join(map(str, d[1])) + '\n')

      plotfilefp.flush()
      plotfilefp.close()

#############################################################################
class Librat_drivers:
  """ setup the scene"""
  def main(self, args):
    me=self.__class__.__name__ +'::'+VLAB.me()

    # initialize with defaults

    # sun angles for eg UK (lat 51' 0" 0' 0") 1 Jun 2013 10:30 am
    # sun_position -lat 51 -lon 0 -date 01:05:2013 -t 10:30:00
    # va, sz, sa = 0, 34, 141
    q = {
            'vz' : 0,
            'va' : 0,
            'sz' : 34,
            'sa' : 141,
      'nsamples' : 1000,
         'xdims' : (0, 300, 30),
         'ydims' : (0, 300, 30),
             'z' : 0,
       'angFile' : 'angles.dat',
      'lookFile' : 'look.dat',
    }

    # fill in: lookFile, angFile, xdums, ydims, z, sz, sa, nsamples...
    for a in args:
      if a == 'angles':
        q['angFile']  = args[a]
      elif a == 'look':
        q['lookFile'] = args[a]
      elif a == 'n':
        q['nsamples'] = args[a]
      else:
        q[a] = args[a]

    if 'msi' in q:
      # across whole swath
      # zmin, vzmax, vzstep = -22, 54, 4
      q['angles'] = (0, 0, 0, 0)
      q['angles'] = (q['vz'], q['va'], q['sz'], q['sa'])
    elif 'olci' in q:
      q['vzmin'], q['vzmax'], q['vzstep'] = -22, 54, 4
      q['vzz'] = VLAB.frange(q['vzmin'], q['vzmax']+q['vzstep'], q['vzstep'])
      q['angles'] = [[0. for col in range(4)] for row in range(len(q['vzz']))]
      q['angles'][:,0] = q['vzz']
      q['angles'][:,1] = q['vz']
      q['angles'][:,2] = q['sz']
      q['angles'][:,3] = q['sa']
    elif 'around' in q:
      q['vamin'], q['vamax'], q['vastep'] =  0, 360, 10
      q['vaa'] = VLAB.frange(q['vamin'], q['vamax']+q['vastep'], q['vastep'])
      q['vz'], q['sz'], q['sa'] = 70., 30, 90
      q['angles'] = [[0. for col in range(4)] for row in range(len(q['vaa']))]
      q['angles'][:,0] = q['vz']
      q['angles'][:,1] = q['vaa']
      q['angles'][:,2] = q['sz']
      q['angles'][:,3] = q['sa']
    elif 'multiple' in q:
      # this one for multiple vz, va angles
      q['vzmin'], q['vzmax'], q['vzstep'] = -70,  70, 10
      q['vamin'], q['vamax'], q['vastep'] =   0, 340, 20
      q['vzz'] = VLAB.frange(q['vzmin'], q['vzmax']+q['vzstep'], q['vzstep'])
      q['vaa'] = VLAB.frange(q['vamin'], q['vamax']+q['vastep'], q['vastep'])

      q['angles'] = [[0. for col in range(4)] for row in range(len(q['vzz'])*len(q['vzz']))]

      for n, va in enumerate(q['vaa']):
        q['angles'][n*len(q['vzz']):(n+1)*len(q['vzz']),0] = q['vzz']
        q['angles'][n*len(q['vzz']):(n+1)*len(q['vzz']),1] = q['vz']
        q['angles'][n*len(q['vzz']):(n+1)*len(q['vzz']),2] = q['sz']
        q['angles'][n*len(q['vzz']):(n+1)*len(q['vzz']),3] = q['sa']
    elif 'random' in q:
      # random n samples of COSINE-WEIGHTED view and illum angles
      # In order that we get the required number, nsamples,  AND can disregard all those with
      # vz > 70 and vz < -70 (not valid for RPV) do twice as many as you need, discard
      # all those outside the angle range and then take the first n.
      nn = q['nsamples']*2

      if VLAB.CONF_RND_REPRODUCE:
        rState = VLAB.rndInit(17)
      else:
        rState = VLAB.rndInit()

      # view angles first
      u1    = [VLAB.rndNextFloat(rState) for i in range(nn)]
      r     = [math.sqrt(i)        for x in u1]
      theta = [2.*math.pi*x for x in [VLAB.rndNextFloat(rState) for i in range(nn)]]


      x  = [r[i] * math.cos(theta[i]) for i in range(nn)]
      y  = [r[i] * math.sin(theta[i]) for i in range(nn)]
      z  = [math.sqrt(1 - u1[i])      for i in range(nn)]

      q['vz'] = [math.acos(z[i])      for i in range(nn)]
      q['va'] = [math.atan(y[i]/x[i]) for i in range(nn)]

      # if x -ve then vz -ve
      for i in range(nn):
        if (x[i]<0):
          q['vz'][i] = q['vz'][i] * -1.

      # sun angles
      u1    = [VLAB.rndNextFloat(rState) for i in range(nn)]
      r     = [math.sqrt(i)              for x in u1]
      theta = [2.*math.pi*x for x in [VLAB.rndNextFloat(rState) for i in range(nn)]]

      x  = [r[i] * math.cos(theta[i]) for i in range(nn)]
      y  = [r[i] * math.sin(theta[i]) for i in range(nn)]
      z  = [math.sqrt(1 - u1[i])      for i in range(nn)]

      #np.savetxt('rpv.angles.test.plot',np.transpose([x,y,z]),fmt='%.6f')
      q['sz'] = [math.acos(z[i])      for i in range(nn)]
      q['sa'] = [math.atan(y[i]/x[i]) for i in range(nn)]

      # if x -ve then sz -ve
      for i in range(nn):
        if (x[i]<0):
          q['sz'][i] = q['sz'][i] * -1.

      vzdeg = [VLAB.r2d(q['vz'][i]) for i in range(nn)]
      vadeg = [VLAB.r2d(q['va'][i]) for i in range(nn)]
      szdeg = [VLAB.r2d(q['sz'][i]) for i in range(nn)]
      sadeg = [VLAB.r2d(q['sa'][i]) for i in range(nn)]
      # zip does a transpose
      angles = zip(vzdeg, vadeg, szdeg, sadeg)

      angles2 = list()
      for row in angles:
        if (row[0] > -70.) and (row[0] <= 70.) and (row[2] > -70.) and (row[2] <= 70.):
          angles2.append(row)

      # now take the first n
      q['angles'] = angles2[0:q['nsamples']]

      # x = np.sin(np.deg2rad(angles[:,0])) * np.cos(np.deg2rad(angles[:,1]))
      # y = np.sin(np.deg2rad(angles[:,0])) * np.sin(np.deg2rad(angles[:,1]))
      # z = np.cos(np.deg2rad(angles[:,0]))

      #np.savetxt('rpv.angles.test',angles,fmt='%.6f')
      #np.savetxt('rpv.angles.test.plot',np.transpose([x,y,z]),fmt='%.6f')
    else:
      # default
      q['vzmin'], q['vzmax'], q['vzstep'] = -70, 70, 5
      q['vzz'] = VLAB.frange(q['vzmin'],q['vzmax']+q['vzstep'],q['vzstep'])
      q['angles'] = [[0. for col in range(4)] for row in range(len(q['vzz']))]
      q['angles'][:,0] = q['vzz']
      q['angles'][:,1] = q['va']
      q['angles'][:,2] = q['sz']
      q['angles'][:,3] = q['sa']

    VLAB.savetxt(q['angFile'], q['angles'], fmt='%.1f')

    if 'lookFile' in q:
      x = VLAB.frange(q['xdims'][0],q['xdims'][1],q['xdims'][2])
      y = VLAB.frange(q['ydims'][0],q['ydims'][1],q['ydims'][2])
      looktrans = [[0. for col in range(len(x)*len(y))] for row in range(3)]
      for n, yy in enumerate(y):
        looktrans[0][n*len(x):(n+1)*len(x)] = x
        looktrans[1][n*len(x):(n+1)*len(x)] = [yy for i in range(len(y))]
        looktrans[2][n*len(x):(n+1)*len(x)] = [0. for i in range(len(x))]
      # zip does a transpose
      look = zip(*looktrans)
      VLAB.savetxt(q['lookFile'],look,fmt='%.1f')

#############################################################################
class Librat_plot:
  """plot results"""
  def main(self, args):
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('======> %s' % me )
    for a in args:
      VLAB.logger.info("%s -> %s" % (a, args[a]))

    q = {
         'spec' : 'SPECTRAL_TEST/result.laegeren.obj_vz_0.0_va_0.0_sz_34.0_sa_141.0_xyz_150.0_150.0_700.0_wb_wb.full_spectrum.dat.direct',
       'wbfile' : 'wb.full_spectrum.dat',
         'root' : 'OLCI_brdf.scene/result.laegeren.obj',
         'spec' : False,
         'brdf' : True
    }

    for a in args:
      if a == 'angles':
        q['angfile'] = args[a]
      else:
        q[a] = args[a]

    if q['spec']:
      self.spec_plot(q['spec'], q['wbfile'])

    if q['brdf']:
      self.brdf_plot(q['root'], q['angfile'], q['wbfile'])

  def spec_plot(self,spec,wbspec):
    wbspec_contents = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, wbspec), transpose=True)

    wb = wbspec_contents[1];
    op1 = spec + '.plot.png'

    data = wbspec_contents
    refl = data[1:,].sum(axis=1)
    # TODO Implement rest of function
    raise Exception("spec_plot()")
  def brdf_plot(self,root,angfile,wbfile):
    opdat = LIBRAT.SDIR + '/' + root + '.brdf.dat'
    opplot = LIBRAT.SDIR + '/' + root + '.brdf.png'
    ang = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, angfile), transpose=True)
    wb = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, wbfile), transpose=True)[1]

    result = [[0. for i in range(len(wb))] for i in range(len(ang[0]))]

    ff = [root.split('/')[0] + '/' + f for f in VLAB.listdir(root.split('/')[0])
          if f.endswith('.direct') and f.startswith(root.split('/')[1])]
    for f in ff:
      fsplit = f.split('_')
      vz = f.split('_')[fsplit.index('vz') + 1]
      va = f.split('_')[fsplit.index('va') + 1]
      sz = f.split('_')[fsplit.index('sz') + 1]
      sa = f.split('_')[fsplit.index('sa') + 1]

      data = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, f), transpose=True)
      refl = [reduce(lambda x, y : x + y, i) for i in data[1:]]
      result[zip(*VLAB.awhere(
        VLAB.treemap(lambda x : x == float(vz), ang)))[1][0]] = refl
      # TODO test and cleanup (brdf_plot)
      # coords = VLAB.treemap(lambda x : x == float(vz), ang[0])
      # coords = map(lambda x, y : x & y, coords,
      #              VLAB.treemap(lambda x : x == float(va), ang[1]))
      # coords = map(lambda x, y : x & y, coords,
      #              VLAB.treemap(lambda x : x == float(sz), ang[2]))
      # coords = map(lambda x, y : x & y, coords,
      #              VLAB.treemap(lambda x : x == float(sa), ang[3]))
      coords = VLAB.makemaskeq(ang[0], float(vz))
      coords = VLAB.maskand(coords, VLAB.makemaskeq(ang[1], float(va)))
      coords = VLAB.maskand(coords, VLAB.makemaskeq(ang[2], float(sz)))
      coords = VLAB.maskand(coords, VLAB.makemaskeq(ang[3], float(sa)))
      coords = VLAB.awhere(coords)
      result[coords[0][0]] = refl
    dataset = VLAB.make_dataset()
    for b in [3, 7]:
      VLAB.plot(dataset, ang[0], [i[b] for i in result],
           "waveband: %.1f" % wb [b])
    chart = VLAB.make_chart("plot.py out", "view zenith angle (deg.)",
                            u"\u03A1", dataset)
    sys.stderr.write('%s: plotting brdf to: %s\n'%(sys.argv[0], opplot))
    sys.stderr.write('%s: saving brdf data to: %s\n'%(sys.argv[0], opdat))
    outdata = [[0. for i in xrange(len(result[0]) + len(ang))]
               for j in xrange(len(result))]
    VLAB.replacerectinarray(outdata, zip(*ang), 0, 0, len(result), len(ang))
    # FIXME: assign result instead of zip(*ang) (see plot.py:88)
    VLAB.replacerectinarray(outdata, result, 0, len(ang), len(result), len(result[0]) + len(ang))
    VLAB.savetxt(opdat, outdata, fmt='%.4f')
    VLAB.save_chart(chart, opplot)

class Librat_rpv_invert:
  """rpv inversion for integration with libmodtran"""
  def main(self, args):
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('======> %s' % me)
    for a in args:
      VLAB.logger.info("%s -> %s" % (a, args[a]))

    # defaults
    q = {
	'dataf' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat',
	'wbfile' : 'wb.MSI.dat',
	'wbNum' : 3, # 665 nm in this case
	'verbose' : 1,
	'plot' : 1,
	'show' : 0
    }

    for a in args:
      if a == 'wb':
        q['wbfile'] = args[a]
      else:
        q[a] = args[a]

    wb   = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, q['wbfile']), transpose=True)[1]
    data = VLAB.valuesfromfile('%s/%s' % (LIBRAT.SDIR, q['dataf']),  transpose=True)

    # check shape of 2 data files i.e. that there are same no. of wbs on each line of datafile ( + 4 angles)
    if len(wb) != len(data) - 4:
      sys.stderr.write('%s: no of wavebands different in brdf file %s and wb file %s\n'%(sys.argv[0],q['dataf'],q['wbfile']))
      raise Exception('%s: no of wavebands different in brdf file %s and wb file %s\n'%(sys.argv[0],q['dataf'],q['wbfile']))

    rho0, k, bigtet, rhoc = 0.03, 1.2, 0.1, 0.2

    if q['three']:
      params = [rho0, k, bigtet]
    else:
      params = [rho0, k, bigtet, rhoc]

    if args['paramfile']:
      opdat = args['paramfile']
    else:
      opdat = q['dataf'] + '.params.dat'

    if q['verbose']: sys.stderr.write('%s: saving params to %s\n'%(sys.argv[0], opdat))

    # create the file if it doesn't exist
    dfp = VLAB.openFileIfNotExists('%s/%s' % (LIBRAT.SDIR, opdat))
    if dfp != None:
      dfp.close()

    # open the previously created file
    opfp = open('%s/%s' % (LIBRAT.SDIR,  opdat), 'w')

    if q['three']:
      opfp.write('# wb rho0 k bigtet\n')
    else:
      opfp.write('# wb rho0 k bigtet rhoc\n')

    ymin, ymax = (0, 0.25)
    xmin, xmax = (-75., 75)

    for wbNum, band in enumerate(wb):
      if q['verbose']: sys.stderr.write('%s: doing band %i (%f)\n'%(sys.argv[0], wbNum, band))

      invdata = [[0. for i in xrange(len(data))] for i in xrange(5)]
      invdata[0:4] = data[0:4]
      invdata[4] = data[4 + q['wbNum']]

      # invdata = [[0. for i in xrange(15)] for i in xrange(5)]
      # invdata[0:4] = [i[0:15] for i in data[0:4]]
      # invdata[4] = data[4 + q['wbNum']][0:15]

      porig = params
      p_est = params
      p_est = VLAB.Minimize_minimize(self.obj, params, args=(invdata,))[0]
      if q['three']:
        p_est = VLAB.min_l_bfgs_b(self.obj, p_est, args=(invdata,), bounds=((0., None), (0., None), (None, None)))
      else:
        p_est = VLAB.min_l_bfgs_b(self.obj, p_est, args=(invdata,), bounds=((0., None), (0., None),(None, None), (None, None)))
      r = self.rpv(p_est, invdata)
      rmse = math.sqrt(reduce(lambda x, y : x + y, map(lambda x : x ** 2, VLAB.suba(r, invdata[4]))))
      if q['three']:
        opfp.write('%.1f %.8f %.8f %.8f\n' % (band, p_est[0], p_est[1], p_est[2]))
      else:
        opfp.write('%.1f %.8f %.8f %.8f %.8f\n' % (band, p_est[0], p_est[1], p_est[2], p_est[3]))
      if q['plot']:
        if q['plotfile']:
          opplot = VLAB.path.join(LIBRAT.SDIR, q['plotfile'] + '.inv.wb.' + str(wbNum) + '.png')
        else:
          opplot = VLAB.path.join(LIBRAT.SDIR, q['dataf'] + '.inv.wb.' + str(wbNum) + '.png')

        if q['verbose']: sys.stderr.write('%s: plotting to %s\n' % (sys.argv[0], opplot))

        dataset = VLAB.make_dataset()
        VLAB.plot(dataset, invdata[0], invdata[4], "original")
        VLAB.plot(dataset, invdata[0], r, "inverted")
        chart = VLAB.make_chart("", "vza (deg)", u"\u03A1", dataset)
        dataset = None
        VLAB.save_chart(chart, opplot)
        chart = None
        # TODO test plot code
    opfp.close()
  def obj(self, p, x):
    fwd = self.rpv(p, x)
    obs = x[4]
    sse = reduce(lambda x, y : x + y,
                 map(lambda x : x ** 2,
                     VLAB.suba(obs, fwd)))
    return sse
  def rpv(self, params, data):
    if len(params) == 4:
      rho0, k, bigtet, rhoc = params
    else:
      rhoc = 1
      rho0, k, bigtet = params
    cosv = VLAB.cosa(data[0])
    coss = VLAB.fabsa(VLAB.cosa(data[2]))
    cosv = VLAB.replaceinarray(cosv, lambda x : x == 0, 1e-20)
    coss = VLAB.replaceinarray(coss, lambda x : x == 0, 1e-20)
    sins = VLAB.sqrta(map(lambda x : 1. - x ** 2, coss))
    sinv = VLAB.sqrta(map(lambda x : 1. - x ** 2, cosv))
    relphi = VLAB.suba(map(VLAB.d2r, data[1]),
                       map(VLAB.d2r, data[3]))
    relphi = map(lambda x : {True : 2 * math.pi - x, False : x}
                 [x > math.pi], relphi)
    cosp = map(lambda x : -1. * x, VLAB.cosa(relphi))
    tans = VLAB.diva(sins, coss)
    tanv = VLAB.diva(sinv, cosv)
    csmllg = VLAB.adda(VLAB.mula(coss, cosv),
                       VLAB.mula(VLAB.mula(sins, sinv), cosp))
    bigg = VLAB.sqrta(
      VLAB.suba(
        VLAB.adda(map(lambda x : x ** 2, tans), map(lambda x : x ** 2, tanv)),
        VLAB.mula(VLAB.mula(map(lambda x : x * 2.0, tans), tanv), cosp)))
    bgthsq = bigtet ** 2
    expon = k - 1.0
    if expon != 0.0:
      f1 = VLAB.mula(VLAB.powa(VLAB.mula(coss, cosv), expon),
                     VLAB.powa(VLAB.adda(coss, cosv), expon))
    else:
      fi = [1. for i in coss]
    denom = VLAB.powa(map(lambda x : 1. + bgthsq + 2. * bigtet * x, csmllg),
                      1.5)
    f2 = list(denom)
    f2 = map(lambda x : {True : (lambda : (1.0 - bgthsq) * 1e20),
                         False : (lambda : (1.0 - bgthsq) / x)}[x == 0.](), relphi)
    f3 = map(lambda x : 1.0 + (1 - rhoc) / (1. + x), bigg)
    return map(lambda x : rho0 * x, VLAB.mula(VLAB.mula(f1, f2), f3))

#############################################################################

class LIBRAT:
  if VLAB.osName().startswith('Windows'):
    SDIR=VLAB.expandEnv('%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_scenes')
  else:
    SDIR=VLAB.expandEnv('$HOME/.beam/beam-vlab/auxdata/librat_scenes')

  """Integration glue for calling external LIBRAT programs"""
  def __init__(self):
    me=self.__class__.__name__ +'::'+VLAB.me()
    VLAB.logger.info('%s: constructor completed...' % me)

  def doProcessingTests(self, pm, args):
 
    """do processing tests for LIBRAT processor"""
    me=self.__class__.__name__ +'::'+VLAB.me()

    if (pm != None):
      pm.beginTask("Computing...", 10)
    # ensure at least 1 second to ensure progress popup feedback
    time.sleep(1)

    cmd = {
    'linux' : {
      'cwd'     : '$HOME/.beam/beam-vlab/auxdata/librat_lin64/src/start',
      'exe'     : '$HOME/.beam/beam-vlab/auxdata/librat_lin64/src/start/start',
      'cmdline' : [
        '-sensor_wavebands', 'wavebands.dat', '-m', '100',
        '-sun_position', '0', '0', '10', 'test.obj'],
      'stdin'   : '$HOME/.beam/beam-vlab/auxdata/librat_lin64/src/start/starttest.ip',
      'stdout'  : None,
      'stderr'  : None,
      'env'     : {
        'LD_LIBRARY_PATH' : '$HOME/.beam/beam-vlab/auxdata/librat_lin64/src/lib/',
        'BPMS'  : '$HOME/.beam/beam-vlab/auxdata/librat_lin64/'
      }},
    'windows'   : {
      'cwd'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_windows/src/start',
      'exe'     : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_windows/src/start/ratstart.exe',
      'cmdline' : [
        '-sensor_wavebands', 'wavebands.dat', '-m', '100',
        '-sun_position', '0', '0', '10', 'test.obj'],
      'stdin'   : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_windows/src/start/starttest.ip',
      'stdout'  : None,
      'stderr'  : None,
      'env'     : {
        'BPMS'  : '%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_windows'
     }}
    }
    VLAB.doExec(cmd)

    # 
    # see https://github.com/netceteragroup/esa-beam/blob/master/beam-3dveglab-vlab/src/main/scenes/librat_scenes/libradtran.README
    #
    drivers      = Librat_drivers()
    dobrdf       = Librat_dobrdf()
    plot         = Librat_plot()
    rpv_invert   = Librat_rpv_invert()
    dolibradtran = Librat_dolibradtran()

    args = {
       'random' : True,
       'n'      : 1000,
       'angles' : 'angles.rpv.2.dat',
    }
    # drivers.main(args)

    # RAMI test case
    args = {
            'v' : True,
         'nice' : 19,
          'obj' : 'HET01_DIS_UNI_NIR_20.obj',
         'hips' : True,
           'wb' : 'wb.MSI.dat',
        'ideal' : (300., 300.),
         'look' : (150., 150., 710.),
          'rpp' : 4,
      'npixels' : 10000,
         'boom' : 786000,
       'angles' : 'angles.rpv.2.dat',
        'opdir' : 'rpv.rami'
    }
    # dobrdf.main(args)

    # LAEGEREN test case
    args = {
            'v' : True,
         'nice' : 19,
          'obj' : 'laegeren.obj.lai.1',
         'hips' : True,
           'wb' : 'wb.MSI.dat',
        'ideal' : (300., 300.),
         'look' : (150., 150., 710.),
          'rpp' : 4,
      'npixels' : 10000,
         'boom' : 786000 ,
       'angles' : 'angles.rpv.2.dat',
        'opdir' : 'rpv.laegeren'
    }
    # dobrdf.main(args)

    # FULL SPECTRUM RAMI test case
    args = {
            'v' : True,
         'nice' : 19,
          'obj' : 'HET01_DIS_UNI_NIR_20.obj',
           'wb' : 'wb.full_spectrum.1nm.dat',
        'ideal' : (80., 80.),
         'look' : (0., 0., 0.),
          'rpp' : 4, 
      'npixels' : 10000,
         'boom' : 786000, 
       'angles' : 'angle.rpv.cosDOM.dat',
        'opdir' : 'dart.rpv.rami'
    }
    # dobrdf.main(args)

    # FULL SPECTRUM LAEGEREN test case
    args = {
            'v' : True,
         'nice' : 19,
          'obj' : 'laegeren.obj.lai.1',
           'wb' : 'wb.full_spectrum.1nm.dat',
        'ideal' : (300., 300),
         'look' : (150., 150., 710),
          'rpp' : 4,
      'npixels' : 10000,
         'boom' : 786000,
       'angles' : 'angle.rpv.cosDOM.dat',
        'opdir' : 'dart.rpv.laegeren'
    }
    # dobrdf.main(args)

    args = {
         'brdf' : True,
       'angles' : 'angles.rpv.2.dat',
       'wbfile' : 'wb.MSI.dat',
        'bands' : (4, 7),
         'root' : 'rpv.rami/result.HET01_DIS_UNI_NIR_20.obj'
    }
    # plot.main(args)

    # RAMI
    args = {
         'brdf' : True,
       'angles' : 'angles.rpv.2.dat', 
       'wbfile' : 'wb.MSI.dat',
         'root' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj',
        'bands' : (4, 7),
            'v' : True
    }
    # plot.main(args)

    # LAEGEREN
    args = {
         'brdf' : True, 
       'angles' : 'angles.rpv.2.dat',
       'wbfile' : 'wb.MSI.dat',
         'root' : 'rpv.laegeren/result.laegeren.obj.lai.1',
        'bands' : (4, 7),
            'v' : True
    }
    # plot.main(args)

    # RAMI
    args = {
         'brdf' : True,
       'wbfile' : 'wb.full_spectrum.1nm.dat',
       'angles' : 'angle.rpv.cosDOM.dat',
         'root' : 'dart.rpv.rami/result.HET01_DIS_UNI_NIR_20.obj',
        'bands' : (250, 450),
            'v' : True
    }
    # plot.main(args)

    # RAMI
    args = {
        'three' : True,
            'v' : True,
         'plot' : True,
        'dataf' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat',
    'paramfile' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat.3params.dat',
     'plotfile' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj.brdf.3params',
    }
    # rpv_invert.main(args)

    # LAEGEREN
    args = {
        'three' : True,
            'v' : True,
         'plot' : True,
        'dataf' : 'rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat',
    'paramfile' : 'rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat.3params.dat',
     'plotfile' : 'rpv.laegeren/result.laegeren.obj.lai.1.brdf.3params'
    }
    # rpv_invert.main(args)

    args = {
        'three' : True,
            'v' : True,
         'plot' : True,
        'dataf' : 'dart.rpv.rami/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat',
    'paramfile' : 'dart.rpv.rami/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat.3params.dat',
     'plotfile' : 'dart.rpv.rami/result.HET01_DIS_UNI_NIR_20.obj.brdf.3params',
    }
    # rpv_invert.main(args)

    args = {
         'three': True,
            'v' : True,
         'plot' : True,
        'dataf' : 'dart.rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat',
    'paramfile' : 'dart.rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat.3params.dat',
     'plotfile' : 'dart.rpv.laegeren/result.laegeren.obj.lai.1.brdf.3params'
    }
    # rpv_invert.main(args)

    args = {
         'v' : True,
     'opdir' : 'rami.TOA',
    'angles' : 'angles.MSI.dat',
    'wbfile' : 'wb.MSI.dat',
       'rpv' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat.3params.dat' ,
      'plot' : 'rami.TOA/rpv.rami.libradtran.dat.all'
    }
    # dolibradtran.main(args)

    args = {
         'v' : True,
     'opdir' : 'laegeren.TOA',
    'angles' : 'angles.MSI.dat',
    'wbfile' : 'wb.MSI.dat',
       'rpv' : 'rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat.3params.dat',
      'plot' : 'laegeren.TOA/rpv.laegeren.libradtran.dat.all'
    }
    # dolibradtran.main(args)

    args = {
         'v' : True,
    'angles' : 'angle.rpv.cosDOM.dat',
    'wbfile' : 'wb.full_spectrum.1nm.dat', 
     'opdir' : 'dart.rami.TOA',
       'rpv' : 'dart.rpv.rami/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat.3params.dat',
      'plot' : 'dart.rami.TOA/rpv.rami.libradtran.dat.all'
    }
    # dolibradtran.main(args)

    args = {
          'v': True,
    'angles' : 'angle.rpv.cosDOM.dat',
    'wbfile' : 'wb.full_spectrum.1nm.dat',
     'opdir' : 'dart.laegeren.TOA',
       'rpv' : 'dart.rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat.3params.dat', 
      'plot' : 'dart.laegeren.TOA/rpv.laegeren.libradtran.dat.all'
    }
    # dolibradtran.main(args)

    args = {
          'v' : True,
      'opdir' : 'rami.TOA',
        'rpv' : 'rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat.3params.dat',
       'plot' : 'rami.TOA/rpv.rami.libradtran.dat.all',
        'lat' : 50,
        'lon' : 0,
       'time' : '2013 0601 12 00 00'
    }
    # dolibradtran.main(args)

    args = {
          'v': True, 
     'opdir' : 'laegeren.TOA.date',
       'rpv' : 'rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat.3params.dat',
      'plot' : 'laegeren.TOA.date/rpv.laegeren.libradtran.dat.all',
       'lat' : 50, 
       'lon' : 0,
      'time' : '2013 06 01 12 00 00'
    }
    # dolibradtran.main(args)


    # DHP Simulation
    args = {
         'v' : True,
      'nice' : 19,
       'obj' : 'laegeren.obj.lai.1',
      'hips' : True,
        'wb' : 'wb.image.dat',
       'dhp' : True,
'samplingPattern' : 'circular',
  'lookFile' : 'dhp.locations.ondem.dat',
    'angles' : 'angles.dhp.dat',
       'fov' : False,
       'rpp' : 8,
   'npixels' : 4000000,
     'opdir' : 'DHP_TEST'
    }
    # dobrdf.main(args)
    VLAB.logger.info('%s: Done...' % me)

  def doProcessing(self, pm, args):
    """do processing for LIBRAT processor"""
    me=self.__class__.__name__ +'::'+VLAB.me()

    VLAB.logger.info('%s: doProcessing()' % (me))
    # defaults
    q = {
       'angles' : 'angle.rpv.cosDOM.dat',
        'bands' : (250, 450),
         'boom' : 786000, 
         'brdf' : True,
        'dataf' : 'dart.rpv.laegeren/result.laegeren.obj.lai.1.brdf.dat',
          'dhp' : True,
          'fov' : False,
         'hips' : False,
        'image' : True,
        'ideal' : (300., 300.),
          'lat' : 50, 
          'lon' : 0,
         'look' : (150., 150., 710.),
     'lookFile' : 'dhp.locations.ondem.dat',
            'n' : 1000,
         'nice' : 19,
      'npixels' : 10000,
          'obj' : 'DEFAULT_OBJ',
        'opdir' : 'DEFAULT_DIR',
    'paramfile' : 'dart.rpv.rami/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat.3params.dat',
         'plot' : 'dart.rami.TOA/rpv.rami.libradtran.dat.all',
     'plotfile' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj.brdf.3params',
         'plot' : 'rami.TOA/rpv.rami.libradtran.dat.all',
       'random' : True,
         'root' : 'rpv.rami/result.HET01_DIS_UNI_NIR_20.obj',
          'rpp' : 1, 
          'rpv' : 'dart.rpv.rami/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat.3params.dat',
          'rpv' : 'rpv.rami.2/result.HET01_DIS_UNI_NIR_20.obj.brdf.dat.3params.dat' ,
'samplingPattern' : 'circular',
        'three' : True,
         'time' : '2013 06 01 12 00 00',
            'v' : True,
       'wbfile' : 'wb.full_spectrum.1nm.dat',
           'wb' : 'wb.full_spectrum.1nm.dat'
    }

    VLAB.logger.info('%s: overwriting defaults' % (me))

    # overwrite defaults
    for a in args:
      if a == 'Sensor':
        q['sensor'] = args[a]
        if args[a] == VLAB.K_SENTINEL2:
          q['wb'] = 'wb.MSI.dat'
          q['wbfile'] = 'wb.MSI.dat'
          q['angfile'] = 'angles.MSI.dat'
        elif args[a] == VLAB.K_SENTINEL3:
          q['wb'] = 'wb.OLCI.dat'
          q['wbfile'] = 'wb.OLCI.dat'
          q['angfile'] = 'angles.OLCI.dat'
        elif args[a] == VLAB.K_MODIS:
          q['wb'] = 'wb.MODIS.dat'
          q['wbfile'] = 'wb.MODIS.dat'
          q['angfile'] = 'angles.MODIS.dat'
        elif args[a] == VLAB.K_MERIS:
          q['wb'] = 'wb.MERIS.dat'
          q['wbfile'] = 'wb.MERIS.dat'
          q['angfile'] = 'angles.MERIS.dat'
        elif args[a] == VLAB.K_LANDSAT_OLI:
          q['wb'] = 'wb.LANDSAT.OLI.dat'
          q['wbfile'] = 'wb.LANDSAT.OLI.dat'
          q['angfile'] = 'angles.LANDSAT.dat'
        elif args[a] == VLAB.K_LANDSAT_ETM:
          q['wb'] = 'wb.LANDSAT.ETM.dat'
          q['wbfile'] = 'wb.LANDSAT.ETM.dat'
          q['angfile'] = 'angles.LANDSAT.dat'
        else:
          q['wb'] = 'wb.full_spectrum.1nm.dat'
          q['wbfile'] = 'wb.full_spectrum.1nm.dat'
          q['angfile'] = 'angles.rami.dat'
      elif a == 'OutputDirectory':
        q['opdir'] = args[a]
      elif a == '3dScene':
        if args[a] == VLAB.K_RAMI:
          q['obj'] = 'HET01_DIS_UNI_NIR_20.obj'
          q['lat'] = 0
          q['lon'] = 0
        elif args[a] == VLAB.K_LAEGERN:
          q['obj'] = 'laegeren.obj'
          q['lat'] = '47.4817'
          q['lon'] = '-8.3947'
        elif args[a] == VLAB.K_THARANDT:
          q['obj'] = 'HET01_DIS_UNI_NIR_20.obj' # FIXME: define when file is available
          q['lat'] = '50.9833'
          q['lon'] = '-13.5808'
        elif args[a] == VLAB.K_USER_DEFINED:
          q['obj'] = 'UserDefined.obj'
          # What about q['lat'] and q['lon']?
      elif a == 'Bands':
        q['bands'] = [int(i) for i in tuple(args[a].split(", "))]
      elif a == 'DHP_ImageFile':
        q['dhp'] = args[a] == 'Yes'
      elif a == 'ImageFile':
        q['image'] = args[a] == 'Yes'
      elif a == 'AsciiFile':
        q['ascii'] = args[a] == 'Yes'
      elif a == 'ViewingAzimuth':
        q['va'] = args[a]
      elif a == 'ViewingZenith':
        q['vz'] = args[a]
      elif a == 'DHP_Zenith':
        q['dhp_vz'] = args[a]
      elif a == 'DHP_Azimuth':
        q['dhp_va'] = args[a]
      elif a == 'IlluminationAzimuth':
        q['sz'] = args[a]
      elif a == 'IlluminationZenith':
        q['sa'] = args[a]
      else:
        q[a] = args[a]

    q['dataf']     = '%s/result.%s.brdf.dat' %(q['opdir'], q['obj'])
    q['paramfile'] = '%s/result.%s.brdf.dat.3params.dat' %(q['opdir'], q['obj'])
    q['plot']      = '%s/rpv.%s.dat.all' %(q['opdir'],q['obj'])
    q['plotfile']  = '%s/result.%s.brdf.3params' %(q['opdir'],q['obj'])
    q['root']      = '%s/result.%s' %(q['opdir'],q['obj'])
    q['rpv']       = '%s/result.%s.brdf.dat.3params.dat' %(q['opdir'],q['obj'])

    if 'dhp' in q:
      q['vz'] = q['dhp_vz']
      q['va'] = q['dhp_va']

    for a in ['dataf', 'paramfile','plot','plotfile','root', 'rpv']:
      VLAB.logger.info ('%s: q["%s"] is %s' % (me, a, q[a]))

    if VLAB.osName().startswith('Windows'):
      fullobjpath = '%s/%s' % (VLAB.expandEnv('%HOMEDRIVE%%HOMEPATH%/.beam/beam-vlab/auxdata/librat_scenes'), args['OutputDirectory'])
    else:
      fullobjpath = '%s/%s' % (VLAB.expandEnv('$HOME/.beam/beam-vlab/auxdata/librat_scenes'), args['OutputDirectory'])
    
    VLAB.logger.info ('%s: ensuring "%s" exists' %(me, fullobjpath))
    if not VLAB.path.exists(fullobjpath):
      VLAB.mkDirPath(fullobjpath)

    VLAB.logger.info('%s: instantiating objects' % (me))
    drivers      = Librat_drivers()
    dobrdf       = Librat_dobrdf()
    plot         = Librat_plot()
    rpv_invert   = Librat_rpv_invert()
    dolibradtran = Librat_dolibradtran()

    if (pm != None):
      pm.beginTask("Computing BRF...", 10)
    # ensure at least 1 second to ensure progress popup feedback
    time.sleep(1)

    # not needed because we are using the dart angles
    # drivers.main()

    VLAB.logger.info('%s: calling dobrdf.main()' % (me))
    dobrdf.main(q)
    VLAB.logger.info('%s: calling plot.main()' % (me))
    plot.main(q)
    VLAB.logger.info('%s: calling rpv.main()' % (me))
    rpv_invert.main(q)
    VLAB.logger.info('%s: calling dolibradtran.main()' % (me))
    dolibradtran.main(q)

    VLAB.logger.info('%s: Done...' % me)

#### LIBRAT end ##############################################################
  
#### MAIN DISPATCH ####################################################

if not sys.platform.startswith('java'):
  VLAB.logger.info('%s: Mode: not jython and not BEAM' % VLAB.me)
  vlab = VLAB()
  vlab.selftests()
else:
  from java.lang import System
  beamVer = System.getProperty('beam.version')
  if beamVer == None:
    if System.getProperty('vlab.fakebeam') != None:
      VLAB.logger.info('%s: Mode: jython and fake BEAM' % VLAB.me)
      vlab = VLAB()
      vlab.fakebeam()
    else:
      VLAB.logger.info('%s: Mode: jython and not BEAM' % VLAB.me)
      vlab = VLAB()
      vlab.selftests()
  else:
    VLAB.logger.info('%s: Mode: jython and BEAM v%s' % (VLAB.me, beamVer))

### BEAM-only code ####################################################

    #
    # BEAM-specific code from here to end of file...
    #

    from java              import awt
    from java              import io
    from java              import lang
    from java.io           import BufferedReader
    from java.io           import BufferedWriter
    from java.io           import File
    from java.io           import FileInputStream
    from java.io           import FileOutputStream
    from java.io           import FileReader
    from java.io           import FileWriter
    from java.io           import InputStreamReader
    from java.io           import IOException
    from java.io           import OutputStreamWriter
    from java.lang         import IllegalArgumentException
    from java.lang         import Integer
    from java.lang         import ProcessBuilder
    from java.lang         import System
    from java.lang         import Thread
    from java.util         import ArrayList
    from java.util         import Arrays
    from java.util         import HashMap
    from java.util         import Map
    from java.util         import Vector
    from javax             import swing

    from com.bc.ceres.core                       import ProgressMonitor
    from com.netcetera.vlab                      import IVLabProcessor
    from com.netcetera.vlab                      import IVLabProcessorUi
    from com.netcetera.vlab                      import VLabProcessor
    from com.netcetera.vlab                      import VLabUi

    from org.esa.beam.dataio.dimap               import DimapProductConstants;
    from org.esa.beam.framework.dataio           import ProductSubsetDef
    from org.esa.beam.framework.dataio           import ProductWriter
    from org.esa.beam.framework.datamodel        import PixelPos
    from org.esa.beam.framework.datamodel        import Product
    from org.esa.beam.framework.datamodel        import RasterDataNode
    from org.esa.beam.framework.dataop.dem       import ElevationModelDescriptor
    from org.esa.beam.framework.dataop.dem       import ElevationModelRegistry
    from org.esa.beam.framework.dataop.maptransf import MapInfo
    from org.esa.beam.framework.dataop.resamp    import ResamplingFactory
    from org.esa.beam.framework.param            import Parameter
    from org.esa.beam.framework.param            import ParamProperties
    from org.esa.beam.framework.param            import ParamValidateException
    from org.esa.beam.framework.param.validators import NumberValidator
    from org.esa.beam.framework.processor        import DefaultRequestElementFactory
    from org.esa.beam.framework.processor        import Processor
    from org.esa.beam.framework.processor        import ProcessorConstants
    from org.esa.beam.framework.processor        import ProcessorException
    from org.esa.beam.framework.processor        import ProcessorUtils
    from org.esa.beam.framework.processor        import ProductRef
    from org.esa.beam.framework.processor        import Request
    from org.esa.beam.framework.processor        import RequestElementFactory
    from org.esa.beam.framework.processor        import RequestElementFactoryException
    from org.esa.beam.framework.processor.ui     import ProcessorUI
    from org.esa.beam.framework.ui               import GridBagUtils
    from org.esa.beam.framework.ui               import UIUtils
    from org.esa.beam.util                       import Guardian

    ##
    ## BEAM-defined Processor implementation
    ##
    class VLabImpl(IVLabProcessor):
      """Implements the BEAM Processor interface"""
  
      def __init__(self):
        me=self.__class__.__name__ +'::'+VLAB.me()
        VLAB.logger.info('%s: constructor completed...' % me)
  
      def getName(self):
        return VLAB.PROCESSOR_NAME
  
      def getSymbolicName(self):
        return VLAB.PROCESSOR_SNAME
  
      def getVersion(self):
        return VLAB.VERSION_STRING
  
      def getCopyrightInformation(self):
        return VLAB.COPYRIGHT_INFO
  
      def getUITitle(self):
        return VLAB.UI_TITLE
  
      def _getP(self, r, k):
        me=self.__class__.__name__ +'::'+VLAB.me()
        return r.getParameter(k).getValueAsText()
  
      def process(self, pm, req):
        VLAB.logger.info('inside process...')
        me=self.__class__.__name__ +'::'+VLAB.me()
        ProcessorUtils.setProcessorLoggingHandler(VLAB.DEFAULT_LOG_PREFIX, 
          req, self.getName(), self.getVersion(), self.getCopyrightInformation())
        #VLAB.log("Parameter list:")
        #for i in range(req.getNumParameters()):
        #  VLAB.log(req.getParameterAt(i).getName() + " = " + req.getParameterAt(i).getValueAsText())
  
        VLAB.logger.info('%s: %s' % (me, ProcessorConstants.LOG_MSG_START_REQUEST))
        pm.beginTask('Running %s...' % VLAB.PROCESSOR_NAME, 10)
        VLAB.logger.info('%s: after pm.beginTask' % me)
  
        # ensure at least 1 second to ensure progress popup feedback
        time.sleep(1)
        VLAB.logger.info('%s: after sleep' % me)
  
        processor = self._getP(req, VLAB.P_RTProcessor)
        VLAB.logger.info('%s: processor is <%s>' % (me, processor))
        if processor == VLAB.K_DART:
          rtProcessor = DART()
        elif processor == VLAB.K_LIBRAT:
          rtProcessor = LIBRAT()
        elif processor == VLAB.K_DUMMY:
          rtProcessor = DUMMY()
        else:
          raise RuntimeError('unknown processor: <%s>' % processor)
  
        pm.beginTask("Computing top of canopy BRF...", 10)
        params = {}
        for i in VLAB.plst:
          val = self._getP(req, i)
          VLAB.logger.info('%s: val for %s is %s' % (me, i, val))
          params[i] = val

        VLAB.logger.info('%s: calling doProcessing' % me)
        rtProcessor.doProcessing(pm, params)

        VLAB.logger.info('%s : %s' % (me, ProcessorConstants.LOG_MSG_FINISHED_REQUEST))
        pm.done()
  
    ##
    ## BEAM-defined UI Implementation
    ##
    class VLabUiImpl(IVLabProcessorUi):
      """Implements the BEAM ProcessorUI interface"""
      def __init__(self):
        self._reqElemFac     = VLabRequestElementFactory() 
        self._defaultFactory = DefaultRequestElementFactory.getInstance()
        self._requestFile    = File('')
        self.pmap            = {} 
  
        me=self.__class__.__name__ +'::'+VLAB.me()
        VLAB.logger.info('%s: constructor completed...' % me)
  
      def getGuiComponent(self):
        self._paramOutputProduct = self._reqElemFac.createDefaultOutputProductParameter()
        v = 2; h = 2;
        guiComponent = swing.JPanel(awt.BorderLayout())
        tabbedPane = swing.JTabbedPane()
        for tabGroups in VLAB.model:
          for tabName in tabGroups:
            tab = swing.JPanel()
            tab.layout = swing.BoxLayout(tab, swing.BoxLayout.PAGE_AXIS)
            tabbedPane.addTab(tabName, tab)
            for group in tabGroups[tabName]:
              tab.add(swing.JLabel(''))
              p = swing.JPanel()
              p.layout = awt.GridLayout(0, 4)
              p.layout.vgap = v; p.layout.hgap = h
              for groupName in group:
                p.border = swing.BorderFactory.createTitledBorder(groupName)
                for groupTuple in group[groupName]:
                  if len(groupTuple) == 4:
                    (lbl, nm, typ, vals) = groupTuple
                    if type(vals) == tuple:
                      dflt = vals[0]
                    else:
                      dflt = vals
                    props = self._defaultFactory.createStringParamProperties()
                    props.setLabel(lbl)
                    props.setDefaultValue(dflt)
                    self._reqElemFac.pMap[nm] = props
                    if type(vals) == tuple:
                      self.pmap[nm] = Parameter(nm, dflt)
                      (self.pmap[nm]).getProperties().setValueSet(vals)
                      (self.pmap[nm]).getProperties().setValueSetBound(True)
                      (self.pmap[nm]).getProperties().setDefaultValue(dflt)
                      (self.pmap[nm]).getProperties().setLabel(lbl)
                    else:
                      self.pmap[nm] = self._reqElemFac.createParameter(nm, dflt)
                    if type(typ) == tuple:
                      self.pmap[nm].getEditor().setEnabled(typ[1])
                    #p.add((self.pmap[nm]).getEditor().getLabelComponent())
                    p.add(swing.JLabel(lbl+':', swing.SwingConstants.RIGHT))
                    p.add((self.pmap[nm]).getEditor().getComponent())
                  else:
                    p.add(swing.JLabel(''))
                    p.add(swing.JLabel(''))
              tab.add(p)
            # hack
            for i in range(50):
              tab.add(swing.Box.createVerticalGlue())
        guiComponent.add(tabbedPane, awt.BorderLayout.NORTH)
        guiComponent.add(swing.JLabel(''), awt.BorderLayout.CENTER)
        VLabUi.setWindowSize(800, 800)
        return guiComponent
  
      def setRequests(self, requests):
        if (not requests.isEmpty()):
          for i in range(requests.size()):
            request = requests.elementAt(i)
            if (request == None):
              continue;
            if (str(VLAB.REQUEST_TYPE) == request.getType()):
              self._requestFile = request.getFile()
              outputProductAt = request.getOutputProductAt(0)
              if (outputProductAt != None):
                self._paramOutputProduct.setValueAsText(outputProductAt.getFilePath(), None)
              # update parameters
              for nm in VLAB.plst:
                (self.pmap[nm]).setValue(request.getParameter(nm).getValue())
  
    #          prefixParam = request.getParameter(ProcessorConstants.LOG_PREFIX_PARAM_NAME)
    #          if (prefixParam != None):
    #            self._logPrefixParameter.setValue(prefixParam.getValue(), None)
    #          logOutputParam = request.getParameter(ProcessorConstants.LOG_TO_OUTPUT_PARAM_NAME)
    #          if (logOutputParam != None):
    #            self._logToOutputParameter.setValue(logOutputParam.getValue(), None)
              break;
        else:
          self.setDefaultRequests()
  
      def setDefaultRequests(self):
        self.setDefaultRequest()
  
      def getRequests(self):
        requests = Vector()
        request = Request()
        request.setFile(self._requestFile)
        request.setType(VLAB.REQUEST_TYPE)
        outputFile = self._paramOutputProduct.getValueAsText()
        request.addOutputProduct(ProcessorUtils.createProductRef(outputFile, DimapProductConstants.DIMAP_FORMAT_NAME));
        for nm in VLAB.plst:
          request.addParameter(self._reqElemFac.createParameter(nm, (self.pmap[nm]).getValueAsText()))
        requests.add(request)
        return requests
  
      def setDefaultRequest(self):
        self._requestFile = None
        outputProductFile = self._paramOutputProduct.getValue()
        if (outputProductFile != None and outputProductFile.getParentFile() != None):
          parentFile = outputProductFile.getParentFile()
          self._paramOutputProduct.setValue(File(parentFile, VLAB.D_PRODNAME), None)
        else:
          self._paramOutputProduct.setDefaultValue()
        (self.pmap[VLAB.P_3dScene]).setDefaultValue()
  
    ##
    ## BEAM-defined mangement of processing parameters (aka request element)
    ##
    class VLabRequestElementFactory(RequestElementFactory):
      """Implements the BEAM RequestElementFactory interface"""
      def __init__(self):
        self._defaultFactory = DefaultRequestElementFactory.getInstance()
        self.pMap = {}
        me=self.__class__.__name__ +'::'+VLAB.me()
        VLAB.logger.info('%s : constructor completed...' % me)
        
      def getInstance(self):
        return VLabRequestElementFactory()
  
      def createParameter(self, name, value):
        Guardian.assertNotNullOrEmpty("name", name)
        try:
          param = self.createParamWithDefaultValueSet(name)
          if (value != None):
            param.setValueAsText(value, None)
        except IllegalArgumentException, e:
          VLAB.logger.info(e.getMessage)
          raise RequestElementFactoryException(e.getMessage())
        return param
        
      def createDefaultInputProductParameter(self):
        return self._defaultFactory.createDefaultInputProductParameter()
       
      def createDefaultLogPatternParameter(self, prefix):
        return self._defaultFactory.createDefaultLogPatternParameter(prefix)
        
      def createDefaultOutputProductParameter(self):
        defaultOutputProductParameter = self._defaultFactory.createDefaultOutputProductParameter()
        properties = defaultOutputProductParameter.getProperties()
        defaultValue = properties.getDefaultValue()
        if (isinstance(defaultValue, File)):
          properties.setDefaultValue(File(defaultValue, VLAB.D_PRODNAME))
        defaultOutputProductParameter.setDefaultValue()
        return defaultOutputProductParameter
        
      def createInputProductRef(self, fileN, fileFmt, typeId):
        try:
          return self._defaultFactory.createInputProductRef(fileN, fileFmt, typeId)
        except RequestElementFactoryException, e:
          raise RequestElementFactoryException(e.getMessage())
        
      def createLogToOutputParameter(self, value):
        try:
          return self._defaultFactory.createLogToOutputParameter(value)
        except ParamValidateException, e:
          raise ParamValidateException(e.getMessage())
                
      def createOutputProductRef(self, fileN, fileFmt, typeId):
        try:
          return self._defaultFactory.createOutputProductRef(fileN, fileFmt, typeId)
        except RequestElementFactoryException, e:
          raise RequestElementFactoryException(e.getMessage())
        
      def createParamWithDefaultValueSet(self, paramName):
        paramProps = self.getParamInfo(paramName)
        param = Parameter(paramName, paramProps.createCopy())
        param.setDefaultValue()
        return param
      
      def getParamInfo(self, parameterName):
        paramProps = self.pMap[parameterName]
        if (paramProps == None):
          if (parameterName.endsWith(VLAB.P_EXPRESSION)):
            VLAB.lof('to be implemented!!!')
          elif (parameterName.endsWith(VLAB.P_CONDITION)):
            VLAB.lof('to be implemented!!!')
          elif (parameterName.endsWith(VLAB.P_OUTPUT)):
            VLAB.lof('to be implemented!!!')
  
        if (paramProps == None):
          raise IllegalArgumentException("Invalid parameter name '" + parameterName + "'.")
        return paramProps
  
      def createStringParamProperties(self):
        return self._defaultFactory.createStringParamProperties()