Index: doc/guide/orientation/orientation.rst
===================================================================
--- doc/guide/orientation/orientation.rst	(revision 82592da3721b7e2098cf6346cf224814189fadb3)
+++ doc/guide/orientation/orientation.rst	(revision 5fb0634b45bba451ca824eecf40346b1ed5305a7)
@@ -4,16 +4,21 @@
 ==================
 
-With two dimensional small angle diffraction data SasView will calculate
+With two dimensional small angle diffraction data sasmodels will calculate
 scattering from oriented particles, applicable for example to shear flow
 or orientation in a magnetic field.
 
 In general we first need to define the reference orientation
-of the particles with respect to the incoming neutron or X-ray beam. This
-is done using three angles: $\theta$ and $\phi$ define the orientation of
-the axis of the particle, angle $\Psi$ is defined as the orientation of
-the major axis of the particle cross section with respect to its starting
-position along the beam direction. The figures below are for an elliptical
-cross section cylinder, but may be applied analogously to other shapes of
-particle.
+of the particle's $a$-$b$-$c$ axes with respect to the incoming
+neutron or X-ray beam. This is done using three angles: $\theta$ and $\phi$
+define the orientation of the $c$-axis of the particle, and angle $\Psi$ is
+defined as the orientation of the major axis of the particle cross section
+with respect to its starting position along the beam direction (or
+equivalently, as rotation about the $c$ axis). There is an unavoidable
+ambiguity when $c$ is aligned with $z$ in that $\phi$ and $\Psi$ both
+serve to rotate the particle about $c$, but this symmetry is destroyed
+when $\theta$ is not a multiple of 180.
+
+The figures below are for an elliptical cross section cylinder, but may
+be applied analogously to other shapes of particle.
 
 .. note::
@@ -29,8 +34,8 @@
 
     Definition of angles for oriented elliptical cylinder, where axis_ratio
-    b/a is shown >1, Note that rotation $\theta$, initially in the $x$-$z$
+    b/a is shown >1. Note that rotation $\theta$, initially in the $x$-$z$
     plane, is carried out first, then rotation $\phi$ about the $z$-axis,
     finally rotation $\Psi$ is around the axis of the cylinder. The neutron
-    or X-ray beam is along the $z$ axis.
+    or X-ray beam is along the $-z$ axis.
 
 .. figure::
@@ -40,57 +45,154 @@
     with $\Psi$ = 0.
 
-Having established the mean direction of the particle we can then apply
-angular orientation distributions. This is done by a numerical integration
-over a range of angles in a similar way to particle size dispersity.
-In the current version of sasview the orientational dispersity is defined
-with respect to the axes of the particle.
+Having established the mean direction of the particle (the view) we can then
+apply angular orientation distributions (jitter). This is done by a numerical
+integration over a range of angles in a similar way to particle size
+dispersity. The orientation dispersity is defined with respect to the
+$a$-$b$-$c$ axes of the particle, with roll angle $\Psi$ about the $c$-axis,
+yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the
+$a$-axis.
+
+More formally, starting with axes $a$-$b$-$c$ of the particle aligned
+with axes $x$-$y$-$z$ of the laboratory frame, the orientation dispersity
+is applied first, using the
+`Tait-Bryan <https://en.wikipedia.org/wiki/Euler_angles#Conventions_2>`_
+$x$-$y'$-$z''$ convention with angles $\Delta\phi$-$\Delta\theta$-$\Delta\Psi$.
+The reference orientation then follows, using the
+`Euler angles <https://en.wikipedia.org/wiki/Euler_angles#Conventions>`_
+$z$-$y'$-$z''$ with angles $\phi$-$\theta$-$\Psi$.  This is implemented
+using rotation matrices as
+
+.. math::
+
+    R = R_z(\phi)\, R_y(\theta)\, R_z(\Psi)\,
+        R_x(\Delta\phi)\, R_y(\Delta\theta)\, R_z(\Delta\Psi)
+
+To transform detector $(q_x, q_y)$ values into $(q_a, q_b, q_c)$ for the
+shape in its canonical orientation, use
+
+.. math::
+
+    [q_a, q_b, q_c]^T = R^{-1} \, [q_x, q_y, 0]^T
+
+
+The inverse rotation is easily calculated by rotating the opposite directions
+in the reverse order, so
+
+.. math::
+
+    R^{-1} = R_z(-\Delta\Psi)\, R_y(-\Delta\theta)\, R_x(-\Delta\phi)\,
+             R_z(-\Psi)\, R_y(-\theta)\, R_z(-\phi)
+
 
 The $\theta$ and $\phi$ orientation parameters for the cylinder only appear
-when fitting 2d data. On introducing "Orientational Distribution" in
-the angles, "distribution of theta" and "distribution of phi" parameters will
+when fitting 2d data. On introducing "Orientational Distribution" in the
+angles, "distribution of theta" and "distribution of phi" parameters will
 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$
-of the cylinder, the $b$ and $a$ axes of the cylinder cross section. (When
-$\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the
-instrument.) The third orientation distribution, in $\Psi$, is about the $c$
-axis of the particle. Some experimentation may be required to understand the
-2d patterns fully. A number of different shapes of distribution are
-available, as described for polydispersity, see :ref:`polydispersityhelp` .
+of the cylinder, which correspond to the $b$ and $a$ axes of the cylinder
+cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and
+$X$ axes of the instrument.) The third orientation distribution, in $\Psi$,
+is about the $c$ axis of the particle. Some experimentation may be required
+to understand the 2d patterns fully. A number of different shapes of
+distribution are available, as described for size dispersity, see
+:ref:`polydispersityhelp`.
 
