Put the datafiles alongside the instrument, add remesh tools and a RE…

…ADME.md

mcstas-comps/examples/Tests_polarization/Test_Pol_Tabled/README.md

@@ -0,0 +1,8 @@
+# `Test_Pol_Tabled.instr` with inputs and tools
+The instrument uses tabulated field-data for input, examples are:
+* A flipping field `flipfield.dat` 
+* A constant field `constfield.dat`
+Further, a couple of remeshing tools are provided, allowing to resample an existing field-file to 'regular' binning, both with positional inputs of `infile` for the input file and `xsiz`, `ysiz`, `zsiz` for the resampling dimensions
+* `remesh.m` function for Matlab users
+* `remesh.py` script for Python users
+When resampling for use with 'regular' interpolation (e.g. on GPU), esure that `xsiz`, `ysiz`, `zsiz` are equal. 

mcstas-comps/examples/Tests_polarization/Test_Pol_Tabled/remesh.m

@@ -0,0 +1,59 @@
+function remesh(infile, xsiz, ysiz, zsiz)
+%REMESH Resample magnetic vector field data onto a new grid.
+% REMESH(INFILE, XSIZ, YSIZ, ZSIZ) resamples the magnetic vector field
+% data loaded from INFILE onto a new grid defined by the dimensions XSIZ,
+% YSIZ, and ZSIZ. The resampled data is saved to a new file with a name
+% constructed based on the dimensions and the input filename.
+%
+% Input arguments:
+% - infile: Input file containing magnetic vector field data
+% - xsiz: Number of points along the x-axis in the new grid
+% - ysiz: Number of points along the y-axis in the new grid
+% - zsiz: Number of points along the z-axis in the new grid
+
+% Check if all input arguments are provided
+if nargin < 4
+ error('I need precisely 4 inputs: infile, xsiz, ysiz, zsiz');
+end
+
+% Load data from the input file
+data = load(infile);
+
+% Extract coordinates and magnetic field components
+XYZ = data(:, [1 2 3]); % Coordinates (x, y, z)
+BXYZ = data(:, [4 5 6]); % Magnetic field components (Bx, By, Bz)
+
+% Create a scattered interpolant for interpolation
+BInterp = scatteredInterpolant(XYZ, BXYZ);
+
+% Determine the range of coordinates
+xmin = min(XYZ(:, 1));
+xmax = max(XYZ(:, 1));
+ymin = min(XYZ(:, 2));
+ymax = max(XYZ(:, 2));
+zmin = min(XYZ(:, 3));
+zmax = max(XYZ(:, 3));
+
+% Generate new grid coordinates
+Xnew = linspace(xmin, xmax, xsiz);
+Ynew = linspace(ymin, ymax, ysiz);
+Znew = linspace(zmin, zmax, zsiz);
+
+% Create meshgrid for new coordinates
+[Xmesh, Ymesh, Zmesh] = meshgrid(Xnew, Ynew, Znew);
+
+% Flatten meshgrid arrays
+XYZnew = [Xmesh(:), Ymesh(:), Zmesh(:)];
+
+% Interpolate magnetic field at new points
+BXYZnew = BInterp(XYZnew);
+
+% Construct output filename
+outfile = sprintf('remeshed_%dx%dx%d_%s', xsiz, ysiz, zsiz, infile);
+
+% Combine coordinates and interpolated magnetic field
+Table = [XYZnew, BXYZnew];
+
+% Save data to output file in ASCII format
+save(outfile, 'Table', '-ascii');
+

mcstas-comps/examples/Tests_polarization/Test_Pol_Tabled/remesh.py

@@ -0,0 +1,56 @@
+import argparse
+import numpy as np
+from scipy.interpolate import griddata
+
+def remesh(infile, xsiz, ysiz, zsiz):
+ data = np.loadtxt(infile)
+
+ XYZ = data[:, :3]
+ BXYZ = data[:, 3:]
+
+ xmin, xmax = np.min(XYZ[:,0]), np.max(XYZ[:,0])
+ ymin, ymax = np.min(XYZ[:,1]), np.max(XYZ[:,1])
+ zmin, zmax = np.min(XYZ[:,2]), np.max(XYZ[:,2])
+
+ Xnew = np.linspace(xmin, xmax, xsiz)
+ Ynew = np.linspace(ymin, ymax, ysiz)
+ Znew = np.linspace(zmin, zmax, zsiz)
+
+ Xmesh, Ymesh, Zmesh = np.meshgrid(Xnew, Ynew, Znew, indexing='ij')
+
+ XYZnew = np.column_stack([Xmesh.ravel(), Ymesh.ravel(), Zmesh.ravel()])
+
+ # Interpolate the magnetic field at the new points
+ BInterp = griddata(XYZ, BXYZ, XYZnew, method='linear')
+
+ outfile = f"remeshed_{xsiz}x{ysiz}x{zsiz}_{infile}"
+
+ Table = np.column_stack([XYZnew, BInterp])
+
+ # Calculate the product of xsiz, ysiz, and zsiz as an integer
+ N = int(xsiz * ysiz * zsiz)
+
+ # Add header to the output file
+ header = f"# {N} ({xsiz} x {ysiz} x {zsiz})\n"
+
+ # Format Table as string with scientific notation and 2 decimal places
+ Table_str = "\n".join([" ".join([f"{val:.2e}" for val in row]) for row in Table])
+
+ with open(outfile, 'w') as f:
+ f.write(header)
+ f.write(Table_str)
+
+def main():
+ parser = argparse.ArgumentParser(description='Remesh magnetic vector field data.')
+ parser.add_argument('infile', type=str, help='Input file containing magnetic vector field data')
+ parser.add_argument('xsiz', type=int, help='Number of points along the x-axis in the new grid')
+ parser.add_argument('ysiz', type=int, help='Number of points along the y-axis in the new grid')
+ parser.add_argument('zsiz', type=int, help='Number of points along the z-axis in the new grid')
+
+ args = parser.parse_args()
+
+ remesh(args.infile, args.xsiz, args.ysiz, args.zsiz)
+
+if __name__ == "__main__":
+ main()
+