diff --git a/docs/sphinx/source/output/diagnostics.rst b/docs/sphinx/source/output/diagnostics.rst
index 48874b679..485f0debb 100644
--- a/docs/sphinx/source/output/diagnostics.rst
+++ b/docs/sphinx/source/output/diagnostics.rst
@@ -10,3 +10,8 @@ where the in a multiprocessor run, the number refers to the log of a specific th
 
 Inspecting these logs is important for understanding whether a Python run is successful,
 and for understanding why if failed if that is the case. 
+
+Additionally in this directory one will find a file **whatever.disk.diag** if the system contained a disk.  This file contains information
+about the structure of the disk as a function of disk radius. It also indicates how many photons hit the disk with how much total energy
+as a function of radius.  It provides several estimates of the effective temperature of the radiation that hits the disk, and determines 
+how the disk temperature of the disk would so that it would reradiate the sum of the viscous energy and the energy associated with illumination.
diff --git a/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst b/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst
index 5b3744935..34c24f0be 100644
--- a/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst
+++ b/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst
@@ -14,8 +14,11 @@
    .. autosummary::
    
       add_gtot
+      combine_collisions
       plot_xsec
+      read_collisions
       read_phot
+      redo_collisions
       redo_lines
       redo_phot
       reweight
diff --git a/source/Makefile b/source/Makefile
index e7b246956..33efedaf7 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -147,7 +147,7 @@ LDFLAGS+= -L$(LIB) -lm -lgsl -lgslcblas $(CUDA_LIBS)
 #Note that version should be a single string without spaces.
 
 
-VERSION = 88d
+VERSION = 88f
 
 
 
@@ -264,6 +264,7 @@ rhf_objects = rad_hydro_files.o $(python_objects)
 modify_wind_objects = modify_wind.o $(python_objects)
 inspect_wind_objects = inspect_wind.o $(python_objects)
 py_optd_obj = py_optd.o py_optd_output.o py_optd_trans.o py_optd_util.o py_optd_extern_init.o $(python_objects)
+windsave2fits_objects = windsave2fits.o $(python_objects)
 
 run_indent:
 	../py_progs/run_indent.py -all_no_headers
@@ -280,6 +281,12 @@ windsave2table: startup $(table_objects) $(CUDA_OBJECTS)
 	mv $@ $(BIN)/windsave2table$(VERSION)
 	@if [ $(INDENT) = yes ] ; then  ../py_progs/run_indent.py -changed ; fi
 
+windsave2fits: startup $(windsave2fits_objects) $(CUDA_OBJECTS)
+	$(CC) $(CFLAGS) $(windsave2fits_objects) $(CUDA_OBJECTS) $(LDFLAGS) -lcfitsio -o windsave2fits 
+	cp $@ $(BIN)
+	mv $@ $(BIN)/windsave2fits$(VERSION)
+	@if [ $(INDENT) = yes ] ; then  ../py_progs/run_indent.py -changed ; fi
+
 rad_hydro_files: startup $(rhf_objects) $(CUDA_OBJECTS)
 	$(CC) $(CFLAGS) $(rhf_objects) $(CUDA_OBJECTS) $(LDFLAGS) -o rad_hydro_files
 	cp $@ $(BIN)
diff --git a/source/communicate_plasma.c b/source/communicate_plasma.c
index 768e46a56..1e84baa98 100644
--- a/source/communicate_plasma.c
+++ b/source/communicate_plasma.c
@@ -986,11 +986,11 @@ reduce_simple_estimators (void)
   ion_helper2 = calloc (sizeof (double), NPLASMA * nions);
   inner_ion_helper = calloc (sizeof (double), NPLASMA * n_inner_tot);
   inner_ion_helper2 = calloc (sizeof (double), NPLASMA * n_inner_tot);
-  /* JM -- added routine to average the qdisk quantities. The 2 is because
-     we only have two doubles to worry about (heat and ave_freq) and
+  /* Routine to average the qdisk quantities. The 3 is because
+     we have three doubles to worry about (emit, heat and ave_freq) and
      two integers (nhit and nphot) */
-  qdisk_helper = calloc (sizeof (double), NRINGS * 2);
-  qdisk_helper2 = calloc (sizeof (double), NRINGS * 2);
+  qdisk_helper = calloc (sizeof (double), NRINGS * 3);
+  qdisk_helper2 = calloc (sizeof (double), NRINGS * 3);
 
   flux_helper = calloc (sizeof (double), NPLASMA * NFLUX_ANGLES * 3);
   flux_helper2 = calloc (sizeof (double), NPLASMA * NFLUX_ANGLES * 3);
@@ -1068,6 +1068,7 @@ reduce_simple_estimators (void)
   {
     qdisk_helper[mpi_i] = qdisk.heat[mpi_i] / np_mpi_global;
     qdisk_helper[mpi_i + NRINGS] = qdisk.ave_freq[mpi_i] / np_mpi_global;
+    qdisk_helper[mpi_i + 2 * NRINGS] = qdisk.emit[mpi_i] / np_mpi_global;
   }
 
 
