Skip to content

Commit

Permalink
fixed direct method bug when using small numbers of points
Browse files Browse the repository at this point in the history
  • Loading branch information
@tomlouden
tomlouden committed Oct 9, 2017
1 parent 455fd17 commit 39dd419
Show file tree
Hide file tree
Showing 8 changed files with 85 additions and 21 deletions.
9 changes: 9 additions & 0 deletions c_src/_web.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -289,6 +289,7 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)

PyArray_Descr *descr = PyArray_DescrFromType(typenum);


/* Parse the input tuple */
if (!PyArg_ParseTuple(args, "iddddiOOOidOOO", &n_layers,&lambda0,&phi0,&p_u1,&p_u2,&bright_type,&bright_obj,&teff_obj,&flux_obj,&n_star,&rp,&lo_2d_obj,&la_2d_obj,&T_2d_obj))
return NULL;
Expand Down Expand Up @@ -396,7 +397,10 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)
}


int nearest = 0;

if(bright_type == 12){
nearest = brightness_params[5];
y1_grid = malloc(sizeof(double) * (int) brightness_params[3]); // dynamic `array (size 4) of pointers to int`
for (int i = 0; i < (int) brightness_params[3]; ++i) {
y1_grid[i] = malloc(sizeof(double) * (int) brightness_params[4]);
Expand All @@ -413,6 +417,7 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)
}