-Earlier versions of SasView had numerical integration issues in some
-circumstances when distributions passed through 90 degrees. The distributions
-in particle coordinates are more robust, but should still be approached with
-care for large ranges of angle.
+Given that the angular dispersion distribution is defined in cartesian space,
+over a cube defined by
+
+.. math::
+
+    [-\Delta \theta, \Delta \theta] \times
+    [-\Delta \phi, \Delta \phi] \times
+    [-\Delta \Psi, \Delta \Psi]
+
+but the orientation is defined over a sphere, we are left with a
+`map projection <https://en.wikipedia.org/wiki/List_of_map_projections>`_
+problem, with different tradeoffs depending on how values in $\Delta\theta$
+and $\Delta\phi$ are translated into latitude/longitude on the sphere.
+
+Sasmodels is using the
+`equirectangular projection <https://en.wikipedia.org/wiki/Equirectangular_projection>`_.
+In this projection, square patches in angular dispersity become wedge-shaped
+patches on the sphere. To correct for the changing point density, there is a
+scale factor of $\sin(\Delta\theta)$ that applies to each point in the
+integral. This is not enough, though. Consider a shape which is tumbling
+freely around the $b$ axis, with $\Delta\theta$ uniform in $[-180, 180]$. At
+$\pm 90$, all points in $\Delta\phi$ map to the pole, so the jitter will have
+a distinct angular preference. If the spin axis is along the beam (which
+will be the case for $\theta=90$ and $\Psi=90$) the scattering pattern
+should be circularly symmetric, but it will go to zero at $q_x = 0$ due to the
+$\sin(\Delta\theta)$ correction. This problem does not appear for a shape
+that is tumbling freely around the $a$ axis, with $\Delta\phi$ uniform in
+$[-180, 180]$, so swap the $a$ and $b$ axes so $\Delta\theta < \Delta\phi$
+and adjust $\Psi$ by 90. This works with the current sasmodels shapes due to
+symmetry.
+
+Alternative projections were considered.
+The `sinusoidal projection <https://en.wikipedia.org/wiki/Sinusoidal_projection>`_
+works by scaling $\Delta\phi$ as $\Delta\theta$ increases, and dropping those
+points outside $[-180, 180]$. The distortions are a little less for middle
+ranges of $\Delta\theta$, but they are still severe for large $\Delta\theta$
+and the model is much harder to explain.
+The `azimuthal equidistance projection <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>`_
+also improves on the equirectangular projection by extending the range of
+reasonable values for the $\Delta\theta$ range, with $\Delta\phi$ forming a
+wedge that cuts to the opposite side of the sphere rather than cutting to the
+pole. This projection has the nice property that distance from the center are
+preserved, and that $\Delta\theta$ and $\Delta\phi$ act the same.
+The `azimuthal equal area projection <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_
+is like the azimuthal equidistance projection, but it preserves area instead
+of distance. It also has the same behaviour for $\Delta\theta$ and $\Delta\phi$.
+The `Guyou projection <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>`_
+has an excellent balance with reasonable distortion in both $\Delta\theta$
+and $\Delta\phi$, as well as preserving small patches. However, it requires
+considerably more computational overhead, and we have not yet derived the
+formula for the distortion correction, measuring the degree of stretch at
+the point $(\Delta\theta, \Delta\phi)$ on the map.
 
 .. note::
-    Note that the form factors for oriented particles are also performing
-    numerical integrations over one or more variables, so care should be taken,
-    especially with very large particles or more extreme aspect ratios. In such 
-    cases results may not be accurate, particularly at very high Q, unless the model
-    has been specifically coded to use limiting forms of the scattering equations.
-    
-    For best numerical results keep the $\theta$ distribution narrower than the $\phi$ 
-    distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may 
-    need to reorder the sizes of the three axes to acheive the desired result. 
-    This is due to the issues of mapping a rectangular distribution onto the 
-    surface of a sphere.
+    Note that the form factors for oriented particles are performing
+    numerical integrations over one or more variables, so care should be
+    taken, especially with very large particles or more extreme aspect
+    ratios. In such cases results may not be accurate, particularly at very
+    high Q, unless the model has been specifically coded to use limiting
+    forms of the scattering equations.
 
-Users can experiment with the values of *Npts* and *Nsigs*, the number of steps 
-used in the integration and the range spanned in number of standard deviations. 
-The standard deviation is entered in units of degrees. For a "rectangular" 
-distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations. 
-The new "uniform" distribution avoids this by letting you directly specify the 
+    For best numerical results keep the $\theta$ distribution narrower than
+    the $\phi$ distribution. Thus for asymmetric particles, such as
+    elliptical_cylinder, you may need to reorder the sizes of the three axes
+    to acheive the desired result. This is due to the issues of mapping a
+    rectanglar distribution onto the surface of a sphere.
+
+Users can experiment with the values of *Npts* and *Nsigs*, the number of steps
+used in the integration and the range spanned in number of standard deviations.
+The standard deviation is entered in units of degrees. For a "rectangular"
+distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.
+The new "uniform" distribution avoids this by letting you directly specify the
 half width.
 
-The angular distributions will be truncated outside of the range -180 to +180 
-degrees, so beware of using saying a broad Gaussian distribution with large value
-of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would 
-give a good integration,as well as becoming rather meaningless. (At some point 
-in the future the actual polydispersity arrays may be made available to the user 
-for inspection.)
+The angular distributions may be truncated outside of the range -180 to +180
+degrees, so beware of using saying a broad Gaussian distribution with large
+value of *Nsigs*, as the array of *Npts* may be truncated to many fewer
+points than would give a good integration,as well as becoming rather
+meaningless. (At some point in the future the actual dispersion arrays may be
+made available to the user for inspection.)
 
 Some more detailed technical notes are provided in the developer section of
 this manual :ref:`orientation_developer` .
 
+This definition of orientation is new to SasView 4.2.  In earlier versions,
+the orientation distribution appeared as a distribution of view angles.
+This led to strange effects when $c$ was aligned with $z$, where changes
+to the $\phi$ angle served only to rotate the shape about $c$, rather than
+having a consistent interpretation as the pitch of the shape relative to
+the flow field defining the reference orientation.  Prior to SasView 4.1,
+the reference orientation was defined using a Tait-Bryan convention, making
+it difficult to control.  Now, rotation in $\theta$ modifies the spacings
+in the refraction pattern, and rotation in $\phi$ rotates it in the detector
+plane.
+
+
 *Document History*
 
 | 2017-11-06 Richard Heenan
+| 2017-12-20 Paul Kienzle
Index: doc/guide/plugin.rst
===================================================================
--- doc/guide/plugin.rst	(revision c654160fb828fbdd130e49cbf46a616ce82875bb)
+++ doc/guide/plugin.rst	(revision 7e6bc45eddf06819ea47d06f521d0c27c582996e)
@@ -538,16 +538,11 @@
 If the scattering is dependent on the orientation of the shape, then you
 will need to include *orientation* parameters *theta*, *phi* and *psi*
-at the end of the parameter table.  Shape orientation uses *a*, *b* and *c*
-axes, corresponding to the *x*, *y* and *z* axes in the laboratory coordinate
-system, with *z* along the beam and *x*-*y* in the detector plane, with *x*
-horizontal and *y* vertical.  The *psi* parameter rotates the shape
-about its *c* axis, the *theta* parameter then rotates the *c* axis toward
-the *x* axis of the detector, then *phi* rotates the shape in the detector
-plane.  (Prior to these rotations, orientation dispersity will be applied
-as roll-pitch-yaw, rotating *c*, then *b* then *a* in the shape coordinate
-system.)  A particular *qx*, *qy* point on the detector, then corresponds
-to *qa*, *qb*, *qc* with respect to the shape.
-
-The oriented C model is called as *Iqabc(qa, qb, qc, par1, par2, ...)* where
+at the end of the parameter table.  As described in the section
+:ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will
+be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its
+canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the
+laboratory frame and beam travelling along $-z$.
+
+The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where
 *par1*, etc. are the parameters to the model.  If the shape is rotationally
 symmetric about *c* then *psi* is not needed, and the model is called
Index: explore/jitter.py
===================================================================
--- explore/jitter.py	(revision ff104792a8c7664b8c2b98af11a2db58f784b8d1)
+++ explore/jitter.py	(revision 8cfb486c14c2f481e9135c412b0ed9506517c141)
@@ -165,4 +165,5 @@
     # constants in kernel_iq.c
     'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
+    'azimuthal_equal_area',
 ]
 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
Index: n_c.bat
===================================================================
--- run_c.bat	(revision 71c5f4dba5152d6c32ed27f1ae53af28d1931ab8)
+++ 	(revision )
@@ -1,3 +1,0 @@
-rm xlate.c
-python sascomp cylinder,_cylpy -double!
-gcc -c xlate.c 2> err.txt
Index: sasmodels/autoc.py
===================================================================
--- sasmodels/autoc.py	(revision 67cc0fff8de4365704eb57c72afe343290c477d4)
+++ sasmodels/autoc.py	(revision 15be191e7e24e07018a77657389eb80b7dfcad55)
@@ -4,7 +4,5 @@
 from __future__ import print_function
 
-import ast
 import inspect
-from functools import reduce
 
 import numpy as np
@@ -95,12 +93,12 @@
                 constants[name] = obj
                 # Claim all constants are declared on line 1
-                snippets.append('#line 1 "%s"'%escaped_filename)
-                snippets.append(define_constant(name, obj))
+                snippets.append('#line 1 "%s"\n'%escaped_filename)
+                snippets.append(py2c.define_constant(name, obj))
             elif isinstance(obj, special.Gauss):
                 for var, value in zip(("N", "Z", "W"), (obj.n, obj.z, obj.w)):
                     var = "GAUSS_"+var
                     constants[var] = value
-                    snippets.append('#line 1 "%s"'%escaped_filename)
-                    snippets.append(define_constant(var, value))
+                    snippets.append('#line 1 "%s"\n'%escaped_filename)
+                    snippets.append(py2c.define_constant(var, value))
                 #libs.append('lib/gauss%d.c'%obj.n)
                 source = (source.replace(name+'.n', 'GAUSS_N')
@@ -121,5 +119,5 @@
 
     # translate source
-    ordered_code = [code[name] for name in ordered_dag(depends) if name in code]
+    ordered_code = [code[name] for name in py2c.ordered_dag(depends) if name in code]
     functions = py2c.translate(ordered_code, constants)
     snippets.extend(functions)
@@ -127,46 +125,4 @@
     # update model info
     info.source = unique_libs
-    info.c_code = "\n".join(snippets)
+    info.c_code = "".join(snippets)
     info.Iq = info.Iqac = info.Iqabc = info.Iqxy = info.form_volume = None
-
-def define_constant(name, value):
-    if isinstance(value, int):
-        parts = ["int ", name, " = ", "%d"%value, ";"]
-    elif isinstance(value, float):
-        parts = ["double ", name, " = ", "%.15g"%value, ";"]
-    else:
-        # extend constant arrays to a multiple of 4; not sure if this
-        # is necessary, but some OpenCL targets broke if the number
-        # of parameters in the parameter table was not a multiple of 4,
-        # so do it for all constant arrays to be safe.
-        if len(value)%4 != 0:
-            value = list(value) + [0.]*(4 - len(value)%4)
-        elements = ["%.15g"%v for v in value]
-        parts = ["double ", name, "[]", " = ",
-                 "{\n   ", ", ".join(elements), "\n};"]
-    return "".join(parts)
-
-
-# Modified from the following:
-#
-#    http://code.activestate.com/recipes/578272-topological-sort/
-#    Copyright (C) 2012 Sam Denton
-#    License: MIT
-def ordered_dag(dag):
-    # type: (Dict[T, Set[T]]) -> Iterator[T]
-    dag = dag.copy()
-
-    # make leaves depend on the empty set
-    leaves = reduce(set.union, dag.values()) - set(dag.keys())
-    dag.update({node: set() for node in leaves})
-    while True:
-        leaves = set(node for node, links in dag.items() if not links)
-        if not leaves:
-            break
-        for node in leaves:
-            yield node
-        dag = {node: (links-leaves)
-               for node, links in dag.items() if node not in leaves}
-    if dag:
-        raise ValueError("Cyclic dependes exists amongst these items:\n%s"
-                            % ", ".join(str(node) for node in dag.keys()))
Index: smodels/codegen.py
===================================================================
--- sasmodels/codegen.py	(revision 59ee4dbe6d75245e5e0e840fdc8f449c1139d3c2)
+++ 	(revision )
@@ -1,686 +1,0 @@
-"""
-    cdegen
-    ~~~~~~~
-
-    Extension to ast that allow ast -> python code generation.
-
-    :copyright: Copyright 2008 by Armin Ronacher.
-    :license: BSD.
-"""
-from __future__ import print_function, division
-
-import ast
-from ast import NodeVisitor
-
-BINOP_SYMBOLS = {}
-BINOP_SYMBOLS[ast.Add] = '+'
-BINOP_SYMBOLS[ast.Sub] = '-'
-BINOP_SYMBOLS[ast.Mult] = '*'
-BINOP_SYMBOLS[ast.Div] = '/'
-BINOP_SYMBOLS[ast.Mod] = '%'
-BINOP_SYMBOLS[ast.Pow] = '**'
-BINOP_SYMBOLS[ast.LShift] = '<<'
-BINOP_SYMBOLS[ast.RShift] = '>>'
-BINOP_SYMBOLS[ast.BitOr] = '|'
-BINOP_SYMBOLS[ast.BitXor] = '^'
-BINOP_SYMBOLS[ast.BitAnd] = '&'
-BINOP_SYMBOLS[ast.FloorDiv] = '//'
-
-BOOLOP_SYMBOLS = {}
-BOOLOP_SYMBOLS[ast.And] = 'and'
-BOOLOP_SYMBOLS[ast.Or] = 'or'
-
-CMPOP_SYMBOLS = {}
-CMPOP_SYMBOLS[ast.Eq] = '=='
-CMPOP_SYMBOLS[ast.NotEq] = '!='
-CMPOP_SYMBOLS[ast.Lt] = '<'
-CMPOP_SYMBOLS[ast.LtE] = '<='
-CMPOP_SYMBOLS[ast.Gt] = '>'
-CMPOP_SYMBOLS[ast.GtE] = '>='
-CMPOP_SYMBOLS[ast.Is] = 'is'
-CMPOP_SYMBOLS[ast.IsNot] = 'is not'
-CMPOP_SYMBOLS[ast.In] = 'in'
-CMPOP_SYMBOLS[ast.NotIn] = 'not in'
-
-UNARYOP_SYMBOLS = {}
-UNARYOP_SYMBOLS[ast.Invert] = '~'
-UNARYOP_SYMBOLS[ast.Not] = 'not'
-UNARYOP_SYMBOLS[ast.UAdd] = '+'
-UNARYOP_SYMBOLS[ast.USub] = '-'
-
-
-def translate(functions, constants):
-    # type: (List[Tuple[str, str, int]], Dict[str, Any]) -> str
-    snippets = []
-    for source, filename, offset in functions:
-        tree = ast.parse(source)
-        snippet = to_source(tree) #, filename, offset)
-        snippets.append(snippet)
-    return "\n".join(snippets)
-
-def to_source(node, indent_with=' ' * 4, add_line_information=False):
-    """This function can convert a node tree back into python sourcecode.
-    This is useful for debugging purposes, especially if you're dealing with
-    custom asts not generated by python itself.
-
-    It could be that the sourcecode is evaluable when the AST itself is not
-    compilable / evaluable.  The reason for this is that the AST contains some
-    more data than regular sourcecode does, which is dropped during
-    conversion.
-
-    Each level of indentation is replaced with `indent_with`.  Per default this
-    parameter is equal to four spaces as suggested by PEP 8, but it might be
-    adjusted to match the application's styleguide.
-
-    If `add_line_information` is set to `True` comments for the line numbers
-    of the nodes are added to the output.  This can be used to spot wrong line
-    number information of statement nodes.
-    """
-    generator = SourceGenerator(indent_with, add_line_information)
-    generator.visit(node)
-
-    return ''.join(generator.result)
-
-class SourceGenerator(NodeVisitor):
-    """This visitor is able to transform a well formed syntax tree into python
-    sourcecode.  For more details have a look at the docstring of the
-    `node_to_source` function.
-    """
-
-    def __init__(self, indent_with, add_line_information=False):
-        self.result = []
-        self.indent_with = indent_with
-        self.add_line_information = add_line_information
-        self.indentation = 0
-        self.new_lines = 0
-
-    def write(self, x):
-        if self.new_lines:
-            if self.result:
-                self.result.append('\n' * self.new_lines)
-            self.result.append(self.indent_with * self.indentation)
-            self.new_lines = 0
-        self.result.append(x)
-
-    def newline(self, node=None, extra=0):
-        self.new_lines = max(self.new_lines, 1 + extra)
-        if node is not None and self.add_line_information:
-            self.write('# line: %s' % node.lineno)
-            self.new_lines = 1
-
-    def body(self, statements):
-        self.new_line = True
-        self.indentation += 1
-        for stmt in statements:
-            self.visit(stmt)
-        self.indentation -= 1
-
-    def body_or_else(self, node):
-        self.body(node.body)
-        if node.orelse:
-            self.newline()
-            self.write('else:')
-            self.body(node.orelse)
-
-    def signature(self, node):
-        want_comma = []
-        def write_comma():
-            if want_comma:
-                self.write(', ')
-            else:
-                want_comma.append(True)
-
-        padding = [None] * (len(node.args) - len(node.defaults))
-        for arg, default in zip(node.args, padding + node.defaults):
-            write_comma()
-            self.visit(arg)
-            if default is not None:
-                self.write('=')
-                self.visit(default)
-        if node.vararg is not None:
-            write_comma()
-            try:
-                self.write('*' + node.vararg.arg)
-            except AttributeError:  # CRUFT: python 2
-                self.write('*' + node.vararg)
-        if node.kwarg is not None:
-            write_comma()
-            try:
-                self.write('**' + node.kwarg.arg)
-            except AttributeError:  # CRUFT: python 2
-                self.write('*' + node.vararg)
-
-    def decorators(self, node):
-        for decorator in node.decorator_list:
-            self.newline(decorator)
-            self.write('@')
-            self.visit(decorator)
-
-    # Statements
-
-    def visit_Assert(self, node):
-        self.newline(node)
-        self.write('assert ')
-        self.visit(node.test)
-        if node.msg is not None:
-           self.write(', ')
-           self.visit(node.msg)
-
-    def visit_Assign(self, node):
-        self.newline(node)
-        for idx, target in enumerate(node.targets):
-            if idx:
-                self.write(', ')
-            self.visit(target)
-        self.write(' = ')
-        self.visit(node.value)
-
-    def visit_AugAssign(self, node):
-        self.newline(node)
-        self.visit(node.target)
-        self.write(' ' + BINOP_SYMBOLS[type(node.op)] + '= ')
-        self.visit(node.value)
-
-    def visit_ImportFrom(self, node):
-        self.newline(node)
-        self.write('from %s%s import ' % ('.' * node.level, node.module))
-        for idx, item in enumerate(node.names):
-            if idx:
-                self.write(', ')
-            self.write(item)
-
-    def visit_Import(self, node):
-        self.newline(node)
-        for item in node.names:
-            self.write('import ')
-            self.visit(item)
-
-    def visit_Expr(self, node):
-        self.newline(node)
-        self.generic_visit(node)
-
-    def visit_FunctionDef(self, node):
-        self.newline(extra=1)
-        self.decorators(node)
-        self.newline(node)
-        self.write('def %s(' % node.name)
-        self.visit(node.args)
-        self.write('):')
-        self.body(node.body)
-
-    def visit_ClassDef(self, node):
-        have_args = []
-        def paren_or_comma():
-            if have_args:
-                self.write(', ')
-            else:
-                have_args.append(True)
-                self.write('(')
-
-        self.newline(extra=2)
-        self.decorators(node)
-        self.newline(node)
-        self.write('class %s' % node.name)
-        for base in node.bases:
-            paren_or_comma()
-            self.visit(base)
-        # XXX: the if here is used to keep this module compatible
-        #      with python 2.6.
-        if hasattr(node, 'keywords'):
-            for keyword in node.keywords:
-                paren_or_comma()
-                self.write(keyword.arg + '=')
-                self.visit(keyword.value)
-            if node.starargs is not None:
-                paren_or_comma()
-                self.write('*')
-                self.visit(node.starargs)
-            if node.kwargs is not None:
-                paren_or_comma()
-                self.write('**')
-                self.visit(node.kwargs)
-        self.write(have_args and '):' or ':')
-        self.body(node.body)
-
-    def visit_If(self, node):
-        self.newline(node)
-        self.write('if ')
-        self.visit(node.test)
-        self.write(':')
-        self.body(node.body)
-        while True:
-            else_ = node.orelse
-            if len(else_) == 0:
-                break
-            elif len(else_) == 1 and isinstance(else_[0], ast.If):
-                node = else_[0]
-                self.newline()
-                self.write('elif ')
-                self.visit(node.test)
-                self.write(':')
-                self.body(node.body)
-            else:
-                self.newline()
-                self.write('else:')
-                self.body(else_)
-                break
-
-    def visit_For(self, node):
-        self.newline(node)
-        self.write('for ')
-        self.visit(node.target)
-        self.write(' in ')
-        self.visit(node.iter)
-        self.write(':')
-        self.body_or_else(node)
-
-    def visit_While(self, node):
-        self.newline(node)
-        self.write('while ')
-        self.visit(node.test)
-        self.write(':')
-        self.body_or_else(node)
-
-    def visit_With(self, node):
-        self.newline(node)
-        self.write('with ')
-        self.visit(node.context_expr)
-        if node.optional_vars is not None:
-            self.write(' as ')
-            self.visit(node.optional_vars)
-        self.write(':')
-        self.body(node.body)
-
-    def visit_Pass(self, node):
-        self.newline(node)
-        self.write('pass')
-
-    def visit_Print(self, node):
-        # XXX: python 2.6 only
-        self.newline(node)
-        self.write('print ')
-        want_comma = False
-        if node.dest is not None:
-            self.write(' >> ')
-            self.visit(node.dest)
-            want_comma = True
-        for value in node.values:
-            if want_comma:
-                self.write(', ')
-            self.visit(value)
-            want_comma = True
-        if not node.nl:
-            self.write(',')
-
-    def visit_Delete(self, node):
-        self.newline(node)
-        self.write('del ')
-        for idx, target in enumerate(node):
-            if idx:
-                self.write(', ')
-            self.visit(target)
-
-    def visit_TryExcept(self, node):
-        self.newline(node)
-        self.write('try:')
-        self.body(node.body)
-        for handler in node.handlers:
-            self.visit(handler)
-
-    def visit_TryFinally(self, node):
-        self.newline(node)
-        self.write('try:')
-        self.body(node.body)
-        self.newline(node)
-        self.write('finally:')
-        self.body(node.finalbody)
-
-    def visit_Global(self, node):
-        self.newline(node)
-        self.write('global ' + ', '.join(node.names))
-
-    def visit_Nonlocal(self, node):
-        self.newline(node)
-        self.write('nonlocal ' + ', '.join(node.names))
-
-    def visit_Return(self, node):
-        self.newline(node)
-        if node.value is None:
-            self.write('return')
-        else:
-            self.write('return ')
-            self.visit(node.value)
-
-    def visit_Break(self, node):
-        self.newline(node)
-        self.write('break')
-
-    def visit_Continue(self, node):
-        self.newline(node)
-        self.write('continue')
-
-    def visit_Raise(self, node):
-        # XXX: Python 2.6 / 3.0 compatibility
-        self.newline(node)
-        self.write('raise')
-        if hasattr(node, 'exc') and node.exc is not None:
-            self.write(' ')
-            self.visit(node.exc)
-            if node.cause is not None:
-                self.write(' from ')
-                self.visit(node.cause)
-        elif hasattr(node, 'type') and node.type is not None:
-            self.visit(node.type)
-            if node.inst is not None:
-                self.write(', ')
-                self.visit(node.inst)
-            if node.tback is not None:
-                self.write(', ')
-                self.visit(node.tback)
-
-    # Expressions
-
-    def visit_Attribute(self, node):
-        self.visit(node.value)
-        self.write('.' + node.attr)
-
-    def visit_Call(self, node):
-        want_comma = []
-        def write_comma():
-            if want_comma:
-                self.write(', ')
-            else:
-                want_comma.append(True)
-
-        self.visit(node.func)
-        self.write('(')
-        for arg in node.args:
-            write_comma()
-            self.visit(arg)
-        for keyword in node.keywords:
-            write_comma()
-            self.write(keyword.arg + '=')
-            self.visit(keyword.value)
-        if getattr(node, 'starargs', None) is not None:
-            write_comma()
-            self.write('*')
-            self.visit(node.starargs)
-        if getattr(node, 'kwargs', None) is not None:
-            write_comma()
-            self.write('**')
-            self.visit(node.kwargs)
-        self.write(')')
-
-    def visit_Name(self, node):
-        self.write(node.id)
-
-    def visit_Str(self, node):
-        self.write(repr(node.s))
-
-    def visit_Bytes(self, node):
-        self.write(repr(node.s))
-
-    def visit_Num(self, node):
-        self.write(repr(node.n))
-
-    def visit_Tuple(self, node):
-        self.write('(')
-        idx = -1
-        for idx, item in enumerate(node.elts):
-            if idx:
-                self.write(', ')
-            self.visit(item)
-        self.write(idx and ')' or ',)')
-
-    def sequence_visit(left, right):
-        def visit(self, node):
-            self.write(left)
-            for idx, item in enumerate(node.elts):
-                if idx:
-                    self.write(', ')
-                self.visit(item)
-            self.write(right)
-        return visit
-
-    visit_List = sequence_visit('[', ']')
-    visit_Set = sequence_visit('{', '}')
-    del sequence_visit
-
-    def visit_Dict(self, node):
-        self.write('{')
-        for idx, (key, value) in enumerate(zip(node.keys, node.values)):
-            if idx:
-                self.write(', ')
-            self.visit(key)
-            self.write(': ')
-            self.visit(value)
-        self.write('}')
-
-    def visit_BinOp(self, node):
-        self.visit(node.left)
-        self.write(' %s ' % BINOP_SYMBOLS[type(node.op)])
-        self.visit(node.right)
-
-    def visit_BoolOp(self, node):
-        self.write('(')
-        for idx, value in enumerate(node.values):
-            if idx:
-                self.write(' %s ' % BOOLOP_SYMBOLS[type(node.op)])
-            self.visit(value)
-        self.write(')')
-
-    def visit_Compare(self, node):
-        self.write('(')
-        self.visit(node.left)
-        for op, right in zip(node.ops, node.comparators):
-            self.write(' %s ' % CMPOP_SYMBOLS[type(op)])
-            self.visit(right)
-        self.write(')')
-
-    def visit_UnaryOp(self, node):
-        self.write('(')
-        op = UNARYOP_SYMBOLS[type(node.op)]
-        self.write(op)
-        if op == 'not':
-            self.write(' ')
-        self.visit(node.operand)
-        self.write(')')
-
-    def visit_Subscript(self, node):
-        self.visit(node.value)
-        self.write('[')
-        self.visit(node.slice)
-        self.write(']')
-
-    def visit_Slice(self, node):
-        if node.lower is not None:
-            self.visit(node.lower)
-        self.write(':')
-        if node.upper is not None:
-            self.visit(node.upper)
-        if node.step is not None:
-            self.write(':')
-            if not (isinstance(node.step, ast.Name) and node.step.id == 'None'):
-                self.visit(node.step)
-
-    def visit_ExtSlice(self, node):
-        for idx, item in node.dims:
-            if idx:
-                self.write(', ')
-            self.visit(item)
-
-    def visit_Yield(self, node):
-        self.write('yield ')
-        self.visit(node.value)
-
-    def visit_Lambda(self, node):
-        self.write('lambda ')
-        self.visit(node.args)
-        self.write(': ')
-        self.visit(node.body)
-
-    def visit_Ellipsis(self, node):
-        self.write('Ellipsis')
-
-    def generator_visit(left, right):
-        def visit(self, node):
-            self.write(left)
-            self.visit(node.elt)
-            for comprehension in node.generators:
-                self.visit(comprehension)
-            self.write(right)
-        return visit
-
-    visit_ListComp = generator_visit('[', ']')
-    visit_GeneratorExp = generator_visit('(', ')')
-    visit_SetComp = generator_visit('{', '}')
-    del generator_visit
-
-    def visit_DictComp(self, node):
-        self.write('{')
-        self.visit(node.key)
-        self.write(': ')
-        self.visit(node.value)
-        for comprehension in node.generators:
-            self.visit(comprehension)
-        self.write('}')
-
-    def visit_IfExp(self, node):
-        self.visit(node.body)
-        self.write(' if ')
-        self.visit(node.test)
-        self.write(' else ')
-        self.visit(node.orelse)
-
-    def visit_Starred(self, node):
-        self.write('*')
-        self.visit(node.value)
-
-    def visit_Repr(self, node):
-        # XXX: python 2.6 only
-        self.write('`')
-        self.visit(node.value)
-        self.write('`')
-
-    # Helper Nodes
-
-    def visit_alias(self, node):
-        self.write(node.name)
-        if node.asname is not None:
-            self.write(' as ' + node.asname)
-
-    def visit_comprehension(self, node):
-        self.write(' for ')
-        self.visit(node.target)
-        self.write(' in ')
-        self.visit(node.iter)
-        if node.ifs:
-            for if_ in node.ifs:
-                self.write(' if ')
-                self.visit(if_)
-
-    def visit_excepthandler(self, node):
-        self.newline(node)
-        self.write('except')
-        if node.type is not None:
-            self.write(' ')
-            self.visit(node.type)
-            if node.name is not None:
-                self.write(' as ')
-                self.visit(node.name)
-        self.write(':')
-        self.body(node.body)
-
-    def visit_arguments(self, node):
-        self.signature(node)
-
-# ===== inspect.getclosurevars backport begin =====
-
-# copied from python 3
-from inspect import ismethod, isfunction, ismodule
-from collections import namedtuple
-builtins = __builtins__
-ClosureVars = namedtuple('ClosureVars', 'nonlocals globals builtins unbound')
-def py2_getclosurevars(func):
-    """
-    Get the mapping of free variables to their current values.
-
-    Returns a named tuple of dicts mapping the current nonlocal, global
-    and builtin references as seen by the body of the function. A final
-    set of unbound names that could not be resolved is also provided.
-    """
-
-    if ismethod(func):
-        func = func.__func__
-
-    if not isfunction(func):
-        raise TypeError("'{!r}' is not a Python function".format(func))
-
-    code = func.__code__
-    # Nonlocal references are named in co_freevars and resolved
-    # by looking them up in __closure__ by positional index
-    if func.__closure__ is None:
-        nonlocal_vars = {}
-    else:
-        nonlocal_vars = {
-            var : cell.cell_contents
-            for var, cell in zip(code.co_freevars, func.__closure__)
-       }
-
-    # Global and builtin references are named in co_names and resolved
-    # by looking them up in __globals__ or __builtins__
-    global_ns = func.__globals__
-    builtin_ns = global_ns.get("__builtins__", builtins.__dict__)
-    if ismodule(builtin_ns):
-        builtin_ns = builtin_ns.__dict__
-    global_vars = {}
-    builtin_vars = {}
-    unbound_names = set()
-    for name in code.co_names:
-        if name in ("None", "True", "False"):
-            # Because these used to be builtins instead of keywords, they
-            # may still show up as name references. We ignore them.
-            continue
-        try:
-            global_vars[name] = global_ns[name]
-        except KeyError:
-            try:
-                builtin_vars[name] = builtin_ns[name]
-            except KeyError:
-                unbound_names.add(name)
-
-    return ClosureVars(nonlocal_vars, global_vars,
-                       builtin_vars, unbound_names)
-
-import inspect
-if not hasattr(inspect, 'getclosurevars'):
-    inspect.getclosurevars = py2_getclosurevars
-
-# ===== inspect.getclosurevars backport end ======
-
-def print_function(f=None):
-    """
-    Print out the code for the function
-    """
-    # Include some comments to see if they get printed
-    import ast
-    import inspect
-    print("function:", f.__code__.co_name)
-    print("closure:", inspect.getclosurevars(f))
-    print("locals:", f.__code__.co_nlocals, f.__code__.co_names)
-    if f is not None:
-        tree = ast.parse(inspect.getsource(f))
-        print(to_source(tree))
-
-from math import sin
-outside = '40 C'
-def _hello(*args, **kw):
-    x = sin(y) + cos(z.real)
-    print("world", outside)
-    a, b = x, y
-    [a, b] = (x, y)
-    (a, b) = [x, y]
-
-if __name__ == "__main__":
-    #print_function(print_function)
-    print_function(_hello)
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 67cc0fff8de4365704eb57c72afe343290c477d4)
+++ sasmodels/modelinfo.py	(revision 5ab99b7733c7a08b8ef5bf95cdbd369b3b0a046a)
@@ -794,5 +794,4 @@
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
-    info.c_code = getattr(kernel_module, 'c_code', None)
     info.source = getattr(kernel_module, 'source', [])
     info.c_code = getattr(kernel_module, 'c_code', None)
