added model selection

c_src/_web.c

@@ -261,17 +261,26 @@ 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;
-    double xi,T_n,delta_T,lambda0,phi0,p_u1,p_u2;
+    int n_layers,n_1,n_2,bright_type;
+    double lambda0,phi0,p_u1,p_u2;
+    PyObject *bright_obj;
 
     /* Parse the input tuple */
-    if (!PyArg_ParseTuple(args, "iddddddd", &n_layers,&xi,&T_n,&delta_T,&lambda0,&phi0,&p_u1,&p_u2))
+    if (!PyArg_ParseTuple(args, "iddddiO", &n_layers,&lambda0,&phi0,&p_u1,&p_u2,&bright_type,&bright_obj))
         return NULL;
 
     /* Call the external C function to compute the area. */
     double **planet_struct = generate_planet(n_layers);
 
-    map_model(planet_struct,n_layers,xi,T_n,delta_T,lambda0,phi0,p_u1,p_u2);
+    PyObject *bright_array = PyArray_FROM_OTF(bright_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (bright_array == NULL) {
+        Py_XDECREF(bright_array);
+        return NULL;
+    }
+    /* Get pointers to the data as C-types. */
+    double *brightness_params    = (double*)PyArray_DATA(bright_array);
+
+    map_model(planet_struct,n_layers,lambda0,phi0,p_u1,p_u2,bright_type,brightness_params);
 
     /* Build the output tuple */
 
@@ -361,16 +370,25 @@ static PyObject *web_separation_of_centers(PyObject *self, PyObject *args)
 
 static PyObject *web_lightcurve(PyObject *self, PyObject *args)
 {
-    int n_layers;
-    double tc,per,a,inc,ecc,omega,a_rs,rp,xi,T_n,delta_T,p_u1,p_u2,T_s;
-    PyObject *t_obj;
+    int n_layers, bright_type;
+    double tc,per,a,inc,ecc,omega,a_rs,rp,p_u1,p_u2;
+    PyObject *t_obj,*bright_obj;
 
     /* Parse the input tuple */
-    if (!PyArg_ParseTuple(args, "iOdddddddddddddd", &n_layers,&t_obj,&tc,&per,&a,&inc,&ecc,&omega,&a_rs,&rp,&xi,&T_n,&delta_T,&p_u1,&p_u2,&T_s))
+    if (!PyArg_ParseTuple(args, "iOddddddddddiO", &n_layers,&t_obj,&tc,&per,&a,&inc,&ecc,&omega,&a_rs,&rp,&p_u1,&p_u2,&bright_type,&bright_obj))
         return NULL;
 
+    PyObject *bright_array = PyArray_FROM_OTF(bright_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (bright_array == NULL) {
+        Py_XDECREF(bright_array);
+        return NULL;
+    }
+    /* Get pointers to the data as C-types. */
+    double *brightness_params    = (double*)PyArray_DATA(bright_array);
+
     PyObject *t_array = PyArray_FROM_OTF(t_obj, NPY_DOUBLE, NPY_IN_ARRAY);
 
+
     /* If that didn't work, throw an exception. */
     if (t_array == NULL) {
         Py_XDECREF(t_array);
@@ -379,12 +397,11 @@ static PyObject *web_lightcurve(PyObject *self, PyObject *args)
 
     /* How many data points are there? */
     int N = (int)PyArray_DIM(t_array, 0);
-
     /* Get pointers to the data as C-types. */
     double *t2    = (double*)PyArray_DATA(t_array);
 
     /* Call the external C function to compute the area. */
-    double *output = lightcurve(n_layers,N,t2,tc,per,a,inc,ecc,omega,a_rs,rp,xi,T_n,delta_T,p_u1,p_u2,T_s);
+    double *output = lightcurve(n_layers,N,t2,tc,per,a,inc,ecc,omega,a_rs,rp,p_u1,p_u2,bright_type,brightness_params);
 
     PyObject *pylist = Convert_Big_Array(output,N);
 

c_src/generate.c

@@ -6,7 +6,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 
-void map_model(double **planet,int n_layers,double xi, double T_n, double delta_T,double lambda0, double phi0, double u1, double u2){
+void map_model(double **planet,int n_layers,double lambda0, double phi0, double u1, double u2,int brightness_model,double *brightness_params){
     double point_T,mu;
 
     double R = 1.0;
@@ -21,10 +21,47 @@ void map_model(double **planet,int n_layers,double xi, double T_n, double delta_
     double lo = -1*old_coords[0];
     free(old_coords);
 
-    point_T = zhang_2016(la,lo,xi,T_n,delta_T);
+    if(brightness_model == 0){
+        double p_t_bright = brightness_params[0];
+        planet[0][16] = p_t_bright/M_PI;
+        planet[0][17] = 0.0;
+    }
+    if(brightness_model == 1){
+        double point_T = brightness_params[1];
+        planet[0][17] = point_T;
+        planet[0][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+    }
+    if(brightness_model == 2){
+        double p_day = brightness_params[0];
+        double p_night = brightness_params[1];
+
+        double p_t_bright = p_night;
+        if((-M_PI/2.0 <= lo) && (lo <= M_PI/2.0)){
+            p_t_bright = p_day;
+        }
+        planet[0][16] = p_t_bright/M_PI;
+        planet[0][17] = 0.0;
+    }
+
+    if(brightness_model == 3){
+        double p_day = brightness_params[1];
+        double p_night = brightness_params[2];
 
-    planet[0][17] = point_T;
-    planet[0][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+        double point_T = p_night;
+        if((-M_PI/2.0 <= lo) && (lo <= M_PI/2.0)){
+            point_T = p_day;
+        }
+        planet[0][17] = point_T;
+        planet[0][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+    }
+    if(brightness_model == 4){
+        double xi =brightness_params[1];
+        double T_n =brightness_params[2];
+        double delta_T =brightness_params[3];
+        double point_T = zhang_2016(la,lo,xi,T_n,delta_T);
+        planet[0][17] = point_T;
+        planet[0][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+    }
 
     for (int k = 1; k < pow(n_layers,2); ++k) {
         double R_mid = (planet[k][13] + planet[k][14])/2.0;
@@ -40,10 +77,48 @@ void map_model(double **planet,int n_layers,double xi, double T_n, double delta_
 
         free(coords);
 
-        point_T = zhang_2016(la,lo,xi,T_n,delta_T);
+        if(brightness_model == 0){
+            double p_t_bright = brightness_params[0];
+            planet[k][16] = p_t_bright/M_PI;
+            planet[k][17] = 0.0;
+        }
+        if(brightness_model == 1){
+            double point_T = brightness_params[0];
+            planet[k][17] = point_T;
+            planet[k][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+        }
+        if(brightness_model == 2){
+            double p_day = brightness_params[0];
+            double p_night = brightness_params[1];
+
+            double p_t_bright = p_night;
+            if((-M_PI/2.0 <= lo) && (lo <= M_PI/2.0)){
+                p_t_bright = p_day;
+            }
 
-        planet[k][17] = point_T;
-        planet[k][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+            planet[k][16] = p_t_bright/M_PI;
+            planet[k][17] = 0.0;
+        }
+        if(brightness_model == 3){
+            double p_day = brightness_params[1];
+            double p_night = brightness_params[2];
+
+            double point_T = p_night;
+            if((-M_PI/2.0 <= lo) && (lo <= M_PI/2.0)){
+                point_T = p_day;
+            }
+            planet[k][17] = point_T;
+            planet[k][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+        }
+
+        if(brightness_model == 4){
+            double xi =brightness_params[0];
+            double T_n =brightness_params[1];
+            double delta_T =brightness_params[2];
+            double point_T = zhang_2016(la,lo,xi,T_n,delta_T);
+            planet[k][17] = point_T;
+            planet[k][16] = bb_flux(l1,l2,point_T,n_bb_seg);
+        }
         // limb darkening
 
         mu = sqrt(1 - pow(R_mid,2));

c_src/generate.h

@@ -1,3 +1,3 @@
 double **generate_planet(int n_layers);
-void map_model(double **planet,int n_layers,double xi, double T_n, double delta_T,double lambda0, double phi0, double u1, double u2);
+void map_model(double **planet,int n_layers,double lambda0, double phi0, double u1, double u2,int brightness_model,double *brightness_params);
 double zhang_2016(double lat, double lon, double zeta, double T_n, double delta_T);

c_src/web.c

@@ -6,11 +6,12 @@
 #include <stdlib.h>
 #include <stdio.h>
 
-double *lightcurve(int n_layers, int n_points, double *t, double tc, double per, double a, double inc, double ecc, double omega, double a_rs, double rp,double xi,double T_n,double delta_T,double u1, double u2,double star_T){
+double *lightcurve(int n_layers, int n_points, double *t, double tc, double per, double a, double inc, double ecc, double omega, double a_rs, double rp,double u1, double u2,int brightness_model,double *brightness_params){
     int n,j;
     double phase,lambda0,phi0;
     double *coords;
     double p_blocked, p_bright,phase_z,phase_dz,phase_dt;
+    double star_bright;
 
     double r2 = 1.0/rp; //invert planet radius ratio - planets always have radius 1 in this code
 
@@ -29,9 +30,16 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
 
     double transit_z = transit_coords[3];
 
-    double star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
+    if(brightness_model == 0){
+        star_bright = 1.0;
+    }
 
-    star_bright = star_bright*M_PI*pow(r2,2);
+    // brightness model 1 is the Xi 2016 model, requires a stellar temperature
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4){
+        double star_T =brightness_params[0];
+        star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
+        star_bright = star_bright*M_PI*pow(r2,2);
+    }
 
     free(coords);
 
@@ -67,7 +75,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
         phi0 = tan(coords[1]/coords[2]);
 
 
-        map_model(planet,n_layers,xi,T_n,delta_T,lambda0,phi0,u1,u2);
+        map_model(planet,n_layers,lambda0,phi0,u1,u2,brightness_model,brightness_params);
 
         p_bright = 0.0;
         for (j = 0; j < pow(n_layers,2); j++) {

c_src/web.h

@@ -1,2 +1,2 @@
-double *lightcurve(int n_layers, int n_points, double *t, double tc, double per, double a, double inc, double ecc, double omega, double r_s, double r2,double xi,double T_n,double delta_T,double u1, double u2,double T_s);
+double *lightcurve(int n_layers, int n_points, double *t, double tc, double per, double a, double inc, double ecc, double omega, double r_s, double r2,double u1, double u2,int brightness_model,double* brightness_params);
 double *call_blocked(int n_layers, int n_points, double *x2, double *y2, double r2);

spiderman/params.py

@@ -1,16 +1,64 @@
 class ModelParams(object):
-    def __init__(self):
-		self.t0= None
-		self.per= None
-		self.a_abs= None
-		self.inc= None
-		self.ecc= None
-		self.w= None
-		self.a= None
-		self.rp= None
-		self.xi= None
-		self.T_n= None
-		self.delta_T= None
-		self.p_u1= None
-		self.p_u2= None
-		self.T_s= None
+
+	def __init__(self,brightness_model='xi'):
+
+		self.n_layers = 5			# The default resolution for the grid
+
+		self.t0= None				# The time of central **PRIMARY** transit [jd]
+		self.per= None				# Orbital period of the planet [days]
+		self.a_abs= None			# Absolute value of the semi major axis [AU]
+		self.inc= None				# Inclination of the planetary orbit (90 is face on) [degrees]
+		self.ecc= None				# Eccentricity
+		self.w= None				# Longitude of periastron [degrees]
+		self.a= None				# Semi major axis, scaled by stellar radius [-]
+		self.rp= None				# Planet radius as a fraction of stellar radius [-]
+		self.p_u1= None				# **PLANETARY** limb darkening coefficients [-]
+		self.p_u2= None				# **PLANETARY** limb darkening coefficients [-]
+
+		if brightness_model == 'uniform brightness':
+			self.n_layers = 1		# The default resolution for the grid
+
+			self.brightness_type= 0	# Integer model identifier
+			self.pb= None			# Relative planet brightness (Star is 1)
+
+		elif brightness_model == 'uniform temperature':
+			self.n_layers = 1		# The default resolution for the grid
+
+			self.brightness_type= 1	# Integer model identifier
+			self.pb= None			# Relative planet brightness (Star is 1)
+			self.T_s= None			# **STELLAR** effective temperature
+
+		elif brightness_model == 'two brightness':
+			self.brightness_type= 2	# Integer model identifier
+			self.pb_d= None			# Relative planet brightness (Star is 1)
+			self.pb_n= None			# Relative planet brightness (Star is 1)
+
+		elif brightness_model == 'two temperature':
+			self.brightness_type= 3	# Integer model identifier
+			self.pb_d= None			# Relative planet brightness (Star is 1)
+			self.pb_n= None			# Relative planet brightness (Star is 1)
+			self.T_s= None			# **STELLAR** effective temperature
+
+		elif brightness_model == 'xi':
+			self.brightness_type= 4	# Integer model identifier
+			self.xi= None			# Ratio between radiative and advective timescales
+			self.T_n= None			# Radiative solution temperature on night side
+			self.delta_T= None		# Day/Night side difference between radiative-only temperature
+			self.T_s= None			# **STELLAR** effective temperature
+
+		else:
+			print('Brightness model "'+str(brightness_model)+'" not recognised!')
+			quit()
+
+	def format_bright_params(self):
+		if (self.brightness_type == 0):
+			brightness_params = [self.pb]
+		if (self.brightness_type == 1):
+			brightness_params = [self.T_s,self.pb]
+		if (self.brightness_type == 2):
+			brightness_params = [self.pb_d,self.pb_n]
+		if (self.brightness_type == 3):
+			brightness_params = [self.T_s,self.pb_d,self.pb_n]
+		if (self.brightness_type == 4):
+			brightness_params = [self.T_s,self.xi,self.T_n,self.delta_T]
+		return brightness_params

spiderman/web.py

@@ -33,5 +33,6 @@ def zhang_2016(lat,lon,xi,T_n,delta_T):
 def separation_of_centers(t,tc,per,a,inc,ecc,omega,a_rs,r2):
 	return _web.separation_of_centers(t,tc,per,a,inc,ecc,omega,a_rs,r2)
 
-def lightcurve(t,sp,n_layers=5):
-	return _web.lightcurve(n_layers,t,sp.t0,sp.per,sp.a_abs,sp.inc,sp.ecc,sp.w,sp.a,sp.rp,sp.xi,sp.T_n,sp.delta_T,sp.p_u1,sp.p_u2,sp.T_s)
+def lightcurve(t,sp):
+	brightness_params = sp.format_bright_params()
+	return _web.lightcurve(sp.n_layers,t,sp.t0,sp.per,sp.a_abs,sp.inc,sp.ecc,sp.w,sp.a,sp.rp,sp.p_u1,sp.p_u2,sp.brightness_type,brightness_params)