ivs.observations.barycentric

1 # -*- coding: utf-8 -*- 2 """ 3 Compute the influence of earth's orbit around the sun, and absolute stellar velocities 4 """ 5 #from numpy import * 6 import numpy as np 7 from ivs.units.conversions import convert 8

9 -def baryvel(dje, deq=0):

10 """ 11 Calculates heliocentric and barycentric velocity components of Earth. 12 13 This function takes into account the Earth-Moon motion, and is useful for 14 radial velocity work to an accuracy of ~1 m/s. 15 16 dvel_hel, dvel_bary = baryvel(dje, deq) 17 18 Inputs: 19 @param dje: Julian ephemeris date. 20 @type dje: array 21 @param deq: epoch of mean equinox of dvelh and dvelb. If deq = 0 then deq is assumed to be equal to dje 22 @type deq: float 23 @return: heliocentric and barycentric velocity components in km/s (The 3-vectors dvelh and dvelb are given in a right-handed coordinate system with the +X axis toward the Vernal Equinox, and +Z axis toward the celestial pole.) 24 @rtype: array (2X3) 25 26 Functions called: 27 premat() -- computes precession matrix 28 29 NOTES: 30 Algorithm taken from FORTRAN program of Stumpff (1980, A&A Suppl, 41,1) 31 Stumpf claimed an accuracy of 42 cm/s for the velocity. A comparison with 32 the JPL FORTRAN planetary ephemeris program PLEPH found agreement to within 33 about 65 cm/s between 1986 and 1994 34 35 REVISION HISTORY: 36 Jeff Valenti, U.C. Berkeley Translated BARVEL.FOR to IDL. 37 W. Landsman, Cleaned up program sent by Chris McCarthy (SfSU) June 1994 38 Converted to IDL V5.0 W. Landsman September 1997 39 Added /JPL keyword W. Landsman July 2001 40 Documentation update W. Landsman Dec 2005 41 Converted to Python S. Koposov 2009-2010 42 Adapted for IvS library 43 """ 44 45 #Define constants 46 dc2pi = 2 * np.pi 47 cc2pi = 2 * np.pi 48 dc1 = 1.0e0 49 dcto = 2415020.0e0 50 dcjul = 36525.0e0 #days in Julian year 51 dcbes = 0.313e0 52 dctrop = 365.24219572e0 #days in tropical year (...572 insig) 53 dc1900 = 1900.0e0 54 au = 1.4959787e8 55 56 #Constants dcfel(i,k) of fast changing elements. 57 dcfel = np.array([1.7400353e00, 6.2833195099091e02, 5.2796e-6, 6.2565836e00, 6.2830194572674e02, -2.6180e-6, 4.7199666e00, 8.3997091449254e03, -1.9780e-5, 1.9636505e-1, 8.4334662911720e03, -5.6044e-5, 4.1547339e00, 5.2993466764997e01, 5.8845e-6, 4.6524223e00, 2.1354275911213e01, 5.6797e-6, 4.2620486e00, 7.5025342197656e00, 5.5317e-6, 1.4740694e00, 3.8377331909193e00, 5.6093e-6]) 58 dcfel = np.reshape(dcfel, (8, 3)) 59 60 #constants dceps and ccsel(i,k) of slowly changing elements. 61 dceps = np.array([4.093198e-1, -2.271110e-4, -2.860401e-8]) 62 ccsel = np.array([1.675104e-2, -4.179579e-5, -1.260516e-7, 2.220221e-1, 2.809917e-2, 1.852532e-5, 1.589963e00, 3.418075e-2, 1.430200e-5, 2.994089e00, 2.590824e-2, 4.155840e-6, 8.155457e-1, 2.486352e-2, 6.836840e-6, 1.735614e00, 1.763719e-2, 6.370440e-6, 1.968564e00, 1.524020e-2, -2.517152e-6, 1.282417e00, 8.703393e-3, 2.289292e-5, 2.280820e00, 1.918010e-2, 4.484520e-6, 4.833473e-2, 1.641773e-4, -4.654200e-7, 5.589232e-2, -3.455092e-4, -7.388560e-7, 4.634443e-2, -2.658234e-5, 7.757000e-8, 8.997041e-3, 6.329728e-6, -1.939256e-9, 2.284178e-2, -9.941590e-5, 6.787400e-8, 4.350267e-2, -6.839749e-5, -2.714956e-7, 1.348204e-2, 1.091504e-5, 6.903760e-7, 3.106570e-2, -1.665665e-4, -1.590188e-7]) 63 ccsel = np.reshape(ccsel, (17, 3)) 64 65 #Constants of the arguments of the short-period perturbations. 66 dcargs = np.array([5.0974222e0, -7.8604195454652e2, 3.9584962e0, -5.7533848094674e2, 1.6338070e0, -1.1506769618935e3, 2.5487111e0, -3.9302097727326e2, 4.9255514e0, -5.8849265665348e2, 1.3363463e0, -5.5076098609303e2, 1.6072053e0, -5.2237501616674e2, 1.3629480e0, -1.1790629318198e3, 5.5657014e0, -1.0977134971135e3, 5.0708205e0, -1.5774000881978e2, 3.9318944e0, 5.2963464780000e1, 4.8989497e0, 3.9809289073258e1, 1.3097446e0, 7.7540959633708e1, 3.5147141e0, 7.9618578146517e1, 3.5413158e0, -5.4868336758022e2]) 67 dcargs = np.reshape(dcargs, (15, 2)) 68 69 #Amplitudes ccamps(n,k) of the short-period perturbations. 70 ccamps = np.array([-2.279594e-5, 1.407414e-5, 8.273188e-6, 1.340565e-5, -2.490817e-7, -3.494537e-5, 2.860401e-7, 1.289448e-7, 1.627237e-5, -1.823138e-7, 6.593466e-7, 1.322572e-5, 9.258695e-6, -4.674248e-7, -3.646275e-7, 1.140767e-5, -2.049792e-5, -4.747930e-6, -2.638763e-6, -1.245408e-7, 9.516893e-6, -2.748894e-6, -1.319381e-6, -4.549908e-6, -1.864821e-7, 7.310990e-6, -1.924710e-6, -8.772849e-7, -3.334143e-6, -1.745256e-7, -2.603449e-6, 7.359472e-6, 3.168357e-6, 1.119056e-6, -1.655307e-7, -3.228859e-6, 1.308997e-7, 1.013137e-7, 2.403899e-6, -3.736225e-7, 3.442177e-7, 2.671323e-6, 1.832858e-6, -2.394688e-7, -3.478444e-7, 8.702406e-6, -8.421214e-6, -1.372341e-6, -1.455234e-6, -4.998479e-8, -1.488378e-6, -1.251789e-5, 5.226868e-7, -2.049301e-7, 0.e0, -8.043059e-6, -2.991300e-6, 1.473654e-7, -3.154542e-7, 0.e0, 3.699128e-6, -3.316126e-6, 2.901257e-7, 3.407826e-7, 0.e0, 2.550120e-6, -1.241123e-6, 9.901116e-8, 2.210482e-7, 0.e0, -6.351059e-7, 2.341650e-6, 1.061492e-6, 2.878231e-7, 0.e0]) 71 ccamps = np.reshape(ccamps, (15, 5)) 72 73 #Constants csec3 and ccsec(n,k) of the secular perturbations in longitude. 74 ccsec3 = -7.757020e-8 75 ccsec = np.array([1.289600e-6, 5.550147e-1, 2.076942e00, 3.102810e-5, 4.035027e00, 3.525565e-1, 9.124190e-6, 9.990265e-1, 2.622706e00, 9.793240e-7, 5.508259e00, 1.559103e01]) 76 ccsec = np.reshape(ccsec, (4, 3)) 77 78 #Sidereal rates. 79 dcsld = 1.990987e-7 #sidereal rate in longitude 80 ccsgd = 1.990969e-7 #sidereal rate in mean anomaly 81 82 #Constants used in the calculation of the lunar contribution. 83 cckm = 3.122140e-5 84 ccmld = 2.661699e-6 85 ccfdi = 2.399485e-7 86 87 #Constants dcargm(i,k) of the arguments of the perturbations of the motion 88 # of the moon. 89 dcargm = np.array([5.1679830e0, 8.3286911095275e3, 5.4913150e0, -7.2140632838100e3, 5.9598530e0, 1.5542754389685e4]) 90 dcargm = np.reshape(dcargm, (3, 2)) 91 92 #Amplitudes ccampm(n,k) of the perturbations of the moon. 93 ccampm = np.array([1.097594e-1, 2.896773e-7, 5.450474e-2, 1.438491e-7, -2.223581e-2, 5.083103e-8, 1.002548e-2, -2.291823e-8, 1.148966e-2, 5.658888e-8, 8.249439e-3, 4.063015e-8]) 94 ccampm = np.reshape(ccampm, (3, 4)) 95 96 #ccpamv(k)=a*m*dl,dt (planets), dc1mme=1-mass(earth+moon) 97 ccpamv = np.array([8.326827e-11, 1.843484e-11, 1.988712e-12, 1.881276e-12]) 98 dc1mme = 0.99999696e0 99 100 #Time arguments. 101 dt = (dje - dcto) / dcjul 102 tvec = np.array([1e0, dt, dt * dt]) 103 104 #Values of all elements for the instant dje. 105 temp = (np.transpose(np.dot(np.transpose(tvec), np.transpose(dcfel)))) % dc2pi 106 dml = temp[0] 107 forbel = temp[1:8] 108 g = forbel[0] #old fortran equivalence 109 110 deps = (tvec * dceps).sum() % dc2pi 111 sorbel = (np.transpose(np.dot(np.transpose(tvec), np.transpose(ccsel)))) % dc2pi 112 e = sorbel[0] #old fortran equivalence 113 114 #Secular perturbations in longitude. 115 dummy = np.cos(2.0) 116 sn = np.sin((np.transpose(np.dot(np.transpose(tvec[0:2]), np.transpose(ccsec[:,1:3])))) % cc2pi) 117 118 #Periodic perturbations of the emb (earth-moon barycenter). 119 pertl = (ccsec[:,0] * sn).sum() + dt * ccsec3 * sn[2] 120 pertld = 0.0 121 pertr = 0.0 122 pertrd = 0.0 123 for k in range(0, 15): 124 a = (dcargs[k,0] + dt * dcargs[k,1]) % dc2pi 125 cosa = np.cos(a) 126 sina = np.sin(a) 127 pertl = pertl + ccamps[k,0] * cosa + ccamps[k,1] * sina 128 pertr = pertr + ccamps[k,2] * cosa + ccamps[k,3] * sina 129 if k < 11: 130 pertld = pertld + (ccamps[k,1] * cosa - ccamps[k,0] * sina) * ccamps[k,4] 131 pertrd = pertrd + (ccamps[k,3] * cosa - ccamps[k,2] * sina) * ccamps[k,4] 132 133 #Elliptic part of the motion of the emb. 134 phi = (e * e / 4e0) * (((8e0 / e) - e) * np.sin(g) + 5 * np.sin(2 * g) + (13 / 3e0) * e * np.sin(3 * g)) 135 f = g + phi 136 sinf = np.sin(f) 137 cosf = np.cos(f) 138 dpsi = (dc1 - e * e) / (dc1 + e * cosf) 139 phid = 2 * e * ccsgd * ((1 + 1.5 * e * e) * cosf + e * (1.25 - 0.5 * sinf * sinf)) 140 psid = ccsgd * e * sinf / np.sqrt(dc1 - e * e) 141 142 #Perturbed heliocentric motion of the emb. 143 d1pdro = dc1 + pertr 144 drd = d1pdro * (psid + dpsi * pertrd) 145 drld = d1pdro * dpsi * (dcsld + phid + pertld) 146 dtl = (dml + phi + pertl) % dc2pi 147 dsinls = np.sin(dtl) 148 dcosls = np.cos(dtl) 149 dxhd = drd * dcosls - drld * dsinls 150 dyhd = drd * dsinls + drld * dcosls 151 152 #Influence of eccentricity, evection and variation on the geocentric 153 # motion of the moon. 154 pertl = 0.0 155 pertld = 0.0 156 pertp = 0.0 157 pertpd = 0.0 158 for k in range(0, 3): 159 a = (dcargm[k,0] + dt * dcargm[k,1]) % dc2pi 160 sina = np.sin(a) 161 cosa = np.cos(a) 162 pertl = pertl + ccampm[k,0] * sina 163 pertld = pertld + ccampm[k,1] * cosa 164 pertp = pertp + ccampm[k,2] * cosa 165 pertpd = pertpd - ccampm[k,3] * sina 166 167 #Heliocentric motion of the earth. 168 tl = forbel[1] + pertl 169 sinlm = np.sin(tl) 170 coslm = np.cos(tl) 171 sigma = cckm / (1.0 + pertp) 172 a = sigma * (ccmld + pertld) 173 b = sigma * pertpd 174 dxhd = dxhd + a * sinlm + b * coslm 175 dyhd = dyhd - a * coslm + b * sinlm 176 dzhd = -sigma * ccfdi * np.cos(forbel[2]) 177 178 #Barycentric motion of the earth. 179 dxbd = dxhd * dc1mme 180 dybd = dyhd * dc1mme 181 dzbd = dzhd * dc1mme 182 for k in range(0, 4): 183 plon = forbel[k + 3] 184 pomg = sorbel[k + 1] 185 pecc = sorbel[k + 9] 186 tl = (plon + 2.0 * pecc * np.sin(plon - pomg)) % cc2pi 187 dxbd = dxbd + ccpamv[k] * (np.sin(tl) + pecc * np.sin(pomg)) 188 dybd = dybd - ccpamv[k] * (np.cos(tl) + pecc * np.cos(pomg)) 189 dzbd = dzbd - ccpamv[k] * sorbel[k + 13] * np.cos(plon - sorbel[k + 5]) 190 191 #Transition to mean equator of date. 192 dcosep = np.cos(deps) 193 dsinep = np.sin(deps) 194 dyahd = dcosep * dyhd - dsinep * dzhd 195 dzahd = dsinep * dyhd + dcosep * dzhd 196 dyabd = dcosep * dybd - dsinep * dzbd 197 dzabd = dsinep * dybd + dcosep * dzbd 198 199 #Epoch of mean equinox (deq) of zero implies that we should use 200 # Julian ephemeris date (dje) as epoch of mean equinox. 201 if deq == 0: 202 dvelh = au * (np.array([dxhd, dyahd, dzahd])) 203 dvelb = au * (np.array([dxbd, dyabd, dzabd])) 204 return (dvelh,dvelb) 205 206 #General precession from epoch dje to deq. 207 deqdat = (dje - dcto - dcbes) / dctrop + dc1900 208 prema = premat(deqdat, deq, fk4=True) 209 dvelh = au * (np.transpose(np.dot(np.transpose(prema), np.transpose(np.array([dxhd, dyahd, dzahd]))))) 210 dvelb = au * (np.transpose(np.dot(np.transpose(prema), np.transpose(np.array([dxbd, dyabd, dzabd]))))) 211 212 return (dvelh, dvelb)

