ivs.units.conversions

1 # -*- coding: utf-8 -*- 2 """ 3 Convert one unit (and uncertainty) to another. 4 5 Contents: 6 7 1. B{The Python module}: basic usage of the Python module 8 2. B{The Terminal tool}: basic usage of the terminal tool 9 3. B{Fundamental constants and base unit systems} 10 - Changing base unit system 11 - Changing the values of the fundamental constants 12 4. B{Calculating with units} 13 5. B{Dealing with and interpreting unit strings} 14 15 Much of the unit conversions has been tested using Appendix B of 16 U{http://physics.nist.gov/cuu/pdf/sp811.pdf}, and examples within, which is the 17 international document describing the SI standard. As a matter of fact, this 18 module should be fully compatible with the SI unit conventions, except for the 19 notation of units (brackets are not allowed). This module extends the SI unit 20 conventions to be more flexible and intuitive, and also allows to ditch the SI 21 unit convention alltogether. 22 23 24 25 Some of the many possibilities include (see L{convert} for an extensive set of 26 examples): 27 28 1. Conversions between equal-type units: meter to nano-lightyears, erg/s 29 to W, cy/d to muHz, but also erg/s/cm2/A to W/m2/mum, sr to deg2, etc... 30 2. Conversions between unequal-type units: angstrom to km/s via the speed 31 of light, F(lambda) to F(nu), F(nu) to lambdaF(lambda)/sr, meter to 32 cycles/arcsec (interferometry), etc... 33 3. Nonlinear conversions: vegamag to erg/s/cm2/A or Jy, Celcius to 34 Fahrenheit or Kelvin, calender date to (any kind of modified) Julian Day, ... 35 4. Conversions of magnitude to flux amplitudes via 'Amag' and 'ppt' or 'ampl' 36 5. Conversions of magnitude colors to flux ratios via 'mag_color' and 'flux_ratio' 37 6. Coordinate transformations (equatorial to galactic etc.) 38 7. Logarithmic conversions, e.g. from logTeff to Teff via '[K]' and 'K' 39 8. Inclusion of uncertainties, both in input values and/or reference values 40 when converting between unequal-type units, automatically recognised when 41 giving two positional argument (value, error) instead of one (value). 42 9. Currency exchange rates. First, you need to call L{set_exchange_rates} 43 for this to work, since the latest currency definitions and rates need to 44 queried from the European Central Bank (automatically done when using the 45 terminal tool). 46 10. Computations with units. 47 48 49 B{Warning 1:} frequency units are technically given in cycles per time (cy). 50 This means that if you want to convert e.g. muHz to d-1 (or 1/d), you need to 51 ask for cy/d. There is a general ambiguity concerning the unit 'cycles': it is 52 not an official SI unit, but is needed to make the distinction between e.g. rad/s 53 and Hz. In true SI-style, the official "unit" of rad/s is actually just s-1, since 54 radians are not really a unit. This gives room for confusion! To make everything 55 even more confusing, there is another unit which is equal to the reciprocal 56 second, namely the Becquerel. This is used for stochastic or non-recurrent 57 phenomena. Basically, the problem is:: 58 59 rad/s == 1/s 60 1/s == Hz 61 62 but:: 63 64 rad/s != Hz 65 66 67 If you have any doubts, set the logger to display debug messages (e.g. 68 C{logger.setLevel(10)}, see L{ivs.aux.loggers}). It will sometimes tell you if 69 certain assumptions are made. 70 71 72 B{Warning 2:} there is some ambiguity between units. For example, C{as} can be 73 interpreted as 'arcsecond', but also as 'attosecond'. In case of ambiguity, 74 the basic units will be preferred (in this case 'arcsecond'), instead of the one 75 with prefixes. This is a list of non-exhaustive ambiguities (the preffered one 76 in italic): 77 78 - C{as}: I{arcsecond} vs attosecond 79 - C{am}: I{arcminute} vs attominute 80 - C{min}: I{minute} vs milli-inch 81 - C{yd}: I{yard} vs yoctoday 82 - C{ch}: I{chain} vs centihour 83 84 Sometimes the case-sensitivity of the conversions comes in handy. There is no 85 ambiguity between: 86 87 - C{foe} (fifty-one-ergs) and C{fOe} (femto Oersted) 88 89 B{Warning 3:} the unit name of angstrom is AA, ampere is A. 90 91 B{Warning 4:} Most of the imperial units are UK/Canada. If you need US, prefix 92 the unit with C{US}: E.g. The gallon (C{gal}) is the international imperial 93 gallon, C{USgal} is the US gallon. 94 95 B{Note 1:} Photometric passbands are given as a string C{"SYSTEM.FILTER"}. For 96 a list of defined systems and filters, see L{ivs.sed.filters}. 97 98 Section 1. The Python module 99 ============================ 100 101 The main function L{convert} (see link for a full list of examples) does all the 102 work and is called via 103 104 >>> result = convert('km','m',1.) 105 106 or when the error is known 107 108 >>> result,e_result = convert('km','m',1.,0.1) 109 110 Be B{careful} when mixing nonlinear conversions (e.g. magnitude to flux) with 111 linear conversions (e.g. Jy to W/m2/m). 112 113 When your favorite conversion is not implemented, there are five places 114 where you can add information: 115 116 1. C{_scalings}: if your favorite prefix (e.g. Tera, nano...) is not 117 available 118 2. C{_aliases}: if your unit is available but not under the name you are 119 used to. 120 3. C{_factors}: if your unit is not available. 121 4. C{_switch}: if your units are available, but conversions from one to 122 another is not straightforward and extra infromation is needed (e.g. to go 123 from angstrom to km/s in a spectrum, a reference wavelength is needed). 124 5. C{_convention}: if your favorite base unit system is not defined (SI,cgs...). 125 126 If you need to add a linear factor, just give the factor in SI units, and the 127 SI base units it consists of (see C{_factors} for examples). If you need to add 128 a nonlinear factor, you need to give a class definition (see the examples). 129 130 Section 2. The Terminal tool 131 ============================ 132 133 For help and a list of all defined units and abbreviations, do:: 134 135 $:> python conversions.py --help 136 ===================================== | ===================================== 137 = Units of absorbed dose = | = Units of acceleration = 138 ===================================== | ===================================== 139 Gy = gray | Gal = Gal 140 ===================================== | ===================================== 141 = Units of angle = | = Units of area = 142 ===================================== | ===================================== 143 am = arcminute | a = are 144 as = arcsecond | ac = acre (international) 145 cy = cycle | b = barn 146 deg = degree | ===================================== 147 rad = radian | = Units of coordinate = 148 rpm = revolutions per minute | ===================================== 149 sr = sterradian | complex_coord = <own unit> 150 ===================================== | deg_coord = degrees 151 = Units of catalytic activity = | ecliptic = ecliptic 152 ===================================== | equatorial = equatorial 153 kat = katal | galactic = galactic 154 ===================================== | rad_coord = radians 155 = Units of currency = | ===================================== 156 ===================================== | = Units of dose equivalent = 157 EUR = EURO | ===================================== 158 ===================================== | Sv = sievert 159 = Units of dynamic viscosity = | rem = rem 160 ===================================== | ===================================== 161 P = poise | = Units of electric capacitance = 162 ===================================== | ===================================== 163 = Units of electric charge = | F = Farad 164 ===================================== | ===================================== 165 C = Coulomb | = Units of electric conductance = 166 ===================================== | ===================================== 167 = Units of electric current = | S = Siemens 168 ===================================== | ===================================== 169 A = Ampere | = Units of electric potential difference = 170 Bi = biot | ===================================== 171 Gi = Ampere | V = Volt 172 ===================================== | ===================================== 173 = Units of electric resistance = | = Units of energy = 174 ===================================== | ===================================== 175 O = Ohm | Cal = large calorie (international table) 176 ===================================== | J = Joule 177 = Units of energy/power = | cal = small calorie (international table) 178 ===================================== | eV = electron volt 179 Lsol = Solar luminosity | erg = ergon 180 ===================================== | foe = (ten to the) fifty one ergs 181 = Units of flux density = | ===================================== 182 ===================================== | = Units of flux = 183 Jy = Jansky | ===================================== 184 ===================================== | ABmag = AB magnitude 185 = Units of frequency = | Amag = amplitude in magnitude 186 ===================================== | STmag = ST magnitude 187 hz = Hertz | ampl = fractional amplitude 188 ===================================== | flux_ratio = flux ratio 189 = Units of inductance = | mag = magnitude 190 ===================================== | mag_color = color 191 H = Henry | pph = amplitude in parts per hundred 192 ===================================== | ppm = amplitude in parts per million 193 = Units of length = | ppt = amplitude in parts per thousand 194 ===================================== | vegamag = Vega magnitude 195 AA = angstrom | ===================================== 196 AU = astronomical unit | = Units of force = 197 Rearth = Earth radius | ===================================== 198 Rjup = Jupiter radius | N = Newton 199 Rsol = Solar radius | dyn = dyne 200 USft = foot (US) | ===================================== 201 USmi = mile (US) | = Units of illuminance = 202 a0 = Bohr radius | ===================================== 203 bs = beard second | lx = lux 204 ch = chain | ph = phot 205 ell = ell | sb = stilb 206 fathom = fathom | ===================================== 207 ft = foot (international) | = Units of kynamatic viscosity = 208 fur = furlong | ===================================== 209 in = inch (international) | St = stokes 210 ly = light year | ===================================== 211 m = meter | = Units of luminous flux = 212 mi = mile (international) | ===================================== 213 nami = nautical mile | lm = lumen 214 pc = parsec | ===================================== 215 perch = pole | = Units of magnetic field strength = 216 pole = perch | ===================================== 217 potrzebie = potrzebie | G = Gauss 218 rd = rod | T = Tesla 219 smoot = smooth | ===================================== 220 yd = yard (international) | = Units of magnetizing field = 221 ===================================== | ===================================== 222 = Units of m/s = | Oe = Oersted 223 ===================================== | ===================================== 224 cc = Speed of light | = Units of power = 225 ===================================== | ===================================== 226 = Units of magnetic flux = | W = Watt 227 ===================================== | dp = Donkeypower 228 Mx = Maxwell | hp = Horsepower 229 Wb = Weber | ===================================== 230 ===================================== | = Units of roentgen = 231 = Units of mass = | ===================================== 232 ===================================== | R = Roentgen 233 Mearth = Earth mass | ===================================== 234 Mjup = Jupiter mass | = Units of temperature = 235 Mlun = Lunar mass | ===================================== 236 Msol = Solar mass | Cel = Celcius 237 carat = carat | Far = Fahrenheit 238 firkin = firkin | K = Kelvin 239 g = gram | Tsol = Solar temperature 240 gr = gram | ===================================== 241 lb = pound | = Units of velocity = 242 mol = molar mass | ===================================== 243 ounce = ounce | knot = nautical mile per hour 244 st = stone | mph = miles per hour 245 ton = gram | 246 u = atomic mass | 247 ===================================== | 248 = Units of pressure = | 249 ===================================== | 250 Pa = Pascal | 251 at = atmosphere (technical) | 252 atm = atmosphere (standard) | 253 ba = barye | 254 bar = baros | 255 mmHg = millimeter of mercury | 256 psi = pound per square inch | 257 torr = Torricelli | 258 ===================================== | 259 = Units of second = | 260 ===================================== | 261 fortnight = fortnight | 262 ===================================== | 263 = Units of time = | 264 ===================================== | 265 Bq = Becquerel | 266 CD = calender day | 267 Ci = Curie | 268 JD = Julian day | 269 MJD = modified Julian day | 270 cr = century | 271 d = day | 272 h = hour | 273 j = jiffy | 274 min = minute | 275 mo = month | 276 s = second | 277 sidereal = sidereal day | 278 wk = week | 279 yr = year | 280 ===================================== | 281 = Units of volume = | 282 ===================================== | 283 USgal = gallon (US) | 284 USgi = gill (Canadian and UK imperial)| 285 bbl = barrel | 286 bu = bushel | 287 gal = gallon (Canadian and UK imperial)| 288 gi = gill (Canadian and UK imperial)| 289 l = liter | 290 ngogn = 1000 cubic potrzebies | 291 292 293 Usage: conversions.py --from=<unit> --to=<unit> [options] value [error] 294 295 Options: 296 -h, --help show this help message and exit 297 --from=_FROM units to convert from 298 --to=_TO units to convert to 299 300 Extra quantities when changing base units (e.g. Flambda to Fnu): 301 -w WAVE, --wave=WAVE 302 wavelength with units (e.g. used to convert Flambda to 303 Fnu) 304 -f FREQ, --freq=FREQ 305 frequency (e.g. used to convert Flambda to Fnu) 306 -p PHOTBAND, --photband=PHOTBAND 307 photometric passband 308 309 To convert from one unit to another, do:: 310 311 $:> python conversions.py --from=nRsol/h --to cm/s=1.2345 312 1.2345 nRsol/h = 0.0238501 cm/s 313 314 In fact, the C{to} parameter is optional (despite it not being a positional 315 argument, C{from} is not optional). The script will assume you want to convert 316 to SI:: 317 318 $:> python conversions.py --from=nRsol/h 1.2345 319 1.2345 nRsol/h = 0.000238501 m1 s-1 320 321 It is also possible to compute with uncertainties, and you can give extra 322 keyword arguments if extra information is needed for conversions:: 323 324 $:> python conversions.py --from=mag --to=erg/s/cm2/AA --photband=GENEVA.U 7.84 0.02 325 7.84 +/- 0.02 mag = 4.12191e-12 +/- 7.59283e-14 erg/s/cm2/AA 326 327 If you want to do coordinate transformations, e.g. from fractional radians to 328 degrees/arcminutes/arcseconds, you can do:: 329 330 $:> python conversions.py --from=rad_coord --to=equatorial 5.412303,0.123 331 (5.412303, 0.123) rad_coord = 20:40:24.51,7:02:50.6 equ 332 333 Section 3. Fundamental constants and base unit systems 334 ====================================================== 335 336 Although fundamental constants are supposed to be constants, there are two 337 reasons why one might to change their value: 338 339 1. You want to work in a different base unit system (e.g. cgs instead of SI) 340 2. You disagree with some of the literature values and want to use your own. 341 342 Section 3.1. Changing base unit system 343 -------------------------------------- 344 345 A solution to the first option might be to define all constants in SI, and define 346 them again in cgs, adding a postfix C{_cgs} to the constants' name. This is done 347 in the module L{ivs.units.constants} to give the user access to values of constants 348 in common-used systems. However, this way of working does not offer any flexibility, 349 and one has to redefine all the values manually for any given convention. Also, you 350 may want to work in cgs by default, which means that typing the postfix every 351 time is cumbersome. For this purpose, you can redefine all constants in a 352 different base system with a simple command: 353 354 >>> print constants.Msol,constants.GG 355 1.988547e+30 6.67384e-11 356 >>> set_convention(units='cgs') 357 ('SI', 'standard', 'rad') 358 >>> print constants.Msol,constants.GG 359 1.988547e+33 6.67384e-08 360 361 Remark that the units are also changed accordingly: 362 363 >>> print constants.Msol_units 364 g1 365 366 You can go crazy with this: 367 368 >>> set_convention(units='imperial') 369 ('cgs', 'standard', 'rad') 370 >>> print(constants.GG,constants.GG_units) 371 (3.9594319112470405e-11, 'yd3 lb-1 s-2') 372 373 Resetting to the default SI can be done by calling L{set_convention} without any 374 arguments, or simply 375 376 >>> set_convention(units='SI') 377 ('imperial', 'standard', 'rad') 378 379 The function L{set_convention} returns the current settings, so you can 380 remember the old settings and make a temporary switch: 381 382 >>> old_settings = set_convention(units='cgs') 383 >>> set_convention(*old_settings) 384 ('cgs', 'standard', 'rad') 385 386 B{Warning:} Changing the value of the constants affects B{all} modules, where 387 an import statement C{from ivs.units import constants} is made. It will B{not} 388 change the values in modules where the import is done via 389 C{from ivs.units.constants import *}. If you want to get the value of a 390 fundamental constant regardless of the preference of base system, you might want 391 to call 392 393 >>> Msol = get_constant('Msol','cgs') 394 395 To change the behaviour of the interpretation of radians and cycle, you can 396 specify the keyword C{frequency}. See L{set_convention} for more information. 397 398 Section 3.2. Changing the values of the fundamental constants 399 ------------------------------------------------------------- 400 401 If you disagree with some of the literature values, you can also redefine the 402 values of the fundamental constants. For example, to use the values defined in 403 the stellar evolution code B{MESA}, simply do 404 405 >>> set_convention(units='cgs',values='mesa') 406 ('SI', 'standard', 'rad') 407 408 You can query for info on the current convention without changing it: 409 410 >>> get_convention() 411 ('cgs', 'mesa', 'rad') 412 413 But for the sake of the examples, we'll switch back to the default SI... 414 415 >>> reset() 416 417 418 Section 4. Calculating with units 419 ================================= 420 421 There exists a class L{Unit}, which allows you to calculate with values with 422 units (optionally including uncertainties). The class L{Unit} offers high 423 flexibility when it comes to initialization, so that a user can create a L{Unit} 424 according to his/her own preferences. See L{Unit.__init__} for all the options. 425 426 The purpose of existence of the class L{Unit} is calculation: you can simply 427 create a couple of units, and multiply, add, divide etc.. with Python's standard 428 operators ('*','+','/',...). The standard rules apply: you can multiply and divide 429 basically any combination of parameters, but adding and subtracting requires Units 430 to have the same dimensions! 431 432 Aside from the basic operators, a L{Unit} also understands the most commonly 433 used operations from numpy (C{np.sin, np.cos, np.sqrt}...). If a function is 434 not implemented, you can easily add it yourself following the existing examples 435 in the code. A L{Unit} does B{not} understand the functions from the C{math} 436 module, only the C{numpy} implementation! 437 438 Better yet, you can also put L{Unit}s in a numpy array. Just make a list of units 439 and make it an array as you would do with normal floats (no need to specify the 440 C{dtype}). You can safely mix units with and without uncertainties, and mix units 441 with different dimensions. You have access to most of numpy functions such as 442 C{mean}, C{sum} etc. Of course, when you mix different dimensions, calling C{sum} 443 will result in an error! 444 445 B{Example 1:} Calculate the inclination of star, given a measured 446 C{vsini=150+/-11} km/s, a rotation period of C{2+/-1} days and a radius of 447 C{12} solar radii (no errors). First we give in our data: 448 449 >>> vsini = Unit((150.,11),'km/s') 450 >>> P = Unit((2.,1.),'d') 451 >>> R = Unit(12.,'Rsol') 452 453 Now we can calculate the C{sini}: 454 455 >>> sini = vsini * P / (2*np.pi*R) 456 457 And take the arcsine to recover the inclination angle in radians. Because we 458 are working with the L{Unit} class, we can also immediately retrieve it in 459 degrees: 460 461 >>> print np.arcsin(sini) 462 0.517004704077+/-0.287337185083 rad 463 >>> print np.arcsin(sini).convert('deg') 464 29.622187532+/-16.4632080024 deg 465 466 467 B{Example 2:} The following is an exam question on the Biophysics exam (1st year 468 Bachelor) about error propagation. 469 470 Question: Suppose there is a party to celebrate the end of the Biophysics exam. 471 You want to invite 4 persons, and you want to know if 1 liter of champagne is 472 enough to fill 5 glasses. The glasses are cylinders with a circumference of 473 15.5+/-0.5cm, and a height of 10.0+/-0.5cm. Calculate the volume of one glass 474 and its uncertainty. Can you pour 5 glasses of champagne from the 1 liter 475 bottle? 476 477 Answer: 478 479 >>> r = Unit( (15.5/(2*np.pi), 0.5/(2*np.pi)), 'cm') 480 >>> h = Unit( (10.0,0.5), 'cm') 481 >>> V = np.pi*r**2*h 482 >>> print(V) 483 0.000191184875389+/-1.56051027314e-05 m3 484 >>> print((5*V).convert('liter')) 485 0.955924376946+/-0.0780255136569 liter 486 487 It is not sufficient within about 1 sigma. 488 489 Section 5, Dealing with and interpreting string units 490 ===================================================== 491 492 Suppose you want to make automated plots with a human readable name and a LaTeX 493 string of the unit. Then you can use L{get_type} and L{unit2texlabel}: 494 495 >>> print(unit2texlabel('erg/s/mum2/AA',full=True)) 496 Flux Density [erg s$^{-1}$ $\mu$m$^{-2}$ $\AA$$^{-1}$] 497 >>> print(unit2texlabel('P',full=True)) 498 Dynamic Viscosity [P] 499 500 or 501 502 >>> print(unit2texlabel('erg/s/mum2/AA')) 503 erg s$^{-1}$ $\mu$m$^{-2}$ $\AA$$^{-1}$ 504 505 506 """ 507 #-- standard libraries 508 import itertools 509 import functools 510 import collections 511 import inspect 512 import re 513 import os 514 import sys 515 import logging 516 import urllib 517 import numpy as np 518 import pytz 519 import datetime 520 521 #-- optional libraries: WARNING: when these modules are not installed, the 522 # module's use is restricted 523 try: import ephem 524 except ImportError: print("Unable to load pyephem, stellar coordinate transformations unavailable") 525 526 #-- from IVS repository 527 from ivs.units import constants 528 from ivs.units.uncertainties import unumpy,AffineScalarFunc,ufloat 529 from ivs.units.uncertainties.unumpy import log10,log,exp,sqrt 530 from ivs.units.uncertainties.unumpy import sin,cos,tan 531 from ivs.units.uncertainties.unumpy import arcsin,arccos,arctan 532 from ivs.sed import filters 533 from ivs.inout import ascii 534 from ivs.aux import loggers 535 from ivs.aux.decorators import memoized 536 537 logger = logging.getLogger("UNITS.CONV") 538 logger.addHandler(loggers.NullHandler())