if(bright_type == 13){
nearest = brightness_params[2];
y1_grid = malloc(sizeof(double) * (int) brightness_params[0]); // dynamic `array (size 4) of pointers to int`
for (int i = 0; i < (int) brightness_params[0]; ++i) {
y1_grid[i] = malloc(sizeof(double) * (int) brightness_params[1]);
Expand Down Expand Up @@ -851,7 +856,10 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)

//NEED TO UPDATE THIS WITH CORRECT STAR BRIGHTNESS VALUES OR REFLECTION MODELS WON'T BE CORRECT!//

int nearest = 0;

if(bright_type == 12){
nearest = brightness_params[5];
y1_grid = malloc(sizeof(double) * (int) brightness_params[3]); // dynamic `array (size 4) of pointers to int`
for (int i = 0; i < (int) brightness_params[3]; ++i) {
y1_grid[i] = malloc(sizeof(double) * (int) brightness_params[4]);
Expand All @@ -868,6 +876,7 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)
}

if(bright_type == 13){
nearest = brightness_params[2];
y1_grid = malloc(sizeof(double) * (int) brightness_params[0]); // dynamic `array (size 4) of pointers to int`
for (int i = 0; i < (int) brightness_params[0]; ++i) {
y1_grid[i] = malloc(sizeof(double) * (int) brightness_params[1]);
Expand Down
16 changes: 12 additions & 4 deletions c_src/generate.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,6 @@ 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]);
Expand All @@ -77,6 +76,7 @@ void map_model(double **planet,int n_layers,double lambda0, double phi0, double
mu = sqrt(1 - pow(R_mid,2));

planet[k][16] = planet[k][16]*(1 - u1*(1-mu) - u2*(pow(1-mu,2)));

}
free(center_coords);

Expand Down Expand Up @@ -216,23 +216,25 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
int *out2;
out2 = malloc(sizeof(int) * 2);

int nlo, nla;
int nlo, nla, nearest;

if(brightness_model == 12){
nlo = brightness_params[3];
nla = brightness_params[4];
nearest = brightness_params[5];
}
if(brightness_model == 13){
nlo = brightness_params[0];
nla = brightness_params[1];
nearest = brightness_params[2];
}

find_top_two(lo_2d, lo*180.0/M_PI, (int) nlo, out1);

// printf("before second %i\n",(int) brightness_params[4]);
find_top_two(la_2d, la*180.0/M_PI, (int) nla, out2);

// printf("positions %i %i %i %i\n",out1[0],out1[1],out2[0],out2[1]);
// printf("positions %i %i %i %i %i %i\n",out1[0],out1[1],out2[0],out2[1],nlo,nla);
// printf("positions %f %f %f %f\n",lo_2d[out1[0]],lo_2d[out1[1]],la_2d[out2[0]],la_2d[out2[1]]);

float *ansy, *ansy1, *ansy2;
Expand Down Expand Up @@ -268,8 +270,15 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
y12[4] = y12_grid[out1[0]][out2[1]];


// printf("%f %f %f %f %f %f\n",lo_2d[out1[0]],lo_2d[out1[1]], la_2d[out2[0]], la_2d[out2[1]], lo*180.0/M_PI, la*180.0/M_PI);
bcuint(y, y1, y2, y12, lo_2d[out1[0]],lo_2d[out1[1]], la_2d[out2[0]], la_2d[out2[1]], lo*180.0/M_PI, la*180.0/M_PI, ansy, ansy1, ansy2);

if (nearest == 1) {
int c_lo = find_minimum(lo_2d, lo*180.0/M_PI, (int) nlo);
int c_la = find_minimum(la_2d, la*180.0/M_PI, (int) nla);
*ansy = T_2d[c_lo][c_la];
}

if(brightness_model == 12){
point_T = *ansy;
point_b = bb_interp(point_T, bb_g);
Expand All @@ -278,7 +287,6 @@ double *call_map_model(double la,double lo,double lambda0, double phi0,int brigh
point_b = *ansy;
}


free(out1);
free(out2);

Expand Down
18 changes: 14 additions & 4 deletions c_src/util.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,22 +31,32 @@ int find_minimum(double *a, double v, int n) {
}

void find_top_two(double *a, double v, int n, int *out) {
int c, index1,index2;
int c, index1,index2,c_start;
double min;

min = pow((a[0] - v),2);
index1 = 0;
index2 = 0;

for (c = 1; c < n; c++) {
c_start = 0;
for (c = c_start; c < n; c++) {
if ( pow((a[c] - v),2) < min) {
index1 = c;
min = pow((a[c] - v),2);
}
}

min = pow((a[0] - v),2);
for (c = 1; c < n; c++) {
if (index1 == 0) {
min = pow((a[1] - v),2);
c_start = 2;
index2 =1;
}
else{
min = pow((a[0] - v),2);
c_start = 1;
}

for (c = c_start; c < n; c++) {
if ( pow((a[c] - v),2) < min) {
if(c != index1){
if(a[c] != a[index1]){
Expand Down
4 changes: 4 additions & 0 deletions c_src/web.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,10 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
// printf("bb_g init 2 %f\n",bb_g[0][1]);
free(coords);

int nearest = 0;

if(brightness_model == 12){
nearest = brightness_params[5];
y1_grid = malloc(sizeof(double) * (int) brightness_params[3]); // dynamic `array (size 4) of pointers to int`
for (int i = 0; i < (int) brightness_params[3]; ++i) {
y1_grid[i] = malloc(sizeof(double) * (int) brightness_params[4]);
Expand All @@ -91,6 +94,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
}

if(brightness_model == 13){
nearest = brightness_params[2];
y1_grid = malloc(sizeof(double) * (int) brightness_params[0]); // dynamic `array (size 4) of pointers to int`
for (int i = 0; i < (int) brightness_params[0]; ++i) {
y1_grid[i] = malloc(sizeof(double) * (int) brightness_params[1]);
Expand Down
25 changes: 21 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',thermal=False):
def __init__(self,brightness_model='zhang',thermal=False, nearest=None):

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

Expand All @@ -23,6 +23,8 @@ def __init__(self,brightness_model='zhang',thermal=False):
self.eclipse = True # specifies whether to include the drop in flux due to the eclipse
self.filter = False # Can use an external response file.
self.grid = [[],[],[[]]] # needed in case the "direct read" method is wanted
self.nearest = nearest # used for choosing which interpolation model to use - default is spline


if brightness_model == 'uniform brightness':
self.n_layers = 1 # The default resolution for the grid
Expand Down Expand Up @@ -248,14 +250,21 @@ def format_bright_params(self):
print('Brightness parameters incorrectly assigned')
print('should be',brightness_param_names)
quit()

elif (self.brightness_type == 12):
brightness_param_names = ['T_s','l1','l2','grid']
n_lo = len(self.grid[0])
n_la = len(self.grid[1])

# chooses whether to use a 2d spline or a nearest neighbor method. default is spline.
if self.nearest == None:
nearest = 0
elif self.nearest == True:
nearest = 1
elif self.nearest == False:
nearest = 0

try:
brightness_params = [self.T_s,self.l1,self.l2,n_lo,n_la]
brightness_params = [self.T_s,self.l1,self.l2,n_lo,n_la,nearest]
except:
print('Brightness parameters incorrectly assigned')
print('should be',brightness_param_names)
Expand All @@ -266,8 +275,16 @@ def format_bright_params(self):
n_lo = len(self.grid[0])
n_la = len(self.grid[1])

# chooses whether to use a 2d spline or a nearest neighbor method. default is spline.
if self.nearest == None:
nearest = 0
elif self.nearest == True:
nearest = 1
elif self.nearest == False:
nearest = 0

try:
brightness_params = [n_lo,n_la]
brightness_params = [n_lo,n_la,nearest]
except:
print('Brightness parameters incorrectly assigned')
print('should be',brightness_param_names)
Expand Down
3 changes: 1 addition & 2 deletions spiderman/plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ def plot_quad(spider_params,min_temp=False,max_temp=False,temp_map=False,min_bri

dp = ((max_temp-min_temp)*min_bright)


spider_params.plot_planet(0,use_phase=True,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)
spider_params.plot_planet(0.25,use_phase=True,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)
spider_params.plot_planet(0.5,use_phase=True,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)
Expand Down Expand Up @@ -535,10 +536,8 @@ def plot_planet(spider_params,t,ax=False,min_temp=False,max_temp=False,temp_map=
else:
new_ax = False


planet = sp.generate_planet(spider_params,t)


if temp_map == True:
b_i = 17
else:
Expand Down
25 changes: 21 additions & 4 deletions spiderman/read_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,10 +47,27 @@ def format_grid(longitude,latitude,temp):

LO, LA = np.unique(longitude),np.unique(latitude)
T = np.reshape(temp,(len(LO),len(LA)))
T = np.vstack((T[-5],T[-4],T[-3],T[-2],T[-1],T,T[0],T[1],T[2],T[3],T[4]))
LO = np.hstack((LO[-5:]-360,LO,LO[0:5]+360))
T = np.vstack((T[:,4],T[:,3],T[:,2],T[:,1],T[:,0],T.T,T[:,-1],T[:,-2],T[:,-3],T[:,-4],T[:,-5])).T
LA = np.hstack((-180 - LA[0:5][::-1],LA,180-LA[-5:][::-1]))

for i in range(0,len(LO)):
if LO[i] > 180:
LO[i] = LO[i] - 360
if LO[i] < -180:
LO[i] = LO[i] + 360

if len(LO) > 100:
T = np.vstack((T[-5],T[-4],T[-3],T[-2],T[-1],T,T[0],T[1],T[2],T[3],T[4]))
LO = np.hstack((LO[-5:]-360,LO,LO[0:5]+360))
else:
T = np.vstack((T[-1],T,T[0]))
LO = np.hstack((LO[-1:]-360,LO,LO[0:1]+360))

if len(LA) > 100:
T = np.vstack((T[:,4],T[:,3],T[:,2],T[:,1],T[:,0],T.T,T[:,-1],T[:,-2],T[:,-3],T[:,-4],T[:,-5])).T
LA = np.hstack((-180 - LA[0:5][::-1],LA,180-LA[-5:][::-1]))
else:
T = np.vstack((T[:,0],T.T,T[:,-1])).T
LA = np.hstack((-180 - LA[0:1][::-1],LA,180-LA[-1:][::-1]))

grid = np.array([LO,LA,T])

return grid
6 changes: 3 additions & 3 deletions spiderman/web.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,13 +57,13 @@ def generate_planet(spider_params,t,use_phase=False,stellar_grid=False,logg=4.5)
# 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,teffs,totals,len(totals),spider_params.rp,spider_params.grid[0],spider_params.grid[1],spider_params.grid[2]))
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,spider_params.grid[0].copy(),spider_params.grid[1].copy(),spider_params.grid[2].copy()))

def call_map_model(spider_params,la,lo):
# DEPRECATED - NOT CURRENTLY USED BY HIGHER LEVEL FUNCTIONS, NEEDS UPDATING

brightness_params = spider_params.format_bright_params()
return np.array(_web.call_map_model(la,lo,spider_params.brightness_type,brightness_params,spider_params.grid[0],spider_params.grid[1],spider_params.grid[2]))
return np.array(_web.call_map_model(la,lo,spider_params.brightness_type,brightness_params,spider_params.grid[0].copy(),spider_params.grid[1].copy(),spider_params.grid[2].copy()))

def blocked(n_layers,x2,y2,r2):
return _web.blocked(n_layers,x2,y2,r2)
Expand Down Expand Up @@ -104,7 +104,7 @@ def lightcurve(t,spider_params,stellar_grid=False,logg=4.5,use_phase=False):

n_wvls = len(filter[0])

out = _web.lightcurve(spider_params.n_layers,t,spider_params.t0,spider_params.per,spider_params.a_abs,spider_params.inc,spider_params.ecc,spider_params.w,spider_params.a,spider_params.rp,spider_params.p_u1,spider_params.p_u2,spider_params.brightness_type,brightness_params,teffs,totals,len(totals), spider_params.eclipse, filter[0], filter[1], n_wvls,use_filter,spider_params.grid[0],spider_params.grid[1],spider_params.grid[2])
out = _web.lightcurve(spider_params.n_layers,t,spider_params.t0,spider_params.per,spider_params.a_abs,spider_params.inc,spider_params.ecc,spider_params.w,spider_params.a,spider_params.rp,spider_params.p_u1,spider_params.p_u2,spider_params.brightness_type,brightness_params,teffs,totals,len(totals), spider_params.eclipse, filter[0], filter[1], n_wvls,use_filter,spider_params.grid[0].copy(),spider_params.grid[1].copy(),spider_params.grid[2].copy())

return np.array(out)

Expand Down

0 comments on commit 39dd419

Please sign in to comment.