213

214 -def bprecess(ra0, dec0, mu_radec=None, parallax=None, rad_vel=None, epoch=None):

215 """ 216 Precess positions from J2000.0 (FK5) to B1950.0 (FK4) 217 218 Calculates the mean place of a star at B1950.0 on the FK4 system from the mean 219 place at J2000.0 on the FK5 system. 220 221 Input: 222 @param ra0: J2000 right ascension in degrees 223 @type ra0: N array 224 @param dec0: J2000 declination in degrees 225 @type dec0: N array 226 @keyword mu_radec: array containing the proper motion in seconds of arc per tropical *century* in right ascension and declination. 227 @type mu_radec: 2xN array 228 @keyword parallax: array giving stellar parallax (seconds of arc) 229 @type parallax: N array 230 @keyword rad_vel: array giving radial velocity in km/s 231 @type rad_vel: N array 232 @keyword epoch: scalar giving epoch of original observations, default 2000. This keyword value is only used if mu_radec is not given 233 @type epoch: 234 @return: The corresponding B1950 right ascension and declination in *degrees*. 235 @rtype: N array, N array 236 237 NOTES: 238 The algorithm is taken from the Explanatory Supplement to the Astronomical Almanac 1992, page 186. Also see Aoki et al (1983), A&A, 128,263 239 240 BPRECESS distinguishes between the following two cases: 241 (1) The proper motion is known and non-zero 242 (2) the proper motion is unknown or known to be exactly zero (i.e. 243 extragalactic radio sources). In this case, the reverse of the 244 algorithm in Appendix 2 of Aoki et al. (1983) is used to ensure that the 245 output proper motion is exactly zero. Better precision can be achieved 246 in this case by inputting the epoch of the original observations. 247 248 REVISION HISTORY: 249 Written (W. Landsman, october 1992) 250 Vectorized (W. Landsman, february, 1994) 251 Treat case where proper motion not known or exactly zero (november 1994) 252 Handling of arrays larger than 32767 (Lars L. Christensen, march 1995) 253 Converted to IDL V5.0 (W. Landsman, september 1997) 254 Fixed bug where A term not initialized for vector input (W. Landsman, february 2000) 255 Converted to python (Sergey Koposov, july 2010) 256 """ 257 scal = True 258 if isinstance(ra0, np.ndarray): 259 ra = ra0 260 dec = dec0 261 n = ra.size 262 scal = False 263 else: 264 n = 1 265 ra = np.array([ra0]) 266 dec = np.array([dec0]) 267 268 if rad_vel is None: 269 rad_vel = np.zeros(n) 270 else: 271 if not isinstance(rad_vel, np.ndarray): 272 rad_vel = np.array([rad_vel],dtype=float) 273 if rad_vel.size != n: 274 raise Exception('ERROR - RAD_VEL keyword vector must be of the same length as RA and DEC') 275 276 if (mu_radec is not None): 277 if (np.array(mu_radec).size != 2 * n): 278 raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,' + strtrim(n, 2) + ')') 279 mu_radec = mu_radec * 1. 280 281 if parallax is None: 282 parallax = np.zeros(n) 283 else: 284 if not isinstance(parallax, np.ndarray): 285 parallax = np.array([parallax],dtype=float) 286 287 if epoch is None: 288 epoch = 2000.0e0 289 290 radeg = 180.e0 / np.pi 291 sec_to_radian = lambda x : np.deg2rad(x/3600.) 292 293 m = np.array([np.array([+0.9999256795e0, -0.0111814828e0, -0.0048590040e0, -0.000551e0, -0.238560e0, +0.435730e0]), 294 np.array([+0.0111814828e0, +0.9999374849e0, -0.0000271557e0, +0.238509e0, -0.002667e0, -0.008541e0]), 295 np.array([+0.0048590039e0, -0.0000271771e0, +0.9999881946e0, -0.435614e0, +0.012254e0, +0.002117e0]), 296 np.array([-0.00000242389840e0, +0.00000002710544e0, +0.00000001177742e0, +0.99990432e0, -0.01118145e0, -0.00485852e0]), 297 np.array([-0.00000002710544e0, -0.00000242392702e0, +0.00000000006585e0, +0.01118145e0, +0.99991613e0, -0.00002716e0]), 298 np.array([-0.00000001177742e0, +0.00000000006585e0, -0.00000242404995e0, +0.00485852e0, -0.00002717e0, +0.99996684e0])]) 299 300 a_dot = 1e-3 * np.array([1.244e0, -1.579e0, -0.660e0]) #in arc seconds per century 301 302 ra_rad = np.deg2rad(ra) 303 dec_rad = np.deg2rad(dec) 304 cosra = np.cos(ra_rad) 305 sinra = np.sin(ra_rad) 306 cosdec = np.cos(dec_rad) 307 sindec = np.sin(dec_rad) 308 dec_1950 = dec * 0. 309 ra_1950 = ra * 0. 310 311 for i in range(n): 312 # Following statement moved inside loop in Feb 2000. 313 a = 1e-6 * np.array([-1.62557e0, -0.31919e0, -0.13843e0]) #in radians 314 r0 = np.array([cosra[i] * cosdec[i], sinra[i] * cosdec[i], sindec[i]]) 315 if (mu_radec is not None): 316 mu_a = mu_radec[i,0] 317 mu_d = mu_radec[i,1] 318 r0_dot = np.array([-mu_a * sinra[i] * cosdec[i] - mu_d * cosra[i] * sindec[i], mu_a * cosra[i] * cosdec[i] - mu_d * sinra[i] * sindec[i], mu_d * cosdec[i]]) + 21.095e0 * rad_vel[i] * parallax[i] * r0 319 else: 320 r0_dot = np.array([0.0e0, 0.0e0, 0.0e0]) 321 322 r_0 = np.concatenate((r0, r0_dot)) 323 r_1 = np.transpose(np.dot(np.transpose(m), np.transpose(r_0))) 324 325 # Include the effects of the E-terms of aberration to form r and r_dot. 326 r1 = r_1[0:3] 327 r1_dot = r_1[3:6] 328 if mu_radec is None: 329 r1 = r1 + sec_to_radian ( r1_dot * (epoch - 1950.0e0) / 100. ) 330 a = a + sec_to_radian ( a_dot * (epoch - 1950.0e0) / 100. ) 331 332 x1 = r_1[0] 333 y1 = r_1[1] 334 z1 = r_1[2] 335 rmag = np.sqrt(x1 ** 2 + y1 ** 2 + z1 ** 2) 336 s1 = r1 / rmag 337 s1_dot = r1_dot / rmag 338 s = s1 339 for j in np.arange(0, 3): 340 r = s1 + a - ((s * a).sum()) * s 341 s = r / rmag 342 343 x = r[0] 344 y = r[1] 345 z = r[2] 346 r2 = x ** 2 + y ** 2 + z ** 2 347 rmag = np.sqrt(r2) 348 349 if mu_radec is not None: 350 r_dot = s1_dot + a_dot - ((s * a_dot).sum()) * s 351 x_dot = r_dot[0] 352 y_dot = r_dot[1] 353 z_dot = r_dot[2] 354 mu_radec[i,0] = (x * y_dot - y * x_dot) / (x ** 2 + y ** 2) 355 mu_radec[i,1] = (z_dot * (x ** 2 + y ** 2) - z * (x * x_dot + y * y_dot)) / (r2 * np.sqrt(x ** 2 + y ** 2)) 356 357 dec_1950[i] = np.arcsin(z / rmag) 358 ra_1950[i] = np.arctan2(y, x) 359 360 if parallax[i] > 0.: 361 rad_vel[i] = (x * x_dot + y * y_dot + z * z_dot) / (21.095 * parallax[i] * rmag) 362 parallax[i] = parallax[i] / rmag 363 364 neg = (ra_1950 < 0) 365 if neg.any() > 0: 366 ra_1950[neg] = ra_1950[neg] + 2.e0 * np.pi 367 368 ra_1950 = np.rad2deg(ra_1950) 369 dec_1950 = np.rad2deg(dec_1950) 370 371 # Make output scalar if input was scalar 372 if scal: 373 return ra_1950[0],dec_1950[0] 374 else: 375 return ra_1950, dec_1950