539 540 #{ Main functions 541 542 -def convert(_from,_to,*args,**kwargs):

543 """ 544 Convert one unit to another. 545 546 Basic explanation 547 ================= 548 549 The unit strings C{_from} and C{_to} should by default be given in the form 550 551 C{erg s-1 cm-2 AA-1} 552 553 Common alternatives are also accepted (see below). 554 555 Square brackets '[]' denote a logarithmic value. 556 557 If one positional argument is given, it can be either a scalar, numpy array 558 or C{uncertainties} object. The function will also return one argument of 559 the same type. 560 561 If two positional arguments are given, the second argument is assumed to be 562 the uncertainty on the first (scalar or numpy array). The function will also 563 return two arguments. 564 565 Basic examples: 566 567 >>> convert('km','cm',1.) 568 100000.0 569 >>> convert('m/s','km/h',1,0.1) 570 (3.5999999999999996, 0.36) 571 572 Keyword arguments can give extra information, for example when converting 573 from Flambda to Fnu, and should be tuples (float(,error),'unit'): 574 575 >>> convert('AA','km/s',4553,0.1,wave=(4552.,0.1,'AA')) 576 (65.85950307557613, 9.314963362464114) 577 578 Extra 579 ===== 580 581 The unit strings C{_from} and C{_to} should by default be given in the form 582 583 C{erg s-1 cm-2 AA-1} 584 585 Common alternatives are also accepted, but don't drive this too far: 586 587 C{erg/s/cm2/AA} 588 589 The crasiest you're allowed to go is something like 590 591 >>> print(convert('10mW m-2/nm','erg s-1 cm-2 AA-1',1.)) 592 1.0 593 594 But there is a limit on the interpretation of this prefactor also. Floats 595 will probably not work, and exponentials require exactly two digits. 596 597 Parentheses are in no circumstances accepted. Some common aliases are also 598 resolved (for a full list, see dictionary C{_aliases}): 599 600 C{erg/s/cm2/angstrom} 601 602 You don't really have to spell both units if either the 'from' or 'to' units 603 is consistently within one convention (SI, cgs, solar...). But of course you 604 have to give at least one!: 605 606 >>> convert('kg','cgs',1.) 607 1000.0 608 >>> convert('g','SI',1.) 609 0.001 610 >>> convert('SI','g',1.) 611 1000.0 612 613 B{WARNINGS}: 614 1. the conversion involving sr and pixels is B{not tested}. 615 2. the conversion involving magnitudes is calibrated but not fully tested 616 3. non-integer powers are not functioning yet 617 618 Examples: 619 620 B{Spectra}: 621 622 >>> convert('AA','km/s',4553.,wave=(4552.,'AA')) 623 65.85950307557613 624 >>> convert('AA','km/s',4553.,wave=(4552.,0.1,'AA')) 625 (65.85950307557613, 6.587397133195861) 626 >>> convert('nm','m/s',455.3,wave=(0.4552,'mum')) 627 65859.50307564587 628 >>> convert('km/s','AA',65.859503075576129,wave=(4552.,'AA')) 629 4553.0 630 >>> convert('nm','Ghz',1000.) 631 299792.4579999999 632 >>> convert('km h-1','nRsol s-1',1.) 633 0.3993883287866966 634 >>> print(convert('erg s-1 cm-2 AA-1','SI',1.)) 635 10000000.0 636 637 B{Fluxes}: 638 639 >>> print(convert('erg/s/cm2/AA','Jy',1e-10,wave=(10000.,'angstrom'))) 640 333.564095198 641 >>> print(convert('erg/s/cm2/AA','Jy',1e-10,freq=(constants.cc/1e-6,'hz'))) 642 333.564095198 643 >>> print(convert('erg/s/cm2/AA','Jy',1e-10,freq=(constants.cc,'Mhz'))) 644 333.564095198 645 >>> print(convert('Jy','erg/s/cm2/AA',333.56409519815202,wave=(10000.,'AA'))) 646 1e-10 647 >>> print(convert('Jy','erg/s/cm2/AA',333.56409519815202,freq=(constants.cc,'Mhz'))) 648 1e-10 649 >>> convert('W/m2/mum','erg/s/cm2/AA',1e-10,wave=(10000.,'AA')) 650 1.0000000000000003e-11 651 >>> print(convert('Jy','W/m2/Hz',1.)) 652 1e-26 653 >>> print convert('W/m2/Hz','Jy',1.) 654 1e+26 655 >>> print convert('Jy','erg/cm2/s/Hz',1.) 656 1e-23 657 >>> print convert('erg/cm2/s/Hz','Jy',1.) 658 1e+23 659 >>> print(convert('Jy','erg/s/cm2',1.,wave=(2.,'micron'))) 660 1.49896229e-09 661 >>> print(convert('erg/s/cm2','Jy',1.,wave=(2.,'micron'))) 662 667128190.396 663 664 #>>> print(convert('Jy','erg/s/cm2/micron/sr',1.,wave=(2.,'micron'),ang_diam=(3.,'mas'))) 665 #4511059.82981 666 #>>> print(convert('Jy','erg/s/cm2/micron/sr',1.,wave=(2.,'micron'),pix=(3.,'mas'))) 667 #3542978.10531 668 #>>> print(convert('erg/s/cm2/micron/sr','Jy',1.,wave=(2.,'micron'),ang_diam=(3.,'mas'))) 669 #2.21677396826e-07 670 671 >>> print(convert('Jy','erg/s/cm2/micron',1.,wave=(2,'micron'))) 672 7.49481145e-10 673 >>> print(convert('10mW m-2 nm-1','erg s-1 cm-2 AA-1',1.)) 674 1.0 675 676 #>>> print convert('Jy','erg s-1 cm-2 micron-1 sr-1',1.,ang_diam=(2.,'mas'),wave=(1.,'micron')) 677 #40599538.4683 678 679 B{Angles}: 680 681 >>> convert('sr','deg2',1.) 682 3282.806350011744 683 684 >>> ang_diam = 3.21 # mas 685 >>> scale = convert('mas','sr',ang_diam/2.) 686 >>> print(ang_diam,scale) 687 (3.21, 6.054800067947964e-17) 688 >>> ang_diam = 2*convert('sr','mas',scale) 689 >>> print(ang_diam,scale) 690 (3.2100000000000004, 6.054800067947964e-17) 691 692 693 B{Magnitudes and amplitudes}: 694 695 >>> print(convert('ABmag','Jy',0.,photband='SDSS.U')) 696 3767.03798984 697 >>> print(convert('Jy','erg cm-2 s-1 AA-1',3630.7805477,wave=(1.,'micron'))) 698 1.08848062485e-09 699 >>> print(convert('ABmag','erg cm-2 s-1 AA-1',0.,wave=(1.,'micron'),photband='SDSS.G')) 700 4.97510278172e-09 701 >>> print(convert('erg cm-2 s-1 AA-1','ABmag',1e-8,wave=(1.,'micron'),photband='SDSS.G')) 702 -0.757994856607 703 >>> print(convert('ppm','muAmag',1.)) 704 1.0857356618 705 >>> print(convert('mAmag','ppt',1.,0.1)) 706 (0.9214583192957981, 0.09218827316735488) 707 >>> print(convert('mag_color','flux_ratio',0.599,0.004,photband='GENEVA.U-B')) 708 (1.1391327795013377, 0.004196720251233046) 709 710 B{Frequency analysis}: 711 712 >>> convert('cy/d','muHz',1.) 713 11.574074074074074 714 >>> convert('muhz','cy/d',11.574074074074074) 715 1.0 716 717 B{Interferometry}: 718 719 >>> convert('m/rad','cy/arcsec',85.,wave=(2.2,'micron')) 720 187.3143767923207 721 >>> convert('cm/rad','cy/arcmin',8500.,wave=(2200,'nm'))/60. 722 187.3143767923207 723 >>> convert('cy/arcsec','m/rad',187.,wave=(2.2,'mum')) 724 84.85734129005544 725 >>> convert('cyc/arcsec','m/rad',187.,wave=(1,'mum')) 726 38.571518768207014 727 >>> convert('cycles/arcsec','m',187.,freq=(300000.,'Ghz')) 728 38.54483473437971 729 >>> convert('cycles/mas','m',0.187,freq=(300000.,'Ghz')) 730 38.5448347343797 731 732 B{Temperature}: 733 734 >>> print(convert('Far','K',123.)) 735 323.705555556 736 >>> print(convert('kFar','kK',0.123)) 737 0.323705555556 738 >>> print(convert('K','Far',323.7)) 739 122.99 740 >>> print(convert('Cel','K',10.)) 741 283.15 742 >>> print(convert('Cel','Far',10.)) 743 50.0 744 >>> print(convert('dCel','kFar',100.)) 745 0.05 746 747 B{Time and Dates}: 748 749 >>> convert('sidereal d','d',1.) 750 1.0027379093 751 >>> convert('JD','CD',2446257.81458) 752 (1985.0, 7.0, 11.314580000005662) 753 >>> convert('CD','JD',(1985,7,11.31)) 754 2446257.81 755 >>> convert('CD','JD',(1985,7,11,7,31,59)) 756 2446257.813877315 757 >>> convert('MJD','CD',0.,jtype='corot') 758 (2000.0, 1.0, 1.5) 759 >>> convert('JD','MJD',2400000.5,jtype='mjd') 760 0.0 761 >>> convert('MJD','CD',0,jtype='mjd') 762 (1858.0, 11.0, 17.0) 763 764 B{Coordinates}: When converting coordinates with pyephem, you get pyephem 765 coordinates back. They have built-in string represenations and float 766 conversions: 767 768 >>> x,y = convert('equatorial','galactic',('17:45:40.4','-29:00:28.1'),epoch='2000') 769 >>> print (x,y) 770 (6.282224277178722, -0.000825178833899176) 771 >>> print x,y 772 359:56:41.8 -0:02:50.2 773 >>> x,y = convert('galactic','equatorial',('00:00:00.00','00:00:00.0'),epoch='2000') 774 >>> print (x,y) 775 (4.64964430303663, -0.5050315085342665) 776 >>> print x,y 777 17:45:37.20 -28:56:10.2 778 779 It is also possible to immediately convert to radians or degrees in floats 780 (this can be useful for plotting): 781 782 >>> convert('equatorial','deg_coord',('17:45:40.4','-29:00:28.1')) 783 (266.41833333333335, -29.00780555555556) 784 >>> convert('equatorial','rad_coord',('17:45:40.4','-29:00:28.1')) 785 (4.649877104342426, -0.5062817157227474) 786 787 B{Magnetism and Electricity}: The stored energy in a magnet, called magnet 788 performance or maximum energy product (often abbreviated BHmax), is 789 typically measured in units of megagauss-oersteds (MGOe). One MGOe is 790 approximately equal to 7957.74715 J/m3 (wikipedia): 791 792 >>> convert('MG Oe','J/m3',1.) 793 7957.747154594768 794 795 796 @param _from: units to convert from 797 @type _from: str 798 @param _to: units to convert to 799 @type _to: str 800 @keyword unpack: set to True if you don't want 'uncertainty objects'. If True 801 and uncertainties are given, they will be returned as a tuple (value, error) 802 instead of uncertainty object. Set to False probably only for internal uses 803 @type unpack: boolean, defaults to True 804 @return: converted value 805 @rtype: float 806 """ 807 #-- remember if user wants to unpack the results to have no trace of 808 # uncertainties, or wants to get uncertainty objects back 809 unpack = kwargs.pop('unpack',True) 810 811 #-- get the input arguments: if only one is given, it is either an 812 # C{uncertainty} from the C{uncertainties} package, or it is just a float 813 if len(args)==1: 814 start_value = args[0] 815 # if two arguments are given, we assume the first is the actual value and 816 # the second is the error on the value 817 elif len(args)==2: 818 start_value = unumpy.uarray([args[0],args[1]]) 819 else: 820 raise ValueError('illegal input') 821 822 #-- (un)logarithmicize (denoted by '[]') 823 m_in = re.search(r'\[(.*)\]',_from) 824 m_out = re.search(r'\[(.*)\]',_to) 825 826 if m_in is not None: 827 _from = m_in.group(1) 828 start_value = 10**start_value 829 if m_out is not None: 830 _to = m_out.group(1) 831 832 #-- It is possible the user gave a convention for either the from or to 833 # units (but not both!) 834 #-- break down the from and to units to their basic elements 835 if _from in _conventions: 836 _from = change_convention(_from,_to) 837 elif _to in _conventions: 838 _to = change_convention(_to,_from) 839 fac_from,uni_from = breakdown(_from) 840 fac_to,uni_to = breakdown(_to) 841 842 #-- convert the kwargs to SI units if they are tuples (make a distinction 843 # when uncertainties are given) 844 if uni_from!=uni_to and is_basic_unit(uni_from,'length') and not ('wave' in kwargs):# or 'freq' in kwargs_SI or 'photband' in kwargs_SI): 845 kwargs['wave'] = (start_value,_from) 846 logger.warning('Assumed input value to serve also for "wave" key') 847 elif uni_from!=uni_to and is_type(uni_from,'frequency') and not ('freq' in kwargs):# or 'freq' in kwargs_SI or 'photband' in kwargs_SI): 848 kwargs['freq'] = (start_value,_from) 849 logger.warning('Assumed input value to serve also for "freq" key') 850 kwargs_SI = {} 851 for key in kwargs: 852 if isinstance(kwargs[key],tuple): 853 kwargs_SI[key] = convert(kwargs[key][-1],'SI',*kwargs[key][:-1],unpack=False) 854 else: 855 kwargs_SI[key] = kwargs[key] 856 #-- add some default values if necessary 857 logger.debug('Convert %s to %s, fac_from-start_value %s/%s'%(uni_from,uni_to,fac_from,start_value)) 858 859 #-- conversion is easy if same units 860 ret_value = 1. 861 862 if uni_from==uni_to: 863 #-- if nonlinear conversions from or to: 864 if isinstance(fac_from,NonLinearConverter): 865 ret_value *= fac_from(start_value,**kwargs_SI) 866 else: 867 try: 868 ret_value *= fac_from*start_value 869 except TypeError: 870 raise TypeError('Cannot multiply value with a float; probably argument is a tuple (value,error), please expand with *(value,error)') 871 872 #-- otherwise a little bit more complicated 873 else: 874 #-- first check where the unit differences are 875 uni_from_ = uni_from.split() 876 uni_to_ = uni_to.split() 877 only_from_c,only_to_c = sorted(list(set(uni_from_) - set(uni_to_))),sorted(list(set(uni_to_) - set(uni_from_))) 878 only_from_c,only_to_c = [list(components(i))[1:] for i in only_from_c],[list(components(i))[1:] for i in only_to_c] 879 #-- push them all bach to the left side (change sign of right hand side components) 880 left_over = " ".join(['%s%d'%(i,j) for i,j in only_from_c]) 881 left_over+= " "+" ".join(['%s%d'%(i,-j) for i,j in only_to_c]) 882 left_over = breakdown(left_over)[1] 883 #-- but be sure to convert everything to SI units so that the switch 884 # can be interpreted. 885 left_over = [change_convention('SI',ilo) for ilo in left_over.split()] 886 only_from = "".join(left_over) 887 only_to = '' 888 889 #-- first we remove any differences concerning (ster)radians 890 # we recently added fac_from* to all these things, maybe this needs to 891 # change? 892 #if 'rad2' in only_from: 893 # start_value = fac_from*_switch['rad2_to_'](start_value,**kwargs_SI) 894 # only_from = only_from.replace('rad2','') 895 # logger.debug('Switching to /sr') 896 # fac_from = 1. 897 #elif 'rad-2' in only_from: 898 # start_value = fac_from*_switch['rad-2_to_'](start_value,**kwargs_SI) 899 # only_from = only_from.replace('rad-2','') 900 # logger.debug('Switching from /sr') 901 # fac_from = 1. 902 #elif 'rad1' in only_from: 903 # start_value = fac_from*_switch['rad1_to_'](start_value,**kwargs_SI) 904 # only_from = only_from.replace('rad1','') 905 # logger.debug('Switching to /rad') 906 # fac_from = 1. 907 #elif 'rad-1' in only_from: 908 # start_value = fac_from*_switch['rad-1_to_'](start_value,**kwargs_SI) 909 # only_from = only_from.replace('rad-1','') 910 # logger.debug('Switching from /rad') 911 # fac_from = 1. 912 913 #-- then we do what is left over (if anything is left over) 914 if only_from or only_to: 915 logger.debug("Convert %s to %s"%(only_from,only_to)) 916 917 #-- nonlinear conversions need a little tweak 918 try: 919 key = '%s_to_%s'%(only_from,only_to) 920 logger.debug('Switching from {} to {} via {:s}'.format(only_from,only_to,_switch[key].__name__)) 921 if isinstance(fac_from,NonLinearConverter): 922 ret_value *= _switch[key](fac_from(start_value,**kwargs_SI),**kwargs_SI) 923 #-- linear conversions are easy 924 else: 925 logger.debug('fac_from=%s, start_value=%s with kwargs %s'%(fac_from,start_value,kwargs_SI)) 926 ret_value *= _switch[key](fac_from*start_value,**kwargs_SI) 927 except KeyError: 928 #-- try to be smart an reverse the units: 929 if not (Unit(1.,uni_from)*Unit(1.,uni_to))[1]: 930 ret_value *= period2freq(fac_from*start_value,**kwargs_SI) 931 logger.warning('It is assumed that the "from" unit is the inverse of the "to" unit') 932 else: 933 logger.critical('cannot convert %s to %s: no %s definition in dict _switch'%(_from,_to,key)) 934 raise 935 else: 936 ret_value *= start_value 937 #-- final step: convert to ... (again distinction between linear and 938 # nonlinear converters) 939 if isinstance(fac_to,NonLinearConverter): 940 ret_value = fac_to(ret_value,inv=True,**kwargs_SI) 941 else: 942 ret_value /= fac_to 943 944 #-- logarithmicize 945 if m_out is not None: 946 ret_value = log10(ret_value) 947 948 #-- unpack the uncertainties if: 949 # 1. the input was not given as an uncertainty 950 # 2. the input was without uncertainties, but extra keywords had uncertainties 951 # 3. the input was with uncertainties (float or array) and unpack==True 952 unpack_case1 = len(args)==2 953 unpack_case2 = len(args)==1 and isinstance(ret_value,AffineScalarFunc) 954 unpack_case3 = len(args)==1 and isinstance(ret_value,np.ndarray) and isinstance(ret_value[0],AffineScalarFunc) 955 if unpack and (unpack_case1 or unpack_case2): 956 ret_value = unumpy.nominal_values(ret_value),unumpy.std_devs(ret_value) 957 #-- convert to real floats if real floats were given 958 if not ret_value[0].shape: 959 ret_value = np.asscalar(ret_value[0]),np.asscalar(ret_value[1]) 960 961 return ret_value

962

963 964 -def nconvert(_froms,_tos,*args,**kwargs):

965 """ 966 Convert a list/array/tuple of values with different units to other units. 967 968 This silently catches some exceptions and replaces the value with nan! 969 970 @rtype: array 971 """ 972 if len(args)==1: 973 ret_value = np.zeros((len(args[0]))) 974 elif len(args)==2: 975 ret_value = np.zeros((len(args[0]),2)) 976 if isinstance(_tos,str): 977 _tos = [_tos for i in _froms] 978 elif isinstance(_froms,str): 979 _froms = [_froms for i in _tos] 980 981 for i,(_from,_to) in enumerate(list(zip(_froms,_tos))): 982 myargs = [iarg[i] for iarg in args] 983 mykwargs = {} 984 for key in kwargs: 985 if not isinstance(kwargs[key],str) and hasattr(kwargs[key],'__iter__'): 986 mykwargs[key] = kwargs[key][i] 987 else: 988 mykwargs[key] = kwargs[key] 989 try: 990 ret_value[i] = convert(_from,_to,*myargs,**mykwargs) 991 except ValueError: #no calibration 992 ret_value[i] = np.nan 993 994 if len(args)==2: 995 ret_value = ret_value.T 996 return ret_value

