Merge pull request #9 from tomlouden/speed_night_day

increased speed in models that require mesh, fixed square plots

c_src/_web.c

@@ -588,7 +588,7 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)
     int n_temps=100;
     int n_bb_seg=20;
 
-    if(bright_type == 1 || bright_type == 3 || bright_type == 4){
+    if(bright_type == 1 || bright_type == 3 || bright_type == 4 || bright_type == 6 || bright_type == 8 || bright_type == 10 || bright_type == 11){
         double l1 = brightness_params[1];
         double l2 = brightness_params[2];
         bb_g = bb_grid(l1, l2, T_start, T_end,n_temps,n_bb_seg);
@@ -601,7 +601,7 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)
 
     //NEED TO UPDATE THIS WITH CORRECT STAR BRIGHTNESS VALUES!//
 
-    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);
+    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,0.0,0.0);
 
     /* Build the output tuple */
 

c_src/brightness_maps.c

@@ -8,11 +8,14 @@
     #define M_PI 3.14159265358979323846
 #endif
 
-double Hotspot_b(double la, double lo,double la0, double lo0,double p_b,double spot_b,double size, int make_grid ,double theta1, double theta2, double r1, double r2, double lambda0, double phi0){
+double Hotspot_b(double la, double lo,double la0, double lo0,double p_b,double spot_b,double size, int make_grid ,double theta1, double theta2, double r1, double r2, double lambda0, double phi0,double la_cen,double lo_cen){
 
     double r_mid;
     double theta_mid;
 
+    la0 = la0*M_PI/180;
+    lo0 = lo0*M_PI/180;
+
     if(r1 == 0){
         r_mid = 0;
         theta_mid = 0;
@@ -25,6 +28,39 @@ double Hotspot_b(double la, double lo,double la0, double lo0,double p_b,double s
 
     if(make_grid != 0){
 
+//        printf("%f %f\n",la_cen,lo_cen);
+
+        double d1  = great_circle(la_cen,lo_cen,lambda0,phi0,r1,theta1);
+//        printf("CENTER TO INNER LINE %f\n",d1);
+        double d2  = great_circle(la_cen,lo_cen,lambda0,phi0,r2,theta1);
+//        printf("CENTER TO OUTER LINE %f\n",d2);
+//        printf("%f %f %f %f\n",la_cen,lo_cen,la0,lo0);
+        double d3 = 180.0*acos( sin(la_cen)*sin(la0) + cos(lo_cen - lo0)*cos(la_cen)*cos(la0))/M_PI;
+//        double d4 = sqrt( pow(d3-size,2) );
+
+        double r_sum1 = sqrt( pow(size + d1,2) );
+        double r_sum2 = sqrt( pow(size + d2,2) );
+
+        double r_diff1 = sqrt( pow(size - d1,2) );
+        double r_diff2 = sqrt( pow(size - d2,2) );
+
+        if((d3 < r_diff2) & (d2 < size)){
+//            printf("CIRCLE 2 INSIDE SPOT\n");
+            return spot_b;
+        }
+
+        if((d3 < r_diff1) & (size<d1)){
+//            printf("SPOT INSIDE CIRCLE 1\n");
+//            printf("%f %f %f\n",d4,d1,size);
+            return p_b;
+        }
+
+        if(d3 > r_sum2){
+//            printf("SPOT too far to intersect 1\n");
+            return p_b;
+        }
+
+
 /*        if(r1 !=0){
             double d1  = great_circle(la0,lo0,lambda0,phi0,r1,theta1);
             double d2  = great_circle(la0,lo0,lambda0,phi0,r2,theta1);
@@ -71,14 +107,35 @@ double Hotspot_b(double la, double lo,double la0, double lo0,double p_b,double s
 
         total_b = total_b / pow(grid_len,2);
 
+        if((d3 < r_diff2) & (d2 < size)){
+            if((spot_b-total_b) > 0.0){
+            printf("CIRCLE 2 INSIDE SPOT %f %f %f\n",spot_b,total_b,spot_b-total_b);
+            }
+//            return spot_b;
+        }
+
+        if((d3 < r_diff1) & (size<d1)){
+            if((p_b-total_b) > 0.0){
+            printf("SPOT INSIDE CIRCLE 1 %f\n",p_b-total_b);
+            }
+//            printf("%f %f %f\n",d4,d1,size);
+//            return p_b;
+        }
+
+        if(d3 > r_sum2){
+            if((p_b-total_b) > 0.0){
+                printf("SPOT too far to intersect 1 %f\n",p_b-total_b);
+                }
+//            return p_b;
+        }
+
+
         return total_b;
 
     }
 
-    printf("not making grid\n");
+//    printf("not making grid\n");
 
-    la0 = la0*M_PI/180;
-    lo0 = lo0*M_PI/180;
     double dist = acos( sin(la)*sin(la0) + cos(lo - lo0)*cos(la)*cos(la0));
     dist = dist*180/M_PI;
     if(dist < size){
@@ -92,15 +149,17 @@ double Hotspot_b(double la, double lo,double la0, double lo0,double p_b,double s
 
 double great_circle(double la0,double lo0,double lambda0,double phi0,double r,double theta){
 
-    la0 = la0*M_PI/180;
-    lo0 = lo0*M_PI/180;
+//    la0 = la0*M_PI/180;
+//    lo0 = lo0*M_PI/180;
 
     double x = r*cos(theta);
     double y = r*sin(theta);
+//    printf("%f %f\n",x,y);
     double *coords = cart_to_ortho(1.0, x, y, lambda0, phi0);
     double la = coords[1];
     double lo = -1*coords[0];
 
+//    printf("%f %f %f %f\n",la0,lo0,la,lo);
 
     double dist = acos( sin(la)*sin(la0) + cos(lo - lo0)*cos(la)*cos(la0));
     free(coords);

c_src/brightness_maps.h

@@ -1,4 +1,4 @@
-double Hotspot_b(double la, double lo,double la0, double lo0,double p_b,double spot_b,double size, int make_grid ,double theta1, double theta2, double r1, double r2, double lambda0, double phi0);
+double Hotspot_b(double la, double lo,double la0, double lo0,double p_b,double spot_b,double size, int make_grid ,double theta1, double theta2, double r1, double r2, double lambda0, double phi0,double la_cen,double lo_cen);
 double Hotspot_T(double la, double lo,double la0, double lo0,double p_T,double spot_T,double size, int make_grid ,double theta1, double theta2, double r1, double r2, double lambda0, double phi0);
 double Uniform_b(double p_bright);
 double Uniform_T(double T_bright);

c_src/generate.c

@@ -22,7 +22,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|| brightness_model == 10) {
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6|| brightness_model == 8|| brightness_model == 10|| brightness_model == 11) {
         l1 = brightness_params[1];
         l2 = brightness_params[2];
     }
@@ -31,6 +31,13 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
         make_grid = brightness_params[3];
     }
 
+    double *center_coords = cart_to_ortho(1.0, 0, 0, lambda0, phi0);
+
+    double la_cen = center_coords[1];
+    // sign change to longitude - we're looking at the planet from the 
+    // other side than in the simulations
+    double lo_cen = -1*center_coords[0];
+
     for (int k = 0; k < pow(n_layers,2); ++k) {
 
         if(k == 0){
@@ -53,7 +60,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],star_bright);
+        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,la_cen,lo_cen);
 
         planet[k][16] = vals[0]; 
         planet[k][17] = vals[1];
@@ -66,10 +73,11 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
 
         planet[k][16] = planet[k][16]*(1 - u1*(1-mu) - u2*(pow(1-mu,2)));
     }
+    free(center_coords);
 
 }
 
-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 *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 la_cen,double lo_cen){
     double point_T,mu,p_t_bright,point_b;
     double l1,l2;
     double *output;
@@ -92,6 +100,7 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
         point_T = Uniform_T(brightness_params[3]);
         point_b = bb_interp(point_T, bb_g);
     }
+
 //    if(brightness_model == 2){
 //        double p_day = brightness_params[0];
 //        double p_night = brightness_params[1];
@@ -129,8 +138,9 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
         double spot_b = brightness_params[4];
         double size = brightness_params[5];
 //        printf("%f %f %f %f %f\n",la0,lo0,p_b,spot_b,size);
-        point_b = Hotspot_b(la, lo, la0,lo0,p_b,spot_b,size,make_grid ,theta1,theta2,r1,r2,lambda0,phi0);
+        point_b = Hotspot_b(la, lo, la0,lo0,p_b,spot_b,size,make_grid ,theta1,theta2,r1,r2,lambda0,phi0,la_cen,lo_cen);
      }
+
     if(brightness_model == 3 || brightness_model == 8){
         double la0 = brightness_params[4];
         double lo0 = brightness_params[5];
@@ -139,8 +149,8 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
         double size = brightness_params[8];
         double b1 = bb_interp(p_T, bb_g); 
         double b2 = bb_interp(spot_T, bb_g); 
-        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);
+        point_T = Hotspot_b(la, lo, la0,lo0,p_T,spot_T,size,make_grid ,theta1,theta2,r1,r2,lambda0,phi0,la_cen,lo_cen);
+        point_b = Hotspot_b(la, lo, la0,lo0,b1,b2,size,make_grid ,theta1,theta2,r1,r2,lambda0,phi0,la_cen,lo_cen);
     }
     if(brightness_model == 9){
         double albedo = brightness_params[0];
@@ -162,6 +172,23 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
         point_b = point_b + lambertian(la,lo,insol,albedo);
     }
 
+    if(brightness_model == 11){
+        double xi =brightness_params[3];
+        double T_n =brightness_params[4];
+        double delta_T =brightness_params[5];
+        double clouds = brightness_params[6];
+
+        point_T = zhang_2016(la,lo,xi,T_n,delta_T);
+
+        point_b = bb_interp(point_T, bb_g);
+
+        if(pow(lo,2) > pow(M_PI/2.0,2)){
+            point_b = point_b - clouds;
+        }
+
+//        printf("%f %f %f\n",ars, star_bright,insol);
+    }
+
     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 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);
+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 la_cen, double lo_cen);

c_src/web.c

@@ -45,7 +45,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
     star_surface_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 || brightness_model == 6 || brightness_model == 8|| brightness_model == 10){
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6 || brightness_model == 8|| brightness_model == 10|| brightness_model == 11){
         double l1 = brightness_params[1];
         double l2 = brightness_params[2];
         double star_T =brightness_params[0];
@@ -125,7 +125,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|| brightness_model == 10){
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4  || brightness_model == 8|| brightness_model == 10 || brightness_model == 11){
         for (int j = 0; j < 4; ++j) {
           free(bb_g[j]);
         }

spiderman/params.py

@@ -98,6 +98,10 @@ def __init__(self,brightness_model='xi'):
 			self.brightness_type= 10 # Integer model identifer
 			self.thermal= True			# Is this a thermal distribution?
 
+		elif brightness_model == 'clouds':
+			self.brightness_type= 11 # Integer model identifer
+			self.thermal= True			# Is this a thermal distribution?
+
 		else:
 			print('Brightness model "'+str(brightness_model)+'" not recognised!')
 			quit()
@@ -209,6 +213,14 @@ def format_bright_params(self):
 				print('should be',brightness_param_names)
 				quit()
 
+		elif (self.brightness_type == 11):
+			brightness_param_names = ['T_s','l1','l2','xi','T_n','delta_T','cloud']
+			try:
+				brightness_params = [self.T_s,self.l1,self.l2,self.xi,self.T_n,self.delta_T,self.clouds]
+			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')
@@ -226,8 +238,8 @@ def calc_substellar(self,t):
 		self.lambda0 = substellar[0]
 		self.phi0 = substellar[1]
 
-	def square_plot(self,ax=False,min_temp=False,max_temp=False,temp_map=False,min_bright=0.2,show_cax=True,mycmap=plt.cm.inferno,theme='black',show_axes=False,nla=100,nlo=100):
-		return splt.square_plot(self,ax=False,min_temp=False,max_temp=False,temp_map=False,min_bright=0.2,scale_planet=1.0,planet_cen=[0.0,0.0],mycmap=plt.get_cmap('inferno'),show_cax=True,theme='black',show_axes=False,nla=nla,nlo=nlo)
+	def square_plot(self,ax=False,min_temp=False,max_temp=False,temp_map=False,min_bright=0.2,show_cax=True,scale_planet=1.0,planet_cen=[0.0,0.0],mycmap=plt.get_cmap('inferno'),theme='white',show_axes=True,nla=100,nlo=100):
+		return splt.square_plot(self,ax=ax,min_temp=min_temp,max_temp=max_temp,temp_map=temp_map,min_bright=min_bright,scale_planet=scale_planet,planet_cen=planet_cen,mycmap=mycmap,show_cax=show_cax,theme=theme,show_axes=show_axes,nla=nla,nlo=nlo)
 
 	def plot_system(self,t,ax=False,min_temp=False,max_temp=False,temp_map=False,min_bright=0.2,use_phase=False,show_cax=True,mycmap=plt.cm.inferno,theme='black',show_axes=False):
 		if use_phase == True: