Merge pull request #4 from tomlouden/interpolated_temperatures

Interpolated temperatures

c_src/_web.c

@@ -436,12 +436,12 @@ static PyObject *web_calc_substellar(PyObject *self, PyObject *args)
 
 static PyObject *web_lightcurve(PyObject *self, PyObject *args)
 {
-    int n_layers, bright_type;
+    int n_layers, bright_type, n_star;
     double tc,per,a,inc,ecc,omega,a_rs,rp,p_u1,p_u2;
-    PyObject *t_obj,*bright_obj;
+    PyObject *t_obj,*bright_obj,*teff_obj,*flux_obj;
 
     /* Parse the input tuple */
-    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))
+    if (!PyArg_ParseTuple(args, "iOddddddddddiOOOi", &n_layers,&t_obj,&tc,&per,&a,&inc,&ecc,&omega,&a_rs,&rp,&p_u1,&p_u2,&bright_type,&bright_obj,&teff_obj,&flux_obj,&n_star))
         return NULL;
 
     PyObject *bright_array = PyArray_FROM_OTF(bright_obj, NPY_DOUBLE, NPY_IN_ARRAY);
@@ -453,8 +453,6 @@ static PyObject *web_lightcurve(PyObject *self, PyObject *args)
     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);
         return NULL;
@@ -465,13 +463,30 @@ static PyObject *web_lightcurve(PyObject *self, PyObject *args)
     /* Get pointers to the data as C-types. */
     double *t2    = (double*)PyArray_DATA(t_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 *star_teff    = (double*)PyArray_DATA(teff_array);
+
+    PyObject *flux_array = PyArray_FROM_OTF(flux_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (t_array == NULL) {
+        Py_XDECREF(flux_array);
+        return NULL;
+    }
+    double *star_flux    = (double*)PyArray_DATA(flux_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,p_u1,p_u2,bright_type,brightness_params);
+
+    double *output = lightcurve(n_layers,N,t2,tc,per,a,inc,ecc,omega,a_rs,rp,p_u1,p_u2,bright_type,brightness_params,star_teff,star_flux,n_star);
 
     PyObject *pylist = Convert_Big_Array(output,N);
 
     /* Clean up. */
     Py_DECREF(t_array);
+    Py_DECREF(teff_array);
+    Py_DECREF(flux_array);
 
     return pylist;
 }

c_src/brightness_maps.c

@@ -67,12 +67,12 @@ 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;
+    double theta = (M_PI/2.0) - (lat+a[1]);
+    double phi = M_PI + (lon+a[2]);
 
     int orders = a[0];
 
-    int k = 1;
+    int k = 3;
 
     x_vec[0] = cos(theta);
 
@@ -82,14 +82,14 @@ double spherical(double lat, double lon, double *a){
       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 );
+        free(fx2_vec);
 
         val = val + a[k]*cos(m*phi)*fx2;
 
         k = k +1;
       }
     }
 
-    return val;
+    return pow(val,2);
 
 }

c_src/generate.c

@@ -20,7 +20,7 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
         l2 = brightness_params[2];
     }
 
-    for (int k = 1; k < pow(n_layers,2); ++k) {
+    for (int k = 0; k < pow(n_layers,2); ++k) {
 
         if(k == 0){
             R_mid = 0;
@@ -78,6 +78,7 @@ 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);
             planet[k][16] = bb_interp(point_T, bb_g);
+//            printf("%f %f %f\n",point_T,bb_flux(l1,l2,point_T,n_bb_seg),planet[k][16]);
         }
         if(brightness_model == 5){
             double p_t_bright = spherical(la,lo,brightness_params);

c_src/orthographic.c

@@ -36,6 +36,8 @@ double *cart_to_ortho(double R, double x, double y, double lambda0, double phi0)
 
 	lambda = lambda0 + atan2(x*sin(c),(rho*cos(c)*cos(phi0) - y*sin(c)*sin(phi0) ) );
 
+	lambda = lambda - M_PI*2*floor(lambda/(M_PI*2));
+
 	if(lambda>M_PI){
 		lambda = lambda-2*M_PI;
 	}

c_src/web.c

@@ -8,22 +8,23 @@
 #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 u1, double u2,int brightness_model,double *brightness_params){
+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,double *stellar_teffs,double *stellar_fluxes,int nstars){
     int n,j;
     double phase,lambda0,phi0;
     double *coords;
     double p_blocked, p_bright,phase_z,phase_dz,phase_dt;
     double star_bright;
     double **bb_g;
+    double *ypp;
 
     double r2 = 1.0/rp; //invert planet radius ratio - planets always have radius 1 in this code
 
     double c = 299792458.0;
 
     int n_bb_seg = 20;
-    double T_start =500;
+    double T_start =0;
     double T_end =10000;
-    int n_temps=32;
+    int n_temps=100;
 
     double *output = malloc(sizeof(double) * n_points);
 
@@ -35,16 +36,20 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
 
     double transit_z = transit_coords[3];
 
-    if(brightness_model == 0){
-        star_bright = 1.0;
-    }
+    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){
         double l1 = brightness_params[1];
         double l2 = brightness_params[2];
         double star_T =brightness_params[0];
-        star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
+        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);
+
         star_bright = star_bright*M_PI*pow(r2,2);
 
     // also requires the precomputation of the blackbody interpolation grid
@@ -108,10 +113,8 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
         for (int j = 0; j < 4; ++j) {
           free(bb_g[j]);
         }
+        free(bb_g);
     }
