diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp
index 8dd18b5db9..e0be6f3c15 100644
--- a/mcstas-comps/contrib/Phonon_BvK_PG.comp
+++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp
@@ -38,7 +38,7 @@
 * yheight: [m]        Height of sample in y direction
 * xwidth: [m]         Horiz. dimension of sample if slap
 * zdepth: [m]         Depth of sample if slap
-
+*
 * sigma_abs: [barns]  Absorption cross section at 2200 m/s per atom
 * sigma_inc: [barns]  Incoherent scattering cross section per atom
 * DW: [1]             Debye-Waller factor; scalar value (no q-dependence)
@@ -71,7 +71,6 @@ DEFINE COMPONENT Phonon_BvK_PG
 SETTING PARAMETERS (hh=0,kk=0,ll=0,radius=0, xwidth=0, yheight=0, zdepth=0, sigma_abs=0,sigma_inc=0,DW=1,T,
 focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0,target_x=0, target_y=0, target_z=0, int target_index=0, int mode_input, int e_steps_low, int e_steps_high, int verbose_input=0, int dispersion=0)
 NOACC
-// The component is currently "NOACC" only.
 
 SHARE
 %{
@@ -1121,12 +1120,6 @@ MCDISPLAY
   } else if (xwidth && yheight) { /* box/rectangle */
     box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0);
   }
-  //  circle("xz", 0,  yheight/2.0, 0, radius);
-  //  circle("xz", 0, -yheight/2.0, 0, radius);
-  //  line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
-  //  line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
-  //  line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
-  //  line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
 %}
 
 END
