Changeset 3f8584a2 in sasmodels
 Timestamp:
 Jul 27, 2016 3:53:07 AM (8 years ago)
 Branches:
 master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 c047acf
 Parents:
 56547a8
 Location:
 sasmodels/models/lib
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/models/lib/sas_J0.c
r0278e3f r3f8584a2 53 53 * except YP, YQ which are designed for absolute error. */ 54 54 55 #if FLOAT_SIZE>4 56 //Cephes double precission 55 57 56 58 constant double PPJ0[8] = { … … 144 146 }; 145 147 146 constant double MOJ0[8] = { 148 double j0(double x); 149 double j0(double x) 150 { 151 double w, z, p, q, xn; 152 153 //const double TWOOPI = 6.36619772367581343075535E1; 154 const double SQ2OPI = 7.9788456080286535587989E1; 155 const double PIO4 = 7.85398163397448309616E1; 156 157 const double DR1 = 5.78318596294678452118E0; 158 const double DR2 = 3.04712623436620863991E1; 159 160 161 if( x < 0 ) 162 x = x; 163 164 if( x <= 5.0 ) { 165 z = x * x; 166 if( x < 1.0e5 ) 167 return( 1.0  z/4.0 ); 168 169 p = (z  DR1) * (z  DR2); 170 p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 ); 171 return( p ); 172 } 173 174 w = 5.0/x; 175 q = 25.0/(x*x); 176 p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 ); 177 q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 ); 178 xn = x  PIO4; 179 180 double sn, cn; 181 SINCOS(xn, sn, cn); 182 p = p * cn  w * q * sn; 183 184 return( p * SQ2OPI / sqrt(x) ); 185 } 186 #else 187 188 constant float MOJ0[8] = { 147 189 6.838999669318810E002, 148 190 1.864949361379502E001, … … 155 197 }; 156 198 157 constant doublePHJ0[8] = {199 constant float PHJ0[8] = { 158 200 3.242077816988247E+001, 159 201 3.630592630518434E+001, … … 166 208 }; 167 209 168 constant doubleJPJ0[8] = {210 constant float JPJ0[8] = { 169 211 6.068350350393235E008, 170 212 6.388945720783375E006, … … 177 219 }; 178 220 179 double sas_J0(double x); 180 double sas_J0(double x) 221 //Cephes single precission 222 float j0f(float x); 223 float j0f(float x) 181 224 { 182 183 //Cephes single precission 184 #if FLOAT_SIZE>4 185 double w, z, p, q, xn; 186 187 //const double TWOOPI = 6.36619772367581343075535E1; 188 const double SQ2OPI = 7.9788456080286535587989E1; 189 const double PIO4 = 7.85398163397448309616E1; 190 191 const double DR1 = 5.78318596294678452118E0; 192 const double DR2 = 3.04712623436620863991E1; 193 194 195 if( x < 0 ) 196 x = x; 197 198 if( x <= 5.0 ) { 199 z = x * x; 200 if( x < 1.0e5 ) 201 return( 1.0  z/4.0 ); 202 203 p = (z  DR1) * (z  DR2); 204 p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 ); 205 return( p ); 206 } 207 208 w = 5.0/x; 209 q = 25.0/(x*x); 210 p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 ); 211 q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 ); 212 xn = x  PIO4; 213 214 double sn, cn; 215 SINCOS(xn, sn, cn); 216 p = p * cn  w * q * sn; 217 218 return( p * SQ2OPI / sqrt(x) ); 219 //Cephes single precission 220 #else 221 double xx, w, z, p, q, xn; 225 float xx, w, z, p, q, xn; 222 226 223 227 //const double YZ1 = 0.43221455686510834878; 224 228 //const double YZ2 = 22.401876406482861405; 225 229 //const double YZ3 = 64.130620282338755553; 226 const doubleDR1 = 5.78318596294678452118;227 const doublePIO4F = 0.7853981633974483096;230 const float DR1 = 5.78318596294678452118; 231 const float PIO4F = 0.7853981633974483096; 228 232 229 233 if( x < 0 ) 230 234 xx = x; 231 235 else 232 236 xx = x; 233 237 234 238 if( x <= 2.0 ) { 235 236 237 238 239 240 241 239 z = xx * xx; 240 if( x < 1.0e3 ) 241 return( 1.0  0.25*z ); 242 243 p = (zDR1) * polevl( z, JPJ0, 4); 244 return( p ); 245 } 242 246 243 247 q = 1.0/x; … … 249 253 p = p * cos(xn + xx); 250 254 return(p); 255 } 251 256 #endif 252 257 253 } 254 255 258 #if FLOAT_SIZE>4 259 #define sas_J0 j0 260 #else 261 #define sas_J0 j0f 262 #endif 
sasmodels/models/lib/sas_J1.c
r0278e3f r3f8584a2 40 40 */ 41 41 42 #if FLOAT_SIZE>4 43 44 //Cephes double pression function 42 45 constant double RPJ1[8] = { 43 46 8.99971225705559398224E8, … … 102 105 0.0 }; 103 106 104 constant double JPJ1[8] = { 107 double j1(double x); 108 double j1(double x) 109 { 110 111 double w, z, p, q, xn; 112 113 const double Z1 = 1.46819706421238932572E1; 114 const double Z2 = 4.92184563216946036703E1; 115 const double THPIO4 = 2.35619449019234492885; 116 const double SQ2OPI = 0.79788456080286535588; 117 118 w = x; 119 if( x < 0 ) 120 w = x; 121 122 if( w <= 5.0 ) { 123 z = x * x; 124 w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 125 w = w * x * (z  Z1) * (z  Z2); 126 return( w ); 127 } 128 129 w = 5.0/x; 130 z = w * w; 131 132 p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 133 q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 134 135 xn = x  THPIO4; 136 137 double sn, cn; 138 SINCOS(xn, sn, cn); 139 p = p * cn  w * q * sn; 140 141 return( p * SQ2OPI / sqrt(x) ); 142 } 143 144 #else 145 //Single precission version of cephes 146 constant float JPJ1[8] = { 105 147 4.878788132172128E009, 106 148 6.009061827883699E007, … … 113 155 }; 114 156 115 constant doubleMO1J1[8] = {157 constant float MO1J1[8] = { 116 158 6.913942741265801E002, 117 159 2.284801500053359E001, … … 124 166 }; 125 167 126 constant doublePH1J1[8] = {168 constant float PH1J1[8] = { 127 169 4.497014141919556E+001, 128 170 5.073465654089319E+001, … … 135 177 }; 136 178 137 double sas_J1(doublex);138 double sas_J1(doublex)179 float j1f(float x); 180 float j1f(float x) 139 181 { 140 182 141 //Cephes double pression function 142 #if FLOAT_SIZE>4 143 144 double w, z, p, q, xn; 145 146 const double Z1 = 1.46819706421238932572E1; 147 const double Z2 = 4.92184563216946036703E1; 148 const double THPIO4 = 2.35619449019234492885; 149 const double SQ2OPI = 0.79788456080286535588; 150 151 w = x; 152 if( x < 0 ) 153 w = x; 154 155 if( w <= 5.0 ) 156 { 157 z = x * x; 158 w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 159 w = w * x * (z  Z1) * (z  Z2); 160 return( w ); 161 } 162 163 w = 5.0/x; 164 z = w * w; 165 166 p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 167 q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 168 169 xn = x  THPIO4; 170 171 double sn, cn; 172 SINCOS(xn, sn, cn); 173 p = p * cn  w * q * sn; 174 175 return( p * SQ2OPI / sqrt(x) ); 176 177 178 //Single precission version of cephes 179 #else 180 double xx, w, z, p, q, xn; 181 182 const double Z1 = 1.46819706421238932572E1; 183 const double THPIO4F = 2.35619449019234492885; /* 3*pi/4 */ 183 float xx, w, z, p, q, xn; 184 185 const float Z1 = 1.46819706421238932572E1; 186 const float THPIO4F = 2.35619449019234492885; /* 3*pi/4 */ 184 187 185 188 186 189 xx = x; 187 190 if( xx < 0 ) 188 xx = x; 189 190 if( xx <= 2.0 ) 191 { 192 z = xx * xx; 193 p = (zZ1) * xx * polevl( z, JPJ1, 4 ); 194 return( p ); 195 } 191 xx = x; 192 193 if( xx <= 2.0 ) { 194 z = xx * xx; 195 p = (zZ1) * xx * polevl( z, JPJ1, 4 ); 196 return( p ); 197 } 196 198 197 199 q = 1.0/x; … … 204 206 205 207 return(p); 208 } 206 209 #endif 207 } 210 211 #if FLOAT_SIZE>4 212 #define sas_J1 j1 213 #else 214 #define sas_J1 j1f 215 #endif 208 216 209 217 //Finally J1c function that equals 2*J1(x)/x 
sasmodels/models/lib/sas_JN.c
r4937980 r3f8584a2 48 48 */ 49 49 50 double sas_JN( int n, double x ); 51 52 double sas_JN( int n, double x ) { 53 50 #if FLOAT_SIZE > 4 51 52 double jn( int n, double x ); 53 double jn( int n, double x ) { 54 55 // PAK: seems to be machine epsilon/2 54 56 const double MACHEP = 1.11022302462515654042E16; 55 57 double pkm2, pkm1, pk, xk, r, ans; … … 57 59 58 60 if( n < 0 ) { 59 60 61 62 63 64 } 65 else66 sign = 1; 61 n = n; 62 if( (n & 1) == 0 ) /* 1**n */ 63 sign = 1; 64 else 65 sign = 1; 66 } else { 67 sign = 1; 68 } 67 69 68 70 if( x < 0.0 ) { 69 70 71 72 71 if( n & 1 ) 72 sign = sign; 73 x = x; 74 } 73 75 74 76 if( n == 0 ) 75 return( sign * sas_J0(x) );77 return( sign * j0(x) ); 76 78 if( n == 1 ) 77 return( sign * sas_J1(x) );79 return( sign * j1(x) ); 78 80 if( n == 2 ) 79 return( sign * (2.0 * sas_J1(x) / x  sas_J0(x)) );81 return( sign * (2.0 * j1(x) / x  j0(x)) ); 80 82 81 83 if( x < MACHEP ) 82 return( 0.0 ); 83 84 #if FLOAT_SIZE > 4 85 k = 53; 86 #else 87 k = 24; 88 #endif 89 84 return( 0.0 ); 85 86 k = 53; 90 87 pk = 2 * (n + k); 91 88 ans = pk; … … 93 90 94 91 do { 95 96 97 92 pk = 2.0; 93 ans = pk  (xk/ans); 94 } while( k > 0 ); 98 95 99 96 /* backward recurrence */ … … 101 98 pk = 1.0; 102 99 103 #if FLOAT_SIZE > 4 104 ans = x/ans; 105 pkm1 = 1.0/ans; 106 107 k = n1; 108 r = 2 * k; 109 110 do { 111 pkm2 = (pkm1 * r  pk * x) / x; 112 pk = pkm1; 113 pkm1 = pkm2; 114 r = 2.0; 115 } while( k > 0 ); 116 117 if( fabs(pk) > fabs(pkm1) ) 118 ans = sas_J1(x)/pk; 100 ans = x/ans; 101 pkm1 = 1.0/ans; 102 103 k = n1; 104 r = 2 * k; 105 106 do { 107 pkm2 = (pkm1 * r  pk * x) / x; 108 pk = pkm1; 109 pkm1 = pkm2; 110 r = 2.0; 111 } while( k > 0 ); 112 113 if( fabs(pk) > fabs(pkm1) ) 114 ans = j1(x)/pk; 115 else 116 ans = j0(x)/pkm1; 117 118 return( sign * ans ); 119 } 120 121 #else 122 123 float jnf(int n, float x) 124 { 125 // PAK: seems to be machine epsilon/2 126 const double MACHEP = 5.9604645e08; 127 float pkm2, pkm1, pk, xk, r, ans; 128 int k, sign; 129 130 if( n < 0 ) { 131 n = n; 132 if( (n & 1) == 0 ) /* 1**n */ 133 sign = 1; 119 134 else 120 ans = sas_J0(x)/pkm1; 121 122 return( sign * ans ); 123 124 #else 125 const double xinv = 1.0/x; 126 pkm1 = ans * xinv; 127 k = n1; 128 r = (float )(2 * k); 129 130 do { 131 pkm2 = (pkm1 * r  pk * x) * xinv; 132 pk = pkm1; 133 pkm1 = pkm2; 134 r = 2.0; 135 } 136 while( k > 0 ); 137 138 r = pk; 139 if( r < 0 ) 140 r = r; 141 ans = pkm1; 142 if( ans < 0 ) 143 ans = ans; 144 145 if( r > ans ) /* if( fabs(pk) > fabs(pkm1) ) */ 146 ans = sign * sas_J1(x)/pk; 147 else 148 ans = sign * sas_J0(x)/pkm1; 149 return( ans ); 150 #endif 135 sign = 1; 136 } else { 137 sign = 1; 138 } 139 140 if( x < 0.0 ) { 141 if( n & 1 ) 142 sign = sign; 143 x = x; 144 } 145 146 if( n == 0 ) 147 return( sign * j0f(x) ); 148 if( n == 1 ) 149 return( sign * j1f(x) ); 150 if( n == 2 ) 151 return( sign * (2.0 * j1f(x) / x  j0f(x)) ); 152 153 if( x < MACHEP ) 154 return( 0.0 ); 155 156 k = 24; 157 pk = 2 * (n + k); 158 ans = pk; 159 xk = x * x; 160 161 do { 162 pk = 2.0; 163 ans = pk  (xk/ans); 164 } while( k > 0 ); 165 166 /* backward recurrence */ 167 168 pk = 1.0; 169 170 const float xinv = 1.0/x; 171 pkm1 = ans * xinv; 172 k = n1; 173 r = (float )(2 * k); 174 175 do { 176 pkm2 = (pkm1 * r  pk * x) * xinv; 177 pk = pkm1; 178 pkm1 = pkm2; 179 r = 2.0; 180 } while( k > 0 ); 181 182 r = pk; 183 if( r < 0 ) 184 r = r; 185 ans = pkm1; 186 if( ans < 0 ) 187 ans = ans; 188 189 if( r > ans ) /* if( fabs(pk) > fabs(pkm1) ) */ 190 ans = sign * j1f(x)/pk; 191 else 192 ans = sign * j0f(x)/pkm1; 193 return( ans ); 151 194 } 152 195 #endif 196 197 #if FLOAT_SIZE>4 198 #define sas_JN jn 199 #else 200 #define sas_JN jnf 201 #endif
Note: See TracChangeset
for help on using the changeset viewer.