Coverage for python / lsst / pipe / tasks / ssp / photfit.py: 28%

1import numpy as np 

2from scipy.interpolate import CubicSpline 

3from scipy.optimize import leastsq 

4import warnings 

5 

6# Constants 

7 

8A = [3.332, 1.862] 

9B = [0.631, 1.218] 

10C = [0.986, 0.238] 

11 

12# values taken from sbpy for convenience 

13 

14alpha_12 = np.deg2rad([7.5, 30.0, 60, 90, 120, 150]) 

15 

16phi_1_sp = [7.5e-1, 3.3486016e-1, 1.3410560e-1, 5.1104756e-2, 2.1465687e-2, 3.6396989e-3] 

17phi_1_derivs = [-1.9098593, -9.1328612e-2] 

18 

19phi_2_sp = [9.25e-1, 6.2884169e-1, 3.1755495e-1, 1.2716367e-1, 2.2373903e-2, 1.6505689e-4] 

20phi_2_derivs = [-5.7295780e-1, -8.6573138e-8] 

21 

22alpha_3 = np.deg2rad([0.0, 0.3, 1.0, 2.0, 4.0, 8.0, 12.0, 20.0, 30.0]) 

23 

24phi_3_sp = [ 

25 1.0, 

26 8.3381185e-1, 

27 5.7735424e-1, 

28 4.2144772e-1, 

29 2.3174230e-1, 

30 1.0348178e-1, 

31 6.1733473e-2, 

32 1.6107006e-2, 

33 0.0, 

34] 

35 

36phi_3_derivs = [-1.0630097, 0] 

37 

38 

39phi_1 = CubicSpline(alpha_12, phi_1_sp, bc_type=((1, phi_1_derivs[0]), (1, phi_1_derivs[1]))) 

40phi_2 = CubicSpline(alpha_12, phi_2_sp, bc_type=((1, phi_2_derivs[0]), (1, phi_2_derivs[1]))) 

41phi_3 = CubicSpline(alpha_3, phi_3_sp, bc_type=((1, phi_3_derivs[0]), (1, phi_3_derivs[1]))) 

42 

43 

44def HG_model(phase, params): 

45 sin_a = np.sin(phase) 

46 tan_ah = np.tan(phase / 2) 

47 

48 W = np.exp(-90.56 * tan_ah * tan_ah) 

49 scale_sina = sin_a / (0.119 + 1.341 * sin_a - 0.754 * sin_a * sin_a) 

50 

51 phi_1_S = 1 - C[0] * scale_sina 

52 phi_2_S = 1 - C[1] * scale_sina 

53 

54 phi_1_L = np.exp(-A[0] * np.power(tan_ah, B[0])) 

55 phi_2_L = np.exp(-A[1] * np.power(tan_ah, B[1])) 

56 

57 phi_1 = W * phi_1_S + (1 - W) * phi_1_L 

58 phi_2 = W * phi_2_S + (1 - W) * phi_2_L 

59 return params[0] - 2.5 * np.log10((1 - params[1]) * phi_1 + (params[1]) * phi_2) 

60 

61 

62def HG1G2_model(phase, params): 

63 phi_1_ev = phi_1(phase) 

64 phi_2_ev = phi_2(phase) 

65 phi_3_ev = phi_3(phase) 

66 

67 msk = phase < 7.5 * np.pi / 180 

68 

69 phi_1_ev[msk] = 1 - 6 * phase[msk] / np.pi 

70 phi_2_ev[msk] = 1 - 9 * phase[msk] / (5 * np.pi) 

71 

72 phi_3_ev[phase > np.pi / 6] = 0 

73 

74 return params[0] - 2.5 * np.log10( 

75 params[1] * phi_1_ev + params[2] * phi_2_ev + (1 - params[1] - params[2]) * phi_3_ev 

76 ) 

77 

78 

79def HG12_model(phase, params): 

80 if params[1] >= 0.2: 

81 G1 = +0.9529 * params[1] + 0.02162 

82 G2 = -0.6125 * params[1] + 0.5572 

83 else: 

84 G1 = +0.7527 * params[1] + 0.06164 

85 G2 = -0.9612 * params[1] + 0.6270 

86 

87 return HG1G2_model(phase, [params[0], G1, G2]) 

88 

89 

90def HG12star_model(phase, params): 

91 G1 = 0 + params[1] * 0.84293649 

92 G2 = 0.53513350 - params[1] * 0.53513350 