997

998 999 -def change_convention(to_,units,origin=None):

1000 """ 1001 Change units from one convention to another. 1002 1003 Example usage: 1004 1005 >>> units1 = 'kg m2 s-2 K-1 mol-1' 1006 >>> print change_convention('cgs',units1) 1007 K-1 g1 cm2 mol-1 s-2 1008 >>> print change_convention('sol','g cm2 s-2 K-1 mol-1') 1009 Tsol-1 Msol1 Rsol2 mol-1 s-2 1010 1011 @param to_: convention name to change to 1012 @type to_: str 1013 @param units: units to change 1014 @type units: str 1015 @param origin: the name of the original convention 1016 @type origin: str 1017 @return: units in new convention 1018 @rtype: str 1019 """ 1020 if origin is None: 1021 origin = constants._current_convention 1022 #-- break down units in base units, and breakdown that string in 1023 # whole units and non-alpha digits (we'll weave them back in after 1024 # translation) 1025 factor,units = breakdown(units) 1026 new_units = re.findall(r'[a-z]+',units, re.I) 1027 powers = re.findall(r'[0-9\W]+',units, re.I) 1028 #-- make the translation dictionary 1029 translator = {} 1030 for key in sorted(_conventions[origin].keys()): 1031 translator[_conventions[origin][key]] = _conventions[to_][key] 1032 #-- translate 1033 new_units = [unit in translator and translator[unit] or unit for unit in new_units] 1034 #-- weave them back in 1035 new_units = "".join(["".join([i,j]) for i,j in zip(new_units,powers)]) 1036 return new_units

1037

1038 -def set_convention(units='SI',values='standard',frequency='rad'):

1039 """ 1040 Consistently change the values of the fundamental constants to the paradigm 1041 of other system units or programs. 1042 1043 This can be important in pulsation codes, where e.g. the value of the 1044 gravitational constant and the value of the solar radius can have siginficant 1045 influences on the frequencies. 1046 1047 Setting C{frequency} to C{rad} gives you the following behaviour: 1048 1049 >>> old_settings = set_convention(frequency='rad') 1050 >>> convert('s-1','rad/s',1.0) 1051 1.0 1052 >>> convert('cy/s','rad/s',1.0) 1053 6.283185307179586 1054 >>> convert('Hz','rad/s',1.0) 1055 6.283185307179586 1056 >>> convert('','mas',constants.au/(1000*constants.pc)) 1057 0.99999999999623 1058 1059 Setting C{frequency} to C{cy} gives you the following behaviour: 1060 1061 >>> old_settings = set_convention(frequency='cy') 1062 >>> convert('s-1','rad/s',1.0) 1063 6.283185307179586 1064 >>> convert('cy/s','rad/s',1.0) 1065 6.283185307179586 1066 >>> convert('Hz','rad/s',1.0) 1067 6.283185307179586 1068 >>> convert('','mas',constants.au/(1000*constants.pc)) 1069 6.2831853071558985 1070 >>> old_settings = set_convention(*old_settings) 1071 1072 @param units: name of the new units base convention 1073 @type units: str 1074 @param values: name of the value set for fundamental constants 1075 @type values: str 1076 @return: name of old convention 1077 @rtype: str 1078 """ 1079 to_return = constants._current_convention,\ 1080 constants._current_values,\ 1081 constants._current_frequency 1082 values = values.lower() 1083 if to_return==(units,values,frequency): 1084 logger.info('No need to change convention or values') 1085 return to_return 1086 if to_return[2]!=frequency and 'cy' in frequency.lower(): 1087 _switch['rad1_to_'] = per_cy 1088 _switch['rad-1_to_'] = times_cy 1089 constants._current_frequency = frequency.lower() 1090 logger.debug('Changed frequency convention to {0}'.format(frequency)) 1091 elif to_return[2]!=frequency and 'rad' in frequency.lower(): 1092 _switch['rad1_to_'] = do_nothing 1093 _switch['rad-1_to_'] = do_nothing 1094 constants._current_frequency = frequency.lower() 1095 logger.debug('Changed frequency convention to {0}'.format(frequency)) 1096 1097 if to_return[:2]==(units,values): 1098 return to_return 1099 #-- first reload the constants to their original values (i.e. SI based) 1100 #reload(constants) 1101 #-- then, where possible, replace all the constants with the value from 1102 # the other convention: 1103 cvars = dir(constants) 1104 cvars = [i for i in cvars if i+'_units' in cvars] 1105 for const in cvars: 1106 #-- break down the old units in their basic components, but replace 1107 # the value of the constant from those of the C{values} set if 1108 # possible 1109 old_units = getattr(constants,const+'_units') 1110 const_ = const+'_'+values 1111 if hasattr(constants,const_): 1112 #logger.info('Override {0} value to {1}'.format(const,values)) 1113 old_value = getattr(constants,const_) 1114 #-- remember, the old values is always in SI 1115 if not to_return[1]=='SI': 1116 SI_units = change_convention('SI',old_units) 1117 old_value = convert(SI_units,old_units,old_value) 1118 else: 1119 old_value = getattr(constants,const) 1120 factor,old_units = breakdown(old_units) 1121 #-- replace old base units with new conventions's base units 1122 new_units = change_convention(units,old_units) 1123 #-- convert the value from the old to the new convenction 1124 new_value = convert(old_units,new_units,old_value) 1125 #-- and attach to the constants module 1126 setattr(constants,const,new_value) 1127 setattr(constants,const+'_units',new_units) 1128 #-- convert the _factors in this module to the new convention. First make 1129 # new list of factors, only to later overwrite the original one as a whole. 1130 # We do this because e.g. while converting everything to yd, we would 1131 # overwrite the definition of 'm' with that of 'yd', so that we would lose 1132 # the conversion information on the fly 1133 new_factors = {} 1134 for fac in _factors: 1135 _from = _factors[fac][1] 1136 #-- some quantities have no units. 1137 if not _from: 1138 new_factors[fac] = _factors[fac] 1139 continue 1140 _to = change_convention(units,_from) 1141 #-- if this unit is one of the base units of any system, make sure not 1142 # to set the '1' at the end, or we'll get in trouble... 1143 try: 1144 for conv in _conventions: 1145 for typ in _conventions[conv]: 1146 if _conventions[conv][typ]==_to[:-1] and _to[-1]=='1': 1147 _to = _to[:-1] 1148 raise StopIteration 1149 except StopIteration: 1150 pass 1151 if _from==_to[:-1] and _to[-1]=='1': 1152 _to = _to[:-1] 1153 #-- some stuff cannot be converted 1154 try: 1155 new_factors[fac] = (convert(_from,_to,_factors[fac][0]),_to,_factors[fac][2],_factors[fac][3]) 1156 except ValueError: 1157 new_factors[fac] = _factors[fac] 1158 continue 1159 except TypeError: 1160 new_factors[fac] = _factors[fac] 1161 continue 1162 for fac in _factors: 1163 _factors[fac] = new_factors[fac] 1164 #-- convert the names of combinations of basic units: 1165 for name in list(_names.keys()): 1166 new_name = change_convention(units,name) 1167 _names[new_name] = _names.pop(name) 1168 #-- convert the switches in this module to the new convention 1169 #for switch in _switch: 1170 1171 constants._current_convention = units 1172 constants._current_values = values 1173 #-- when we set everything back to SI, make sure we have no rounding errors: 1174 if units=='SI' and values=='standard' and frequency=='rad': 1175 reload(constants) 1176 logger.warning('Reloading of constants') 1177 logger.info('Changed convention to {0} with values from {1} set'.format(units,values)) 1178 return to_return

1179

1180 -def reset():

1181 """ 1182 Resets all values and conventions to the original values. 1183 """ 1184 set_convention()

1185

1186 -def get_convention():

1187 """ 1188 Returns convention and name of set of values. 1189 1190 >>> get_convention() 1191 ('SI', 'standard', 'rad') 1192 1193 @rtype: str,str 1194 @return: convention, values 1195 """ 1196 return constants._current_convention,constants._current_values,constants._current_frequency

1197

1198 1199 -def get_constant(constant_name,units='SI',value='standard'):

1200 """ 1201 Convenience function to retrieve the value of a constant in a particular 1202 system. 1203 1204 >>> Msol = get_constant('Msol','SI') 1205 >>> Msol = get_constant('Msol','kg') 1206 1207 @param constant_name: name of the constant 1208 @type constant_name: str 1209 @param units: name of the unit base system 1210 @type units: str 1211 @param value: name of the parameter set the get the value from 1212 @type value: str 1213 @return: value of the constant in the unit base system 1214 @rtype: float 1215 """ 1216 value = value.lower() 1217 #-- see if we need (and have) a different value than the standard one 1218 const_ = constant_name+'_'+value 1219 old_units = getattr(constants,constant_name+'_units') 1220 if hasattr(constants,const_): 1221 old_value = getattr(constants,const_) 1222 old_units = change_convention('SI',old_units) 1223 else: 1224 old_value = getattr(constants,constant_name) 1225 new_value = convert(old_units,units,old_value) 1226 return new_value

1227

1228 -def get_constants(units='SI',values='standard'):

1229 """ 1230 Convenience function to retrieve the value of all constants in a particular 1231 system. 1232 1233 This can be helpful when you want to attach a copy of the fundamental constants 1234 to some module or class instances, and be B{sure} these values never change. 1235 1236 Yes, I know, there's a lot that can go wrong with values that are supposed 1237 to be CONSTANT! 1238 1239 @rtype: dict 1240 """ 1241 cvars = dir(constants) 1242 cvars = [i for i in cvars if i+'_units' in cvars] 1243 myconstants = {} 1244 for cvar in cvars: 1245 myconstants[cvar] = get_constant(cvar,units,values) 1246 return myconstants

1247

1248 1249 1250 1251 1252 #} 1253 #{ Conversions basics and helper functions 1254 1255 1256 -def solve_aliases(unit):

1257 """ 1258 Resolve simple aliases in a unit's name. 1259 1260 Resolves aliases and replaces division signs with negative power. 1261 1262 @param unit: unit (e.g. erg s-1 angstrom-1) 1263 @type unit: str 1264 @return: aliases-resolved unit (e.g. erg s-1 AA-1) 1265 @rtype: str 1266 """ 1267 #-- resolve aliases 1268 for alias in _aliases: 1269 unit = unit.replace(alias[0],alias[1]) 1270 1271 #-- replace slash-forward with negative powers 1272 if '/' in unit: 1273 unit_ = [uni.split('/') for uni in unit.split()] 1274 for i,uni in enumerate(unit_): 1275 for j,after_div in enumerate(uni[1:]): 1276 if not after_div[-1].isdigit(): after_div += '1' 1277 #m = re.search(r'(\d*)(.+?)(-{0,1}\d+)',after_div) 1278 m = re.search(r'(\d*)(.+?)(-{0,1}[\d\.]+)',after_div) 1279 if m is not None: 1280 factor,basis,power = m.group(1),m.group(2),m.group(3) 1281 if not '.' in power: power = int(power) 1282 else: power = float(power) 1283 if factor: factor = float(factor) 1284 else: factor = 1. 1285 else: 1286 factor,basis,power = 1.,after_div,1 1287 uni[1+j] = '%s%g'%(basis,-power) 1288 if factor!=1: uni[1+j] = '%d%s'%(factor,uni[1+j]) 1289 ravelled = [] 1290 for uni in unit_: 1291 ravelled += uni 1292 unit = " ".join(ravelled) 1293 return unit

1294

1295 1296 -def components(unit):

1297 """ 1298 Decompose a unit into: factor, SI base units, power. 1299 1300 You probably want to call L{solve_aliases} first. 1301 1302 Examples: 1303 1304 >>> factor1,units1,power1 = components('m') 1305 >>> factor2,units2,power2 = components('g2') 1306 >>> factor3,units3,power3 = components('hg3') 1307 >>> factor4,units4,power4 = components('Mg4') 1308 >>> factor5,units5,power5 = components('mm') 1309 >>> factor6,units6,power6 = components('W3') 1310 >>> factor7,units7,power7 = components('s-2') 1311 1312 >>> print "{0}, {1}, {2}".format(factor1,units1,power1) 1313 1.0, m, 1 1314 >>> print "{0}, {1}, {2}".format(factor2,units2,power2) 1315 0.001, kg, 2 1316 >>> print "{0}, {1}, {2}".format(factor3,units3,power3) 1317 0.1, kg, 3 1318 >>> print "{0}, {1}, {2}".format(factor4,units4,power4) 1319 1000.0, kg, 4 1320 >>> print "{0}, {1}, {2}".format(factor5,units5,power5) 1321 0.001, m, 1 1322 >>> print "{0}, {1}, {2}".format(factor6,units6,power6) 1323 1.0, m2 kg1 s-3, 3 1324 >>> print "{0}, {1}, {2}".format(factor7,units7,power7) 1325 1.0, s, -2 1326 1327 @param unit: unit name 1328 @type unit: str 1329 @return: 3-tuple with factor, SI base unit and power 1330 @rtype: (float,str,int) 1331 """ 1332 factor = 1. 1333 #-- manually check if there is a prefactor of the form '10-14' or '10e-14' 1334 # and expand it if necessary 1335 m = re.search('\d\d[-+]\d\d',unit[:5]) 1336 if m is not None: 1337 factor *= float(m.group(0)[:2])**(float(m.group(0)[2]+'1')*float(m.group(0)[3:5])) 1338 unit = unit[5:] 1339 m = re.search('\d\d[[eE][-+]\d\d',unit[:6]) 1340 if m is not None: 1341 factor *= float(m.group(0)[:2])**(float(m.group(0)[3]+'1')*float(m.group(0)[4:6])) 1342 unit = unit[6:] 1343 1344 if not unit[-1].isdigit(): unit += '1' 1345 #-- decompose unit in base name and power 1346 #m = re.search(r'(\d*)(.+?)(-{0,1}\d+)',unit) 1347 m = re.search(r'(\d*)(.+?)(-{0,1}[\d\.]+)',unit) 1348 if m is not None: 1349 factor_,basis,power = m.group(1),m.group(2),m.group(3) 1350 #-- try to make power an integer, otherwise make it a float 1351 if not '.' in power: power = int(power) 1352 else: power = float(power) 1353 if factor_: 1354 factor *= float(factor_) 1355 else: 1356 factor,basis,power = 1.,unit,1 1357 #-- decompose the base name (which can be a composition of a prefix 1358 # (e.g., 'mu') and a unit name (e.g. 'm')) into prefix and unit name 1359 #-- check if basis is part of _factors dictionary. If not, find the 1360 # combination of _scalings and basis which is inside the dictionary! 1361 1362 for scale in _scalings: 1363 scale_unit,base_unit = basis[:len(scale)],basis[len(scale):] 1364 if scale_unit==scale and base_unit in _factors and not basis in _factors: 1365 #if basis in _factors: 1366 # raise ValueError,'ambiguity between %s and %s-%s'%(basis,scale_unit,base_unit) 1367 factor *= _scalings[scale] 1368 basis = base_unit 1369 break 1370 #-- if we didn't find any scalings, check if the 'raw' unit is already 1371 # a base unit 1372 else: 1373 if not basis in _factors: 1374 raise ValueError('Unknown unit %s'%(basis)) 1375 1376 #-- switch from base units to SI units 1377 if hasattr(_factors[basis][0],'__call__'): 1378 factor = factor*_factors[basis][0]() 1379 else: 1380 factor *= _factors[basis][0] 1381 basis = _factors[basis][1] 1382 1383 return factor,basis,power

1384

1385 1386 1387 -def breakdown(unit):

1388 """ 1389 Decompose a unit into SI base units containing powers. 1390 1391 Examples: 1392 1393 >>> factor1,units1 = breakdown('erg s-1 W2 kg2 cm-2') 1394 >>> factor2,units2 = breakdown('erg s-1 cm-2 AA-1') 1395 >>> factor3,units3 = breakdown('W m-3') 1396 1397 1398 >>> print "{0}, {1}".format(factor1,units1) 1399 0.001, kg5 m4 s-9 1400 >>> print "{0}, {1}".format(factor2,units2) 1401 10000000.0, kg1 m-1 s-3 1402 >>> print "{0}, {1}".format(factor3,units3) 1403 1.0, kg1 m-1 s-3 1404 1405 @param unit: unit's name 1406 @type unit: str 1407 @return: 2-tuple factor, unit's base name 1408 @rtype: (float,str) 1409 """ 1410 #-- solve aliases 1411 unit = solve_aliases(unit) 1412 #-- break down in basic units 1413 units = unit.split() 1414 total_factor = 1. 1415 total_units = [] 1416 total_power = [] 1417 for unit in units: 1418 factor,basis,power = components(unit) 1419 total_factor = total_factor*factor**power 1420 basis = basis.split() 1421 for base in basis: 1422 factor_,basis_,power_ = components(base) 1423 try: 1424 total_factor = total_factor*factor_ 1425 except TypeError: 1426 pass 1427 if basis_ in total_units: 1428 index = total_units.index(basis_) 1429 total_power[index] += power_*power 1430 else: 1431 total_units.append(basis_) 1432 total_power.append(power_*power) 1433 1434 #-- make sure to return a sorted version 1435 total_units = sorted(['%s%s'%(i,j) for i,j in zip(total_units,total_power) if j!=0]) 1436 return total_factor," ".join(total_units)

1437

1438 -def compress(unit,ignore_factor=False):

1439 """ 1440 Compress basic units to more human-readable type. 1441 1442 If you are only interested in a shorter form of the units, regardless 1443 of the values, you need to set C{ignore_factor=True}. 1444 1445 For example: erg/s/cm2/AA can be written shorter as W1 m-3, but 1 erg/s/cm2/AA 1446 is not equal to 1 W1 m-3. Thus: 1447 1448 >>> print(compress('erg/s/cm2/AA',ignore_factor=True)) 1449 W1 m-3 1450 >>> print(compress('erg/s/cm2/AA',ignore_factor=False)) 1451 erg/s/cm2/AA 1452 1453 You cannot write the latter shorter without converting the units! 1454 1455 For the same reasoning, compressing non-basic units will not work out: 1456 1457 >>> print(compress('Lsol')) 1458 Lsol 1459 >>> print(compress('Lsol',ignore_factor=True)) 1460 W1 1461 1462 So actually, this function is really designed to compress a combination 1463 of basic units: 1464 1465 >>> print(compress('kg1 m-1 s-3')) 1466 W1 m-3 1467 1468 @param unit: combination of units to compress 1469 @type unit: str 1470 @param ignore_factor: if True, then a `sloppy compress' will be performed, 1471 that is, it is not guaranteed that the multiplication factor equals 1 (see 1472 examples) 1473 @type ignore_factor: bool 1474 @rtype: str 1475 @return: shorter version of original units 1476 1477 """ 1478 #-- if input is not a basic unit (or is not trivially expandable as 1479 # a basic unit, do nothing) 1480 breakdown_units = breakdown(unit) 1481 if not ignore_factor and breakdown_units[0]!=1.: 1482 return unit 1483 #-- else, look for ways to simplify it 1484 orig_units = Unit(1.,breakdown_units[1]) 1485 for fac in _factors: 1486 try: 1487 if not np.allclose(_factors[fac][0],1.): continue 1488 y = Unit(1.,fac+'1') 1489 left_over = (orig_units/y) 1490 value,left_over = left_over[0],left_over[1] 1491 except: 1492 continue 1493 new_unit = ' '.join([y[1],left_over]) 1494 if len(new_unit)<len(orig_units[1]): 1495 orig_units = Unit(1.,new_unit) 1496 1497 return orig_units[1]

1498

1499 -def split_units_powers(unit):

