ivs.roche.local

1 """ 2 Derive local quantities such as effective temperature and surface gravity for Roche potentials. 3 """ 4 import numpy as np 5 import pylab as pl 6 from numpy import pi,cos,sin,sqrt 7 from scipy.optimize import newton 8 try: 9 from scipy.spatial import Delaunay 10 except ImportError: 11 print 'No Delaunay' 12 13 from ivs.sed import limbdark 14 from ivs.coordinates import vectors 15 16 17 #{ Gridding functions 18

19 -def get_grid(*args,**kwargs):

20 """ 21 Construct a coordinate grid 22 23 If you give two resolutions, the first is for theta, the second for phi 24 25 @param args: one or two integers indicating number of grid points in theta 26 and phi direction 27 @type args: integer 28 @keyword gtype: grid type ('spher' or 'delaunay') 29 @type gtype: str 30 @return: theta,phi(,grid) 31 """ 32 gtype = kwargs.get('gtype','spherical') 33 full = kwargs.get('full',False) 34 35 if 'spher' in gtype.lower(): 36 if len(args)==1: res1 = res2 = args[0] # same resolution for both coordinates 37 else: res1,res2 = args # different resolution 38 39 dtheta = pi/res1 # step size coordinate 1 40 dphi = pi/res2 # step size coordinate 2 41 42 #-- full grid or only one quadrant 43 if full: 44 theta0,thetan = dtheta, pi-dtheta 45 phi0,phin = 0,2*pi-2*dphi 46 else: 47 theta0,thetan = dtheta/4,pi/2-dtheta/4 48 phi0,phin = 0,pi-dphi 49 50 theta,phi = np.mgrid[theta0:thetan:res1*1j,phi0:phin:res2*1j] 51 return theta,phi 52 53 elif 'cil' in gtype.lower(): 54 if len(args)==1: res1 = res2 = args[0] # same resolution for both coordinates 55 else: res1,res2 = args # different resolution 56 57 dsintheta = 2./res1 # step size sin(theta) 58 dphi = pi/res2 # step size coordinate 2 59 60 #-- full grid or only one quadrant 61 if full: 62 sintheta0,sinthetan = -1+dsintheta, 1.-dsintheta 63 phi0,phin = 0,2*pi-dphi 64 else: 65 sintheta0,sinthetan = dsintheta,1-dsintheta 66 phi0,phin = 0,pi-dphi 67 68 sintheta,phi = np.mgrid[sintheta0:sinthetan:res1*1j,phi0:phin:res2*1j] 69 return np.arccos(sintheta),phi 70 71 elif 'delaunay' in gtype.lower(): 72 #-- resolution doesn't matter that much anymore 73 if len(args)==1: res1 = res2 = args[0] # same resolution for both coordinates 74 else: res1,res2 = args # different resolution 75 76 u,v = np.random.uniform(size=res1*res2),np.random.uniform(size=res1*res2) 77 phi = 2*pi*u 78 theta = np.arccos(2*v-1) 79 80 x = sin(theta)*sin(phi) 81 y = sin(theta)*cos(phi) 82 z = cos(theta) 83 84 points = np.array([x,y,z]).T 85 grid = Delaunay(points) 86 centers = np.zeros((len(grid.convex_hull),3)) 87 for i,indices in enumerate(grid.convex_hull): 88 centers[i] = [x[indices].sum()/3,y[indices].sum()/3,z[indices].sum()/3] 89 theta,phi = np.arccos(centers[:,2]),np.arctan2(centers[:,1],centers[:,0])+pi 90 91 return theta,phi,grid 92 93 elif 'tri' in gtype.lower(): 94 #-- resolution doesn't matter that much anymore 95 if len(args)==1: res1 = res2 = args[0] # same resolution for both coordinates 96 else: res1,res2 = args # different resolution 97 98 u,v = np.random.uniform(size=res1*res2),np.random.uniform(size=res1*res2) 99 phi = 2*pi*u 100 theta = np.arccos(2*v-1) 101 return theta,phi

102 103 104 105 106 107 108

109 -def stitch_grid(theta,phi,*quant,**kwargs):