@@ -810,6 +809,4 @@
     info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore
     # Default single and opencl to True for C models.  Python models have callable Iq.
-    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
-    info.single = getattr(kernel_module, 'single', not callable(info.Iq))
     info.random = getattr(kernel_module, 'random', None)
 
@@ -827,4 +824,9 @@
         except Exception as exc:
             logger.warn(str(exc) + " while converting %s from C to python"%name)
+
+    # Needs to come after autoc.convert since the Iq symbol may have been
+    # converted from python to C
+    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
+    info.single = getattr(kernel_module, 'single', not callable(info.Iq))
 
     if callable(info.Iq) and parameters.has_2d:
Index: sasmodels/models/_cylpy.py
===================================================================
--- sasmodels/models/_cylpy.py	(revision 67cc0fff8de4365704eb57c72afe343290c477d4)
+++ sasmodels/models/_cylpy.py	(revision c01ed3e2bc363908dd28daf057e445c91561293d)
@@ -140,4 +140,5 @@
 py2c = True
 
+# TODO: "#define INVALID (expr)" is not supported
 def invalid(v):
     return v.radius < 0 or v.length < 0
@@ -206,5 +207,5 @@
             phi_pd=10, phi_pd_n=5)
 
-qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
+qx, qy = 0.2 * cos(2.5), 0.2 * sin(2.5)
 # After redefinition of angles, find new tests values.  Was 10 10 in old coords
 tests = [
Index: sasmodels/py2c.py
===================================================================
--- sasmodels/py2c.py	(revision 0bd0877de0c83cef75a40dbe39227fa18c60fdb6)
+++ sasmodels/py2c.py	(revision 0bd0877de0c83cef75a40dbe39227fa18c60fdb6)
@@ -1,81 +1,149 @@
-"""
-    py2c
-    ~~~~
-
-    Convert simple numeric python code into C code.
-
-    The translate() function works on
-
-    Variables definition in C
-    -------------------------
-    Defining variables within the translate function is a bit of a guess work,
-    using following rules:
-    *   By default, a variable is a 'double'.
-    *   Variable in a for loop is an int.
-    *   Variable that is references with brackets is an array of doubles. The
-        variable within the brackets is integer. For example, in the
-        reference 'var1[var2]', var1 is a double array, and var2 is an integer.
-    *   Assignment to an argument makes that argument an array, and the index
-        in that assignment is 0.
-        For example, the following python code::
-            def func(arg1, arg2):
-                arg2 = 17.
-        is translated to the following C code::
-            double func(double arg1)
-            {
-                arg2[0] = 17.0;
-            }
-        For example, the following python code is translated to the
-        following C code::
-
-            def func(arg1, arg2):          double func(double arg1) {
-                arg2 = 17.                      arg2[0] = 17.0;
-                                            }
-    *   All functions are defined as double, even if there is no
-        return statement.
-
-Based on codegen.py:
-
+r"""
+py2c
+~~~~
+
+Convert simple numeric python code into C code.
+
+This code is intended to translate direct algorithms for scientific code
+(mostly if statements and for loops operating on double precision values)
+into C code. Unlike projects like numba, cython, pypy and nuitka, the
+:func:`translate` function returns the corresponding C which can then be
+compiled with tinycc or sent to the GPU using CUDA or OpenCL.
+
+There is special handling certain constructs, such as *for i in range* and
+small integer powers.
+
+**TODO: make a nice list of supported constructs***
+
+Imports are not supported, but they are at least ignored so that properly
+constructed code can be run via python or translated to C without change.
+
+Most other python constructs are **not** supported:
+* classes
+* builtin types (dict, set, list)
+* exceptions
+* with context
+* del
+* yield
+* async
+* list slicing
+* multiple return values
+* "is/is not", "in/not in" conditionals
+
+There is limited support for list and list comprehensions, so long as they
+can be represented by a fixed array whose size is known at compile time, and
+they are small enough to be stored on the stack.
+
+Variables definition in C
+-------------------------
+Defining variables within the translate function is a bit of a guess work,
+using following rules:
+*   By default, a variable is a 'double'.
+*   Variable in a for loop is an int.
+*   Variable that is references with brackets is an array of doubles. The
+    variable within the brackets is integer. For example, in the
+    reference 'var1[var2]', var1 is a double array, and var2 is an integer.
+*   Assignment to an argument makes that argument an array, and the index
+    in that assignment is 0.
+    For example, the following python code::
+        def func(arg1, arg2):
+            arg2 = 17.
+    is translated to the following C code::
+        double func(double arg1)
+        {
+            arg2[0] = 17.0;
+        }
+    For example, the following python code is translated to the
+    following C code::
+
+        def func(arg1, arg2):          double func(double arg1) {
+            arg2 = 17.                      arg2[0] = 17.0;
+                                        }
+*   All functions are defined as double, even if there is no
+    return statement.
+
+Debugging
+---------
+
+*print* is partially supported using a simple regular expression. This
+requires a stylized form. Be sure to use print as a function instead of
+the print statement. If you are including substition variables, use the
+% string substitution style. Include parentheses around the substitution
+tuple, even if there is only one item; do not include the final comma even
+if it is a single item (yes, it won't be a tuple, but it makes the regexp
+much simpler). Keep the item on a single line. Here are three forms that work::
+
+    print("x") => printf("x\n");
+    print("x %g"%(a)) => printf("x %g\n", a);
+    print("x %g %g %g"%(a, b, c)) => printf("x %g %g %g\n", a, b, c);
+
+You can generate *main* using the *if __name__ == "__main__":* construct.
+This does a simple substitution with "def main():" before translation and
+a substitution with "int main(int argc, double *argv[])" after translation.
+The result is that the content of the *if* block becomes the content of *main*.
+Along with the print statement, you can run and test a translation standalone
+using::
+
+    python py2c.py source.py
+    cc source.c
+    ./a.out
+
+Known issues
+------------
+The following constructs may cause problems:
+
+* implicit arrays: possible namespace collision for variable "vec#"
+* swap fails: "x,y = y,x" will set x==y
+* top-level statements: code outside a function body causes errors
+* line number skew: each statement should be tagged with its own #line
+  to avoid skew as comments are skipped and loop bodies are wrapped with
+  braces, etc.
+
+References
+----------
+
+Based on a variant of codegen.py:
+
+    https://github.com/andreif/codegen
     :copyright: Copyright 2008 by Armin Ronacher.
     :license: BSD.
 """
-"""
-Update Notes
-============
-11/22/2017, O.E.   Each 'visit_*' method is to build a C statement string. It
-                    shold insert 4 blanks per indentation level.
-                    The 'body' method will combine all the strings, by adding
-                    the 'current_statement' to the c_proc string list
-   11/2017, OE: variables, argument definition implemented.
-   Note: An argument is considered an array if it is the target of an
-        assignment. In that case it is translated to <var>[0]
-11/27/2017, OE: 'pow' basicly working
-  /12/2017, OE: Multiple assignment: a1,a2,...,an=b1,b2,...bn implemented
-  /12/2017, OE: Power function, including special cases of
-                square(x)(pow(x,2)) and cube(x)(pow(x,3)), implemented in
-                translate_power, called from visit_BinOp
-12/07/2017, OE: Translation of integer division, '\\' in python, implemented
-                in translate_integer_divide, called from visit_BinOp
-12/07/2017, OE: C variable definition handled in 'define_c_vars'
-              : Python integer division, '//', translated to C in
-                'translate_integer_divide'
-12/15/2017, OE: Precedence maintained by writing opening and closing
-                parenthesesm '(',')', in procedure 'visit_BinOp'.
-12/18/2017, OE: Added call to 'add_current_line()' at the beginning
-                of visit_Return
-01/4/2018, O.E. Added functions 'get_files_names()', , 'print_usage()'. The get_files_names()
-                function retrieves the in/out file names form the command line, and returns
-                true/false if the number of parameters is valid. In case of no input
-                parameters, usage is prompt and the program terminates.
-01/4/2018, O.E. 'translate(functions, constants=None)' returns string, instaed of a list
-01/04/2017 O.E. Fixed bug in 'visit_If': visiting node.orelse in case else exists.
-01/05/2018 O.E. New class C_Vector defined. It holds name, type, size and initial values.
-                Assignment to that class in 'sequence_visit'.
-                Reading and inserting declaration in 'insert_c_vars'
-
-
-"""
+
+# Update Notes
+# ============
+# 2017-11-22, OE: Each 'visit_*' method is to build a C statement string. It
+#                 shold insert 4 blanks per indentation level. The 'body'
+#                 method will combine all the strings, by adding the
+#                 'current_statement' to the c_proc string list
+# 2017-11-22, OE: variables, argument definition implemented.  Note: An
+#                 argument is considered an array if it is the target of an
+#                 assignment. In that case it is translated to <var>[0]
+# 2017-11-27, OE: 'pow' basicly working
+# 2017-12-07, OE: Multiple assignment: a1,a2,...,an=b1,b2,...bn implemented
+# 2017-12-07, OE: Power function, including special cases of
+#                 square(x)(pow(x,2)) and cube(x)(pow(x,3)), implemented in
+#                 translate_power, called from visit_BinOp
+# 2017-12-07, OE: Translation of integer division, '\\' in python, implemented
+#                 in translate_integer_divide, called from visit_BinOp
+# 2017-12-07, OE: C variable definition handled in 'define_c_vars'
+#               : Python integer division, '//', translated to C in
+#                 'translate_integer_divide'
+# 2017-12-15, OE: Precedence maintained by writing opening and closing
+#                 parenthesesm '(',')', in procedure 'visit_BinOp'.
+# 2017-12-18, OE: Added call to 'add_current_line()' at the beginning
+#                 of visit_Return
+# 2018-01-03, PK: Update interface for use in sasmodels
+# 2018-01-03, PK: support "expr if cond else expr" syntax
+# 2018-01-03, PK: x//y => (int)((x)/(y)) and x/y => ((double)(x)/(double)(y))
+# 2018-01-03, PK: True/False => true/false
+# 2018-01-03, PK: f(x) was introducing an extra semicolon
+# 2018-01-03, PK: simplistic print function, for debugging
+# 2018-01-03, PK: while expr: ... => while (expr) { ... }
+# 2018-01-04, OE: Fixed bug in 'visit_If': visiting node.orelse in case else exists.
+
+from __future__ import print_function
+
+import sys
 import ast
-import sys
 from ast import NodeVisitor
 
@@ -117,14 +185,5 @@
 
 
-def to_source(tree, constants=None, fname=None, lineno=0):
-    """
-    This function can convert a syntax tree into C sourcecode.
-    """
-    generator = SourceGenerator(constants=constants, fname=fname, lineno=lineno)
-    generator.visit(tree)
-#    c_code = "\n".join(generator.c_proc)
-    c_code = "".join(generator.c_proc)
-    return (c_code, generator.warnings)
-
+# TODO: should not allow eval of arbitrary python
 def isevaluable(s):
     try:
@@ -133,41 +192,4 @@
     except Exception:
         return False
-
-class C_Vector():
-    def __init__(self):
-        self.size = 0
-        self.type = "double"
-        self.values = []
-        self.assign_line = -1
-        self.name = ""
-
-    def get_leading_space(self, indent_with="    ", indent_level = 1):
-        leading_space = ""
-        for n in range(indent_level):
-            leading_space += indent_with
-        return leading_space
-        
-    def item_declare_string(self, indent_with="    ", indent_level = 1):
-        declare_string = self.name + "[" + str(len(self.values)) + "]"
-        return declare_string
-    
-    def declare_string(self, indent_with="    ", indent_level = 1):
-        declare_string = self.get_leading_space(indent_with, indent_level)
-        declare_string += self.type + " " + item_declare_string(indent_with, indent_level)
-#        self.name + "[" + str(len(self.values)) + "];\n"
-        return declare_string
-    
-    def get_assign(self, indent_with="    ", indent_level = 1):
-        assign_loop = []
-        c_string = self.get_leading_space(indent_with, indent_level - 1)
-        c_string += "for (n=0 ; n < " + str(len(self.values)) + " ; n++) {"
-        assign_loop.append(c_string)
-        for n,value in enumerate(self.values):
-            c_string = self.get_leading_space(indent_with, indent_level)
-            c_string += self.name + "[" + str(n) + "] = " + str(value) + ";"
-            assign_loop.append(c_string)
-        c_string = self.get_leading_space(indent_with, indent_level - 1) + "}"
-        assign_loop.append(c_string)
-        return (assign_loop)
 
 class SourceGenerator(NodeVisitor):
@@ -203,4 +225,5 @@
         self.c_vectors = []
         self.c_constants = constants if constants is not None else {}
+        self.in_expr = False
         self.in_subref = False
         self.in_subscript = False
@@ -209,8 +232,4 @@
         self.is_sequence = False
         self.visited_args = False
-        self.inside_if = False
-        self.C_Vectors = []
-        self.assign_target = ""
-        self.assign_c_vector = []
 
     def write_python(self, x):
@@ -248,33 +267,8 @@
             self.add_unique_var(arg.id)
 
-    def write_print(self, node):
-        self.write_c ("printf (")
-        c_str = self.current_statement
-        self.current_statement = ""
-        self.visit(node.args[0])
-        print_params = self.current_statement # save print agruments
-        self.current_statement = c_str
-        print_params = print_params[1:(len(print_params) - 1)] # remove parentheses
-        prcnt = print_params.rfind('%')
-        if prcnt >= 0:
-            format_string = print_params[0:prcnt-1].strip()
-            format_string = format_string[1:len(format_string) - 1]
-            var_string = print_params[prcnt+1:len(print_params)].strip()
-        else:
-            format_string = print_params
-            var_string = ""
-        format_string= format_string[:len(format_string)]
-        self.write_c('"' + format_string + '\\n"')
-        if (len(var_string) > 0):
-            self.write_c(', ' + var_string)
-        self.write_c(");")
-        self.add_current_line()
-
-
-
     def newline(self, node=None, extra=0):
         self.new_lines = max(self.new_lines, 1 + extra)
         if node is not None and self.add_line_information:
-            self.write_c('# line: %s' % node.lineno)
+            self.write_c('// line: %s' % node.lineno)
             self.new_lines = 1
         if self.current_statement:
@@ -297,4 +291,6 @@
         self.body(node.body)
         if node.orelse:
+            self.unsupported(node, "for...else/while...else not supported")
+
             self.newline()
             self.write_c('else:')
@@ -326,6 +322,6 @@
                 except AttributeError:
                     arg_name = arg.id
-                w_str = ("Default Parameters are unknown to C: '%s = %s" \
-                        % (arg_name, str(default.n)))
+                w_str = ("C does not support default parameters: %s=%s"
+                         % (arg_name, str(default.n)))
                 self.warnings.append(w_str)
 
@@ -366,11 +362,12 @@
 
     def add_semi_colon(self):
-        semi_pos = self.current_statement.find(';')
-        if semi_pos > 0.0:
-            self.current_statement = self.current_statement.replace(';', '')
+        #semi_pos = self.current_statement.find(';')
+        #if semi_pos >= 0:
+        #    self.current_statement = self.current_statement.replace(';', '')
         self.write_c(';')
 
     def visit_Assign(self, node):
         self.add_current_line()
+        self.in_expr = True
         for idx, target in enumerate(node.targets): # multi assign, as in 'a = b = c = 7'
             if idx:
@@ -378,6 +375,4 @@
             self.define_c_vars(target)
             self.visit(target)
-            if hasattr(target,'id'):
-                self.assign_target = target.id
         # Capture assigned tuple names, if any
         targets = self.tuples[:]
@@ -387,30 +382,24 @@
         self.visited_args = False
         self.visit(node.value)
-        if len(self.assign_c_vector):
-            for c_statement in self.assign_c_vector:
-                self.add_c_line (c_statement)
-#                self.c_proc.append(c_statement)
-            self.assign_c_vector.clear()
-        else:
+        self.add_semi_colon()
+        self.add_current_line()
+        # Assign tuples to tuples, if any
+        # TODO: doesn't handle swap:  a,b = b,a
+        for target, item in zip(targets, self.tuples):
+            self.visit(target)
+            self.write_c(' = ')
+            self.visit(item)
             self.add_semi_colon()
             self.add_current_line()
-        # Assign tuples to tuples, if any
-        # TODO: doesn't handle swap:  a,b = b,a
-            for target, item in zip(targets, self.tuples):
-                self.visit(target)
-                self.write_c(' = ')
-                self.visit(item)
-                self.add_semi_colon()
-                self.add_current_line()
-            if self.is_sequence and not self.visited_args:
-                for target in node.targets:
-                    if hasattr(target, 'id'):
-                        if target.id in self.c_vars and target.id not in self.c_dcl_pointers:
-                            if target.id not in self.c_dcl_pointers:
-                                self.c_dcl_pointers.append(target.id)
-                                if target.id in self.c_vars:
-                                    self.c_vars.remove(target.id)
+        if self.is_sequence and not self.visited_args:
+            for target in node.targets:
+                if hasattr(target, 'id'):
+                    if target.id in self.c_vars and target.id not in self.c_dcl_pointers:
+                        if target.id not in self.c_dcl_pointers:
+                            self.c_dcl_pointers.append(target.id)
+                            if target.id in self.c_vars:
+                                self.c_vars.remove(target.id)
         self.current_statement = ''
-        self.assign_target = ''
+        self.in_expr = False
 
     def visit_AugAssign(self, node):
@@ -418,8 +407,10 @@
             if node.target.id not in self.arguments:
                 self.c_vars.append(node.target.id)
+        self.in_expr = True
         self.visit(node.target)
         self.write_c(' ' + BINOP_SYMBOLS[type(node.op)] + '= ')
         self.visit(node.value)
         self.add_semi_colon()
+        self.in_expr = False
         self.add_current_line()
 
@@ -441,6 +432,8 @@
 
     def visit_Expr(self, node):
+        #self.in_expr = True
         self.newline(node)
         self.generic_visit(node)
+        #self.in_expr = False
 
     def write_c_pointers(self, start_var):
@@ -469,15 +462,4 @@
             have_decls = True
             start_var += 1
-        if len(self.C_Vectors) > 0:
-            vec_declare = []
-            for c_vec in self.C_Vectors:
-                item_declare = c_vec.item_declare_string()
-                if item_declare not in vec_declare:
-                    vec_declare.append (item_declare)
-                if c_vec.name in self.c_vars:
-                    self.c_vars.remove(c_vec.name)
-            all_declare = "    double " + ", ".join (vec_declare) + ";\n"
-            self.c_proc.insert(start_var, all_declare)
-            start_var += 1
 
         if self.c_vars:
@@ -519,5 +501,9 @@
         self.current_function = node.name
 
+        # remember the location of the next warning that will be inserted
+        # so that we can stuff the function name ahead of the warning list
+        # if any warnings are generated by the function.
         warning_index = len(self.warnings)
+
         self.newline(extra=1)
         self.decorators(node)
@@ -533,12 +519,8 @@
         self.insert_signature()
         self.insert_c_vars(start_vars)
-#        self.c_pointers = []
-        self.c_pointers.clear()
-        self.C_Vectors.clear()
+        del self.c_pointers[:]
         self.current_function = ""
-        if (warning_index != len(self.warnings)):
+        if warning_index != len(self.warnings):
             self.warnings.insert(warning_index, "Warning in function '" + node.name + "':")
-            self.warnings.append("Note: C compilation will fail")
-
 
     def visit_ClassDef(self, node):
@@ -576,8 +558,9 @@
 
     def visit_If(self, node):
-        self.add_current_line()
-        self.inside_if = True
+
         self.write_c('if ')
+        self.in_expr = True
         self.visit(node.test)
+        self.in_expr = False
         self.write_c(' {')
         self.body(node.body)
@@ -592,5 +575,7 @@
                 #self.newline()
                 self.write_c('else if ')
+                self.in_expr = True
                 self.visit(node.test)
+                self.in_expr = False
                 self.write_c(' {')
                 self.body(node.body)
@@ -601,8 +586,7 @@
                 self.newline()
                 self.write_c('else {')
-                self.body(node.orelse)
+                self.body(else_)
                 self.add_c_line('}')
                 break
-        self.inside_if = False
 
     def get_for_range(self, node):
@@ -631,7 +615,7 @@
         return start, stop, step
 
-    def add_c_int_var (self, iterator):
-        if iterator not in self.c_int_vars:
-            self.c_int_vars.append(iterator)
+    def add_c_int_var(self, name):
+        if name not in self.c_int_vars:
+            self.c_int_vars.append(name)
 
     def visit_For(self, node):
@@ -647,7 +631,7 @@
                     iterator = self.current_statement
                     self.current_statement = ''
-                    self.add_c_int_var (iterator)
+                    self.add_c_int_var(iterator)
                     start, stop, step = self.get_for_range(node)
-                    self.write_c("for(" + iterator + "=" + str(start) +
+                    self.write_c("for (" + iterator + "=" + str(start) +
                                  " ; " + iterator + " < " + str(stop) +
                                  " ; " + iterator + " += " + str(step) + ") {")
@@ -664,5 +648,5 @@
             self.write_c(':')
             # report the error
-            self.unsupported("unsupported " + self.current_statement)
+            self.unsupported(node, "unsupported " + self.current_statement)
 
     def visit_While(self, node):
@@ -670,5 +654,5 @@
         self.write_c('while ')
         self.visit(node.test)
-        self.write_c(' { ')
+        self.write_c(' {')
         self.body_or_else(node)
         self.write_c('}')
@@ -677,4 +661,5 @@
     def visit_With(self, node):
         self.unsupported(node)
+
         self.newline(node)
         self.write_python('with ')
@@ -691,6 +676,6 @@
 
     def visit_Print(self, node):
-        # TODO: print support would be nice, though hard to do
         self.unsupported(node)
+
         # CRUFT: python 2.6 only
         self.newline(node)
@@ -711,4 +696,5 @@
     def visit_Delete(self, node):
         self.unsupported(node)
+
         self.newline(node)
         self.write_python('del ')
@@ -720,4 +706,5 @@
     def visit_TryExcept(self, node):
         self.unsupported(node)
+
         self.newline(node)
         self.write_python('try:')
@@ -728,4 +715,5 @@
     def visit_TryFinally(self, node):
         self.unsupported(node)
+
         self.newline(node)
         self.write_python('try:')
@@ -737,4 +725,5 @@
     def visit_Global(self, node):
         self.unsupported(node)
+
         self.newline(node)
         self.write_python('global ' + ', '.join(node.names))
@@ -746,11 +735,12 @@
     def visit_Return(self, node):
         self.add_current_line()
+        self.in_expr = True
         if node.value is None:
             self.write_c('return')
         else:
-            self.write_c('return(')
+            self.write_c('return ')
             self.visit(node.value)
-        self.write_c(')')
         self.add_semi_colon()
+        self.in_expr = False
         self.add_c_line(self.current_statement)
         self.current_statement = ''
@@ -766,4 +756,5 @@
     def visit_Raise(self, node):
         self.unsupported(node)
+
         # CRUFT: Python 2.6 / 3.0 compatibility
         self.newline(node)
@@ -788,4 +779,5 @@
     def visit_Attribute(self, node):
         self.unsupported(node, "attribute reference a.b not supported")
+
         self.visit(node.value)
         self.write_python('.' + node.attr)
@@ -807,7 +799,4 @@
             elif node.func.id == "SINCOS":
                 self.write_sincos(node)
-                return
-            elif node.func.id == "print":
-                self.write_print(node)
                 return
             else:
@@ -835,8 +824,23 @@
                 self.visit(node.kwargs)
         self.write_c(')')
-        if (self.inside_if == False):
-            self.write_c(';')
+        if not self.in_expr:
+            self.add_semi_colon()
+
+    TRANSLATE_CONSTANTS = {
+        # python 2 uses normal name references through vist_Name
+        'True': 'true',
+        'False': 'false',
+        'None': 'NULL',  # "None" will probably fail for other reasons
+        # python 3 uses NameConstant
+        True: 'true',
+        False: 'false',
+        None: 'NULL',  # "None" will probably fail for other reasons
+        }
 
     def visit_Name(self, node):
+        translation = self.TRANSLATE_CONSTANTS.get(node.id, None)
+        if translation:
+            self.write_c(translation)
+            return
         self.write_c(node.id)
         if node.id in self.c_pointers and not self.in_subref:
@@ -853,13 +857,28 @@
                 name not in self.c_constants and not name.isdigit()):
             if self.in_subscript:
-                self.add_c_int_var (self, node.id)
+                self.add_c_int_var(node.id)
             else:
                 self.c_vars.append(node.id)
 
+    def visit_NameConstant(self, node):
+        translation = self.TRANSLATE_CONSTANTS.get(node.value, None)
+        if translation is not None:
+            self.write_c(translation)
+        else:
+            self.unsupported(node, "don't know how to translate %r"%node.value)
+
     def visit_Str(self, node):
-        self.write_c(repr(node.s))
+        s = node.s
+        s = s.replace('\\','\\\\').replace('"','\\"').replace('\n','\\n')
+        self.write_c('"')
+        self.write_c(s)
+        self.write_c('"')
 
     def visit_Bytes(self, node):
-        self.write_c(repr(node.s))
+        s = node.s
+        s = s.replace('\\','\\\\').replace('"','\\"').replace('\n','\\n')
+        self.write_c('"')
+        self.write_c(s)
+        self.write_c('"')
 
     def visit_Num(self, node):
@@ -876,5 +895,4 @@
         def visit(self, node):
             self.is_sequence = True
-            c_vec = C_Vector()
             s = ""
             for idx, item in enumerate(node.elts):
@@ -882,28 +900,11 @@
                     s += ', '
                 if hasattr(item, 'id'):
-                    c_vec.values.append(item.id)
-#                    s += item.id
+                    s += item.id
                 elif hasattr(item, 'n'):
-                    c_vec.values.append(item.n)
-#                    s += str(item.n)
-                else:
-                    temp = self.current_statement
-                    self.current_statement = ""
-                    self.visit(item)
-                    c_vec.values.append(self.current_statement)
-#                    s += self.current_statement
-                    self.current_statement = temp
-            if (self.assign_target):
-                self.add_c_int_var ('n')
-                c_vec.name = self.assign_target
-                self.C_Vectors.append(c_vec)
-                assign_strings = c_vec.get_assign()
-                for assign_str in assign_strings:
-                    self.assign_c_vector.append(assign_str)
-#                    self.add_c_line(assign_str)
-#            if s:
-#                self.c_vectors.append(s)
-#                vec_name = "vec"  + str(len(self.c_vectors))
-#                self.write_c(vec_name)
+                    s += str(item.n)
+            if s:
+                self.c_vectors.append(s)
+                vec_name = "vec"  + str(len(self.c_vectors))
+                self.write_c(vec_name)
         return visit
 
@@ -914,4 +915,5 @@
     def visit_Dict(self, node):
         self.unsupported(node)
+
         self.write_python('{')
         for idx, (key, value) in enumerate(zip(node.keys, node.values)):
@@ -964,12 +966,19 @@
         if is_negative_exp:
             self.write_c(")")
-        self.write_c(" ")
+        #self.write_c(" ")
 
     def translate_integer_divide(self, node):
-        self.write_c("(int)(")
+        self.write_c("(int)((")
         self.visit(node.left)
-        self.write_c(") /(int)(")
+        self.write_c(")/(")
         self.visit(node.right)
-        self.write_c(")")
+        self.write_c("))")
+
+    def translate_float_divide(self, node):
+        self.write_c("((double)(")
+        self.visit(node.left)
+        self.write_c(")/(double)(")
+        self.visit(node.right)
+        self.write_c("))")
 
     def visit_BinOp(self, node):
@@ -979,4 +988,6 @@
         elif '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.FloorDiv]:
             self.translate_integer_divide(node)
+        elif '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.Div]:
+            self.translate_float_divide(node)
         else:
             self.visit(node.left)
@@ -1043,4 +1054,5 @@
     def visit_Yield(self, node):
         self.unsupported(node)
+
         self.write_python('yield ')
         self.visit(node.value)
@@ -1048,4 +1060,5 @@
     def visit_Lambda(self, node):
         self.unsupported(node)
+
         self.write_python('lambda ')
         self.visit(node.args)
@@ -1055,4 +1068,5 @@
     def visit_Ellipsis(self, node):
         self.unsupported(node)
+
         self.write_python('Ellipsis')
 
@@ -1075,4 +1089,5 @@
     def visit_DictComp(self, node):
         self.unsupported(node)
+
         self.write_python('{')
         self.visit(node.key)
@@ -1083,23 +1098,12 @@
         self.write_python('}')
 
-    def visit_NameConstant (self, node):
-        if (hasattr (node, "value")):
-            if (node.value == True):
-                val = "1"
-            elif node.value == False:
-                val = "0"
-            else:
-                val = ""
-            if (len(val) > 0):
-                self.write_c("(" + val + ")")
-
     def visit_IfExp(self, node):
-        self.write_c('(')
+        self.write_c('((')
         self.visit(node.test)
-        self.write_c(' ? ')
+        self.write_c(')?(')
         self.visit(node.body)
-        self.write_c(' : ')
+        self.write_c('):(')
         self.visit(node.orelse)
-        self.write_c(');')
+        self.write_c('))')
 
     def visit_Starred(self, node):
@@ -1117,4 +1121,5 @@
     def visit_alias(self, node):
         self.unsupported(node)
+
         self.write_python(node.name)
         if node.asname is not None:
@@ -1179,7 +1184,7 @@
     const = "constant "  # OpenCL needs globals to be constant
     if isinstance(value, int):
-        parts = [const + "int ", name, " = ", "%d"%value, ";"]
+        parts = [const + "int ", name, " = ", "%d"%value, ";\n"]
     elif isinstance(value, float):
-        parts = [const + "double ", name, " = ", "%.15g"%value, ";"]
+        parts = [const + "double ", name, " = ", "%.15g"%value, ";\n"]
     else:
         try:
@@ -1195,5 +1200,5 @@
         elements = ["%.15g"%v for v in value]
         parts = [const + "double ", name, "[]", " = ",
-                 "{\n   ", ", ".join(elements), "\n};"]
+                 "{\n   ", ", ".join(elements), "\n};\n"]
 
     return "".join(parts)
@@ -1239,44 +1244,53 @@
                          % ", ".join(str(node) for node in dag.keys()))
 
+import re
+PRINT_ARGS = re.compile(r'print[(]"(?P<template>[^"]*)" *% *[(](?P<args>[^\n]*)[)] *[)] *\n')
+SUBST_ARGS = r'printf("\g<template>\\n", \g<args>)\n'
+PRINT_STR = re.compile(r'print[(]"(?P<template>[^"]*)" *[)] *\n')
+SUBST_STR = r'printf("\g<template>\n")'
 def translate(functions, constants=None):
-    # type: (List[(str, str, int)], Dict[str, any]) -> List[str]
+    # type: (Sequence[(str, str, int)], Dict[str, any]) -> List[str]
     """
-    Convert a set of functions
+    Convert a list of functions to a list of C code strings.
+
+    Returns list of corresponding code snippets (with trailing lines in
+    each block) and a list of warnings generated by the translator.
+
+    A function is given by the tuple (source, filename, line number).
+
+    Global constants are given in a dictionary of {name: value}.  The
+    constants are used for name space resolution and type inferencing.
+    Constants are not translated by this code. Instead, call
+    :func:`define_constant` with name and value, and maybe block_size
+    if arrays need to be padded to the next block boundary.
+
+    Function prototypes are not generated. Use :func:`ordered_dag`
+    to list the functions in reverse order of dependency before calling
+    translate. [Maybe a future revision will return the function prototypes
+    so that a suitable "*.h" file can be generated.
     """
     snippets = []
-    snippets.append("#include <math.h>")
-    snippets.append("")
     warnings = []
     for source, fname, lineno in functions:
-        line_directive = '#line %d "%s"'%(lineno, fname.replace('\\', '\\\\'))
+        line_directive = '#line %d "%s"\n'%(lineno, fname.replace('\\', '\\\\'))
         snippets.append(line_directive)
+        # Replace simple print function calls with printf statements
+        source = PRINT_ARGS.sub(SUBST_ARGS, source)
+        source = PRINT_STR.sub(SUBST_STR, source)
         tree = ast.parse(source)
-        c_code,w = to_source(tree, constants=constants, fname=fname, lineno=lineno)
-        if w:
-            warnings.append ("\n".join(w))
+        generator = SourceGenerator(constants=constants, fname=fname, lineno=lineno)
+        generator.visit(tree)
+        c_code = "".join(generator.c_proc)
         snippets.append(c_code)
-    return "\n".join(snippets),warnings
-#    return snippets
-
-def print_usage ():
-        print("""\
-            Usage: python py2c.py <infile> [<outfile>]
-            if outfile is omitted, output file is '<infile>.c'
-        """)
-
-def get_files_names ():
-    import os
-    fname_in = fname_out = ""
-    valid_params =  len(sys.argv) > 1
-    if (valid_params == False):
-        print_usage()
-    else:
-        fname_in = sys.argv[1]
-        if len(sys.argv) == 2:
-            fname_base = os.path.splitext(fname_in)[0]
-            fname_out = str(fname_base) + '.c'
-        else:
-            fname_out = sys.argv[2]
-    return valid_params, fname_in,fname_out
+        warnings.extend(generator.warnings)
+    return snippets, warnings
+
+def to_source(tree, constants=None, fname=None, lineno=0):
+    """
+    This function can convert a syntax tree into C sourcecode.
+    """
+    c_code = "".join(generator.c_proc)
+    return c_code
+
 
 C_HEADER = """
@@ -1299,13 +1313,27 @@
 }
 """
+
+USAGE = """\
+Usage: python py2c.py <infile> [<outfile>]
+
+if outfile is omitted, output file is '<infile>.c'
+"""
+
 def main():
-    print("Parsing...using Python" + sys.version)
-    valid_params, fname_in, fname_out = get_files_names ()
-    if (valid_params == False):
-        print("Input parameters error.\nExiting")
+    import os
+    #print("Parsing...using Python" + sys.version)
+    if len(sys.argv) == 1:
+        print(USAGE)
         return
 
+    fname_in = sys.argv[1]
+    if len(sys.argv) == 2:
+        fname_base = os.path.splitext(fname_in)[0]
+        fname_out = str(fname_base) + '.c'
+    else:
+        fname_out = sys.argv[2]
+
     with open(fname_in, "r") as python_file:
-            code = python_file.read()
+        code = python_file.read()
     name = "gauss"
     code = (code
@@ -1315,17 +1343,17 @@
             .replace('if __name__ == "__main__"', "def main()")
            )
-    translation,warnings = translate([(code, fname_in, 1)])
-    translation = translation.replace("double main()", "int main(int argc, char *argv[])")
+
+    translation, warnings = translate([(code, fname_in, 1)])
+    c_code = "".join(translation)
+    c_code = c_code.replace("double main()", "int main(int argc, char *argv[])")
 
     with open(fname_out, "w") as file_out:
         file_out.write(C_HEADER)
-        file_out.write(str(translation))
+        file_out.write(c_code)
+
     if warnings:
-        print("".join(str(warnings[0])))
-    print("...Done")
+        print("\n".join(warnings))
+    #print("...Done")
 
 if __name__ == "__main__":
-    try:
-        main()
-    except Exception as excp:
-        print ("Error:\n" + str(excp.args))
+    main()
Index: varsall.bat
===================================================================
--- vcvarsall.bat	(revision 71c5f4dba5152d6c32ed27f1ae53af28d1931ab8)
+++ 	(revision )
@@ -1,57 +1,0 @@
-@echo off
-
-if "%1" == "" goto x86
-if not "%2" == "" goto usage
-
-if /i %1 == x86       goto x86
-if /i %1 == amd64     goto amd64
-if /i %1 == x64       goto amd64
-if /i %1 == x86_amd64 goto x86_amd64
-goto usage
-
-:x86
-echo Setting environment for using Microsoft Visual Studio 2008 x86 tools.
-set VCINSTALLDIR=%~dp0VC\
-set WindowsSdkDir=%~dp0WinSDK\
-if not exist "%VCINSTALLDIR%bin\cl.exe" goto missing
-set PATH=%VCINSTALLDIR%Bin;%WindowsSdkDir%Bin;%PATH%
-set INCLUDE=%VCINSTALLDIR%Include;%WindowsSdkDir%Include;%INCLUDE%
-set LIB=%VCINSTALLDIR%Lib;%WindowsSdkDir%Lib;%LIB%
-set LIBPATH=%VCINSTALLDIR%Lib;%WindowsSdkDir%Lib;%LIBPATH%
-goto :eof
-
-:amd64
-echo Setting environment for using Microsoft Visual Studio 2008 x64 tools.
-set VCINSTALLDIR=%~dp0VC\
-set WindowsSdkDir=%~dp0WinSDK\
-if not exist "%VCINSTALLDIR%Bin\amd64\cl.exe" goto missing
-set PATH=%VCINSTALLDIR%Bin\amd64;%WindowsSdkDir%Bin\x64;%WindowsSdkDir%Bin;%PATH%
-set INCLUDE=%VCINSTALLDIR%Include;%WindowsSdkDir%Include;%INCLUDE%
-set LIB=%VCINSTALLDIR%Lib\amd64;%WindowsSdkDir%Lib\x64;%LIB%
-set LIBPATH=%VCINSTALLDIR%Lib\amd64;%WindowsSdkDir%Lib\x64;%LIBPATH%
-goto :eof
-
-:x86_amd64
-echo Setting environment for using Microsoft Visual Studio 2008 x64 cross tools.
-set VCINSTALLDIR=%~dp0VC\
-set WindowsSdkDir=%~dp0WinSDK\
-if not exist "%VCINSTALLDIR%Bin\x86_amd64\cl.exe" goto missing
-set PATH=%VCINSTALLDIR%Bin\x86_amd64;%WindowsSdkDir%Bin;%PATH%
-set INCLUDE=%VCINSTALLDIR%Include;%WindowsSdkDir%Include;%INCLUDE%
-set LIB=%VCINSTALLDIR%Lib\amd64;%WindowsSdkDir%Lib\x64;%LIB%
-set LIBPATH=%VCINSTALLDIR%Lib\amd64;%WindowsSdkDir%Lib\x64;%LIBPATH%
-goto :eof
-
-:usage
-echo Error in script usage. The correct usage is:
-echo     %0 [option]
-echo where [option] is: x86 ^| amd64 ^| x86_amd64
-echo:
-echo For example:
-echo     %0 x86_ia64
-goto :eof
-
-:missing
-echo The specified configuration type is missing.  The tools for the
-echo configuration might not be installed.
-goto :eof
