Regarding AST:

Simon Krughoff points out that the null map has 1400 lines of code (perhaps only 700 of which are actual code). This is intimidating. However, AST does provide a mechanism whereby a function written in C can be wrapped as a mapping. The primary documented limitation is that the mapping does not become part of AST, so any transform based on it is not portable. This seems tolerable to me, since we can provide our extensions as part of the LSST software stack.

I propose to implement a radial polynomial (with the center as a parameter) as an AST extension, to get an idea of what is involved.

Another limitation of registering a C function as an AST mapping is that parameters can only be passed as a single string. Rather than write a radial polynomial, am tempted to write vector <-> norm, unit vector, at which point a radial polynomial is a sequence of:

• apply offset
• vector -> norm, unit vector
• apply a 1-d polynomial to the norm
• norm, unit vector -> vector
• remove offset

An more traditional alternative to vector <-> norm, unit vector is x,y <-> r,theta, but that is less general and likely slower.

Comment from David Berry:

The way to create a radial polynomial mapping would be to take a copy of a similar existing class (I'd suggest PcdMap since it is already pretty close to what you want) and just go through it to change the few lines that implement the transformation, the accessors for parameters, etc. Starting with a case-sensitive global edit of "PcdMap" to "RadialMap" would get you a long way. You would then need to grep for pcdmap throughout all the AST source files and duplicate the required occurrences for RadialMap (just files globals.c and loader.c, unless you also want to create a Fortran API for the new mapping class)

Whilst it is true that unitmap.c contains 1400 lines, in fact only 300 are genuine lines of code (and full commenting of what is going on is surely a good thing rather than a bad thing?). Of those 300 lines of code, the vast majority are close to boiler-plate stuff that is very nearly the same for each class (mainly to do with defining required methods such as attribute accessors, constructors, deleters, etc - all very similar to what you need to do to create a new class in python in fact).

My estimate is that it would take me an hour, plus or minus 15 mins, to create a radial polynomial mapping. Clearly it would take someone else longer, but the fact remains that creating new mappings is a remarkably simple job in AST for anyone who is even reasonably familiar with its inwards. I'll admit some other tasks are harder (like creating a new class of Frame for instance), but creating new mappings seems to be the most popular requirement I'm hearing about.

OK. I"ll try to make a map instead

For my first mapping I wrote UnitNormMap, a map outputs unit vector and norm relative to a specified center. This can then be used with a 1-dimensional PolyMap to make a radial polynomial map, which can be used to model optical distortion.

I modified a copy of astShiftMap. This required creating the actual .c and .h files, and then updating the following files based on searching for ShiftMap and shiftmap:

• Ast.c (the python interface)
• ast.h
• globals.c
• globals.h
• loader.c
• setup.py
• test.py

The changes are trivial in globals.c, globals.h, loader.c, and setup.py but the others require some care.

Additional changes are needed to support FORTRAN. The manual might also want an update, though it is possible that the reference section is automatically generated.

Furthermore, I modified a fork of PyAST, which contains a distribution of AST, but not its unit tests. PyAST has its own tests, but leaves most of the testing of mappings to AST. I added extensive testing for the new mapping in python, but did not attempt to add C tests.

The changes are here: https://github.com/r-owen/starlink-pyast.git on branch unitnormmap. The new map is being considered for inclusion in AST.