110 """ 111 Stitch a grid together that was originally defined on 1 quadrant. 112 113 We add the three other quandrants. 114 """ 115 seamless = kwargs.get('seamless',False) 116 vtype = kwargs.get('vtype',['scalar' for i in quant]) 117 gtype = kwargs.get('gtype','spher') 118 ravel = kwargs.get('ravel',False) 119 120 if gtype == 'spher': 121 #-- basic coordinates 122 alltheta = np.vstack([np.hstack([theta,theta]),np.hstack([theta+pi/2,theta+pi/2])]) 123 allphi = np.vstack([np.hstack([phi,phi+pi]),np.hstack([phi,phi+pi])]) 124 125 #-- all other quantities 126 allquan = [] 127 for i,iquant in enumerate(quant): 128 #-- if they are scalar values, they do not change direction 129 top1,top2 = iquant,iquant[:,::-1] 130 bot1,bot2 = iquant[::-1],iquant[::-1][:,::-1] 131 if vtype[i]=='scalar': 132 allquan.append(np.vstack([np.hstack([top1,top2]),np.hstack([bot1,bot2])])) 133 #-- vector components do change direction 134 elif vtype[i]=='x': 135 allquan.append(np.vstack([np.hstack([top1,top2]),np.hstack([bot1,bot2])])) 136 elif vtype[i]=='y': 137 allquan.append(np.vstack([np.hstack([top1,-top2]),np.hstack([bot1,-bot2])])) 138 elif vtype[i]=='z': 139 allquan.append(np.vstack([np.hstack([top1,top2]),np.hstack([-bot1,-bot2])])) 140 elif vtype[i]=='vx': 141 allquan.append(np.vstack([np.hstack([top1,-top2]),np.hstack([bot1,-bot2])])) 142 elif vtype[i]=='vy': 143 allquan.append(np.vstack([np.hstack([top1,top2]),np.hstack([bot1,bot2])])) 144 elif vtype[i]=='vz': 145 allquan.append(np.vstack([np.hstack([top1,top2]),np.hstack([-bot1,-bot2])])) 146 147 out = [alltheta,allphi]+allquan 148 149 #-- for plotting reasons, remove the vertical seam and the the bottom hole. 150 if seamless: 151 #-- vertical seam 152 out = [np.column_stack([i,i[:,0]]) for i in out] 153 out[1][:,-1] += 2*pi 154 #-- fill bottom hole 155 out = [np.vstack([i,i[0]]) for i in out] 156 out[0][-1] += pi 157 158 #-- ravel to 1d arrays if asked for 159 if ravel: 160 out = [i.ravel() for i in out] 161 else: 162 out = [theta,phi] + list(quant) 163 164 return out

165 166 #}

167 -def surface_normals(r,phi,theta,grid,gtype='spher'):

168 """ 169 Numerically compute surface normals of a grid (in absence of analytical alternative). 170 171 Also computes the surface elements, making L{surface_elements} obsolete. 172 """ 173 if gtype=='spher': 174 raise NotImplementedError 175 elif gtype=='delaunay': 176 raise NotImplementedError 177 elif gtype=='triangular': 178 #-- compute the angle between the surface normal and the radius vector 179 x,y,z = vectors.spher2cart_coord(r,phi,theta) 180 181 centers = np.zeros((len(grid.convex_hull),3)) 182 normals = np.zeros((len(grid.convex_hull),3)) 183 sizes = np.zeros(len(grid.convex_hull)) 184 185 #vertx,verty,vertz = points.T 186 187 #-- compute centers,normals and sizes 188 for i,indices in enumerate(grid.convex_hull): 189 #-- center is triangle's barycenter 190 centers[i] = [x[indices].sum()/3,y[indices].sum()/3,z[indices].sum()/3] 191 #-- size is size of triangle 192 a = sqrt((x[indices[0]]-x[indices[1]])**2 + (y[indices[0]]-y[indices[1]])**2 + (z[indices[0]]-z[indices[1]])**2) 193 b = sqrt((x[indices[0]]-x[indices[2]])**2 + (y[indices[0]]-y[indices[2]])**2 + (z[indices[0]]-z[indices[2]])**2) 194 c = sqrt((x[indices[1]]-x[indices[2]])**2 + (y[indices[1]]-y[indices[2]])**2 + (z[indices[1]]-z[indices[2]])**2) 195 s = 0.5*(a+b+c) 196 sizes[i] = sqrt( s*(s-a)*(s-b)*(s-c)) 197 #-- normal is cross product of two sides 198 side1 = [x[indices[1]]-x[indices[0]],y[indices[1]]-y[indices[0]],z[indices[1]]-z[indices[0]]] 199 side2 = [x[indices[2]]-x[indices[0]],y[indices[2]]-y[indices[0]],z[indices[2]]-z[indices[0]]] 200 normals[i] = np.cross(side1,side2) 201 202 #-- make sure the normal is pointed outwards 203 normal_r,normal_phi,normal_theta = vectors.cart2spher(centers.T,normals.T) 204 normal_r = np.abs(normal_r) 205 centers_sph = vectors.cart2spher_coord(*centers.T) 206 normals = np.array(vectors.spher2cart(centers_sph,(normal_r,normal_phi,normal_theta))) 207 208 #-- normalise and compute angles 209 normals_T = normals.T 210 normals = normals_T / vectors.norm(normals_T) 211 #cos_gamma = vectors.cos_angle(a,normals) 212 print centers.shape,sizes.shape,normals.shape 213 return centers, sizes, normals#, cos_gamma