@@ -1080,7 +1081,7 @@ reduce_simple_estimators (void)
   MPI_Allreduce (flux_helper, flux_helper2, NPLASMA * 3 * NFLUX_ANGLES, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
   MPI_Allreduce (ion_helper, ion_helper2, NPLASMA * nions, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
   MPI_Allreduce (inner_ion_helper, inner_ion_helper2, NPLASMA * n_inner_tot, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-  MPI_Allreduce (qdisk_helper, qdisk_helper2, 2 * NRINGS, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce (qdisk_helper, qdisk_helper2, 3 * NRINGS, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
   /* Unpacking stuff */
   for (mpi_i = 0; mpi_i < NPLASMA; mpi_i++)
@@ -1157,6 +1158,7 @@ reduce_simple_estimators (void)
   {
     qdisk.heat[mpi_i] = qdisk_helper2[mpi_i];
     qdisk.ave_freq[mpi_i] = qdisk_helper2[mpi_i + NRINGS];
+    qdisk.emit[mpi_i] = qdisk_helper2[mpi_i + 2 * NRINGS];
   }
 
   /* now we've done all the doubles so we can free their helper arrays */
diff --git a/source/disk_init.c b/source/disk_init.c
index b87121e39..4b41a29ad 100644
--- a/source/disk_init.c
+++ b/source/disk_init.c
@@ -161,8 +161,8 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_extract, ftot)
     }
     else if (spectype == SPECTYPE_BB_FCOL)
     {
-      fcol = disk_colour_correction(t);
-      emit = emittance_bb (freqmin, freqmax, fcol * t) / pow(fcol, 4.0);
+      fcol = disk_colour_correction (t);
+      emit = emittance_bb (freqmin, freqmax, fcol * t) / pow (fcol, 4.0);
     }
     else
     {
@@ -220,8 +220,8 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_extract, ftot)
     }
     else if (spectype == SPECTYPE_BB_FCOL)
     {
-      fcol = disk_colour_correction(t);
-      emit = emittance_bb (freqmin, freqmax, fcol * t) / pow(fcol, 4.0);
+      fcol = disk_colour_correction (t);
+      emit = emittance_bb (freqmin, freqmax, fcol * t) / pow (fcol, 4.0);
     }
     else
     {
@@ -354,13 +354,72 @@ qdisk_init (rmin, rmax, m, mdot)
 }
 
 
+/**********************************************************/
+/**
+ * @brief      Reinitialize portions a structure (qdisk) for recording information about 
+ *      photons/energy impinging
+ * 	on the disk.
+ *
+ * @return     Always return zero
+ *
+ *
+ * ###Notes###
+ *
+ * The structure qdisk records information about how many photons hit the disk
+ * and what the impact of these photons would be if these
+ * photons are used to modify the temperatures structure
+ * of the disk
+ *
+ * This routine records information about the 
+ * photons which were emitted from the disk 
+ * and zeros various arrays that will record information
+ * about photons that will hit the disk
+ * 
+ *
+ **********************************************************/
+
+int
+qdisk_reinit (p)
+     PhotPtr p;
+{
+  int nphot, i, n;
+  double rho;
+  struct photon pp;
+  for (n = 0; n < NRINGS; n++)
+  {
+    qdisk.emit[n] = qdisk.nhit[n] = qdisk.heat[n] = qdisk.nphot[n] = qdisk.w[n] = qdisk.ave_freq[n] = 0;
+  }
+
+  for (nphot = 0; nphot < NPHOT; nphot++)
+  {
+    stuff_phot (&p[nphot], &pp);
+    if (pp.origin_orig == PTYPE_DISK)
+    {
+      rho = sqrt (pp.x[0] * pp.x[0] + pp.x[1] * pp.x[1]);
+
+      i = 0;
+      while (rho > qdisk.r[i] && i < NRINGS - 1)
+        i++;
+      i--;                      /* So that the heating refers to the heating between i and i+1 */
+
+      qdisk.emit[i] += pp.w;
+      qdisk.nphot[i] += 1;
+
+
+    }
+  }
+
+  return (0);
+}
+
+
 /**********************************************************/
 /**
  * @brief      Save the information in the qdisk structure to a file
  *
  * @param [in] char *  diskfile   Name of the file whihc is writteen
- * @param [in] double  ztot   The total luminosity of the disk as
- *                      calculate over multiple cycles
+ * @param [in] int  ichoice   A switch that just says whether the weighted
+ *       frequency has been converted to a real frequency (0 = No, 1=Yes)
  * @return     Always returns 0
  *
  * The routine reformats the data about disk heating which has
@@ -368,55 +427,91 @@ qdisk_init (rmin, rmax, m, mdot)
  *
  * ###Notes###
  *
- * The data concerning heating by the disk is built up during
- * a the ionization cycles
  *
  * The file that is produced should be readable as an astropy
  * table
  *
+ * Here is a definition of the columnts
+ *
+ * * r            the radius of the wirng
+ * * v            the velocity of the ring
+ * * zdisk        the zheight ot the rign
+ * * t_disk       the temperature derived from viscous heating or the input grid
+ * * g            the log of the gravity
+ * * emit         the luminosity of photons emitted by the ring in this cycle
+ * * nemit        the number photons emitted by the ring in this cycle
+ * * heat         the total luminosity of photons that hit the ring
+ * * nhit         the number of photon bundles hitting a ring
+ * * ehit/emit    the ratio of the energy of photons that hit the ring/to those
+ *              to those that are emitted by the ring
+ * * t_heat       the temperature calculated from the energy of the
+ *              photons hitting a ring
+ * * t_freq       The temperature calculated from the weighted frequency of 
+ *              photons hitting the ring
+ * * W_freq       Assuming t_freq for the temperature, this gives ther
+ *              ratio of the energy flux to LTE
+ * * t_rerad      The temperature required to radiate away both the
+ *              the viscous heating and the heating from photons
+ *              in the next cycle
+* 
  **********************************************************/
 
 int
-qdisk_save (diskfile, ztot)
+qdisk_save (diskfile, ichoice)
      char *diskfile;
-     double ztot;
+     int ichoice;
 {
   FILE *qptr;
   int n;
   double area, theat, ttot;
+  double ratio[NRINGS];
   qptr = fopen (diskfile, "w");
   fprintf (qptr,
-           "         r          v      zdisk    t_disk         g      heat  nhit nhit/emit    t_heat  t_irrad    W_irrad     t_tot\n");
-//  fprintf (qptr,
-//           "---------- ---------- ---------- --------- --------- --------- ----- --------- --------- --------- --------- ---------\n");
+           "         r          v      zdisk    t_disk         g       emit      nemit      heat       nhit ehit/emit    t_heat    t_freq     W_freq   t_rerad\n");
 
   for (n = 0; n < NRINGS; n++)
   {
     if (n < NRINGS - 1)
     {
+      // The factor of 2 comes from the fact that the disk has two sides
       area = (2. * PI * (qdisk.r[n + 1] * qdisk.r[n + 1] - qdisk.r[n] * qdisk.r[n]));
     }
     else
     {
       area = 0;
     }
-    theat = qdisk.heat[n] / area;
-    theat = pow (theat / STEFAN_BOLTZMANN, 0.25);
-    //theat is temperature if no internal energy production
     if (qdisk.nhit[n] > 0)
     {
 
-      qdisk.ave_freq[n] /= qdisk.heat[n];
+      /* During photon transport ave_freq is actually the w x frequency */
+      if (ichoice == 0)
+        qdisk.ave_freq[n] /= qdisk.heat[n];
+
       qdisk.t_hit[n] = PLANCK * qdisk.ave_freq[n] / (BOLTZMANN * 3.832);
       //Basic conversion from freq to T
       qdisk.w[n] = qdisk.heat[n] / (4. * PI * STEFAN_BOLTZMANN * area * qdisk.t_hit[n] * qdisk.t_hit[n] * qdisk.t_hit[n] * qdisk.t_hit[n]);
+      ratio[n] = 99.;
+      if (qdisk.emit[n] > 0.0)
+        ratio[n] = qdisk.heat[n] / qdisk.emit[n];
+      theat = qdisk.heat[n] / area;
+      theat = pow (theat / STEFAN_BOLTZMANN, 0.25);
+      //theat is temperature if no internal energy production
+
+    }
+    else
+    {
+      qdisk.ave_freq[n] = 0.0;
+      qdisk.t_hit[n] = 0.0;
+      qdisk.w[n] = 0.0;
+      ratio[n] = 0.0;
+      theat = 0;
     }
     ttot = pow (qdisk.t[n], 4) + pow (theat, 4);
     ttot = pow (ttot, 0.25);
     fprintf (qptr,
-             "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %5d %8.3e %8.3e %8.3e %8.3e %8.3e\n",
+             "%9.4e %9.4e %0.4e %8.3e %8.3e  %8.3e %10d %8.3e %10d %8.3e %8.3e %8.3e %10.3e %8.3e\n",
              qdisk.r[n], qdisk.v[n], zdisk (qdisk.r[n]), qdisk.t[n], qdisk.g[n],
-             qdisk.heat[n], qdisk.nhit[n], qdisk.heat[n] * NRINGS / ztot, theat, qdisk.t_hit[n], qdisk.w[n], ttot);
+             qdisk.emit[n], qdisk.nphot[n], qdisk.heat[n], qdisk.nhit[n], ratio[n], theat, qdisk.t_hit[n], qdisk.w[n], ttot);
   }
 
   fclose (qptr);
diff --git a/source/disk_photon_gen.c b/source/disk_photon_gen.c
index 9d5f12cd1..826f3d51d 100644
--- a/source/disk_photon_gen.c
+++ b/source/disk_photon_gen.c
@@ -189,7 +189,7 @@ photo_gen_disk (p, weight, f1, f2, spectype, istart, nphot)
     }
     else if (spectype == SPECTYPE_BB_FCOL)
     {
-      fcol = disk_colour_correction(disk.t[nring]);
+      fcol = disk_colour_correction (disk.t[nring]);
       p[i].freq = planck (fcol * disk.t[nring], freqmin, freqmax);
     }
     else if (spectype == SPECTYPE_UNIFORM)
