Changeset 37538

Timestamp:

Oct 30, 2014, 9:25:31 AM (12 years ago)

Author:

eugene

Message:

vsh polynomials kind of defined...

Location:

branches/eam_branches/ipp-20140904/Ohana/src/libohana

Files:

• Makefile (modified) (2 diffs)
• include/ohana.h (modified) (1 diff)
• src/spherical.c (modified) (3 diffs)

branches/eam_branches/ipp-20140904/Ohana/src/libohana/Makefile

@@ -35 +35 @@
 $(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 +47 @@
 $(SRC)/gprint.$(ARCH).o  \
 $(SRC)/version.$(ARCH).o
+
+# $(SRC)/legendre.$(ARCH).o      \
+# $(SRC)/gsl_utils.$(ARCH).o     \
 
 TYPETEST = \

branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h

@@ -447 +447 @@
 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

branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c

@@ -1 +1 @@
 # 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 +84 @@
 
   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 +111 @@
       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++;
     }