Chris,
I was not able to get an astrophyics example, but I do have a reasonable
one
that performs a spectral analysis of the output of an analog to digital
converter, something radio astronomers are known to do.
I am including the code, but it requires a rather large data file to run, which
I will not include. The code uses my 'engfmt' library from pypi to perform
conversion to SI form. In this example, there is no need for conversion from SI
form.
#!/usr/bin/env python3
import numpy as np
from numpy.fft import fft, fftfreq, fftshift
import matplotlib as mpl
mpl.use('SVG')
from matplotlib.ticker import FuncFormatter
import matplotlib.pyplot as pl
from engfmt import Quantity, set_preferences
set_preferences(spacer=' ')
def mag(spectrum):
return np.absolute(spectrum)
def freq_fmt(val, pos):
return Quantity(val, 'Hz').to_eng()
def volt_fmt(val, pos):
return Quantity(val, 'V').to_eng()
freq_formatter = FuncFormatter(freq_fmt)
volt_formatter = FuncFormatter(volt_fmt)
data = np.fromfile('delta-sigma.smpl', sep=' ')
time, wave = data.reshape((2, len(data)//2), order='F')
timestep = time[1] - time[0]
nonperiodicity = wave[-1] - wave[0]
period = timestep * len(time)
print('timestep = {}'.format(Quantity(timestep, 's')))
print('nonperiodicity = {}'.format(Quantity(nonperiodicity, 'V')))
print('timepoints = {}'.format(len(time)))
print('freq resolution = {}'.format(Quantity(1/period, 'Hz')))
window = np.kaiser(len(time), 11)/0.37
# beta=11 corresponds to alpha=3.5 (beta = pi*alpha), /.4 # 0.4 is the
# processing gain with alpha=3.5 is 0.37
#window = 1
windowed = window*wave
spectrum = 2*fftshift(fft(windowed))/len(time)
freq = fftshift(fftfreq(len(wave), timestep))
fig = pl.figure()
ax = fig.add_subplot(111)
ax.plot(freq, mag(spectrum))
ax.set_yscale('log')
ax.xaxis.set_major_formatter(freq_formatter)
ax.yaxis.set_major_formatter(volt_formatter)
pl.savefig('spectrum.svg')
ax.set_xlim((0, 1e6))
pl.savefig('spectrum-zoomed.svg')
When run, this program prints the following diagnostics to stdout:
timestep = 20 ns
nonperiodicity = 2.3 pV
timepoints = 27994
freq resolution = 1.7861 kHz
It also generates two SVG files. I have converted one to PNG and attached it.
A few comments:
1. The data in the input file ('delta-sigma.smpl') has low dynamic range and is
machine generated, and not really meant for direct human consumption. As
such, it does not benefit from using SI scale factors. But there are
certainly cases where the data has both high dynamic range and is intended
for people to examine it directly. In those cases it would be very helpful
if
NumPy was able to directly read the file. As the language exists today,
I would need to read the file myself, manually convert it, and feed the
result to NumPy.
2. Many of these numbers that are output do have high dynamic range and are
intended to be consumed directly by humans. These benefit from using SI
scale
factors. For example, the 'freq resolution' can vary from Hz to MHz and
'nonperiodicity' can vary from fV to mV.
3. Extra effort was expended to make the axis labels on the graph use SI scale
factors so as to make the results 'publication quality'.
My hope is that if Python accepted SI literals directly, then both NumPy nad
MatPlotLib would also be extended to accept/use these formats directly,
eliminating the need for me to do the conversions and manage the axes.
-Ken
On Mon, Aug 29, 2016 at 06:02:29PM +1000, Chris Angelico wrote:
> On Mon, Aug 29, 2016 at 5:07 PM, Ken Kundert
> <python-id...@shalmirane.com> wrote:
> >
> > I talked to astrophysicist about your comments, and what she said was:
> > 1. She would love it if Python had built in support for real numbers with SI
> > scale factors
> > 2. I told her about my library for reading and writing numbers with SI scale
> > factors, and she was much less enthusiastic because using it would
> > require
> > convincing the rest of the group, which would be too much effort.
> > 3. She was amused by the "kilo pico speed of light" comment, but she was
> > adamant
> > that the fact that you, or some system administrator, does not understand
> > what kpc means has absolutely no affect on her desired to use SI scale
> > factors. Her comment: I did not write it for him.
> > 4. She pointed out that the software she writes and uses is intended either
> > for
> > herself of other astrophysicists. No system administrators involved.
>
> So can you share some of her code, and show how the ability to scale
> unitless numbers would help it?
>
> ChrisA
> _______________________________________________
> Python-ideas mailing list
> Python-ideas@python.org
> https://mail.python.org/mailman/listinfo/python-ideas
> Code of Conduct: http://python.org/psf/codeofconduct/
_______________________________________________
Python-ideas mailing list
Python-ideas@python.org
https://mail.python.org/mailman/listinfo/python-ideas
Code of Conduct: http://python.org/psf/codeofconduct/
