Object-oriented implementations of the MPDATA advection equation solver in C++, Python and Fortran

Three object-oriented implementations of a prototype solver of the advection equation are introduced. The presented programs are based on Blitz++ (C++), NumPy (Python), and Fortran's built-in array containers. The solvers include an implementation of the Multidimensional Positive-Definite Advective Transport Algorithm (MPDATA). The introduced codes exemplify how the application of object-oriented programming (OOP) techniques allows to reproduce the mathematical notation used in the literature within the program code. A discussion on the tradeoffs of the programming language choice is presented. The main angles of comparison are code brevity and syntax clarity (and hence maintainability and auditability) as well as performance. In the case of Python, a significant performance gain is observed when switching from the standard interpreter (CPython) to the PyPy implementation of Python. Entire source code of all three implementations is embedded in the text and is licensed under the terms of the GNU GPL license

Extended dust emission and atomic hydrogen, a reservoir of diffuse H_2 in NGC 1068

We report on sensitive sub-mm imaging observations of the prototype Seyfert~2/starburst galaxy NGC 1068 at 850 $\mu$m and 450 $\mu$m using the Submillimetre Common-User Bolometer Array (SCUBA) on the James Clerk Maxwell Telescope (JCMT). We find clear evidence of dust emission associated with the extended HI component which together with the very faint $^{12}$CO J=1--0 emission give a gas-to-dust ratio of $\rm M_{\rm gas}/M_{\rm dust} \sim 70-150$. This contrasts with the larger ratio $\rm M_{\rm gas}/M_{\rm dust}\sim 330$ estimated within a galactocentric radius of $r\leq 1.36$ kpc, where the gas is mostly molecular and starburst activity occurs. The large gas-to-dust ratio found for the starburst region is attributed to a systematic overestimate of the molecular gas mass in starburst environments when the luminosity of the $^{12}$CO J=1--0 line and a standard galactic conversion factor is used. On the other hand sub-mm imaging proves to be a more powerful tool than conventional CO imaging for revealing the properties of the diffuse $\rm H_2$ that coexists with HI. This molecular gas phase is characterized by low densities ($\rm n(H_2)<10^3$ cm$^{-3}$), very faint emission from sub-thermally excited CO, and contains more mass than HI, namely $\rm M(H_2)/M(HI)\sim 5$.Comment: Accepted for publication in the Astrophysical Journal Letter

Runtime-Flexible Multi-dimensional Arrays and Views for C++98 and C++0x

Multi-dimensional arrays are among the most fundamental and most useful data structures of all. In C++, excellent template libraries exist for arrays whose dimension is fixed at runtime. Arrays whose dimension can change at runtime have been implemented in C. However, a generic object-oriented C++ implementation of runtime-flexible arrays has so far been missing. In this article, we discuss our new implementation called Marray, a package of class templates that fills this gap. Marray is based on views as an underlying concept. This concept brings some of the flexibility known from script languages such as R and MATLAB to C++. Marray is free both for commercial and non-commercial use and is publicly available from www.andres.sc/marrayComment: Free source code availabl

Array operators using multiple dispatch: a design methodology for array implementations in dynamic languages

Arrays are such a rich and fundamental data type that they tend to be built into a language, either in the compiler or in a large low-level library. Defining this functionality at the user level instead provides greater flexibility for application domains not envisioned by the language designer. Only a few languages, such as C++ and Haskell, provide the necessary power to define $n$-dimensional arrays, but these systems rely on compile-time abstraction, sacrificing some flexibility. In contrast, dynamic languages make it straightforward for the user to define any behavior they might want, but at the possible expense of performance. As part of the Julia language project, we have developed an approach that yields a novel trade-off between flexibility and compile-time analysis. The core abstraction we use is multiple dispatch. We have come to believe that while multiple dispatch has not been especially popular in most kinds of programming, technical computing is its killer application. By expressing key functions such as array indexing using multi-method signatures, a surprising range of behaviors can be obtained, in a way that is both relatively easy to write and amenable to compiler analysis. The compact factoring of concerns provided by these methods makes it easier for user-defined types to behave consistently with types in the standard library.Comment: 6 pages, 2 figures, workshop paper for the ARRAY '14 workshop, June 11, 2014, Edinburgh, United Kingdo

Primary Beam Shape Calibration from Mosaicked, Interferometric Observations

Image quality in mosaicked observations from interferometric radio telescopes is strongly dependent on the accuracy with which the antenna primary beam is calibrated. The next generation of radio telescope arrays such as the Allen Telescope Array (ATA) and the Square Kilometer Array (SKA) have key science goals that involve making large mosaicked observations filled with bright point sources. We present a new method for calibrating the shape of the telescope's mean primary beam that uses the multiple redundant observations of these bright sources in the mosaic. The method has an analytical solution for simple Gaussian beam shapes but can also be applied to more complex beam shapes through $\chi^2$ minimization. One major benefit of this simple, conceptually clean method is that it makes use of the science data for calibration purposes, thus saving telescope time and improving accuracy through simultaneous calibration and observation. We apply the method both to 1.43 GHz data taken during the ATA Twenty Centimeter Survey (ATATS) and to 3.14 GHz data taken during the ATA's Pi Gigahertz Sky Survey (PiGSS). We find that the beam's calculated full width at half maximum (FWHM) values are consistent with the theoretical values, the values measured by several independent methods, and the values from the simulation we use to demonstrate the effectiveness of our method on data from future telescopes such as the expanded ATA and the SKA. These results are preliminary, and can be expanded upon by fitting more complex beam shapes. We also investigate, by way of a simulation, the dependence of the accuracy of the telescope's FWHM on antenna number. We find that the uncertainty returned by our fitting method is inversely proportional to the number of antennas in the array.Comment: Accepted by PASP. 8 pages, 8 figure