376 377 378 379 ## COLE START HERE

380 -def convolve(image, psf, ft_psf=None, ft_image=None, no_ft=None, correlate=None, auto_correlation=None):

381 """ 382 Convolution of an image with a Point Spread Function (PSF) 383 The default is to compute the convolution using a product of Fourier transforms (for speed). 384 385 Inputs: 386 @param image: 2-D array (matrix) to be convolved with psf 387 @type image: array 388 @param psf: the Point Spread Function, (size < or = to size of image). 389 @type psf array 390 @keyword ft_psf: passes out/in the Fourier transform of the PSF, (so that it can be re-used the next time function is called). 391 @type ft_psc: bool 392 @param ft_image: passes out/in the Fourier transform of image. 393 @type ft_image: bool 394 @param correlate: uses the conjugate of the Fourier transform of PSF to compute the cross-correlation of image and PSF (equivalent to IDL function convol() with NO rotation of PSF) 395 @type correlate: bool 396 @param auto_corr: computes the auto-correlation function of image using FFT. 397 @type auto_corr: bool 398 399 METHOD: 400 When using FFT, PSF is centered & expanded to size of image. 401 402 HISTORY: 403 Written (Frank Varosi, 1992) 404 Appropriate precision type for result depending on input image (Markus Hundertmark, february 2006) 405 Fix the bug causing the recomputation of FFT(psf) and/or FFT(image) (Sergey Koposov, december 2006) 406 """ 407 from numpy.fft import fft2, ifft2 408 n_params = 2 409 psf_ft = ft_psf 410 imft = ft_image 411 noft = no_ft 412 auto = auto_correlation 413 sp = np.array(np.shape(psf_ft)) 414 sif = np.array(np.shape(imft)) 415 sim = np.array(np.shape(image)) 416 sc = sim / 2 417 npix = np.array(image, copy=0).size 418 419 if image.ndim!=2 or noft!=None: 420 if (auto is not None): 421 message("auto-correlation only for images with FFT", inf=True) 422 return image 423 else: 424 if (correlate is not None): 425 return convol(image, psf) 426 else: 427 return convol(image, rotate(psf, 2)) 428 429 if imft==None or (imft.ndim!=2) or imft.shape!=im.shape: #add the type check 430 imft = ifft2(image) 431 if (auto is not None): 432 return np.roll(np.roll(npix * real(fft2(imft * conjugate(imft))), sc[0], 0),sc[1],1) 433 434 if (ft_psf==None or ft_psf.ndim!=2 or ft_psf.shape!=image.shape or ft_psf.dtype!=image.dtype): 435 sp = np.array(shape(psf)) 436 loc = np.maximum((sc - sp / 2), 0) #center PSF in new array, 437 s = np.maximum((sp / 2 - sc), 0) #handle all cases: smaller or bigger 438 l = np.minimum((s + sim - 1), (sp - 1)) 439 psf_ft = np.conjugate(image) * 0 #initialise with correct size+type according to logic of conj and set values to 0 (type of ft_psf is conserved) 440 psf_ft[loc[1]:loc[1]+l[1]-s[1]+1,loc[0]:loc[0]+l[0]-s[0]+1] = psf[s[1]:(l[1])+1,s[0]:(l[0])+1] 441 psf_ft = ifft2(psf_ft) 442 443 if (correlate is not None): 444 conv = npix * np.real(fft2(imft * np.conjugate(psf_ft))) 445 else: 446 conv = npix * np.real(fft2(imft * psf_ft)) 447 448 sc = sc + (sim % 2) #shift correction for odd size images. 449 return np.roll(np.roll(conv, sc[0],0), sc[1],1)