diff --git a/source/estimators_simple.c b/source/estimators_simple.c
index 27359ba05..c85f55f98 100644
--- a/source/estimators_simple.c
+++ b/source/estimators_simple.c
@@ -16,6 +16,7 @@
 
 #include "atomic.h"
 #include "python.h"
+#include <gsl/gsl_integration.h>
 
 
 /**********************************************************/
@@ -471,6 +472,83 @@ update_force_estimators (xplasma, p, phot_mid, ds, w_ave, ndom, z, frac_ff, frac
 }
 
 
+typedef struct
+{
+  double T;                     // Temperature in Kelvin
+} planck_params;
+
+// Planck spectral radiance function
+double
+planck_spectral_radiance (double nu, void *params)
+{
+  planck_params *p = (planck_params *) params;
+  double T = p->T;
+  double exponent = PLANCK * nu / (BOLTZMANN * T);
+  return (2 * PLANCK * pow (nu, 3) / pow (VLIGHT, 2)) / (exp (exponent) - 1);
+}
+
+// Function to calculate nu * J(nu)
+double
+nu_times_radiance (double nu, void *params)
+{
+  return nu * planck_spectral_radiance (nu, params);
+}
+
+// Function to calculate mean frequency over a specified range
+double
+mean_frequency_nu_range (double T, double nu_min, double nu_max)
+{
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc (1000);
+
+  double numerator, denominator, error;
+  planck_params params = { T };
+
+  gsl_function F;
+  F.function = &nu_times_radiance;
+  F.params = &params;
+
+  gsl_integration_qags (&F, nu_min, nu_max, 0, 1e-7, 1000, w, &numerator, &error);
+
+  F.function = &planck_spectral_radiance;
+
+  gsl_integration_qags (&F, nu_min, nu_max, 0, 1e-7, 1000, w, &denominator, &error);
+
+  gsl_integration_workspace_free (w);
+
+  return numerator / denominator;
+}
+
+// Function to estimate temperature given a mean frequency and a frequency range
+double
+estimate_temperature_from_mean_frequency (double mean_nu_target, double nu_min, double nu_max, double initial_guess)
+{
+  double T_guess = initial_guess;
+  double tol = 1e-6;            // Convergence tolerance
+  int max_iter = 100;           // Maximum number of iterations
+  double mean_nu;
+
+
+  for (int i = 0; i < max_iter; i++)
+  {
+    // Calculate the mean frequency for the current temperature guess
+    mean_nu = mean_frequency_nu_range (T_guess, nu_min, nu_max);
+
+    // Adjust the temperature guess
+    double T_new = T_guess * (mean_nu_target / mean_nu);
+
+    // Check for convergence
+    if (fabs (T_new - T_guess) < tol)
+    {
+      return T_new;
+    }
+
+    T_guess = T_new;
+  }
+
+  return T_guess;               // Return the estimated temperature
+}
+
+
 /**********************************************************/
 /**
  * @brief Normalise simple estimators and cell spectra
@@ -496,6 +574,10 @@ update_force_estimators (xplasma, p, phot_mid, ds, w_ave, ndom, z, frac_ff, frac
  *
  **********************************************************/
 