1500 """ 1501 Split unit in a list of units and a list of powers. 1502 1503 >>> y = 'AA A kg-2 F-11 Wb-1' 1504 >>> split_units_powers(y) 1505 (['AA', 'A', 'kg', 'F', 'Wb'], ['1', '1', '-2', '-11', '-1']) 1506 >>> y = 'angstrom A1/kg2 F-11/Wb' 1507 >>> split_units_powers(y) 1508 (['AA', 'A', 'kg', 'F', 'Wb'], ['1', '1', '-2', '-11', '-1']) 1509 >>> y = '10angstrom A1/kg2 F-11/Wb' 1510 >>> split_units_powers(y) 1511 (['10AA', 'A', 'kg', 'F', 'Wb'], ['1', '1', '-2', '-11', '-1']) 1512 >>> y = '10angstrom A1/kg2 5F-11/Wb' 1513 >>> split_units_powers(y) 1514 (['10AA', 'A', 'kg', '5F', 'Wb'], ['1', '1', '-2', '-11', '-1']) 1515 >>> y = '10angstrom A0.5/kg2 5F-11/100Wb2.5' 1516 >>> split_units_powers(y) 1517 (['10AA', 'A', 'kg', '5F', '100Wb'], ['1', '0.5', '-2', '-11', '-2.5']) 1518 1519 @param unit: unit 1520 @type unit: str 1521 @rtype: list,list 1522 @return: list of unit names, list of powers 1523 """ 1524 #-- compile regular expressions 1525 units = re.compile('\A[\d].[a-z]+|[a-z]+|\s\d+[a-z]+',re.IGNORECASE) 1526 powers = re.compile('[-\.\d]+[\s]|[-\.\d]+\Z') 1527 #-- solve aliases 1528 unit = solve_aliases(unit) 1529 #-- and make sure every unit has at least one power: 1530 unit = ' '.join([comp[-1].isdigit() and comp or comp+'1' for comp in unit.split()]) 1531 #-- now split and remove spaces 1532 units = [i.strip() for i in units.findall(unit)] 1533 powers = [i.strip() for i in powers.findall(unit)] 1534 return units,powers

1535

1536 1537 -def is_basic_unit(unit,type):

1538 """ 1539 Check if a unit is represents one of the basic quantities (length, mass...) 1540 1541 >>> is_basic_unit('K','temperature') 1542 True 1543 >>> is_basic_unit('m','length') 1544 True 1545 >>> is_basic_unit('cm','length') 1546 True 1547 >>> is_basic_unit('km','length') # not a basic unit in any type of convention 1548 False 1549 1550 @parameter unit: unit to check 1551 @type unit: str 1552 @parameter type: type to check 1553 @type type: str 1554 @rtype: bool 1555 """ 1556 for conv in _conventions: 1557 if unit==_conventions[conv][type] or (unit[:-1]==_conventions[conv][type] and unit[-1]=='1'): 1558 return True 1559 else: 1560 return False

1561

1562 -def is_type(unit,type):

1563 """ 1564 Check if a unit is of a certain type (mass, length, force...) 1565 1566 >>> is_type('K','temperature') 1567 True 1568 >>> is_type('m','length') 1569 True 1570 >>> is_type('cm','length') 1571 True 1572 >>> is_type('ly','length') 1573 True 1574 >>> is_type('ly','time') 1575 False 1576 >>> is_type('W','power') 1577 True 1578 >>> is_type('kg m2/s3','power') 1579 True 1580 >>> is_type('kg/s2/m','pressure') 1581 True 1582 >>> is_type('W','force') 1583 False 1584 1585 @parameter unit: unit to check 1586 @type unit: str 1587 @parameter type: type to check 1588 @type type: str 1589 @rtype: bool 1590 """ 1591 #-- solve aliases and breakdown in basic SI units 1592 unit = solve_aliases(unit) 1593 comps = breakdown(unit)[1] 1594 for fac in _factors: 1595 this_unit = change_convention(constants._current_convention,_factors[fac][1]) 1596 if comps==this_unit and type in _factors[fac][2].split('/'): 1597 return True 1598 #-- via names 1599 if comps in _names and _names[comps]==type: 1600 return True 1601 1602 return False

1603

1604 -def get_type(unit):

1605 """ 1606 Return human readable name of the unit 1607 1608 @rtype: str 1609 """ 1610 comps = breakdown(unit)[1] 1611 if comps in _names: 1612 return _names[comps] 1613 for i in _factors: 1614 test_against = breakdown(i)[1] 1615 if test_against==comps or unit==test_against: 1616 return _factors[i][2] 1617 return unit

1618

1619 -def round_arbitrary(x, base=5):

1620 """ 1621 Round to an arbitrary base. 1622 1623 Example usage: 1624 1625 >>> round_arbitrary(1.24,0.25) 1626 1.25 1627 >>> round_arbitrary(1.37,0.75) 1628 1.5 1629 1630 1631 @param x: number to round 1632 @type x: float 1633 @param base: base to round to 1634 @type base: integer or float 1635 @return: rounded number 1636 @rtype: float 1637 """ 1638 return base * round(float(x)/base)

1639

1640 -def unit2texlabel(unit,full=False):

1641 """ 1642 Convert a unit string to a nicely formatted TeX label. 1643 1644 For fluxes, also the Flambda, lambda Flambda, or Fnu, or nuFnu is returned. 1645 1646 @param unit: unit name 1647 @type unit: str 1648 @param full: return "name [unit]" or only "unit" 1649 @type full: bool 1650 @rtype: str 1651 """ 1652 unit = solve_aliases(unit) 1653 fac,base = breakdown(unit) 1654 names,powers = split_units_powers(unit) 1655 powers = [(power!='1' and power or ' ') for power in powers] 1656 powers = [(power[0]!=' ' and '$^{{{0}}}$'.format(power) or power) for power in powers] 1657 unit_ = ' '.join([(name+power) for name,power in zip(names,powers)]) 1658 1659 translate = {'kg1 m-1 s-3' :r'$F_\lambda$ [{0}]'.format(unit_), 1660 'cy-1 kg1 s-2':r'$F_\nu$ [{0}]'.format(unit_), 1661 'kg1 s-3':r'$\lambda F_\lambda$ [{0}]'.format(unit_), 1662 } 1663 #translate = {} 1664 #-- translate 1665 if base in translate: 1666 label = translate[base] 1667 else: 1668 label = unit_ 1669 #-- make TeX 1670 label = label.replace('AA','$\AA$') 1671 label = label.replace('mu','$\mu$') 1672 label = label.replace('sol','$_\odot$') 1673 1674 if full: 1675 label = '{0} [{1}]'.format(get_type(unit).title(),label.strip()) 1676 1677 return label

1678

1679 1680 1681 1682 1683 -def get_help():

1684 """ 1685 Return a string with a list and explanation of all defined units. 1686 1687 @rtype: str 1688 """ 1689 try: 1690 set_exchange_rates() 1691 except IOError: 1692 logger.warning('Unable to connect to ecb.europa.eu') 1693 help_text = {} 1694 for fac in sorted(_factors.keys()): 1695 if _factors[fac][2] not in help_text: 1696 help_text[_factors[fac][2]] = [] 1697 help_text[_factors[fac][2]].append('%-15s = %-30s'%(fac,_factors[fac][3])) 1698 text = [[],[]] 1699 bar = '%-48s'%('================='+20*'=') 1700 for i,key in enumerate(sorted(help_text.keys())): 1701 text[i%2] += [bar,'%-48s'%("= Units of %-20s ="%(key)),bar] 1702 text[i%2] += help_text[key] 1703 out = '' 1704 #for i,j in itertools.zip_longest(*text,fillvalue=''): # for Python 3 1705 for i,j in itertools.izip_longest(*text,fillvalue=''): 1706 out += '%s| %s\n'%(i,j) 1707 1708 return out

1709

1710 -def units2sphinx():

1711 set_exchange_rates() 1712 text = [] 1713 divider = "+{0:s}+{1:s}+{2:s}+{3:s}+{4:s}+".format(22*'-',52*'-',22*'-',32*'-',52*'-') 1714 text.append(divider) 1715 text.append("| {0:20s} | {1:50} | {2:20s} | {3:30s} | {4:50s} |".format('name','factor','units','type','description')) 1716 text.append(divider) 1717 for key in _factors: 1718 text.append("| {0:20s} | {1:50} | {2:20s} | {3:30s} | {4:50s} |".format(key,*_factors[key])) 1719 text.append(divider) 1720 return "\n".join(text)

1721

1722 -def derive_wrapper(fctn):

1723 @functools.wraps(fctn) 1724 def func(*args,**kwargs): 1725 argspec = inspect.getargspec(fctn) 1726 print argspec 1727 result = fctn(*args,**kwargs) 1728 return result

1729 return func 1730 #}

1731 #{ Linear change-of-base conversions 1732 1733 -def distance2velocity(arg,**kwargs):

1734 """ 1735 Switch from distance to velocity via a reference wavelength. 1736 1737 @param arg: distance (SI, m) 1738 @type arg: float 1739 @keyword wave: reference wavelength (SI, m) 1740 @type wave: float 1741 @return: velocity (SI, m/s) 1742 @rtype: float 1743 """ 1744 if 'wave' in kwargs: 1745 wave = kwargs['wave'] 1746 velocity = (arg-wave) / wave * constants.cc 1747 else: 1748 raise ValueError('reference wavelength (wave) not given') 1749 return velocity

1750

1751 -def velocity2distance(arg,**kwargs):

1752 """ 1753 Switch from velocity to distance via a reference wavelength. 1754 1755 @param arg: velocity (SI, m/s) 1756 @type arg: float 1757 @keyword wave: reference wavelength (SI, m) 1758 @type wave: float 1759 @return: distance (SI, m) 1760 @rtype: float 1761 """ 1762 if 'wave' in kwargs: 1763 wave = kwargs['wave'] 1764 distance = wave / constants.cc * arg + wave 1765 else: 1766 raise ValueError('reference wavelength (wave) not given') 1767 return distance

1768

1769 -def fnu2flambda(arg,**kwargs):

1770 """ 1771 Switch from Fnu to Flambda via a reference wavelength. 1772 1773 Flambda and Fnu are spectral irradiance in wavelength and frequency, 1774 respectively 1775 1776 @param arg: spectral irradiance (SI,W/m2/Hz) 1777 @type arg: float 1778 @keyword photband: photometric passband 1779 @type photband: str ('SYSTEM.FILTER') 1780 @keyword wave: reference wavelength (SI, m) 1781 @type wave: float 1782 @keyword freq: reference frequency (SI, Hz) 1783 @type freq: float 1784 @return: spectral irradiance (SI, W/m2/m) 1785 @rtype: float 1786 """ 1787 if 'photband' in kwargs: 1788 lameff = filters.eff_wave(kwargs['photband']) 1789 lameff = convert('AA','m',lameff) 1790 kwargs['wave'] = lameff 1791 if 'wave' in kwargs: 1792 wave = kwargs['wave'] 1793 flambda = constants.cc/wave**2 * arg 1794 elif 'freq' in kwargs: 1795 freq = kwargs['freq']/(2*np.pi) 1796 flambda = freq**2/constants.cc * arg 1797 else: 1798 raise ValueError('reference wave/freq not given') 1799 #if constants._current_frequency=='rad': 1800 return flambda*2*np.pi

1801 #else:

1802 # return flambda 1803 1804 -def flambda2fnu(arg,**kwargs):

1805 """ 1806 Switch from Flambda to Fnu via a reference wavelength. 1807 1808 Flambda and Fnu are spectral irradiance in wavelength and frequency, 1809 respectively 1810 1811 @param arg: spectral irradiance (SI, W/m2/m) 1812 @type arg: float 1813 @keyword photband: photometric passband 1814 @type photband: str ('SYSTEM.FILTER') 1815 @keyword wave: reference wavelength (SI, m) 1816 @type wave: float 1817 @keyword freq: reference frequency (SI, Hz) 1818 @type freq: float 1819 @return: spectral irradiance (SI,W/m2/Hz) 1820 @rtype: float 1821 """ 1822 if 'photband' in kwargs: 1823 lameff = filters.eff_wave(kwargs['photband']) 1824 lameff = convert('AA','m',lameff) 1825 kwargs['wave'] = lameff 1826 if 'wave' in kwargs: 1827 wave = kwargs['wave'] 1828 fnu = wave**2/constants.cc * arg 1829 elif 'freq' in kwargs: 1830 freq = kwargs['freq']/(2*np.pi) 1831 fnu = constants.cc/freq**2 * arg 1832 else: 1833 raise ValueError('reference wave/freq not given') 1834 #if constants._current_frequency=='rad': 1835 return fnu/(2*np.pi)

1836 #else:

1837 # return fnu 1838 1839 -def fnu2nufnu(arg,**kwargs):

1840 """ 1841 Switch from Fnu to nuFnu via a reference wavelength. 1842 1843 Flambda and Fnu are spectral irradiance in wavelength and frequency, 1844 respectively 1845 1846 @param arg: spectral irradiance (SI,W/m2/Hz) 1847 @type arg: float 1848 @keyword photband: photometric passband 1849 @type photband: str ('SYSTEM.FILTER') 1850 @keyword wave: reference wavelength (SI, m) 1851 @type wave: float 1852 @keyword freq: reference frequency (SI, Hz) 1853 @type freq: float 1854 @return: spectral irradiance (SI, W/m2/m) 1855 @rtype: float 1856 """ 1857 if 'photband' in kwargs: 1858 lameff = filters.eff_wave(kwargs['photband']) 1859 lameff = convert('AA','m',lameff) 1860 kwargs['wave'] = lameff 1861 if 'wave' in kwargs: 1862 wave = kwargs['wave'] 1863 fnu = constants.cc/wave * arg 1864 elif 'freq' in kwargs: 1865 freq = kwargs['freq'] 1866 fnu = freq * arg 1867 else: 1868 raise ValueError('reference wave/freq not given') 1869 return fnu*(2*np.pi)

1870

1871 -def nufnu2fnu(arg,**kwargs):

1872 """ 1873 Switch from nuFnu to Fnu via a reference wavelength. 1874 1875 Flambda and Fnu are spectral irradiance in wavelength and frequency, 1876 respectively 1877 1878 @param arg: spectral irradiance (SI,W/m2/Hz) 1879 @type arg: float 1880 @keyword photband: photometric passband 1881 @type photband: str ('SYSTEM.FILTER') 1882 @keyword wave: reference wavelength (SI, m) 1883 @type wave: float 1884 @keyword freq: reference frequency (SI, Hz) 1885 @type freq: float 1886 @return: spectral irradiance (SI, W/m2/m) 1887 @rtype: float 1888 """ 1889 if 'photband' in kwargs: 1890 lameff = filters.eff_wave(kwargs['photband']) 1891 lameff = convert('AA','m',lameff) 1892 kwargs['wave'] = lameff 1893 if 'wave' in kwargs: 1894 wave = kwargs['wave'] 1895 fnu = wave/constants.cc * arg 1896 elif 'freq' in kwargs: 1897 freq = kwargs['freq'] 1898 fnu = arg / freq 1899 else: 1900 raise ValueError('reference wave/freq not given') 1901 return fnu/(2*np.pi)

1902

1903 -def flam2lamflam(arg,**kwargs):

1904 """ 1905 Switch from lamFlam to Flam via a reference wavelength. 1906 1907 Flambda and Fnu are spectral irradiance in wavelength and frequency, 1908 respectively 1909 1910 @param arg: spectral irradiance (SI,W/m2/Hz) 1911 @type arg: float 1912 @keyword photband: photometric passband 1913 @type photband: str ('SYSTEM.FILTER') 1914 @keyword wave: reference wavelength (SI, m) 1915 @type wave: float 1916 @keyword freq: reference frequency (SI, Hz) 1917 @type freq: float 1918 @return: spectral irradiance (SI, W/m2/m) 1919 @rtype: float 1920 """ 1921 if 'photband' in kwargs: 1922 lameff = filters.eff_wave(kwargs['photband']) 1923 lameff = convert('AA','m',lameff) 1924 kwargs['wave'] = lameff 1925 if 'wave' in kwargs: 1926 wave = kwargs['wave'] 1927 lamflam = wave * arg 1928 elif 'freq' in kwargs: 1929 freq = kwargs['freq'] 1930 lamflam = constants.cc/freq * arg 1931 else: 1932 raise ValueError('reference wave/freq not given') 1933 return lamflam

1934

1935 -def lamflam2flam(arg,**kwargs):

1936 """ 1937 Switch from lamFlam to Flam via a reference wavelength. 1938 1939 Flambda and Fnu are spectral irradiance in wavelength and frequency, 1940 respectively 1941 1942 @param arg: spectral irradiance (SI,W/m2/Hz) 1943 @type arg: float 1944 @keyword photband: photometric passband 1945 @type photband: str ('SYSTEM.FILTER') 1946 @keyword wave: reference wavelength (SI, m) 1947 @type wave: float 1948 @keyword freq: reference frequency (SI, Hz) 1949 @type freq: float 1950 @return: spectral irradiance (SI, W/m2/m) 1951 @rtype: float 1952 """ 1953 if 'photband' in kwargs: 1954 lameff = filters.eff_wave(kwargs['photband']) 1955 lameff = convert('AA','m',lameff) 1956 kwargs['wave'] = lameff 1957 if 'wave' in kwargs: 1958 wave = kwargs['wave'] 1959 flam = arg / wave 1960 elif 'freq' in kwargs: 1961 freq = kwargs['freq'] 1962 flam = arg / (cc/freq) 1963 else: 1964 raise ValueError('reference wave/freq not given') 1965 return flam

1966

1967 -def distance2spatialfreq(arg,**kwargs):

1968 """ 1969 Switch from distance to spatial frequency via a reference wavelength. 1970 1971 @param arg: distance (SI, m) 1972 @type arg: float 1973 @keyword photband: photometric passband 1974 @type photband: str ('SYSTEM.FILTER') 1975 @keyword wave: reference wavelength (SI, m) 1976 @type wave: float 1977 @keyword freq: reference frequency (SI, Hz) 1978 @type freq: float 1979 @return: spatial frequency (SI, cy/as) 1980 @rtype: float 1981 """ 1982 if 'photband' in kwargs: 1983 lameff = filters.eff_wave(kwargs['photband']) 1984 lameff = convert('AA','m',lameff) 1985 kwargs['wave'] = lameff 1986 if 'wave' in kwargs: 1987 spatfreq = 2*np.pi*arg/kwargs['wave'] 1988 elif 'freq' in kwargs: 1989 spatfreq = 2*np.pi*arg/(constants.cc/kwargs['freq']) 1990 else: 1991 raise ValueError('reference wave/freq not given') 1992 return spatfreq

1993

1994 -def spatialfreq2distance(arg,**kwargs):

1995 """ 1996 Switch from spatial frequency to distance via a reference wavelength. 1997 1998 @param arg: spatial frequency (SI, cy/as) 1999 @type arg: float 2000 @keyword photband: photometric passband 2001 @type photband: str ('SYSTEM.FILTER') 2002 @keyword wave: reference wavelength (SI, m) 2003 @type wave: float 2004 @keyword freq: reference frequency (SI, Hz) 2005 @type freq: float 2006 @return: distance (SI, m) 2007 @rtype: float 2008 """ 2009 if 'photband' in kwargs: 2010 lameff = filters.eff_wave(kwargs['photband']) 2011 lameff = convert('AA','m',lameff) 2012 kwargs['wave'] = lameff 2013 if 'wave' in kwargs: 2014 distance = kwargs['wave']*arg/(2*np.pi) 2015 elif 'freq' in kwargs: 2016 distance = constants.cc/kwargs['freq']*arg/(2*np.pi) 2017 else: 2018 raise ValueError('reference wave/freq not given') 2019 return distance

2020

2021 -def per_sr(arg,**kwargs):

2022 """ 2023 Switch from [Q] to [Q]/sr 2024 2025 @param arg: some SI unit 2026 @type arg: float 2027 @return: some SI unit per steradian 2028 @rtype: float 2029 """ 2030 if 'ang_diam' in kwargs: 2031 radius = kwargs['ang_diam']/2. 2032 surface = np.pi*radius**2 2033 elif 'radius' in kwargs: 2034 radius = kwargs['radius'] 2035 surface = np.pi*radius**2 2036 elif 'pix' in kwargs: 2037 pix = kwargs['pix'] 2038 surface = pix**2 2039 else: 2040 raise ValueError('angular size (ang_diam/radius) not given') 2041 Qsr = arg/surface 2042 return Qsr

2043

2044 -def times_sr(arg,**kwargs):

2045 """ 2046 Switch from [Q]/sr to [Q] 2047 2048 @param arg: some SI unit per steradian 2049 @type arg: float 2050 @return: some SI unit 2051 @rtype: float 2052 """ 2053 if 'ang_diam' in kwargs: 2054 radius = kwargs['ang_diam']/2. 2055 surface = np.pi*radius**2 2056 elif 'radius' in kwargs: 2057 radius = kwargs['radius'] 2058 surface = np.pi*radius**2 2059 elif 'pix' in kwargs: 2060 pix = kwargs['pix'] 2061 surface = pix**2 2062 else: 2063 raise ValueError('angular size (ang_diam/radius) not given') 2064 Q = arg*surface 2065 return Q

2066

2067 -def per_cy(arg,**kwargs):

2068 """ 2069 Switch from radians/s to cycle/s 2070 2071 @rtype: float 2072 """ 2073 return arg / (2*np.pi)

2074

2075 -def times_cy(arg,**kwargs):