450

451 -def cv_coord(a,b,c,fr=None,to=None,degr=False):

452 #import numpy 453 if degr: 454 degrad = np.deg2rad 455 raddeg = np.rad2deg 456 else: 457 degrad = lambda x: x 458 raddeg = lambda x: x 459 if fr=='sph': 460 cosa = np.cos(degrad(a)) 461 sina = np.sin(degrad(a)) 462 cosb = np.cos(degrad(b)) 463 sinb = np.sin(degrad(b)) 464 x=c*cosa*cosb 465 y=c*sina*cosb 466 z=c*sinb 467 elif fr=='rect': 468 x=a 469 y=b 470 z=c 471 elif fr is None: 472 raise Exception('You must specify the input coordinate system') 473 else: 474 raise Exception('Unknown input coordinate system') 475 if to=='rect': 476 return (x,y,z) 477 elif to=='sph': 478 ra = raddeg(np.arctan2(y,x)) 479 dec = raddeg(np.arctan2(z,np.sqrt(x**2+y**2))) 480 rad = np.sqrt(x**2+y**2+z**2) 481 return (ra,dec,rad) 482 elif to is None: 483 raise Exception('You must specify the output coordinate system') 484 else: 485 raise Exception('Unknown output coordinate system')

486

487 -def daycnv(xjd):

488 """ 489 Converts Julian dates to Gregorian calendar dates 490 491 @param xjd: Julian date, positive double precision scalar or vector 492 @type: array 493 494 OUTPUTS: 495 YR: Year (Integer) 496 MN: Month (Integer) 497 DAY: Day (Integer) 498 HR: Hours and fractional hours (Real). If XJD is a vector, then YR,MN,DAY and HR will be vectors of the same length. 499 500 REVISION HISTORY: 501 Converted to IDL from Yeoman's Comet Ephemeris Generator, B. Pfarr, STX, 6/16/88 502 Converted to IDL V5.0 W. Landsman September 1997 503 """ 504 # Adjustment needed because Julian day starts at noon, calendar day at midnight 505 506 jd = np.array(xjd).astype(int) #Truncate to integral day 507 frac = np.array(xjd).astype(float) - jd + 0.5 #Fractional part of calendar day 508 after_noon = (frac >= 1.0) 509 510 if after_noon.any(): #Is it really the next calendar day? 511 if frac.ndim>0: # proper array 512 frac[after_noon] = frac[after_noon] - 1.0 513 jd[after_noon] = jd[after_noon] + 1 514 else: # scalar 515 frac = frac - 1.0 516 jd = jd + 1 517 hr = frac * 24.0 518 l = jd + 68569 519 n = 4 * l / 146097 520 l = l - (146097 * n + 3) / 4 521 yr = 4000 * (l + 1) / 1461001 522 l = l - 1461 * yr / 4 + 31 #1461 = 365.25 * 4 523 mn = 80 * l / 2447 524 day = l - 2447 * mn / 80 525 l = mn / 11 526 mn = mn + 2 - 12 * l 527 yr = 100. * (n - 49) + yr + l 528 return (yr, mn, day, hr)

