diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp
new file mode 100644
index 000000000..4ba9ae927
--- /dev/null
+++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp
@@ -0,0 +1,1295 @@
+/*******************************************************************************
+*
+*  McStas, neutron ray-tracing package
+*  Copyright(C) 2000 Risoe National Laboratory.
+*
+* Component: Phonon_BvK_PG 
+*
+* %I
+* Written by: Kim Lefmann
+* Date: 06.12.2025
+* Origin: NBI, KU
+*
+* A sample for phonon scattering based on cross section expressions from the Boothroyd textbook
+* based upon the algorithm from component Phonon_simple (with expressions from Squires, Ch. 3)
+* Using the Born-von Karman force constnts and the normal mode formalism
+*
+* 
+* %D
+* Single-cylinder or slap shaped shape.
+* Absorption included.
+* No multiple scattering.
+* No incoherent scattering emitted, but incoherent is present in attenuation
+* No attenuation from coherent scattering. No Bragg scattering.
+* Specialized for Pyrolytic Graphite(PG)
+*
+* Algorithm:
+* 0. Always perform the scattering if possible (otherwise ABSORB)
+* 1. Choose direction within a focusing solid angle
+* 2. Select a phonon mode at random (alternatively mode number fixed by user)
+* 3. Calculate the zeros of (E_i-E_f-hbar omega(kappa)) as a function of k_f
+* 4. Choose one value of k_f (there always at least one possible; see e.g. Squires)
+* 5. Perform the correct weight transformation
+* For details: see article A. F. Davidsen et al, manuscript for J Neutron Res 2026
+*
+* %P
+* INPUT PARAMETERS:
+* radius: [m]         Outer radius of sample in (x,z) plane
+* 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)
+* T: [K]              Temperature
+* focus_r: [m]        Radius of sphere containing target.
+* focus_xw: [m]       horiz. dimension of a rectangular area
+* focus_yh: [m]       vert.  dimension of a rectangular area
+* focus_aw: [deg]     horiz. angular dimension of a rectangular area
+* focus_ah: [deg]     vert.  angular dimension of a rectangular area
+* target_x: [m]       position of target to focus at . Transverse coordinate
+* target_y: [m]       position of target to focus at. Vertical coordinate
+* target_z: [m]       position of target to focus at. Straight ahead.
+* target_index: [1]   relative index of component to focus at, e.g. next is +1 
+* mode_input: [1]     order of the phonon mode to be scattered (unphysical). In the range 0 through 11 . If mode_input == 12, then the mode is choosen at random
+* verbose_input: [1]        toggles verbose mode (verbose>1) or eigenvector mode (verbose==1)
+* dispersion: [1]   0: do nothing special; 1: calculate and output the dispersion; 2: calculate and output the dispersion while varying v_f
+*
+* OUTPUT PARAMETERS:
+* V_rho: [AA^-3]      Unit cell density
+* V_my_s: [m^-1]      Attenuation factor due to incoherent scattering
+* V_my_a_v: [m^-1]    Attenuation factor due to absorbtion
+*
+* %L
+* The test/example instrument <a href="../examples/Test_Phonon.instr">Test_Phonon.instr</a>.
+*
+* %E
+******************************************************************************/
+
+DEFINE COMPONENT Phonon_BvK_PG
+DEFINITION PARAMETERS ()
+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)
+OUTPUT PARAMETERS (q_x, q_y, q_z)
+/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
+
+SHARE
+%{
+
+#ifndef PHONON_SIMPLE
+#define PHONON_SIMPLE $Revision$   // KL comment 040125: what is this used for?
+#pragma acc routine
+
+// ---------------- THINGS THAT COULD BE GENERAL FOR McStas --------------
+
+#define T2E (1/11.605)   // Kelvin to meV converter
+#define pi    3.14159265358979323846264338327950288419716939937510L
+#define sqrt3 1.73205080756887729352744634150587236694280525381038L
+#define THz2meV 4.13566 // THz to meV converter
+
+#define Da2kg 1.66054e-27 //Dalton to kg converter
+#define dyn2N 1e-3 // dyn/cm to N/m converter
+#define hbar 1.0546E-34
+#define mn (1.0087*Da2kg)
+
+double nbose(double omega, double T)  /* Give it another name ?? */
+  {
+    double nb;
+    nb= (omega>0) ? 1+1/(exp(omega/(T*T2E))-1) : 1/(exp(-omega/(T*T2E))-1);
+    return nb;
+  }
+#undef T2E
+
+/* Routine types from Numerical Recipies book */
+#define UNUSED (-1.11e30)
+#define MAXRIDD 60
+
+void fatalerror_cpu(char *s)
+  {
+    fprintf(stderr,"%s \n",s);
+    exit(1);
+  }
+
+  void nrerror(const char *error_text){
+  printf("Numerical Recipes run-time error ...\n");
+  printf("%s\n",error_text);
+  printf("... now exiting to system.\n");
+  exit(1);
+}
+
+#pragma acc routine
+void fatalerror(char *s)
+  {
+  #ifndef OPENACC
+    fatalerror_cpu(s);
+  #endif
+  }
+
+  #pragma acc routine
+
+void bubblesort (int size , double* inputarray , int* index)
+//this function modifies the inputarray by bubble sorting, and also modifies the index array as the
+//index after the sorting. The index array should be initialized as index[] = {0,1,2,3,4,5,6,7...}
+{
+    int tempindex = 0;
+    int i;
+    int j;
+    for (i = 0 ; i < size-1 ;  i++)
+    {
+        for (j = 0 ; j < size-1 ; j++)
+        {
+            if (inputarray[index[j]] > inputarray[index[j+1]])
+            {
+                tempindex = index[j];
+                index[j] = index[j+1];
+                index[j+1] = tempindex;
+            }
+        }
+     }
+    return;
+}
+
+
+// ---------------- END GENERAL FOR McStas ----------------------------------
+
+// ---------------- THINGS THAT ARE SPECIFIC FOR THE PG COMPONENT ------------
+
+#define SMALL_NUMBER 1E-6
+#define V_HIGH 8000   // Highest velocity to search, corresponding to 320 meV, far above the to of the PG phonon band (202 meV along main axes)
+#define MATRIX_TEST 0 // Change to 1 to write matrices to disk
+#define TIME_EIGENSYSTEMS 1
+
+#define DIM 12     // dimensionality of the phonon eigenvectors (4 atoms in the unit cell times 3 dimensions)
+#define DV 0.001   // Velocity change used for numerical derivative [m/s]  CHECK if this should be smaller
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <complex.h>
+#include "src/eigen.h"
+   
+extern FILE *matrix_test_file, *parms_test_file;
+extern int matrix_test_count;
+int mode, verbose, shape;
+
+#define a_latt 2.461 //PG lattice constant in AA; Nicklow1972 gives 2.45 (WHICH PAGE?)
+#define c_latt 6.708 // PG lattice constant in AA; Nicklow1972 gives 6.70
+#define b_length 6.646 // PG scattering legth in fm
+
+double M = 12.011 * Da2kg; //atomic mass of C
+
+#define K_l1  3.62e+5      // Force constant longitudinal nn1 - in (a,b) plane; here units of dyn
+#define K_t1  1.99e+5      // Force constant transverse nn1 - in (a,b) plane
+#define K_l2  1.33e+5      // Force constant longitudinal nn2 - in (a,b) plane
+#define K_t2 -0.520e+5     // Force constant transverse nn2 - in (a,b) plane
+#define K_l3 -0.037e+5     // Force constant longitudinal nn3 - in (a,b) plane
+#define K_t3  0.288e+5     // Force constant transverse nn3 - in (a,b) plane
+#define K_l4  0.058e+5     // Force constant longitudinal nn4 - along c
+#define K_t4  0.0077e+5    // Force constant transverse nn4 - along c
+
+    //Lattice basis
+
+    // PG is hexagonal close packed with 4 atoms per cell
+    // Coordinates from Nicklow1972, Fig. 7: Atoms A, B, C, D
+    double Delta[4*3] = {0 , 0 , 0 , a_latt/(2*sqrt3) , a_latt/2 , 0 , 0 , 0 , c_latt/2 ,  -a_latt/(2*sqrt3) , a_latt/2 , c_latt/2};
+
+    // Lattice vectors from paper fig. 7
+    double avec[3] = {a_latt*sqrt3/2 , a_latt*0.5 , 0}; // THIS IS NEVER USED ??!!
+    double bvec[3] = {a_latt*(-sqrt3/2) , a_latt*0.5 , 0};
+    double cvec[3] = {0 , 0 , c_latt};
+
+    // reciprocal lattice vectors
+    double astar = 4*PI/(sqrt3*a_latt);  // length of reciprocal lattice vector a*
+    double cstar = 2*PI/c_latt;  // length of reciprocal lattice vector c*
+
+    // Rotation matrices
+    // m(12*i+j) = M(i,j)
+    double Rot120[3][3] = {{-0.5 , sqrt3/2 , 0} , {-sqrt3/2 , -0.5 , 0} , {0 , 0 , 1}};
+    double Rot120_flat[9] = {-0.5 , sqrt3/2 , 0 , -sqrt3/2 , -0.5 , 0 , 0 , 0 , 1};
+    double Rot60[3][3] = {{0.5 , sqrt3/2 , 0} , {-sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}};
+    double Rot120rev[3][3] = {{-0.5 , -sqrt3/2 , 0} , {sqrt3/2 , -0.5 , 0} , {0 , 0 , 1}}; //rotate -120
+    double Rot60rev[3][3] = {{0.5 , -sqrt3/2 , 0} , {sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate -60
+    double Rot180[3][3] = {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}}; // rotate 180 or -180
+//    double Rot5[3][3] = {{0.5 , sqrt3/2 , 0} , {-sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate 60
+//    double Rot5rev[3][3] = {{0.5 , -sqrt3/2 , 0} , {sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate -60
+
+    //1st neighbour
+
+//    double r_j1[3][3] = {{-a_latt/sqrt3 , 0 , 0} , {a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0}}; // This holds from atoms A and D
+//    double r_j2[3][3] = {{a_latt/sqrt3 , 0 , 0} , {-a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0}}; // This holds from atoms B and C
+
+    double Phi_nn1[3][3] = {{K_l1*dyn2N , 0 , 0} , {0 , K_t1*dyn2N , 0} , {0 , 0 , K_t1*dyn2N}};
+    double Phi1[3][3] = {{K_l1*dyn2N , 0 , 0} , {0 , K_t1*dyn2N , 0} , {0 , 0 , K_t1*dyn2N}}; //Phi1 = Phi_nn1
+    double Phi2[3][3];
+    double Phi3[3][3];
+
+    //2nd neighbour
+//    double r_j11[9][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0}};//r_j1 + r_add
+//    double r_j21[9][3] = {{-a_latt*(-1/sqrt3) , -a_latt*0 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*(-0.5) , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0}};//r_j2 + r_add
+
+    double Phi_nn2[3][3] = {{K_t2*dyn2N , 0 , 0} , {0 , K_l2*dyn2N , 0} , {0 , 0 , K_t2*dyn2N}};
+    double Phi4[3][3] = {{K_t2*dyn2N , 0 , 0} , {0 , K_l2*dyn2N , 0} , {0 , 0 , K_t2*dyn2N}};
+    double Phi5[3][3];
+    double Phi6[3][3];
+    double Phi7[3][3];
+    double Phi8[3][3];
+    double Phi9[3][3];
+
+    //3rd neighbour
+//    double r_j12[12][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {2*a_latt/sqrt3 , 0 , 0} , {-a_latt/sqrt3 , a_latt , 0} , {-a_latt/sqrt3 , -a_latt , 0}};//r_j11 + r_add2
+    double r_j2[12][3] = {{a_latt/sqrt3 , 0 , 0} , {-a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {-2*a_latt/sqrt3 , 0 , 0} , {a_latt/sqrt3 , -a_latt , 0} , {a_latt/sqrt3 , a_latt , 0}};//r_j21 - r_add2
+
+    double Phi_nn3[3][3] = {{K_l3*dyn2N , 0 , 0} , {0 , K_t3*dyn2N , 0} , {0 , 0 , K_t3*dyn2N}};
+
+    // c2 = 2/sqrt(6);
+    // c1 = 1/sqrt(6);
+    // c0 = 1/sqrt(2);
+    // c3 = 1/sqrt(3);
+
+    double Phi10[3][3] = {{K_l3*dyn2N , 0 , 0} , {0 , K_t3*dyn2N , 0} , {0 , 0 , K_t3*dyn2N}}; //Phi_nn3
+    double Phi11[3][3];
+    double Phi12[3][3];
+
+    //4th neighbour
+    double r_j13[14][3] = {{-a_latt/sqrt3 , 0 , 0} , {a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {2*a_latt/sqrt3 , 0 , 0} , {-a_latt/sqrt3 , a_latt , 0} , {-a_latt/sqrt3 , -a_latt , 0} , {0 , 0 , c_latt/2} , {0 , 0 , -c_latt/2}};//r_j12 + r_add3 Sublattice B and D do not have this coupling
+
+    double Phi_nn4[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}};
+
+    double Phi13[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; //Phi_nn4;
+    double Phi14[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; //Phi_nn4;
+
+// ----------------- END SPECIFIC FOR THE PG COPONENT --------------------------
+
+
+
+/* TODO: Det er meget langsomt at bruge pointere. Erstat med flad memory. */
+// m(12*i+j) = M(i,j)
+//3*3 matrices multiplication
+void matrixmultiplication( double matrix1[3][3], double matrix2[3][3], double result[3][3])
+    {
+        // multiplies matrix1 by matrix2 and places the result in result
+//      double matrix1[3][3]={{*(array1), *(array1+1), *(array1+2)}, {*(array1+3), *(array1+4), *(array1+5)}, {*(array1+6), *(array1+7), *(array1+8)}};
+//      double matrix2[3][3]={{*(array2), *(array2+1), *(array2+2)}, {*(array2+3), *(array2+4), *(array2+5)}, {*(array2+6), *(array2+7), *(array2+8)}};
+
+      result[0][0]=matrix2[0][0]*matrix1[0][0]+matrix2[1][0]*matrix1[0][1]+matrix2[2][0]*matrix1[0][2];
+      result[0][1]=matrix2[0][1]*matrix1[0][0]+matrix2[1][1]*matrix1[0][1]+matrix2[2][1]*matrix1[0][2];
+      result[0][2]=matrix2[0][2]*matrix1[0][0]+matrix2[1][2]*matrix1[0][1]+matrix2[2][2]*matrix1[0][2];
+      result[1][0]=matrix2[0][0]*matrix1[1][0]+matrix2[1][0]*matrix1[1][1]+matrix2[2][0]*matrix1[1][2];
+      result[1][1]=matrix2[0][1]*matrix1[1][0]+matrix2[1][1]*matrix1[1][1]+matrix2[2][1]*matrix1[1][2];
+      result[1][2]=matrix2[0][2]*matrix1[1][0]+matrix2[1][2]*matrix1[1][1]+matrix2[2][2]*matrix1[1][2];
+      result[2][0]=matrix2[0][0]*matrix1[2][0]+matrix2[1][0]*matrix1[2][1]+matrix2[2][0]*matrix1[2][2];
+      result[2][1]=matrix2[0][1]*matrix1[2][0]+matrix2[1][1]*matrix1[2][1]+matrix2[2][1]*matrix1[2][2];
+      result[2][2]=matrix2[0][2]*matrix1[2][0]+matrix2[1][2]*matrix1[2][1]+matrix2[2][2]*matrix1[2][2];
+
+      return;
+    }
+
+// TODO: (James Avery) Benchmark against pointer versions: matrix3mupltiplication or matrixmultiplication
+// m(12*i+j) = M(i,j)
+void matrix3x3_multiply3(const double A[9], const double B[9], const double C[9], double result[9])
+{
+  // R_{i,j} =  \sum_{k,l=1}^3 A_{i,k}*B_{k,l}*C_{l,j}
+  for(int i=0;i<3;i++)
+    for(int j=0;j<3;j++){
+      double ij_sum = 0;
+      for(int k=0;k<3;k++)
+	for(int l=0;l<3;l++)
+	  ij_sum += A[i*3+k]*B[k*3+l]*C[l*3+j];
+      result[i*3+j] = ij_sum;
+    }
+}
+
+/* TODO: (James Avery) It is slow to use points. Replace with flat memory. */
+  // m(12*i+j) = M(i,j)
+  // Counter argument: (Kim Lefmann) pointers are now only used in the initialization. In the TRACE code, memory is flat
+  void matrix3multiplication(double matrix1[3][3], double matrix2[3][3], double matrix3[3][3], double result[3][3])
+    {
+        // Performs the multiplication (matrix1*matrix2)*matrix3 and placed the result in result
+//#define TEST_MULTIPLY
+        double temp[3][3];
+        matrixmultiplication(matrix1,matrix2,temp);
+        matrixmultiplication(temp,matrix3,result);
+#ifdef TEST_MULTIPLY
+        printf("Matrix3multiply called with \n");
+        printf("matrix1 = (");
+        for (int i=0; i<3; i++)
+        {
+            printf("( %g %g %g), ",matrix1[i][0], matrix1[i][1], matrix1[i][2]);
+        }
+        printf(")\n Matrix2 = (");
+        for (int i=0; i<3; i++)
+        {
+            printf("( %g %g %g), ",matrix2[i][0], matrix2[i][1], matrix2[i][2]);
+        }
+        printf(")\n Matrix3 = (");
+                for (int i=0; i<3; i++)
+        {
+            printf("( %g %g %g), ",matrix3[i][0], matrix3[i][1], matrix3[i][2]);
+        }
+        printf(")\n Temp = (");
+                for (int i=0; i<3; i++)
+        {
+            printf("( %g %g %g), ",temp[i][0], temp[i][1], temp[i][2]);
+        }
+        printf(")\n Result = (");
+                 for (int i=0; i<3; i++)
+        {
+            printf("( %g %g %g), ",result[i][0], result[i][1], result[i][2]);
+        }
+        printf(")\n");
+#endif // TEST_MULTIPLY
+
+        return;
+    }
+    
+
+/* TODO: (James Avery) Replace by standard complex.h operation */
+complex complexexp_qdotr (double* q, double* rj)
+{
+//    double qrjtest
+    double qrj = *(q)**(rj)+*(q+1)**(rj+1)+*(q+2)**(rj+2);
+//    printf(" q= (%g %g %g), r = (%g %g %g) qrj = %g \n",q[0],q[1],q[2],rj[0],rj[1],rj[2],qrj);
+//     return (cexp(qrj));    // KL comment 030125 WHY does this not run?
+    return (cos(qrj)+I*sin(qrj));
+}
+
+
+void initialize_omega_q()
+{
+    // Make the matrix algebra that finalizes the force constant set-up
+
+    matrix3multiplication(Rot120,Phi_nn1,Rot120rev,Phi2); //Phi2=Rot120*Phi_nn1*Rot120^(-1)
+    matrix3multiplication(Rot120rev,Phi_nn1,Rot120,Phi3); //Phi3=Rot120^{-1}*Phi_nn1*Rot120
+    matrix3multiplication(Rot60,Phi_nn2,Rot60rev,Phi5);   //Phi5=Rot60*Phi_nn2*Rot60^{-1}
+    matrix3multiplication(Rot120,Phi_nn2,Rot120rev,Phi6); //Phi6=Rot120*Phi_nn2*Rot120^(-1);
+    matrix3multiplication(Rot180,Phi_nn2,Rot180,Phi7);    //Phi7=Rot180*Phi_nn2*Rot180;
+    matrix3multiplication(Rot120rev,Phi_nn2,Rot120,Phi8); //Phi8=Rot120^{-1}*Phi_nn2*Rot120;
+    matrix3multiplication(Rot60rev,Phi_nn2,Rot60,Phi9);   //Phi9=Rot60^{-1}*Phi_nn2*Rot60;
+    matrix3multiplication(Rot120,Phi_nn3,Rot120rev,Phi11);//Phi11=Rot120*Phi_nn3*Rot120^(-1);
+    matrix3multiplication(Rot120rev,Phi_nn3,Rot120,Phi12);//Phi12=Rot120^{-1}*Phi_nn3*Rot120;
+
+    return;
+}
+
+//---------------------------- START CENTRAL FUNCTION OMEGA-Q ----------------------------
+
+    double omega_q(complex* parms, complex* eigenvectorgood) {
+    
+    // _______________ DETERMINE Q ____________________________
+    
+      double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z;
+      double q, qx, qy, qz, Jq, hw_phonon, hw_neutron;
+      double h,k,l,hk_inplane, hk_outofplane;
+      int disp;
+
+      for(int ii=0;ii<8;ii++) if(isnan(creal(parms[ii])) || isnan(cimag(parms[ii]))){
+	  fprintf(stderr,"omega_q: Array parms contains a nan at position %d.\n",ii);
+	  fprintf(stderr,"parms = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g + i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n");
+	  abort();
+	}
+
+      vf=parms[0];
+      vi=parms[1];
+      vv_x=parms[2];
+      vv_y=parms[3];
+      vv_z=parms[4];
+      vi_x=parms[5];
+      vi_y=parms[6];
+      vi_z=parms[7];
+      h = parms[9];
+      k = parms[10];
+      l = parms[11];
+      disp = (int)parms[12];
+ 
+      qx=V2K*(vi_x-vf*vv_x);
+      qy=V2K*(vi_y-vf*vv_y);
+      qz=V2K*(vi_z-vf*vv_z);
+    
+      if (disp==1) /* calculate dispersion directly */
+      {
+          printf("Dispersion calculation: ");
+          qx = h*astar;
+          qy = k*astar; //correct only for qy=0 !
+          qz = l*cstar;
+      }
+      
+    double qvec[3]={qx,qy,qz};
+    
+    if (verbose>=4)
+        printf("Enter omega_q; q = (%g, %g, %g) \n",qx,qy,qz);
+    
+    // ________________ END DETERMINE Q ________________-
+    
+
+    complex Phi_diag[3*3];
+    int i;
+    int j;
+    for (i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            Phi_diag[i*3+j] = Phi1[i][j]+Phi2[i][j]+Phi3[i][j]+Phi4[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0]))+Phi5[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0]))+Phi6[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0]))+Phi7[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0]))+Phi8[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0]))+Phi9[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0]))+Phi10[i][j]+Phi11[i][j]+Phi12[i][j];
+        }
+    }
+    
+    complex Phi_diag1[3*3];
+    for (i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            Phi_diag1[i*3+j] = Phi13[i][j]+Phi14[i][j];
+        }
+    }
+    
+    complex Phi_offdiag1[3*3];
+    complex Phi_offdiag2[3*3];
+    for (i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            Phi_offdiag1[i*3+j] =-Phi1[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[0][0]))-Phi2[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[1][0]))-Phi3[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[2][0]))-Phi10[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[9][0]))-Phi11[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[10][0]))-Phi12[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[11][0]));
+            Phi_offdiag2[i*3+j] = conj(Phi_offdiag1[i*3+j]);
+        }
+    }
+    
+    complex Phi_6Dim[6*6];
+    for (i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            Phi_6Dim[i*6+j] = Phi_diag[i*3+j]+Phi_diag1[i*3+j];
+            Phi_6Dim[i*6+j+3] = Phi_offdiag1[i*3+j];
+            Phi_6Dim[(i+3)*6+j] = Phi_offdiag2[i*3+j];
+            Phi_6Dim[(i+3)*6+j+3] = Phi_diag[i*3+j];  // Because Phi_diag1 does not appear for the B,C, sublattices
+        }
+    }
+    
+    complex Phi_AC[3*3];
+    for (i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            Phi_AC[i*3+j] = -Phi13[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[12][0]))-Phi14[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[13][0]));
+        }
+    }
+    
+//    complex Zero_3Dim[3][3] = {{0 , 0 , 0},{0 , 0 , 0},{0 , 0 , 0}};
+
+    complex Phi_6Dim_offdiag[6*6];
+    for (i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            Phi_6Dim_offdiag[i*6+j] = Phi_AC[i*3+j];
+            Phi_6Dim_offdiag[i*6+j+3] = (complex)0.0;
+            Phi_6Dim_offdiag[(i+3)*6+j] = (complex)0.0;
+            Phi_6Dim_offdiag[(i+3)*6+j+3] = (complex)0.0;
+        }
+    }
+
+    complex Phi_12Dim[12*12];
+    for (i = 0; i < 6; i++) {
+        for (j = 0; j < 6; j++) {
+            Phi_12Dim[i*DIM+j] = Phi_6Dim[i*6+j];
+            Phi_12Dim[i*DIM+j+6] = Phi_6Dim_offdiag[i*6+j];
+            Phi_12Dim[(i+6)*DIM+j] = conj(Phi_6Dim_offdiag[i*6+j]);
+            Phi_12Dim[(i+6)*DIM+j+6] = Phi_6Dim[i*6+j];
+        }
+    }
+    
+#ifdef TEST_MATRICES
+    printf("Phi_12Dim=\n");
+    for (i = 0; i < 12; i++) {
+      for (j = 0; j< 12; j++) {
+	if (j == 0) {
+	  printf("{");
+	}
+	printf("%g+(%g)i  ", creal(Phi_12Dim[i*DIM+j]) , cimag((Phi_12Dim[i*DIM+j])));
+	if (j == 11) {
+	  printf("}\n");
+	}
+      }
+    }
+    
+    printf("Phi_6Dim=\n");
+    for (i = 0; i < 6; i++) {
+      for (j = 0; j< 6; j++) {
+	if (j == 0) {
+	  printf("{");
+	}
+	printf("%g+(%g)i  ", creal(Phi_6Dim[i*6+j]) , cimag((Phi_6Dim[i*6+j])));
+	if (j == 5) {
+	  printf("}\n");
+	}
+      }
+    }
+
+    printf("Phi_diag=\n");
+    for (i = 0; i < 3; i++) {
+      for (j = 0; j< 3; j++) {
+	if (j == 0) {
+	  printf("{");
+	}
+	printf("%g+(%g)i  ", creal(Phi_diag[i*3+j]) , cimag((Phi_diag[i*3+j])));
+	if (j == 2) {
+	  printf("}\n");
+	}
+      }
+    }
+    printf("Phi2={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}}\n",Phi2[0][0],Phi2[0][1],Phi2[0][2],Phi2[1][0],Phi2[1][1],Phi2[1][2],Phi2[2][0],Phi2[2][1],Phi2[2][2]);
+    printf("Phi3={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi3[0][0],Phi3[0][1],Phi3[0][2],Phi3[1][0],Phi3[1][1],Phi3[1][2],Phi3[2][0],Phi3[2][1],Phi3[2][2]);
+    printf("Phi5={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi5[0][0],Phi5[0][1],Phi5[0][2],Phi5[1][0],Phi5[1][1],Phi5[1][2],Phi5[2][0],Phi5[2][1],Phi5[2][2]);
+    printf("Phi6={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi6[0][0],Phi6[0][1],Phi6[0][2],Phi6[1][0],Phi6[1][1],Phi6[1][2],Phi6[2][0],Phi6[2][1],Phi6[2][2]);
+    printf("Phi7={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi7[0][0],Phi7[0][1],Phi7[0][2],Phi7[1][0],Phi7[1][1],Phi7[1][2],Phi7[2][0],Phi7[2][1],Phi7[2][2]);
+    printf("Phi8={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi8[0][0],Phi8[0][1],Phi8[0][2],Phi8[1][0],Phi8[1][1],Phi8[1][2],Phi8[2][0],Phi8[2][1],Phi8[2][2]);
+    printf("Phi9={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi9[0][0],Phi9[0][1],Phi9[0][2],Phi9[1][0],Phi9[1][1],Phi9[1][2],Phi9[2][0],Phi9[2][1],Phi9[2][2]);
+    printf("Phi11={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi11[0][0],Phi11[0][1],Phi11[0][2],Phi11[1][0],Phi11[1][1],Phi11[1][2],Phi11[2][0],Phi11[2][1],Phi11[2][2]);
+    printf("Phi12={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi12[0][0],Phi12[0][1],Phi12[0][2],Phi12[1][0],Phi12[1][1],Phi12[1][2],Phi12[2][0],Phi12[2][1],Phi12[2][2]);
+    printf("Phi14={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi14[0][0],Phi14[0][1],Phi14[0][2],Phi14[1][0],Phi14[1][1],Phi14[1][2],Phi14[2][0],Phi14[2][1],Phi14[2][2]);
+
+    printf(" 1-exp(q.r_j13[3]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0])));
+    printf(" 1-exp(q.r_j13[4]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0])));
+    printf(" 1-exp(q.r_j13[5]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0])));
+    printf(" 1-exp(q.r_j13[6]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0])));
+    printf(" 1-exp(q.r_j13[7]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0])));
+    printf(" 1-exp(q.r_j13[8]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0])));
+    
+#endif // TEST_MATRICES
+    
+  double eigenvalue[DIM]; 
+  complex Q[DIM*DIM];
+  complex matrix[DIM*DIM];
+  int index[DIM];
+  
+  
+  for (i=0; i<DIM; i++)   // transforms from 2D matrix to 1D pseudo-matrix
+    for (j=i; j<DIM; j++){
+      matrix[i*DIM+j] = Phi_12Dim[i*DIM+j];
+      matrix[j*DIM+i] = conj(matrix[i*DIM+j]); // Forces the pseudo-matrix Hermitean
+    }  // m(12*i+j) = M(i,j)
+
+#if MATRIX_TEST
+    fwrite(matrix,sizeof(complex),DIM*DIM,matrix_test_file);
+    fwrite(parms, sizeof(complex),8,    parms_test_file);
+#endif  
+  
+ eigensystem_hermitian(matrix, DIM, eigenvalue, Q); /* Produces unitary transformation matrix Q: eigenvectors are columns */
+
+   for (i=0; i<DIM; i++)  index[i]=i;
+
+  if (verbose >= 6)
+  {
+      printf("Before bubblesort\n");
+      for (i=0; i<DIM; i++)
+      {
+          printf("i = %i eigenvalue[i] = %g index[i] %i \n",i,eigenvalue[i],index[i]);
+      }
+  }
+
+  // Sort eigenvalues from lowest to highest
+  bubblesort(DIM, eigenvalue, index);  //mergesort(eigenvalue,DIM,index);
+  if (verbose >= 5)
+  {
+      printf("After bubblesort\n");
+      for (i=0; i<DIM; i++)
+      {
+          printf("i = %i eigenvalue[i] = %g index[i] %i \n",i,eigenvalue[i],index[i]);
+      }
+  }   
+
+  double eigenvaluegood[DIM];
+  
+    for (j = 0; j < DIM ; j++){
+      eigenvectorgood[j] = Q[index[mode]*DIM + j]; /* Eigenvectors are columns of Q */ // No, they are rows!
+    }
+
+      if (verbose>=7)
+      {
+          printf("Q = [");
+          for (i = 0; i < DIM ; i++)
+          {
+              printf("\n (");
+              for (j=0; j<DIM; j++)
+                  printf(" %g +i %g, ",creal(Q[i*DIM+j]),cimag(Q[i*DIM+j]));
+              printf(")");
+          }
+          printf("]\n");
+      }
+
+      for (i=0; i<DIM; i++)
+        eigenvaluegood[i] = sqrt(creal(eigenvalue[index[i]]))/sqrt(M)/(2*PI*1E12)*THz2meV; // convert to units of THz
+	
+    if (verbose>=5)
+    {
+        printf("Diagonalization done, eigenvalues (energies) are: \n (");
+        for (i=0; i<DIM; i++)
+            printf(" %g, ",eigenvaluegood[i]);
+        printf(" )\n");
+        printf("mode = %i , eigenvectorgood[mode] = (",mode);
+        for (j=0; j<DIM; j++)
+            printf(" %g +i %g, ",creal(eigenvectorgood[j]),cimag(eigenvectorgood[j]));
+        printf(" )\n");
+    }
+        
+  // ___________________ RETURN RESULTS TO CALLING FUNCTION __________________
+  
+  hw_neutron = fabs(VS2E*(vi*vi-vf*vf)); // neutron energy transfer 
+  // fabs used to find both positive and negative solutions, controlled by findroots()  
+  hw_phonon = eigenvaluegood[mode];  // phonon energy
+
+  if(verbose >= 4){
+    printf("Return from omega_q \n");
+    printf("hw_neutron = %g (vi = %g, vf = %g)\n",hw_neutron, vi, vf);
+    printf("hw_phonon  = %g = eigenvaluegood[%d] = %g = eigenvalues[%d] = %g\n",
+	    hw_phonon,mode,eigenvaluegood[mode],index[mode], eigenvalue[index[mode]]);
+//    printf("old parms  = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g+i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n");
+  }
+  
+  parms[8] = hw_phonon;
+  
+  if (disp==1) {
+      printf("mode = %d, (h,k,l) = (%g,%g,%g), hw_q = %g",mode,h,k,l,hw_phonon);
+      printf(" polarization: ( ");
+      for (j = 0; j < DIM; j++)
+      {
+        printf("%g + i %g ,",creal(eigenvectorgood[j]),cimag(eigenvectorgood[j]));
+      }
+      printf(")\n");
+  }
+  if (disp==2) {
+      hk_inplane = qx/astar;
+      hk_outofplane = qy/astar;   // Derivation assume astar along x and 60 degrees between astar and bstar:
+                                  // h_ip = h + k/2  h_op = sqrt(3)*k/2  =>  k = 2/sqrt(3)*h_op   h = h_ip - h_op/sqrt(3)
+      k = 2/sqrt(3)*hk_outofplane;
+      h = hk_inplane - hk_outofplane/sqrt(3);
+      l = qz/cstar;
+      printf("mode = %d, vf = %g, (h,k,l) = (%g,%g,%g), hw_q = %g",mode,vf,h,k,l,hw_phonon);
+      printf(" polarization: (%g + i %g , %g + i %g , %g + i %g) \n",creal(eigenvectorgood[0]),cimag(eigenvectorgood[0]),creal(eigenvectorgood[1]),cimag(eigenvectorgood[1]), creal(eigenvectorgood[2]),cimag(eigenvectorgood[2]));
+  }
+
+  return (hw_phonon - hw_neutron);
+}
+
+//--------------------- END CENTRAL FUNCTION OMEGA-Q ------------------------
+
+// ------------------START OF THE GENERAL ROOT-FINDING CODE--------------------
+
+double last_fh, last_x2;
+double zridd(double (*func)(complex*,complex*), double x1, double x2, complex *parms, complex* vector, double xacc)
+    {
+      int j;
+      double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew;
+
+      //      fprintf(stderr,"i zridd(%g,%g, parms,%g)\n",x1,x2,xacc);
+      parms[0]=x1;
+      //      fprintf(stderr,"zridd fl = omega_q\n");      
+      if (xacc>0)
+      {
+          fl=(*func)(parms,vector);
+      }
+      else
+      {
+          fl = last_fh;
+          xacc = -xacc;
+          if (fabs(x1-last_x2)>xacc)   // KL comment 030125 WHY does the library routine zridd use complex abs() function ?? changed to fabs()
+              printf("*** error in zridd *** x1: %g last_x2: %g \n",x1,last_x2);
+      }
+      //      fprintf(stderr,"fl = %g\n",fl);                  
+      parms[0]=x2;
+      //      fprintf(stderr,"zridd fh = omega_q; parms[0] = x2 = %g+i%g = %g\n",creal(parms[0]), cimag(parms[0]),x2);
+      last_fh=fh=(*func)(parms,vector);
+      last_x2 = x2;
+      //      fprintf(stderr,"fh = %g\n",fh);            
+      if (fl*fh >= 0)
+      {
+        if (fl==0) return x1;
+        if (fh==0) return x2;
+        return UNUSED;
+      }
+      else
+      {
+        xl=x1;
+        xh=x2;
+//        printf("zridd sign change: v_low : %g v_high: %g \n",xl,xh);
+        ans=UNUSED;
+        for (j=1; j<MAXRIDD; j++)
+        {
+          xm=0.5*(xl+xh);
+          parms[0]=xm;
+	  //	  fprintf(stderr,"zridd fm = omega_q\n");	  
+          fm=(*func)(parms,vector);
+          s=sqrt(fm*fm-fl*fh);
+          if (s == 0.0)
+            return ans;
+          xnew=xm+(xm-xl)*((fl >= fh ? 1.0 : -1.0)*fm/s);
+          if (fabs(xnew-ans) <= xacc)
+            return ans;
+          ans=xnew;
+          parms[0]=ans;
+	  //	  fprintf(stderr,"s = %g;  fl, fm, fh = %g,%g,%g; fnew = %g; fm*fm-fl*fh = %g\n",
+	  //		  s, fl, fm,fh, fnew, fm*fm-fl*fh);
+	  //fprintf(stderr,"zridd fnew = omega_q\n");
+          fnew=(*func)(parms,vector);
+          if (fnew == 0.0) return ans;
+          if (fabs(fm)*SIGN(fnew) != fm)
+          {
+            xl=xm;
+            fl=fm;
+            xh=ans;
+            fh=fnew;
+          }
+          else
+            if (fabs(fl)*SIGN(fnew) != fl)
+            {
+              xh=ans;
+              fh=fnew;
+            }
+            else
+              if(fabs(fh)*SIGN(fnew) != fh)
+              {
+                xl=ans;
+                fl=fnew;
+              }
+              else
+                fatalerror("never get here in zridd");
+          if (fabs(xh-xl) <= xacc)
+            return ans;
+        }
+        fatalerror("zridd exceeded maximum iterations");
+      }
+      return 0.0;  /* Never get here */
+    }
+
+#pragma acc routine 
+double zridd_gpu(double x1, double x2, complex *parms, complex* vector, double xacc)
+    {
+      int j;
+      double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew;
+
+      parms[0]=x1;
+      //      fprintf(stderr,"fl=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0]));                  
+      fl=omega_q(parms,vector);
+      parms[0]=x2;
+      //      fprintf(stderr,"fh=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0]));                  
+      fh=omega_q(parms,vector);
+      if (fl*fh >= 0)
+      {
+        if (fl==0) return x1;
+        if (fh==0) return x2;
+        return UNUSED;
+      }
+      else
+      {
+        xl=x1;
+        xh=x2;
+        ans=UNUSED;
+        for (j=1; j<MAXRIDD; j++)
+        {
+          xm=0.5*(xl+xh);
+          parms[0]=xm;
+	  //	  fprintf(stderr,"fm=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0]));            	  
+          fm=omega_q(parms,vector);
+          s=sqrt(fm*fm-fl*fh);
+          if (s == 0.0)
+            return ans;
+          xnew=xm+(xm-xl)*((fl >= fh ? 1.0 : -1.0)*fm/s);
+          if (fabs(xnew-ans) <= xacc)
+            return ans;
+          ans=xnew;
+          parms[0]=ans;
+	  //	  fprintf(stderr,"fnew=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0]));            	  
+          fnew=omega_q(parms,vector);
+          if (fnew == 0.0) return ans;
+          if (fabs(fm)*SIGN(fnew) != fm)
+          {
+            xl=xm;
+            fl=fm;
+            xh=ans;
+            fh=fnew;
+          }
+          else
+            if (fabs(fl)*SIGN(fnew) != fl)
+            {
+              xh=ans;
+              fh=fnew;
+            }
+            else
+              if(fabs(fh)*SIGN(fnew) != fh)
+              {
+                xl=ans;
+                fl=fnew;
+              }
+              else
+                fatalerror("never get here in zridd");
+          if (fabs(xh-xl) <= xacc)
+            return ans;
+        }
+        fatalerror("zridd exceeded maximum iterations");
+      }
+      return 0.0;  /* Never get here */
+    }
+ 
+#define ROOTACC 1e-8
+
+  int findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(complex*,complex*), complex *parms, complex* vector, int low_steps, int high_steps)
+  // *f  is here omega_q
+  // *parms is here p_call
+    {
+      double root,range;
+      int i;
+
+     if (verbose==2)
+          printf("Enter findroots() \n");
+     
+ // First, find roots in energy loss side, if any
+     range = brack_mid-brack_low;
+     for (i=0; i<(low_steps-1); i++) 
+     {
+        if (i==0)
+           root = zridd(f, brack_low+range*i/low_steps,
+                   brack_low+range*(i+1)/low_steps,
+                   (complex *)parms, (complex *)vector, ROOTACC);
+        else
+           root = zridd(f, brack_low+range*i/low_steps,
+                   brack_low+range*(i+1)/low_steps,
+                   (complex *)parms, (complex *)vector, -ROOTACC);  // Re-use the previous fh value
+      if (root != UNUSED)
+      {
+        list[(*index)++]=root;
+        if (verbose >=4)
+            printf("--- findroots returned a root on the energy loss side: vf = %g \n",root);
+      }
+     }
+   
+     range = (brack_mid-brack_low)/(double)low_steps;
+     for (i=0; i<low_steps; i++)  // Small steps close to zero energy transfer
+     {
+         if (low_steps==1)  // then the previous loop did not execute
+            root = zridd(f, brack_low+range*(low_steps-1+i/(double)low_steps),
+                   brack_low+range*(low_steps-1+(i+1)/(double)low_steps),
+                   (complex *)parms, (complex *)vector, ROOTACC);
+        else
+            root = zridd(f, brack_low+range*(low_steps-1+i/(double)low_steps),
+                   brack_low+range*(low_steps-1+(i+1)/(double)low_steps),
+                   (complex *)parms, (complex *)vector, -ROOTACC);  // Re-use the previous fh value
+      if (root != UNUSED)
+      {
+        list[(*index)++]=root;
+        if (verbose >=4)
+            printf("--- findroots returned a root on the energy loss side: vf = %g \n",root);
+      }
+     }
+
+ // Second, find roots in energy gain side, there is always some
+     range = (brack_high-brack_mid)/(double)high_steps;
+     for (i=0; i<high_steps; i++) // Close to zero: small steps
+     {
+         if (i==0)
+            root = zridd(f, brack_mid+range*i/(double)high_steps, 
+                   brack_mid+range*(i+1)/(double)high_steps, 
+                   (complex *)parms, (complex *)vector, ROOTACC);
+        else
+            root = zridd(f, brack_mid+range*i/(double)high_steps, 
+                   brack_mid+range*(i+1)/(double)high_steps, 
+                   (complex *)parms, (complex *)vector, -ROOTACC);  // Re-use the previous fh value
+      if (root != UNUSED)
+      {
+        list[(*index)++]=root;
+        if (verbose >= 4)
+            printf("*** findroots returned a root on the energy gain side: vf = %g \n",root);
+      }
+     }
+     
+     range = brack_high-brack_mid;
+     for (i=1; i<high_steps; i++)  // Larger steps away from zero
+     {
+      root = zridd(f, brack_mid+range*i/(double)high_steps, 
+                   brack_mid+range*(i+1)/(double)high_steps, 
+                   (complex *)parms, (complex *)vector, -ROOTACC);  // Re-use the previous fh value - there was always one
+      if (root != UNUSED)
+      {
+        list[(*index)++]=root;
+        if (verbose >= 4)
+            printf("*** findroots returned a root on the energy gain side: vf = %g \n",root);
+      }
+     }
+     
+    }
+  
+#pragma acc routine 
+  int findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, double *parms, complex* vector, int low_steps, int high_steps)
+    {
+      double root,range=brack_mid-brack_low;
+      int i;
+
+     for (i=0; i<low_steps; i++)
+     {
+       root = zridd_gpu(brack_low+range*i/(int)low_steps,
+                   brack_low+range*(i+1)/(int)low_steps,
+                   (complex *)parms, (complex *)vector, ROOTACC);
+      if (root != UNUSED)
+      {
+        list[(*index)++]=root;
+      }
+     } 
+    for (i=0; i<high_steps; i++)
+     {
+      root = zridd_gpu(brack_mid+range*i/(int)high_steps,
+                   brack_high+range*(i+1)/(int)high_steps,
+                   (complex *)parms, (complex *)vector, ROOTACC);
+       if (root != UNUSED)
+      {
+        list[(*index)++]=root;
+      }
+     }
+    }
+  
+#undef UNUSED
+#undef MAXRIDD
+#endif
+%}
+
+DECLARE
+%{
+  double V_rho;
+  double V_my_s;
+  double V_my_a_v;
+  complex **Matrix;
+  int i,j;
+  double q[3];    /* Scattering vector */
+  double qx;
+  double qy;
+  double qz;
+  double q_x;
+  double q_y;
+  double q_z;
+  complex eigenvectormode[DIM];
+  complex p_call[15];  /* Parameter list to transfer to omega_q; elsewhere known as parms */
+%}
+
+INITIALIZE
+%{
+  int i,j;
+
+  if (xwidth && yheight && zdepth)  shape=1; /* box */
+  else if (radius > 0 &&  yheight)  shape=0; /* cylinder */
+
+
+  V_rho = 1/(a_latt*a_latt*c_latt*sqrt3/2);  // (unit: Å^-3)
+  V_my_s = (V_rho * 100 * sigma_inc);
+  V_my_a_v = (V_rho * 100 * sigma_abs * 2200);
+
+  /* now compute target coords if a component index is supplied */
+  if (!target_index && !target_x && !target_y && !target_z) target_index=1;
+  if (target_index){
+    Coords ToTarget;
+    ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index), POS_A_CURRENT_COMP);
+    ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
+    coords_get(ToTarget, &target_x, &target_y, &target_z);
+  }
+  if (!(target_x || target_y || target_z)) {
+    printf("Phonon_BkV_PG: %s: The target is not defined. Using direct beam (Z-axis).\n",
+      NAME_CURRENT_COMP);
+    target_z=1;
+  }
+  initialize_omega_q();
+  verbose = verbose_input;
+
+    // m(12*i+j) = M(i,j)
+    // OK for: Rot120
+  for (i=0; i<3; i++)
+    for (j=0; j<3; j++)
+    {
+      printf("i,j: %i, %i  ROT120: %g  ROT120FLAT; %g \n",i,j,Rot120[i][j],Rot120_flat[3*i+j]);
+    }
+  %}
+
+TRACE
+%{
+  double t0, t1, t2, t3;    /* Entry/exit times for shape */
+  double v_i, v_f;          /* Neutron velocities: initial, final */
+  double vx_i, vy_i, vz_i;  /* Neutron initial velocity vector */
+  double dt0, dt;           /* Flight times through sample */
+  double l_full;            /* Flight path length for non-scattered neutron */
+  double l_i, l_o;          /* Flight path lenght in/out for scattered neutron */
+  double my_a_i;            /* Initial attenuation factor */
+  double my_a_f;            /* Final attenuation factor */
+  double solid_angle;       /* Solid angle of target as seen from scattering point */
+  double aim_x=0, aim_y=0, aim_z=1;   /* Position of target relative to scattering point */
+  double omega;             /* Energy transfer */
+  double qsquare;           /* Square of the scattering vector */
+  double bose_factor;       /* Calculated value of the Bose factor */
+  int nf, index;            /* Number of allowed final velocities */
+  double vf_list[20];       /* List of allowed final velocities */
+  double J_factor, J0;          /* Jacobian from delta fnc.s in cross section; elastic J-factor */
+  double f1, f2;            /* Probed values of omega_q minus omega */
+  double p_att,p_MC,p_ph1,p_ph2,p_ph3,p_J, p_tot;    /* Temporary multipliers */
+  complex Gn;               /* Phonon structure factor */
+  complex q_dot_e;          /* q dot e */
+  double q_dot_r;           /* q dot r */
+  double qh,qk,ql,qhkx,qhky,qh_Bragg,qk_Bragg,ql_Bragg,qh_res,qk_res,qh_add,qk_add;    /* components of q */
+  double dist_00, dist_01, dist_10, dist_11, dist_min; /* used to calculate nearest Bragg point */
+  int i,j, intersect;
+  
+    if (shape == 0)
+      intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight);
+    else if (shape == 1)
+      intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
+
+  if(intersect)
+  {
+    //    fprintf(stderr,"intersect\n");      
+    if(t0 < 0)
+      ABSORB; /* Neutron came from the sample or begins inside */
+
+    if (verbose>=2) 
+        printf("neutron entered Phonon_BvK_PG \n");
+    /* Neutron enters at t=t0. */
+    dt0 = t3-t0;                /* Time in sample */
+    v_i = sqrt(vx*vx + vy*vy + vz*vz);
+    l_full = v_i * dt0;   /* Length of path through sample if not scattered */
+    dt = rand01()*dt0;    /* Time of scattering (relative to t0) */
+    l_i = v_i*dt;                 /* Penetration in sample at scattering */
+    vx_i=vx;
+    vy_i=vy;
+    vz_i=vz;
+    PROP_DT(dt+t0);             /* Point of scattering */
+
+    aim_x = target_x-x;         /* Vector pointing at target (e.g. analyzer) */
+    aim_y = target_y-y;
+    aim_z = target_z-z;
+
+    if(focus_aw && focus_ah) {
+      randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
+        aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
+    } else if(focus_xw && focus_yh) {
+      randvec_target_rect(&vx, &vy, &vz, &solid_angle,
+        aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP);
+    } else {
+      randvec_target_sphere(&vx,&vy,&vz,&solid_angle,aim_x,aim_y,aim_z, focus_r);
+    }
+    NORM(vx, vy, vz);
+    nf=0;
+    if (mode_input == 12)
+        mode = round(12*rand01()-0.5); // Select the phonon mode
+    else
+        mode = mode_input;             // Use the given mode
+    if ((mode < 0) || (mode > 11))     // Undefined mode specification
+    {
+        printf("mode = %d ",mode);
+        nrerror("illegal value of mode");
+    }
+
+      p_call[0]=-1;   // in the end this will be v_f
+      p_call[1]=v_i;
+      p_call[2]=vx;
+      p_call[3]=vy;
+      p_call[4]=vz;
+      p_call[5]=vx_i;
+      p_call[6]=vy_i;
+      p_call[7]=vz_i;
+      p_call[9]=hh;
+      p_call[10]=kk;
+      p_call[11]=ll;
+      p_call[12] = (double)dispersion;
+      
+      if (dispersion==1)
+      {
+          f1=omega_q(p_call,eigenvectormode);
+          p_call[12]=0.0;
+      }
+
+/*      if (dispersion==2)
+      {
+        for(i=0; i<e_steps_low; i++)  // Loop over v_f values on the energy loss side
+        {
+          v_f = i*v_i/e_steps_low;
+          p_call[0] = v_f;
+          f1=omega_q(p_call,eigenvectormode);
+          qx=V2K*(vx_i-v_f*vx);
+          qy=V2K*(vy_i-v_f*vy);
+          qz=V2K*(vz_i-v_f*vz);
+          h=
+
+          printf("v_f = %g , omega_q = %g \n",v_f,f1);
+        }
+          p_call[12]=0.0;
+      }
+*/
+    #ifndef OPENACC
+    if (verbose >=2) 
+        printf("Call findroots \n");      
+    findroots(0, v_i, v_i+V_HIGH, vf_list, &nf, omega_q, p_call, eigenvectormode, e_steps_low, e_steps_high);
+    if (verbose >=2)
+        printf( "Findroots returned %d roots \n",nf);
+    #else
+    findroots_gpu(0, v_i, v_i+V_HIGH, vf_list, &nf, p_call, eigenvectormode, e_steps_low, e_steps_high);
+    #endif
+
+    if (verbose>=3)
+    {
+        printf("Findroots done, mode %i , last phonon energy %g \n",mode,creal(p_call[8]));
+    }
+    
+    // ________________ ALL FROM HERE IS CALCULATION OF INTENSITY - SHOULD BE REVISITED ___________________
+   
+        index=(int)floor(rand01()*nf); // Select random solution
+        v_f=vf_list[index];
+        p_call[0]=v_f;      // transfer choice of v_f to omega_q
+
+        f1=omega_q(p_call,eigenvectormode);
+        if (verbose >= 2)
+        {
+            printf("Chosen solution is %i; v_f = %g, hw = %g, energy match is %g; eigenvector is: (", index, v_f, creal(p_call[8]), f1);
+            for (j=0; j<DIM; j++)
+                printf(" %g +i %g, ",creal(eigenvectormode[j]),cimag(eigenvectormode[j]));
+            printf(" )\n");
+            printf("--------------------------------------------\n");
+        }
+      p_call[0]=v_f-DV;
+      //      fprintf(stderr,"f1=omega_q(%g+i%g)\n",creal(p_call[0]),cimag(p_call[0]));            
+      f1=omega_q(p_call,eigenvectormode);
+      //      fprintf(stderr,"retur fra f1=omega_q\n");            
+      p_call[0]=v_f+DV;
+      //      fprintf(stderr,"f2=omega_q(%g+i%g)\n",creal(p_call[0]),cimag(p_call[0]));                  
+      f2=omega_q(p_call,eigenvectormode);
+      //      fprintf(stderr,"retur fra f2=omega_q\n");                  
+      J_factor = 1.6022E-22*1E-10*fabs(f2-f1)/(2*DV*V2K);  // unit: meV Å converted to SI units
+      omega=VS2E*(v_i*v_i-v_f*v_f);  // unit: meV
+      vx *= v_f;
+      vy *= v_f;
+      vz *= v_f;
+
+      q_x=q[0]=V2K*(vx_i-vx);
+      q_y=q[1]=V2K*(vy_i-vy);
+      q_z=q[2]=V2K*(vz_i-vz);
+
+      ql = q_z/cstar;
+      qhkx = q_x/astar;     // These are the Cartesian coordinates
+      qhky = q_y/astar;     // These are the Cartesian coordinates
+      qk = 2*qhky/sqrt3;    // transform to hexagonal coordinates
+      qh = qhkx-qhky/sqrt3; // transform to hexagonal coordinates
+
+      // Now find out which Bragg peak is closest in hexagonal plane
+
+      ql_Bragg = round(ql); // We know what to do along the c-axis
+
+      qh_Bragg = floor(qh); // Divide into integer and residual part
+      qh_res = qh-qh_Bragg;
+      qk_Bragg = floor(qk);
+      qk_res = qk-qk_Bragg;
+
+      dist_00 = qh_res*qh_res + qk_res*qk_res + qh_res * qk_res; // Square distance in hexagonal coordinates to (0,0)
+      dist_min = dist_00;
+      qh_add=0;
+      qk_add=0;
+
+      dist_01 = qh_res*qh_res + (qk_res-1)*(qk_res-1) + qh_res * (qk_res-1); // Square distance in hexagonal coordinates to (0,1)
+      if (dist_01 < dist_min)
+      {
+        dist_min = dist_01;
+        qh_add=0;
+        qk_add=1;
+      }
+      dist_10 = (qh_res-1)*(qh_res-1) + qk_res*qk_res + (qh_res-1) * qk_res; // Square distance in hexagonal coordinates to (1,0)
+      if (dist_10 < dist_min)
+      {
+        dist_min = dist_10;
+        qh_add=1;
+        qk_add=0;
+      }
+      dist_11 = (qh_res-1)*(qh_res-1) + (qk_res-1)*(qk_res-1) + (qh_res-1) * (qk_res-1); // Square distance in hexagonal coordinates to (1,1)
+      if (dist_11 < dist_min)
+      {
+        dist_min = dist_11;
+        qh_add=1;
+        qk_add=1;
+      }
+
+      qh_Bragg += qh_add;  // Potentially correct to the right BZ
+      qk_Bragg += qk_add;
+
+        if (verbose >= 3)
+        {
+            printf("q transforms in r.l.u to (h, k, l) = (%g, %g, %g). Nearest Bragg point: (%g, %g, %g)\n",qh,qk,ql,qh_Bragg,qk_Bragg,ql_Bragg);
+            printf("--------------------------------------------\n");
+        }
+
+
+      if (dispersion==1)
+          printf("Dispersion found:  (h,k,l) = (%g, %g, %g), hw = %g \n",q[0]/astar,q[1]/astar,q[2]/cstar,omega);
+      
+      qsquare=1;
+      Gn = 0;
+
+      double q_Bragg[3];  // Bragg point in Å-1; Cartesian coordinates
+      q_Bragg[0]=qh_Bragg*astar+qk_Bragg*astar/2.0;
+      q_Bragg[1]=qk_Bragg*astar*2.0/sqrt3;
+      q_Bragg[2]=ql_Bragg*cstar;
+
+      for (i=0; i<4; i++) // Calculate phonon structure factor; loop over atoms in unit cell
+      {
+          q_dot_e = (complex)0;
+          q_dot_r = 0;
+          for (j=0; j<3; j++) // loop over x,y,z
+          {
+              q_dot_e += (complex)q[j]*eigenvectormode[3*i+j];
+              q_dot_r += q_Bragg[j]*Delta[3*i+j];
+          }
+          if (verbose >= 7)
+          {
+              printf("i = %i , q dot r = %g , q dot e = (%g + i %g) Fn contrib = (%g + i %g)",i,q_dot_r,creal(q_dot_e),cimag(q_dot_e),creal(b_length*q_dot_e*cexp(I*q_dot_r)),cimag(b_length*q_dot_e*cexp(I*q_dot_r)));
+          }
+          Gn += b_length*q_dot_e*cexp(I*q_dot_r); // Unit fm Å-1  ; for the general lattice, this should be divided by sqrt(mass[j])
+      }
+      if (verbose >= 7) 
+          printf("\n");
+
+    if (shape == 0)
+      intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight);
+    else if (shape == 1)
+      intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
+
+      if(!intersect)
+	{
+	  printf("FATAL ERROR: Did not hit sample surface from inside.\n");
+	  exit(1);
+	}
+	
+	if (verbose>=2)
+        printf("neutron left Phonon_BvK_PG, p_init = %g,",p);
+
+      dt = t3;
+      l_o = v_f*dt;
+      my_a_i = V_my_a_v/v_i;
+      my_a_f = V_my_a_v/v_f;
+      bose_factor=nbose(omega,T);
+      J0 = hbar*hbar*v_f*V2K*1E10/mn;  // The J-factor for a flat mode
+
+
+      p_att = exp(-(V_my_s*(l_i+l_o)+my_a_i*l_i+my_a_f*l_o)); /* Absorption factor (units: 1) */
+      p_MC = nf*solid_angle*l_full;     /* Focusing factors; assume random choice of n_f possibilities (units: m) */
+      if (mode_input == 12)
+          p_MC *= 12;  //The factor 12 is due to MC choice between the 12 modes
+      p_ph1 = (v_f/v_i)*DW*bose_factor*V_rho*1E30;   /* Cross section factor 1 (units: m^-3) */
+      p_ph2 = hbar*hbar/(2*(fabs(omega)*1.6022E-22));  /* Jacobian of delta functions in cross section; and constants. Units: J s^2 */
+      p_ph3 = (Gn*conj(Gn))*1E-10/M;  /* Structure factor.   (units: kg^-1) */
+      p_J = J0/J_factor;
+      p_tot = p_att*p_MC*p_ph1*p_ph2*p_ph3*p_J;
+      p *= p_tot;
+
+      if (verbose>=2) {
+       printf("nf = %d, solid_angle = %g, l_full = %g \n",nf,solid_angle,l_full);
+        printf("(p in SI units): p_att = %g, p_MC= %g, p_ph1 = %g, p_ph2 = %g, p_ph3 = %g, p_J = %g, bose_factor = %g, V_rho = %g, J_factor/J0 = %g, G_factor = %g +i %g, conj(G) = %g + I %g, |G|^2 = %g + I %g, ",p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, bose_factor, V_rho, J_factor/J0, creal(Gn), cimag(Gn), creal(conj(Gn)), cimag(conj(Gn)), creal(Gn*conj(Gn)), cimag(Gn*conj(Gn)) );
+        printf(" p_tot = %g, p_final = %g\n",p_tot, p) ;
+        printf("q is (h,k,l) = (%g, %g, %g) \n",q_x/astar,q_y/astar,q_z/cstar);
+      }
+  } /* else transmit: Neutron did not hit the sample */
+%}
+
+MCDISPLAY
+%{
+      if (radius > 0 &&  yheight) {	/* cylinder */
+          cylinder(0,0,0,radius,yheight,0, 0, 1, 0);
+      }
+      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/ESS/KVASIR/KVASIRscript.ipynb b/mcstas-comps/examples/ESS/KVASIR/KVASIRscript.ipynb
new file mode 100755
index 000000000..b2caf7bbe
--- /dev/null
+++ b/mcstas-comps/examples/ESS/KVASIR/KVASIRscript.ipynb
@@ -0,0 +1,414 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import mcstasscript as ms # type: ignore\n",
+    "\n",
+    "configurator = ms.Configurator()\n",
+    "configurator.set_mcrun_path(\"/Applications/McStas-3.2.app/Contents/Resources/mcstas/3.2/bin\")\n",
+    "configurator.set_mcstas_path(\"/Applications/McStas-3.2.app/Contents/Resources/mcstas/3.2\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "KVASIR = ms.McStas_instr(\"KVASIR\")\n",
+    "#EXAMPLE: interated intensity on detector when running vanadium as sample: 0.0268478\n",
+    "t  = KVASIR.add_parameter(\"t\", value=0) # Wedge no. (330 total to cover 0-170deg)\n",
+    "\n",
+    "l_analyzers =[2.5035035035035036, 2.5105105105105103, 2.5175175175175175, 2.5245245245245247, 2.5305305305305303, 2.5375375375375375, 2.5435435435435436, 2.5505505505505504, 2.5565565565565564, 2.5625625625625625, 2.5695695695695697, 2.5755755755755754, 2.5815815815815815, 2.5875875875875876, 2.5925925925925926, 2.5985985985985987, 2.6046046046046047, 2.6096096096096097, 2.6156156156156154, 2.6206206206206204, 2.6266266266266265, 2.6316316316316315, 2.6366366366366365, 2.6416416416416415, 2.6466466466466465, 2.6516516516516515, 2.6566566566566565, 2.6616616616616615, 2.6656656656656654, 2.6706706706706704, 2.674674674674675, 2.67967967967968, 2.6836836836836837, 2.6886886886886887, 2.6926926926926926, 2.6966966966966965, 2.7007007007007005, 2.704704704704705, 2.7087087087087087, 2.7127127127127126, 2.7157157157157155, 2.71971971971972, 2.7237237237237237, 2.7267267267267266, 2.7307307307307305, 2.7337337337337337, 2.736736736736737, 2.7407407407407405, 2.7437437437437437, 2.746746746746747, 2.74974974974975, 2.7527527527527527, 2.754754754754755, 2.7577577577577577, 2.7607607607607605, 2.7637637637637638, 2.7657657657657655, 2.7687687687687688, 2.7707707707707705, 2.7727727727727727, 2.7757757757757755, 2.7777777777777777, 2.77977977977978, 2.781781781781782, 2.7837837837837838, 2.7857857857857855, 2.7877877877877877, 2.78978978978979, 2.7907907907907905, 2.7927927927927927, 2.793793793793794, 2.7957957957957955, 2.796796796796797, 2.798798798798799, 2.7997997997998, 2.8008008008008005, 2.801801801801802, 2.8028028028028027, 2.804804804804805, 2.804804804804805] #array containg nummerical solution for distance from sample to n'th analyzer in wedge\n",
+    "\n",
+    "N_analyzer =60                                                      #1 total number of analyzers\n",
+    "width = 0.0127                                                      #m width of He3 tubes, height of analyzer crystals\n",
+    "dt = 2*np.arcsin(width/(2*2.5))*180/np.pi                           #deg vertical angle between each analyzer crystal seen from sample position\n",
+    "theta =  np.linspace(-dt/2,-dt*N_analyzer-dt/2, N_analyzer)         #deg array containg rotations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# PARAMETERS\n",
+    "width= KVASIR.add_declare_var(\"double\", \"width\", value=0.0127)                              #m width of He3 tubes, height of analyzer crystals\n",
+    "d_PG= KVASIR.add_declare_var(\"double\", \"d_PG\", value=3.355)                                 #Å lattice spacing og pyrolytic graphite\n",
+    "u = KVASIR.add_declare_var(\"double\",\"u\", value=1e-5)                                        #1 small number\n",
+    "h_D = KVASIR.add_declare_var(\"double\", \"h_D\", value=0.4)                                    #m heigt of detector \n",
+    "\n",
+    "Ef = KVASIR.add_declare_var(\"double\",\"Ef\", value=1.95)                                      #meV focusing energy                       \n",
+    "d_AD = KVASIR.add_declare_var(\"double\",\"d_AD\", value=1.25)                                  #m distance from analyzer to detactor\n",
+    "d_SA = KVASIR.add_declare_var(\"double\",\"d_SA\", value=2.5)                                   #m distance from sample to central analyzer\n",
+    "d_GS = KVASIR.add_declare_var(\"double\",\"d_GS\", value=0.58)                                  #m distance from guide end to sample position\n",
+    "\n",
+    "lamda_f = KVASIR.add_declare_var(\"double\", \"lambda_f\")                                      #Å final wavelenght\n",
+    "KVASIR.append_initialize(\"lambda_f=1/(0.11056*sqrt(Ef));\")\n",
+    "\n",
+    "#ANALYZER CALCULATION\n",
+    "tA = KVASIR.add_declare_var(\"double\", \"tA\")                                                 #deg Bragg angle of analyzer\n",
+    "KVASIR.append_initialize(\"tA=asin(lambda_f/(2*d_PG))*RAD2DEG;\")\n",
+    "\n",
+    "dtS = KVASIR.add_declare_var(\"double\", \"dtS\", value=[*theta], array = len(theta))           #deg array containg rotations\n",
+    "\n",
+    "len = KVASIR.add_declare_var(\"double\", \"len\", value=l_analyzers, array = len(l_analyzers))  #array containg nummerical solution for distance from sample to n'th analyzer in wedge\n",
+    "\n",
+    "#DETECTOR CALCULATION\n",
+    "d_SDz = KVASIR.add_declare_var(\"double\", \"d_SDz\")                                           #m vetical distance from scattering plane to center of detector\n",
+    "KVASIR.append_initialize(\"d_SDz=d_SA+d_AD*cos(2*tA*DEG2RAD);\")\n",
+    "\n",
+    "d_SDy = KVASIR.add_declare_var(\"double\", \"d_SDy\")                                           #m distance from sample to detector in scattering plane\n",
+    "KVASIR.append_initialize(\"d_SDy=d_AD*sin(2*tA*DEG2RAD);\")\n",
+    "\n",
+    "tDy= KVASIR.add_declare_var(\"double\", \"tDy\")                                                #deg horizontal angle between wedges                                          \n",
+    "KVASIR.append_initialize(\"tDy=asin(width/(d_SDz))*RAD2DEG;\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## SOURCE OPTION 1: DUMMY SOURCE OPTION\n",
+    "# source = KVASIR.add_component(f\"source\", \"Source_Maxwell_3\", AT=[0,0,0], RELATIVE=\"ABSOLUTE\")\n",
+    "# source.xwidth=0.01\n",
+    "# source.yheight=0.01\n",
+    "# source.Lmin= \"lambda_f-D_lambda/2\"\n",
+    "# source.Lmax= \"lambda_f+D_lambda/2\"\n",
+    "# source.dist= d_SA\n",
+    "# source.focus_yh= f\"width*2*{N_analyzer}\"\n",
+    "# source.focus_xw=f\"width*3\"\n",
+    "# source.T1=300\n",
+    "# source.T2=300\n",
+    "# source.T3=300\n",
+    "# source.I1=1E15\n",
+    "# source.I2=1E15\n",
+    "# source.I3=1E15\n",
+    "\n",
+    "## SOURCE OPTION 2: IMPORT MCPL FILE FROM PRIMARY SPECTROMETER\n",
+    "source = KVASIR.add_component(f\"source\", \"MCPL_input\", AT=[0,0,0], RELATIVE=\"ABSOLUTE\")\n",
+    "source.filename= '\"Virtual_output_endOfGuide_narrow_100mus.mcpl\"' #EXAMPLE MCPL: BIFROST GUIDE, 100mus PSC opening time, wavelength band is centered at 1.95meV and narrowed to speed up simulations\n",
+    "source.pos_smear = 0.005\n",
+    "source.dir_smear = 0.05\n",
+    "\n",
+    "#ARM DEFINES SAMPLE POSITION (0.58m is distance from BIFROST guide opening to sample position)\n",
+    "Arm_sample = KVASIR.add_component(f\"Arm_sample\", \"Arm\", AT=[0,0,d_GS], ROTATED = [0,f\"tDy*t\",0], RELATIVE=\"ABSOLUTE\")\n",
+    "\n",
+    "## SAMPLE OPTION 1: VANADIUM\n",
+    "sample = KVASIR.add_component(f\"sample\", \"Incoherent\", AT=[0,0,0], RELATIVE=\"Arm_sample\")\n",
+    "sample.radius=0.005\n",
+    "sample.yheight=0.01\n",
+    "sample.thickness=0.001\n",
+    "sample.target_index = 2\n",
+    "sample.focus_ah =0.5\n",
+    "sample.focus_aw = 20 \n",
+    "\n",
+    "## SAMPLE OPTION 1: POWDER SAMPLE\n",
+    "# sample = KVASIR.add_component(f\"sample\", \"PowderN\", AT=[0,0,0], RELATIVE=\"Arm_sample\")\n",
+    "# sample.radius=0.005\n",
+    "# sample.yheight=0.01\n",
+    "# sample.thickness=0.001\n",
+    "# sample.reflections=\"Na2Ca3Al2F14.laz\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for n in range(0,N_analyzer,2):                                                                  #loops placing analyzers (every second analyzer is placed to avoid shading)\n",
+    "\n",
+    "    az = KVASIR.add_declare_var(\"double\", f\"az_{n}\")                                             #m z-component of n'th analyzer relative to sample position\n",
+    "    KVASIR.append_initialize(f\"az_{n} = len[{n}]*cos(dtS[{n}]*DEG2RAD)-d_SDz;\")\n",
+    "    \n",
+    "    ay = KVASIR.add_declare_var(\"double\", f\"ay_{n}\")                                             #m y-component of n'th analyzer relative to sample position\n",
+    "    KVASIR.append_initialize(f\"ay_{n} = len[{n}]*sin(dtS[{n}]*DEG2RAD)-d_SDy;\")\n",
+    "\n",
+    "    bz = KVASIR.add_declare_var(\"double\", f\"bz_{n}\")                                             #m z-component of detector relative to n'th analyzer\n",
+    "    KVASIR.append_initialize(f\"bz_{n} = len[{n}]*cos(dtS[{n}]*DEG2RAD);\")                        \n",
+    "    \n",
+    "    by = KVASIR.add_declare_var(\"double\", f\"by_{n}\")                                             #m y-component of detector relative to n'th analyzer\n",
+    "    KVASIR.append_initialize(f\"by_{n} = len[{n}]*sin(dtS[{n}]*DEG2RAD);\")\n",
+    "    \n",
+    "    dtA = KVASIR.add_declare_var(\"double\", f\"dtA_{n}\")                                           #deg scattering angle of n'th analyzer\n",
+    "    KVASIR.append_initialize(f\"dtA_{n} = acos( (az_{n}*bz_{n} + ay_{n}*by_{n}) / (len[{n}]* sqrt(az_{n}*az_{n} + ay_{n}*ay_{n}) ) )*RAD2DEG;\")\n",
+    "\n",
+    "    # GENERARING UPPER ANALYZER WEDGE\n",
+    "    Arm = KVASIR.add_component(f\"Arm_{n}\", \"Arm\", AT=[0,0,0], ROTATED=[f\"-dtS[{n}]\",0,0], RELATIVE=\"Arm_sample\") \n",
+    "    Analyser = KVASIR.add_component(f\"Analyzer_{n}\",\"Monochromator_flat\", AT = [0,0,f\"len[{n}]\"], ROTATED=[f\"dtA_{n}/2\",90,0], RELATIVE=f\"Arm_{n}\")\n",
+    "    Analyser.zwidth = f\"2*len[{n}]*sin(tDy/2*DEG2RAD)\"\n",
+    "    Analyser.yheight = width\n",
+    "    Analyser.mosaich = 90\n",
+    "    Analyser.mosaicv = 90\n",
+    "    Analyser.Q = 1.87278\n",
+    "    Analyser.r0 = 0.8\n",
+    "\n",
+    "    ## GENERARING LOWER ANALYZER WEDGE\n",
+    "    # Arm = KVASIR.add_component(f\"Armm_{n}\", \"Arm\", AT=[0,0,0], ROTATED=[f\"dtS[{n}]\",0,0], RELATIVE=\"Arm_sample\")\n",
+    "    # Analyser = KVASIR.add_component(f\"Analyzerm_{n}\",\"Monochromator_flat\", AT = [0,0,f\"len[{n}]\"], ROTATED=[f\"-dtA_{n}/2\",90,0], RELATIVE=f\"Armm_{n}\" )\n",
+    "    # Analyser.zwidth = f\"2*len[{n}]*sin(tDy/2*DEG2RAD)\"\n",
+    "    # Analyser.yheight = width\n",
+    "    # Analyser.mosaich = 90\n",
+    "    # Analyser.mosaicv = 90\n",
+    "    # Analyser.Q = 1.87278\n",
+    "    # Analyser.r0 = 0.8"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N_detectors  = 1#24#230                                                                         #Number of detectors, centered around wedge (24 is recomended for single wedge to capture mosaic spread)\n",
+    "\n",
+    "for n in range(N_detectors):\n",
+    "    ## UPPER ToF MONITOR SECTION\n",
+    "    KVASIR.add_component(f\"arm_detector_{n}\", \"Arm\", AT = [0,0,0], RELATIVE=\"Arm_sample\", ROTATED=[0,f\"tDy*{n}-tDy*{N_detectors/2}\",0])\n",
+    "    det_He3 = KVASIR.add_component(f\"det_He3_{n}_ToF\", \"Monitor_nD\", AT = [0,\"d_SDy\", d_SDz], AT_RELATIVE=\"PREVIOUS\", ROTATED=[0,f\"tDy*{n}-tDy*{N_detectors/2}\",0], ROTATED_RELATIVE=\"Arm_sample\")\n",
+    "    det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=80, t, limits = [0.26,0.28], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    det_He3.xwidth = width\n",
+    "    det_He3.yheight = \"h_D\"\n",
+    "    det_He3.filename = f'\"Detector_{n}_ToF.dat\"'\n",
+    "\n",
+    "    ## UPPER ENERGY MONITOR SECTION\n",
+    "    # det_He3 = KVASIR.add_component(f\"det_He3_{n}\", \"Monitor_nD\", AT = [0,0, \"-u\"], AT_RELATIVE=\"PREVIOUS\")\n",
+    "    # det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=80, energy, limits = [1.9,2.1], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    # det_He3.xwidth = width\n",
+    "    # det_He3.yheight = \"h_D\"\n",
+    "    # det_He3.filename = f'\"Detector_{n}.dat\"'\n",
+    "    # # det_He3.set_WHEN(\"wedge == 1\")\n",
+    "\n",
+    "\n",
+    "\n",
+    "    ## LOWER ToF MONITOR SECTION\n",
+    "    # KVASIR.add_component(f\"arm_detector_{n}_lower\", \"Arm\", AT = [0,0,0], RELATIVE=\"Arm_sample\", ROTATED=[0,f\"tDy*{n}\",0])\n",
+    "    # det_He3 = KVASIR.add_component(f\"det_He3_{n}_lower\", \"Monitor_nD\", AT = [0,\"-d_SDy\", d_SDz], AT_RELATIVE=\"PREVIOUS\", ROTATED=[0,f\"tDy*{n}\",0], ROTATED_RELATIVE=\"Arm_sample\")\n",
+    "    # det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=300, t, limits = [0.27,0.278], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    # det_He3.xwidth = width\n",
+    "    # det_He3.yheight = \"h_D\"\n",
+    "    # det_He3.filename = f'\"Detector_{n}_ToF_lower.dat\"'\n",
+    "\n",
+    "    ## LOWER ENERGY MONITOR SECTION\n",
+    "    # det_He3 = KVASIR.add_component(f\"det_He3_{n}_ToF_lower\", \"Monitor_nD\", AT = [0,0, \"-u\"], AT_RELATIVE=\"PREVIOUS\")\n",
+    "    # det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=300, energy, limits = [1.9,2.1], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    # det_He3.xwidth = width\n",
+    "    # det_He3.yheight = \"h_D\"\n",
+    "    # det_He3.filename = f'\"Detector_{n}_lower.dat\"'\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#KVASIR.show_components()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "---- Found 5 places in McStas output with keyword 'error'. \n",
+      "\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "--------------------------------------------------------------------------\n",
+      "Primary job  terminated normally, but 1 process returned\n",
+      "a non-zero exit code. Per user-direction, the job has been aborted.\n",
+      "--------------------------------------------------------------------------\n",
+      "--------------------------------------------------------------------------\n",
+      "mpirun detected that one or more processes exited with non-zero status, thus causing\n",
+      "the job to be terminated. The first process to do so was:\n",
+      "\n",
+      "  Process name: [[57467,1],5]\n",
+      "----------------------------------------------------------------------\n",
+      "\n",
+      "loading system configuration\n",
+      "loading user configuration from /Users/amaliedavidsen/.mcstas/3.2/mccode_config.json\n",
+      "INFO: Using directory: \"/Volumes/T7/KVASIR/KVASIR_for_repo/KVASIR\"\n",
+      "INFO: Regenerating c-file: KVASIR.c\n",
+      "CFLAGS= -Wl,-rpath,CMD(mcpl-config --show libdir) -LCMD(mcpl-config --show libdir) -lmcpl -ICMD(mcpl-config --show includedir)\n",
+      "      \n",
+      "-----------------------------------------------------------\n",
+      "\n",
+      "Generating single GPU kernel or single CPU section layout: \n",
+      "\n",
+      "-----------------------------------------------------------\n",
+      "\n",
+      "Generating GPU/CPU -DFUNNEL layout:\n",
+      "\n",
+      "-----------------------------------------------------------\n",
+      "INFO: Recompiling: ./KVASIR.out\n",
+      "./KVASIR.c:3393:1: warning: non-void function does not return a value in all control paths [-Wreturn-type]\n",
+      " 3393 | } /* siminfo_init */\n",
+      "      | ^\n",
+      "./KVASIR.c:14678:37: warning: format specifies type 'unsigned long' but the argument has type 'unsigned long long' [-Wformat]\n",
+      " 14674 |     MPI_MASTER(\n",
+      "       |     ~~~~~~~~~~~\n",
+      " 14675 | #endif\n",
+      "       | ~~~~~~\n",
+      " 14676 |            fprintf(stdout, \"\\n\\n Warning: You are using MCPL_input with a repeat_count of %lu:\\n - Minimum neutron count requested is %lu x %lu <= %lu\",\n",
+      "       |            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      "       |                                                                                                                                                    %llu\n",
+      " 14677 |                (long unsigned)repeat_count,(long unsigned)nparticles,\n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14678 |                (long unsigned)repeat_count,(long unsigned)repeat_cnt*nparticles); \n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14679 | #if defined (USE_MPI)\n",
+      "       | ~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14680 |   fprintf(stdout, \" x %i MPI nodes = %lu neutrons total\\n\",\n",
+      "       |   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14681 |     mpi_node_count,(long unsigned)mpi_node_count*repeat_cnt*nparticles);\n",
+      "       |     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14682 |      );\n",
+      "       |      ~\n",
+      "./KVASIR.c:456:5: note: expanded from macro 'MPI_MASTER'\n",
+      "  456 |   { statement; } \\\n",
+      "      |     ^~~~~~~~~\n",
+      "./KVASIR.c:14681:20: warning: format specifies type 'unsigned long' but the argument has type 'unsigned long long' [-Wformat]\n",
+      " 14674 |     MPI_MASTER(\n",
+      "       |     ~~~~~~~~~~~\n",
+      " 14675 | #endif\n",
+      "       | ~~~~~~\n",
+      " 14676 |            fprintf(stdout, \"\\n\\n Warning: You are using MCPL_input with a repeat_count of %lu:\\n - Minimum neutron count requested is %lu x %lu <= %lu\",\n",
+      "       |            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14677 |                (long unsigned)repeat_count,(long unsigned)nparticles,\n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14678 |                (long unsigned)repeat_count,(long unsigned)repeat_cnt*nparticles); \n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14679 | #if defined (USE_MPI)\n",
+      "       | ~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14680 |   fprintf(stdout, \" x %i MPI nodes = %lu neutrons total\\n\",\n",
+      "       |   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      "       |                                      %llu\n",
+      " 14681 |     mpi_node_count,(long unsigned)mpi_node_count*repeat_cnt*nparticles);\n",
+      "       |     ~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14682 |      );\n",
+      "       |      ~\n",
+      "./KVASIR.c:456:5: note: expanded from macro 'MPI_MASTER'\n",
+      "  456 |   { statement; } \\\n",
+      "      |     ^~~~~~~~~\n",
+      "./KVASIR.c:14773:36: warning: format specifies type 'unsigned long' but the argument has type 'unsigned long long' [-Wformat]\n",
+      " 14771 |       fprintf(stdout, \"\\n\\n Warning: You are using MCPL_input with a repeat_count of %lu:\\n - Neutron count requested is %lu x %lu <= %lu\",\n",
+      "       |                                                                                                                                       ~~~\n",
+      "       |                                                                                                                                       %llu\n",
+      " 14772 |               (long unsigned)repeat_count,(long unsigned)read_neutrons,\n",
+      " 14773 |               (long unsigned)repeat_count,(long unsigned)repeat_cnt*read_neutrons);\n",
+      "       |                                           ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      "4 warnings generated.\n",
+      "[cctools-port]: generating fake signature for './KVASIR.out'\n",
+      "INFO: ===\n",
+      "Simulation 'KVASIR' (KVASIR.instr): running on 6 nodes (master is 'Amalies-MacBook-Air.local', MPI version 3.1).\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "--------------------------------------------------------------------------\n",
+      "Primary job  terminated normally, but 1 process returned\n",
+      "a non-zero exit code. Per user-direction, the job has been aborted.\n",
+      "--------------------------------------------------------------------------\n",
+      "--------------------------------------------------------------------------\n",
+      "mpirun detected that one or more processes exited with non-zero status, thus causing\n",
+      "the job to be terminated. The first process to do so was:\n",
+      "\n",
+      "  Process name: [[57467,1],5]\n",
+      "  Exit code:    1\n",
+      "--------------------------------------------------------------------------\n",
+      "\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "KVASIR.set_parameters(t=50)\n",
+    "KVASIR.settings(ncount=1e6, mpi=6)\n",
+    "data = KVASIR.backengine()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "No data to plot\n"
+     ]
+    }
+   ],
+   "source": [
+    "ms.make_sub_plot(data, log=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#KVASIR.show_instrument(format=\"window\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "mcstas-dev",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/mcstas-comps/examples/ESS/KVASIR/KVASIRscript_submit.ipynb b/mcstas-comps/examples/ESS/KVASIR/KVASIRscript_submit.ipynb
new file mode 100755
index 000000000..70055f6d6
--- /dev/null
+++ b/mcstas-comps/examples/ESS/KVASIR/KVASIRscript_submit.ipynb
@@ -0,0 +1,413 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import mcstasscript as ms # type: ignore\n",
+    "\n",
+    "configurator = ms.Configurator()\n",
+    "configurator.set_mcrun_path(\"/Applications/McStas-3.2.app/Contents/Resources/mcstas/3.2/bin\")\n",
+    "configurator.set_mcstas_path(\"/Applications/McStas-3.2.app/Contents/Resources/mcstas/3.2\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "KVASIR = ms.McStas_instr(\"KVASIR\")\n",
+    "t  = KVASIR.add_parameter(\"t\", value=0) # Wedge no. (330 total to cover 0-170deg)\n",
+    "\n",
+    "l_analyzers =[2.5035035035035036, 2.5105105105105103, 2.5175175175175175, 2.5245245245245247, 2.5305305305305303, 2.5375375375375375, 2.5435435435435436, 2.5505505505505504, 2.5565565565565564, 2.5625625625625625, 2.5695695695695697, 2.5755755755755754, 2.5815815815815815, 2.5875875875875876, 2.5925925925925926, 2.5985985985985987, 2.6046046046046047, 2.6096096096096097, 2.6156156156156154, 2.6206206206206204, 2.6266266266266265, 2.6316316316316315, 2.6366366366366365, 2.6416416416416415, 2.6466466466466465, 2.6516516516516515, 2.6566566566566565, 2.6616616616616615, 2.6656656656656654, 2.6706706706706704, 2.674674674674675, 2.67967967967968, 2.6836836836836837, 2.6886886886886887, 2.6926926926926926, 2.6966966966966965, 2.7007007007007005, 2.704704704704705, 2.7087087087087087, 2.7127127127127126, 2.7157157157157155, 2.71971971971972, 2.7237237237237237, 2.7267267267267266, 2.7307307307307305, 2.7337337337337337, 2.736736736736737, 2.7407407407407405, 2.7437437437437437, 2.746746746746747, 2.74974974974975, 2.7527527527527527, 2.754754754754755, 2.7577577577577577, 2.7607607607607605, 2.7637637637637638, 2.7657657657657655, 2.7687687687687688, 2.7707707707707705, 2.7727727727727727, 2.7757757757757755, 2.7777777777777777, 2.77977977977978, 2.781781781781782, 2.7837837837837838, 2.7857857857857855, 2.7877877877877877, 2.78978978978979, 2.7907907907907905, 2.7927927927927927, 2.793793793793794, 2.7957957957957955, 2.796796796796797, 2.798798798798799, 2.7997997997998, 2.8008008008008005, 2.801801801801802, 2.8028028028028027, 2.804804804804805, 2.804804804804805] #array containg nummerical solution for distance from sample to n'th analyzer in wedge\n",
+    "\n",
+    "N_analyzer =60                                                      #1 total number of analyzers\n",
+    "width = 0.0127                                                      #m width of He3 tubes, height of analyzer crystals\n",
+    "dt = 2*np.arcsin(width/(2*2.5))*180/np.pi                           #deg vertical angle between each analyzer crystal seen from sample position\n",
+    "theta =  np.linspace(-dt/2,-dt*N_analyzer-dt/2, N_analyzer)         #deg array containg rotations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# PARAMETERS\n",
+    "width= KVASIR.add_declare_var(\"double\", \"width\", value=0.0127)                              #m width of He3 tubes, height of analyzer crystals\n",
+    "d_PG= KVASIR.add_declare_var(\"double\", \"d_PG\", value=3.355)                                 #Å lattice spacing og pyrolytic graphite\n",
+    "u = KVASIR.add_declare_var(\"double\",\"u\", value=1e-5)                                        #1 small number\n",
+    "h_D = KVASIR.add_declare_var(\"double\", \"h_D\", value=0.4)                                    #m heigt of detector \n",
+    "\n",
+    "Ef = KVASIR.add_declare_var(\"double\",\"Ef\", value=1.95)                                      #meV focusing energy                       \n",
+    "d_AD = KVASIR.add_declare_var(\"double\",\"d_AD\", value=1.25)                                  #m distance from analyzer to detactor\n",
+    "d_SA = KVASIR.add_declare_var(\"double\",\"d_SA\", value=2.5)                                   #m distance from sample to central analyzer\n",
+    "d_GS = KVASIR.add_declare_var(\"double\",\"d_GS\", value=0.58)                                  #m distance from guide end to sample position\n",
+    "\n",
+    "lamda_f = KVASIR.add_declare_var(\"double\", \"lambda_f\")                                      #Å final wavelenght\n",
+    "KVASIR.append_initialize(\"lambda_f=1/(0.11056*sqrt(Ef));\")\n",
+    "\n",
+    "#ANALYZER CALCULATION\n",
+    "tA = KVASIR.add_declare_var(\"double\", \"tA\")                                                 #deg Bragg angle of analyzer\n",
+    "KVASIR.append_initialize(\"tA=asin(lambda_f/(2*d_PG))*RAD2DEG;\")\n",
+    "\n",
+    "dtS = KVASIR.add_declare_var(\"double\", \"dtS\", value=[*theta], array = len(theta))           #deg array containg rotations\n",
+    "\n",
+    "len = KVASIR.add_declare_var(\"double\", \"len\", value=l_analyzers, array = len(l_analyzers))  #array containg nummerical solution for distance from sample to n'th analyzer in wedge\n",
+    "\n",
+    "#DETECTOR CALCULATION\n",
+    "d_SDz = KVASIR.add_declare_var(\"double\", \"d_SDz\")                                           #m vetical distance from scattering plane to center of detector\n",
+    "KVASIR.append_initialize(\"d_SDz=d_SA+d_AD*cos(2*tA*DEG2RAD);\")\n",
+    "\n",
+    "d_SDy = KVASIR.add_declare_var(\"double\", \"d_SDy\")                                           #m distance from sample to detector in scattering plane\n",
+    "KVASIR.append_initialize(\"d_SDy=d_AD*sin(2*tA*DEG2RAD);\")\n",
+    "\n",
+    "tDy= KVASIR.add_declare_var(\"double\", \"tDy\")                                                #deg horizontal angle between wedges                                          \n",
+    "KVASIR.append_initialize(\"tDy=asin(width/(d_SDz))*RAD2DEG;\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## SOURCE OPTION 1: DUMMY SOURCE OPTION\n",
+    "# source = KVASIR.add_component(f\"source\", \"Source_Maxwell_3\", AT=[0,0,0], RELATIVE=\"ABSOLUTE\")\n",
+    "# source.xwidth=0.01\n",
+    "# source.yheight=0.01\n",
+    "# source.Lmin= \"lambda_f-D_lambda/2\"\n",
+    "# source.Lmax= \"lambda_f+D_lambda/2\"\n",
+    "# source.dist= d_SA\n",
+    "# source.focus_yh= f\"width*2*{N_analyzer}\"\n",
+    "# source.focus_xw=f\"width*3\"\n",
+    "# source.T1=300\n",
+    "# source.T2=300\n",
+    "# source.T3=300\n",
+    "# source.I1=1E15\n",
+    "# source.I2=1E15\n",
+    "# source.I3=1E15\n",
+    "\n",
+    "## SOURCE OPTION 2: IMPORT MCPL FILE FROM PRIMARY SPECTROMETER\n",
+    "source = KVASIR.add_component(f\"source\", \"MCPL_input\", AT=[0,0,0], RELATIVE=\"ABSOLUTE\")\n",
+    "source.filename= '\"Virtual_output_endOfGuide_narrow_100mus.mcpl\"' #EXAMPLE MCPL: BIFROST GUIDE, 100mus PSC opening time, wavelength band is centered at 1.95meV and narrowed to speed up simulations\n",
+    "source.pos_smear = 0.005\n",
+    "source.dir_smear = 0.05\n",
+    "\n",
+    "#ARM DEFINES SAMPLE POSITION (0.58m is distance from BIFROST guide opening to sample position)\n",
+    "Arm_sample = KVASIR.add_component(f\"Arm_sample\", \"Arm\", AT=[0,0,d_GS], ROTATED = [0,f\"tDy*t\",0], RELATIVE=\"ABSOLUTE\")\n",
+    "\n",
+    "## SAMPLE OPTION 1: VANADIUM\n",
+    "sample = KVASIR.add_component(f\"sample\", \"Incoherent\", AT=[0,0,0], RELATIVE=\"Arm_sample\")\n",
+    "sample.radius=0.005\n",
+    "sample.yheight=0.01\n",
+    "sample.thickness=0.001\n",
+    "sample.target_index = 2\n",
+    "sample.focus_ah =0.5\n",
+    "sample.focus_aw = 20 \n",
+    "\n",
+    "## SAMPLE OPTION 1: POWDER SAMPLE\n",
+    "# sample = KVASIR.add_component(f\"sample\", \"PowderN\", AT=[0,0,0], RELATIVE=\"Arm_sample\")\n",
+    "# sample.radius=0.005\n",
+    "# sample.yheight=0.01\n",
+    "# sample.thickness=0.001\n",
+    "# sample.reflections=\"Na2Ca3Al2F14.laz\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for n in range(0,N_analyzer,2):                                                                  #loops placing analyzers (every second analyzer is placed to avoid shading)\n",
+    "\n",
+    "    az = KVASIR.add_declare_var(\"double\", f\"az_{n}\")                                             #m z-component of n'th analyzer relative to sample position\n",
+    "    KVASIR.append_initialize(f\"az_{n} = len[{n}]*cos(dtS[{n}]*DEG2RAD)-d_SDz;\")\n",
+    "    \n",
+    "    ay = KVASIR.add_declare_var(\"double\", f\"ay_{n}\")                                             #m y-component of n'th analyzer relative to sample position\n",
+    "    KVASIR.append_initialize(f\"ay_{n} = len[{n}]*sin(dtS[{n}]*DEG2RAD)-d_SDy;\")\n",
+    "\n",
+    "    bz = KVASIR.add_declare_var(\"double\", f\"bz_{n}\")                                             #m z-component of detector relative to n'th analyzer\n",
+    "    KVASIR.append_initialize(f\"bz_{n} = len[{n}]*cos(dtS[{n}]*DEG2RAD);\")                        \n",
+    "    \n",
+    "    by = KVASIR.add_declare_var(\"double\", f\"by_{n}\")                                             #m y-component of detector relative to n'th analyzer\n",
+    "    KVASIR.append_initialize(f\"by_{n} = len[{n}]*sin(dtS[{n}]*DEG2RAD);\")\n",
+    "    \n",
+    "    dtA = KVASIR.add_declare_var(\"double\", f\"dtA_{n}\")                                           #deg scattering angle of n'th analyzer\n",
+    "    KVASIR.append_initialize(f\"dtA_{n} = acos( (az_{n}*bz_{n} + ay_{n}*by_{n}) / (len[{n}]* sqrt(az_{n}*az_{n} + ay_{n}*ay_{n}) ) )*RAD2DEG;\")\n",
+    "\n",
+    "    # GENERARING UPPER ANALYZER WEDGE\n",
+    "    Arm = KVASIR.add_component(f\"Arm_{n}\", \"Arm\", AT=[0,0,0], ROTATED=[f\"-dtS[{n}]\",0,0], RELATIVE=\"Arm_sample\") \n",
+    "    Analyser = KVASIR.add_component(f\"Analyzer_{n}\",\"Monochromator_flat\", AT = [0,0,f\"len[{n}]\"], ROTATED=[f\"dtA_{n}/2\",90,0], RELATIVE=f\"Arm_{n}\")\n",
+    "    Analyser.zwidth = f\"2*len[{n}]*sin(tDy/2*DEG2RAD)\"\n",
+    "    Analyser.yheight = width\n",
+    "    Analyser.mosaich = 90\n",
+    "    Analyser.mosaicv = 90\n",
+    "    Analyser.Q = 1.87278\n",
+    "    Analyser.r0 = 0.8\n",
+    "\n",
+    "    ## GENERARING LOWER ANALYZER WEDGE\n",
+    "    # Arm = KVASIR.add_component(f\"Armm_{n}\", \"Arm\", AT=[0,0,0], ROTATED=[f\"dtS[{n}]\",0,0], RELATIVE=\"Arm_sample\")\n",
+    "    # Analyser = KVASIR.add_component(f\"Analyzerm_{n}\",\"Monochromator_flat\", AT = [0,0,f\"len[{n}]\"], ROTATED=[f\"-dtA_{n}/2\",90,0], RELATIVE=f\"Armm_{n}\" )\n",
+    "    # Analyser.zwidth = f\"2*len[{n}]*sin(tDy/2*DEG2RAD)\"\n",
+    "    # Analyser.yheight = width\n",
+    "    # Analyser.mosaich = 90\n",
+    "    # Analyser.mosaicv = 90\n",
+    "    # Analyser.Q = 1.87278\n",
+    "    # Analyser.r0 = 0.8"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N_detectors  = 1#24#230                                                                         #Number of detectors, centered around wedge (24 is recomended for single wedge to capture mosaic spread)\n",
+    "\n",
+    "for n in range(N_detectors):\n",
+    "    ## UPPER ToF MONITOR SECTION\n",
+    "    KVASIR.add_component(f\"arm_detector_{n}\", \"Arm\", AT = [0,0,0], RELATIVE=\"Arm_sample\", ROTATED=[0,f\"tDy*{n}-tDy*{N_detectors/2}\",0])\n",
+    "    det_He3 = KVASIR.add_component(f\"det_He3_{n}_ToF\", \"Monitor_nD\", AT = [0,\"d_SDy\", d_SDz], AT_RELATIVE=\"PREVIOUS\", ROTATED=[0,f\"tDy*{n}-tDy*{N_detectors/2}\",0], ROTATED_RELATIVE=\"Arm_sample\")\n",
+    "    det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=80, t, limits = [0.26,0.28], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    det_He3.xwidth = width\n",
+    "    det_He3.yheight = \"h_D\"\n",
+    "    det_He3.filename = f'\"Detector_{n}_ToF.dat\"'\n",
+    "\n",
+    "    ## UPPER ENERGY MONITOR SECTION\n",
+    "    # det_He3 = KVASIR.add_component(f\"det_He3_{n}\", \"Monitor_nD\", AT = [0,0, \"-u\"], AT_RELATIVE=\"PREVIOUS\")\n",
+    "    # det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=80, energy, limits = [1.9,2.1], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    # det_He3.xwidth = width\n",
+    "    # det_He3.yheight = \"h_D\"\n",
+    "    # det_He3.filename = f'\"Detector_{n}.dat\"'\n",
+    "    # # det_He3.set_WHEN(\"wedge == 1\")\n",
+    "\n",
+    "\n",
+    "\n",
+    "    ## LOWER ToF MONITOR SECTION\n",
+    "    # KVASIR.add_component(f\"arm_detector_{n}_lower\", \"Arm\", AT = [0,0,0], RELATIVE=\"Arm_sample\", ROTATED=[0,f\"tDy*{n}\",0])\n",
+    "    # det_He3 = KVASIR.add_component(f\"det_He3_{n}_lower\", \"Monitor_nD\", AT = [0,\"-d_SDy\", d_SDz], AT_RELATIVE=\"PREVIOUS\", ROTATED=[0,f\"tDy*{n}\",0], ROTATED_RELATIVE=\"Arm_sample\")\n",
+    "    # det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=300, t, limits = [0.27,0.278], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    # det_He3.xwidth = width\n",
+    "    # det_He3.yheight = \"h_D\"\n",
+    "    # det_He3.filename = f'\"Detector_{n}_ToF_lower.dat\"'\n",
+    "\n",
+    "    ## LOWER ENERGY MONITOR SECTION\n",
+    "    # det_He3 = KVASIR.add_component(f\"det_He3_{n}_ToF_lower\", \"Monitor_nD\", AT = [0,0, \"-u\"], AT_RELATIVE=\"PREVIOUS\")\n",
+    "    # det_He3.options = '\"cylinder, y, limits=[-0.2,0.2], bins=300, energy, limits = [1.9,2.1], bins=500\"' # t, limits = [0.0069,0.00695] energy, limits = [1.82,1.84]\n",
+    "    # det_He3.xwidth = width\n",
+    "    # det_He3.yheight = \"h_D\"\n",
+    "    # det_He3.filename = f'\"Detector_{n}_lower.dat\"'\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#KVASIR.show_components()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "---- Found 5 places in McStas output with keyword 'error'. \n",
+      "\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "----------------------------------------------------------------------\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "--------------------------------------------------------------------------\n",
+      "Primary job  terminated normally, but 1 process returned\n",
+      "a non-zero exit code. Per user-direction, the job has been aborted.\n",
+      "--------------------------------------------------------------------------\n",
+      "--------------------------------------------------------------------------\n",
+      "mpirun detected that one or more processes exited with non-zero status, thus causing\n",
+      "the job to be terminated. The first process to do so was:\n",
+      "\n",
+      "  Process name: [[57467,1],5]\n",
+      "----------------------------------------------------------------------\n",
+      "\n",
+      "loading system configuration\n",
+      "loading user configuration from /Users/amaliedavidsen/.mcstas/3.2/mccode_config.json\n",
+      "INFO: Using directory: \"/Volumes/T7/KVASIR/KVASIR_for_repo/KVASIR\"\n",
+      "INFO: Regenerating c-file: KVASIR.c\n",
+      "CFLAGS= -Wl,-rpath,CMD(mcpl-config --show libdir) -LCMD(mcpl-config --show libdir) -lmcpl -ICMD(mcpl-config --show includedir)\n",
+      "      \n",
+      "-----------------------------------------------------------\n",
+      "\n",
+      "Generating single GPU kernel or single CPU section layout: \n",
+      "\n",
+      "-----------------------------------------------------------\n",
+      "\n",
+      "Generating GPU/CPU -DFUNNEL layout:\n",
+      "\n",
+      "-----------------------------------------------------------\n",
+      "INFO: Recompiling: ./KVASIR.out\n",
+      "./KVASIR.c:3393:1: warning: non-void function does not return a value in all control paths [-Wreturn-type]\n",
+      " 3393 | } /* siminfo_init */\n",
+      "      | ^\n",
+      "./KVASIR.c:14678:37: warning: format specifies type 'unsigned long' but the argument has type 'unsigned long long' [-Wformat]\n",
+      " 14674 |     MPI_MASTER(\n",
+      "       |     ~~~~~~~~~~~\n",
+      " 14675 | #endif\n",
+      "       | ~~~~~~\n",
+      " 14676 |            fprintf(stdout, \"\\n\\n Warning: You are using MCPL_input with a repeat_count of %lu:\\n - Minimum neutron count requested is %lu x %lu <= %lu\",\n",
+      "       |            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      "       |                                                                                                                                                    %llu\n",
+      " 14677 |                (long unsigned)repeat_count,(long unsigned)nparticles,\n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14678 |                (long unsigned)repeat_count,(long unsigned)repeat_cnt*nparticles); \n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14679 | #if defined (USE_MPI)\n",
+      "       | ~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14680 |   fprintf(stdout, \" x %i MPI nodes = %lu neutrons total\\n\",\n",
+      "       |   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14681 |     mpi_node_count,(long unsigned)mpi_node_count*repeat_cnt*nparticles);\n",
+      "       |     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14682 |      );\n",
+      "       |      ~\n",
+      "./KVASIR.c:456:5: note: expanded from macro 'MPI_MASTER'\n",
+      "  456 |   { statement; } \\\n",
+      "      |     ^~~~~~~~~\n",
+      "./KVASIR.c:14681:20: warning: format specifies type 'unsigned long' but the argument has type 'unsigned long long' [-Wformat]\n",
+      " 14674 |     MPI_MASTER(\n",
+      "       |     ~~~~~~~~~~~\n",
+      " 14675 | #endif\n",
+      "       | ~~~~~~\n",
+      " 14676 |            fprintf(stdout, \"\\n\\n Warning: You are using MCPL_input with a repeat_count of %lu:\\n - Minimum neutron count requested is %lu x %lu <= %lu\",\n",
+      "       |            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14677 |                (long unsigned)repeat_count,(long unsigned)nparticles,\n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14678 |                (long unsigned)repeat_count,(long unsigned)repeat_cnt*nparticles); \n",
+      "       |                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14679 | #if defined (USE_MPI)\n",
+      "       | ~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14680 |   fprintf(stdout, \" x %i MPI nodes = %lu neutrons total\\n\",\n",
+      "       |   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      "       |                                      %llu\n",
+      " 14681 |     mpi_node_count,(long unsigned)mpi_node_count*repeat_cnt*nparticles);\n",
+      "       |     ~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      " 14682 |      );\n",
+      "       |      ~\n",
+      "./KVASIR.c:456:5: note: expanded from macro 'MPI_MASTER'\n",
+      "  456 |   { statement; } \\\n",
+      "      |     ^~~~~~~~~\n",
+      "./KVASIR.c:14773:36: warning: format specifies type 'unsigned long' but the argument has type 'unsigned long long' [-Wformat]\n",
+      " 14771 |       fprintf(stdout, \"\\n\\n Warning: You are using MCPL_input with a repeat_count of %lu:\\n - Neutron count requested is %lu x %lu <= %lu\",\n",
+      "       |                                                                                                                                       ~~~\n",
+      "       |                                                                                                                                       %llu\n",
+      " 14772 |               (long unsigned)repeat_count,(long unsigned)read_neutrons,\n",
+      " 14773 |               (long unsigned)repeat_count,(long unsigned)repeat_cnt*read_neutrons);\n",
+      "       |                                           ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
+      "4 warnings generated.\n",
+      "[cctools-port]: generating fake signature for './KVASIR.out'\n",
+      "INFO: ===\n",
+      "Simulation 'KVASIR' (KVASIR.instr): running on 6 nodes (master is 'Amalies-MacBook-Air.local', MPI version 3.1).\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "MCPL \u001b[91mERROR\u001b[0m: Unable to open file!\n",
+      "--------------------------------------------------------------------------\n",
+      "Primary job  terminated normally, but 1 process returned\n",
+      "a non-zero exit code. Per user-direction, the job has been aborted.\n",
+      "--------------------------------------------------------------------------\n",
+      "--------------------------------------------------------------------------\n",
+      "mpirun detected that one or more processes exited with non-zero status, thus causing\n",
+      "the job to be terminated. The first process to do so was:\n",
+      "\n",
+      "  Process name: [[57467,1],5]\n",
+      "  Exit code:    1\n",
+      "--------------------------------------------------------------------------\n",
+      "\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "KVASIR.set_parameters(t=50)\n",
+    "KVASIR.settings(ncount=1e6, mpi=6)\n",
+    "data = KVASIR.backengine()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "No data to plot\n"
+     ]
+    }
+   ],
+   "source": [
+    "ms.make_sub_plot(data, log=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#KVASIR.show_instrument(format=\"window\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "mcstas-dev",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/mcstas-comps/examples/ESS/KVASIR/Virtual_output_example.mcpl b/mcstas-comps/examples/ESS/KVASIR/Virtual_output_example.mcpl
new file mode 100755
index 000000000..dca479c18
Binary files /dev/null and b/mcstas-comps/examples/ESS/KVASIR/Virtual_output_example.mcpl differ