214 215 #{ Derivation of local quantities 216

217 -def surface_elements((r,mygrid),(surfnormal_x,surfnormal_y,surfnormal_z),gtype='spher'):

218 """ 219 Compute surface area of elements in a grid. 220 221 theta,phi must be generated like mgrid(theta_range,phi_range) 222 223 usually, the surfnormals are acquired via differentiation of a gravity potential, 224 and is then equal to the *negative* of the local surface gravity. 225 """ 226 theta,phi = mygrid[:2] 227 if gtype=='spher': 228 #-- compute the grid size at each location 229 dtheta = theta[1:]-theta[:-1] 230 dtheta = np.vstack([dtheta,dtheta[-1]]) 231 232 dphi = phi[:,1:]-phi[:,:-1] 233 dphi = np.column_stack([dphi,dphi[:,-1]]) 234 #-- compute the angle between the surface normal and the radius vector 235 x,y,z = vectors.spher2cart_coord(r,phi,theta) 236 237 a = np.array([x,y,z]) 238 b = np.array([surfnormal_x,surfnormal_y,surfnormal_z]) 239 240 cos_gamma = vectors.cos_angle(a,b) 241 242 return r**2 * sin(theta) * dtheta * dphi / cos_gamma, cos_gamma 243 244 elif gtype=='delaunay': 245 #-- compute the angle between the surface normal and the radius vector 246 x,y,z = vectors.spher2cart_coord(r,phi,theta) 247 a = np.array([x,y,z]) 248 b = np.array([surfnormal_x,surfnormal_y,surfnormal_z]) 249 cos_gamma = vectors.cos_angle(a,b) 250 251 delaunay_grid = mygrid[2] 252 253 sizes = np.zeros(len(delaunay_grid.convex_hull)) 254 points = delaunay_grid.points 255 vertx,verty,vertz = points.T 256 257 #from enthought.mayavi import mlab 258 #mlab.figure() 259 #mlab.triangular_mesh(vertx,verty,vertz,delaunay_grid.convex_hull,scalars=np.ones_like(vertx),colormap='gray',representation='wireframe') 260 #mlab.points3d(x/r,y/r,z/r,scale_factor=0.02) 261 262 centers = np.zeros((len(delaunay_grid.convex_hull),3)) 263 for i,indices in enumerate(delaunay_grid.convex_hull): 264 #centers[i] = [vertx[indices].sum()/3,verty[indices].sum()/3,vertz[indices].sum()/3] 265 a = sqrt((vertx[indices[0]]-vertx[indices[1]])**2 + (verty[indices[0]]-verty[indices[1]])**2 + (vertz[indices[0]]-vertz[indices[1]])**2) 266 b = sqrt((vertx[indices[0]]-vertx[indices[2]])**2 + (verty[indices[0]]-verty[indices[2]])**2 + (vertz[indices[0]]-vertz[indices[2]])**2) 267 c = sqrt((vertx[indices[1]]-vertx[indices[2]])**2 + (verty[indices[1]]-verty[indices[2]])**2 + (vertz[indices[1]]-vertz[indices[2]])**2) 268 s = 0.5*(a+b+c) 269 sizes[i] = sqrt( s*(s-a)*(s-b)*(s-c)) 270 271 #theta,phi = np.arccos(centers[:,2]),np.arctan2(centers[:,1],centers[:,0])+pi 272 #mlab.points3d(centers[:,0],centers[:,1],centers[:,2],sizes,scale_factor=0.05,scale_mode='none',colormap='RdBu') 273 #mlab.show() 274 275 #pl.show() 276 277 return sizes*r**2, cos_gamma