-    free(bb_g);
-
-
 
     return output;
 }

c_src/web.h

@@ -1,4 +1,4 @@
-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 *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* stellar_teff,double* stellar_flux,int nstars);
 double *call_blocked(int n_layers, int n_points, double *x2, double *y2, double r2);
 double calc_phase(double t, double t0, double per);
 double *calc_substellar(double phase, double *coords);

spiderman/__init__.py

@@ -1,6 +1,28 @@
+from os.path import expanduser,join
+home = expanduser("~")
+param_file = join(home,'.spidermanrc')
+
+class RcParams(object):
+	def __init__(self):
+		self.dict = {}
+		try:
+			print 'looking for spidermanrc file at '+param_file
+			for line in open(param_file):
+				key = line.split(':')[0].replace(' ','')
+				cval = line.split(':')[1].replace(' ','')
+				self.dict[key] = cval
+			RcParams.read = True
+		except:
+			RcParams.read = False
+			print 'no spidermanrc file detected'
+	def __getitem__(self, key):
+	    return self.dict[key]
+rcParams = RcParams()
+
 from spiderman.params import *
 from spiderman.web import *
 from spiderman.plot import *
 from spiderman.test import *
+from spiderman.stellar_grid import *
 
-__all__ = ["web","params","_web","plot","test"]
+__all__ = ["web","params","_web","plot","test","stellar_grid"]

spiderman/params.py

@@ -26,13 +26,15 @@ def __init__(self,brightness_model='xi'):
 
 			self.brightness_type= 0	# Integer model identifier
 			self.pb= None			# Relative planet brightness (Star is 1)
+			self.thermal= False			# Is this a thermal distribution?
 
 		elif brightness_model == 'uniform temperature':
 			self.n_layers = 1		# The default resolution for the grid
 
 			self.brightness_type= 1	# Integer model identifier
 			self.T_p= None			# Relative planet brightness (Star is 1)
 			self.T_s= None			# **STELLAR** effective temperature
+			self.thermal= True			# Is this a thermal distribution?
 
 		elif brightness_model == 'two brightness':
 			self.brightness_type= 2	# Integer model identifier
@@ -44,17 +46,20 @@ def __init__(self,brightness_model='xi'):
 			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
+			self.thermal= True			# Is this a thermal distribution?
 
 		elif brightness_model == 'zhang':
 			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
+			self.thermal= True			# Is this a thermal distribution?
 
 		elif brightness_model == 'spherical':
 			self.brightness_type= 5	# Integer model identifier
 			self.a= None			# Ratio between radiative and advective timescales
+			self.thermal= False			# Is this a thermal distribution?
 
 		else:
 			print('Brightness model "'+str(brightness_model)+'" not recognised!')
@@ -102,9 +107,9 @@ def format_bright_params(self):
 				print('should be',brightness_param_names)
 				quit()
 		if (self.brightness_type == 5):
-			brightness_param_names = ['orders','sph']
+			brightness_param_names = ['orders','sph','la_o','lo_o']
 			try:
-				brightness_params = [self.orders] + self.sph
+				brightness_params = [self.orders,self.la_o,self.lo_o] + self.sph
 			except:
 				print('Brightness parameters incorrectly assigned')
 				print('should be',brightness_param_names)