529

530 -def euler(ai, bi, select=1, fk4=False):

531 """ 532 Transform between Galactic, celestial, and ecliptic coordinates. 533 534 Input: 535 @param ai: Input longitude in degrees 536 @type ai: float 537 @param bi: Input latitude in degrees 538 @type bi: float 539 @keyword select: Integer (1-6) specifying type of coordinate transformation. 540 @type select: int 1-6 541 @keyword fk4: set the fk4, otherwise use J2000 542 @type fk4: bool 543 544 SELECT From To | SELECT From To 545 1 RA-Dec (2000) Galactic | 4 Ecliptic RA-Dec 546 2 Galactic RA-DEC | 5 Ecliptic Galactic 547 3 RA-Dec Ecliptic | 6 Galactic Ecliptic 548 """ 549 550 twopi = 2.0e0 * np.pi 551 fourpi = 4.0e0 * np.pi 552 553 if fk4: 554 equinox = '(B1950)' 555 psi = np.array ([0.57595865315e0, 4.9261918136e0, 0.00000000000e0, 0.0000000000e0, 0.11129056012e0, 4.7005372834e0]) 556 stheta = np.array ([0.88781538514e0, -0.88781538514e0, 0.39788119938e0, -0.39788119938e0, 0.86766174755e0, -0.86766174755e0]) 557 ctheta = np.array([0.46019978478e0, 0.46019978478e0, 0.91743694670e0, 0.91743694670e0, 0.49715499774e0, 0.49715499774e0]) 558 phi = np.array([4.9261918136e0, 0.57595865315e0, 0.0000000000e0, 0.00000000000e0, 4.7005372834e0, 0.11129056012e0]) 559 else: 560 equinox = '(J2000)' 561 psi = np.array([0.57477043300e0, 4.9368292465e0, 0.00000000000e0, 0.0000000000e0, 0.11142137093e0, 4.71279419371e0]) 562 stheta = np.array([0.88998808748e0, -0.88998808748e0, 0.39777715593e0, -0.39777715593e0, 0.86766622025e0, -0.86766622025e0]) 563 ctheta = np.array([0.45598377618e0, 0.45598377618e0, 0.91748206207e0, 0.91748206207e0, 0.49714719172e0, 0.49714719172e0]) 564 phi = np.array([4.9368292465e0, 0.57477043300e0, 0.0000000000e0, 0.00000000000e0, 4.71279419371e0, 0.11142137093e0]) 565 566 i = select - 1 567 a = np.deg2rad(ai) - phi[i] 568 b = np.deg2rad(bi) 569 sb = np.sin(b) 570 cb = np.cos(b) 571 cbsa = cb * np.sin(a) 572 b = -stheta[i] * cbsa + ctheta[i] * sb 573 bo = np.rad2deg(arcsin(minimum(b, 1.0))) 574 del b 575 a = np.arctan2(ctheta[i] * cbsa + stheta[i] * sb, cb * np.cos(a)) 576 del cb, cbsa, sb 577 ao = np.rad2deg(((a + psi[i] + fourpi) % twopi) ) 578 return (ao, bo)

579

580 -def gal_uvw(distance=None, lsr=None, ra=None, dec=None, pmra=None, pmdec=None, vrad=None, plx=None):

581 """ 582 Calculate the Galactic space velocity (U,V,W) of star 583 584 Calculates the Galactic space velocity U, V, W of star given its 585 (1) coordinates, (2) proper motion, (3) distance (or parallax), and (4) radial velocity. 586 587 OUTPUT PARAMETERS: 588 U - Velocity (km/s) positive toward the Galactic *anti*center 589 V - Velocity (km/s) positive in the direction of Galactic rotation 590 W - Velocity (km/s) positive toward the North Galactic Pole 591 INPUT KEYWORDS: 592 User must supply a position, proper motion,radial velocity and distance 593 (or parallax). Either scalars or vectors can be supplied. 594 (1) Position: 595 RA - Right Ascension in *Degrees* 596 Dec - Declination in *Degrees* 597 (2) Proper Motion 598 PMRA = Proper motion in RA in arc units (typically milli-arcseconds/yr) 599 PMDEC = Proper motion in Declination (typically mas/yr) 600 (3) Radial Velocity 601 VRAD = radial velocity in km/s 602 (4) Distance or Parallax 603 DISTANCE - distance in parsecs 604 or 605 PLX - parallax with same distance units as proper motion measurements 606 typically milliarcseconds (mas) 607 608 REVISION HISTORY: 609 Written, W. Landsman December 2000 610 fix the bug occuring if the input arrays are longer than 32767 611 and update the Sun velocity Sergey Koposov June 2008 612 vectorization of the loop -- performance on large arrays 613 is now 10 times higher Sergey Koposov December 2008 614 """ 615 616 n_params = 3 617 618 if ra is None or dec is None: 619 raise Exception('ERROR - The RA, Dec (J2000) position keywords must be supplied (degrees)') 620 if plx is None and distance is None: 621 raise Exception('ERROR - Either a parallax or distance must be specified') 622 if distance is not None: 623 if np.any(distance == 0): 624 raise Exception('ERROR - All distances must be > 0') 625 plx = 1e3 / distance 626 if plx is not None and any(plx==0): 627 raise Exception('ERROR - Parallaxes must be > 0') 628 629 cosd = np.cos(np.deg2rad(dec)) 630 sind = np.sin(np.deg2rad(dec)) 631 cosa = np.cos(np.deg2rad(ra)) 632 sina = np.sin(np.deg2rad(ra)) 633 k = 4.74047 634 a_g = np.array([[0.0548755604, +0.4941094279, -0.8676661490], [0.8734370902, -0.4448296300, -0.1980763734], [0.4838350155, 0.7469822445, +0.4559837762]]) 635 vec1 = vrad 636 vec2 = k * pmra / plx 637 vec3 = k * pmdec / plx 638 u = (a_g[0,0] * cosa * cosd + a_g[1,0] * sina * cosd + a_g[2,0] * sind) * vec1 + (-a_g[0,0] * sina + a_g[1,0] * cosa) * vec2 + (-a_g[0,0] * cosa * sind - a_g[1,0] * sina * sind + a_g[2,0] * cosd) * vec3 639 v = (a_g[0,1] * cosa * cosd + a_g[1,1] * sina * cosd + a_g[2,1] * sind) * vec1 + (-a_g[0,1] * sina + a_g[1,1] * cosa) * vec2 + (-a_g[0,1] * cosa * sind - a_g[1,1] * sina * sind + a_g[2,1] * cosd) * vec3 640 w = (a_g[0,2] * cosa * cosd + a_g[1,2] * sina * cosd + a_g[2,2] * sind) * vec1 + (-a_g[0,2] * sina + a_g[1,2] * cosa) * vec2 + (-a_g[0,2] * cosa * sind - a_g[1,2] * sina * sind + a_g[2,2] * cosd) * vec3 641 lsr_vel = np.array([-10.00, 5.25, 7.17]) 642 if (lsr is not None): 643 u = u + lsr_vel[0] 644 v = v + lsr_vel[1] 645 w = w + lsr_vel[2] 646 return (u,v,w)

