maths_fns.pseudo

31 """Return n! 32 33 @param n: The number to return the factorial of. 34 @type n: int 35 @return: n! 36 @rtype: int 37 """ 38 39 # A list of factorials. 40 fact = [ 41 1, 42 1, 43 2, 44 6, 45 24, 46 120, 47 720, 48 5040, 49 40320, 50 362880, 51 3628800, 52 39916800, 53 479001600, 54 6227020800, 55 87178291200, 56 1307674368000, 57 20922789888000, 58 355687428096000, 59 6402373705728000, 60 121645100408832000, 61 2432902008176640000, 62 51090942171709440000, 63 1124000727777607680000, 64 25852016738884976640000, 65 620448401733239439360000, 66 15511210043330985984000000, 67 403291461126605635584000000, 68 10888869450418352160768000000, 69 304888344611713860501504000000, 70 8841761993739701954543616000000, 71 265252859812191058636308480000000, 72 8222838654177922817725562880000000, 73 263130836933693530167218012160000000, 74 8683317618811886495518194401280000000, 75 295232799039604140847618609643520000000, 76 10333147966386144929666651337523200000000, 77 371993326789901217467999448150835200000000, 78 13763753091226345046315979581580902400000000, 79 523022617466601111760007224100074291200000000, 80 20397882081197443358640281739902897356800000000 81 ] 82 83 # Return the value. 84 return fact[n]

