diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index 1d28b0a593efa00f77ab2ee27627d08aa84eb4e9..25b6880f7e69e83fb1101863d91bd95f85df8c58 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -1599,8 +1599,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
       }
     } 
   } // i168 loop
-  last_configuration++;
-  oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
+  oindex = (jindex - 1) * (num_configs + 1) + num_configs;
   output->vec_sphere_sizes[oindex] = sze;
   output->vec_sphere_ref_indices[oindex] = entn;
   // label 160
@@ -2081,6 +2080,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 #ifdef USE_NVTX
   nvtxRangePop();
 #endif
+  if (jer == 0) output->vec_jxi[jindex - 1] = jxi488;
   interval_end = chrono::high_resolution_clock::now();
   elapsed = interval_end - interval_start;
   message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index e121e95fbd78493cc5ef34122921f047b1b8b634..4d6e61a2fb6273f9e665e038098d6f76ad0a868b 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -4058,7 +4058,7 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
   int result = 0;
   FILE *p_outfile = fopen(output.c_str(), "w");
   if (p_outfile != NULL) {
-    if (vec_xi[0] == 1) {
+    if (vec_jxi[0] == 1) {
       fprintf(p_outfile, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
 #ifdef USE_ILP64
       fprintf(
@@ -4102,7 +4102,10 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
       if (idfc < 0) {
 	fprintf( p_outfile, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n \n", vec_vk[0]);
       }
-      for (int jxi = 0; jxi < xi_block_size; jxi++) {
+      // End preamble writing
+    }
+    // Wavelength loop
+    for (int jxi = 0; jxi < xi_block_size; jxi++) {
 	int done_dirs = 0;
 	double alamb = 2.0 * 3.141592653589793238 / vec_vk[jxi];
 	fprintf(p_outfile, "========== JXI =%3d ====================\n", vec_jxi[jxi]);
@@ -4119,20 +4122,22 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
 	  break;
 	}
 	for (int i168 = 1; i168 <= configurations; i168++) {
-	  if (vec_sphere_ref_indices[i168 - 1] == cc0) {
-	    fprintf(p_outfile, "  SPHERE N.%2d: SIZE=%15.7lE\n", i168, vec_sphere_sizes[i168 - 1]);
+	  int cindex = jxi * (configurations + 1) + i168 - 1;
+	  if (vec_sphere_ref_indices[cindex] == cc0) {
+	    fprintf(p_outfile, "  SPHERE N.%2d: SIZE=%15.7lE\n", i168, vec_sphere_sizes[cindex]);
 	  } else {
 	    fprintf(
 	      p_outfile, "  SPHERE N.%2d: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
-	      i168, vec_sphere_sizes[i168 - 1], real(vec_sphere_ref_indices[i168 - 1]),
-	      imag(vec_sphere_ref_indices[i168 - 1])
+	      i168, vec_sphere_sizes[cindex], real(vec_sphere_ref_indices[cindex]),
+	      imag(vec_sphere_ref_indices[cindex])
 	    );
 	  }
 	} // i168 configuration loop
+	int cindex = jxi * (configurations + 1) + configurations;
 	fprintf(
 	  p_outfile, "  EXT. SPHERE: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
-	  vec_sphere_sizes[configurations], real(vec_sphere_ref_indices[configurations]),
-	  imag(vec_sphere_ref_indices[configurations])
+	  vec_sphere_sizes[cindex], real(vec_sphere_ref_indices[cindex]),
+	  imag(vec_sphere_ref_indices[cindex])
 	);
 	fprintf(p_outfile, "     ENSEMBLE AVERAGE, MODE%2d\n", iavm);
 	if (inpol == 0) fprintf(p_outfile, "   LIN -1\n");
@@ -4465,8 +4470,7 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
 	    } // jths loop
 	  } // jph loop
 	} // jth loop
-      } // jxi wavelength loop
-    }
+    } // jxi wavelength loop
     fclose(p_outfile);
   } else {
     result = -1;
diff --git a/src/testing/test_inclusion_outputs.cpp b/src/testing/test_inclusion_outputs.cpp
index e546e8989069dafe5ce35dc4174dec5fcdf83fc1..ac93471acc09b4b03d09e294b08a69b6cab02884 100644
--- a/src/testing/test_inclusion_outputs.cpp
+++ b/src/testing/test_inclusion_outputs.cpp
@@ -27,23 +27,21 @@ int test_inclusion_devel();
 
 int main() {
   int result = 0;
-  // result += test_inclusion_hdf5_output(); // 1 if failed
+  result += test_inclusion_hdf5_output(); // 1 if failed
   result += test_inclusion_devel(); // 10 if failed
   return result;
 }
 
 int test_inclusion_hdf5_output() {
   int result = 0;
-  /*
   try {
-    const string hdf5_file = "../../test_data/cluster/c_OCLU_24.hd5";
-    ClusterOutputInfo *oi = new ClusterOutputInfo(hdf5_file);
-    oi->write("c_OCLU_24", "LEGACY");
+    const string hdf5_file = "../../test_data/inclusion/c_OINCLU.hd5";
+    InclusionOutputInfo *oi = new InclusionOutputInfo(hdf5_file);
+    oi->write("c_OINCLU", "LEGACY");
     delete oi;
   } catch (const exception& ex) {
     result = 1;
   }
-  */
   return result;
 }
 
diff --git a/test_data/inclusion/c_OINCLU.hd5 b/test_data/inclusion/c_OINCLU.hd5
index b136f7b2ec732e52174c4a69a88e45c7876394cf..c32b085708391b2755357c05f1fb32a63f5772a0 100644
Binary files a/test_data/inclusion/c_OINCLU.hd5 and b/test_data/inclusion/c_OINCLU.hd5 differ