Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c	(revision 37536)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c	(revision 37537)
@@ -1,5 +1,24 @@
 # include "ohana.h"
 
-void VSH () {
+typedef struct {
+  double *dR_B;
+  double *dR_E;
+  double *dD_B;
+  double *dD_E;
+  int lmax;
+  int Nterms;
+} VSHterms;
+
+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 
+
+}
+
+void getVSHterms (int lmax, double R, double D, VSHterms *VlmSet) {
 
   // given a point on the sky R,D
@@ -15,9 +34,12 @@
   // dD^E = \sum_{l,m} e_{l,m} dD^E_{l,m} (l,m,R,D)
   
+  // dR^E_{l,m} is Valeri's Elm, item[0], dD^E_{l,m} is Valeri's Elm, item[1]
+  // dR^B_{l,m} is Valeri's Blm, item[0], dD^B_{l,m} is Valeri's Blm, item[1]
+
   // NOTE : as implemented by Valeri, the coeffs b_{l,m} are the same for both R and D
   // components; ditto for e_{l,m}
 
-  // there are thus for X_l,m terms that need to be generated, for dR^B, dR^E, dD^B, dD^E.  Valeri calls
-  // these S_lm_lat and S_lm_lon
+  // there are thus 4 X_l,m terms that need to be generated, for dR^B, dR^E, dD^B, dD^E.
+  // Valeri calls these S_lm_lat and S_lm_lon
 
   // dR^B_{l,m} = +S^D_{l,m}(l,m,R,D)
@@ -27,11 +49,63 @@
   // dD^E_{l,m} = +S^D_{l,m}(l,m,R,D)
 
-  // S^R_{l,m} : 
+  // S^R_{l,m} == Slm_lng
   // (m >  0) : -m S_{l,m}(l,-m,R,D) 
   // (m <= 0) : -m S_{l,m}(l,-m,R,D)
   // NOTE: these are in fact the same...
+  // NOTE: sending in -m means we get sin() for (m > 0) and cos() for (m <= 0)
 
-  // S_^D_{l,m} :
-  // P_{l,m}(D) 
+  // S_{l,m}(l,m,R,D) :
+  // aP_{l,m} = gsl_Plm(l,|m|,cos(90-D))
+
+  // (m >  0) : S_{l,m} = aP_{l,m}*cos(m*R)    = gsl_Plm(l,|m|,cos(90-D))*cos(|m|*R)  [cos(x) == cos(-x)]
+  // (m <  0) : S_{l,m} = aP_{l,m}*sin(|m|*R)  = gsl_Plm(l,|m|,cos(90-D))*sin(|m|*R)
+  // (m == 0) : S_{l,m} = aP_{l,m}             = gsl_Plm(l,|m|,cos(90-D))
+
+  // S^R_{l,m} = -m S_{l,m}(l,-m,R,D)
+  // (m >  0) : S^R_{l,m} = -m aP_lm*sin(|m|*R)
+  // (m <  0) : S^R_{l,m} = -m aP_lm*cos(|m|*R)
+  // (m == 0) : S^R_{l,m} = 0
+
+  // S_^D_{l,m} == Slm_lat:
+  //
+  // aP_{l,  m}(D) = gsl_Plm(l,  |m|,cos(90-D))
+  // aP_{l+1,m}(D) = gsl_Plm(l+1,|m|,cos(90-D))
+
+  // dPlm_dD = ((l+1)*sin(D)*aP_{l,m} - (l+1-|m|)*aP_{l+1,m}) / cos(D)
+  // dSlm_dD
+  // (m == 0) : dPlm_dD
+  // (m > 0) : dPlm_dD*cos( m *R)
+  // (m < 0) : dPlm_dD*sin(|m|*R)
+
+  double Rrad = R * RAD_DEG;
+  double Drad = D * RAD_DEG;
+
+  int l, m;
+  int n = 0;
+  for (l = 1; l < lmax; l++) {
+
+    for (m = -l; m <= l; m++) {
+    
+      int m_abs = abs(m);
+      double aP_lm = gsl_sf_legendre_Plm(l, m_abs, cos(M_PI/2 - Drad)); // alt-P_lm = P_lm (l, |m|, cos(90 - D))
+
+      // note that for m == 0, S_R_lm -> 0 (due to m*)
+      double S_R_lm = (m > 0) ? -m * aP_lm * sin (m_abs * Rrad) : -m * aP_lm * cos (m_abs * Rrad) ; 
+      
+      double aP_lp1m = gsl_sf_legendre_Plm(l + 1, m_abs, cos(M_PI/2 - Drad));
+      // reuse for aP_lm above?
+      
+      double dPlm_dD = ((l+1) * sin(Drad) * aP_lm - (l+1-m_abs)*aP_lp1m) / cos(D);
+      // XXX : this diverges @ D = 90.0??
+
+      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;
+      n++;
+    }
+  }
 }
+