88 """The 2D pseudo-elliptic cosine function. 89 90 This is defined as:: 91 92 int_-pi^pi (1 - cos(theta_max)).dphi, 93 94 where:: 95 96 theta_max = 1 / sqrt(cos**2(phi)/x**2 + sin**2(phi)/y**2). 97 98 The function is approximated as:: 99 100 _inf 101 \ (-1)**n 102 2*pi*x*y > ----------------- fn(x, y), 103 /___ 4**n * (2n + 2)! 104 n=0 105 106 and where the individual functions are:: 107 108 n: 0; fn(x,y) = 1 109 n: 1; fn(x,y) = 2*a + 2*b 110 n: 2; fn(x,y) = 6*a^2 + 4*a*b + 6*b^2 111 n: 3; fn(x,y) = 20*a^3 + 12*a^2*b + 12*a*b^2 + 20*b^3 112 n: 4; fn(x,y) = 70*a^4 + 40*a^3*b + 36*a^2*b^2 + 40*a*b^3 + 70*b^4 113 n: 5; fn(x,y) = 252*a^5 + 140*a^4*b + 120*a^3*b^2 + 120*a^2*b^3 + 140*a*b^4 + 252*b^5 114 n: 6; fn(x,y) = 924*a^6 + 504*a^5*b + 420*a^4*b^2 + 400*a^3*b^3 + 420*a^2*b^4 + 504*a*b^5 + 924*b^6 115 n: 7; fn(x,y) = 3432*a^7 + 1848*a^6*b + 1512*a^5*b^2 + 1400*a^4*b^3 + 1400*a^3*b^4 + 1512*a^2*b^5 + 1848*a*b^6 + 3432*b^7 116 n: 8; fn(x,y) = 12870*a^8 + 6864*a^7*b + 5544*a^6*b^2 + 5040*a^5*b^3 + 4900*a^4*b^4 + 5040*a^3*b^5 + 5544*a^2*b^6 + 6864*a*b^7 + 12870*b^8 117 n: 9; fn(x,y) = 48620*a^9 + 25740*a^8*b + 20592*a^7*b^2 + 18480*a^6*b^3 + 17640*a^5*b^4 + 17640*a^4*b^5 + 18480*a^3*b^6 + 20592*a^2*b^7 + 25740*a*b^8 + 48620*b^9 118 n: 10; fn(x,y) = 184756*a^10 + 97240*a^9*b + 77220*a^8*b^2 + 68640*a^7*b^3 + 64680*a^6*b^4 + 63504*a^5*b^5 + 64680*a^4*b^6 + 68640*a^3*b^7 + 77220*a^2*b^8 + 97240*a*b^9 + 184756*b^10 119 n: 11; fn(x,y) = 705432*a^11 + 369512*a^10*b + 291720*a^9*b^2 + 257400*a^8*b^3 + 240240*a^7*b^4 + 232848*a^6*b^5 + 232848*a^5*b^6 + 240240*a^4*b^7 + 257400*a^3*b^8 + 291720*a^2*b^9 + 369512*a*b^10 + 705432*b^11 120 121 122 @param x: The theta_x cone angle (in rad). 123 @type x: float 124 @param y: The theta_y cone angle (in rad). 125 @type y: float 126 @return: The approximate function value. 127 @rtype: float 128 """ 129 130 # Init. 131 fn = 0.0 132 fact = 1.0 133 134 # The x powers. 135 a = x**2 136 a2 = a**2 137 a3 = a**3 138 a4 = a**4 139 a5 = a**5 140 a6 = a**6 141 a7 = a**7 142 a8 = a**8 143 a9 = a**9 144 a10 = a**10 145 146 # The y powers. 147 b = y**2 148 b2 = b**2 149 b3 = b**3 150 b4 = b**4 151 b5 = b**5 152 b6 = b**6 153 b7 = b**7 154 b8 = b**8 155 b9 = b**9 156 b10 = b**10 157 158 # Term 0. 159 f0 = 1.0 160 fn = fn + f0 / factorial(2) 161 162 # Term 1. 163 f1 = 2*a + 2*b 164 fn = fn - f1 / 4.0 / factorial(4) 165 166 # Term 2. 167 f2 = 6*a2 + 4*a*b + 6*b2 168 fn = fn + f2 / 16.0 / factorial(6) 169 170 # Term 3. 171 f3 = 20*a3 + 12*a2*b + 12*a*b2 + 20*b3 172 fn = fn - f3 / 64.0 / factorial(8) 173 174 # Term 4. 175 f4 = 70*a4 + 40*a3*b + 36*a2*b2 + 40*a*b3 + 70*b4 176 fn = fn + f4 / 256.0 / factorial(10) 177 178 # Term 5. 179 f5 = 252*a5 + 140*a4*b + 120*a3*b2 + 120*a2*b3 + 140*a*b4 + 252*b5 180 fn = fn - f5 / 1024.0 / factorial(12) 181 182 # Term 6. 183 f6 = 924*a6 + 504*a5*b + 420*a4*b2 + 400*a3*b3 + 420*a2*b4 + 504*a*b5 + 924*b6 184 fn = fn + f6 / 4096.0 / factorial(14) 185 186 # Term 7. 187 f7 = 3432*a7 + 1848*a6*b + 1512*a5*b2 + 1400*a4*b3 + 1400*a3*b4 + 1512*a2*b5 + 1848*a*b6 + 3432*b7 188 fn = fn - f7 / 16384.0 / factorial(16) 189 190 # Term 8. 191 f8 = 12870*a8 + 6864*a7*b + 5544*a6*b2 + 5040*a5*b3 + 4900*a4*b4 + 5040*a3*b5 + 5544*a2*b6 + 6864*a*b7 + 12870*b8 192 fn = fn + f8 / 65536.0 / factorial(18) 193 194 # Term 9. 195 f9 = 48620*a9 + 25740*a8*b + 20592*a7*b2 + 18480*a6*b3 + 17640*a5*b4 + 17640*a4*b5 + 18480*a3*b6 + 20592*a2*b7 + 25740*a*b8 + 48620*b9 196 fn = fn - f9 / 262144.0 / factorial(20) 197 198 # Term 10. 199 f10 = 184756*a10 + 97240*a9*b + 77220*a8*b2 + 68640*a7*b3 + 64680*a6*b4 + 63504*a5*b5 + 64680*a4*b6 + 68640*a3*b7 + 77220*a2*b8 + 97240*a*b9 + 184756*b10 200 fn = fn + f10 / 1048576.0 / factorial(22) 201 202 # Common factor. 203 fn = 2.0 * pi * x * y * fn 204 205 # Return the approximate function value. 206 return fn

Source Code for Module maths_fns.pseudo_ellipse