Package gavo :: Package stc :: Module units
[frames] | no frames]

Source Code for Module gavo.stc.units

  1  """ 
  2  Definition and conversion of units in STC 
  3   
  4  For every physical quantity we deal with, there is a standard unit defined: 
  5   
  6          - angles: deg  (we way want to use rad here) 
  7          - distances: m 
  8          - time: s 
  9          - frequencies: Hz 
 10          - wavelength/energy: m 
 11   
 12  We keep dictionaries of conversion factors to those units.  Since turning 
 13  things around, these factors give "how many bases are in the unit", e.g. 
 14  km -> 1000. 
 15   
 16  The main interface are functions returning converter functions.  Pass 
 17  a value in fromUnit to them and receive a value in toUnit.  Simple factors 
 18  unfortunately don't cut it here since conversion from wavelength to 
 19  frequency needs division of the value. 
 20  """ 
 21   
 22  #c Copyright 2008-2019, the GAVO project 
 23  #c 
 24  #c This program is free software, covered by the GNU GPL.  See the 
 25  #c COPYING file in the source distribution. 
 26   
 27   
 28  import math 
 29   
 30  from gavo.utils import memoized 
 31  from gavo.stc import common 
 32   
 33   
 34  toRad=math.pi/180. 
 35  oneAU = 1.49597870691e11   # IAU 
 36  onePc = oneAU/2/math.tan(0.5/3600*toRad) 
 37  lightspeed = 2.99792458e8  # SI 
 38  planckConstant = 4.13566733e-15  # CODATA 2008, in eV s 
 39  julianYear = 365.25*24*3600 
