Changeset 37630

Timestamp:

Nov 18, 2014, 8:22:08 PM (12 years ago)

Author:

eugene

Message:

add spherical harmonic fitting functions

Location:

branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data

Files:

• Makefile (modified) (1 diff)
• init.c (modified) (2 diffs)
• shapply.c (added)
• shdot.c (added)
• shfit.c (modified) (3 diffs)

branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile

@@ -164 +164 @@
 $(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     \

branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c

@@ -152 +152 @@
 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 +325 @@
   {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"},

branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c

@@ -3 +3 @@
 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 +34 @@
   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 +54 @@
 
     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;
   }