Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/Makefile
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/Makefile	(revision 37537)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/Makefile	(revision 37538)
@@ -35,6 +35,5 @@
 $(SRC)/time.$(ARCH).o		 \
 $(SRC)/gaussj.$(ARCH).o		 \
-$(SRC)/legendre.$(ARCH).o	 \
-$(SRC)/gsl_utils.$(ARCH).o	 \
+$(SRC)/spherical.$(ARCH).o       \
 $(SRC)/vstats.$(ARCH).o		 \
 $(SRC)/config.$(ARCH).o		 \
@@ -48,4 +47,7 @@
 $(SRC)/gprint.$(ARCH).o	 \
 $(SRC)/version.$(ARCH).o
+
+# $(SRC)/legendre.$(ARCH).o	 \
+# $(SRC)/gsl_utils.$(ARCH).o	 \
 
 TYPETEST = \
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h	(revision 37537)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h	(revision 37538)
@@ -447,3 +447,15 @@
 double legendre_Plm_sphere(int l, int m, double x, double *err);
 
-# endif
+typedef struct {
+  double *dR_B;
+  double *dR_E;
+  double *dD_B;
+  double *dD_E;
+  int lmax;
+  int Nterms;
+} VSHterms;
+
+VSHterms *InitVSHterms (int lmax);
+void getVSHterms (VSHterms *terms, double R, double D);
+
+# endif
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c	(revision 37537)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c	(revision 37538)
@@ -1,24 +1,26 @@
 # include "ohana.h"
-
-typedef struct {
-  double *dR_B;
-  double *dR_E;
-  double *dD_B;
-  double *dD_E;
-  int lmax;
-  int Nterms;
-} VSHterms;
+# include <gsl/gsl_sf_legendre.h>
 
 VSHterms *InitVSHterms (int lmax) {
 
   // 0,0  -- N = 1, S = 1
-  // 1,-1; 1,0; 1,1 -- N = 3, S = 4
-  // 2,-2; 2,-1; 2,0; 2,1; 2,2 -- N = 5, S = 9
-  // 3,-3; 3,-2; 3,-1; 3,0; 3,1; 3,2, 3,3 -- N = 7, S = 16 
-  int 
+  // 1,-1; 1,0; 1,1 -- N = 3, S = 4 (l+1)^2
+  // 2,-2; 2,-1; 2,0; 2,1; 2,2 -- N = 5, S = 9 (l+1)^2
+  // 3,-3; 3,-2; 3,-1; 3,0; 3,1; 3,2, 3,3 -- N = 7, S = 16 (l+1)^2
 
+  VSHterms *terms = NULL;
+  ALLOCATE (terms, VSHterms, 1);
+  terms->lmax = lmax;
+  terms->Nterms = SQ(lmax+1);
+
+  ALLOCATE (terms->dR_B, double, terms->Nterms);
+  ALLOCATE (terms->dR_E, double, terms->Nterms);
+  ALLOCATE (terms->dD_B, double, terms->Nterms);
+  ALLOCATE (terms->dD_E, double, terms->Nterms);
+  
+  return terms;
 }
 
-void getVSHterms (int lmax, double R, double D, VSHterms *VlmSet) {
+void getVSHterms (VSHterms *terms, double R, double D) {
 
   // given a point on the sky R,D
@@ -82,6 +84,14 @@
 
   int l, m;
+
+  // Vlm(l=0) = 0, so assign directly (we ignore for fitting)
   int n = 0;
-  for (l = 1; l < lmax; l++) {
+  terms->dR_B[n] = 0.0;
+  terms->dD_B[n] = 0.0;
+  terms->dR_E[n] = 0.0;
+  terms->dD_E[n] = 0.0;
+  n++;
+
+  for (l = 1; l < terms->lmax; l++) {
 
     for (m = -l; m <= l; m++) {
@@ -101,8 +111,8 @@
       double S_D_lm = (m < 0) ? dPlm_dD * sin(m_abs * Rrad) : dPlm_dD * cos(m_abs * Rrad);
   
-      VlmSet->dR_B[n] = +S_D_lm;
-      VlmSet->dD_B[n] = -S_R_lm / cos(Drad);
-      VlmSet->dR_E[n] = +S_R_lm / cos(Drad);
-      VlmSet->dR_B[n] = +S_D_lm;
+      terms->dR_B[n] = +S_D_lm;
+      terms->dD_B[n] = -S_R_lm / cos(Drad);
+      terms->dR_E[n] = +S_R_lm / cos(Drad);
+      terms->dD_E[n] = +S_D_lm;
       n++;
     }