added combined model

c_src/_web.c

@@ -2,6 +2,7 @@
 #include <numpy/arrayobject.h>
 #include "heron.h"
 #include "segment.h"
+#include "spline.h"
 #include "areas.h"
 #include "intersection.h"
 #include "generate.h"
@@ -270,12 +271,12 @@ static PyObject *web_line_intersect(PyObject *self, PyObject *args)
 
 static PyObject *web_generate_planet(PyObject *self, PyObject *args)
 {
-    int n_layers,n_1,n_2,bright_type;
-    double lambda0,phi0,p_u1,p_u2;
-    PyObject *bright_obj;
+    int n_layers,n_1,n_2,bright_type,n_star;
+    double lambda0,phi0,p_u1,p_u2,rp;
+    PyObject *bright_obj,*teff_obj,*flux_obj;
 
     /* Parse the input tuple */
-    if (!PyArg_ParseTuple(args, "iddddiO", &n_layers,&lambda0,&phi0,&p_u1,&p_u2,&bright_type,&bright_obj))
+    if (!PyArg_ParseTuple(args, "iddddiOOOid", &n_layers,&lambda0,&phi0,&p_u1,&p_u2,&bright_type,&bright_obj,&teff_obj,&flux_obj,&n_star,&rp))
         return NULL;
 
     /* Call the external C function to compute the area. */
@@ -290,21 +291,54 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)
     /* Get pointers to the data as C-types. */
     double *brightness_params    = (double*)PyArray_DATA(bright_array);
 
