Update Phonon_simple fixing and modernising

Make sure that vf_list has length of 100 and allow the user to choose the amount of scan steps actually taken, if they need to optimise the instrument. This counteracts overflow of vf_list.

Also change from parms array to parms struct in order to improve readability.

Author: Lomholy

mcstas-comps/samples/Phonon_simple.comp

@@ -66,7 +66,9 @@
 DEFINE COMPONENT Phonon_simple
 
 SETTING PARAMETERS (radius,yheight,sigma_abs,sigma_inc,a,b,M,c,DW,T,
-target_x=0, target_y=0, target_z=0, int target_index=0,focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0, gap=0)
+              target_x=0, target_y=0, target_z=0, int target_index=0,
+              focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0, 
+              gap=0, int e_steps_low=50, int e_steps_high=50)
 
 /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
 SHARE
@@ -75,6 +77,24 @@ SHARE
 #define PHONON_SIMPLE $Revision$
 #define T2E (1/11.605)   /* Kelvin to meV */
 
+	struct neutron_and_phonon_params {
+		// Statically allocate vectors that are always 3 
+    double vf; // Final velocity size
+    double vi; // Initial velocity size
+    double vv_x; // vv is the unit vector of the final velocity vector
+    double vv_y;
+    double vv_z;
+    double vi_x; // vi is the initial velocity vector
+    double vi_y;
+    double vi_z;
+    double a_; // d spacing of the cubic lattice
+    double c_; // Speed of sound in the material
+    double gap_; // optional spin gap
+    double ah; // Half of a 
+    int e_steps_high_;
+    int e_steps_low_;
+	};
+
 #pragma acc routine 
 double nbose(double omega, double T)  /* Other name ?? */
   {
@@ -103,26 +123,26 @@ void fatalerror(char *s)
   }
 
   #pragma acc routine
-  double omega_q(double* parms)
+  double omega_q(struct neutron_and_phonon_params *parms)
     {
       /* dispersion in units of meV  */
       double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z;
       double q, qx, qy, qz, Jq, res_phonon, res_neutron;
       double ah, a, c;
       double gap;
 
-      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];
-      a   =parms[8];
-      c   =parms[9];
-      gap =parms[10];
-      ah=a/2.0;
+      vf=parms->vf;
+      vi=parms->vi;
+      vv_x=parms->vv_x;
+      vv_y=parms->vv_y;
+      vv_z=parms->vv_z;
+      vi_x=parms->vi_x;
+      vi_y=parms->vi_y;
+      vi_z=parms->vi_z;
+      a   =parms->a_;
+      c   =parms->c_;
+      gap =parms->gap_;
+      ah  =parms->ah;
 
       qx=V2K*(vi_x-vf*vv_x);
       qy=V2K*(vi_y-vf*vv_y);
@@ -141,15 +161,15 @@ void fatalerror(char *s)
     }
 
 