2076 """ 2077 Switch from cycle/s to radians/s 2078 2079 @rtype: float 2080 """ 2081 return 2*np.pi*arg

2082

2083 -def period2freq(arg,**kwargs):

2084 """ 2085 Convert period to frequency and back. 2086 2087 @rtype: float 2088 """ 2089 return 1./arg

2090

2091 2092 2093 -def do_nothing(arg,**kwargs):

2094 logger.warning('Experimental: probably just dropped the "cy" unit, please check results') 2095 return arg

2096

2097 2098 -def do_sqrt(arg,**kwargs):

2099 return sqrt(arg)

2100

2101 -def do_quad(arg,**kwargs):

2102 return arg**2

2103

2104 -def Hz2_to_ss(change,frequency):

2105 """ 2106 Convert Frequency shift in Hz2 to period change in seconds per second. 2107 2108 Frequency in Hz! 2109 """ 2110 #-- to second per second 2111 second_per_second = abs(1/frequency - 1/(frequency-change)) 2112 return second_per_second

2113

2114 #} 2115 2116 2117 #{ Stellar calibrations 2118 2119 -def derive_luminosity(radius,temperature,units=None):

2120 """ 2121 Convert radius and effective temperature to stellar luminosity. 2122 2123 Units given to radius and temperature must be understandable by C{Unit}. 2124 2125 Stellar luminosity is returned, by default in Solar units. 2126 2127 I recommend to give the input as follows, since this is most clear also 2128 from a programmatical point of view: 2129 2130 >>> print(derive_luminosity((1.,'Rsol'),(5777,'K'),units='erg/s')) 2131 3.83916979644e+33 erg/s 2132 2133 The function is pretty smart, however: it knows what the output units 2134 should be, so you could ask to put them in the 'standard' units of some 2135 convention: 2136 2137 >>> print(derive_luminosity((1.,'Rsol'),(5777,'K'),units='cgs')) 2138 3.83916979644e+33 g1 cm2 s-3 2139 2140 If you give nothing, everything will be interpreted in the current 2141 convention's units: 2142 2143 >>> print(derive_luminosity(7e8,5777.)) 2144 3.88892117726e+26 W1 2145 2146 You are also free to give units (you can do this in a number of ways), and 2147 if you index the return value, you can get the value ([0], float or 2148 uncertainty) or unit ([1],string): 2149 2150 >>> derive_luminosity((1.,'Rsol'),(1.,0.1,'Tsol'),units='Lsol')[0] 2151 1.0000048246799305+/-0.4000019298719723 2152 2153 >>> derive_luminosity((1.,'Rsol'),((1.,0.1),'Tsol'),units='Lsol')[0] 2154 1.0000048246799305+/-0.4000019298719723 2155 2156 Finally, you can simply work with Units: 2157 2158 >>> radius = Unit(1.,'Rsol') 2159 >>> temperature = Unit(5777.,'K') 2160 >>> print(derive_luminosity(radius,temperature,units='foe/s')) 2161 3.83916979644e-18 foe/s 2162 2163 Everything should be independent of the current convention, unless it 2164 needs to be: 2165 2166 >>> old = set_convention('cgs') 2167 >>> derive_luminosity((1.,'Rsol'),(1.,0.1,'Tsol'),units='Lsol')[0] 2168 1.0000048246799305+/-0.4000019298719722 2169 >>> print(derive_luminosity(7e10,5777.)) 2170 3.88892117726e+33 erg1 s-1 2171 >>> sol = set_convention('sol') 2172 >>> print(derive_luminosity(1.,1.)) 2173 3.9982620567e-22 Msol1 Rsol2 s-3 2174 >>> print(derive_luminosity(1.,1.,units='Lsol')) 2175 1.00000482468 Lsol 2176 >>> old = set_convention(*old) 2177 2178 @param radius: (radius(, error), units) 2179 @type radius: 2 or 3 tuple, Unit or float (current convention) 2180 @param temperature: (effective temperature(, error), units) 2181 @type temperature: 2 or 3 tuple, Unit or float (current convention) 2182 @return: luminosity 2183 @rtype: Unit 2184 """ 2185 if units is None: 2186 units = constants._current_convention 2187 radius = Unit(radius,unit='length') 2188 temperature = Unit(temperature,unit='temperature') 2189 sigma = Unit('sigma') 2190 2191 lumi = 4*np.pi*sigma*radius**2*temperature**4 2192 2193 return lumi.convert(units)

2194

2195 2196 -def derive_radius(luminosity,temperature, units='m'):

2197 """ 2198 Convert luminosity and effective temperature to stellar radius. 2199 2200 Units given to luminosity and temperature must be understandable by C{convert}. 2201 2202 Stellar radius is returned in SI units. 2203 2204 Example usage: 2205 2206 >>> print(derive_radius((3.9,'[Lsol]'),(3.72,'[K]'),units='Rsol')) 2207 108.091293736 Rsol 2208 >>> print(derive_radius(3.8e26,5777.,units='Rjup')) 2209 9.89759670363 Rjup 2210 2211 @param luminosity: (Luminosity(, error), units) 2212 @type luminosity: 2 or 3 tuple 2213 @param temperature: (effective temperature(, error), units) 2214 @type temperature: 2 or 3 tuple, Unit or float (current convention) 2215 @return: radius 2216 @rtype: Unit 2217 """ 2218 if units is None: 2219 units = constants._current_convention 2220 luminosity = Unit(luminosity,unit='power') 2221 temperature = Unit(temperature,unit='temperature') 2222 sigma = Unit('sigma') 2223 2224 R = np.sqrt(luminosity / (temperature**4)/(4*np.pi*sigma)) 2225 2226 return R.convert(units)

2227

2228 -def derive_radius_slo(numax,Deltanu0,teff,unit='Rsol'):

2229 """ 2230 Derive stellar radius from solar-like oscillations diagnostics. 2231 2232 Large separation is frequency spacing between modes of different radial order. 2233 2234 Small separation is spacing between modes of even and odd angular degree but 2235 same radial order. 2236 2237 @param numax: (numax(, error), units) 2238 @type numax: 2 or 3 tuple 2239 @param Deltanu0: (large separation(, error), units) 2240 @type Deltanu0: 2 or 3 tuple 2241 @param teff: (effective temperature(, error), units) 2242 @type teff: 2 or 3 tuple 2243 @return: radius (and error) in whatever units 2244 @rtype: 1- or 2-tuple 2245 """ 2246 numax_sol = convert(constants.numax_sol[-1],'mHz',*constants.numax_sol[:-1],unpack=False) 2247 Deltanu0_sol = convert(constants.Deltanu0_sol[-1],'muHz',*constants.Deltanu0_sol[:-1],unpack=False) 2248 #-- take care of temperature 2249 if len(teff)==3: 2250 teff = (unumpy.uarray([teff[0],teff[1]]),teff[2]) 2251 teff = convert(teff[-1],'5777K',*teff[:-1],unpack=False) 2252 #-- take care of numax 2253 if len(numax)==3: 2254 numax = (unumpy.uarray([numax[0],numax[1]]),numax[2]) 2255 numax = convert(numax[-1],'mHz',*numax[:-1],unpack=False) 2256 #-- take care of large separation 2257 if len(Deltanu0)==3: 2258 Deltanu0 = (unumpy.uarray([Deltanu0[0],Deltanu0[1]]),Deltanu0[2]) 2259 Deltanu0 = convert(Deltanu0[-1],'muHz',*Deltanu0[:-1],unpack=False) 2260 R = sqrt(teff)/numax_sol * numax/Deltanu0**2 * Deltanu0_sol**2 2261 return convert('Rsol',unit,R)

2262

2263 #@derive_wrapper 2264 -def derive_radius_rot(veq_sini,rot_period,incl=(90.,'deg'),unit='Rsol'):

2265 """ 2266 Derive radius from rotation period and inclination angle. 2267 """ 2268 incl = Unit(incl) 2269 radius = veq_sini*rot_period/(2*np.pi)/np.sin(incl) 2270 return radius

2271

2272 2273 2274 2275 -def derive_logg(mass,radius, unit='[cm/s2]'):

2276 """ 2277 Convert mass and radius to stellar surface gravity. 2278 2279 Units given to mass and radius must be understandable by C{convert}. 2280 2281 Logarithm of surface gravity is returned in CGS units. 2282 2283 @param mass: (mass(, error), units) 2284 @type mass: 2 or 3 tuple 2285 @param radius: (radius(, error), units) 2286 @type radius: 2 or 3 tuple 2287 @return: log g (and error) in CGS units 2288 @rtype: 1- or 2-tuple 2289 """ 2290 #-- take care of mass 2291 if len(mass)==3: 2292 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2293 M = convert(mass[-1],'g',*mass[:-1],unpack=False) 2294 #-- take care of radius 2295 if len(radius)==3: 2296 radius = (unumpy.uarray([radius[0],radius[1]]),radius[2]) 2297 R = convert(radius[-1],'cm',*radius[:-1],unpack=False) 2298 #-- calculate surface gravity in logarithmic CGS units 2299 logg = log10(constants.GG_cgs*M / (R**2)) 2300 logg = convert('[cm/s2]',unit,logg) 2301 return logg

2302

2303 2304 -def derive_logg_zg(mass, zg, unit='cm s-2', **kwargs):

2305 """ 2306 Convert mass and gravitational redshift to stellar surface gravity. Provide the gravitational 2307 redshift as a velocity. 2308 2309 Units given to mass and zg must be understandable by C{convert}. 2310 2311 Logarithm of surface gravity is returned in CGS units unless otherwhise stated in unit. 2312 2313 >>> print derive_logg_zg((0.47, 'Msol'), (2.57, 'km/s')) 2314 5.97849814386 2315 2316 @param mass: (mass(, error), units) 2317 @type mass: 2 or 3 tuple 2318 @param zg: (zg(, error), units) 2319 @type zg: 2 or 3 tuple 2320 @param unit: unit of logg 2321 @type unit: str 2322 @return: log g (and error) 2323 @rtype: 1- or 2-tuple 2324 """ 2325 2326 #-- take care of mass 2327 if len(mass)==3: 2328 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2329 M = convert(mass[-1],'g',*mass[:-1],unpack=False) 2330 #-- take care of zg 2331 if len(zg)==3: 2332 zg = (unumpy.uarray([zg[0],zg[1]]),zg[2]) 2333 gr = convert(zg[-1],'cm s-1',*zg[:-1],unpack=False, **kwargs) 2334 2335 logg = np.log10( gr**2 * constants.cc_cgs**2 / (constants.GG_cgs * M) ) 2336 return convert('cm s-2', unit, logg)

2337

2338 -def derive_logg_slo(teff,numax, unit='[cm/s2]'):

2339 """ 2340 Derive stellar surface gravity from solar-like oscillations diagnostics. 2341 2342 Units given to teff and numax must be understandable by C{convert}. 2343 2344 Logarithm of surface gravity is returned in CGS units. 2345 2346 @param teff: (effective temperature(, error), units) 2347 @type teff: 2 or 3 tuple 2348 @param numax: (numax(, error), units) 2349 @type numax: 2 or 3 tuple 2350 @return: log g (and error) in CGS units 2351 @rtype: 1- or 2-tuple 2352 """ 2353 numax_sol = convert(constants.numax_sol[-1],'mHz',*constants.numax_sol[:-1],unpack=False) 2354 #-- take care of temperature 2355 if len(teff)==3: 2356 teff = (unumpy.uarray([teff[0],teff[1]]),teff[2]) 2357 teff = convert(teff[-1],'5777K',*teff[:-1],unpack=False) 2358 #-- take care of numax 2359 if len(numax)==3: 2360 numax = (unumpy.uarray([numax[0],numax[1]]),numax[2]) 2361 numax = convert(numax[-1],'mHz',*numax[:-1],unpack=False) 2362 #-- calculate surface gravity in logarithmic CGS units 2363 GG = convert(constants.GG_units,'Rsol3 Msol-1 s-2',constants.GG) 2364 surf_grav = GG*sqrt(teff)*numax / numax_sol 2365 logg = convert('Rsol s-2',unit,surf_grav) 2366 return logg

2367

2368 -def derive_mass(surface_gravity,radius,unit='kg'):

2369 """ 2370 Convert surface gravity and radius to stellar mass. 2371 """ 2372 #-- take care of logg 2373 if len(surface_gravity)==3: 2374 surface_gravity = (unumpy.uarray([surface_gravity[0],surface_gravity[1]]),surface_gravity[2]) 2375 grav = convert(surface_gravity[-1],'m/s2',*surface_gravity[:-1],unpack=False) 2376 #-- take care of radius 2377 if len(radius)==3: 2378 radius = (unumpy.uarray([radius[0],radius[1]]),radius[2]) 2379 R = convert(radius[-1],'m',*radius[:-1],unpack=False) 2380 #-- calculate mass in SI 2381 M = grav*R**2/constants.GG 2382 return convert('kg',unit,M)

2383

2384 -def derive_zg(mass, logg, unit='cm s-1', **kwargs):

2385 """ 2386 Convert mass and stellar surface gravity to gravitational redshift. 2387 2388 Units given to mass and logg must be understandable by C{convert}. 2389 2390 Gravitational redshift is returned in CGS units unless otherwhise stated in unit. 2391 2392 >>> print derive_zg((0.47, 'Msol'), (5.98, 'cm s-2'), unit='km s-1') 2393 2.57444756874 2394 2395 @param mass: (mass(, error), units) 2396 @type mass: 2 or 3 tuple 2397 @param logg: (logg(, error), units) 2398 @type logg: 2 or 3 tuple 2399 @param unit: unit of zg 2400 @type unit: str 2401 @return: log g (and error) 2402 @rtype: 1- or 2-tuple 2403 """ 2404 #-- take care of mass 2405 if len(mass)==3: 2406 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2407 M = convert(mass[-1],'g',*mass[:-1],unpack=False) 2408 #-- take care of logg 2409 if len(logg)==3: 2410 logg = (unumpy.uarray([logg[0],logg[1]]),logg[2]) 2411 g = convert(logg[-1],'cm s-2',*logg[:-1],unpack=False, **kwargs) 2412 2413 zg = 10**g / constants.cc_cgs * np.sqrt(constants.GG_cgs * M / 10**g) 2414 return convert('cm s-1', unit, zg)

2415

2416 -def derive_numax(mass,radius,temperature,unit='mHz'):

2417 """ 2418 Derive the predicted nu_max according to Kjeldsen and Bedding (1995). 2419 2420 Example: compute the predicted numax for the Sun in mHz 2421 >>> print derive_numax((1.,'Msol'),(1.,'Rsol'),(5777.,'K')) 2422 3.05 2423 """ 2424 #-- take care of mass 2425 if len(mass)==3: 2426 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2427 M = convert(mass[-1],'Msol',*mass[:-1],unpack=False) 2428 #-- take care of radius 2429 if len(radius)==3: 2430 radius = (unumpy.uarray([radius[0],radius[1]]),radius[2]) 2431 R = convert(radius[-1],'Rsol',*radius[:-1],unpack=False) 2432 #-- take care of effective temperature 2433 if len(temperature)==3: 2434 temperature = (unumpy.uarray([temperature[0],temperature[1]]),temperature[2]) 2435 teff = convert(temperature[-1],'5777K',*temperature[:-1],unpack=False) 2436 #-- predict nu_max 2437 nu_max = M/R**2/sqrt(teff)*3.05 2438 return convert('mHz',unit,nu_max)

2439

2440 -def derive_nmax(mass,radius,temperature):

2441 """ 2442 Derive the predicted n_max according to Kjeldsen and Bedding (1995). 2443 2444 Example: compute the predicted numax for the Sun in mHz 2445 >>> print derive_nmax((1.,'Msol'),(1.,'Rsol'),(5777.,'K')) 2446 21.0 2447 """ 2448 #-- take care of mass 2449 if len(mass)==3: 2450 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2451 M = convert(mass[-1],'Msol',*mass[:-1]) 2452 #-- take care of radius 2453 if len(radius)==3: 2454 radius = (unumpy.uarray([radius[0],radius[1]]),radius[2]) 2455 R = convert(radius[-1],'Rsol',*radius[:-1]) 2456 #-- take care of effective temperature 2457 if len(temperature)==3: 2458 temperature = unumpy.uarray([temperature[0],temperature[1]]) 2459 teff = convert(temperature[-1],'5777K',*temperature[:-1]) 2460 #-- predict n_max 2461 n_max = sqrt(M/teff/R)*22.6 - 1.6 2462 return n_max

2463

2464 -def derive_Deltanu0(mass,radius,unit='mHz'):

2465 """ 2466 Derive the predicted large spacing according to Kjeldsen and Bedding (1995). 2467 2468 Example: compute the predicted large spacing for the Sun in mHz 2469 >>> print derive_Deltanu0((1.,'Msol'),(1.,'Rsol'),unit='muHz') 2470 134.9 2471 """ 2472 #-- take care of mass 2473 if len(mass)==3: 2474 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2475 M = convert(mass[-1],'Msol',*mass[:-1]) 2476 #-- take care of radius 2477 if len(radius)==3: 2478 radius = (unumpy.uarray([radius[0],radius[1]]),radius[2]) 2479 R = convert(radius[-1],'Rsol',*radius[:-1]) 2480 #-- predict large spacing 2481 Deltanu0 = sqrt(M)*R**(-1.5)*134.9 2482 return convert('muHz',unit,Deltanu0)

2483

2484 -def derive_ampllum(luminosity,mass,temperature,wavelength,unit='ppm'):

2485 """ 2486 Derive the luminosity amplitude around nu_max of solar-like oscillations. 2487 2488 See Kjeldsen and Bedding (1995). 2489 2490 >>> print derive_ampllum((1.,'Lsol'),(1,'Msol'),(5777.,'K'),(550,'nm')) 2491 (4.7, 0.3) 2492 """ 2493 #-- take care of mass 2494 if len(mass)==3: 2495 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2496 M = convert(mass[-1],'Msol',*mass[:-1]) 2497 #-- take care of effective temperature 2498 if len(temperature)==3: 2499 temperature = unumpy.uarray([temperature[0],temperature[1]]) 2500 teff = convert(temperature[-1],'5777K',*temperature[:-1]) 2501 #-- take care of luminosity 2502 if len(luminosity)==3: 2503 luminosity = unumpy.uarray([luminosity[0],luminosity[1]]) 2504 lumi = convert(luminosity[-1],'Lsol',*luminosity[:-1]) 2505 #-- take care of wavelength 2506 if len(wavelength)==3: 2507 wavelength = unumpy.uarray([wavelength[0],wavelength[1]]) 2508 wave = convert(wavelength[-1],'550nm',*wavelength[:-1]) 2509 ampllum = lumi / wave / teff**2 / M * ufloat((4.7,0.3)) 2510 return convert('ppm',unit,ampllum)

2511

2512 -def derive_ampllum_from_velo(velo,temperature,wavelength,unit='ppm'):

2513 """ 2514 Derive the luminosity amplitude around nu_max of solar-like oscillations. 2515 2516 See Kjeldsen and Bedding (1995). 2517 2518 >>> print derive_ampllum_from_velo((1.,'m/s'),(5777.,'K'),(550,'nm')) 2519 (20.1, 0.05) 2520 """ 2521 #-- take care of velocity 2522 if len(velo)==3: 2523 velo = (unumpy.uarray([velo[0],velo[1]]),velo[2]) 2524 velo = convert(velo[-1],'m/s',*velo[:-1]) 2525 #-- take care of effective temperature 2526 if len(temperature)==3: 2527 temperature = unumpy.uarray([temperature[0],temperature[1]]) 2528 teff = convert(temperature[-1],'5777K',*temperature[:-1]) 2529 #-- take care of wavelength 2530 if len(wavelength)==3: 2531 wavelength = unumpy.uarray([wavelength[0],wavelength[1]]) 2532 wave = convert(wavelength[-1],'550nm',*wavelength[:-1]) 2533 ampllum = velo / wave / teff**2 * ufloat((20.1,0.05)) #not sure about error 2534 return convert('ppm',unit,ampllum)

2535

2536 -def derive_amplvel(luminosity,mass,unit='cm/s'):

2537 """ 2538 Derive the luminosity amplitude around nu_max of solar-like oscillations. 2539 2540 See Kjeldsen and Bedding (1995). 2541 2542 >>> print derive_amplvel((1.,'Lsol'),(1,'Msol'),unit='cm/s') 2543 (23.4, 1.4) 2544 """ 2545 #-- take care of mass 2546 if len(mass)==3: 2547 mass = (unumpy.uarray([mass[0],mass[1]]),mass[2]) 2548 M = convert(mass[-1],'Msol',*mass[:-1]) 2549 #-- take care of luminosity 2550 if len(luminosity)==3: 2551 luminosity = unumpy.uarray([luminosity[0],luminosity[1]]) 2552 lumi = convert(luminosity[-1],'Lsol',*luminosity[:-1]) 2553 amplvel = lumi / M * ufloat((23.4,1.4)) 2554 return convert('cm/s',unit,amplvel)