93 

94 return HG1G2_model(phase, [params[0], G1, G2]) 

95 

96 

97def chi2(params, mag, phase, mag_err, model): 

98 pred = model(phase, params) 

99 return (mag - pred) / mag_err 

100 

101 

102def fit(mag, phase, sigma, model=HG12_model, params=[0.1]): 

103 phase = np.deg2rad(phase) 

104 

105 sol = leastsq(chi2, [mag[0]] + params, (mag, phase, sigma, model), full_output=True) 

106 

107 return sol 

108 

109 

110def fitHG12(mag, magSigma, phaseAngle, tdist, rdist): 

111 # This uses the HG12 function from Muinonen et al (2010) 

112 nobsv = len(mag) 

113 

114 # ensure these are plain ndarrays 

115 (mag, magSigma, phaseAngle, tdist, rdist) = map(np.asarray, (mag, magSigma, phaseAngle, tdist, rdist)) 

116 

117 # correct the mag to 1AU distance 

118 dmag = -5.0 * np.log10(tdist * rdist) 

119 mag = mag + dmag 

120 # phaseAngle = np.deg2rad(phaseAngle) 

121 

122 # fit, suppressing warnings 

123 with warnings.catch_warnings(): 

124 warnings.simplefilter("ignore") 

125 

126 sol = fit(mag, phaseAngle, magSigma) 

127 if sol[-1] in [1, 2, 3] and sol[1] is not None: # least squares solver flags 

128 chi2 = np.sum(sol[2]["fvec"] ** 2) 

129 H = sol[0][0] 

130 G = sol[0][1] 

131 H_err = np.sqrt(sol[1][0, 0]) 

132 G_err = np.sqrt(sol[1][1, 1]) 

133 HG_cov = sol[1][0, 1] 

134 

135 return H, G, H_err, G_err, HG_cov, chi2 / (nobsv - 2), nobsv 

136 else: 

137 # fit failed 

138 return (np.nan,) * 6 + (0,) 

139 

140 

141#################### 

142 

143 

144def phase_angle_deg(r_obj_sun, r_obs_sun): 

145 """ 

146 Compute phase angle (Sun–Object–Observer) in degrees. 

147 

148 Parameters 

149 ---------- 

150 r_obj_sun : array, shape (3,) or (3, N) 

151 Object position vector wrt Sun (Sun → object). 

152 r_obs_sun : array, shape (3,) or (3, N) 

153 Observer position vector wrt Sun (Sun → observer). 

154 

155 Returns 

156 ------- 

157 float or ndarray 

158 Phase angle(s) in degrees, in [0, 180]. 

159 """ 

160 r_obj_sun = np.asarray(r_obj_sun) 

161 r_obs_sun = np.asarray(r_obs_sun) 

162 

163 # Vectors at the object 

164 v_sun = -r_obj_sun # object → Sun 

165 v_obs = r_obs_sun - r_obj_sun # object → observer 

166 

167 # Dot products and norms along axis 0 

168 dot = np.sum(v_sun * v_obs, axis=0) 

169 norm_sun = np.linalg.norm(v_sun, axis=0) 

170 norm_obs = np.linalg.norm(v_obs, axis=0) 

171 

172 cosang = dot / (norm_sun * norm_obs) 

173 cosang = np.clip(cosang, -1.0, 1.0) 

174 

175 return np.degrees(np.arccos(cosang)) 

176 

177 

178def hg_V_mag(H, G, r, delta, phase_deg): 

179 """ 

180 Compute apparent V magnitude from the IAU H–G system. 

181 

182 Parameters 

183 ---------- 

184 H : float or ndarray 

185 Absolute magnitude (V-band). 

186 G : float or ndarray 

187 Slope parameter. 

188 r : float or ndarray 

189 Heliocentric distance in AU. 

190 delta : float or ndarray 

191 Observer distance (Δ) in AU. 

192 phase_deg : float or ndarray 

193 Phase angle in degrees. 

194 """ 

195 a = np.radians(phase_deg) / 2.0 

196 

197 # Phase functions 

198 phi1 = np.exp(-3.33 * np.tan(a) ** 0.63) 

199 phi2 = np.exp(-1.87 * np.tan(a) ** 1.22) 

200 

201 phi = (1 - G) * phi1 + G * phi2 

202 

203 return H + 5 * np.log10(r * delta) - 2.5 * np.log10(phi)