Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile	(revision 37629)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile	(revision 37630)
@@ -164,4 +164,6 @@
 $(SRC)/shterms.$(ARCH).o	   \
 $(SRC)/shfit.$(ARCH).o	   \
+$(SRC)/shdot.$(ARCH).o	   \
+$(SRC)/shapply.$(ARCH).o	   \
 $(SRC)/wd.$(ARCH).o		   \
 $(SRC)/write_vectors.$(ARCH).o	   \
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c	(revision 37629)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c	(revision 37630)
@@ -152,4 +152,6 @@
 int shterms          PROTO((int, char **));
 int shfit            PROTO((int, char **));
+int shdot            PROTO((int, char **));
+int shapply          PROTO((int, char **));
 int wd               PROTO((int, char **));
 int write_vectors    PROTO((int, char **));
@@ -323,4 +325,6 @@
   {1, "shterms",      shterms,          "Spherical Harmonics terms"},
   {1, "shfit",        shfit,            "Spherical Harmonics fits"},
+  {1, "shdot",        shdot,            "Spherical Harmonics dot product"},
+  {1, "shapply",      shapply,          "Spherical Harmonics fit application"},
   {1, "wd",           wd,               "write an image to a file"},
   {1, "write",        write_vectors,    "write vectors to datafile"},
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shapply.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shapply.c	(revision 37630)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shapply.c	(revision 37630)
@@ -0,0 +1,63 @@
+# include "data.h"
+
+int shapply (int argc, char **argv) {
+  
+  int i, j;
+  Vector *Rvec, *Dvec, *Frvec, *Fivec, *Vrvec, *Vivec;
+
+  if (argc != 8) {
+    gprint (GP_ERR, "USAGE: shapply R D lmax Fr Fi Vr Vi\n");
+    gprint (GP_ERR, "  apply the Ylm values (Fr,Fi) at the given R,D values to get the fit at those positions\n");
+    return (FALSE);
+  }
+
+  if ((Rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Dvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if (Rvec->Nelements != Dvec->Nelements) {
+    gprint (GP_ERR, "R and D sizes do not match\n");
+    return FALSE;
+  }
+
+  int lmax = atoi(argv[3]);
+
+  if ((Frvec = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fivec = SelectVector (argv[5], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if (Frvec->Nelements != Fivec->Nelements) { 
+    gprint (GP_ERR, "Fr and Fi sizes do not match\n");
+    return FALSE;
+  }
+
+  if ((Vrvec = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Vivec = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  SHterms *terms = SHtermsInit (lmax);
+  if (Frvec->Nelements != terms->Nterms) {
+    gprint (GP_ERR, "Lmax (%d) terms (%d) does not match size of Fr,Fi vectors (%d) \n", lmax, terms->Nterms, Frvec->Nelements);
+    return FALSE;
+  }
+
+  ResetVector (Vrvec, OPIHI_FLT, Rvec->Nelements);
+  ResetVector (Vivec, OPIHI_FLT, Rvec->Nelements);
+
+  // measure the dot product \sum(F_i * Ylm_i)
+  for (i = 0; i < Rvec->Nelements; i++) {
+
+    SHtermsForRD (terms, Rvec->elements.Flt[i], Dvec->elements.Flt[i]); 
+
+    double *Fr = Frvec->elements.Flt;
+    double *Fi = Fivec->elements.Flt;
+
+    double Vr = 0.0;
+    double Vi = 0.0;
+    for (j = 0; j < terms->Nterms; j++) {
+      Vr += Fr[j] * terms->Fr[j];
+      Vi += Fi[j] * terms->Fi[j];
+    }
+    Vrvec->elements.Flt[i] = Vr;
+    Vivec->elements.Flt[i] = Vi;
+  }
+
+  SHtermsFree (terms);
+
+  return (TRUE);
+}
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shdot.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shdot.c	(revision 37630)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shdot.c	(revision 37630)
@@ -0,0 +1,69 @@
+# include "data.h"
+
+int shdot (int argc, char **argv) {
+  
+  int i, j;
+  Vector *Rvec, *Dvec, *Fvec, *Lvec, *Mvec, *Frvec, *Fivec;
+  double *Fr, *Fi;
+
+  if (argc != 9) {
+    gprint (GP_ERR, "USAGE: shdot R D value lmax Lout Mout Fr Fi\n");
+    gprint (GP_ERR, "  find the dot product of the scalar field value (at points R,D) to spherical harmonics up to the given lmax\n");
+    return (FALSE);
+  }
+
+  if ((Rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Dvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  int lmax = atoi(argv[4]);
+
+  if ((Lvec  = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Mvec  = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Frvec = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fivec = SelectVector (argv[8], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  SHterms *terms = SHtermsInit (lmax);
+
+  ALLOCATE_ZERO (Fr, double, terms->Nterms);
+  ALLOCATE_ZERO (Fi, double, terms->Nterms);
+
+  ResetVector (Frvec, OPIHI_FLT, terms->Nterms);
+  ResetVector (Fivec, OPIHI_FLT, terms->Nterms);
+  ResetVector (Lvec,  OPIHI_FLT, terms->Nterms);
+  ResetVector (Mvec,  OPIHI_FLT, terms->Nterms);
+
+  // measure the dot product \sum(F_i * Ylm_i)
+  for (i = 0; i < Rvec->Nelements; i++) {
+
+    SHtermsForRD (terms, Rvec->elements.Flt[i], Dvec->elements.Flt[i]); 
+
+    double Fv = Fvec->elements.Flt[i];
+
+    for (j = 0; j < terms->Nterms; j++) {
+      Fr[j] += Fv * terms->Fr[j];
+      Fi[j] += Fv * terms->Fi[j];
+    }
+  }
+
+  // for (j = 0; j < terms->Nterms; j++) {
+  //   fprintf (stderr, "%d : %d %d : %f %f\n", j, terms->l[j], terms->m[j], Fr[j], Fi[j]);
+  // }
+
+  for (j = 0; j < terms->Nterms; j++) {
+    Fr[j] /= Fvec->Nelements;
+    Fi[j] /= Fvec->Nelements;
+  }
+
+  for (j = 0; j < terms->Nterms; j++) {
+    // fprintf (stderr, "%d : %d %d : %f %f\n", j, terms->l[j], terms->m[j], 4*M_PI*Fr[j], 4*M_PI*Fi[j]);
+    Lvec[0].elements.Flt[j]  = terms->l[j];
+    Mvec[0].elements.Flt[j]  = terms->m[j];
+    Frvec[0].elements.Flt[j] = 4*M_PI*Fr[j];
+    Fivec[0].elements.Flt[j] = 4*M_PI*Fi[j];
+  }
+
+  SHtermsFree (terms);
+
+  return (TRUE);
+}
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c	(revision 37629)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c	(revision 37630)
@@ -3,5 +3,5 @@
 int shfit (int argc, char **argv) {
   
-  int i, j;
+  int i, j, k;
   Vector *Rvec, *Dvec, *Fvec, *Lvec, *Mvec, *Frvec, *Fivec;
   double *Fr, *Fi;
@@ -34,5 +34,17 @@
   ResetVector (Mvec,  OPIHI_FLT, terms->Nterms);
 
-  // measure the dot product \sum(F_i * Ylm_i)
+  double **Ar, **br, **Ai, **bi;
+  ALLOCATE (Ar, double *, terms->Nterms);
+  ALLOCATE (Ai, double *, terms->Nterms);
+  ALLOCATE (br, double *, terms->Nterms);
+  ALLOCATE (bi, double *, terms->Nterms);
+  for (i = 0; i < terms->Nterms; i++) {
+    ALLOCATE_ZERO (Ar[i], double, terms->Nterms);
+    ALLOCATE_ZERO (Ai[i], double, terms->Nterms);
+    ALLOCATE_ZERO (br[i], double, 1);
+    ALLOCATE_ZERO (bi[i], double, 1);
+  }
+
+  // measure the dot product \sum(F_i * Ylm_i) and the cross terms (\sum(Y_lm * Y_jk))
   for (i = 0; i < Rvec->Nelements; i++) {
 
@@ -42,16 +54,86 @@
 
     for (j = 0; j < terms->Nterms; j++) {
-      Fr[j] += Fv * terms->Fr[j];
-      Fi[j] += Fv * terms->Fi[j];
+      br[j][0] += Fv * terms->Fr[j];
+      bi[j][0] += Fv * terms->Fi[j];
+    }
+
+    for (j = 0; j < terms->Nterms; j++) {
+      for (k = j; k < terms->Nterms; k++) {
+	Ar[j][k] += terms->Fr[j] * terms->Fr[k];
+	Ai[j][k] += terms->Fi[j] * terms->Fi[k];
+      }	
     }
   }
 
-  // for (j = 0; j < terms->Nterms; j++) {
-  //   fprintf (stderr, "%d : %d %d : %f %f\n", j, terms->l[j], terms->m[j], Fr[j], Fi[j]);
-  // }
+  for (j = 1; j < terms->Nterms; j++) {
+    for (k = 0; k < j; k++) {
+      Ar[j][k] = Ar[k][j];
+      Ai[j][k] = Ai[k][j];
+    }	
+  }
 
   for (j = 0; j < terms->Nterms; j++) {
-    Fr[j] /= Fvec->Nelements;
-    Fi[j] /= Fvec->Nelements;
+    if (terms->m[j] == 0) {
+      Ai[j][j] = 1.0; // Im(Y_l0) == 0.0
+      bi[j][0] = 0.0;
+    } else {
+      if (terms->m[j] < 0) {
+	for (k = 0; k < terms->Nterms; k++) {
+	  Ar[j][k] = 0.0;
+	  Ar[k][j] = 0.0;
+	  Ai[j][k] = 0.0;
+	  Ai[k][j] = 0.0;
+	}
+	Ar[j][j] = 1.0;
+	Ai[j][j] = 1.0;
+	br[j][0] = 0.0;
+	bi[j][0] = 0.0;
+	// Re(Y_lm) = Re(Y_l(-m)) [m even], Re(Y_lm) = Re(-Y_l(-m)) [m odd]
+      }
+    }
+  }
+
+  if (0) {
+    fprintf (stderr, "--- Ar --- : br \n");
+    for (i = 0; i < terms->Nterms; i++) {
+      for (j = 0; j < terms->Nterms; j++) {
+	fprintf (stderr, "%10.6f ", Ar[i][j]);
+      }
+      fprintf (stderr, " : %10.6f\n", br[i][0]);
+    }
+
+    fprintf (stderr, "--- Ai --- : bi \n");
+    for (i = 0; i < terms->Nterms; i++) {
+      for (j = 0; j < terms->Nterms; j++) {
+	fprintf (stderr, "%10.6f ", Ai[i][j]);
+      }
+      fprintf (stderr, " : %10.6f\n", bi[i][0]);
+    }
+  }
+
+  fprintf (stderr, "--- Ar --- : br \n");
+  for (j = 0; j < terms->Nterms; j++) {
+    fprintf (stderr, "%10.6f ", Ar[j][j]);
+    fprintf (stderr, " : %10.6f\n", br[j][0]);
+  }
+
+  fprintf (stderr, "--- Ai --- : bi \n");
+  for (i = 0; i < terms->Nterms; i++) {
+    fprintf (stderr, "%10.6f : ", Ai[i][i]);
+    fprintf (stderr, " : %10.6f\n", bi[i][0]);
+  }
+
+  if (!dgaussjordan (Ar, br, terms->Nterms, 1)) {
+    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
+    return FALSE;
+  }
+  if (!dgaussjordan (Ai, bi, terms->Nterms, 1)) {
+    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
+    return FALSE;
+  }
+
+  for (j = 0; j < terms->Nterms; j++) {
+    Fr[j] = br[j][0] / (float) Frvec[0].Nelements;
+    Fi[j] = bi[j][0] / (float) Fivec[0].Nelements;
   }