diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 033e111cad90872665c1ccb4fb7e191f000bdd37..38c452348df0eb94c241fe797880f17ac18a23c5 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -548,6 +548,7 @@ int sphere_jxi488_cycle(
   Logger *logger = new Logger(LOG_INFO);
   int oindex = 0;
   int jxi = jxi488 + 1;
+  int jxindex = jxi - oi->first_xi;
   int nsph = gconf->number_of_spheres;
   int l_max = gconf->l_max;
   int in_pol = gconf->in_pol;
@@ -570,14 +571,14 @@ int sphere_jxi488_cycle(
   if (idfc >= 0) {
     vk = xi * wn;
     vkarg = vk;
-    oi->vec_vk[jxi488] = vk;
-    oi->vec_xi[jxi488] = xi;
+    oi->vec_vk[jxindex] = vk;
+    oi->vec_xi[jxindex] = xi;
   } else { // IDFC < 0
     vk = sconf->xip * wn;
     vkarg = xi * vk;
     sqsfi = 1.0 / (xi * xi);
-    oi->vec_vk[jxi488] = vk;
-    oi->vec_xi[jxi488] = xi;
+    oi->vec_vk[jxindex] = vk;
+    oi->vec_xi[jxindex] = xi;
   }
   vtppoanp->append_line(VirtualBinaryLine(vk));
   double thsca = (gconf->isam > 1) ? sa->ths - sa->th : 0.0;
@@ -607,7 +608,7 @@ int sphere_jxi488_cycle(
     if (nsh % 2 == 0) sid->c1->dc0[ici] = exdc;
     dme(l_max, i, npnt, npntts, vkarg, exdc, exri, sid->c1, jer, lcalc, sid->arg);
     if (jer != 0) {
-      oi->vec_ier[jxi] = 1;
+      oi->vec_ier[jxindex] = 1;
       oi->lcalc = lcalc;
       return jer;
     }
@@ -635,7 +636,7 @@ int sphere_jxi488_cycle(
     int i = i170 + 1;
     if (sid->c1->iog[i170] >= i) {
       last_configuration++;
-      oindex = jxi488 * configurations + last_configuration - 1;
+      oindex = jxindex * configurations + last_configuration - 1;
       double albeds = sid->c1->sscs[i170] / sid->c1->sexs[i170];
       sid->c1->sqscs[i170] *= sqsfi;
       sid->c1->sqabs[i170] *= sqsfi;
@@ -698,15 +699,15 @@ int sphere_jxi488_cycle(
     } // End if iog[i170] >= i
   } // i170 loop
   if (nsph != 1) {
-    oi->vec_fsat[jxi488] = sid->tfsas;
+    oi->vec_fsat[jxindex] = sid->tfsas;
     double csch = 2.0 * vk * sqsfi / sid->c1->gcs;
     dcomplex s0 = sid->tfsas * exri;
     double qschu = csch * imag(s0);
     double pschu = csch * real(s0);
     double s0mag = cs0 * cabs(s0);
-    oi->vec_qschut[jxi488] = qschu;
-    oi->vec_pschut[jxi488] = pschu;
-    oi->vec_s0magt[jxi488] = s0mag;
+    oi->vec_qschut[jxindex] = qschu;
+    oi->vec_pschut[jxindex] = pschu;
+    oi->vec_s0magt[jxindex] = s0mag;
   }
   double th = sa->th;
   int done_dirs = 0;
@@ -721,12 +722,11 @@ int sphere_jxi488_cycle(
     for (int jph484 = 0; jph484 < sa->nph; jph484++) {
       int jph = jph484 + 1;
       bool goto182 = (sa->nk == 1) && (jxi > 1);
-      double cost, sint, cosp, sinp;
       if (!goto182) {
-	upvmp(th, ph, 0, cost, sint, cosp, sinp, sid->u, sid->upmp, sid->unmp);
+	upvmp(th, ph, 0, sid->cost, sid->sint, sid->cosp, sid->sinp, sid->u, sid->upmp, sid->unmp);
       }
       if (gconf->isam >= 0) {
-	wmamp(0, cost, sint, cosp, sinp, in_pol, l_max, 0, nsph, sid->argi, sid->u, sid->upmp, sid->unmp, sid->c1);
+	wmamp(0, sid->cost, sid->sint, sid->cosp, sid->sinp, in_pol, l_max, 0, nsph, sid->argi, sid->u, sid->upmp, sid->unmp, sid->c1);
 	for (int i183 = 0; i183 < nsph; i183++) {
 	  double rapr = sid->c1->sexs[i183] - sid->gaps[i183];
 	  frx = rapr * sid->u[0];
@@ -752,16 +752,15 @@ int sphere_jxi488_cycle(
 	double phs = sa->phs;
 	for (int jphs480 = 0; jphs480 < sa->nphs; jphs480++) {
 	  int jphs = jphs480 + 1;
-	  double costs, sints, cosps, sinps;
 	  if (gconf->isam >= 1) {
 	    phs = ph + phsph;
 	    if (phs >= 360.0) phs -= 360.0;
 	  }
 	  bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 1));
 	  if (!goto190) {
-	    upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, sid->us, sid->upsmp, sid->unsmp);
+	    upvmp(ths, phs, icspnv, sid->costs, sid->sints, sid->cosps, sid->sinps, sid->us, sid->upsmp, sid->unsmp);
 	    if (gconf->isam >= 0) {
-	      wmamp(2, costs, sints, cosps, sinps, in_pol, l_max, 0, nsph, sid->args, sid->us, sid->upsmp, sid->unsmp, sid->c1);
+	      wmamp(2, sid->costs, sid->sints, sid->cosps, sid->sinps, in_pol, l_max, 0, nsph, sid->args, sid->us, sid->upsmp, sid->unsmp, sid->c1);
 	    }
 	  }
 	  if (nkks != 0 || jxi == 1) {
@@ -772,10 +771,10 @@ int sphere_jxi488_cycle(
 	    );
 	    if (gconf->isam < 0) {
 	      wmasp(
-		cost, sint, cosp, sinp, costs, sints, cosps, sinps,
-		sid->u, sid->up, sid->un, sid->us, sid->ups, sid->uns,
-		isq, ibf, in_pol, l_max, 0, nsph, sid->argi, sid->args,
-		sid->c1
+		sid->cost, sid->sint, sid->cosp, sid->sinp, sid->costs,
+		sid->sints, sid->cosps, sid->sinps, sid->u, sid->up,
+		sid->un, sid->us, sid->ups, sid->uns, isq, ibf, in_pol,
+		l_max, 0, nsph, sid->argi, sid->args, sid->c1
 	      );
 	    }
 	    for (int i193 = 0; i193 < 3; i193++) {
@@ -823,24 +822,24 @@ int sphere_jxi488_cycle(
 	  last_configuration = 0;
 	  for (int ns226 = 0; ns226 < nsph; ns226++) {
 	    int ns = ns226 + 1;
-	    oindex = jxi488 * nsph * ndirs + nsph * done_dirs + ns226;
+	    oindex = jxindex * nsph * ndirs + nsph * done_dirs + ns226;
 	    oi->vec_dir_sas11[oindex] = sid->c1->sas[ns226][0][0];
 	    oi->vec_dir_sas21[oindex] = sid->c1->sas[ns226][1][0];
 	    oi->vec_dir_sas12[oindex] = sid->c1->sas[ns226][0][1];
 	    oi->vec_dir_sas22[oindex] = sid->c1->sas[ns226][1][1];
-	    oi->vec_dir_fx[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frx;
-	    oi->vec_dir_fy[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = fry;
-	    oi->vec_dir_fz[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frz;
+	    oi->vec_dir_fx[jxindex * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frx;
+	    oi->vec_dir_fy[jxindex * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = fry;
+	    oi->vec_dir_fz[jxindex * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frz;
 	    for (int i225 = 0; i225 < 16; i225++) sid->c1->vint[i225] = sid->c1->vints[ns226][i225];
 	    mmulc(sid->c1->vint, sid->cmullr, sid->cmul);
 	    for (int imul = 0; imul < 4; imul++) {
-	      int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;
+	      int muls_index = 16 * jxindex * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;
 	      for (int jmul = 0; jmul < 4; jmul++) {
 		oi->vec_dir_muls[muls_index + jmul] = sid->cmul[imul][jmul];
 	      }
 	    }
 	    for (int imul = 0; imul < 4; imul++) {
-	      int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;
+	      int muls_index = 16 * jxindex * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;
 	      for (int jmul = 0; jmul < 4; jmul++) {
 		oi->vec_dir_mulslr[muls_index + jmul] = sid->cmullr[imul][jmul];
 	      }
@@ -867,7 +866,7 @@ int sphere_jxi488_cycle(
     } // jph484 loop on elevation
     th += sa->thstp;
   } // jth486 loop on azimuth
-  oi->vec_jxi[jxi488] = jxi;
+  oi->vec_jxi[jxindex] = jxi;
   logger->log("INFO: finished scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");
   return jer;
 }
@@ -886,6 +885,14 @@ SphereIterationData::SphereIterationData(
   c1 = new ParticleDescriptorSphere(gconf, sconf);
   argi = new double[1]();
   args = new double[1]();
+  cost = 0.0;
+  sint = 0.0;
+  cosp = 0.0;
+  sinp = 0.0;
+  costs = 0.0;
+  sints = 0.0;
+  cosps = 0.0;
+  sinps = 0.0;
   scan = 0.0;
   cfmp = 0.0;
   cfsp = 0.0;
@@ -959,6 +966,14 @@ SphereIterationData::SphereIterationData(const SphereIterationData &rhs) {
   args = new double[1];
   argi[0] = rhs.argi[0];
   args[0] = rhs.args[0];
+  cost = rhs.cost;
+  sint = rhs.sint;
+  cosp = rhs.cosp;
+  sinp = rhs.sinp;
+  costs = rhs.costs;
+  sints = rhs.sints;
+  cosps = rhs.cosps;
+  sinps = rhs.sinps;
   scan = rhs.scan;
   cfmp = rhs.cfmp;
   cfsp = rhs.cfsp;
@@ -1095,6 +1110,14 @@ SphereIterationData::SphereIterationData(const mixMPI *mpidata, const int device
   MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(&s0, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&cost, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sint, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&cosp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sinp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&costs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sints, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&cosps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sinps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
   MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
   MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
   MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -1173,6 +1196,14 @@ int SphereIterationData::mpibcast(const mixMPI *mpidata) {
   MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(&s0, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&cost, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sint, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&cosp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sinp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&costs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sints, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&cosps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&sinps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
   MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
   MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
   MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);