40 41 -def makeConverterMaker(label, conversions):
42 """returns a conversion function that converts between any of the units 43 mentioned in the dict conversions. 44 """ 45 def getConverter(fromUnit, toUnit, reverse=False): 46 if fromUnit not in conversions or toUnit not in conversions: 47 raise common.STCUnitError("One of '%s' or '%s' is no valid %s unit"%( 48 fromUnit, toUnit, label)) 49 fact = conversions[fromUnit]/conversions[toUnit] 50 if reverse: 51 fact = 1/fact 52 def convert(val): 53 return fact*val
54 return convert 55 return getConverter 56 57 58 # Factors like "one kilometer is 1e3 meters" 59 distFactors = { 60 "m": 1., 61 "km": 1e3, 62 "mm": 1e-3, 63 "AU": oneAU, 64 "pc": onePc, 65 "kpc": (1e3*onePc), 66 "Mpc": (1e6*onePc), 67 "lyr": (lightspeed*julianYear), 68 } 69 getDistConv = makeConverterMaker("distance", distFactors) 70 71 72 angleFactors = { 73 "deg": 1., 74 "rad": 1/toRad, 75 "h": 360./24., 76 "arcmin": 1/60., 77 "arcsec": 1/3600., 78 } 79 getAngleConv = makeConverterMaker("angle", angleFactors) 80 81 82 timeFactors = { 83 "s": 1., 84 "h": 3600., 85 "d": (3600.*24), 86 "a": julianYear, 87 "yr": julianYear, 88 "cy": (julianYear*100), 89 } 90 getTimeConv = makeConverterMaker("time", timeFactors) 91 92 93 # spectral units have the additional intricacy that a factor is not 94 # enough when wavelength needs to be converted to a frequency. 95 freqFactors = { 96 "Hz": 1., 97 "kHz": 1e3, 98 "MHz": 1e6, 99 "GHz": 1e9, 100 "eV": 1/planckConstant, 101 "keV": 1/planckConstant*1e3, 102 "MeV": 1/planckConstant*1e6, 103 "GeV": 1/planckConstant*1e9, 104 "TeV": 1/planckConstant*1e12, 105 } 106 getFreqConv = makeConverterMaker("frequency", freqFactors) 107 108 wlFactors = { 109 "m": 1., 110 "mm": 1e-3, 111 "um": 1e-6, 112 "nm": 1e-9, 113 "Angstrom": 1e-10, 114 } 115 getWlConv = makeConverterMaker("wavelength", wlFactors)
116 117 -def getSpectralConv(fromUnit, toUnit, reverse=False):
118 if fromUnit in wlFactors: 119 if toUnit in wlFactors: 120 conv = getWlConv(fromUnit, toUnit, reverse) 121 else: # toUnit is freq 122 fromFunc = getWlConv(fromUnit, "m") 123 toFunc = getFreqConv("Hz", toUnit, reverse) 124 def conv(val): 125 return toFunc(lightspeed/fromFunc(val))
126 else: # fromUnit is freq 127 if toUnit in freqFactors: 128 conv = getFreqConv(fromUnit, toUnit, reverse) 129 else: # toUnit is wl 130 fromFunc = getFreqConv(fromUnit, "Hz", reverse) 131 toFunc = getWlConv("m", toUnit) 132 def conv(val): 133 return toFunc(lightspeed/fromFunc(val)) 134 return conv 135 136 137 distUnits = set(distFactors) 138 angleUnits = set(angleFactors) 139 timeUnits = set(timeFactors) 140 spectralUnits = set(wlFactors) | set(freqFactors) 141 142 systems = [(distUnits, getDistConv), (angleUnits, getAngleConv), 143 (timeUnits, getTimeConv), (spectralUnits, getSpectralConv)]
144 145 146 @memoized 147 -def getBasicConverter(fromUnit, toUnit, reverse=False):
148 """returns a function converting fromUnit values to toUnit values. 149 """ 150 for units, factory in systems: 151 if fromUnit in units and toUnit in units: 152 return factory(fromUnit, toUnit, reverse) 153 raise common.STCUnitError("No known conversion from '%s' to '%s'"%( 154 fromUnit, toUnit))
155 156 157 # the maximal parallax distance as parallax. This is used in the parallax 158 # converters to avoid divisions by zero. 159 maxDistance = 1e7
160 161 162 @memoized 163 -def getParallaxConverter(fromUnit, toUnit, reverse=False):
164 """returns a function converting distances to/from parallaxes. 165 """ 166 if fromUnit not in angleUnits: 167 fromUnit, toUnit, reverse = toUnit, fromUnit, not reverse 168 if fromUnit not in angleUnits: 169 raise common.STCUnitError("No spatial conversion between %s and %s"%( 170 fromUnit, toUnit)) 171 # first convert angular unit to arcsec, then invert, yielding pc, 172 # and convert that to distance unit 173 angularConv = getBasicConverter(fromUnit, "arcsec", reverse) 174 distanceConv = getBasicConverter("pc", toUnit, reverse) 175 if reverse: 176 def conv(val): #noflake: local function 177 res = distanceConv(val) 178 if res>maxDistance: 179 return 0. 180 else: 181 return angularConv(1./res)
182 else: 183 def conv(val): #noflake: local function 184 res = angularConv(val) 185 if res<1/maxDistance: 186 return distanceConv(maxDistance) 187 else: 188 return distanceConv(1./res) 189 return conv 190
191 192 @memoized 193 -def getRedshiftConverter(spaceUnit, timeUnit, toSpace, toTime, 194 reverse=False):
195 """returns a function converting redshifts in spaceUnit/timeUnit to 196 toSpace/toTime. 197 198 This will actually work for any unit of the form unit1/unit2 as long 199 as unit2 is purely multiplicative. 200 """ 201 spaceFun = getBasicConverter(spaceUnit, toSpace, reverse) 202 timeFun = getBasicConverter(timeUnit, toTime, not reverse) 203 def convert(val): 204 return spaceFun(timeFun(val))
205 return convert 206
207 208 -def _expandUnits(fromUnits, toUnits):
209 """makes sure fromUnits and toUnits have the same length. 210 211 This is a helper for vector converters. 212 """ 213 if isinstance(toUnits, basestring): 214 toUnits = (toUnits,)*len(fromUnits) 215 if len(fromUnits)!=len(toUnits): 216 raise common.STCUnitError( 217 "Values in %s cannot be converted to values in %s"%(fromUnits, toUnits)) 218 return toUnits
219
220 @memoized 221 -def getVectorConverter(fromUnits, toUnits, reverse=False):
222 """returns a function converting from fromUnits to toUnits. 223 224 fromUnits is a tuple, toUnits is a tuple of which only the first item 225 is interpreted. This first item must be a tuple or a single string; in the 226 latter case, all components are supposed to be of that unit. 227 228 ToUnits may be shorter than fromUnits. In this case, the additional 229 fromUnits are ignored. This is mainly for cases in which geometries go 230 with SPHER3 positions. 231 232 The resulting functions accepts sequences of len(toUnits) and returns 233 tuples of the same length. 234 235 As a special service for Geometries, these spatial converters have 236 additional attributes fromUnit and toUnit giving what transformation 237 they implement. 238 """ 239 if not fromUnits: # no base unit given, we give up 240 def convert(val): #noflake: local function 241 return val
242 return convert 243 244 toUnits = _expandUnits(fromUnits, toUnits) 245 convs = [] 246 convs.append(getBasicConverter(fromUnits[0], toUnits[0], reverse)) 247 if len(toUnits)>1: 248 convs.append(getBasicConverter(fromUnits[1], toUnits[1], reverse)) 249 if len(toUnits)>2: 250 try: 251 convs.append(getBasicConverter(fromUnits[2], toUnits[2], reverse)) 252 except common.STCUnitError: # try parallax for the last unit only 253 convs.append(getParallaxConverter(fromUnits[2], toUnits[2], reverse)) 254 255 def convert(val): #noflake: local function 256 if not hasattr(val, "__iter__"): 257 # someone sneaked in a scalar. Sigh 258 val = (val,) 259 return tuple(f(c) for f, c in zip(convs, val)) 260 if reverse: 261 convert.fromUnit, convert.toUnit = toUnits, fromUnits 262 else: 263 convert.fromUnit, convert.toUnit = fromUnits, toUnits 264 return convert 265
266 267 @memoized 268 -def getVelocityConverter(fromSpaceUnits, fromTimeUnits, toSpace, toTime, 269 reverse=False):
270 """returns a function converting from fromSpaceUnits/fromTimeUnits to 271 toSpace/toTime. 272 273 fromXUnits is a tuple, toX may either be a tuple of length fromXUnits or a a 274 single string like in getVectorUnits. 275 276 The resulting functions accepts sequences of proper length and returns 277 tuples. 278 """ 279 toSpace = _expandUnits(fromSpaceUnits, toSpace) 280 toTime = _expandUnits(fromTimeUnits, toTime) 281 convs = tuple(getRedshiftConverter(fs, ft, ts, tt, reverse) 282 for fs, ft, ts, tt in zip( 283 fromSpaceUnits, fromTimeUnits, toSpace, toTime)) 284 def convert(val): 285 return tuple(f(c) for f, c in zip(convs, val))
286 return convert 287
288 289 @memoized 290 -def getUnitConverter(fromCoo, toCoo):
291 """returns a pair of unit info and a conversion function to take fromCoo 292 to the units of toCoo. 293 294 toCoo may be None, in which case the unit of fromCoo is returned together 295 with an identity function. If the units already match, (None, None) is 296 returned. 297 298 The unit info is a dictionary suitable for change(). 299 """ 300 if toCoo is None or toCoo.getUnitArgs() is None: 301 return fromCoo.getUnitArgs(), None 302 if fromCoo.getUnitArgs() is None: 303 return toCoo.getUnitArgs(), None 304 if fromCoo.getUnitArgs()==toCoo.getUnitArgs(): 305 return None, None 306 return toCoo.getUnitArgs(), toCoo.getUnitConverter( 307 fromCoo.getUnitArgs())
308
309 310 -def iterUnitAdapted(baseSTC, sysSTC, attName, dependentName):
311 """iterates over all keys that need to be changed to adapt units in baseSTC's 312 attName facet to conform to what sysSTC gives. 313 314 If something in baseSTC is not specified in sysSTC, it is ignored here 315 (i.e., it will remain unchanged if the result is used in a change). 316 317 Since units are only on coordinates, and areas inherit these units, 318 they are transformed as well, and their name is given by dependentName. 319 See also conform.conformUnits. 320 """ 321 coo = getattr(baseSTC, attName) 322 if coo is None: 323 return 324 overrides, conv = getUnitConverter(coo, getattr(sysSTC, attName)) 325 if conv is None: # units are already ok 326 return 327 overrides.update(coo.iterTransformed(conv)) 328 yield attName, coo.change(**overrides) 329 areas = getattr(baseSTC, dependentName) 330 if areas: 331 transformed = [] 332 for a in areas: 333 transformed.append(a.adaptValuesWith(conv)) 334 yield dependentName, tuple(transformed)
335