647

648 -def helcorr(ra2000, dec2000, jd, obs_long=None, obs_lat=None, obs_alt=None, debug=False):

649 """ 650 calculates heliocentric Julian date, baricentric and heliocentric radial 651 velocity corrections from: 652 653 INPUT: 654 @param ra2000: Right ascension of object for epoch 2000.0 (hours) 655 @type ra2000: float 656 @param dec2000: declination of object for epoch 2000.0 (degrees) 657 @type dec2000: float 658 @param jd: julian date for the middle of exposure 659 @type jd: float 660 @keyword obs_long: longitude of observatory (degrees, western direction is positive) 661 @type obs_long: float 662 @keyword obs_lat: latitude of observatory (degrees) 663 @type obs_lat: float 664 @keyword obs_alt: altitude of observatory (meters) 665 @type obs_alt: float 666 667 Algorithms used are taken from the IRAF task noao.astutils.rvcorrect 668 and some procedures of the IDL Astrolib are used as well. 669 Accuracy is about 0.5 seconds in time and about 1 m/s in velocity. 670 671 History: 672 written by Peter Mittermayer, Nov 8,2003 673 2005-January-13 Kudryavtsev Made more accurate calculation of the sideral time. 674 Conformity with MIDAS compute/barycorr is checked. 675 2005-June-20 Kochukhov Included precession of RA2000 and DEC2000 to current epoch 676 """ 677 _radeg = 180.0 / np.pi 678 679 # If observatory not given, set La Palma 680 if obs_long is None: 681 obs_long = 17.878333 682 if obs_lat is None: 683 obs_lat = 28.762222 684 if obs_alt is None: 685 obs_alt = 2333. 686 687 #covert JD to Gregorian calendar date 688 xjd = jd*1.0 689 jd = jd-2400000.0 690 year, month, dayraw = convert("JD", "CD", xjd) 691 day = dayraw-dayraw%24. 692 ut = 24.*(dayraw%24) 693 694 #current epoch 695 epoch = year + month / 12. + day / 365. 696 697 #precess ra2000 and dec2000 to current epoch 698 ra,dec=precess(ra2000*15., dec2000, 2000.0, epoch) 699 700 #calculate heliocentric julian date 701 hjd = np.array(helio_jd(jd, ra, dec)).astype(float) 702 703 #DIURNAL VELOCITY (see IRAF task noao.astutil.rvcorrect) 704 dlat = -(11. * 60. + 32.743) * np.sin(2 * obs_lat / _radeg) + 1.1633 * np.sin(4 * obs_lat / _radeg) - 0.0026 * np.sin(6 * obs_lat / _radeg) 705 lat = obs_lat + dlat / 3600 706 707 #calculate distance of observer from earth center 708 r = 6378160.0 * (0.998327073 + 0.001676438 * np.cos(2 * lat / _radeg) - 0.00000351 * np.cos(4 * lat / _radeg) + 0.000000008 * np.cos(6 * lat / _radeg)) + obs_alt 709 710 #calculate rotational velocity (perpendicular to the radius vector) in km/s 711 #23.934469591229 is the siderial day in hours for 1986 712 v = 2. * np.pi * (r / 1000.) / (23.934469591229 * 3600.) 713 714 #calculating local mean siderial time (see astronomical almanach) 715 tu = (jd - 51545.0) / 36525 716 gmst = 6.697374558 + ut + (236.555367908 * (jd - 51545.0) + 0.093104 * tu ** 2 - 6.2e-6 * tu ** 3) / 3600 717 lmst = (gmst - obs_long / 15) % 24 718 719 #projection of rotational velocity along the line of sight 720 vdiurnal = v * np.cos(lat / _radeg) * np.cos(dec / _radeg) * np.sin((ra - lmst * 15) / _radeg) 721 722 #BARICENTRIC and HELIOCENTRIC VELOCITIES 723 vh,vb=baryvel(xjd, 0) 724 725 #project to line of sight 726 vbar = vb[0] * np.cos(dec / _radeg) * np.cos(ra / _radeg) + vb[1] * np.cos(dec / _radeg) * np.sin(ra / _radeg) + vb[2] * np.sin(dec / _radeg) 727 vhel = vh[0] * np.cos(dec / _radeg) * np.cos(ra / _radeg) + vh[1] * np.cos(dec / _radeg) * np.sin(ra / _radeg) + vh[2] * np.sin(dec / _radeg) 728 corr = (vdiurnal + vbar) #using baricentric velocity for correction 729 730 return (corr, hjd)

731

732 -def helio_jd(date, ra, dec, b1950=False, time_diff=False):

733 """ 734 Convert geocentric (reduced) Julian date to heliocentric Julian date 735 736 This procedure correct for the extra light travel time between the Earth and the Sun. 737 An online calculator for this quantity is available at http://www.physics.sfasu.edu/astro/javascript/hjd.html 738 739 Input: 740 @param date: Reduced Julian date (= JD - 2400000) 741 @type date: float 742 @param ra: J2000 right ascension in DEGREES (unless B1950 keyword is set) 743 @type ra: float 744 @param dec: J2000 declination in DEGREES (unless B1950 keyword is set) 745 @type dec: float 746 @keyword b1950: if set, then input coordinates are assumed to be in equinox B1950 coordinates 747 @type b1950: bool 748 @keyword time_diff: if set, then HELIO_JD() returns the time difference (heliocentric JD - geocentric JD ) in seconds 749 @type time_diff: bool 750 @return: heliocentric reduced Julian date. If time_diff is True, then helio_jd() instead returns the time difference in seconds between the geocentric and heliocentric Julian date. 751 """ 752 #Because XYZ uses default B1950 coordinates, we'll convert everything to B1950 753 if not b1950: 754 ra1, dec1 = bprecess(ra, dec) 755 else: 756 ra1 = ra 757 dec1 = dec 758 759 delta_t = (np.array(date).astype(float) - 33282.42345905e0) / 36525.0e0 760 epsilon_sec = np.poly1d([44.836e0, -46.8495, -0.00429, 0.00181][::-1])(delta_t) 761 epsilon = np.deg2rad(23.433333e0 + epsilon_sec / 3600.0e0) 762 ra1 = np.deg2rad(ra1) 763 dec1 = np.deg2rad(dec1) 764 x, y, z, tmp, tmp, tmp = xyz(date) 765 766 #Find extra distance light must travel in AU, multiply by 1.49598e13 cm/AU, 767 #and divide by the speed of light, and multiply by 86400 second/year 768 time = -499.00522e0 * (np.cos(dec1) * np.cos(ra1) * x + (np.tan(epsilon) * np.sin(dec1) + np.cos(dec1) * np.sin(ra1)) * y) 769 if time_diff: 770 return time 771 else: 772 return np.array(date).astype(float) + time / 86400.0e0

773

774 -def precess(ra0, dec0, equinox1, equinox2, doprint=False, fk4=False, radian=False):

