def D_Hplus1(self, x, y, dunits='km', doseunits='Sv'):
"""Returns dose rate at x, y at 1 hour after burst. This value includes dose rate from all activity that WILL be deposited at that location, not just that that has arrived by H+1 hr."""
rx, ry = self.translation * (convert_units(x, dunits, 'mi'), convert_units(y, dunits, 'mi')) * ~Affine.rotation(-self.wd + 270)
f_x = self.yld * 2e6 * self.phi(rx) * self.g(rx) * self.ff
s_y = np.sqrt(self.s_02 + ((8 * abs(rx + 2 * self.s_x) * self.s_02) / self.L) + (2 * (self.s_x * self.T_c * self.s_h * self.shear)**2 / self.L_2) + (((rx + 2 * self.s_x) * self.L_0 * self.T_c * self.s_h * self.shear)**2 / self.L**4))
a_2 = 1 / (1 + ((0.001 * self.H_c * self.wind) / self.s_0) * (1 - norm.cdf(2 * x / self.wind)))
f_y = np.exp(-0.5 * (ry / (a_2 * s_y))**2) / (2.5066282746310002 * s_y)
return convert_units(f_x * f_y, 'Roentgen', doseunits)
评论列表
文章目录