Convert binary .RWD8 format from Avaspec spectrometers into a text file and graphical plots

convert_rwd8_to_txt.py

#!/usr/bin/python3  
#-*- coding: utf-8 -*-

## Import common moduli
import sys 
import matplotlib.pyplot as plt
import numpy as np
import struct

"""
Converts binary .RWD8 format from Avaspec spectrometers into a text file and graphical plots

Takes into account the calibration for dark, as well as integration time
"""

DEAD_PIXELS = [338, 240, 156]


def hsv_to_rgb(h, s=1, v=1, red_norm=.9, green_norm=.8): ## adapted from https://docs.python.org/2/library/colorsys.html with perceptual coeffs for red and green
    if s == 0.0: return v*red_norm, v*green_norm, v
    i = int(h*6.0) # XXX assume int() truncates!
    f = (h*6.0) - i
    p = v*(1.0 - s)
    q = v*(1.0 - s*f)
    t = v*(1.0 - s*(1.0-f))
    i = i%6
    if i == 0: return np.array([v*red_norm, t*green_norm, p])
    if i == 1: return np.array([q*red_norm, v*green_norm, p])
    if i == 2: return np.array([p*red_norm, v*green_norm, t])
    if i == 3: return np.array([p*red_norm, q*green_norm, v])
    if i == 4: return np.array([t*red_norm, p*green_norm, v])
    if i == 5: return np.array([v*red_norm, p*green_norm, q])

def rgb_palette(n_colors, red_norm=.7, green_norm=.5): # todo stretch hue around orange-yellow-green a bit?
    return np.array([hsv_to_rgb(i,1,1,red_norm, green_norm) for i in np.linspace(0,1-1/n_colors,n_colors)])

def remove_dead_pixels(arr):
    for dp in DEAD_PIXELS:
        arr = np.delete(arr, dp)
    return arr


fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 6))
figl, axl = plt.subplots(nrows=1, ncols=1, figsize=(12, 6))
fige, axe = plt.subplots(nrows=1, ncols=1, figsize=(12, 6))
figle, axle = plt.subplots(nrows=1, ncols=1, figsize=(12, 6))

for fn,c in zip(sys.argv[1:], rgb_palette(len(sys.argv)-1)):
    with open(fn, 'rb') as f:
        label = fn.replace(".RWD8","")
        raw = f.read()
        int_t = struct.unpack(f'f', raw[int('5d',16):int('5d',16)+4])[0] / 1000 # seconds
        if not int_t: continue # probably wrong file?
        print(f"converting & plotting {fn} (integration time = {int_t:.4f} s)") # in hex, it is 0x2148:0x6147

        wl = np.array(struct.unpack(f'{2048}f', raw[int('148',16):int('148',16)+2048*4]))
        y_raw = np.array(struct.unpack(f'{2048}f', raw[8520:16712]))
        y_dark = np.array(struct.unpack(f'{2048}f', raw[16712:24904]))
        wl = remove_dead_pixels(wl)
        y = remove_dead_pixels(y_raw - y_dark) / int_t
        np.savetxt(label+'_nm.txt', np.vstack([wl, y]).T, header='# wavelength_(nm) counts_per_second')
        np.savetxt(label+'_eV.txt', np.vstack([1240/wl, y*wl**2]).T, header='# energy_(eV) counts_per_second')

        convol = 2**-np.linspace(-2,2,25)**2; ysm = np.convolve(y,convol/np.sum(convol), mode='same') ## simple smoothing
        for a in (ax,axl):
            a.plot(wl, ysm, label=label, alpha=.9, lw=1, c=c)
        for a in (axe,axle):
            a.plot(1240/wl, ysm*wl**2 / 400**2, label=label, alpha=.9,lw=1, c=c)

for a in [ax, axe, axl, axle]:
    a.set_ylabel("intensity (a. u.)")
    if a in (ax,axl):
        a.set_xlabel("wavelength (nm)")
    else:
        a.set_xlabel("energy (eV)")
    a.grid()
    a.legend(prop={'size':10}, loc='upper right')
    a.set_ylim(ymin=np.mean(np.abs(y[1:]-y[:-1])) / 10)

fig.savefig("DL_spectra.png", bbox_inches='tight')

fige.savefig("DL_spectra_eV.png", bbox_inches='tight')

axl.set_yscale("log")
figl.savefig("DL_spectra_logarithmic.png", bbox_inches='tight')

axle.set_yscale("log")
figle.savefig("DL_spectra_logarithmic_eV.png", bbox_inches='tight')

""" # tip: visually seek for valid floats in any byte array
for idx in ('5a','5f','5c','5e'):
    print('    ', idx, ':', struct.unpack(f'f', raw[int(idx,16):int(idx,16)+4]))
    print('        ', idx, ':', struct.unpack(f'f', raw[int(idx,16)+1:int(idx,16)+5]))
"""