2555

2556 -def derive_galactic_uvw(ra, dec, pmra, pmdec, d, vrad, lsr=False, unit='km s-1'):

2557 """ 2558 Calculate the Galactic space velocity (U,V,W) of a star based on the propper motion, 2559 location, radial velocity and distance. (U,V,W) are returned in km/s unless stated 2560 otherwhise in unit. 2561 2562 Follows the general outline of Johnson & Soderblom (1987, AJ, 93,864) 2563 except that U is positive outward toward the Galactic *anti*center, and 2564 the J2000 transformation matrix to Galactic coordinates is taken from 2565 the introduction to the Hipparcos catalog. 2566 2567 Uses solar motion from Coskunoglu et al. 2011 MNRAS: 2568 (U,V,W)_Sun = (-8.5, 13.38, 6.49) 2569 2570 Example for HD 6755: 2571 2572 >>> derive_galactic_uvw((17.42943586, 'deg'), (61.54727506, 'deg'), (628.42, 'mas yr-1'), (76.65, 'mas yr-1'), (139, 'pc'), (-321.4, 'km s-1'), lsr=True) 2573 (142.66328352779027, -483.55149105148121, 93.216106970932813) 2574 2575 @param ra: Right assension of the star 2576 @param dec: Declination of the star 2577 @param pmra: propper motion in RA 2578 @param pmdec: propper motion in DEC 2579 @param d: distance 2580 @param vrad: radial velocity 2581 @param lsr: if True, correct for solar motion 2582 @param unit: units for U,V and W 2583 2584 @return: (U,V,W) 2585 """ 2586 #-- take care of ra 2587 if len(ra)==3: 2588 ra = (unumpy.uarray([ra[0],ra[1]]),ra[2]) 2589 ra = convert(ra[-1],'deg',*ra[:-1]) 2590 #-- take care of dec 2591 if len(dec)==3: 2592 dec = (unumpy.uarray([dec[0],dec[1]]),dec[2]) 2593 dec = convert(dec[-1],'deg',*dec[:-1]) 2594 #-- take care of pmra 2595 if len(pmra)==3: 2596 pmra = (unumpy.uarray([pmra[0],pmra[1]]),pmra[2]) 2597 pmra = convert(pmra[-1],'mas yr-1',*pmra[:-1]) 2598 #-- take care of pmdec 2599 if len(pmdec)==3: 2600 pmdec = (unumpy.uarray([pmdec[0],pmdec[1]]),pmdec[2]) 2601 pmdec = convert(pmdec[-1],'mas yr-1',*pmdec[:-1]) 2602 #-- take care of d 2603 if len(d)==3: 2604 d = (unumpy.uarray([d[0],d[1]]),d[2]) 2605 d = convert(d[-1],'pc',*d[:-1]) 2606 #-- take care of pmdec 2607 if len(vrad)==3: 2608 vrad = (unumpy.uarray([vrad[0],vrad[1]]),vrad[2]) 2609 vrad = convert(vrad[-1],'km s-1',*vrad[:-1]) 2610 2611 plx = 1e3 / d # parallax in mas 2612 2613 cosd = np.cos(np.radians(dec)) 2614 sind = np.sin(np.radians(dec)) 2615 cosa = np.cos(np.radians(ra)) 2616 sina = np.sin(np.radians(ra)) 2617 2618 k = 4.74047 #Equivalent of 1 A.U/yr in km/s 2619 A_G = np.array( [ [ 0.0548755604, +0.4941094279, -0.8676661490], 2620 [ 0.8734370902, -0.4448296300, -0.1980763734], 2621 [ 0.4838350155, 0.7469822445, +0.4559837762] ]).T 2622 2623 vec1 = vrad 2624 vec2 = k * pmra / plx 2625 vec3 = k * pmdec / plx 2626 2627 u = ( A_G[0,0]*cosa*cosd+A_G[0,1]*sina*cosd+A_G[0,2]*sind)*vec1+ \ 2628 (-A_G[0,0]*sina +A_G[0,1]*cosa )*vec2+ \ 2629 (-A_G[0,0]*cosa*sind-A_G[0,1]*sina*sind+A_G[0,2]*cosd)*vec3 2630 v = ( A_G[1,0]*cosa*cosd+A_G[1,1]*sina*cosd+A_G[1,2]*sind)*vec1+ \ 2631 (-A_G[1,0]*sina +A_G[1,1]*cosa )*vec2+ \ 2632 (-A_G[1,0]*cosa*sind-A_G[1,1]*sina*sind+A_G[1,2]*cosd)*vec3 2633 w = ( A_G[2,0]*cosa*cosd+A_G[2,1]*sina*cosd+A_G[2,2]*sind)*vec1+ \ 2634 (-A_G[2,0]*sina +A_G[2,1]*cosa )*vec2+ \ 2635 (-A_G[2,0]*cosa*sind-A_G[2,1]*sina*sind+A_G[2,2]*cosd)*vec3 2636 2637 if lsr: 2638 #lsr_vel=[-10.0,5.2,7.2] # Dehnen & Binney (1998) 2639 #lsr_vel=[-11.1, 12.24, 7.25] # Schonrich et al. (2010) 2640 lsr_vel=[-8.5,13.38,6.49] # Coskunoglu et al. (2011) 2641 u = u+lsr_vel[0] 2642 v = v+lsr_vel[1] 2643 w = w+lsr_vel[2] 2644 2645 return convert('km s-1',unit,u), convert('km s-1',unit,v), convert('km s-1',unit,w)

2646

2647 -def convert_Z_FeH(Z=None, FeH=None, Zsun=0.0122):

2648 """ 2649 Convert between Z and [Fe/H] using the conversion of Bertelli et al. (1994). 2650 Provide one of the 2 arguments, and the other will be returned. 2651 2652 log10(Z/Zsun) = 0.977 [Fe/H] 2653 """ 2654 2655 if Z != None: 2656 return np.log10(Z/Zsun) / 0.977 2657 else: 2658 return 10**(0.977*FeH + np.log10(Zsun))

2659

2660 2661 #} 2662 2663 2664 2665 #{ Nonlinear change-of-base functions 2666 2667 -class NonLinearConverter():

2668 """ 2669 Base class for nonlinear conversions 2670 2671 This class keeps track of prefix-factors and powers. 2672 2673 To have a real nonlinear converter, you need to define the C{__call__} 2674 attribute. 2675 """

2676 - def __init__(self,prefix=1.,power=1.):

2677 self.prefix = prefix 2678 self.power = power

2679 - def __rmul__(self,other):

2680 if type(other)==type(5) or type(other)==type(5.): 2681 return self.__class__(prefix=self.prefix*other)

2682 - def __div__(self,other):

2683 if type(other)==type(5) or type(other)==type(5.): 2684 return self.__class__(prefix=self.prefix*other)

2685 - def __pow__(self,other):

2686 if type(other)==type(5) or type(other)==type(5.): 2687 return self.__class__(prefix=self.prefix,power=self.power+other)

2688

2689 -class Fahrenheit(NonLinearConverter):

2690 """ 2691 Convert Fahrenheit to Kelvin and back 2692 """

2693 - def __call__(self,a,inv=False):

2694 if not inv: return (a*self.prefix+459.67)*5./9. 2695 else: return (a*9./5.-459.67)/self.prefix

2696

2697 -class Celcius(NonLinearConverter):

2698 """ 2699 Convert Celcius to Kelvin and back 2700 """

2701 - def __call__(self,a,inv=False):

2702 if not inv: return a*self.prefix+273.15 2703 else: return (a-273.15)/self.prefix

2704

2705 2706 -class AmplMag(NonLinearConverter):

2707 """ 2708 Convert a Vega magnitude to W/m2/m (Flambda) and back 2709 """

2710 - def __call__(self,meas,inv=False):

2711 #-- this part should include something where the zero-flux is retrieved 2712 if not inv: return 10**(meas*self.prefix/2.5) - 1. 2713 else: return (2.5*log10(1.+meas))/self.prefix

2714

2715 -class VegaMag(NonLinearConverter):

2716 """ 2717 Convert a Vega magnitude to W/m2/m (Flambda) and back 2718 """

2719 - def __call__(self,meas,photband=None,inv=False,**kwargs):

2720 #-- this part should include something where the zero-flux is retrieved 2721 zp = filters.get_info() 2722 match = zp['photband']==photband.upper() 2723 if sum(match)==0: raise ValueError("No calibrations for %s"%(photband)) 2724 F0 = convert(zp['Flam0_units'][match][0],'W/m3',zp['Flam0'][match][0]) 2725 mag0 = float(zp['vegamag'][match][0]) 2726 if not inv: return 10**(-(meas-mag0)/2.5)*F0 2727 else: return -2.5*log10(meas/F0)+mag0

2728

2729 -class ABMag(NonLinearConverter):

2730 """ 2731 Convert an AB magnitude to W/m2/Hz (Fnu) and back 2732 """

2733 - def __call__(self,meas,photband=None,inv=False,**kwargs):

2734 zp = filters.get_info() 2735 F0 = convert('W/m2/Hz',constants._current_convention,3.6307805477010024e-23) 2736 match = zp['photband']==photband.upper() 2737 if sum(match)==0: raise ValueError("No calibrations for %s"%(photband)) 2738 mag0 = float(zp['ABmag'][match][0]) 2739 if np.isnan(mag0): mag0 = 0. 2740 if not inv: 2741 try: 2742 return 10**(-(meas-mag0)/2.5)*F0 2743 except OverflowError: 2744 return np.nan 2745 else: return -2.5*log10(meas/F0)

2746

2747 -class STMag(NonLinearConverter):

2748 """ 2749 Convert an ST magnitude to W/m2/m (Flambda) and back 2750 2751 mag = -2.5*log10(F) - 21.10 2752 2753 F0 = 3.6307805477010028e-09 erg/s/cm2/AA 2754 """

2755 - def __call__(self,meas,photband=None,inv=False,**kwargs):

2756 zp = filters.get_info() 2757 F0 = convert('erg/s/cm2/AA',constants._current_convention,3.6307805477010028e-09)#0.036307805477010027 2758 match = zp['photband']==photband.upper() 2759 if sum(match)==0: raise ValueError("No calibrations for %s"%(photband)) 2760 mag0 = float(zp['STmag'][match][0]) 2761 if np.isnan(mag0): mag0 = 0. 2762 if not inv: return 10**(-(meas-mag0)/-2.5)*F0 2763 else: return -2.5*log10(meas/F0)

2764

2765 2766 -class Color(NonLinearConverter):

2767 """ 2768 Convert a color to a flux ratio and back 2769 2770 B-V = -2.5log10(FB) + CB - (-2.5log10(FV) + CV) 2771 B-V = -2.5log10(FB) + CB + 2.5log10(FV) - CV 2772 B-V = -2.5log10(FB/FV) + (CB-CV) 2773 2774 and thus 2775 2776 FB/FV = 10 ** [((B-V) - (CB-CV)) / (-2.5)] 2777 2778 where 2779 2780 CB = 2.5log10[FB(m=0)] 2781 CV = 2.5log10[FV(m=0)] 2782 2783 Stromgren colour indices: 2784 2785 m1 = v - 2b + y 2786 c1 = u - 2v + b 2787 Hbeta = HBN - HBW 2788 """

2789 - def __call__(self,meas,photband=None,inv=False,**kwargs):

2790 #-- we have two types of colours: the stromgren M1/C1 type, and the 2791 # normal Band1 - Band2 type. We need to have conversions back and 2792 # forth: this translates into four cases. 2793 system,band = photband.split('.') 2794 if '-' in band and not inv: 2795 band0,band1 = band.split('-') 2796 f0 = convert('mag','SI',meas,photband='.'.join([system,band0]),unpack=False) 2797 f1 = convert('mag','SI',0.00,photband='.'.join([system,band1])) 2798 return f0/f1 2799 elif '-' in band and inv: 2800 #-- the units don't really matter, we choose SI' 2801 # the flux ratio is converted to color by assuming that the 2802 # denominator flux is equal to one. 2803 band0,band1 = band.split('-') 2804 m0 = convert('W/m3','mag',meas,photband='.'.join([system,band0]),unpack=False) 2805 m1 = convert('W/m3','mag',1.00,photband='.'.join([system,band1])) 2806 return m0-m1 2807 elif photband=='STROMGREN.C1' and not inv: 2808 fu = convert('mag','SI',meas,photband='STROMGREN.U',unpack=False) 2809 fb = convert('mag','SI',0.00,photband='STROMGREN.B') 2810 fv = convert('mag','SI',0.00,photband='STROMGREN.V') 2811 return fu*fb/fv**2 2812 elif photband=='STROMGREN.C1' and inv: 2813 mu = convert('W/m3','mag',meas,photband='STROMGREN.U',unpack=False) 2814 mb = convert('W/m3','mag',1.00,photband='STROMGREN.B') 2815 mv = convert('W/m3','mag',1.00,photband='STROMGREN.V') 2816 return mu-2*mv+mb 2817 elif photband=='STROMGREN.M1' and not inv: 2818 fv = convert('mag','SI',meas,photband='STROMGREN.V',unpack=False) 2819 fy = convert('mag','SI',0.00,photband='STROMGREN.Y') 2820 fb = convert('mag','SI',0.00,photband='STROMGREN.B') 2821 return fv*fy/fb**2 2822 elif photband=='STROMGREN.M1' and inv: 2823 mu = convert('W/m3','mag',meas,photband='STROMGREN.V',unpack=False) 2824 mb = convert('W/m3','mag',1.00,photband='STROMGREN.Y') 2825 mv = convert('W/m3','mag',1.00,photband='STROMGREN.B') 2826 return mv-2*mb+my 2827 else: 2828 raise ValueError("No color calibrations for %s"%(photband))

2829

2830 -class DecibelSPL(NonLinearConverter):

2831 """ 2832 Convert a Decibel to intensity W/m2 and back. 2833 2834 This is only valid for Decibels as a sound Pressure level 2835 """

2836 - def __call__(self,meas,inv=False):

2837 F0 = 1e-12 # W/m2 2838 if not inv: return 10**(-meas)*F0 2839 else: return log10(meas/F0)

2840

2841 -class JulianDay(NonLinearConverter):

2842 """ 2843 Convert a calender date to Julian date and back 2844 """

2845 - def __call__(self,meas,inv=False,**kwargs):

2846 if inv: 2847 #L= meas+68569 2848 #N= 4*L//146097 2849 #L= L-(146097*N+3)//4 2850 #I= 4000*(L+1)//1461001 2851 #L= L-1461*I//4+31 2852 #J= 80*L//2447 2853 #day = L-2447*J//80+0.5 2854 #L= J//11 2855 #month = J+2-12*L 2856 #year = 100*(N-49)+I+L 2857 2858 j = meas + 0.5 + 32044 2859 g = np.floor(j/146097) 2860 dg = np.fmod(j,146097) 2861 c = np.floor((np.floor(dg/36524) + 1)*3/4) 2862 dc = dg - c * 36524 2863 b = np.floor(dc/1461) 2864 db = np.fmod(dc,1461) 2865 a = np.floor((np.floor(db/365) + 1)*3/4) 2866 da = db - a * 365 2867 y = g*400+c*100+b*4+a 2868 m = np.floor((da*5+308)/153)-2 2869 d = da - np.floor((m+4)*153/5) + 124 2870 year = y - 4800 + np.floor((m+2)/12) 2871 month = np.fmod(m+2,12)+1 2872 day = d - 1 2873 return year,month,day 2874 else: 2875 year,month,day = meas[:3] 2876 hour = len(meas)>3 and meas[3] or 0. 2877 mint = len(meas)>4 and meas[4] or 0. 2878 secn = len(meas)>5 and meas[5] or 0. 2879 a = (14 - month)//12 2880 y = year + 4800 - a 2881 m = month + 12*a - 3 2882 jd = day + ((153*m + 2)//5) + 365*y + y//4 - y//100 + y//400 - 32045 2883 jd += hour/24. 2884 jd += mint/24./60. 2885 jd += secn/24./3600. 2886 jd -= 0.5 2887 return jd

2888

2889 -class ModJulianDay(NonLinearConverter):

2890 """ 2891 Convert a Modified Julian Day to Julian Day and back 2892 2893 The CoRoT conversion has been checked with the CoRoT data archive: it is 2894 correct at least to the second (the archive tool's precision). 2895 """ 2896 ZP = {'COROT':2451545., 2897 'HIP':2440000., 2898 'MJD':2400000.5, 2899 'PYEPHEM':2415020.0}

2900 - def __call__(self,meas,inv=False,jtype='MJD'):

2901 if inv: 2902 return meas-self.ZP[jtype.upper()] 2903 else: 2904 return meas+self.ZP[jtype.upper()]

2905

2906 -class GalCoords(NonLinearConverter):

2907 """ 2908 Convert Galactic coords to complex coords and back 2909 """

2910 - def __call__(self,mycoord,inv=False,epoch='2000'):

2911 if inv: 2912 x,y = mycoord.real,mycoord.imag 2913 equ = ephem.Equatorial(x,y,epoch=epoch) 2914 gal = ephem.Galactic(equ,epoch=epoch) 2915 return gal.lon,gal.lat 2916 else: 2917 x,y = mycoord 2918 gal = ephem.Galactic(x,y,epoch=epoch) 2919 equ = ephem.Equatorial(gal,epoch=epoch) 2920 return float(equ.ra) + 1j*float(equ.dec)

2921

2922 -class EquCoords(NonLinearConverter):

2923 """ 2924 Convert Equatorial coords to complex coords and back 2925 """

2926 - def __call__(self,mycoord,inv=False,epoch='2000'):

2927 if inv: 2928 x,y = mycoord.real,mycoord.imag 2929 equ = ephem.Equatorial(x,y,epoch=epoch) 2930 return equ.ra,equ.dec 2931 else: 2932 x,y = mycoord 2933 equ = ephem.Equatorial(x,y,epoch=epoch) 2934 return float(equ.ra) + 1j*float(equ.dec)

2935

2936 -class EclCoords(NonLinearConverter):

2937 """ 2938 Convert Ecliptic coords to complex coords and back 2939 """

2940 - def __call__(self,mycoord,inv=False,epoch='2000'):

2941 if inv: 2942 x,y = mycoord.real,mycoord.imag 2943 equ = ephem.Equatorial(x,y,epoch=epoch) 2944 ecl = ephem.Ecliptic(equ,epoch=epoch) 2945 return ecl.long,ecl.lat 2946 else: 2947 x,y = mycoord 2948 ecl = ephem.Ecliptic(x,y,epoch=epoch) 2949 equ = ephem.Equatorial(ecl,epoch=epoch) 2950 return float(equ.ra) + 1j*float(equ.dec)

2951

2952 -class DegCoords(NonLinearConverter):

2953 """ 2954 Convert Complex coords to degrees coordinates and back 2955 """

2956 - def __call__(self,mycoord,inv=False):

2957 if inv: 2958 x,y = mycoord.real,mycoord.imag 2959 return x/np.pi*180,y/np.pi*180 2960 else: 2961 x,y = mycoord 2962 return x/180.*np.pi + 1j*y/180.*np.pi

2963

2964 -class RadCoords(NonLinearConverter):

2965 """ 2966 Convert Complex coords to degrees coordinates and back 2967 """

2968 - def __call__(self,mycoord,inv=False):

2969 if inv: 2970 x,y = mycoord.real,mycoord.imag 2971 return x,y 2972 else: 2973 x,y = mycoord 2974 return x + 1j*y

2975

2976 #} 2977 #{ Currencies 2978 @memoized 2979 -def set_exchange_rates():

2980 """ 2981 Download currency exchange rates from the European Central Bank. 2982 """ 2983 myurl = 'http://www.ecb.europa.eu/stats/eurofxref/eurofxref-daily.xml' 2984 url = urllib.URLopener() 2985 #url = urllib.request.URLopener() # for Python 3 2986 logger.info('Downloading current exchanges rates from ecb.europa.eu') 2987 filen,msg = url.retrieve(myurl) 2988 ff = open(filen,'r') 2989 for line in ff.readlines(): 2990 if '<Cube currency=' in line: 2991 prefix,curr,interfix,rate,postfix = line.split("'") 2992 _factors[curr] = (1/float(rate),'EUR','currency','<some currency>') 2993 ff.close() 2994 #-- now also retrieve the name of the currencies: 2995 myurl = 'http://www.ecb.europa.eu/stats/exchange/eurofxref/html/index.en.html' 2996 url = urllib.URLopener() 2997 #url = urllib.request.URLopener() # for Python 3 2998 logger.info('Downloading information on currency names from ecb.europa.eu') 2999 filen,msg = url.retrieve(myurl) 3000 ff = open(filen,'r') 3001 gotcurr = False 3002 for line in ff.readlines(): 3003 if gotcurr: 3004 name = line.split('>')[1].split('<')[0] 3005 if curr in _factors: 3006 _factors[curr] = (_factors[curr][0],_factors[curr][1],_factors[curr][2],name) 3007 gotcurr = False 3008 if '<td headers="aa" id="' in line: 3009 curr = line.split('>')[1].split('<')[0] 3010 gotcurr = True 3011 ff.close()