+    PyObject *teff_array = PyArray_FROM_OTF(teff_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (teff_array == NULL) {
+        Py_XDECREF(teff_array);
+        return NULL;
+    }
+    double *stellar_teffs    = (double*)PyArray_DATA(teff_array);
+
+    PyObject *flux_array = PyArray_FROM_OTF(flux_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (flux_array == NULL) {
+        Py_XDECREF(flux_array);
+        return NULL;
+    }
+    double *stellar_fluxes    = (double*)PyArray_DATA(flux_array);
+
+
     /* NEED TO GENERATE THE GRID HERE */
     double **bb_g;
+    double *ypp;
 
-    double T_start =500;
+    double T_start =0;
     double T_end =10000;
-    int n_temps=32;
-    int n_bb_seg=10;
+    int n_temps=100;
+    int n_bb_seg=20;
+
+    double star_bright = 1.0;
+
+    double r2 = 1.0/rp; //invert planet radius ratio - planets always have radius 1 in this code
 
-    if(bright_type == 1 || bright_type == 3 || bright_type == 4 || bright_type == 8){
+    if(bright_type == 1 || bright_type == 3 || bright_type == 4 || bright_type == 8|| bright_type == 10){
         double l1 = brightness_params[1];
         double l2 = brightness_params[2];
+        double star_T =brightness_params[0];
+        double ypval;
+        double yppval;
+
+//        star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
+        ypp = spline_cubic_set( n_star, stellar_teffs, stellar_fluxes, 0, 0, 0, 0 );
+
+        star_bright = spline_cubic_val( n_star, stellar_teffs, stellar_fluxes, ypp, star_T, &ypval, &yppval);
+        free(ypp);
+
+        star_bright = star_bright*M_PI*pow(r2,2);
+
+
         bb_g = bb_grid(l1, l2, T_start, T_end,n_temps,n_bb_seg);
     }
 
-    map_model(planet_struct,n_layers,lambda0,phi0,p_u1,p_u2,bright_type,brightness_params,bb_g);
+    map_model(planet_struct,n_layers,lambda0,phi0,p_u1,p_u2,bright_type,brightness_params,bb_g,star_bright);
 
     /* Build the output tuple */
 
@@ -476,7 +510,7 @@ static PyObject *web_lightcurve(PyObject *self, PyObject *args)
     double *star_teff    = (double*)PyArray_DATA(teff_array);
 
     PyObject *flux_array = PyArray_FROM_OTF(flux_obj, NPY_DOUBLE, NPY_IN_ARRAY);
-    if (t_array == NULL) {
+    if (flux_array == NULL) {
         Py_XDECREF(flux_array);
         return NULL;
     }
@@ -542,10 +576,10 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)
     /* NEED TO GENERATE THE GRID HERE */
     double **bb_g;
 
-    double T_start =500;
-    double T_end =10000;
-    int n_temps=32;
-    int n_bb_seg=10;
+    double T_start =0;
+    double T_end =1000;
+    int n_temps=100;
+    int n_bb_seg=20;
 
     if(bright_type == 1 || bright_type == 3 || bright_type == 4){
         double l1 = brightness_params[1];
@@ -556,7 +590,9 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)
     double lambda0 = 0;
     double phi0 = 0;
 
-    double *vals = call_map_model(la,lo,lambda0,phi0,bright_type,brightness_params,bb_g,0,0.0,0.0,0.0,0.0);
+    double star_bright = 1.0;
+
+    double *vals = call_map_model(la,lo,lambda0,phi0,bright_type,brightness_params,bb_g,0,0.0,0.0,0.0,0.0,star_bright);
 
     /* Build the output tuple */
 

c_src/brightness_maps.c

@@ -239,3 +239,18 @@ double kreidberg_2016(double lat, double lon, double insol, double albedo, doubl
     }
     return T;
 }
+
+double lambertian(double lat, double lon, double insol, double albedo){
+    double b;
+
+//    printf("%f %f \n",insol,albedo);
+
+    if((-M_PI/2.0 <= lon) && (lon <= M_PI/2.0)){
+        b = albedo*insol*cos(lat)*cos(lon)/M_PI;
+    }
+    else{
+        b = 0.0;
+    }
+
+    return b;
+}

c_src/brightness_maps.h

@@ -7,4 +7,5 @@ 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 spherical(double lat, double lon, double *a);
 double kreidberg_2016(double lat, double lon, double insol, double albedo, double redist);
-double great_circle(double la0,double lo0,double lambda0,double phi0,double r,double theta);
+double lambertian(double lat, double lon, double insol, double albedo);
+double great_circle(double la0,double lo0,double lambda0,double phi0,double r,double theta);

c_src/generate.c

@@ -7,7 +7,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 
-void map_model(double **planet,int n_layers,double lambda0, double phi0, double u1, double u2,int brightness_model,double *brightness_params,double **bb_g){
+void map_model(double **planet,int n_layers,double lambda0, double phi0, double u1, double u2,int brightness_model,double *brightness_params,double **bb_g,double star_bright){
     double point_T,mu,p_t_bright;
     double l1,l2,R_mid,theta_mid,mid_x,mid_y;
     double la, lo;
@@ -18,7 +18,7 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
 
 //    printf("bb_g 4 %f\n",bb_g[0][1]);
 
-    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6|| brightness_model == 8) {
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6|| brightness_model == 8|| brightness_model == 10) {
         l1 = brightness_params[1];
         l2 = brightness_params[2];
     }
@@ -49,7 +49,7 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
 
         free(coords);
 
-        double *vals = call_map_model(la,lo,lambda0,phi0,brightness_model,brightness_params,bb_g,make_grid,planet[k][10],planet[k][11],planet[k][13],planet[k][14]);
+        double *vals = call_map_model(la,lo,lambda0,phi0,brightness_model,brightness_params,bb_g,make_grid,planet[k][10],planet[k][11],planet[k][13],planet[k][14],star_bright);
 
         planet[k][16] = vals[0]; 
         planet[k][17] = vals[1];
@@ -65,7 +65,7 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
 
 }
 
-double *call_map_model(double la,double lo,double lambda0, double phi0,int brightness_model,double *brightness_params,double **bb_g,int make_grid,double theta1, double theta2, double r1, double r2){
+double *call_map_model(double la,double lo,double lambda0, double phi0,int brightness_model,double *brightness_params,double **bb_g,int make_grid,double theta1, double theta2, double r1, double r2, double star_bright){
     double point_T,mu,p_t_bright,point_b;
     double l1,l2;
     double *output;
@@ -75,7 +75,7 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
 
 //    printf("bb_g %f\n",bb_g[0][1]);
 
-    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6|| brightness_model == 8) {
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6|| brightness_model == 8|| brightness_model == 10) {
         l1 = brightness_params[1];
         l2 = brightness_params[2];
     }
@@ -137,6 +137,25 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
         point_T = Hotspot_b(la, lo, la0,lo0,p_T,spot_T,size,make_grid ,theta1,theta2,r1,r2,lambda0,phi0);
         point_b = Hotspot_b(la, lo, la0,lo0,b1,b2,size,make_grid ,theta1,theta2,r1,r2,lambda0,phi0);
     }
+    if(brightness_model == 9){
+        double albedo = brightness_params[0];
+        double ars = pow(brightness_params[1],2);
+        double insol = ars*star_bright;
+        point_b = lambertian(la,lo,insol,albedo);
+    }
+    if(brightness_model == 10){
+        double xi =brightness_params[3];
+        double T_n =brightness_params[4];
+        double delta_T =brightness_params[5];
+        point_T = zhang_2016(la,lo,xi,T_n,delta_T);
+        point_b = bb_interp(point_T, bb_g);
+
+        double albedo = brightness_params[6];
+        double ars = pow(brightness_params[7],2);
+        double insol = ars*star_bright;
+//        printf("%f %f %f\n",ars, star_bright,insol);
+        point_b = point_b + lambertian(la,lo,insol,albedo);
+    }
 
     output = malloc(sizeof(double) * 2); // dynamic `array (size 2) of pointers to double`
 

c_src/generate.h

@@ -1,3 +1,3 @@
 double **generate_planet(int n_layers);
-void map_model(double **planet,int n_layers,double lambda0, double phi0, double u1, double u2,int brightness_model,double *brightness_params,double **bb_g);
-double *call_map_model(double la,double lo,double lambda0, double phi0,int brightness_model,double *brightness_params,double **bb_g,int make_grid, double theta1, double theta2, double r1, double r2);
+void map_model(double **planet,int n_layers,double lambda0, double phi0, double u1, double u2,int brightness_model,double *brightness_params,double **bb_g, double star_bright);
+double *call_map_model(double la,double lo,double lambda0, double phi0,int brightness_model,double *brightness_params,double **bb_g,int make_grid, double theta1, double theta2, double r1, double r2, double star_bright);

c_src/web.c

@@ -39,22 +39,24 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
     star_bright = 1.0;
 
     // brightness model 1 is the Xi 2016 model, requires a stellar temperature
-    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6 || brightness_model == 8){
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6 || brightness_model == 8|| brightness_model == 10){
         double l1 = brightness_params[1];
         double l2 = brightness_params[2];
         double star_T =brightness_params[0];
         double ypval;
         double yppval;
 
 //        star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
+
         ypp = spline_cubic_set( nstars, stellar_teffs, stellar_fluxes, 0, 0, 0, 0 );
         star_bright = spline_cubic_val( nstars, stellar_teffs, stellar_fluxes, ypp, star_T, &ypval, &yppval);
         free(ypp);
 
+
         star_bright = star_bright*M_PI*pow(r2,2);
 
     // also requires the precomputation of the blackbody interpolation grid
-        printf("%f %f %f %f %f %f\n",l1,l2,T_start,T_end,n_temps,n_bb_seg);
+//        printf("%f %f %f %f %f %f\n",l1,l2,T_start,T_end,n_temps,n_bb_seg);
         bb_g = bb_grid(l1, l2, T_start, T_end,n_temps,n_bb_seg);
 //        printf("bb_g init %f\n",bb_g[0][1]);
 
@@ -87,7 +89,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
 
 //        printf("bb_g 3 %f\n",bb_g[0][1]);
 
-        map_model(planet,n_layers,lambda0,phi0,u1,u2,brightness_model,brightness_params,bb_g);
+        map_model(planet,n_layers,lambda0,phi0,u1,u2,brightness_model,brightness_params,bb_g,star_bright);
 
         p_bright = 0.0;
         for (j = 0; j < pow(n_layers,2); j++) {
@@ -114,7 +116,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
     free(coords);
     free(transit_coords);
 
-    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4  || brightness_model == 8){
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4  || brightness_model == 8|| brightness_model == 10){
         for (int j = 0; j < 4; ++j) {
           free(bb_g[j]);
         }

spiderman/params.py

@@ -85,6 +85,14 @@ def __init__(self,brightness_model='xi'):
 			self.thermal= True			# Is this a thermal distribution?
 			self.grid_size = 10
 
+		elif brightness_model == 'lambertian':
+			self.brightness_type= 9 # Integer model identifer
+			self.thermal= False			# Is this a thermal distribution?
+
+		elif brightness_model == 'combine':
+			self.brightness_type= 10 # Integer model identifer
+			self.thermal= True			# Is this a thermal distribution?
+
 		else:
 			print('Brightness model "'+str(brightness_model)+'" not recognised!')
 			quit()
@@ -177,6 +185,26 @@ def format_bright_params(self):
 				print('should be',brightness_param_names)
 				quit()
 
+		elif (self.brightness_type == 9):
+			brightness_param_names = ['albedo']
+			ars = 1.0/self.a
+			try:
+				brightness_params = [self.albedo,ars]
+			except:
+				print('Brightness parameters incorrectly assigned')
+				print('should be',brightness_param_names)
+				quit()
+		elif (self.brightness_type == 10):
+			brightness_param_names = ['T_s','l1','l2','xi','T_n','delta_T','albedo']
+			ars = 1.0/self.a
+			try:
+				brightness_params = [self.T_s,self.l1,self.l2,self.xi,self.T_n,self.delta_T,self.albedo,ars]
+			except:
+				print('Brightness parameters incorrectly assigned')
+				print('should be',brightness_param_names)
+				quit()
+
+
 		elif any(b == None for b in brightness_params):
 			print('Brightness parameters incorrectly assigned')
 			print('should be',brightness_param_names)

spiderman/plot.py

@@ -436,6 +436,8 @@ def plot_planet(spider_params,t,ax=False,min_temp=False,max_temp=False,temp_map=
 		min_temp = np.min(temps)
 		max_temp = np.max(temps)
 
+	temps = planet[:,b_i]
+	print('min ',np.min(temps), 'max',np.max(temps))
 
 	dp = ((max_temp-min_temp)*min_bright)
 

spiderman/web.py

@@ -23,12 +23,25 @@ def circle_intersect(x1,y1,r1,x2,y2,r2):
 def line_intersect(x1,y1,x2,y2,r2):
 	return _web.line_intersect(x1,y1,x2,y2,r2)
 
-def generate_planet(spider_params,t,use_phase=False):
+def generate_planet(spider_params,t,use_phase=False,stellar_grid=False):
 	if use_phase == True:
 		t = spider_params.t0 + spider_params.per*t
 	brightness_params = spider_params.format_bright_params()
+
+	if spider_params.thermal == True:
+		if stellar_grid == False:
+			star_grid = sp.stellar_grid.gen_grid(spider_params.l1,spider_params.l2)
+			teffs = star_grid[0]
+			totals = star_grid[1]
+		else:
+			teffs = stellar_grid[0]
+			totals = stellar_grid[1]
+	else:
+		teffs = []
+		totals = []
+
 	spider_params.calc_substellar(t)
-	return np.array(_web.generate_planet(spider_params.n_layers,spider_params.lambda0,spider_params.phi0,spider_params.p_u1,spider_params.p_u2,spider_params.brightness_type,brightness_params))
+	return np.array(_web.generate_planet(spider_params.n_layers,spider_params.lambda0,spider_params.phi0,spider_params.p_u1,spider_params.p_u2,spider_params.brightness_type,brightness_params,teffs,totals,len(totals),spider_params.rp))
 
 def call_map_model(spider_params,la,lo):
 	brightness_params = spider_params.format_bright_params()