1# -*- coding: utf-8 -*- 

2from ._constants import * 

3import numpy as np 

4 

5 

6__all__ = [ 

7 "get_ortho_latitude_lines", 

8 "get_ortho_longitude_lines", 

9 "get_projection", 

10] 

11 

12 

13def RAxisAngle(axis=[0, 1, 0], theta=0): 

14 """ 

15 

16 """ 

17 cost = np.cos(theta) 

18 sint = np.sin(theta) 

19 

20 return np.reshape( 

21 [ 

22 cost + axis[0] * axis[0] * (1 - cost), 

23 axis[0] * axis[1] * (1 - cost) - axis[2] * sint, 

24 axis[0] * axis[2] * (1 - cost) + axis[1] * sint, 

25 axis[1] * axis[0] * (1 - cost) + axis[2] * sint, 

26 cost + axis[1] * axis[1] * (1 - cost), 

27 axis[1] * axis[2] * (1 - cost) - axis[0] * sint, 

28 axis[2] * axis[0] * (1 - cost) - axis[1] * sint, 

29 axis[2] * axis[1] * (1 - cost) + axis[0] * sint, 

30 cost + axis[2] * axis[2] * (1 - cost), 

31 ], 

32 [3, 3], 

33 ) 

34 

35 

36def get_ortho_latitude_lines(inc=np.pi / 2, obl=0, dlat=np.pi / 6, npts=1000): 

37 """ 

38 

39 """ 

40 # Angular quantities 

41 ci = np.cos(inc) 

42 si = np.sin(inc) 

43 co = np.cos(obl) 

44 so = np.sin(obl) 

45 

46 # Latitude lines 

47 res = [] 

48 latlines = np.arange(-np.pi / 2, np.pi / 2, dlat)[1:] 

49 for lat in latlines: 

50 

51 # Figure out the equation of the ellipse 

52 y0 = np.sin(lat) * si 

53 a = np.cos(lat) 

54 b = a * ci 

55 x = np.linspace(-a, a, npts) 

56 y1 = y0 - b * np.sqrt(1 - (x / a) ** 2) 

57 y2 = y0 + b * np.sqrt(1 - (x / a) ** 2) 

58 

59 # Mask lines on the backside 

60 if si != 0: 

61 if inc > np.pi / 2: 

62 ymax = y1[np.argmax(x ** 2 + y1 ** 2)] 

63 y1[y1 < ymax] = np.nan 

64 ymax = y2[np.argmax(x ** 2 + y2 ** 2)] 

65 y2[y2 < ymax] = np.nan 

66 else: 

67 ymax = y1[np.argmax(x ** 2 + y1 ** 2)] 

68 y1[y1 > ymax] = np.nan 

69 ymax = y2[np.argmax(x ** 2 + y2 ** 2)] 

70 y2[y2 > ymax] = np.nan 

71 

72 # Rotate them 

73 for y in (y1, y2): 

74 xr = -x * co - y * so 

75 yr = -x * so + y * co 

76 res.append((xr, yr)) 

77 

78 return res 

79 

80 

81def get_ortho_longitude_lines( 

82 inc=np.pi / 2, obl=0, theta=0, dlon=np.pi / 6, npts=1000 

83): 

84 """ 

85 

86 """ 

87 

88 # Angular quantities 

89 ci = np.cos(inc) 

90 si = np.sin(inc) 

91 co = np.cos(obl) 

92 so = np.sin(obl) 

93 

94 # Are we (essentially) equator-on? 

95 equator_on = (inc > 88 * np.pi / 180) and (inc < 92 * np.pi / 180) 

96 

97 # Longitude grid lines 

98 res = [] 

99 if equator_on: 

100 offsets = np.arange(-np.pi / 2, np.pi / 2, dlon) 

101 else: 

102 offsets = np.arange(0, 2 * np.pi, dlon) 

103 

104 for offset in offsets: 

105 

106 # Super hacky, sorry. This can probably 

107 # be coded up more intelligently. 

108 if equator_on: 

109 sgns = [1] 

110 if np.cos(theta + offset) >= 0: 

111 bsgn = 1 

112 else: 

113 bsgn = -1 

114 else: 

115 bsgn = 1 

116 if np.cos(theta + offset) >= 0: 

117 sgns = np.array([1, -1]) 

118 else: 

119 sgns = np.array([-1, 1]) 

120 

121 for lon, sgn in zip([0, np.pi], sgns): 

122 

123 # Viewed at i = 90 

124 y = np.linspace(-1, 1, npts) 

125 b = bsgn * np.sin(lon - theta - offset) 

126 x = b * np.sqrt(1 - y ** 2) 

127 z = sgn * np.sqrt(np.abs(1 - x ** 2 - y ** 2)) 

128 

129 if equator_on: 

130 

131 pass 

132 

133 else: 

134 

135 # Rotate by the inclination 

136 R = RAxisAngle([1, 0, 0], np.pi / 2 - inc) 

137 v = np.vstack( 

138 (x.reshape(1, -1), y.reshape(1, -1), z.reshape(1, -1)) 

139 ) 

140 x, y, _ = np.dot(R, v) 

141 

142 # Mask lines on the backside 

143 if si != 0: 

144 if inc < np.pi / 2: 

145 imax = np.argmax(x ** 2 + y ** 2) 

146 y[: imax + 1] = np.nan 

147 else: 

148 imax = np.argmax(x ** 2 + y ** 2) 

149 y[imax:] = np.nan 

150 

151 # Rotate by the obliquity 

152 xr = -x * co - y * so 

153 yr = -x * so + y * co 

154 res.append((xr, yr)) 

155 

156 return res 

157 

158 

159def get_projection(projection): 

160 """ 

161 

162 """ 

163 if projection.lower().startswith("rect"): 

164 projection = STARRY_RECTANGULAR_PROJECTION 

165 elif projection.lower().startswith("ortho"): 

166 projection = STARRY_ORTHOGRAPHIC_PROJECTION 

167 else: 

168 raise ValueError("Unknown map projection.") 

169 return projection