775 """ 776 Precess coordinates from EQUINOX1 to EQUINOX2. 777 778 For interactive display, one can use the procedure ASTRO which calls PRECESS 779 or use the /PRINT keyword. The default (RA,DEC) system is FK5 based on epoch 780 J2000.0 but FK4 based on B1950.0 is available via the /FK4 keyword. 781 782 Input: 783 @param ra0: Input right ascension in DEGREES 784 @type ra0: float 785 @param dec0: Input declination in degrees 786 @type dec0: float 787 @param equinox1: first equinox 788 @type equinox1: float 789 @param equinox2: second equinox 790 @type equinox2: float 791 @keyword fk4: If this keyword is set and non-zero, the FK4 (B1950.0) system will be used otherwise FK5 (J2000.0) will be used instead. 792 @type fk4: float 793 @keyword radian: If this keyword is set, input is in radian instead of degrees 794 @type radian 795 796 Restrictions: 797 Accuracy of precession decreases for declination values near 90 degrees. 798 PRECESS should not be used more than 2.5 centuries from 2000 on the FK5 system (1950.0 on the FK4 system). 799 800 Procedure: 801 Algorithm from Computational Spherical Astronomy by Taff (1983), p24. (FK4). 802 FK5 constants from "Astronomical Almanac Explanatory Supplement 1992, page 104 Table 3.211.1. 803 804 Revision history 805 Written, Wayne Landsman, STI Corporation August 1986 806 Correct negative output RA values February 1989 807 Provided FK5 (J2000.0) I. Freedman January 1994 808 Precession Matrix computation now in PREMAT W. Landsman June 1994 809 Work for arrays, not just vectors (W. Landsman, september 2003) 810 Convert to Python (S. Koposov, july 2010) 811 Converted for use at IvS (K. Smolders) 812 """ 813 scal = True 814 if isinstance(ra0, np.ndarray): 815 ra = ra0.copy() 816 dec = dec0.copy() 817 scal = False 818 else: 819 ra = np.array([ra0]) 820 dec = np.array([dec0]) 821 npts = ra.size 822 823 if not radian: 824 ra_rad = np.deg2rad(ra) 825 dec_rad = np.deg2rad(dec) 826 else: 827 ra_rad = ra 828 dec_rad = dec 829 830 a = np.cos(dec_rad) 831 x = np.zeros((npts, 3)) 832 x[:,0] = a * np.cos(ra_rad) 833 x[:,1] = a * np.sin(ra_rad) 834 x[:,2] = np.sin(dec_rad) 835 836 # Use PREMAT function to get precession matrix from Equinox1 to Equinox2 837 r = premat(equinox1, equinox2, fk4=fk4) 838 x2 = np.transpose(np.dot(np.transpose(r), np.transpose(x))) 839 ra_rad = np.zeros(npts) + np.arctan2(x2[:,1], x2[:,0]) 840 dec_rad = np.zeros(npts) + np.arcsin(x2[:,2]) 841 842 if not radian: 843 ra = np.rad2deg(ra_rad) 844 ra = ra + (ra < 0.) * 360.e0 845 dec = np.rad2deg(dec_rad) 846 else: 847 ra = ra_rad 848 dec = dec_rad 849 ra = ra + (ra < 0.) * 2.0e0 * np.pi 850 851 if doprint: 852 print 'Equinox (%.2f): %f,%f' % (equinox2, ra, dec) 853 if scal: 854 ra, dec = ra[0], dec[0] 855 return ra, dec

856 857

858 -def precess_xyz(x, y, z, equinox1, equinox2):

859 """ 860 Precess equatorial geocentric rectangular coordinates. 861 862 Input: 863 @param x: heliocentric rectangular x-coordinates 864 @type x: float 865 @param y: heliocentric rectangular y-coordinates 866 @type y: float 867 @param z: heliocentric rectangular z-coordinates 868 @type z: float 869 @param equinox1: first equinox 870 @type equinox1: float 871 @param equinox2: second equinox 872 @type equinox2: float 873 874 NOTES: 875 The equatorial geocentric rectangular coords are converted to RA and Dec, 876 precessed in the normal way, then changed back to x, y and z using unit 877 vectors. 878 """ 879 #take input coords and convert to ra and dec (in radians) 880 ra = np.arctan2(y, x) 881 _del = np.sqrt(x * x + y * y + z * z) #magnitude of distance to Sun 882 dec = np.arcsin(z / _del) 883 884 # precess the ra and dec 885 ra,dec = precess(ra, dec, equinox1, equinox2, radian=True) 886 887 #convert back to x, y, z 888 xunit = np.cos(ra) * np.cos(dec) 889 yunit = np.sin(ra) * np.cos(dec) 890 zunit = np.sin(dec) 891 x = xunit * _del 892 y = yunit * _del 893 z = zunit * _del 894 return x,y,z

895

896 -def premat(equinox1, equinox2, fk4=False):

897 """ 898 Return the precession matrix needed to go from EQUINOX1 to EQUINOX2. 899 900 This matrix is used by the procedures PRECESS and BARYVEL to precess astronomical coordinates 901 @param equinox1: Original equinox of coordinates, numeric scalar. 902 @type equinox1: float 903 @param equinox2: Equinox of precessed coordinates. 904 @type equinox2: float 905 @return: double precision 3 x 3 precession matrix, used to precess equatorial rectangular coordinates 906 @rtype: 3 by 3 array 907 @keyword fk4: If this keyword is set, the FK4 (B1950.0) system precession angles are used to compute the precession matrix. The default is to use FK5 (J2000.0) precession angles 908 909 Revision history: 910 Written, Wayne Landsman, HSTX Corporation, June 1994 911 Converted to IDL V5.0 W. Landsman September 1997 912 """ 913 deg_to_rad = np.pi / 180.0e0 914 sec_to_rad = deg_to_rad / 3600.e0 915 t = 0.001e0 * (equinox2 - equinox1) 916 917 if not fk4: 918 st = 0.001e0 * (equinox1 - 2000.e0) 919 # Compute 3 rotation angles 920 a = sec_to_rad * t * (23062.181e0 + st * (139.656e0 + 0.0139e0 * st) + t * (30.188e0 - 0.344e0 * st + 17.998e0 * t)) 921 b = sec_to_rad * t * t * (79.280e0 + 0.410e0 * st + 0.205e0 * t) + a 922 c = sec_to_rad * t * (20043.109e0 - st * (85.33e0 + 0.217e0 * st) + t * (-42.665e0 - 0.217e0 * st - 41.833e0 * t)) 923 else: 924 st = 0.001e0 * (equinox1 - 1900.e0) 925 # Compute 3 rotation angles 926 a = sec_to_rad * t * (23042.53e0 + st * (139.75e0 + 0.06e0 * st) + t * (30.23e0 - 0.27e0 * st + 18.0e0 * t)) 927 b = sec_to_rad * t * t * (79.27e0 + 0.66e0 * st + 0.32e0 * t) + a 928 c = sec_to_rad * t * (20046.85e0 - st * (85.33e0 + 0.37e0 * st) + t * (-42.67e0 - 0.37e0 * st - 41.8e0 * t)) 929 930 sina = np.sin(a) 931 sinb = np.sin(b) 932 sinc = np.sin(c) 933 cosa = np.cos(a) 934 cosb = np.cos(b) 935 cosc = np.cos(c) 936 r = np.zeros((3, 3)) 937 r[0,:] = np.array([cosa * cosb * cosc - sina * sinb, sina * cosb + cosa * sinb * cosc, cosa * sinc]) 938 r[1,:] = np.array([-cosa * sinb - sina * cosb * cosc, cosa * cosb - sina * sinb * cosc, -sina * sinc]) 939 r[2,:] = np.array([-cosb * sinc, -sinb * sinc, cosc]) 940 941 return r

