added spherical harmonics model

c_src/_web.c

@@ -10,6 +10,7 @@
 #include "web.h"
 #include "pyutil.h"
 #include "blackbody.h"
+#include "brightness_maps.h"
 #include <stdio.h>
 
 static char module_docstring[] =

c_src/brightness_maps.c

@@ -1,4 +1,8 @@
 #include "math.h"
+#include "harmonics.h"
+#include "legendre_polynomial.h"
+#include <stdlib.h>
+#include <stdio.h>
 
 double Uniform_b(double p_bright){
     return p_bright/M_PI;
@@ -58,4 +62,35 @@ double zhang_2016(double lat, double lon, double xi, double T_n, double delta_T)
     }
 
     return T;
+}
+
+double spherical(double lat, double lon, double *a){
+    double x_vec[1];
+    double fx2;
+    double *fx2_vec;
+    double theta = (M_PI/2.0) - lat;
+    double phi = M_PI + lon;
+
+    int orders = a[0];
+
+    int k = 1;
+
+    x_vec[0] = cos(theta);
+
+    double val = 0.0;
+    for (int l = 0; l < (orders); ++l) {
+
+      for (int m = 0; m < (l+1); ++m) {
+        fx2_vec = pm_polynomial_value(1,l,m,x_vec);
+        fx2 = fx2_vec[l];
+        free ( fx2_vec );
+
+        val = val + a[k]*cos(m*phi)*fx2;
+
+        k = k +1;
+      }
+    }
+
+    return val;
+
 }

c_src/brightness_maps.h

@@ -2,4 +2,5 @@ double Uniform_b(double p_bright);
 double Uniform_T(double T_bright);
 double Two_b(double la, double lo, double p_day, double p_night);
 double Two_T(double la, double lo, double p_day, double p_night);
-double zhang_2016(double lat, double lon, double zeta, double T_n, double delta_T);
+double zhang_2016(double lat, double lon, double zeta, double T_n, double delta_T);
+double spherical(double lat, double lon, double *a);

c_src/generate.c

@@ -59,6 +59,13 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
         planet[0][16] = bb_flux(l1,l2,point_T,n_bb_seg);
     }
 
+    if(brightness_model == 5){
+        double p_t_bright = spherical(la,lo,brightness_params);
+        planet[0][16] = p_t_bright;
+        planet[0][17] = 0.0;
+    }
+
+
     for (int k = 1; k < pow(n_layers,2); ++k) {
 
         if(k == 0){
@@ -115,6 +122,12 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
             planet[k][17] = point_T;
             planet[k][16] = bb_flux(l1,l2,point_T,n_bb_seg);
         }
+        if(brightness_model == 5){
+            double p_t_bright = spherical(la,lo,brightness_params);
+            planet[k][16] = p_t_bright;
+            planet[k][17] = 0.0;
+        }
+
         // limb darkening
 
         mu = sqrt(1 - pow(R_mid,2));