@@ -216,9 +221,79 @@ def plot_quad(self,min_temp=False,max_temp=False,temp_map=False,min_bright=0.2,s
 
 		return fig
 
-	def lightcurve(self,t,use_phase=False):
+	def plot_uncertainty(self,fs,min_temp=False,max_temp=False,temp_map=True,min_bright=0.2,scale_planet=1.0,planet_cen=[0.0,0.0],use_phase=False,show_cax=True,mycmap=plt.cm.viridis_r,theme='black'):
+
+		if theme == 'black':
+			bg = 'black'
+			tc = ("#04d9ff")
+		else:
+			bg = 'white'
+			tc = 'black'
+
+		fig, axs = plt.subplots(2,2,figsize=(8/0.865,8),facecolor=bg)
+
+		# need a "get max" or "get min" function or something similar so the scale can be set properly
+
+		fs = np.array(fs)
+
+		min_temp = np.min(fs)
+		max_temp = np.max(fs)
+
+		dp = ((max_temp-min_temp)*min_bright)
+
+		sp.plot_dist(self,fs[0],ax=axs[0,0],show_cax=False,min_temp=min_temp,max_temp=max_temp,temp_map=temp_map,mycmap=mycmap,theme=theme,min_bright=min_bright)
+		sp.plot_dist(self,fs[1],ax=axs[0,1],show_cax=False,min_temp=min_temp,max_temp=max_temp,temp_map=temp_map,mycmap=mycmap,theme=theme,min_bright=min_bright)
+		sp.plot_dist(self,fs[2],ax=axs[1,0],show_cax=False,min_temp=min_temp,max_temp=max_temp,temp_map=temp_map,mycmap=mycmap,theme=theme,min_bright=min_bright)
+		sp.plot_dist(self,fs[3],ax=axs[1,1],show_cax=False,min_temp=min_temp,max_temp=max_temp,temp_map=temp_map,mycmap=mycmap,theme=theme,min_bright=min_bright)
+
+#		divider = make_axes_locatable(fig)
+
+#		zero_temp = min_temp - dp
+		zero_temp = min_temp
+
+		data = [np.linspace(zero_temp,max_temp,1000)]*2
+		fake, fake_ax = plt.subplots()
+		mycax = fake_ax.imshow(data, interpolation='none', cmap=mycmap)
+		plt.close(fake)
+
+		fig.subplots_adjust(right=0.80)
+		cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
+
+		if show_cax == True:
+#			cax = divider.append_axes("right", size="20%", pad=0.05)
+		#	cbar = plt.colorbar(mycax, cax=cax,ticks=[1100,1300,1500,1700,1900])
+			cbar = plt.colorbar(mycax, cax=cbar_ax)
+			cbar.ax.tick_params(colors=tc)
+
+			cbar.ax.invert_yaxis()
+
+			if temp_map == True:
+				cbar.set_label(r'Temperature precision (\%)',color=tc)  # horizontal colorbar
+			else:
+				cbar.set_label(r'Flux precision (\%)',color=tc)  # horizontal colorbar
+
+		fig.subplots_adjust(wspace=0, hspace=0)
+
+		return fig
+
+	def lightcurve(self,t,use_phase=False,stellar_grid=False):
+
 		brightness_params = self.format_bright_params()
+
+		if self.thermal == True:
+			if stellar_grid == False:
+				star_grid = sp.stellar_grid.gen_grid(self.l1,self.l2)
+				teffs = star_grid[0]
+				totals = star_grid[1]
+			else:
+				teffs = stellar_grid[0]
+				totals = stellar_grid[1]
+		else:
+			teffs = []
+			totals = []
+
 		if use_phase == True:
 			t = self.t0 + self.per*t
-		out = _web.lightcurve(self.n_layers,t,self.t0,self.per,self.a_abs,self.inc,self.ecc,self.w,self.a,self.rp,self.p_u1,self.p_u2,self.brightness_type,brightness_params)
+
+		out = _web.lightcurve(self.n_layers,t,self.t0,self.per,self.a_abs,self.inc,self.ecc,self.w,self.a,self.rp,self.p_u1,self.p_u2,self.brightness_type,brightness_params,teffs,totals,len(totals))
 		return np.array(out)