Source Code for Module gavo.stc.geo

 1  """ 
 2  Coordinate systems for positions on earth. 
 3  """ 
 4   
 5  #c Copyright 2008-2019, the GAVO project 
 6  #c 
 7  #c This program is free software, covered by the GNU GPL.  See the 
 8  #c COPYING file in the source distribution. 
 9   
10   
11  import math 
12   
13  from gavo.utils import DEG 
14   
15   
16 -class WGS84(object): 
17          """the WGS84 reference system. 
18          """ 
19          a = 6378137. 
20          f = 1/298.257223563 
21          GM = 3.986005e14    # m3s-1 
22          J2 = 0.00108263 
23          omega = 7.292115e-5 # rad s-1 
24   
25   
26 -def _getC_S(phi, refSys): 
27          """returns the values of the auxiliary functions C and S. 
28   
29          phi must be in rad. 
30   
31          See Astron. Almanac, Appendix K. 
32          """ 
33          B = (1-refSys.f)**2 
34          C = 1/math.sqrt( 
35                  math.cos(phi)**2 
36                  +B*math.sin(phi)**2) 
37          S = C*B 
38          return C, S 
39   
40   
41 -def geocToGeod(long, phip, rho=1, refSys=WGS84): 
42          """returns geodetic coordinates long, phi for geocentric coordinates. 
43   
44          refSys defaults is the reference system the geodetic coordinates are 
45          expressed in. 
46   
47          This will not work at the poles -- patches welcome. 
48   
49          See Astron. Almanac, Appendix K; we go for the iterative solution 
50          discussed there. 
51          """ 
52          long, phip = long*DEG, phip*DEG 
53          x = refSys.a*rho*math.cos(phip)*math.cos(long) 
54          y = refSys.a*rho*math.cos(phip)*math.sin(long) 
55          z = refSys.a*rho*math.sin(phip) 
56   
57          e2 = 2*refSys.f-refSys.f**2 
58          r = math.sqrt(x**2+y**2) 
59          phi = math.atan2(z, r) 
60   
61          while True: 
62                  phi1 = phi 
63                  C = 1/math.sqrt((1-e2*math.sin(phi1)**2)) 
64                  phi = math.atan2(z+refSys.a*C*e2*math.sin(phi1), r) 
65                  if abs(phi1-phi)<1e-14: # phi is always of order 1 
66                          break 
67          return long/DEG, phi/DEG, r/math.cos(phi)-refSys.a*C 
68   
69   
70 -def geodToGeoc(long, phi, height, refSys=WGS84): 
71          """returns geocentric coordinates lambda, phi', rho for geodetic coordinates. 
72   
73          refSys defaults is the reference system the geodetic coordinates are 
74          expressed in. 
75   
76          height is in meter, long, phi in degrees. 
77   
78          See Astron. Almanac, Appendix K. 
79          """ 
80          long, phi = long*DEG, phi*DEG 
81          C, S = _getC_S(phi, refSys) 
82          rcp = (C+height/refSys.a)*math.cos(phi) 
83          rsp = (S+height/refSys.a)*math.sin(phi) 
84          rho = math.sqrt(rcp**2+rsp**2) 
85          phip = math.atan2(rsp, rcp) 
86          return long/DEG, phip/DEG, rho 
87