3012

3013 3014 3015 3016 3017 #} 3018 3019 #{ Computations with units 3020 -class Unit(object):

3021 """ 3022 Class to calculate with numbers (and uncertainties) containing units. 3023 3024 You can put Units in an array and calculate with them. It is then allowed 3025 to mix Units with and without uncertainties. It is also allowed to mix 3026 Units with different units in one array, though some operations (e.g. sum) 3027 will not be possible. 3028 3029 Initalisation is done via (see L{__init__} for more options): 3030 3031 >>> a = Unit(2.,'m') 3032 >>> b = Unit(4.,'km') 3033 >>> c = Unit(3.,'cm2') 3034 3035 And you can calculate via: 3036 3037 >>> print a*b 3038 8000.0 m2 3039 >>> print a/c 3040 6666.66666667 m-1 3041 >>> print a+b 3042 4002.0 m 3043 3044 B{Example 1:} You want to calculate the equatorial velocity of the Sun: 3045 3046 >>> distance = Unit(2*np.pi,'Rsol') 3047 >>> time = Unit(22.,'d') 3048 >>> print (distance/time) 3049 2299.03495719 m1 s-1 3050 3051 or directly to km/s: 3052 3053 >>> print (distance/time).convert('km/s') 3054 2.29903495719 km/s 3055 3056 and with uncertainties: 3057 3058 >>> distance = Unit((2*np.pi,0.1),'Rsol') 3059 >>> print (distance/time).convert('km/s') 3060 2.29903495719+/-0.0365902777778 km/s 3061 3062 B{Example 2}: The surface gravity of the Sun: 3063 3064 >>> G = Unit('GG') 3065 >>> M = Unit('Msol') 3066 >>> R = Unit('Rsol') 3067 >>> logg = np.log10((G*M/R**2).convert('cgs')) 3068 >>> print logg 3069 4.43830739117 3070 3071 or 3072 3073 >>> G = Unit('GG') 3074 >>> M = Unit(1.,'Msol') 3075 >>> R = Unit(1.,'Rsol') 3076 >>> logg = np.log10((G*M/R**2).convert('cgs')) 3077 >>> print logg 3078 4.43830739117 3079 3080 or 3081 3082 >>> old = set_convention('cgs') 3083 3084 >>> G = Unit('GG') 3085 >>> M = Unit((1.,0.01),'Msol') 3086 >>> R = Unit((1.,0.01),'Rsol') 3087 >>> logg = np.log10(G*M/R**2) 3088 >>> print logg 3089 4.43830739117+/-0.00971111983789 3090 3091 >>> old = set_convention('SI') 3092 3093 B{Example 3}: The speed of light in vacuum: 3094 3095 >>> eps0 = Unit('eps0') 3096 >>> mu0 = Unit('mu0') 3097 >>> cc = Unit('cc') 3098 >>> cc_ = 1./np.sqrt(eps0*mu0) 3099 >>> print eps0 3100 8.85418781762e-12 F m-1 3101 >>> print mu0 3102 1.2566370614e-06 T m A-1 3103 >>> print cc 3104 299792458.0 m s-1 3105 >>> print cc_ 3106 299792458.004 m1 s-1 3107 >>> print cc_/cc 3108 1.00000000001 3109 3110 """

3111 - def __new__(cls,*args,**kwargs):

3112 """ 3113 Create a new Unit instance. 3114 3115 Overriding the default __new__ method makes sure that it is 3116 possible to initialize a Unit with an existing Unit instance. In 3117 that case, the original instance will simply be returned, and the 3118 __init__ functions is smart enough that it will not re-initialize 3119 the class instance. 3120 """ 3121 if not isinstance(args[0],Unit): 3122 return super(Unit, cls).__new__(cls) 3123 else: 3124 return args[0]

3125

3126 - def __init__(self,value,unit=None,**kwargs):

3127 """ 3128 Different ways to initialize a Unit. 3129 3130 Without uncertainties: 3131 3132 >>> x1 = Unit(5.,'m') 3133 >>> x4 = Unit((5.,'m')) 3134 >>> x5 = Unit(5.,'length') 3135 3136 The latter only works with the basic names ('length', 'time', 'temperature', 3137 etc..., and assures that the value is interpreted correctly in the 3138 current convention (SI,cgs...) 3139 3140 With uncertainties: 3141 3142 >>> x2 = Unit((5.,1.),'m') 3143 >>> x3 = Unit((5.,1.,'m')) 3144 3145 Initiating with an existing instance will not do anything (not even 3146 making a copy! Only a new reference to the original object is made): 3147 3148 >>> x6 = Unit(x3) 3149 3150 This is the output from printing the examples above: 3151 3152 >>> print x1 3153 5.0 m 3154 >>> print x2 3155 5.0+/-1.0 m 3156 >>> print x3 3157 5.0+/-1.0 m 3158 >>> print x4 3159 5.0 m 3160 >>> print x5 3161 5.0 m 3162 >>> print x6 3163 5.0+/-1.0 m 3164 """ 3165 if isinstance(value,Unit): 3166 return None 3167 #-- input values 3168 kwargs.setdefault('unpack',False) 3169 #-- handle different input types -- tuples etc.. 3170 if isinstance(value,tuple) and isinstance(value[-1],str): 3171 unit = value[-1] 3172 value = value[:-1] 3173 if len(value)==1: 3174 value = value[0] 3175 #-- set the input values 3176 self.value = value 3177 self.unit = unit 3178 self.kwargs = kwargs 3179 #-- make sure we can calculate with defined constants 3180 if isinstance(self.value,str) and hasattr(constants,self.value): 3181 self.unit = getattr(constants,self.value+'_units') 3182 self.value = getattr(constants,self.value) 3183 elif isinstance(self.value,str) or isinstance(self.value,tuple): 3184 try: 3185 self.value = ufloat(self.value) 3186 except: 3187 self.value = unumpy.uarray(self.value) 3188 3189 #-- values and units to work with 3190 if self.unit is not None: 3191 #-- perhaps someone says the unit is of "type" length. If so, 3192 # take the current conventions' unit of the type 3193 if not self.unit in _factors and self.unit in _conventions[constants._current_convention]: 3194 self.unit = _conventions[constants._current_convention][self.unit] 3195 #-- OK, maybe people are making it really difficult and give their 3196 # unit as 'pressure' or so... then we can still try to do 3197 # something 3198 elif not self.unit in _factors: 3199 for fac in _factors: 3200 if self.unit==_factors[fac][2]: 3201 self.unit = _factors[fac][1] 3202 break 3203 if self.unit and len(self.unit) and self.unit[0]=='[': 3204 self.value = 10**self.value 3205 self.unit = self.unit[1:-1]

3206

3207 - def convert(self,unit,compress=True):

3208 """ 3209 Convert this Unit to different units. 3210 3211 By default, the converted unit will be compressed (i.e. it will 3212 probably not consist of only basic units, but be more readable). 3213 """ 3214 if unit in _conventions: 3215 unit_ = change_convention(unit,self.unit) 3216 else: 3217 unit_ = unit 3218 new_value = convert(self.unit,unit_,self.value,**self.kwargs) 3219 new_unit = Unit(new_value,unit_) 3220 if compress: 3221 new_unit.compress() 3222 return new_unit

3223

3224 - def compress(self):

3225 """ 3226 Compress current unit to a more readable version. 3227 """ 3228 self.unit = compress(self.unit) 3229 return self

3230

3231 - def as_tuple(self):

3232 return self[0],self[1]

3233

3234 - def get_value(self):

3235 """ 3236 Returns (value,error) in case of uncertainties. 3237 """ 3238 val = unumpy.nominal_values(self.value) 3239 err = unumpy.std_devs(self.value) 3240 if not val.shape: val = float(val) 3241 if not err.shape: err = float(err) 3242 return val,err

3243

3244 - def __getitem__(self,key):

3245 """ 3246 Implements indexing, [0] returns value, [1] return unit. 3247 """ 3248 if key==0: 3249 return self.value 3250 elif key==1: 3251 return self.unit

3252

3253 - def __lt__(self,other):

3254 """ 3255 Compare SI-values of Units with eacht other. 3256 """ 3257 return self.convert('SI')[0]<other.convert('SI')[0]

3258

3259 - def __radd__(self,other):

3260 return self.__add__(other)

3261

3262 - def __add__(self,other):

3263 """ 3264 Add a Unit to a Unit. 3265 3266 You can add a non-Unit to a Unit only if the Unit has an empty unit string. 3267 """ 3268 unit1 = breakdown(self.unit)[1] 3269 if not hasattr(other,'unit'): 3270 unit2 = '' 3271 else: 3272 unit2 = breakdown(other.unit)[1] 3273 if unit1!=unit2: 3274 raise ValueError('unequal units %s and %s'%(unit1,unit2)) 3275 elif unit2=='': 3276 return self.value+other 3277 else: 3278 other_value = convert(other.unit,self.unit,other.value,unpack=False) 3279 return Unit(self.value+other_value,self.unit)

3280

3281 - def __sub__(self,other):

3282 """ 3283 Subtract a Unit from a Unit. 3284 """ 3285 return self.__add__(-1*other)

3286

3287 - def __rsub__(self,other):

3288 """ 3289 Subtract a Unit from a Unit. 3290 """ 3291 return (self*-1).__radd__(other)

3292

3293 - def __mul__(self,other):

3294 """ 3295 Multiply a Unit with something. 3296 3297 If `something' is a Unit, simply join the two Unit strings and call 3298 L{breakdown} to collect. 3299 """ 3300 if hasattr(other,'unit'): 3301 unit1 = change_convention(constants._current_convention,self.unit) 3302 unit2 = change_convention(constants._current_convention,other.unit) 3303 value1 = convert(self.unit,unit1,self.value,unpack=False) 3304 value2 = convert(other.unit,unit2,other.value,unpack=False) 3305 new_unit = ' '.join([unit1,unit2]) 3306 fac,new_unit = breakdown(new_unit) 3307 outcome = value1*value2 3308 else: 3309 outcome = other*self.value 3310 new_unit = self.unit 3311 return Unit(outcome,new_unit)

3312

3313 - def __rmul__(self,other):

3314 return self.__mul__(other)

3315

3316 - def __div__(self,other):

3317 """ 3318 Divide two units. 3319 3320 >>> x = Unit(5.,'m15') 3321 >>> y = Unit(2.5,'m6.5') 3322 >>> print(x/y) 3323 2.0 m8.5 3324 """ 3325 if hasattr(other,'unit'): 3326 unit1 = change_convention(constants._current_convention,self.unit) 3327 unit2 = change_convention(constants._current_convention,other.unit) 3328 value1 = convert(self.unit,unit1,self.value,unpack=False) 3329 value2 = convert(other.unit,unit2,other.value,unpack=False) 3330 #-- reverse signs in second units 3331 uni_b_ = '' 3332 isalpha = True 3333 prev_char_min = False 3334 for i,char in enumerate(unit2): 3335 if char=='-': 3336 prev_char_min = True 3337 continue 3338 if isalpha and not char.isalpha() and not prev_char_min: 3339 uni_b_ += '-' 3340 prev_char_min = False 3341 uni_b_ += char 3342 isalpha = char.isalpha() 3343 new_unit = ' '.join([unit1,uni_b_]) 3344 fac,new_unit = breakdown(new_unit) 3345 outcome = value1/value2 3346 else: 3347 outcome = self.value/other 3348 new_unit = self.unit 3349 return Unit(outcome,new_unit)

3350

3351 - def __rdiv__(self,other):

3352 if hasattr(other,'unit'): 3353 unit1 = change_convention(constants._current_convention,self.unit) 3354 unit2 = change_convention(constants._current_convention,other.unit) 3355 value1 = convert(self.unit,unit1,self.value,unpack=False) 3356 value2 = convert(other.unit,unit2,other.value,unpack=False) 3357 #-- reverse signs in second units 3358 uni_b_ = '' 3359 isalpha = True 3360 for i,char in enumerate(unit2): 3361 if char=='-': continue 3362 if isalpha and not char.isalpha(): 3363 uni_b_ += '-' 3364 uni_b_ += char 3365 isalpha = char.isalpha() 3366 new_unit = ' '.join([unit1,uni_b_]) 3367 fac,new_unit = breakdown(new_unit) 3368 outcome = value2/value1 3369 else: 3370 unit1 = change_convention(constants._current_convention,self.unit) 3371 value1 = convert(self.unit,unit1,self.value,unpack=False) 3372 new_unit = '' 3373 isalpha = True 3374 prev_char_min = False 3375 for i,char in enumerate(unit1): 3376 if char=='-': 3377 prev_char_min = True 3378 continue 3379 if isalpha and not char.isalpha() and not prev_char_min: 3380 new_unit += '-' 3381 prev_char_min = False 3382 new_unit += char 3383 isalpha = char.isalpha() 3384 fac,new_unit = breakdown(new_unit) 3385 outcome = other/value1 3386 return Unit(outcome,new_unit)

3387

3388 - def __pow__(self,power):

3389 """ 3390 Raise a unit to a power: 3391 3392 >>> x = Unit(5.,'m') 3393 >>> print(x**0.5) 3394 2.2360679775 m0.5 3395 >>> print(x**-1) 3396 0.2 m-1 3397 >>> print(x**3) 3398 125.0 m3 3399 3400 >>> x = Unit(9.,'m11') 3401 >>> print(x**0.5) 3402 3.0 m5.5 3403 """ 3404 unit1 = change_convention(constants._current_convention,self.unit) 3405 value1 = convert(self.unit,unit1,self.value,unpack=False) 3406 mycomps = [components(u) for u in unit1.split()] 3407 mycomps = [(u[0]**power,u[1],u[2]*power) for u in mycomps] 3408 factor = np.product([u[0] for u in mycomps]) 3409 new_unit = ' '.join(['%s%g'%(u[1],u[2]) for u in mycomps]) 3410 fac,new_unit = breakdown(new_unit) 3411 return Unit(value1**power,new_unit)

3412 #new_unit = ' '.join(power*[self._basic_unit]) 3413 #fac,new_unit = breakdown(new_unit) 3414 #return Unit(self._SI_value**power,new_unit) 3415

3416 - def sin(self): return Unit(sin(self.convert('rad').value),'')

3417 - def cos(self): return Unit(cos(self.convert('rad').value),'')

3418 - def tan(self): return Unit(tan(self.convert('rad').value),'')

3419 - def arcsin(self): return Unit(arcsin(self.value),'rad')

3420 - def arccos(self): return Unit(arccos(self.value),'rad')

3421 - def arctan(self): return Unit(arctan(self.value),'rad')

3422 - def log10(self): return Unit(log10(self.value),'')

3423 - def log(self): return Unit(log(self.value),'')

3424 - def exp(self): return Unit(exp(self.value),'')

3425 - def sqrt(self): return self**0.5

3426

3427 - def __str__(self):

3428 return '{0} {1}'.format(self.value,self.unit)

3429

3430 - def __repr__(self):

3431 return "Unit('{value}','{unit}')".format(value=repr(self.value),unit=self.unit)