278

279 -def temperature(surface_gravity,g_pole,T_pole,beta=1.):

280 """ 281 Calculate local temperature. 282 283 beta is gravity darkening parameter. 284 """ 285 Grav = abs(surface_gravity/g_pole)**beta 286 Teff = Grav**0.25 * T_pole 287 return Teff

288

289 -def intensity(teff,grav,mu=None,photband='OPEN.BOL'):

290 """ 291 Calculate local intensity. 292 293 beta is gravity darkening parameter. 294 """ 295 if mu is None: 296 mu = np.ones_like(teff) 297 if (teff<3500).any() or np.isnan(teff).any(): 298 print 'WARNING: point outside of grid, minimum temperature is 3500K' 299 teff = np.where((teff<3500) | np.isnan(teff),3500,teff) 300 if (grav<0.01).any() or np.isnan(grav).any(): 301 print 'WARNING: point outside of grid, minimum gravity is 0 dex' 302 grav = np.where((np.log10(grav*100)<0.) | np.isnan(grav),0.01,grav) 303 intens = np.array([limbdark.get_itable(teff=iteff,logg=np.log10(igrav*100),absolute=True,mu=imu,photbands=[photband])[0] for iteff,igrav,imu in zip(teff.ravel(),grav.ravel(),mu.ravel())]) 304 return intens.reshape(teff.shape)

305 306

307 -def projected_intensity(teff,gravity,areas,line_of_sight,photband='OPEN.BOL'):

308 """ 309 Compute projected intensity in the line of sight. 310 311 gravity is vector directed inwards in the star 312 line of sight is vector. 313 """ 314 ones = np.ones_like(gravity[0]) 315 losx = line_of_sight[0]*ones 316 losy = line_of_sight[1]*ones 317 losz = line_of_sight[2]*ones 318 angles = vectors.angle(-gravity,np.array([losx,losy,losz])) 319 mus = cos(angles) 320 grav_ = vectors.norm(gravity) 321 #-- intensity is less if we look at the limb 322 intens = intensity(teff,grav_,mu=mus,photband=photband) 323 #-- intensity is less if the surface element area is small (it does not need 324 # to be projected anymore!) 325 return intens*areas*mus,mus

326 327

328 -def project(star,view_long=(0,0,0),view_lat=(pi/2,0,0),photband='OPEN.BOL', 329 only_visible=False,plot_sort=False,scale_factor=1.):