-double zridd(double (*func)(double*), double x1, double x2, double *parms, double xacc)
+double zridd(double (*func)(struct neutron_and_phonon_params*), double x1, double x2, struct neutron_and_phonon_params *parms, double xacc)
     {
       int j;
       double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew;
 
-      parms[0]=x1;
-      fl=(*func)(parms);
-      parms[0]=x2;
-      fh=(*func)(parms);
+      parms->vf=x1;
+      fl=func(parms);
+      parms->vf=x2;
+      fh=func(parms);
       if (fl*fh >= 0)
       {
         if (fl==0) return x1;
@@ -164,17 +184,17 @@ double zridd(double (*func)(double*), double x1, double x2, double *parms, doubl
         for (j=1; j<MAXRIDD; j++)
         {
           xm=0.5*(xl+xh);
-          parms[0]=xm;
-          fm=(*func)(parms);
+          parms->vf=xm;
+          fm=func(parms);
           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;
-          fnew=(*func)(parms);
+          parms->vf=ans;
+          fnew=func(parms);
           if (fnew == 0.0) return ans;
           if (fabs(fm)*SIGN(fnew) != fm)
           {
@@ -206,14 +226,14 @@ double zridd(double (*func)(double*), double x1, double x2, double *parms, doubl
     }
 
 #pragma acc routine 
-double zridd_gpu(double x1, double x2, double *parms, double xacc)
+double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, double xacc)
     {
       int j;
       double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew;
 
-      parms[0]=x1;
+      parms->vf=x1;
       fl=omega_q(parms);
-      parms[0]=x2;
+      parms->vf=x2;
       fh=omega_q(parms);
       if (fl*fh >= 0)
       {
@@ -229,7 +249,7 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc)
         for (j=1; j<MAXRIDD; j++)
         {
           xm=0.5*(xl+xh);
-          parms[0]=xm;
+          parms->vf=xm;
           fm=omega_q(parms);
           s=sqrt(fm*fm-fl*fh);
           if (s == 0.0)
@@ -238,7 +258,7 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc)
           if (fabs(xnew-ans) <= xacc)
             return ans;
           ans=xnew;
-          parms[0]=ans;
+          parms->vf=ans;
           fnew=omega_q(parms);
           if (fnew == 0.0) return ans;
           if (fabs(fm)*SIGN(fnew) != fm)
@@ -273,30 +293,40 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc)
 
 #define ROOTACC 1e-8
 
-  int findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(double*), double *parms)
+  void findroots(double brack_low, double brack_mid, double brack_high, 
+                double *list, int* index, double (*f)(struct neutron_and_phonon_params*), 
+                struct neutron_and_phonon_params *parms)
     {
-      double root,range=brack_mid-brack_low;
-      int i, steps=100;
-
-     for (i=0; i<steps; i++)
-     {
-      root = zridd(f, brack_low+range*i/(int)steps,
-                   brack_low+range*(i+1)/(int)steps,
-                   (double *)parms, ROOTACC);
+      double root;
+      // Energy gain and energy loss spaces are not equally big. We check uniformly
+      // So we use two different ranges
+      double range_low = brack_mid - brack_low; 
+      double range_high = brack_high - brack_mid;
+      // First in energy loss for the neutron
+     for (int i=0; i<parms->e_steps_low_; i++){
+      root = zridd(f, brack_low+range_low*i/ parms->e_steps_low_,
+                   brack_low+range_low*(i+1)/ parms->e_steps_low_,
+                   parms, ROOTACC);
       if (root != UNUSED)
       {
         list[(*index)++]=root;
       }
      }
-      root = zridd(f, brack_mid, brack_high, (double *)parms, ROOTACC);
-      if (root != UNUSED)
-      {
+
+    // Then in energy gain for the neutron
+    for (int i=0; i<parms->e_steps_high_; i++){
+      root = zridd(f, brack_low+range_high*i/ parms->e_steps_high_,
+                   brack_low+range_high*(i+1)/ parms->e_steps_high_,
+                   parms, ROOTACC);
+      if (root != UNUSED){
         list[(*index)++]=root;
       }
+     }
     }
 
 #pragma acc routine 
-  int findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, double *parms)
+  void findroots_gpu(double brack_low, double brack_mid, double brack_high, 
+                    double *list, int* index, struct neutron_and_phonon_params *parms)
     {
       double root,range=brack_mid-brack_low;
       int i, steps=100;
@@ -305,13 +335,13 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc)
      {
        root = zridd_gpu(brack_low+range*i/(int)steps,
                    brack_low+range*(i+1)/(int)steps,
-                   (double *)parms, ROOTACC);
+                   parms, ROOTACC);
       if (root != UNUSED)
       {
         list[(*index)++]=root;
       }
      }
-      root = zridd_gpu(brack_mid, brack_high, (double *)parms, ROOTACC);
+      root = zridd_gpu(brack_mid, brack_high, parms, ROOTACC);
       if (root != UNUSED)
       {
         list[(*index)++]=root;
@@ -329,13 +359,25 @@ DECLARE
   double V_my_s;
   double V_my_a_v;
   double DV;
+  double vf_list[100];// List of allowed final velocities. Has length of scan_steps
+  struct neutron_and_phonon_params parms;
+  
 %}
 INITIALIZE
 %{
+  memset(vf_list, 0, sizeof(vf_list));
   V_rho = 4/(a*a*a);
   V_my_s = (V_rho * 100 * sigma_inc);
   V_my_a_v = (V_rho * 100 * sigma_abs * 2200);
   DV = 0.001;   /* Velocity change used for numerical derivative */
+  
+  // Set constant parameters for parms object
+  parms.a_ = a;
+  parms.c_ = c;
+  parms.gap_ = gap;
+  parms.ah = a/2.0;
+  parms.e_steps_high_ = e_steps_high;
+  parms.e_steps_low_ = e_steps_low;
 
   /* now compute target coords if a component index is supplied */
   if (!target_index && !target_x && !target_y && !target_z) target_index=1;
@@ -368,11 +410,9 @@ TRACE
   double bose_factor;        /* Calculated value of the Bose factor */
   double omega;              /* energy transfer */
   int nf, index;                   /* Number of allowed final velocities */
-  double vf_list[2];             /* List of allowed final velocities */
   double J_factor;            /* Jacobian from delta fnc.s in cross section */
   double f1, f2;            /* probed values of omega_q minus omega */
   double p1,p2,p3,p4,p5;    /* temporary multipliers */
-  double parms[11];
 
   if(cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight))
   {
@@ -405,29 +445,27 @@ TRACE
     }
     NORM(vx, vy, vz);
     nf=0;
-      parms[0]=-1;
-      parms[1]=v_i;
-      parms[2]=vx;
-      parms[3]=vy;
-      parms[4]=vz;
-      parms[5]=vx_i;
-      parms[6]=vy_i;
-      parms[7]=vz_i;
-      parms[8]=a;
-      parms[9]=c;
-      parms[10]=gap;
+    parms.vf = -1;
+    parms.vi = v_i;
+    parms.vv_x = vx;
+    parms.vv_y = vy;
+    parms.vv_z = vz;
+    parms.vi_x = vx_i;
+    parms.vi_y = vy_i;
+    parms.vi_z = vz_i;
+
     #ifndef OPENACC
-    findroots(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, omega_q, parms);
+    findroots(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, omega_q, &parms);
     #else
-    findroots_gpu(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, parms);
+    findroots_gpu(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, &parms);
     #endif
     index=(int)floor(rand01()*nf);
     if (index<2) {
       v_f=vf_list[index];
-      parms[0]=v_f-DV;
-      f1=omega_q(parms);
-      parms[0]=v_f+DV;
-      f2=omega_q(parms);
+      parms.vf=v_f-DV;
+      f1=omega_q(&parms);
+      parms.vf=v_f+DV;
+      f2=omega_q(&parms);
       J_factor = fabs(f2-f1)/(2*DV);
       omega=VS2E*(v_i*v_i-v_f*v_f);
       vx *= v_f;
@@ -460,6 +498,9 @@ TRACE
       ABSORB; // findroots returned junk
     }
   } /* else transmit: Neutron did not hit the sample */
+
+  // Finally, set vf_list to zero again
+  memset(vf_list, 0, nf * sizeof(double));
 %}
 
 MCDISPLAY