+/* switch to choose how to calculate T_r. This could be made to depend on the ionization mode
+   if desired, or it may be that TRUE is the desired behaviour long-term */
+#define BAND_CORRECTED_TRAD FALSE
+
 int
 normalise_simple_estimators (xplasma)
      PlasmaPtr xplasma;
@@ -527,7 +609,19 @@ normalise_simple_estimators (xplasma)
     xplasma->j_scatt /= (4. * PI * invariant_volume_time);
 
     xplasma->t_r_old = xplasma->t_r;    // Store the previous t_r in t_r_old immediately before recalculating
-    radiation_temperature = xplasma->t_r = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832);
+
+    /* the method of calculation of the band corrected radiation temperature depends on 
+       the flag BAND_CORRECTED_TRAD -- see issue #1097 */
+    if (BAND_CORRECTED_TRAD == FALSE)
+    {
+      radiation_temperature = xplasma->t_r = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832);
+    }
+    else 
+    {
+      radiation_temperature = xplasma->t_r =
+        estimate_temperature_from_mean_frequency (xplasma->ave_freq, xband.f1[0], xband.f2[xband.nbands - 1], radiation_temperature);
+    }
+
     xplasma->w =
       PI * xplasma->j / (STEFAN_BOLTZMANN * radiation_temperature * radiation_temperature * radiation_temperature * radiation_temperature);
 
diff --git a/source/ionization.c b/source/ionization.c
index 11d41c838..127cc3e92 100644
--- a/source/ionization.c
+++ b/source/ionization.c
@@ -572,7 +572,8 @@ calc_te (PlasmaPtr xplasma, double tmin, double tmax)
  * at a specifc temperature.  
  *
  * @details
- * This routine is used to estiamate a new temperature for a cell given the
+ * This routine is used in the process of estimating a new temperature for a cell
+ * given the
  * heating of the cell in the previous ionization cycle.  When 0 the heating
  * and cooling are matched.
  *
@@ -583,6 +584,13 @@ calc_te (PlasmaPtr xplasma, double tmin, double tmax)
  *
  *
  * ### Notes ###
+ *
+ * This routine is poorly named; what it does is simply calculate the
+ * sum the heating that occurred during the current cycle, and calculate
+ * the cooling for the input temperature given the current abundances
+ * it returns the difference between these two numbers.  It does update
+ * t_e in the plasma cell.
+ *
  * The abundances of ions in the cell are not modified.  Results are stored
  * in the cell of interest.  This routine is used in connection with a zero
  * finding routine
diff --git a/source/photon_gen.c b/source/photon_gen.c
index ae6e8cb12..2ebb4f228 100644
--- a/source/photon_gen.c
+++ b/source/photon_gen.c
@@ -159,6 +159,10 @@ define_phot (p, f1, f2, nphot_tot, ioniz_or_extract, iwind, freq_sampling)
 
         iphot_start += xband.nphot[n];
       }
+      else
+      {
+        Error ("photon_gen: No photons for band %d\n", n);
+      }
 
     }
 
diff --git a/source/python.h b/source/python.h
index 41aac6daf..4b7424e9a 100644
--- a/source/python.h
+++ b/source/python.h
@@ -755,12 +755,13 @@ struct xdisk
   double t[NRINGS];             /**< The temperature at the middle of the annulus */
   double g[NRINGS];             /**< The gravity at the middle of the annulus */
   double v[NRINGS];             /**< The velocity at the middle of the annulus */
+  double emit[NRINGS];          /**< The radiative energy of the photons emitted from the disk */
   double heat[NRINGS];          /**< The total energy flux of photons hitting each annulus */
   double ave_freq[NRINGS];      /**< The flux weighted average of frequency of photons hitting each annulus */
-  double w[NRINGS];             /**< The radiative weight of the photons that hit the disk */
-  double t_hit[NRINGS];         /**< The effective T of photons hitting the disk */
-  int nphot[NRINGS];            /**< The number of photons created in each annulus */
-  int nhit[NRINGS];             /**< The number of photons which hit each annulus */
+  double w[NRINGS];             /**< The weight relative to a BB of the photons that hit the disk */
+  double t_hit[NRINGS];         /**< The effective T of photons hitting the disk, based on the average frequency */
+  int nphot[NRINGS];        /**< The number of disk photons created in each annulus */
+  int nhit[NRINGS];         /**< The number of photons which hit each annulus */
 };
 extern struct xdisk disk, qdisk;   /**<  disk defines zones in the disk which in a specified frequency band emit equal amounts
                                    of radiation. disk gets reinitialized whenever the frequency interval of interest
@@ -1236,6 +1237,11 @@ extern int size_Jbar_est, size_gamma_est, size_alpha_est;
 #define NEBULARMODE_MATRIX_SPECTRALMODEL 9      /**< matrix solver spectral model */
 #define NEBULARMODE_MATRIX_ESTIMATORS 10        /**<  matrix solver spectral model */
 
+#define NEBULARMODE_MATRIX_MULTISHOT   11    /**<  matrix solver spectral model based on power laws which
+                                          * updates T_e multiple times before arriving at a final
+                                          * solution 
+                                          */
+
 // modes for the wind_luminosity routine
 #define MODE_OBSERVER_FRAME_TIME 0
 #define MODE_CMF_TIME 1