330 """ 331 Project and transform coordinates and vectors to align with the line-of-sight. 332 333 Parameter C{star} should be a record array containing fields 'teff','gravx', 334 'gravy','gravz','areas','vx','vy','vz' 335 336 and either you suply ('r','theta','phi') or ('x','y','z') 337 338 The XY direction is then the line-of-sight, and the YZ plane is the plane 339 of the sky. 340 341 An extra column 'projflux' and 'eyeflux' will be added. Projected flux 342 takes care of limb darkening, and projected surface area. Eye flux only 343 takes care of limbdarkening, and should only be used for plotting reasons. 344 345 view_long[0] of 0 means looking in the XY line, pi means looking in the YX line. 346 view_lat[0] of pi/2 means edge on, 0 or pi is pole-on. 347 348 This function updates all Cartesian coordinates present in the star, but not 349 the polar coordinates! The projected fluxes are added as a field 'projflux' 350 to the returned record array. 351 352 If you set 'only_visible' to True, only the information on the visible parts 353 of the star will be contained. 354 355 If you set 'plot_sort' to True, the arrays will be returned in a sorted order, 356 where the areas at the back come first. This is especially handy for plotting. 357 358 @parameters star: record array containing all necessary information on the 359 star 360 @type star: numpy record array 361 @parameter view_long: longitude viewing angle (radians) and coordinate zeropoint 362 @type view_long: tuple floats (radians,x,y) 363 @parameter view_lat: inclination viewing angle (radians) and coordinate zeropoint 364 @type view_lat: tuple floats (radians,x,z) 365 @parameter photband: photometric passband 366 @type photband: string 367 @parameter only_visible: flag to return only information on visible surface elements 368 @type only_visible: boolean 369 @parameter plot_sort: flag to sort the surface elements from back to front 370 @type plot_sort: boolean 371 """ 372 myshape = star['gravx'].shape 373 gravx,gravy,gravz = np.array([star['gravx'].ravel(),star['gravy'].ravel(),star['gravz'].ravel()]) 374 areas = star['areas'].ravel() 375 teff = star['teff'].ravel() 376 vx,vy,vz = star['vx'].ravel(),star['vy'].ravel(),star['vz'].ravel() 377 #-- if 'x' is not in the star's record array, we assume the polar coordinates 378 # are in there and convert them to Cartesian coordinates 379 if not 'x' in star.dtype.names: 380 x,y,z = vectors.spher2cart_coord(star['r'].ravel(),star['phi'].ravel(),star['theta'].ravel()) 381 else: 382 x,y,z = star['x'].ravel(),star['y'].ravel(),star['z'].ravel(), 383 384 #-- first we rotate in the XY plane (only for surface coordinates is the 385 # coordinate zeropoint important, the rest are vectors!): 386 x,y = vectors.rotate(x,y,view_long[0],x0=view_long[1],y0=view_long[2]) 387 gravx,gravy = vectors.rotate(gravx,gravy,view_long[0]) 388 vx,vy = vectors.rotate(vx,vy,view_long[0]) 389 #-- then we rotate in the YZ plane: 390 if view_lat[0]!=pi/2: 391 rot_i = -(pi/2 - view_lat[0]) 392 x,z = vectors.rotate(x,z,rot_i) 393 gravx,gravz = vectors.rotate(gravx,gravz,rot_i) 394 vx,vz = vectors.rotate(vx,vz,rot_i) 395 #-- ... and project the fluxes in the line of sight, which is now in the XY 396 # direction: 397 view_vector = np.array([1.,0,0])#np.array([-sin(pi/2),0,-cos(pi/2)]) 398 grav_local = np.array([gravx,gravy,gravz]) 399 proj_flux,mus = projected_intensity(teff,grav_local,areas,view_vector,photband=photband) 400 401 #-- we now construct a copy of the star record array with the changed 402 # coordinates 403 new_star = star.copy() 404 new_star['gravx'],new_star['gravy'],new_star['gravz'] = gravx,gravy,gravz 405 new_star['vx'],new_star['vy'],new_star['vz'] = vx,vy,vz 406 if 'x' in star.dtype.names: 407 new_star['x'],new_star['y'],new_star['z'] = x*scale_factor,y*scale_factor,z*scale_factor 408 else: 409 new_star = pl.mlab.rec_append_fields(new_star,'x',x*scale_factor) 410 new_star = pl.mlab.rec_append_fields(new_star,'y',y*scale_factor) 411 new_star = pl.mlab.rec_append_fields(new_star,'z',z*scale_factor) 412 new_star = pl.mlab.rec_append_fields(new_star,'projflux',proj_flux) 413 new_star = pl.mlab.rec_append_fields(new_star,'eyeflux',proj_flux/areas) 414 new_star = pl.mlab.rec_append_fields(new_star,'mu',mus) 415 416 #-- clip visible areas and sort in plotting order if necessary 417 if only_visible: 418 new_star = new_star[-np.isnan(new_star['projflux'])] 419 if plot_sort: 420 new_star = new_star[np.argsort(new_star['x'])] 421 return new_star

422

Source Code for Module ivs.roche.local