942 943

944 -def sphdist (ra1, dec1, ra2, dec2):

945 """ 946 Measures the spherical distance in degrees. The input has to be in degrees too 947 """ 948 dec1_r = np.deg2rad(dec1) 949 dec2_r = no.deg2rad(dec2) 950 return 2 * np.rad2deg(np.arcsin(np.sqrt((np.sin((dec1_r - dec2_r) / 2))**2+np.cos(dec1_r) * np.cos(dec2_r) *(np.sin((np.deg2rad(ra1 - ra2)) / 2))**2)))

951

952 -def xyz(date, equinox=None):

953 """ 954 Calculate geocentric X,Y, and Z and velocity coordinates of the Sun 955 956 Calculates geocentric X,Y, and Z vectors and velocity coordinates 957 (dx, dy and dz) of the Sun. (The positive X axis is directed towards the 958 equinox, the y-axis, towards the point on the equator at right ascension 6h, 959 and the z axis toward the north pole of the equator). Typical position 960 accuracy is <1e-4 AU (15000 km). 961 962 Input: 963 @param date: reduced julian date (=JD - 2400000) 964 @type date: float 965 @keyword equinox: equinox of output. Default is 1950. 966 @type equinox: float 967 968 Output: 969 x,y,z: scalars or vectors giving heliocentric rectangular coordinates 970 (in A.U) for each date supplied. Note that sqrt(x^2 + y^2 + z^2) 971 gives the Earth-Sun distance for the given date. 972 xvel, yvel, zvel: velocity vectors corresponding to X, Y and Z. 973 """ 974 picon = np.pi / 180.0e0 975 t = (date - 15020.0e0) / 36525.0e0 #Relative Julian century from 1900 976 977 # NOTE: longitude arguments below are given in *equinox* of date. 978 # Precess these to equinox 1950 to give everything an even footing. 979 # Compute argument of precession from equinox of date back to 1950 980 pp = (1.396041e0 + 0.000308e0 * (t + 0.5e0)) * (t - 0.499998e0) 981 982 # Compute mean solar longitude, precessed back to 1950 983 el = 279.696678e0 + 36000.76892e0 * t + 0.000303e0 * t * t - pp 984 985 # Compute Mean longitude of the Moon 986 c = 270.434164e0 + 480960.e0 * t + 307.883142e0 * t - 0.001133e0 * t * t - pp 987 988 # Compute longitude of Moon's ascending node 989 n = 259.183275e0 - 1800.e0 * t - 134.142008e0 * t + 0.002078e0 * t * t - pp 990 991 # Compute mean solar anomaly 992 g = 358.475833e0 + 35999.04975e0 * t - 0.00015e0 * t * t 993 994 # Compute the mean jupiter anomaly 995 j = 225.444651e0 + 2880.0e0 * t + 154.906654e0 * t * t 996 997 # Compute mean anomaly of Venus 998 v = 212.603219e0 + 58320.e0 * t + 197.803875e0 * t + 0.001286e0 * t * t 999 1000 # Compute mean anomaly of Mars 1001 m = 319.529425e0 + 19080.e0 * t + 59.8585e0 * t + 0.000181e0 * t * t 1002 1003 # Convert degrees to radians for trig functions 1004 el = el * picon 1005 g = g * picon 1006 j = j * picon 1007 c = c * picon 1008 v = v * picon 1009 n = n * picon 1010 m = m * picon 1011 1012 # Calculate X,Y,Z using trigonometric series 1013 x = 0.999860e0 * np.cos(el) - 0.025127e0 * np.cos(g - el) + 0.008374e0 * np.cos(g + el) + 0.000105e0 * np.cos(g + g + el) + 0.000063e0 * t * np.cos(g - el) + 0.000035e0 * np.cos(g + g - el) - 0.000026e0 * np.sin(g - el - j) - 0.000021e0 * t * np.cos(g + el) + 0.000018e0 * np.sin(2.e0 * g + el - 2.e0 * v) + 0.000017e0 * np.cos(c) - 0.000014e0 * np.cos(c - 2.e0 * el) + 0.000012e0 * np.cos(4.e0 * g + el - 8.e0 * m + 3.e0 * j) - 0.000012e0 * np.cos(4.e0 * g - el - 8.e0 * m + 3.e0 * j) - 0.000012e0 * np.cos(g + el - v) + 0.000011e0 * np.cos(2.e0 * g + el - 2.e0 * v) + 0.000011e0 * np.cos(2.e0 * g - el - 2.e0 * j) 1014 y = 0.917308e0 * np.sin(el) + 0.023053e0 * np.sin(g - el) + 0.007683e0 * np.sin(g + el) + 0.000097e0 * np.sin(g + g + el) - 0.000057e0 * t * np.sin(g - el) - 0.000032e0 * np.sin(g + g - el) - 0.000024e0 * np.cos(g - el - j) - 0.000019e0 * t * np.sin(g + el) - 0.000017e0 * np.cos(2.e0 * g + el - 2.e0 * v) + 0.000016e0 * np.sin(c) + 0.000013e0 * np.sin(c - 2.e0 * el) + 0.000011e0 * np.sin(4.e0 * g + el - 8.e0 * m + 3.e0 * j) + 0.000011e0 * np.sin(4.e0 * g - el - 8.e0 * m + 3.e0 * j) - 0.000011e0 * np.sin(g + el - v) + 0.000010e0 * np.sin(2.e0 * g + el - 2.e0 * v) - 0.000010e0 * np.sin(2.e0 * g - el - 2.e0 * j) 1015 z = 0.397825e0 * np.sin(el) + 0.009998e0 * np.sin(g - el) + 0.003332e0 * np.sin(g + el) + 0.000042e0 * np.sin(g + g + el) - 0.000025e0 * t * np.sin(g - el) - 0.000014e0 * np.sin(g + g - el) - 0.000010e0 * np.cos(g - el - j) 1016 1017 #Precess_to new equator? 1018 if equinox is not None: 1019 x, y, z = precess_xyz(x, y, z, 1950, equinox) 1020 1021 xvel = -0.017200e0 * np.sin(el) - 0.000288e0 * np.sin(g + el) - 0.000005e0 * np.sin(2.e0 * g + el) - 0.000004e0 * np.sin(c) + 0.000003e0 * np.sin(c - 2.e0 * el) + 0.000001e0 * t * np.sin(g + el) - 0.000001e0 * np.sin(2.e0 * g - el) 1022 yvel = 0.015780 * np.cos(el) + 0.000264 * np.cos(g + el) + 0.000005 * np.cos(2.e0 * g + el) + 0.000004 * np.cos(c) + 0.000003 * np.cos(c - 2.e0 * el) - 0.000001 * t * np.cos(g + el) 1023 zvel = 0.006843 * np.cos(el) + 0.000115 * np.cos(g + el) + 0.000002 * np.cos(2.e0 * g + el) + 0.000002 * np.cos(c) + 0.000001 * np.cos(c - 2.e0 * el) 1024 1025 #Precess to new equator? 1026 if equinox is not None: 1027 xvel, yvel, zvel = precess_xyz(xvel, yvel, zvel, 1950, equinox) 1028 1029 return x, y, z, xvel, yvel, zvel

1030

Source Code for Module ivs.observations.barycentric_correction