Skip to content

Commit

Permalink
changed final normalisation of spherical modes and plotting colors
Browse files Browse the repository at this point in the history
  • Loading branch information
@tomlouden
tomlouden committed Oct 3, 2017
1 parent 8d601f1 commit 45724bf
Show file tree
Hide file tree
Showing 8 changed files with 86 additions and 43 deletions.
4 changes: 2 additions & 2 deletions c_src/_web.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -357,7 +357,7 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)
star_surface_bright = star_bright/(M_PI*pow(r2,2));


if(bright_type == 1 || bright_type == 3 || bright_type == 4|| bright_type == 6 || bright_type == 8 || bright_type == 10 || bright_type == 11 || bright_type == 12){
if(bright_type == 1 || bright_type == 3 || bright_type == 4|| bright_type == 6 || bright_type == 8 || bright_type == 10 || bright_type == 11 || bright_type == 12 || bright_type == 14){
l1 = brightness_params[1];
l2 = brightness_params[2];

Expand Down Expand Up @@ -819,7 +819,7 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)
}


if(bright_type == 1 || bright_type == 3 || bright_type == 4 || bright_type == 6 || bright_type == 8 || bright_type == 10 || bright_type == 11 || bright_type == 12){
if(bright_type == 1 || bright_type == 3 || bright_type == 4 || bright_type == 6 || bright_type == 8 || bright_type == 10 || bright_type == 11 || bright_type == 12 || bright_type == 14){
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,use_filter, n_wvls, wvl_grid);
Expand Down
6 changes: 6 additions & 0 deletions c_src/blackbody.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,12 @@ double bb_flux(double l1, double l2, double T,int n_segments, int use_filter, in
// ypp = spline_cubic_set( n_wvls, wvl_g[0], wvl_g[1], 0, 0, 0, 0 );
// }

if(T< 0.0){
printf("CANNOT DO NEGATIVE TEMPERATURES!\n");
exit(0);
}


for (int k = 0; k <n_segments; ++k) {
l_lower = (l1+k*wvl_int);
l_upper = (l1+(k+1)*wvl_int);
Expand Down
14 changes: 9 additions & 5 deletions c_src/brightness_maps.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,17 +259,21 @@ 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 spherical(double lat, double lon, double *a, int therm_flag){
double x_vec[1];
double fx2;
double *fx2_vec;
double theta = (M_PI/2.0) - (lat+a[1]);
double phi = M_PI + (lon+a[2]);
int orders = a[0];
int avoid_neg = 0; // This should probably be a user setting.
double norm;

int k = 3;
if(therm_flag == 1){
k = 6;
}

double theta = (M_PI/2.0) - (lat+a[k-1]);
double phi = M_PI + (lon+a[k-2]);
int orders = a[k-3];

x_vec[0] = cos(theta);

Expand Down Expand Up @@ -301,7 +305,7 @@ double spherical(double lat, double lon, double *a){
if(avoid_neg == 1){
return pow(pow(val,2),0.5)/M_PI;
}
return val/M_PI;
return val;
}

double kreidberg_2016(double lat, double lon, double insol, double albedo, double redist){
Expand Down
2 changes: 1 addition & 1 deletion c_src/brightness_maps.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ 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 spherical(double lat, double lon, double *a);
double spherical(double lat, double lon, double *a, int therm_flag);
double kreidberg_2016(double lat, double lon, double insol, double albedo, double redist);
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);
42 changes: 27 additions & 15 deletions c_src/generate.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,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|| brightness_model == 11|| brightness_model == 12) {
if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6|| brightness_model == 8|| brightness_model == 10|| brightness_model == 11|| brightness_model == 12 || brightness_model == 14) {
l1 = brightness_params[1];
l2 = brightness_params[2];
}
Expand Down Expand Up @@ -65,8 +65,10 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double

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,lo_2d,la_2d,T_2d,y1_grid,y2_grid,y12_grid);


planet[k][16] = vals[0];
planet[k][17] = vals[1];
// printf("%f %f inside",planet[k][16],planet[k][17]);

free(vals);

Expand All @@ -88,19 +90,19 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
double R = 1.0;
int n_bb_seg = 10;

point_T = 0.0;
point_b = 0.0;

// 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|| brightness_model == 10|| brightness_model == 11|| brightness_model == 12) {
l1 = brightness_params[1];
l2 = brightness_params[2];
}


if(brightness_model == 0){
point_b = Uniform_b(brightness_params[0]);
}

if(brightness_model == 1){
else if(brightness_model == 1){
point_T = Uniform_T(brightness_params[3]);
point_b = bb_interp(point_T, bb_g);
}
Expand All @@ -118,26 +120,33 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
// point_b = bb_interp(point_T, bb_g);
// }

if(brightness_model == 4){
else if(brightness_model == 4){
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);
}
if(brightness_model == 5){
point_b = spherical(la,lo,brightness_params);
else if(brightness_model == 5){

point_b = spherical(la,lo,brightness_params,0);
}

if(brightness_model == 6){
else if(brightness_model == 14){

point_T = spherical(la,lo,brightness_params,1);
point_b = bb_interp(point_T, bb_g);
}

else if(brightness_model == 6){
double insol = brightness_params[3];
double albedo = brightness_params[4];
double redist = brightness_params[5];
point_T = kreidberg_2016(la, lo, insol, albedo, redist);
point_b = bb_interp(point_T, bb_g);
}

if(brightness_model == 2 || brightness_model == 7){
else if(brightness_model == 2 || brightness_model == 7){
double la0 = brightness_params[0];
double lo0 = brightness_params[1];
double p_b = brightness_params[2];
Expand All @@ -147,7 +156,7 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
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){
else if(brightness_model == 3 || brightness_model == 8){
double la0 = brightness_params[4];
double lo0 = brightness_params[5];
double p_T = brightness_params[6];
Expand All @@ -158,13 +167,13 @@ 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,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){
else 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){
else if(brightness_model == 10){
double xi =brightness_params[3];
double T_n =brightness_params[4];
double delta_T =brightness_params[5];
Expand All @@ -178,7 +187,7 @@ 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){
else if(brightness_model == 11){
double xi =brightness_params[3];
double T_n =brightness_params[4];
double delta_T =brightness_params[5];
Expand All @@ -193,7 +202,7 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
}
}

if(brightness_model == 12 || brightness_model == 13){
else if(brightness_model == 12 || brightness_model == 13){

/* printf("\n");
printf("%f\n",lo*180.0/M_PI);
Expand Down Expand Up @@ -283,6 +292,9 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
free(ansy2);

}
else{
printf("BRIGHTNESS MODEL NOT RECOGNISED\n");
}

// printf("%f %f %f\n",ars, star_bright,insol);
output = malloc(sizeof(double) * 2); // dynamic `array (size 2) of pointers to double`
Expand Down
4 changes: 2 additions & 2 deletions c_src/web.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,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|| brightness_model == 11|| brightness_model == 12){
if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6 || brightness_model == 8|| brightness_model == 10|| brightness_model == 11|| brightness_model == 12 || brightness_model == 14){
double l1 = brightness_params[1];
double l2 = brightness_params[2];
double star_T =brightness_params[0];
Expand Down Expand Up @@ -169,7 +169,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 == 6 || brightness_model == 8|| brightness_model == 10 || brightness_model == 11|| brightness_model == 12){
if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6 || brightness_model == 8|| brightness_model == 10 || brightness_model == 11|| brightness_model == 12 || brightness_model == 14){
for (int j = 0; j < 4; ++j) {
free(bb_g[j]);
}
Expand Down
24 changes: 20 additions & 4 deletions spiderman/params.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@

class ModelParams(object):

def __init__(self,brightness_model='zhang'):
def __init__(self,brightness_model='zhang',thermal=False):

self.n_layers = 5 # The default resolution for the grid

Expand Down Expand Up @@ -67,9 +67,12 @@ def __init__(self,brightness_model='zhang'):
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?
self.thermal= thermal # Is this a thermal distribution?
if thermal == True:
self.brightness_type= 14 # Integer model identifier
else:
self.brightness_type= 5 # Integer model identifier

elif brightness_model == 'kreidberg':
self.brightness_type= 6 # Integer model identifer
Expand Down Expand Up @@ -164,10 +167,22 @@ def format_bright_params(self):
print('Brightness parameters incorrectly assigned')
print('should be',brightness_param_names)
quit()
elif (self.brightness_type == 5):

elif (self.brightness_type == 5 or self.brightness_type == 14):

brightness_param_names = ['degree','sph','la0','lo0']
if self.brightness_type == 14:
brightness_param_names = ['T_s','l1','l2'] + brightness_param_names


brightness_params = [self.degree,self.la0,self.lo0] + self.sph
if self.brightness_type == 14:
brightness_params = [self.T_s,self.l1,self.l2] + brightness_params

try:
brightness_params = [self.degree,self.la0,self.lo0] + self.sph
if self.brightness_type == 14:
brightness_params = [self.T_s,self.l1,self.l2] + brightness_params
except:
print('Brightness parameters incorrectly assigned')
print('should be',brightness_param_names)
Expand All @@ -178,6 +193,7 @@ def format_bright_params(self):
print('You have not specified the correct number of mode coefficients!')
print('You gave '+str(int(len(self.sph)))+', there should be '+str(int(total_modes)))
quit()

elif (self.brightness_type == 6):
brightness_param_names = ['T_s','l1','l2','insol','albedo','redist']
try:
Expand Down
33 changes: 19 additions & 14 deletions spiderman/plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,8 @@ def plot_quad(spider_params,min_temp=False,max_temp=False,temp_map=False,min_bri
blims2 = spider_params.get_lims(0.25,temp_map=temp_map,use_phase=True)
blims3 = spider_params.get_lims(0.5,temp_map=temp_map,use_phase=True)
blims4 = spider_params.get_lims(0.75,temp_map=temp_map,use_phase=True)
min_temp = np.min(np.array([blims1,blims2,blims3,blims4]))
max_temp = np.max(np.array([blims1,blims2,blims3,blims4]))
min_temp = 0.9*np.min(np.array([blims1,blims2,blims3,blims4]))
max_temp = 1.1*np.max(np.array([blims1,blims2,blims3,blims4]))
if min_temp == 0.0:
min_temp = 1e-19

Expand All @@ -105,18 +105,18 @@ def plot_quad(spider_params,min_temp=False,max_temp=False,temp_map=False,min_bri

# divider = make_axes_locatable(fig)

min_val = (dp + min_temp-min_temp)/(dp + max_temp-min_temp)
max_val = (dp + max_temp-min_temp)/(dp + max_temp-min_temp)
# min_val = (dp + min_temp-min_temp)/(dp + max_temp-min_temp)
# max_val = (dp + max_temp-min_temp)/(dp + max_temp-min_temp)

zero_temp = min_temp - dp
# zero_temp = min_temp - dp

# zero_temp = min_temp
zero_temp = min_temp


if temp_map == True:
data = [np.linspace(zero_temp,max_temp,1000)]*2
else:
data = [np.linspace(zero_temp/max_temp,max_temp/max_temp,1000)]*2
data = [np.linspace(min_temp / (0.5*(min_temp + max_temp)),max_temp / (0.5*(min_temp + max_temp)),1000)]*2
fake, fake_ax = plt.subplots()
mycax = fake_ax.imshow(data, interpolation='none', cmap=mycmap)
plt.close(fake)
Expand Down Expand Up @@ -431,6 +431,9 @@ def plot_square(spider_params,ax=False,min_temp=False,max_temp=False,temp_map=Fa
fluxes += [row]
fluxes = np.array(fluxes)

if temp_map == False:
fluxes = fluxes/np.max(fluxes)

lala, lolo = np.meshgrid(los,las)

plt.plot([0],[0],'x',color=('#0cff0c'),ms=10,mew=2)
Expand Down Expand Up @@ -543,11 +546,11 @@ def plot_planet(spider_params,t,ax=False,min_temp=False,max_temp=False,temp_map=

if min_temp == False:
temps = planet[:,b_i]
min_temp = np.min(temps)
max_temp = np.max(temps)
min_temp = np.min(temps)*0.9
max_temp = np.max(temps)*1.1

temps = planet[:,b_i]
# print('min ',np.min(temps), 'max',np.max(temps))


dp = ((max_temp-min_temp)*min_bright)

Expand All @@ -559,13 +562,15 @@ def plot_planet(spider_params,t,ax=False,min_temp=False,max_temp=False,temp_map=
for j in range (0,len(planet)):


if dp == 0.0:
val = 1
else:
val = (dp + planet[j][b_i]-min_temp)/(dp + max_temp-min_temp)
# if dp == 0.0:
# val = 1
# else:
# val = (dp + planet[j][b_i]-min_temp)/(dp + max_temp-min_temp)

# print(val,(dp + planet[j][b_i]-min_temp),(dp + max_temp-min_temp))

val = (planet[j][b_i]-min_temp)/(max_temp-min_temp)

c = mycmap(val)

n = planet[j]
Expand Down

0 comments on commit 45724bf

Please sign in to comment.