3432 3433 3434 3435 #} 3436 3437 3438 3439 _fluxcalib = os.path.join(os.path.abspath(os.path.dirname(__file__)),'fluxcalib.dat') 3440 #-- basic units which the converter should know about 3441 _factors = collections.OrderedDict([ 3442 # DISTANCE 3443 ('m', ( 1e+00, 'm','length','meter')), # meter 3444 ('AA', ( 1e-10, 'm','length','angstrom')), # Angstrom 3445 ('AU', (constants.au, constants.au_units,'length','astronomical unit')), # astronomical unit 3446 ('pc', (constants.pc, constants.pc_units,'length','parsec')), # parsec 3447 ('ly', (constants.ly, constants.ly_units,'length','light year')), # light year 3448 ('bs', ( 5e-9, 'm','length','beard second')), 3449 ('Rsol', (constants.Rsol, constants.Rsol_units,'length','Solar radius')), # Solar radius 3450 ('Rearth',(constants.Rearth,constants.Rearth_units,'length','Earth radius')), # Earth radius 3451 ('Rjup', (constants.Rjup,constants.Rjup_units,'length','Jupiter radius')), # Jupiter radius 3452 ('ft', (0.3048, 'm','length','foot (international)')), # foot (international) 3453 ('USft', (1200./3937., 'm','length','foot (US)')), # foot (US) 3454 ('fathom',(1.828803657607315,'m','length','fathom')), 3455 ('rd', (5.0292, 'm','length','rod')), 3456 ('perch', (5.0292, 'm','length','pole')), 3457 ('pole', (5.0292, 'm','length','perch')), 3458 ('in', (0.0254, 'm','length','inch (international)')), # inch (international) 3459 ('mi', (1609.344, 'm','length','mile (international)')), # mile (international) 3460 ('USmi', (1609.344, 'm','length','mile (US)')), # US mile or survey mile or statute mile 3461 ('nami', (1852., 'm','length','nautical mile')), # nautical mile 3462 ('a0', (constants.a0, constants.a0_units,'length','Bohr radius')), # Bohr radius 3463 ('ell', (1.143, 'm','length','ell')), # ell 3464 ('yd', (0.9144, 'm','length','yard (international)')), # yard (international) 3465 ('potrzebie',(0.002263348517438173216473,'m','length','potrzebie')), 3466 ('smoot', (1.7018, 'm','length','smooth')), 3467 ('fur', (201.168, 'm','length','furlong')), 3468 ('ch', (20.11684, 'm','length','chain')), # (approximate) chain, based on US survey foot 3469 # VELOCITY 3470 ('mph', (0.44704, 'm s-1','velocity','miles per hour')), 3471 ('knot', (1852./3600., 'm s-1','velocity','nautical mile per hour')), 3472 # MASS 3473 ('g', ( 1e-03, 'kg','mass','gram')), # gram 3474 ('gr', (6.479891e-5, 'kg','mass','grain')), 3475 ('ton', (1e3, 'kg','mass','metric tonne')), # metric ton(ne) 3476 ('u', (1.66053892173e-27,'kg','mass','atomic mass')), # atomic mass unit (NIST 2010) 3477 ('Msol', (constants.Msol, constants.Msol_units,'mass','Solar mass')), # Solar mass 3478 ('Mearth',(constants.Mearth, constants.Mearth_units,'mass','Earth mass')), # Earth mass 3479 ('Mjup', (constants.Mjup, constants.Mjup_units,'mass','Jupiter mass')), # Jupiter mass 3480 ('Mlun', (constants.Mlun, constants.Mlun_units,'mass','Lunar mass')), # Lunar mass 3481 ('lb', (0.45359237, 'kg','mass','pound')), # pound 3482 ('st', (6.35029318, 'kg','mass','stone')), # stone 3483 ('ounce', (0.0283495231, 'kg','mass','ounce')), # ounce 3484 ('mol', (1./constants.NA,'mol','mass','molar mass')), # not really a mass... 3485 ('firkin',(40.8233, 'kg','mass','firkin')), 3486 ('carat', (0.0002, 'kg','mass','carat')), # carat, metric 3487 # TIME 3488 ('s', ( 1e+00, 's','time','second')), # second 3489 ('min', ( 60., 's','time','minute')), # minute 3490 ('h', (3600., 's','time','hour')), # hour 3491 ('d', (86400., 's','time','day')), # day 3492 ('wk', (7*24*3600., 's','time','week')), # week 3493 ('mo', (30*24*3600., 's','time','month')), # month 3494 ('sidereal', (1.0027379093,'','time','sidereal day')), # sidereal 3495 ('yr', (31557600.0,'s','time','Julian year')), # year (1 Julian century = 36525 d, see NIST appendix B) 3496 ('cr', (36525*86400,'s','time','Julian century')), # century 3497 ('hz', (2*np.pi, 'rad s-1','frequency','Hertz')),# Hertz (periodic phenomena) 3498 # ('rhz', (1e+00, 'rad s-1','angular frequency','RadHertz')),# Rad Hertz (periodic phenomena) 3499 ('Bq', (1e+00, 's-1','time','Becquerel')), # Becquerel (stochastic or non-recurrent phenomena) 3500 ('Ci', (3.7e10, 's-1','time','Curie')), 3501 ('R', (2.58e-4, 'C kg-1','roentgen','Roentgen')), 3502 ('JD', (1e+00, 'JD','time','Julian day')), # Julian Day 3503 ('CD', (JulianDay, 'JD','time','calender day')), # Calender Day 3504 ('MJD', (ModJulianDay, 'JD','time','modified Julian day')), # Modified Julian Day 3505 ('j', (1/60., 's','time','jiffy')), # jiffy 3506 ('fortnight',(1209600., 's','second','fortnight')), 3507 # ANGLES 3508 ('rad', (1e+00, 'rad','angle','radian')), # radian 3509 ('cy', (2*np.pi, 'rad','angle','cycle')), # cycle 3510 ('deg', (np.pi/180., 'rad','angle','degree')), # degree 3511 ('am', (np.pi/180./60., 'rad','angle','arcminute')), # arcminute 3512 ('as', (np.pi/180./3600., 'rad','angle','arcsecond')), # arcsecond 3513 ('sr', (1, 'sr','angle','sterradian')), # sterradian #1/39.4784176045 3514 ('rpm', (0.104719755, 'rad s-1','angle','revolutions per minute')),# revolutions per minute 3515 # COORDINATES 3516 ('complex_coord',(1e+00+0*1j, 'complex_coord','coordinate','<own unit>')), # own unit 3517 ('equatorial', (EquCoords, 'complex_coord','coordinate','equatorial')), # Equatorial coordinates 3518 ('galactic', (GalCoords, 'complex_coord','coordinate','galactic')), # Galactic coordinates 3519 ('ecliptic', (EclCoords, 'complex_coord','coordinate','ecliptic')), # Ecliptic coordinates 3520 ('deg_coord', (DegCoords, 'complex_coord','coordinate','degrees')), # Coordinates in degrees 3521 ('rad_coord', (RadCoords, 'complex_coord','coordinate','radians')), # Coordinates in radians 3522 # FORCE 3523 ('N', (1e+00, 'kg m s-2','force','Newton')), # newton 3524 ('dyn', (1e-05, 'kg m s-2','force','dyne')), # dyne 3525 # TEMPERATURE 3526 ('K', (1e+00, 'K','temperature','Kelvin')), # Kelvin 3527 ('Far', (Fahrenheit, 'K','temperature','Fahrenheit')), # Fahrenheit 3528 ('Cel', (Celcius, 'K','temperature','Celcius')), # Celcius 3529 ('Tsol', (constants.Tsol,constants.Tsol_units,'temperature','Solar temperature')), # solar temperature 3530 # ENERGY & POWER 3531 ('J', ( 1e+00, 'kg m2 s-2','energy','Joule')), # Joule 3532 ('W', ( 1e+00, 'kg m2 s-3','power','Watt')), # Watt 3533 ('erg', ( 1e-07, 'kg m2 s-2','energy','ergon')), # ergon 3534 ('foe', ( 1e44, 'kg m2 s-2','energy','(ten to the) fifty one ergs')), 3535 ('eV', (1.60217646e-19,'kg m2 s-2','energy','electron volt')), # electron volt 3536 ('cal', (4.1868, 'kg m2 s-2','energy','small calorie (international table)')),# calorie (International table) 3537 ('Cal', (4186.8, 'kg m2 s-2','energy','large calorie (international table)')),# calorie (International table) 3538 ('Lsol', (constants.Lsol, constants.Lsol_units,'energy/power','Solar luminosity')), # solar luminosity 3539 ('hp', (745.699872, 'kg m2 s-3','power','Horsepower')), # horsepower 3540 ('dp', (250., 'kg m2 s-3','power','Donkeypower')), 3541 ('Gy', (1., 'm2 s-2','absorbed dose','gray')), # absobred dose, specific energy (imparted) or kerma 3542 ('Sv', (1., 'm2 s-2','dose equivalent','sievert')), 3543 ('kat', (1., 'mol s-1','catalytic activity','katal')), 3544 ('rem', (1e-2, 'm2 s-2','dose equivalent','rem')), 3545 ('dB', (DecibelSPL, 'kg s-3','sound intensity','Decibel')), # W/m2 3546 # VELOCITY 3547 ('cc', (constants.cc, constants.cc_units,'m/s','Speed of light')), 3548 # ACCELERATION 3549 ('Gal', (1., 'm s-2','acceleration','Gal')), 3550 # PRESSURE 3551 ('Pa', ( 1e+00, 'kg m-1 s-2','pressure','Pascal')), # Pascal 3552 ('bar', ( 1e+05, 'kg m-1 s-2','pressure','baros')), # baros 3553 ('ba', ( 0.1, 'kg m-1 s-2','pressure','barye')), # Barye 3554 ('at', ( 98066.5, 'kg m-1 s-2','pressure','atmosphere (technical)')), # atmosphere (technical) 3555 ('atm', ( 101325, 'kg m-1 s-2','pressure','atmosphere (standard)')), # atmosphere (standared) 3556 ('torr', ( 133.322, 'kg m-1 s-2','pressure','Torricelli')), # Torricelli 3557 ('psi', ( 6894., 'kg m-1 s-2','pressure','pound per square inch')), # pound per square inch 3558 ('mmHg', (133.322, 'kg m-1 s-2','pressure','millimeter of mercury')), 3559 ('P', ( 0.1, 'kg m-1 s-1','dynamic viscosity','poise')), # CGS 3560 ('St', ( 1e-4, 'm2 s-1', 'kynamatic viscosity','stokes')), # CGS 3561 # MAGNETISM & Electricity 3562 ('A', ( 1., 'A', 'electric current','Ampere')), 3563 ('Gi', ( 0.7957747, 'A', 'electric current','Ampere')), # (approximate) 3564 ('Bi', ( 10., 'A', 'electric current','biot')), 3565 ('C', ( 1., 'A s', 'electric charge','Coulomb')), 3566 ('T', ( 1., 'kg s-2 A-1','magnetic field strength','Tesla')), 3567 ('G', ( 1e-4, 'kg s-2 A-1','magnetic field strength','Gauss')), 3568 ('Oe', (1000/(4*np.pi),'A m-1','magnetizing field','Oersted')), 3569 ('V', ( 1., 'kg m2 s-3 A-1','electric potential difference','Volt')), 3570 ('F', ( 1., 'kg-1 m-2 s4 A2','electric capacitance','Farad')), 3571 ('O', ( 1., 'kg m2 s-3 A-2','electric resistance','Ohm')), 3572 ('S', ( 1., 'kg-1 m-2 s3 A2','electric conductance','Siemens')), 3573 ('Wb', ( 1., 'kg m2 s-2 A-1','magnetic flux','Weber')), 3574 ('Mx', ( 1e-8, 'kg m2 s-2 A-1','magnetic flux','Maxwell')), 3575 ('H', ( 1., 'kg m2 s-2 A-2','inductance','Henry')), 3576 ('lm', ( 1., 'cd sr','luminous flux','lumen')), 3577 ('lx', ( 1., 'cd sr m-2','illuminance','lux')), 3578 ('sb', ( 1e4, 'cd m-2','illuminance','stilb')), 3579 ('ph', ( 1e4, 'cd sr m-2','illuminance','phot')), 3580 # AREA & VOLUME 3581 ('ac', (4046.873, 'm2','area','acre (international)')), # acre (international) 3582 ('a', (100., 'm2','area','are')), # are 3583 ('b', (1e-28, 'm2','area','barn')), # barn 3584 ('l', ( 1e-3, 'm3','volume','liter')), 3585 ('gal', (4.54609e-3, 'm3','volume','gallon (Canadian and UK imperial)')), 3586 ('USgal', (3.785412e-3, 'm3','volume','gallon (US)')), 3587 ('gi', (1.420643e-4, 'm3','volume','gill (Canadian and UK imperial)')), 3588 ('USgi', (1.182941e-4, 'm3','volume','gill (Canadian and UK imperial)')), 3589 ('bbl', (0.1589873, 'm3','volume','barrel')), # (approximate) US barrel) 3590 ('bu', (0.03523907, 'm3','volume','bushel')), # (approximate) US bushel 3591 ('ngogn', (1.1594560721765171e-05,'m3','volume','1000 cubic potrzebies')), 3592 # FLUX 3593 # -- absolute magnitudes 3594 ('Jy', (1e-26/(2*np.pi),'kg s-2 rad-1','flux density','Jansky')), # W/m2/Hz 3595 ('fu', (1e-26/(2*np.pi),'kg s-2 rad-1','flux density','flux unit')), # W/m2/Hz 3596 ('vegamag', (VegaMag, 'kg m-1 s-3','flux','Vega magnitude')), # W/m2/m 3597 ('mag', (VegaMag, 'kg m-1 s-3','flux','magnitude')), # W/m2/m 3598 ('STmag', (STMag, 'kg m-1 s-3','flux','ST magnitude')), # W/m2/m 3599 ('ABmag', (ABMag, 'kg s-2 cy-1','flux','AB magnitude')), # W/m2/Hz 3600 # -- magnitude differences (colors) 3601 ('mag_color',(Color, 'flux_ratio','flux','color')), 3602 ('flux_ratio',(1+00, 'flux_ratio','flux','flux ratio')), 3603 # -- magnitude amplitudes 3604 ('ampl', (1e+00, 'ampl','flux','fractional amplitude')), 3605 ('Amag', (AmplMag, 'ampl','flux','amplitude in magnitude')), 3606 ('pph', (1e-02, 'ampl','flux','amplitude in parts per hundred')), # amplitude 3607 ('ppt', (1e-03, 'ampl','flux','amplitude in parts per thousand')), # amplitude 3608 ('ppm', (1e-06, 'ampl','flux','amplitude in parts per million')), # amplitude 3609 # -- currency 3610 ('EUR', (1e+00, 'EUR','currency','EURO')) 3611 ]) 3612 #-- set of conventions: 3613 _conventions = {'SI': dict(mass='kg',length='m', time='s',temperature='K', 3614 electric_current='A',lum_intens='cd',amount='mol',currency='EUR'), # International standard 3615 'cgs':dict(mass='g', length='cm',time='s',temperature='K', 3616 electric_current='A',lum_intens='cd',amount='mol',currency='EUR'), # Centi-gramme-second 3617 'sol':dict(mass='Msol',length='Rsol',time='s',temperature='Tsol', 3618 electric_current='A',lum_intens='cd',amount='mol',currency='EUR'), # solar 3619 'imperial':dict(mass='lb',length='yd',time='s',temperature='K', 3620 electric_current='A',lum_intens='cd',amount='mol',currency='GBP'), # Imperial (UK/US) system 3621 } 3622 3623 #-- some names of combinations of units 3624 _names = {'m1':'length', 3625 's1':'time', 3626 'kg1':'mass', 3627 'rad':'angle', 3628 'deg':'angle', 3629 'K':'temperature', 3630 'm1 s-2':'acceleration', 3631 'kg1 s-2':'surface tension', 3632 'cy-1 kg1 s-2':'flux density', # W/m2/Hz 3633 'kg1 m-1 s-3':'flux density', # W/m3 3634 'kg1 s-3':'flux', # W/m2 3635 'cy1 s-1':'frequency', 3636 'cy-1 kg1 rad-2 s-2':'specific intensity', # W/m2/Hz/sr 3637 'kg1 m-1 rad-2 s-3':'specific intensity', # W/m3/sr 3638 'kg1 rad-2 s-3':'total intensity', # W/m2/sr 3639 'kg1 m2 s-3':'luminosity', # W or power 3640 'm1 s-1':'velocity', 3641 'kg1 m-1 s-2':'pressure', 3642 'm2':'area', 3643 'm3':'volume', 3644 'kg1 m2 s-2':'energy', 3645 'kg1 m1 s-2':'force', 3646 'A':'electric current', 3647 'A s':'electric charge', 3648 'kg s-2 A-1':'magnetic field strength', 3649 'A m-1':'magnetizing field', 3650 'kg m2 s-3 A-1':'electric potential difference', 3651 'kg-1 m-2 s4 A2':'electric capacitance', 3652 'kg m2 s-3 A-2':'electric resistance', 3653 'kg-1 m-2 s3 A2':'electric conductance', 3654 'kg m2 s-2 A-1':'magnetic flux', 3655 'kg m2 s-2 A-2':'inductance', 3656 'cd':'luminous flux', 3657 'm-2 cd':'illuminance', 3658 } 3659 3660 #-- scaling factors for prefixes 3661 _scalings ={'y': 1e-24, # yocto 3662 'z': 1e-21, # zepto 3663 'a': 1e-18, # atto 3664 'f': 1e-15, # femto 3665 'p': 1e-12, # pico 3666 'n': 1e-09, # nano 3667 'mu': 1e-06, # micro 3668 'u': 1e-06, # micro 3669 'm': 1e-03, # milli 3670 'c': 1e-02, # centi 3671 'd': 1e-01, # deci 3672 'da': 1e+01, # deka 3673 'h': 1e+02, # hecto 3674 'k': 1e+03, # kilo 3675 'M': 1e+06, # mega 3676 'G': 1e+09, # giga 3677 'T': 1e+12, # tera 3678 'P': 1e+15, # peta 3679 'E': 1e+18, # exa 3680 'Z': 1e+21, # zetta 3681 'Y': 1e+24 # yotta 3682 } 3683 3684 #-- some common aliases 3685 _aliases = [('micron','mum'),('au','AU'),('lbs','lb'), 3686 ('micro','mu'),('milli','m'),('kilo','k'),('mega','M'),('giga','G'), 3687 ('nano','n'),('dyne','dyn'), 3688 ('watt','W'),('Watt','W'), 3689 ('Hz','hz'), ('liter','l'),('litre','l'), 3690 ('Tesla','T'),('Ampere','A'),('Coulomb','C'), 3691 ('joule','J'),('Joule','J'), 3692 ('jansky','Jy'),('Jansky','Jy'),('jy','Jy'), 3693 ('arcsec','as'),('arcmin','am'), 3694 ('cycles','cy'),('cycle','cy'),('cyc','cy'), 3695 ('angstrom','AA'),('Angstrom','AA'), 3696 ('inch','in'),('stone','st'), 3697 ('^',''),('**',''), 3698 ('celcius','C'),('fahrenheit','F'),('hr','h'), 3699 ('Vegamag','vegamag'),('mile','mi'), 3700 ('oz','ounce'),('sun','sol'),('_sun','sol'),('_sol','sol'), 3701 ('solMass','Msol'),('solLum','Lsol'), 3702 ('pk','hp'),('mph','mi/h'),('f.u.','fu') 3703 ] 3704 3705 #-- Change-of-base function definitions 3706 _switch = {'s1_to_': distance2velocity, # switch from wavelength to velocity 3707 's-1_to_': velocity2distance, # switch from wavelength to velocity 3708 'm1rad-1_to_': distance2spatialfreq, # for interferometry 3709 'm-1rad1_to_': spatialfreq2distance, # for interferometry 3710 'm1rad1s1_to_': fnu2flambda, # dunno! 3711 'm-1rad-1s-1_to_': flambda2fnu,# dunno! 3712 'm1rad-1s1_to_': fnu2flambda, # switch from Fnu to Flam 3713 'm-1rad1s-1_to_':flambda2fnu, # switch from Flam to Fnu 3714 'rad-1s1_to_': fnu2nufnu, # switch from Fnu to nuFnu 3715 'rad1s-1_to_': nufnu2fnu, # switch from nuFnu to Fnu 3716 'm1_to_': lamflam2flam, # switch from lamFlam to Flam 3717 'm-1_to_': flam2lamflam, # switch from Flam to lamFlam 3718 #'rad2_to_': per_sr, 3719 #'rad-2_to_': times_sr, 3720 'sr1_to_': per_sr, 3721 'sr-1_to_': times_sr, 3722 'rad1_to_': do_nothing,#per_cy, 3723 'rad-1_to_': do_nothing,#times_cy, 3724 'rad-2sr1_to_': do_nothing, 3725 'rad2sr-1_to_': do_nothing, 3726 'rad-1sr1_to_': do_sqrt, 3727 'rad1sr-1_to_': do_quad, 3728 'rad-1s2_to_': period2freq, # same for both, since just inverse 3729 'rad1s-2_to_': period2freq, 3730 #'cy1_to_': do_nothing, 3731 #'cy-1_to_': do_nothing, 3732 'cy2_to_': do_nothing, 3733 'cy-2_to_': do_nothing, 3734 } 3735 3736 3737 if __name__=="__main__": 3738 3739 if not sys.argv[1:]: 3740 import doctest 3741 doctest.testmod() 3742 quit() 3743 from optparse import OptionParser, Option, OptionGroup 3744 import datetime 3745 import copy 3746 3747 logger = loggers.get_basic_logger()

3748 3749 #-- make sure we can parse strings as Python code 3750 - def check_pythoncode(option, opt, value):

3751 try: 3752 return eval(value) 3753 except ValueError: 3754 raise OptionValueError( 3755 "option %s: invalid python code: %r" % (opt, value))

3756

3757 #-- generate a custom help log 3758 - class MyOption (Option):

3759 TYPES = Option.TYPES + ("pythoncode",) 3760 TYPE_CHECKER = copy.copy(Option.TYPE_CHECKER) 3761 TYPE_CHECKER["pythoncode"] = check_pythoncode

3762 usage = "List of available units:\n" + get_help() + "\nUsage: %prog --from=<unit> --to=<unit> [options] value [error]" 3763 3764 #-- define all the input parameters 3765 parser = OptionParser(option_class=MyOption,usage=usage) 3766 parser.add_option('--from',dest='_from',type='str', 3767 help="units to convert from",default='SI') 3768 parser.add_option('--to',dest='_to',type='str', 3769 help="units to convert to",default='SI') 3770 3771 group = OptionGroup(parser,'Extra quantities when changing base units (e.g. Flambda to Fnu)') 3772 group.add_option('--wave','-w',dest='wave',type='str', 3773 help="wavelength with units (e.g. used to convert Flambda to Fnu)",default=None) 3774 group.add_option('--freq','-f',dest='freq',type='str', 3775 help="frequency (e.g. used to convert Flambda to Fnu)",default=None) 3776 group.add_option('--photband','-p',dest='photband',type='str', 3777 help="photometric passband",default=None) 3778 parser.add_option_group(group) 3779 3780 #-- prepare inputs for functions 3781 (options, args) = parser.parse_args() 3782 options = vars(options) 3783 _from = options.pop('_from') 3784 _to = options.pop('_to') 3785 3786 #-- in case of normal floats or floats with errors 3787 if not any([',' in i for i in args]): 3788 args = tuple([float(i) for i in args]) 3789 #-- in case of tuples (like coordinates) 3790 else: 3791 args = tuple((eval(args[0]),)) 3792 #-- remove None types 3793 for option in copy.copy(options): 3794 if options[option] is None: 3795 options.pop(option) 3796 #-- set type correctly 3797 for option in options: 3798 if isinstance(options[option],str) and ',' in options[option]: 3799 entry = options[option].split(',') 3800 options[option] = (float(entry[0]),entry[1]) 3801 3802 #-- check if currencies are asked. If so, download the latest exchange rates 3803 # we cannot exactly know if something is a currency before we download all 3804 # the names, and we don't want to download the names if no currencies are 3805 # given. Therefore, we make an educated guess here: 3806 from_is_probably_currency = (hasattr(_from,'isupper') and _from.isupper() and len(_from)==3) 3807 to_is_probably_currency = (hasattr(_to,'isupper') and _to.isupper() and len(_to)==3) 3808 if from_is_probably_currency or to_is_probably_currency: 3809 set_exchange_rates() 3810 3811 #-- do the conversion 3812 output = convert(_from,_to,*args,**options) 3813 3814 #-- and nicely print to the screen 3815 if _to=='SI': 3816 fac,_to = breakdown(_from) 3817 if _from=='SI': 3818 fac,_from = breakdown(_to) 3819 if isinstance(output,tuple) and len(output)==2 and len(args)==2: 3820 print(("%g +/- %g %s = %g +/- %g %s"%(args[0],args[1],_from,output[0],output[1],_to))) 3821 elif isinstance(output,tuple) and len(output)==2: 3822 print(("%s %s = %s,%s %s"%(args[0],_from,output[0],output[1],_to))) 3823 elif _to.lower()=='cd': 3824 year,month,day = output 3825 year,month = int(year),int(month) 3826 day,fraction = int(day),day-int(day) 3827 hour = fraction*24 3828 hour,fraction = int(hour),hour-int(hour) 3829 minute = fraction*60 3830 minute,fraction = int(minute),minute-int(minute) 3831 second = int(fraction*60) 3832 dt = datetime.datetime(year,month,day,hour,minute,second) 3833 print(("%.10g %s = %s %s (YYYY-MM-DD HH:MM:SS)"%(args[0],_from,dt,_to))) 3834 else: 3835 print(("%.10g %s = %.10g %s"%(args[0],_from,output,_to))) 3836

Source Code for Module ivs.units.conversions