diff --git a/misc/finiterectlat-constant-driving.py b/misc/finiterectlat-constant-driving.py
index 4d3a0ea..4bf2cb9 100755
--- a/misc/finiterectlat-constant-driving.py
+++ b/misc/finiterectlat-constant-driving.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 import math
-from qpms.argproc import ArgParser, make_dict_action, sslice
+from qpms.argproc import ArgParser, make_dict_action, sslice, annotate_pdf_metadata
 figscale=3
 
 ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax', 'single_omega'])
@@ -44,7 +44,7 @@ import numpy as np
 import qpms
 from qpms.cybspec import BaseSpec
 from qpms.cytmatrices import CTMatrix, TMatrixGenerator
-from qpms.qpms_c import Particle
+from qpms.qpms_c import Particle, qpms_library_version
 from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
 from qpms.cycommon import DebugFlags, dbgmsg_enable
 from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar
@@ -225,7 +225,7 @@ for iri in range(ss.nirreps):
             scattered_full[spi, y] += scattered_ir_unpacked[spi, y]
     if a.save_gradually:
         iriout = outfile_tmp + ".%d" % iri
-        np.savez(iriout, iri=iri, meta=vars(a), 
+        np.savez(iriout, iri=iri, meta={**vars(a), 'qpms_version' : qpms.__version__()}, 
 		 omega=omega, wavenumber=wavenumber, nelem=nelem, wavevector=np.array(a.wavevector), phases=phases, 
                  positions = ss.positions[:,:2],
                  scattered_ir_packed = scattered_ir[iri], 
@@ -251,7 +251,7 @@ if not math.isnan(a.ccd_distance):
     logging.info("Far fields done")
 
 outfile = defaultprefix + ".npz" if a.output is None else a.output
-np.savez(outfile, meta=vars(a), 
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, 
 		 omega=omega, wavenumber=wavenumber, nelem=nelem, wavevector=np.array(a.wavevector), phases=phases, 
                  positions = ss.positions[:,:2],
                  scattered_ir_packed = scattered_ir, 
@@ -355,6 +355,7 @@ if a.plot or (a.plot_out is not None):
                     axes[y,gg].yaxis.set_major_formatter(plt.NullFormatter())
         fig.text(0, 0, str(slicepairs[spi]), horizontalalignment='left', verticalalignment='bottom')
         pp.savefig()
+    annotate_pdf_metadata(pp, scriptname="finiterectlat-constant-driving.py")
     pp.close()
 
 exit(0)
diff --git a/misc/finiterectlat-modes.py b/misc/finiterectlat-modes.py
index 39b946a..701d61d 100755
--- a/misc/finiterectlat-modes.py
+++ b/misc/finiterectlat-modes.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 import math
-from qpms.argproc import ArgParser
+from qpms.argproc import ArgParser, annotate_pdf_metadata
 
 
 ap = ArgParser(['rectlattice2d_finite', 'background_analytical', 'single_particle', 'single_lMax', ])
@@ -100,7 +100,7 @@ results['inside_contour'] = inside_ellipse((results['eigval'].real, results['eig
 results['refractive_index_internal'] = [medium(om).n for om in results['eigval']]
 
 outfile = defaultprefix + (('_ir%s_%s.npz' % (str(iri), irname)) if iri is not None else '.npz') if a.output is None else a.output
-np.savez(outfile, meta=vars(a), **results)
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, **results)
 logging.info("Saved to %s" % outfile)
 
 exit(0)
@@ -110,7 +110,7 @@ if a.plot or (a.plot_out is not None):
     import matplotlib
     matplotlib.use('pdf')
     from matplotlib import pyplot as plt
-
+    from matplotlib.backends.backend_pdf import PdfPages
     fig = plt.figure()
     ax = fig.add_subplot(111)
     ax.plot(sinalpha_list, σ_ext*1e12,label='$\sigma_\mathrm{ext}$')
@@ -119,9 +119,11 @@ if a.plot or (a.plot_out is not None):
     ax.legend()
     ax.set_xlabel('$\sin\\alpha$')
     ax.set_ylabel('$\sigma/\mathrm{\mu m^2}$')
-    
+        
     plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
-    fig.savefig(plotfile)
+    with PdfPages(plotfile) as pdf:
+        pdf.savefig(fig)
+        annotate_pdf_metadata(pdf, scriptname='finiterectlat-modes.py')
 
 exit(0)
 
diff --git a/misc/finiterectlat-scatter.py b/misc/finiterectlat-scatter.py
index 98a3f30..82eae39 100755
--- a/misc/finiterectlat-scatter.py
+++ b/misc/finiterectlat-scatter.py
@@ -111,7 +111,7 @@ for i, omega in enumerate(ap.allomegas):
             σ_scat_arr_ir[i, j, iri] = np.vdot(fi,np.dot(translation_matrix, fi)).real/wavenumber**2
         if a.save_gradually:
             iriout = outfile_tmp + ".%d.%d" % (i, iri)
-            np.savez(iriout, omegai=i, iri=iri, meta=vars(a), omega=omega, k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
+            np.savez(iriout, omegai=i, iri=iri, meta={**vars(a), 'qpms_version' : qpms.__version__()}, omega=omega, k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
                      wavenumber=wavenumber, σ_ext_list_ir=σ_ext_arr_ir[i,:,iri], σ_scat_list_ir=σ_scat_list_ir[i,:,iri])
             logging.info("partial results saved to %s"%iriout)
 
@@ -122,7 +122,8 @@ for i, omega in enumerate(ap.allomegas):
 
 
 outfile = defaultprefix + ".npz" if a.output is None else a.output
-np.savez(outfile, meta=vars(a), k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, 
+        k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
         σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr,
         σ_ext_ir=σ_ext_arr_ir,σ_abs_ir=σ_abs_arr_ir,σ_scat_ir=σ_scat_arr_ir, omega=ap.allomegas, wavenumbers=wavenumbers
        )
@@ -231,6 +232,7 @@ if a.plot or (a.plot_out is not None):
 
             pdf.savefig(fig)
             plt.close(fig)
+        annotate_pdf_metadata(pdf, scriptname="finiterectlat-scatter.py")
 
 
 exit(0)
diff --git a/misc/infiniterectlat-k0realfreqsvd.py b/misc/infiniterectlat-k0realfreqsvd.py
index fb85093..ee9e2da 100755
--- a/misc/infiniterectlat-k0realfreqsvd.py
+++ b/misc/infiniterectlat-k0realfreqsvd.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 import math
-from qpms.argproc import ArgParser
+from qpms.argproc import ArgParser, annotate_pdf_metadata
 
 ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax', 'omega_seq'])
 ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
@@ -93,7 +93,7 @@ for i, omega in enumerate(omegas):
         SVs[iri][i] = np.linalg.svd(ImTW_packed, compute_uv = False)
 
 outfile = defaultprefix + ".npz" if a.output is None else a.output
-np.savez(outfile, meta=vars(a), omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
        )
 logging.info("Saved to %s" % outfile)
 
@@ -102,7 +102,8 @@ if a.plot or (a.plot_out is not None):
     import matplotlib
     matplotlib.use('pdf')
     from matplotlib import pyplot as plt
-
+    from matplotlib.backends.backend_pdf import PdfPages
+    
     fig = plt.figure()
     ax = fig.add_subplot(111)
     cc = plt.rcParams['axes.prop_cycle']()
@@ -117,9 +118,11 @@ if a.plot or (a.plot_out is not None):
     ax.set_xlabel('$\hbar \omega / \mathrm{eV}$')
     ax.set_ylabel('Singular values')
     ax.legend()
-    
+
     plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
-    fig.savefig(plotfile)
+    with PdfPages(plotfile) as pdf:
+        pdf.savefig(fig)
+        annotate_pdf_metadata(pdf, scriptname='infiniterectlat-k0realfreqsvd.py')
 
 exit(0)
 
diff --git a/misc/infiniterectlat-scatter.py b/misc/infiniterectlat-scatter.py
index 9caa0de..116d7b8 100755
--- a/misc/infiniterectlat-scatter.py
+++ b/misc/infiniterectlat-scatter.py
@@ -2,7 +2,7 @@
 
 import math
 pi = math.pi
-from qpms.argproc import ArgParser
+from qpms.argproc import ArgParser, annotate_pdf_metadata
 
 ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax', 'omega_seq_real_ng', 'planewave'])
 ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
@@ -106,7 +106,7 @@ with pgsl_ignore_error(15): # avoid gsl crashing on underflow
 σ_abs_arr = σ_ext_arr - σ_scat_arr
 
 outfile = defaultprefix + ".npz" if a.output is None else a.output
-np.savez(outfile, meta=vars(a), dir_sph=dir_sph_list, k_cart = k_cart_arr, omega = ap.allomegas, E_cart = E_cart_list, wavenumbers= wavenumbers, σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr, unitcell_area=ss.unitcell_volume
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, dir_sph=dir_sph_list, k_cart = k_cart_arr, omega = ap.allomegas, E_cart = E_cart_list, wavenumbers= wavenumbers, σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr, unitcell_area=ss.unitcell_volume
        )
 logging.info("Saved to %s" % outfile)
 
@@ -214,6 +214,6 @@ if a.plot or (a.plot_out is not None):
             
             pdf.savefig(fig)
             plt.close(fig)
-
+        annotate_pdf_metadata(pdf)
 exit(0)
 
diff --git a/misc/lat2d_modes.py b/misc/lat2d_modes.py
index f7dad4b..8eea75b 100755
--- a/misc/lat2d_modes.py
+++ b/misc/lat2d_modes.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 import math
-from qpms.argproc import ArgParser, sfloat
+from qpms.argproc import ArgParser, sfloat, annotate_pdf_metadata
 
 ap = ArgParser(['const_real_background', 'lattice2d', 'multi_particle']) # TODO general analytical background
     
@@ -112,7 +112,7 @@ res['inside_contour'] = inside_ellipse((res['eigval'].real, res['eigval'].imag),
 #del res['omega']  If contour points are not needed...
 #del res['ImTW'] # not if dbg=false anyway
 outfile = defaultprefix + ".npz" if a.output is None else a.output
-np.savez(outfile, meta=vars(a), empty_freqs=np.array(empty_freqs), 
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, empty_freqs=np.array(empty_freqs), 
         ss_positions=ss.positions, ss_fullvec_poffsets=ss.fullvec_poffsets, 
         ss_fullvec_psizes=ss.fullvec_psizes,
         ss_bspecs_flat = np.concatenate(ss.bspecs),
@@ -132,6 +132,7 @@ if a.plot or (a.plot_out is not None):
     import matplotlib
     matplotlib.use('pdf')
     from matplotlib import pyplot as plt
+    from matplotlib.backends.backend_pdf import PdfPages
 
     fig = plt.figure()
     ax = fig.add_subplot(111,)
@@ -155,8 +156,11 @@ if a.plot or (a.plot_out is not None):
     ax.set_ylim([ymin-.1*yspan, ymax+.1*yspan])
     ax.set_xlabel('$\hbar \Re \omega / \mathrm{eV}$')
     ax.set_ylabel('$\hbar \Im \omega / \mathrm{meV}$')
+
     plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
-    fig.savefig(plotfile)
+    with PdfPages(plotfile) as pdf:
+        pdf.savefig(fig)
+        annotate_pdf_metadata(pdf, scriptname='lat2d_modes.py')
 
 exit(0)
 
diff --git a/misc/lat2d_realfreqsvd.py b/misc/lat2d_realfreqsvd.py
index 90551b7..b955575 100755
--- a/misc/lat2d_realfreqsvd.py
+++ b/misc/lat2d_realfreqsvd.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 import math
-from qpms.argproc import ArgParser, sfloat
+from qpms.argproc import ArgParser, sfloat, annotate_pdf_metadata
 
 ap = ArgParser(['background', 'lattice2d', 'multi_particle',  'omega_seq'])
 
@@ -88,7 +88,7 @@ for i, omega in enumerate(omegas):
         SVs[iri][i] = np.linalg.svd(ImTW_packed, compute_uv = False)
 
 outfile = defaultprefix + ".npz" if a.output is None else a.output
-np.savez(outfile, meta=vars(a), omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
        )
 logging.info("Saved to %s" % outfile)
 
@@ -97,6 +97,7 @@ if a.plot or (a.plot_out is not None):
     import matplotlib
     matplotlib.use('pdf')
     from matplotlib import pyplot as plt
+    from matplotlib.backends.backend_pdf import PdfPages
 
     fig = plt.figure()
     ax = fig.add_subplot(111)
@@ -118,9 +119,11 @@ if a.plot or (a.plot_out is not None):
     ax.set_xlabel('$\hbar \omega / \mathrm{eV}$')
     ax.set_ylabel('Singular values')
     ax.legend()
-    
+
     plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
-    fig.savefig(plotfile)
+    with PdfPages(plotfile) as pdf:
+        pdf.savefig(fig)
+        annotate_pdf_metadata(pdf, scriptname='lat2d_realfreqsvd.py')
 
 exit(0)
 
diff --git a/misc/rectlat_simple_modes.py b/misc/rectlat_simple_modes.py
index 1f706fa..b39b5e4 100755
--- a/misc/rectlat_simple_modes.py
+++ b/misc/rectlat_simple_modes.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 import math
-from qpms.argproc import ArgParser
+from qpms.argproc import ArgParser, annotate_pdf_metadata
 
 ap = ArgParser(['rectlattice2d', 'const_real_background', 'single_particle', 'single_lMax']) # const_real_background needed for calculation of the diffracted orders
 ap.add_argument("-k", nargs=2, type=float, required=True, help='k vector', metavar=('K_X', 'K_Y'))
@@ -116,7 +116,7 @@ res['refractive_index_internal'] = [emg(om).n for om in res['eigval']]
 #del res['omega']  If contour points are not needed...
 #del res['ImTW'] # not if dbg=false anyway
 outfile = defaultprefix + ".npz" if a.output is None else a.output
-np.savez(outfile, meta=vars(a), empty_freqs=np.array(empty_freqs), 
+np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, empty_freqs=np.array(empty_freqs), 
         ss_positions=ss.positions, ss_fullvec_poffsets=ss.fullvec_poffsets, 
         ss_fullvec_psizes=ss.fullvec_psizes,
         ss_bspecs_flat = np.concatenate(ss.bspecs),
@@ -133,6 +133,7 @@ if a.plot or (a.plot_out is not None):
     import matplotlib
     matplotlib.use('pdf')
     from matplotlib import pyplot as plt
+    from matplotlib.backends.backend_pdf import PdfPages
 
     fig = plt.figure()
     ax = fig.add_subplot(111,)
@@ -156,8 +157,11 @@ if a.plot or (a.plot_out is not None):
     ax.set_ylim([ymin-.1*yspan, ymax+.1*yspan])
     ax.set_xlabel('$\hbar \Re \omega / \mathrm{eV}$')
     ax.set_ylabel('$\hbar \Im \omega / \mathrm{meV}$')
+
     plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
-    fig.savefig(plotfile)
+    with PdfPages(plotfile) as pdf:
+        pdf.savefig(fig)
+        annotate_pdf_metadata(pdf, scriptname="rectlat_simple_modes.py")
 
 exit(0)