diff --git a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr
new file mode 100644
index 0000000000..6cc880f468
--- /dev/null
+++ b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr
@@ -0,0 +1,315 @@
+/*******************************************************************************
+* Instrument: Test_Phonon_BvK
+*
+* %I
+* Written by: Kim Lefmann
+* Date: 20.07.2025
+* Origin: NBI, KU
+* %INSTRUMENT_SITE: Generic
+*
+* Simple thermal triple axis spectrometer for test of Phonon_BvK_PG component
+* Modified from older version by Dorothy Wang, Kim Lefmann
+*
+* %D
+* Simple TAS to test the Phonon_BvK_PG component
+*
+* %Example: -y Detector: e_mon_beforeana_zoom_I=6.90143e-16
+*
+* %P
+* E:       [meV] Energy transfer
+* Ef:      [meV] Final energy
+* Dlambda: [AA] width of wavelength band
+* h:       [rlu] q-value in hexagonal plane
+* l:       [rlu] q-value out of hexagonal plane
+* dA3:     [deg] A3 offset
+* phononmode: [1] choice of a specific phonon mode  
+*
+* %L
+* <reference/HTML link>
+*
+* %E
+*******************************************************************************/
+DEFINE INSTRUMENT Test_Phonon_BvK_PG(E=4.5,Ef=25,Dlambda=0.05,h=0,l=3.5,dA3=90, Temp=2, width = 0.001, coll = 20, int phononmode=0,int E_steps_high=100, int E_steps_low=100, int Verbose=0, int DISP=0)
+
+//SHELL " cc -c eigen.c "
+//DEPENDENCY " eigen.o "
+
+DECLARE
+%{
+  double Gqx,Gqy,Gqz;
+//double Ef=24.8
+double Ei;
+//meV
+double qx,qy,qz;
+double thetaM;
+double twothetaS;
+double thetaA;
+double A3;
+//double pi;
+double QM;
+double alpha;
+double lambda_i;
+double SMALL__NUMBER;
+%}
+
+INITIALIZE
+%{
+// Set monochromator/analyzer Q-value for PG
+QM = 1.8734;
+SMALL__NUMBER = 1e-6;
+//pi = 3.14159265;
+
+double a = 2.461;
+double c =  6.708;
+//PG lattice constants, in AA, same as in Phonon_BvK_PG.comp
+
+double astar = 4*pi/(sqrt(3)*a);
+double cstar = 2*pi/c;
+//PG reciprocal lattice vectors, in 1/AA
+
+//calculate Ei
+Ei=Ef+E;
+
+//calculate ki, kf, lambda_i, q
+double ki = V2K*SE2V*sqrt(Ei);
+double kf = V2K*SE2V*sqrt(Ef);
+lambda_i=2*pi/ki;
+double q = sqrt(h*h*astar*astar+l*l*cstar*cstar);
+
+//calculate 2thetaM and 2thetaA
+thetaM = asin(QM/(2*ki))*RAD2DEG;
+thetaA = asin(QM/(2*kf))*RAD2DEG;
+
+//calculate scattering angle and sample rotation
+twothetaS = acos((q*q-ki*ki-kf*kf)/(-2*ki*kf))*RAD2DEG;
+alpha = acos((kf*kf-ki*ki-q*q)/(-2*ki*q));
+A3=(asin(l*cstar/(q+SMALL__NUMBER))-alpha)*RAD2DEG;
+
+//Printout calculations
+printf("a = %g c = %g astar = %g cstar = %g \n",a,c,astar,cstar);
+printf("ki = %g kf = %g q = %g lambda_i = %g\n",ki,kf,q,lambda_i);
+printf("thetaA = %g thetaM = %g\n",thetaA,thetaM);
+printf("twothetaS = %g \n",twothetaS);
+printf("A3 = %g \n",A3);
+printf("alpha = %g \n",alpha);
+
+%}
+
+TRACE
+
+COMPONENT origin = Progress_bar()
+AT (0, 0, 0) RELATIVE ABSOLUTE
+
+COMPONENT source = Source_Maxwell_3(
+    xwidth=width,
+    yheight=width,
+    Lmin=lambda_i-Dlambda/2, 
+    Lmax=lambda_i+Dlambda/2, 
+    dist=7.5, 
+    focus_yh=width, 
+    focus_xw=width, 
+    T1=300,     T2=300,     T3=300, 
+    I1=1E15,    I2=1E15,    I3=1E15)
+AT (0, 0, 0) RELATIVE PREVIOUS
+
+COMPONENT l_monitor_before_mono = L_monitor(
+    nL=200, 
+    filename="L_before_mono", 
+    Lmin=0, 
+    Lmax=4,   
+    xwidth=0.16, 
+    yheight=0.25, 
+    restore_neutron=1)
+AT (0, 0, 7) RELATIVE PREVIOUS
+
+COMPONENT monochromator_flat = Monochromator_flat(
+     mosaicv=30,
+     mosaich=30,
+     zwidth = width,
+     yheight = width,
+     Q=QM)
+AT (0, 0, 7.5) RELATIVE source
+ROTATED (0, thetaM, 0) RELATIVE source
+
+COMPONENT arm1 = Arm()
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, thetaM, 0) RELATIVE PREVIOUS
+
+COMPONENT collimator_linear1 = Collimator_linear(
+    xwidth=0.1, 
+    yheight=0.2, 
+    length=0.2, 
+    divergence=coll)
+AT (0, 0, 0.6) RELATIVE arm1
+
+COMPONENT l_monitor_before_sample_big = L_monitor(
+    nL=300, 
+    filename="L_before_sample", 
+    Lmin=lambda_i - Dlambda, 
+    Lmax=lambda_i + Dlambda,
+    xwidth = 0.1,
+    yheight = 0.1,
+    restore_neutron=1)
+AT (0, 0, 0.9) RELATIVE arm1
+
+COMPONENT e_monitor_before_sample_big = E_monitor(
+    nE=100, 
+    filename="E_before_sample", 
+    Emin=Ei-3, 
+    Emax=Ei+3, 
+    xwidth = 0.1,
+    yheight = 0.1,
+   restore_neutron=1)
+AT (0, 0, 0.01) RELATIVE PREVIOUS
+
+COMPONENT PSD_monitor_before_sample = PSD_monitor(
+    nx = 200,
+    ny = 200,
+    xwidth = width/2,
+    yheight = width/2,
+    filename="PSD_before_sample", 
+    restore_neutron=1)
+AT (0, 0, 0.998) RELATIVE arm1
+
+COMPONENT arm2 = Arm()
+AT (0, 0, 1) RELATIVE arm1
+ROTATED (0, -A3+dA3, 0) RELATIVE arm1
+
+SPLIT 10
+/*COMPONENT incoherent = Incoherent( 
+    yheight=width/2, 
+    xwidth = width/2,
+   zdepth = width/2000,
+    focus_xw=width, 
+    focus_yh=width, 
+    target_index=3,
+   order = 1, 
+    sigma_abs=0) */
+COMPONENT phonon_bvk_pg = Phonon_BvK_PG(
+//COMPONENT phonon_bvk_pg = Phonon_BvK_PG_eigenvector_Dec2024(
+//    radius=width/2,
+    xwidth=width/2, zdepth=width/20,
+    yheight=width/2, 
+    sigma_abs=0, 
+    sigma_inc=0, 
+    DW=1, 
+    T=Temp, 
+    mode_input=phononmode,
+    target_index=3,
+    focus_xw=width,
+    focus_yh=width,
+    e_steps_low = E_steps_low,
+    e_steps_high = E_steps_high,
+    verbose_input=Verbose,
+    hh = h, kk = 0, ll = l, dispersion = DISP
+)
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, 0, 0) RELATIVE PREVIOUS  // To use the (0,k,l) scattering plane, rotate (0, 0, -60); for (h,k,0) rotate (90, 0, 0)
+
+COMPONENT arm3 = Arm()
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, -twothetaS, 0) RELATIVE arm1
+
+/*COMPONENT PSD_sample_image = PSD_monitor(
+    nx=200, 
+    ny=200, 
+   xwidth = 0.01,
+   yheight = 0.01,
+    filename="PSD_sample_image", 
+    restore_neutron=1)
+AT (0, 0, 1.005) RELATIVE arm1*/
+
+/*COMPONENT collimator_linear2 = Collimator_linear(
+    xwidth=0.1, 
+    yheight=0.2, 
+    length=0.2, 
+    divergence=coll)
+AT (0, 0, 0.5) RELATIVE PREVIOUS
+*/
+
+COMPONENT PSD_monitor_after_sample = PSD_monitor(
+    nx=200, 
+    ny=200, 
+    filename="PSD_after_sample", 
+    restore_neutron=1)
+AT (0, 0, 0.9999) RELATIVE arm3
+
+COMPONENT slit1 = Slit(
+    xwidth=width, 
+    yheight=width)
+AT (0, 0, 1) RELATIVE arm3
+
+COMPONENT e_monitor_beforeana = E_monitor(
+    nE=750,
+    filename="E_before_ana", 
+    xwidth=2*width, 
+    yheight=2*width, 
+    Emin=0, 
+    Emax=75)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+
+COMPONENT e_mon_beforeana_zoom = E_monitor(
+    nE=750,
+    filename="E_before_ana_zoom", 
+    xwidth=2*width, 
+    yheight=2*width, 
+    Emin=Ef*0.9, 
+    Emax=Ef*1.1)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+COMPONENT analyzer = Monochromator_flat(
+    zwidth=4*width,
+    yheight=4*width,
+    mosaicv=30,
+    mosaich=30,
+    Q=QM
+    )
+AT (0, 0, 0.1) RELATIVE PREVIOUS
+ROTATED (0, thetaA, 0) RELATIVE PREVIOUS
+
+COMPONENT arm4 = Arm()
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, thetaA, 0) RELATIVE PREVIOUS
+
+COMPONENT slit2 = Slit(
+    xwidth=width,
+    yheight=width)
+AT (0, 0, 1) RELATIVE PREVIOUS
+EXTEND
+%{
+//    qx = MC_GETPAR(phonon_bvk_pg,q_x);
+//    qy = MC_GETPAR(phonon_bvk_bg,q_y);
+//    qz = MC_GETPAR(phonon_bvk_bg,q_z);
+//    printf("hit; q = (%g, %g, %g)\n",0,0,0);
+%}
+
+COMPONENT Emon_after_analyzer = E_monitor(
+    nE=750,
+    filename="E_after_analyzer", 
+    xwidth=0.05,
+    yheight=0.10, 
+    Emin=0, 
+    Emax=250)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+
+COMPONENT Emon_1storder = E_monitor(
+    nE=50,
+    filename="E_1st_order", 
+    xwidth=width,
+    yheight=width, 
+    Emin=Ef*0.9, 
+    Emax=Ef*1.1)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+
+/*COMPONENT PSD_detector = PSD_monitor(
+    nx=128, 
+    ny=128, 
+    filename="PSD_detector", 
+    xwidth=0.05,
+    yheight=0.10)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+*/
+FINALLY
+%{
+%}
+
+END
diff --git a/mcstas-comps/share/eigen.c b/mcstas-comps/share/eigen.c
new file mode 100644
index 0000000000..67957f47c2
--- /dev/null
+++ b/mcstas-comps/share/eigen.c
@@ -0,0 +1,472 @@
+/************************************************************
+ * McStas Eigensolver library                               *
+ * Contributed by James Emil Avery, 202x                    *
+ * Department of Computer Science, University of Copenhagen *
+ ************************************************************/
+
+/* Original file "eigen.c": */
+
+#ifndef EIGEN_LIB
+#define EIGEN_LIB
+
+#include <stddef.h>
+#include <float.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <stdio.h>
+#include <complex.h>
+#include <time.h>
+
+/* TODO: Tridiagonal input/outpout in (3,m) storage. */
+/* TODO: Symmetric tridiagonal input/outpout in (2,m) storage. */
+/* TODO: Add pivot to handle ill-conditioned matrices */
+// #include "eigen.h" 
+
+#define G PRINTF_G
+#define ROWS 0
+#define COLS 1
+
+/* TOOD: Generic matrix functions to separate file */
+INLINE real_t vector_norm2(const cdouble*x, int n){
+  real_t sum = 0;
+  for(int i=0;i<n;i++) sum += cabs(cmul(x[i],conj(x[i])));
+  return sum;
+}
+
+INLINE real_t vector_norm(const cdouble *x, int n){
+  return sqrt(vector_norm2(x,n));
+}
+
+INLINE void extract_region(cdouble *S, int N,
+			   int i0, int j0, int m, int n,
+			   cdouble *D)
+{
+  for(int i=0;i<m;i++)
+    for(int j=0;j<n;j++)
+      D[i*n+j] = S[(i+i0)*N+(j+j0)];
+}
+
+INLINE void identity_matrix(cdouble *Q, int n)
+{
+  memset(Q,0,n*n*sizeof(cdouble));
+  for(int i=0;i<n;i++) Q[i*n+i] = cplx(1,0);
+}
+
+INLINE void real_identity_matrix(real_t *Q, int n)
+{
+  memset(Q,0,n*n*sizeof(real_t));
+  for(int i=0;i<n;i++) Q[i*n+i] = 1;
+}
+
+/* TODO: Make this work with both real and complex cdoubles (cabs should be replaced by a macro) */
+INLINE real_t max_norm(cdouble *A, int m, int n)
+{
+  real_t mx = 0;
+  for(int i=0;i<m;i++){ 	
+    real_t row_norm = 0;
+    for(int j=0;j<n;j++) row_norm += cabs(A[i*n+j]);
+    mx = fmax(mx,row_norm);
+  }
+  return mx;  
+}
+
+INLINE void matrix_inplace_multiply(cdouble *A, cdouble *B, int m, int n, int q)
+{
+  for(int i=0;i<m;i++){
+    cdouble* row=malloc(q*sizeof(cdouble));
+    for(int j=0;j<q;j++){
+      cdouble sum = cplx(0,0);
+      for(int k=0;k<n;k++) sum = cadd(sum,cmul(A[i*n+k],B[k*q+j]));
+      row[j] = sum;
+    }
+    for(int j=0;j<q;j++) A[i*n+j] = row[j];
+  }
+}
+
+
+/* TODO: Document */
+INLINE void apply_reflection(/*in/out*/cdouble* A, const cdouble *v, int m, int n, cdouble sigma, int transpose)
+{
+  cdouble* vHA=malloc(n*sizeof(cdouble));
+  cdouble conj_sigma = conj(sigma);
+  int stride[2] = {n*(!transpose) + 1*transpose, n*transpose + 1*(!transpose)};
+  
+  for(int j=0;j<n;j++){		
+    cdouble sum = cplx(0,0);
+    for(int k=0;k<m;k++) sum = cadd(sum,cmul(conj(v[k]),A[k*stride[0]+j*stride[1]]));
+    vHA[j] = sum;
+  }
+
+  for(size_t i=0;i<m;i++)       /* A += -2*outer(v,vTA) */
+    for(size_t j=0;j<n;j++){
+      A[i*stride[0]+j*stride[1]] = csub(A[i*stride[0]+j*stride[1]],cmul(cmul(conj_sigma,v[i]),vHA[j])); 
+    }  
+}
+
+/* A: NxN
+   Transform mxn region starting at (i0,j0)
+ */
+INLINE void reflect_region(/*in/out*/cdouble* A, int N, int i0, int j0, int m, int n, const cdouble *v, cdouble sigma, int cols)
+{
+  cdouble* vHA=malloc(n*sizeof(cdouble));
+  cdouble conj_sigma = conj(sigma);
+  int stride[2] = {N*(!cols) + 1*cols, N*cols + 1*(!cols)};
+
+  /* printf("reflect_region(A,%d,\n"
+  	 "                 %d,%d\n"
+  	 "                 %d,%d)\n",N,i0,j0,m,n);
+  printf("stride = [%d,%d]\n",stride[0],stride[1]);
+   */
+  for(int j=0;j<n;j++){		
+    cdouble sum = cplx(0,0);
+    for(int i=0;i<m;i++){
+      size_t IJ = (i+i0)*stride[0] + (j+j0)*stride[1];
+      sum = cadd(sum,cmul(conj(v[i]),A[IJ]));
+    }
+    vHA[j] = sum;
+  }
+
+  for(size_t i=0;i<m;i++)       /* A += -2*outer(v,vTA) */
+    for(size_t j=0;j<n;j++){
+      size_t IJ = (i+i0)*stride[0] + (j+j0)*stride[1];
+      A[IJ] = csub(A[IJ],cmul(conj_sigma,cmul(v[i],vHA[j]))); 
+    }  
+}
+
+/* Apply n 2D reflections to an mxm matrix Q */
+/* TODO: Transform rows instead (after this works) */
+/* TODO: Since we no longer have views, pass view as function parameter */
+void apply_real_reflections(real_t* V, int n, real_t* Q, int m)
+{
+  if(Q != 0){       // Do we want eigenvectors?
+    for(int k=0;k<n;k++){
+      real_t v0 = V[2*k], v1 = V[2*k+1];      
+      // Udrullet:
+      //       apply_reflection(Q({k,k+2},{0,m}), v);
+      for(int l=0;l<m;l++){
+	real_t q0 = Q[k*m+l], q1 = Q[(k+1)*m+l]; /* cdouble *q = &Q[l*m+k] */
+	real_t vTA = q0*v0 + q1*v1;
+	Q[k*m+l]     -= 2*v0*vTA;
+	Q[(k+1)*m+l] -= 2*v1*vTA;
+      }      
+    }  
+  }  
+}
+
+/* INLINE real_t COPYSIGN(real_t to, real_t from) */
+/* { */
+/*   if(fabs(from/to)<1e-14) return to; */
+/*   else return copysign(to,from); */
+/* } */
+
+INLINE void reflection_vector(/*in*/cdouble* a, const real_t anorm,
+			      /*out*/cdouble* v, cdouble* sigma, int n)
+{ // Reflection vector that eliminates *row* a (as opposed to *column*)
+
+  for(int i=0;i<n;i++) v[i] = a[i];
+
+  // In the complex case, straightforward Householder reflection can only transform a Hermitian matrix
+  // to a complex Hermitian tridiagonal matrix.
+  // The following transform reduces to Householder reflection for real matrices, but
+  // in the complex case, it produces a real symmetric tridiagonal matrix.
+  
+  // Using the nomenclature from [LAWN72], Pages 4-5 and setting xi_i = a[i-1]:
+  cdouble x1 = a[0];
+  real_t nu = copysign(anorm,creal(x1)); 
+  cdouble norm_inv = cdiv(cplx(1.0,0),radd(nu,x1));
+  *sigma = rmul(1/nu,radd(nu,x1));
+  v[0] = radd(nu,v[0]);
+  
+  for(size_t i=0;i<n;i++) v[i] = cmul(v[i],norm_inv);
+}
+
+
+
+// Decompose a matrix (complex or real) into A = Q H Q.H(), where H is upper Hessenberg.
+// If A is Hermitian (symmetric in real case), then H is tridiagonal.
+// TODO: Row pivot
+void print_vector(const char* name, cdouble* a, int l)
+{
+  printf("%s = array([",name); for(int i=0;i<l;i++) printf("%.3g + %.3gj%s",creal(a[i]), cimag(a[i]), i+1<l?", ":"])\n");
+}
+
+void print_matrix(const char *name, cdouble *A, int m, int n)
+{
+  printf("%s = array([\n",name);
+  for(int i=0;i<m;i++){
+    printf("\t[");
+    for(int j=0;j<n;j++)
+      printf("%.3g + %.3gj%s",creal(A[i*n+j]), cimag(A[i*n+j]), j+1<n?", ":"]");
+    printf("%s",i+1<m?",\n":"\n])\n");
+  }
+}
+
+/* TODO: Kan jeg håndtere a[0] =~= 0 mere robust? Kan jeg inkludere pivots? */
+void QHQ(/*in/out*/cdouble* A, int n, cdouble* Q)
+{
+  cdouble  sigma;	        // Elementary operation scale (2 for reflection)
+  cdouble* v=malloc(n*sizeof(cdouble));
+  cdouble* vc=malloc(n*sizeof(cdouble));
+  cdouble* a=malloc(n*sizeof(cdouble));    // Reflection vector
+
+  real_t numerical_zero = max_norm(A,n,n)*10*machine_precision;
+  
+  for(int k=0;k<n-1;k++){
+    int l = n-(k+1);		     /* length of kth subdiagonal column */
+    extract_region(A,n,k+1,k,l,1,a); /* a = A[k:(k+1),(k+1):n], kth row postdiagonal. */    
+    real_t anorm = vector_norm(a,l);
+/*     printf("\n\n%d: ",k); print_vector("a",a,l);
+    printf("%d: ",k); print_matrix("A before",A,n,n);    
+ */
+    if(anorm < numerical_zero) continue; /* Already eliminated, don't divide by 0 */
+    
+    reflection_vector(a,anorm,v,&sigma,n);    /* Vector definining elimination operations */
+    for(int i=0;i<l;i++) vc[i] = conj(v[i]);  /* ...and its conjugate. */
+    
+    reflect_region(A,n,/* reg. start:*/k+1,k, /* reg. size:*/l,l+1, v,       sigma,  ROWS);/* TODO: fix */
+    reflect_region(A,n,/* reg. start:*/k+1,k, /* reg. size:*/l,l+1, vc,conj(sigma),  COLS);
+    // printf("%d: ",k); print_matrix("A after",A,n,n);
+
+    if(Q != 0) reflect_region(Q,n, k+1,0, l,n, v, sigma, ROWS); 
+  }
+  free(v);free(vc);free(a);
+}
+
+// TODO: T_QTQ based on Givens rotations (should be possible to do with fewer operations)
+//size_t QTQ_calls = 0;
+void T_QTQ(const int n, const real_t *Din, const real_t *Lin, real_t *Dout, real_t *Lout, real_t *Vout, real_t shift)
+{
+  // Unrolled
+  //  real_t numerical_zero = T.max_norm()*10*machine_precision;
+  // specialized max_norm = max(sum(abs(A),axis=1)) for tridiagonal matrix. 
+  real_t max_norm = 0, numerical_zero = 0;
+  //  for(int i=0;i<n;i++) max_norm = max(max_norm, fabs(Din[i]) + 2*fabs(Lin[i]));
+  //  numerical_zero = 10*max_norm*machine_precision;//TODO: max_norm for symmetric tridiagonal
+  numerical_zero = 100*machine_precision;
+  
+  real_t a[2], v[2];
+  real_t *D=malloc((n+1)*sizeof(real_t));
+  real_t *L=malloc((n+1)*sizeof(real_t));
+  real_t *U=malloc((2*(n+1))*sizeof(real_t));
+
+  for(int i=0;i<n+1;i++){
+    D[i] = Din[i]-shift;		// Diagonal
+    L[i] = 0;			// Zero padding to avoid branching in inner loop
+    U[i] = 0;                   // Zero padding to avoid branching in inner loop
+    U[(n+1)+i] = 0;		// Second upper diagonal for fill-in. U[n+k] = T(k,k+2) is the element two rows above (k+2)st diagonal element.
+    if(i<n-1){
+      L[ i ] = Lin[i];	// First lower subdiagonal. L[k] = T(k+1,k) is element below kth diagonal element.
+      U[ i ] = Lin[i];	// First upper subdiagonal. U[k] = T(k,k+1) is element above (k+1)st diagonal element.
+      Vout[2*i] = 0; Vout[2*i+1] = 0;	// i'th reflection vector. (0,0) yields "no reflection". Must be initialized due to skipped steps.          
+    }
+  }
+
+  for(int k=0;k<n-1;k++)
+    if(fabs(L[k]) > numerical_zero)  // Only process if subdiagonal element is not already zero.
+    {
+      a[0] = D[k]; a[1] = L[k];       // a = T[k:k+2,k] is the vector of nonzeros in kth subdiagonal column.
+      
+      real_t anorm = sqrt(a[0]*a[0] + a[1]*a[1]); 
+
+      // // Udrullet
+      // //    reflection_vector(a,anorm,v);
+      v[0] = D[k]; v[1] = L[k];
+      real_t alpha = -copysign(anorm,a[0]); // Koster ingenting
+      v[0] -= alpha;
+
+      real_t vnorm = sqrt(v[0]*v[0]+v[1]*v[1]);
+      real_t norm_inv = 1/vnorm;               /* Normalize */
+      v[0] *= norm_inv;  v[1] *= norm_inv;
+
+      Vout[2*k] = v[0]; Vout[2*k+1] = v[1];
+      
+      // // Udrullet 
+      // //    apply_reflection(T({k,k+2},{k,k+3}),v);
+      // //      if(k+1<n){			// k=n-1 case handled by padding with zeros
+      real_t   vTA[3] = {D[ k ]*v[0] + L[ k ]*v[1],  // T(k,k  )*v[0] + T(k+1,k  )*v[1]
+      			 U[ k ]*v[0] + D[k+1]*v[1],  // T(k,k+1)*v[0] + T(k+1,k+1)*v[1]
+      			 U[(n+1)+k]*v[0] + U[k+1]*v[1]}; // T(k,k+2)*v[0] + T(k+1,k+2)*v[1]
+
+      D[ k ]     -= 2*v[0]*vTA[0];
+      L[ k ]     -= 2*v[1]*vTA[0];
+      U[ k ]     -= 2*v[0]*vTA[1];
+      D[k+1]     -= 2*v[1]*vTA[1];
+      U[(n+1)+k] -= 2*v[0]*vTA[2];
+      U[k+1]     -= 2*v[1]*vTA[2];
+    }
+
+  // Transform from the right = transform columns of the transpose.
+  {
+    int k = 0;
+    const real_t *v = &Vout[0];
+    real_t   vTA[2] = {D[ k ]*v[0] + U[  k  ]*v[1],  // T(k,k  )*v[0] + T(k,  k+1)*v[1]
+  		          0        + D[ k+1 ]*v[1]}; // T(k+1,k)*v[0] + T(k+1,k+1)*v[1]. Lower subdiagonal is zero at this stage.
+    
+    D[k]       -= 2*v[0]*vTA[0]; // T(k,k)     -= 2*v[0]*vTA[0]
+    U[k]       -= 2*v[1]*vTA[0]; // T(k,k+1)   -= 2*v[1]*vTA[0]
+    L[k]       -= 2*v[0]*vTA[1]; // T(k+1,k)   -= 2*v[0]*vTA[1]
+    D[k+1]     -= 2*v[1]*vTA[1]; // T(k+1,k+1) -= 2*v[1]*vTA[1]        
+  }    
+
+  for(int k=1;k<n-1;k++){
+    const real_t *v = &Vout[2*k];
+
+    real_t   vTA[3] = {U[k-1]*v[0] + U[(n+1)+k-1]*v[1], // T(k-1,k)*v[0] + T(k-1,k+1)*v[1]  
+  		       D[ k ]*v[0] + U[  k  ]*v[1],     // T(k,k  )*v[0] + T(k,  k+1)*v[1]
+  		       L[ k ]*v[0] + D[ k+1 ]*v[1]};    // T(k+1,k)*v[0] + T(k+1,k+1)*v[1]. Lower subdiagonal is zero at this stage
+
+    U[k-1]     -= 2*v[0]*vTA[0];     // T(k-1,k)   -= 2*v[0]*vTA[0]
+    U[(n+1)+(k-1)] -= 2*v[1]*vTA[0]; // T(k-1,k+1) -= 2*v[1]*vTA[0]
+    D[k]       -= 2*v[0]*vTA[1];     // T(k,  k)     -= 2*v[0]*vTA[1]
+    U[k]       -= 2*v[1]*vTA[1];     // T(k,  k+1)   -= 2*v[1]*vTA[1]
+    L[k]       -= 2*v[0]*vTA[2];     // T(k+1,k)   -= 2*v[0]*vTA[2]
+    D[k+1]     -= 2*v[1]*vTA[2];     // T(k+1,k+1) -= 2*v[1]*vTA[2]        
+  } 
+  
+  // Copy working diagonals to output
+  for(int i=0;i<n;i++){
+    Dout[i] = D[i] + shift;	  // Diagonal
+    if(i<n-1){
+      Lout[i] = U[i];  // First lower subdiagonal. L[k] = T(k+1,k) is element below kth diagonal element.
+    }
+  }
+  free(D); free(L); free(U);
+}
+
+
+INLINE real_pair eigvalsh2x2(real_t a, real_t b, real_t c, real_t d){
+  real_pair result;
+  real_t D = sqrt(4*b*c+(a-d)*(a-d));
+
+  result.value[0] = (a+d-D)/2;
+  result.value[1] = (a+d+D)/2;
+
+  return result;
+}
+
+
+int    nth_time = 0;
+real_t eigensystem_cputime = 0;
+#define gershgorin_tolerance 1e4*machine_precision
+#define max_iterations       30
+
+// TODO: Til tridiagonale matricer er Givens-rotation nemmere/hurtigere (kun een sqrt)
+// TODO: Assumes all different eigenvalues. Does this break with multiples?
+// TODO: Stop after max_steps for fixed k. Return max Gershgorin radius as convergence -- or max Rayleigh quotient residual?
+// TODO: Implement implicit QR iteration using Francis' Q theorem/bulge chasing
+real_pair eigensystem_hermitian(const cdouble *A,
+				const int n,
+				real_t *lambdas, cdouble* Q)
+{
+  cdouble *T=malloc(n*n*sizeof(cdouble));	         /* T: Tridiagonal transform, dense representation */
+  real_t *V=malloc((2*(n-1))*sizeof(real_t));
+  real_t *Qhat=malloc(n*n*sizeof(real_t));	 /* V: Reflection vectors, Qhat: Inner transform */
+  memcpy(T, A, n*n*sizeof(cdouble));
+
+  real_t max_error    = gershgorin_tolerance;
+  size_t n_iterations = 0;
+
+#if TIME_EIGENSYSTEMS
+  clock_t start_time = clock();
+  if(((++nth_time) % 1000)==0) printf("%dth call to eigensystem_hermitian, cpu-time used: %gs, on average %gs per call.\n",
+				      nth_time,eigensystem_cputime, eigensystem_cputime/nth_time);  
+#endif  
+  if(Q !=0){ // Do we want to compute eigenvectors?
+    identity_matrix(Q,   n);       // Yes, so initialize Q and Qhat to identity    
+    real_identity_matrix(Qhat,n);  // After initial transform, everything becomes real
+  }
+  
+  // 1. Initial O(n^3) decomposition A = Q T Q.T to tridiagonal form
+  QHQ(T,n,Q);
+
+  /* Copy diagonal + first off-diagonal from dense T. */
+  real_t *D=malloc(n*sizeof(real_t));
+  real_t *L=malloc(n*sizeof(real_t));
+  for(int i=0;i<n;i++){
+    D[i] = creal(T[i*n+i]);
+    L[i] = (i+1<n)? creal(T[i*n+i+1]) : 0;
+  }
+
+  
+  // 2. After tridiagonal decomposition, we can do an eigenvalue
+  //    QR-iteration step in O(n), and an eigenvector QR-iteration
+  //    step in O(n^2).
+  for(int k=n-1;k>=0;k--){
+    // We start by targeting the (n,n)-eigenvalue, and gradually
+    // deflate, working on smaller and smaller submatrices.
+    real_t d  = D[k];		// d = T(k,k)
+    real_t shift = d;
+
+    // The Gershgorin disk radius is defined by just the row-sums of
+    // absolute off-diagonal elements, since T is symmetric. As T is
+    // tridiagonal, only T(k,k-1),T(k,k), and T(k,k+1) are nonzero.
+    // Thus, the k'th Gershgorin radius is just |T(k,k-1)| +
+    // |T(k,k+1)| = |T(k,k-1)| + |T(k+1,k)| = |L[k-1]|+|L[k]|.
+    int i=0;
+
+    real_t GR = (k>0?fabs(L[k-1]):0)+fabs(L[k]); // Initial Gershgorin radius
+    int not_done = 2;
+    while(not_done > 0){
+      i++;   
+      T_QTQ(k+1, D,L, D,L, V, shift);
+      if(Q!=0) apply_real_reflections(V,k,Qhat,n);
+      
+      GR = (k>0?fabs(L[k-1]):0)+(k+1<n?fabs(L[k]):0); // Gershgorin radius
+
+      // Best guess to eigenvalue in position n-1,n-1.
+      if(k>0){
+	real_pair l  = eigvalsh2x2(D[k-1],L[k-1], /* T[(k-1):k, (k-1):k] */
+				   L[k-1],D[k]  );
+	real_t l0 = l.value[0], l1 = l.value[1];
+	
+	shift    = fabs(l0-d) < fabs(l1-d)? l0 : l1; // Pick closest eigenvalue
+      } else
+	shift    = D[k];
+
+      if(GR <= gershgorin_tolerance) not_done--;      
+      if(i>max_iterations && GR > gershgorin_tolerance){
+	fprintf(stderr,"Cannot converge eigenvalue %d to tolerance " G 
+	       " in %d iterations using machine precision " G " (d=" G ", shift=" G ", GR=" G ")\n"
+	       "D[k] = " G ", L[k-1] = " G ", L[k] = " G "\n",
+	       k,gershgorin_tolerance,i,
+	       machine_precision,d,shift,GR,
+	       D[k], (k>0)?L[k-1]:0, (k+1<n)?L[k]:0);
+	
+	max_error = fmax(max_error,GR);
+	break;
+      }
+      n_iterations++;
+    }
+  }
+  for(int k=0;k<n;k++) lambdas[k] = D[k];
+
+  if(Q != 0){
+    cdouble *Q_tmp=malloc(n*n*sizeof(cdouble));
+    for(int i=0;i<n;i++){
+      cdouble *col=malloc(n*sizeof(cdouble));
+      for(int j=0;j<n;j++){
+	cdouble sum = cplx(0,0);
+	for(int k=0;k<n;k++) sum = cadd(sum,rmul(Qhat[i*n+k],Q[k*n+j]));
+	col[j] = sum;
+      }
+      for(int j=0;j<n;j++) Q_tmp[i*n+j] = col[j];
+    }
+    memcpy(Q,Q_tmp,n*n*sizeof(cdouble));
+  }
+
+#if TIME_EIGENSYSTEMS
+  clock_t end_time = clock();
+  eigensystem_cputime += (end_time-start_time)/(real_t)CLOCKS_PER_SEC;
+#endif
+  
+  real_pair result;
+  result.value[0] = max_error;
+  result.value[1] = n_iterations;
+  
+  free(T); free(V); free(Qhat); free(D); free(L);
+  return result;
+}
+
+
+#endif
diff --git a/mcstas-comps/share/eigen.h b/mcstas-comps/share/eigen.h
new file mode 100644
index 0000000000..b18a42cb0e
--- /dev/null
+++ b/mcstas-comps/share/eigen.h
@@ -0,0 +1,71 @@
+/************************************************************
+ * McStas Eigensolver library                               *
+ * Contributed by James Emil Avery, 202x                    *
+ * Department of Computer Science, University of Copenhagen *
+ ************************************************************/
+
+#ifndef EIGEN_LIB_H
+#define EIGEN_LIB_H
+
+/* Original file "types.h": */
+#include <float.h>
+
+// PW: Non-float64 logic suppressed.
+  #define real_t double
+  #define scalar cdouble
+  #define machine_precision DBL_EPSILON
+ #define PRINTF_G "%g"
+
+#ifndef  _MSC_EXTENSIONS
+#define INLINE inline __attribute__((always_inline))
+#else
+#define INLINE inline
+#endif
+
+typedef struct {
+  real_t value[2];
+} real_pair;
+
+
+/* End of "types.h"         */
+/* ------------------------ */
+/* Original file "eigen.h": */
+
+void print_vector(const char *name, cdouble *a, int l);
+void print_matrix(const char *name, cdouble *A, int m, int n);
+
+real_t vector_norm(const cdouble *x, int n);
+void   extract_region(cdouble *S, int N,
+		    int i0, int j0, int m, int n,
+		    cdouble *D);
+double max_norm(cdouble* A, int m, int n);
+void identity_matrix(cdouble *Q, int n);
+void real_identity_matrix(double *Q, int n);
+void matrix_inplace_multiply(cdouble *A, cdouble *B, 
+			     int m, int n, int q);
+
+void reflection_vector(/*in*/cdouble *a, double anorm,
+		       /*out*/cdouble *v, cdouble *sigma, int n);
+void apply_reflection(/*in/out*/cdouble *A, const cdouble *v,
+		      int m, int n, cdouble sigma, int transpose);
+
+void reflect_region(/*in/out*/cdouble *A, int N,
+		    int i0, int j0, int m, int n,
+		    const cdouble *v, cdouble sigma, int cols);
+
+void apply_real_reflections(double *V, int n, double *Q, int m);
+
+
+void QHQ(/*in/out*/cdouble *A, int n, cdouble *Q);
+void T_QTQ(const int n,
+	   const real_t *Din, const real_t *Lin,
+	   real_t *Dout, real_t *Lout, real_t *Vout,
+	   real_t shift);
+
+real_pair eigvalsh2x2(real_t a, real_t b, real_t c, real_t d);
+
+real_pair eigensystem_hermitian(const cdouble *A, int n,
+				real_t *lambdas, cdouble *Q);
+
+/* End of eigen.h */
+#endif