diff --git a/source/run.c b/source/run.c
index cd4956072..93d8c92a6 100644
--- a/source/run.c
+++ b/source/run.c
@@ -48,8 +48,8 @@ int
 calculate_ionization (restart_stat)
      int restart_stat;
 {
-  int n, nn;
-  double zz, z_abs_all, z_abs_all_orig, z_orig[N_ISTAT], z_abs[N_ISTAT], z_else, z_else_orig, ztot;
+  int nn;
+  double zz, z_abs_all, z_abs_all_orig, z_orig[N_ISTAT], z_abs[N_ISTAT], z_else, z_else_orig;
   double radiated[20], radiated_orig[20];
   int nphot_istat[N_ISTAT];
   WindPtr w;
@@ -146,7 +146,6 @@ calculate_ionization (restart_stat)
 
     geo.n_ioniz = 0.0;
     geo.cool_tot_ioniz = 0.0;
-    ztot = 0.0;                 /* ztot is the luminosity of the disk multipled by the number of cycles, which is used by save_disk_heating */
 
     if (!geo.wind_radiation || (geo.wcycle == 0 && geo.run_type != RUN_TYPE_PREVIOUS))
       iwind = -1;               /* Do not generate photons from wind */
@@ -190,11 +189,9 @@ calculate_ionization (restart_stat)
     geo.lum_star_back = 0;
     geo.lum_disk_back = 0;
 
+    /* Prepare qdisk for recording photon pages */
+    qdisk_reinit (p);
 
-    for (n = 0; n < NRINGS; n++)
-    {
-      qdisk.heat[n] = qdisk.nphot[n] = qdisk.w[n] = qdisk.ave_freq[n] = 0;
-    }
 
 
     zz = 0.0;
@@ -205,7 +202,6 @@ calculate_ionization (restart_stat)
 
     Log ("!!python: Total photon luminosity before transphot %18.12e\n", zz);
     Log_flush ();
-    ztot += zz;                 /* Total luminosity in all cycles, used for calculating disk heating */
 
     /* kbf_need determines how many & which bf processes one needs to considere.  It was introduced
      * as a way to speed up the program.  It has to be recalculated evey time one changes
@@ -323,7 +319,14 @@ calculate_ionization (restart_stat)
 #endif
       if (geo.disk_type != DISK_NONE)
       {
-        qdisk_save (files.disk, ztot);
+        qdisk_save (files.disk, 1);
+        if (modes.make_tables)
+        {
+          strcpy (dummy, "");
+          sprintf (dummy, "diag_%.100s/%.100s.%02d.disk.diag", files.root, files.root, geo.wcycle + 1);
+          qdisk_save (dummy, 0);
+
+        }
       }
 #ifdef MPI_ON
     }
diff --git a/source/templates.h b/source/templates.h
index 14892e70c..e90e95a94 100644
--- a/source/templates.h
+++ b/source/templates.h
@@ -163,7 +163,8 @@ double disk_colour_correction(double t);
 /* disk_init.c */
 double disk_init(double rmin, double rmax, double m, double mdot, double freqmin, double freqmax, int ioniz_or_extract, double *ftot);
 int qdisk_init(double rmin, double rmax, double m, double mdot);
-int qdisk_save(char *diskfile, double ztot);
+int qdisk_reinit(PhotPtr p);
+int qdisk_save(char *diskfile, int ichoice);
 int read_non_standard_disk_profile(char *tprofile);
 /* disk_photon_gen.c */
 int photo_gen_disk(PhotPtr p, double weight, double f1, double f2, int spectype, int istart, int nphot);
@@ -200,6 +201,10 @@ double alpha_st_e_integrand(double freq, void *params);
 int update_banded_estimators(PlasmaPtr xplasma, PhotPtr p, double ds, double w_ave, int ndom);
 int update_flux_estimators(PlasmaPtr xplasma, PhotPtr phot_mid, double ds_obs, double w_ave, int ndom);
 int update_force_estimators(PlasmaPtr xplasma, PhotPtr p, PhotPtr phot_mid, double ds, double w_ave, int ndom, double z, double frac_ff, double frac_auger, double frac_tot);
+double planck_spectral_radiance(double nu, void *params);
+double nu_times_radiance(double nu, void *params);
+double mean_frequency_nu_range(double T, double nu_min, double nu_max);
+double estimate_temperature_from_mean_frequency(double mean_nu_target, double nu_min, double nu_max, double initial_guess);
 int normalise_simple_estimators(PlasmaPtr xplasma);
 void update_persistent_directional_flux_estimators(int nplasma, double flux_persist_scale);
 /* extract.c */
@@ -355,6 +360,8 @@ int invert_matrix(double *matrix, double *inverted_matrix, int num_rows);
 /* matrix_ion.c */
 int matrix_ion_populations(PlasmaPtr xplasma, int mode);
 int populate_ion_rate_matrix(double rate_matrix[nions][nions], double pi_rates[nions], double inner_rates[n_inner_tot], double rr_rates[nions], double b_temp[nions], double xne, double nh1, double nh2);
+/* matrix_ion2.c */
+int matrix_ion_populations2(PlasmaPtr xplasma, int mode);
 /* models_extern_init.c */
 /* para_update.c */
 int get_parallel_nrange(int rank, int ntotal, int nproc, int *my_nmin, int *my_nmax);
diff --git a/source/windsave2fits.c b/source/windsave2fits.c
new file mode 100644
index 000000000..f56e2e186
--- /dev/null
+++ b/source/windsave2fits.c
@@ -0,0 +1,700 @@
+
+/***********************************************************/
+/** @file  windsave2fits.c 
+ * @author ksl
+ * @date   August,2024   
+ *
+ * @brief  Routines to write variious portions of  in a wind structure 
+ * to a fits file
+ *
+ *###Notes###
+ * This routine was written primarily to better understand
+ * how well the spectra we use for estimating ionization rates
+ * model the actual cell spectra; but it can be rather
+ * straightforwardly adatpted to create images or table
+ * of other portions of a windsave file, that involve
+ * large amounts of data.
+ *
+ *
+ ***********************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "atomic.h"
+#include "python.h"
+#include "fitsio.h"
+
+
+/* The first portions of the program contain some utility routines 
+   for working with cfitsio */
+
+char inroot[LINELENGTH], outroot[LINELENGTH], model_file[LINELENGTH], folder[LINELENGTH];
+int model_flag, ksl_flag, cmf2obs_flag, obs2cmf_flag;
+
+double line_matom_lum_single (double lum[], PlasmaPtr xplasma, int uplvl);
+int line_matom_lum (int uplvl);
+int create_matom_level_map ();
+
+// Define a structure to hold spectral data
+typedef struct
+{
+  int num_wavelengths;
+  int num_spectra;
+  float **data;                 // 2D array to hold spectra data
+} Spectra;
+
+/* Next functions are utilities that it should be possitle to call repeatedly*/
+
+/**********************************************************/
+/**
+ * @brief      Function to perpare 2d image data from a 2d array
+ *
+ * @details The routine simply returns a flattened version of the 2d array
+ * ### Notes ###
+ *
+ *
+ **********************************************************/
+
+float *
+prepare_image_data (float **data, int width, int height)
+{
+  // Allocate memory for a contiguous 1D array to store the 2D image data
+  float *image_data = calloc (width * height, sizeof (float));
+  if (!image_data)
+  {
+    fprintf (stderr, "Memory allocation error\n");
+    return NULL;
+  }
+
+  // Flatten the 2D data into a 1D array for FITS writing
+  for (int i = 0; i < height; i++)
+  {
+    for (int j = 0; j < width; j++)
+    {
+      image_data[i * width + j] = data[i][j];
+    }
+  }
+
+  return image_data;
+}
+
+
+
+/**********************************************************/
+/**
+ * @brief      Function to write a 2d image data to a fits file with
+ * a ginve extension name
+ *
+ * @details The routine normally returns 0                                
+ * ### Notes ###
+ * This routine is called after prepare_image data, which will
+ * have flattened and converted the original array to a float
+ *
+ *
+ **********************************************************/
+
+int
+write_image_extension (fitsfile *fptr, float *image_data, int width, int height, const char *extname)
+{
+  int status = 0;               // CFITSIO status value MUST be initialized to zero
+  long naxes[2];                // Array to store dimension sizes
+
+  // Set the dimensions of the 2D image: [width, height]
+  naxes[0] = width;             // Number of columns
+  naxes[1] = height;            // Number of rows
+
+  // Create a new image extension of type FLOAT_IMG
+  if (fits_create_img (fptr, FLOAT_IMG, 2, naxes, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Set the extension name using the EXTNAME keyword
+  if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) extname, NULL, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Write the image data to the FITS extension
+  if (fits_write_img (fptr, TFLOAT, 1, width * height, image_data, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  return 0;
+}
+
+
+
+
+/**********************************************************/
+/**
+ * @brief      Function to write a 1d image extension with
+ * a ginve extension name
+ *
+ * @details The routine normally returns 0                                   
+ * ### Notes ###
+ *
+ *
+ **********************************************************/
+
+int
+write_1d_image_extension (fitsfile *fptr, float *data, long length, const char *extname)
+{
+  int status = 0;               // CFITSIO status value MUST be initialized to zero
+  long naxes[1];                // Array to store the dimension size (1D array)
+
+  naxes[0] = length;            // Length of the 1D array
+
+  // Create a new image extension of type FLOAT_IMG
+  if (fits_create_img (fptr, FLOAT_IMG, 1, naxes, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Set the extension name using the EXTNAME keyword
+  if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) extname, NULL, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Write the 1D array data to the FITS image extension
+  if (fits_write_img (fptr, TFLOAT, 1, length, data, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  return 0;
+}
+
+
+/* This routine is not longer used but has been kept for reference
+// Function to write a binary table with mixed data types to a FITS extension
+int
+write_mixed_table_extension (fitsfile *fptr, int num_rows, int num_columns, int *int_data, float *float_data, const char *extname)
+{
+  int status = 0;
+  // Define the names, formats, and units for each column
+  char *ttype[] = { "INT_COL", "FLOAT_COL" };   // Column names
+  char *tform[] = { "J", "E" }; // Formats: 'J' for integer, 'E' for float
+  char *tunit[] = { "", "" };   // Units
+
+  // Create a new binary table extension
+  if (fits_create_tbl (fptr, BINARY_TBL, num_rows, num_columns, ttype, tform, tunit, NULL, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Set the extension name using the EXTNAME keyword
+  if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) extname, NULL, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Write the integer data to the first column
+  if (fits_write_col (fptr, TINT, 1, 1, 1, num_rows, int_data, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Write the float data to the second column
+  if (fits_write_col (fptr, TFLOAT, 2, 1, 1, num_rows, float_data, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  return 0;
+}
+*/
+
+/* This ends the utility routines */
+
+/**********************************************************/
+/**
+ * @brief      write_spectral_model_table extension to the
+ * fits file
+ *
+ * @details The routine normally returns 0                                   
+ * ### Notes ###
+ *
+ *
+ **********************************************************/
+
+int
+write_spectra_model_table (fitsfile *fptr)
+{
+
+  /* We need one row for each cell and band so
+     our table will be NPLASMA*nbands long
+   */
+
+  int i, j, k;
+  int num_rows, num_cols;
+  int nbands;
+  int *ichoice, *iplasma, *iband, *nxtot;
+  float *exp_w, *exp_temp, *pl_log_w, *pl_alpha;
+
+  nbands = geo.nxfreq;
+  num_rows = nbands * NPLASMA;
+  printf ("xtest %d %d\n", nbands, NPLASMA);
+  num_cols = 8;                 // for now
+  char *table_name = "spec_model";
+
+  ichoice = calloc (num_rows, sizeof (int *));
+  iplasma = calloc (num_rows, sizeof (int *));
+  iband = calloc (num_rows, sizeof (int *));
+  exp_w = calloc (num_rows, sizeof (float *));
+  exp_temp = calloc (num_rows, sizeof (float *));
+  pl_log_w = calloc (num_rows, sizeof (float *));
+  pl_alpha = calloc (num_rows, sizeof (float *));
+  nxtot = calloc (num_rows, sizeof (int *));
+
+  k = 0;
+  for (i = 0; i < NPLASMA; i++)
+  {
+    for (j = 0; j < nbands; j++)
+    {
+      iplasma[k] = i;
+      iband[k] = j;
+      ichoice[k] = plasmamain[i].spec_mod_type[j];
+      exp_w[k] = plasmamain[i].exp_w[j];
+      exp_temp[k] = plasmamain[i].exp_temp[j];
+      pl_log_w[k] = plasmamain[i].pl_log_w[j];
+      pl_alpha[k] = plasmamain[i].pl_alpha[j];
+      nxtot[k] = plasmamain[i].nxtot[j];
+      printf ("nxtot %d\n", plasmamain[i].nxtot[j]);
+      k++;
+    }
+  }
+
+
+  int status = 0;
+  // Define the names, formats, and units for each column
+  char *ttype[] = { "nplasma", "band", "spec_mod_type", "exp_w", "exp_temp", "pl_log_w", "pl_alpha", "nxtot" }; // Column names
+  char *tform[] = { "J", "J", "J", "E", "E", "E", "E", "J" };   // Formats: 'J' for integer, 'E' for float
+  char *tunit[] = { "", "", "", "", "", "", "", "" };   // Units
+
+  // Create a new binary table extension
+  if (fits_create_tbl (fptr, BINARY_TBL, num_rows, num_cols, ttype, tform, tunit, NULL, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Set the extension name using the EXTNAME keyword
+  if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) table_name, NULL, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  printf ("num_rows: %d\n", num_rows);
+  for (i = 0; i < num_rows; i++)
+  {
+    printf ("%d %d %d\n", i, iplasma[i], ichoice[i]);
+  }
+
+  // Write the integer data to the first column
+  if (fits_write_col (fptr, TINT, 1, 1, 1, num_rows, iplasma, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+
+  // Write the float data to the second column
+  if (fits_write_col (fptr, TINT, 2, 1, 1, num_rows, iband, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Write the float data to the third  column
+  if (fits_write_col (fptr, TINT, 3, 1, 1, num_rows, ichoice, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+
+  // Write the float data to the fourth column
+  if (fits_write_col (fptr, TFLOAT, 4, 1, 1, num_rows, exp_w, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+
+
+  // Write the float data to the fifth columnmn
+  if (fits_write_col (fptr, TFLOAT, 5, 1, 1, num_rows, exp_temp, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+
+  // Write the float data to the sixth  column
+  if (fits_write_col (fptr, TFLOAT, 6, 1, 1, num_rows, pl_log_w, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+
+  // Write the float data to the seventh column
+  if (fits_write_col (fptr, TFLOAT, 7, 1, 1, num_rows, pl_alpha, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+  // Write the integer data to the eightth column
+  if (fits_write_col (fptr, TINT, 8, 1, 1, num_rows, nxtot, &status))
+  {
+    fits_report_error (stderr, status);
+    return status;
+  }
+
+
+  return 0;
+}
+
+
+
+
+/**********************************************************/
+/**
+ * @brief      write_a fits file that contains information
+ * releveant to how well the cell spectra are modelled
+ * fits file
+ *
+ * @details The routine normally returns 0                                   
+ * ### Notes ###
+ *
+ * This is the main routine, in the sense that everyting is
+ * controlled frrom here
+ *
+ **********************************************************/
+
+
+
+int
+make_spec (inroot)
+     char *inroot;
+{
+
+  fitsfile *fptr;               // Pointer to the FITS file
+  int status = 0;               // CFITSIO status value MUST be initialized to zero
+  char outfile[LINELENGTH];
+
+  /* The obscure exclamation mark at the beginning of the name means yu want to overwrie 
+     and existing file, if it exists */
+
+  sprintf (outfile, "!%s_cellspec.fits", inroot);
+  printf ("outfile %s\n", outfile);
+
+
+
+  // Create a new FITS file
+  //if (fits_create_file (&fptr, "!multi_extension.fits", &status))
+  if (fits_create_file (&fptr, outfile, &status))
+  {
+    fits_report_error (stderr, status); // Print any error message
+    return status;
+  }
+
+// Create an empty primary HDU
+  if (fits_create_img (fptr, FLOAT_IMG, 0, NULL, &status))
+  {
+    fits_report_error (stderr, status);
+  }
+
+  printf ("Hello World %s \n", inroot);
+  printf ("Plasma  %d \n", NPLASMA);
+  printf ("NBINS in spec %d \n", NBINS_IN_CELL_SPEC);
+
+
+  Spectra spectra;
+  spectra.num_spectra = NPLASMA;
+  spectra.num_wavelengths = NBINS_IN_CELL_SPEC;
+  spectra.data = calloc (spectra.num_spectra, sizeof (float *));
+
+  for (int i = 0; i < spectra.num_spectra; i++)
+  {
+    spectra.data[i] = calloc (spectra.num_wavelengths, sizeof (float));
+
+    // Fill with sample data
+    for (int j = 0; j < spectra.num_wavelengths; j++)
+    {
+      // spectra.data[i][j] = (float) (i + j);
+      spectra.data[i][j] = (float) plasmamain[i].cell_spec_flux[j];
+    }
+  }
+
+  // Prepare image data for the first extension
+  float *image_data1 = prepare_image_data (spectra.data, spectra.num_wavelengths, spectra.num_spectra);
+  if (!image_data1)
+  {
+    // Free allocated memory before exiting
+    for (int i = 0; i < spectra.num_spectra; i++)
+    {
+      free (spectra.data[i]);
+    }
+    free (spectra.data);
+    return 1;
+  }
+
+
+  status = write_image_extension (fptr, image_data1, spectra.num_wavelengths, spectra.num_spectra, "Jnu");
+
+  Spectra freq;
+  freq.num_spectra = 1;
+  freq.num_wavelengths = NBINS_IN_CELL_SPEC;
+  freq.data = calloc (freq.num_spectra, sizeof (float *));
+  freq.data[0] = calloc (freq.num_wavelengths, sizeof (float));
+
+  for (int i = 0; i < NBINS_IN_CELL_SPEC; i++)
+  {
+    freq.data[0][i] = (float) pow (10., geo.cell_log_freq_min + i * geo.cell_delta_lfreq);
+  }
+
+
+  float *image_data2 = prepare_image_data (freq.data, NBINS_IN_CELL_SPEC, 1);
+/*
+  double freq[2000][1];
+  int i;
+
+  for (i = 0; i < NBINS_IN_CELL_SPEC; i++)
+  {
+    freq[i][0] = (float) pow (10., geo.cell_log_freq_min + i * geo.cell_delta_lfreq);
+  }
+
+  float *image_data2 = prepare_image_data (freq, NBINS_IN_CELL_SPEC, 1);
+*/
+  // status = write_image_extension (fptr, image_data2, NBINS_IN_CELL_SPEC, 1, "nu");
+  status = write_1d_image_extension (fptr, image_data2, NBINS_IN_CELL_SPEC, "nu");
+
+  /* Elimainate this for now
+
+     // Prepare data for the binary table
+     int num_rows = 5;
+     int int_data[] = { 1, 2, 3, 4, 5 };   // Example integer data
+     float float_data[] = { 0.1, 0.2, 0.3, 0.4, 0.5 };     // Example float data
+
+     write_mixed_table_extension (fptr, num_rows, 2, int_data, float_data, "Table_Ext");
+   */
+
+
+/* Now create a 1d extension hold the frequencies of the broad bands */
+
+  int i;
+  float xfreq[100];
+
+  for (i = 0; i <= geo.nxfreq; i++)
+  {
+    xfreq[i] = geo.xfreq[i];
+  }
+
+  status = write_1d_image_extension (fptr, xfreq, geo.nxfreq + 1, "nu_model");
+
+/* Now write my table */
+
+  write_spectra_model_table (fptr);
+
+  if (fits_close_file (fptr, &status))
+  {
+    fits_report_error (stderr, status);
+  }
+
+  return 0;
+}
+
+/**********************************************************/
+/**
+ * @brief      parses the command line options
+ *
+ * @param [in]  int  argc   the number of command line arguments
+ * @param [in]  char *  argv[]   The command line arguments
+ *
+ *
+ * ###Notes###
+ *
+ * The general purpose of each of the command line options
+ * should be fairly obvious from reading the code.
+ *
+ *
+ * Although this routine uses the standard Log and Error commands
+ * the diag files have not been created yet and so this information
+ * is really simply written to the terminal.
+ * 
+ * WARNING: this has not been updeated. Currently there
+ * routine expects only the rootname of a windsave file
+ *
+ *
+ **********************************************************/
+
+int
+xparse_command_line (argc, argv)
+     int argc;
+     char *argv[];
+{
+  int j = 0;
+  int i;
+  char dummy[LINELENGTH];
+  int mkdir ();
+  char *fgets_rc;
+
+
+  sprintf (outroot, "%s", "new");
+
+  model_flag = ksl_flag = obs2cmf_flag = cmf2obs_flag = 0;
+
+  if (argc == 1)
+  {
+    printf ("Parameter file name (e.g. my_model.pf, or just my_model):");
+    fgets_rc = fgets (dummy, LINELENGTH, stdin);
+    if (!fgets_rc)
+    {
+      printf ("Input rootname is NULL or invalid\n");
+      exit (1);
+    }
+    get_root (inroot, dummy);
+  }
+  else
+  {
+
+    for (i = 1; i < argc; i++)
+    {
+      if (strcmp (argv[i], "-out_root") == 0)
+      {
+        if (sscanf (argv[i + 1], "%s", dummy) != 1)
+        {
+          printf ("python: Expected out_root after -out_root switch\n");
+          exit (0);
+        }
+
+        get_root (outroot, dummy);
+        i++;
+        j = i;
+
+      }
+      if (strcmp (argv[i], "-model_file") == 0)
+      {
+        if (sscanf (argv[i + 1], "%s", dummy) != 1)
+        {
+          printf ("python: Expected a model file containing density, velocity and temperature after -model_file switch\n");
+          exit (0);
+        }
+        get_root (model_file, dummy);
+        i++;
+        j = i;
+        printf ("got a model file %s\n", model_file);
+        model_flag = 1;
+      }
+      else if (strcmp (argv[i], "-ksl") == 0)
+      {
+        printf ("Carrying out a simple hard wired ion modification\n");
+        ksl_flag = 1;
+      }
+      else if (strcmp (argv[i], "--dry-run") == 0)
+      {
+        modes.quit_after_inputs = 1;
+        j = i;
+      }
+      else if (strcmp (argv[i], "-cmf") == 0)
+      {
+        obs2cmf_flag = 1;
+      }
+      else if (strcmp (argv[i], "-obs") == 0)
+      {
+        cmf2obs_flag = 1;
+      }
+      else if (strncmp (argv[i], "-", 1) == 0)
+      {
+        printf ("python: Unknown switch %s\n", argv[i]);
+        exit (0);
+      }
+    }
+
+    /* The last command line variable is always the windsave file */
+
+    if (j + 1 == argc)
+    {
+      printf ("All of the command line has been consumed without specifying a parameter file name, so exiting\n");
+      exit (1);
+    }
+    strcpy (dummy, argv[argc - 1]);
+    get_root (inroot, dummy);
+
+  }
+
+  return (0);
+}
+
+
+
+/**********************************************************/
+/**
+ * @brief      the main routine which carries out the effort
+ *
+ * @param [in]  int  argc   the number of command line arguments
+ * @param [in]  char *  argv[]   The command line arguments
+ *
+ *
+ * ###Notes###
+ *
+ * This routine oversees the effort.  The basic steps are
+ *
+ * - parse the command line to get the names of files
+ * - call the routine to craetate fits file
+ *
+ *
+ **********************************************************/
+
+
+
+int
+main (argc, argv)
+     int argc;
+     char *argv[];
+{
+
+  char infile[LINELENGTH], outfile[LINELENGTH];
+  FILE *fopen ();
+  int mkdir ();
+
+
+  xparse_command_line (argc, argv);
+
+  sprintf (infile, "%.150s.wind_save", inroot);
+  sprintf (outfile, "%.150s.txt", inroot);
+
+  printf ("Reading %s and writing to %s\n", infile, outfile);
+
+  zdom = calloc (MAX_DOM, sizeof (domain_dummy));
+  if (zdom == NULL)
+  {
+    printf ("Unable to allocate memory for domain\n");
+    return EXIT_FAILURE;
+  }
+
+  wind_read (infile);
+
+  make_spec (inroot);
+
+  return (0);
+}
diff --git a/source/windsave2table_sub.c b/source/windsave2table_sub.c
index ba3e72f94..9a5c3b9f7 100644
--- a/source/windsave2table_sub.c
+++ b/source/windsave2table_sub.c
@@ -498,10 +498,10 @@ create_heat_table (ndom, rootname)
   strcpy (column_name[25], "heat_shock");
 
   c[26] = get_one (ndom, "heat_lines_macro");
-  strcpy (column_name[26], "heat_lines_macro");
+  strcpy (column_name[26], "ht_ln_macro");
 
   c[27] = get_one (ndom, "heat_photo_macro");
-  strcpy (column_name[27], "heat_photo_macro");
+  strcpy (column_name[27], "ht_ph_macro");
 
   /* This should be the maximum number above +1 */
   ncols = 28;