Changeset 9c77bcf in sasmodels
- Timestamp:
- Mar 18, 2016 12:49:31 PM (9 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 3a45c2c
- Parents:
- af92b73 (diff), 9f7a852 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/J1.py
raceac39 r95ce773 50 50 model = core.load_model('bessel', dtype=dtype) 51 51 calculator = direct_model.DirectModel(data.empty_data1D(x), model) 52 #ret = calculator(background=0) 53 #print ret 52 54 53 return calculator(background=0) 55 54 -
explore/J1c.py
rcbd37a7 r95ce773 2 2 Show numerical precision of $2 J_1(x)/x$. 3 3 """ 4 import sys; sys.path.insert(0, '..') 4 5 5 6 import numpy as np 6 from sympy.mpmath import mp 7 try: 8 from mpmath import mp 9 except: 10 # CRUFT: mpmath split out into its own package 11 from sympy.mpmath import mp 7 12 #import matplotlib; matplotlib.use('TkAgg') 8 13 import pylab -
sasmodels/models/bessel.py
ra5af4e1 r95ce773 78 78 79 79 Iq = """ 80 return J1 (q);80 return J1c(q); 81 81 """ 82 82 -
sasmodels/models/lib/j0_cephes.c
r094e320 r95ce773 53 53 * except YP, YQ which are designed for absolute error. */ 54 54 55 double j0( double ); 56 double j0(double x) { 55 double J0(double x) { 57 56 58 57 //Cephes single precission … … 60 59 double w, z, p, q, xn; 61 60 62 const double TWOOPI = 6.36619772367581343075535E-1;61 //const double TWOOPI = 6.36619772367581343075535E-1; 63 62 const double SQ2OPI = 7.9788456080286535587989E-1; 64 63 const double PIO4 = 7.85398163397448309616E-1; … … 67 66 const double DR2 = 3.04712623436620863991E1; 68 67 69 const double PP[8] = {70 7.96936729297347051624E-4,71 8.28352392107440799803E-2,72 1.23953371646414299388E0,73 5.44725003058768775090E0,74 8.74716500199817011941E0,75 5.30324038235394892183E0,76 9.99999999999999997821E-1,77 0.078 };79 80 const double PQ[8] = {81 9.24408810558863637013E-4,82 8.56288474354474431428E-2,83 1.25352743901058953537E0,84 5.47097740330417105182E0,85 8.76190883237069594232E0,86 5.30605288235394617618E0,87 1.00000000000000000218E0,88 0.089 };90 91 const double QP[8] = {92 -1.13663838898469149931E-2,93 -1.28252718670509318512E0,94 -1.95539544257735972385E1,95 -9.32060152123768231369E1,96 -1.77681167980488050595E2,97 -1.47077505154951170175E2,98 -5.14105326766599330220E1,99 -6.05014350600728481186E0,100 };101 102 const double QQ[8] = {103 /* 1.00000000000000000000E0,*/104 6.43178256118178023184E1,105 8.56430025976980587198E2,106 3.88240183605401609683E3,107 7.24046774195652478189E3,108 5.93072701187316984827E3,109 2.06209331660327847417E3,110 2.42005740240291393179E2,111 };112 113 const double YP[8] = {114 1.55924367855235737965E4,115 -1.46639295903971606143E7,116 5.43526477051876500413E9,117 -9.82136065717911466409E11,118 8.75906394395366999549E13,119 -3.46628303384729719441E15,120 4.42733268572569800351E16,121 -1.84950800436986690637E16,122 };123 124 const double YQ[7] = {125 /* 1.00000000000000000000E0,*/126 1.04128353664259848412E3,127 6.26107330137134956842E5,128 2.68919633393814121987E8,129 8.64002487103935000337E10,130 2.02979612750105546709E13,131 3.17157752842975028269E15,132 2.50596256172653059228E17,133 };134 135 const double RP[8] = {136 -4.79443220978201773821E9,137 1.95617491946556577543E12,138 -2.49248344360967716204E14,139 9.70862251047306323952E15,140 0.0,141 0.0,142 0.0,143 0.0144 };145 146 const double RQ[8] = {147 /* 1.00000000000000000000E0,*/148 4.99563147152651017219E2,149 1.73785401676374683123E5,150 4.84409658339962045305E7,151 1.11855537045356834862E10,152 2.11277520115489217587E12,153 3.10518229857422583814E14,154 3.18121955943204943306E16,155 1.71086294081043136091E18,156 };157 68 158 69 if( x < 0 ) … … 165 76 166 77 p = (z - DR1) * (z - DR2); 167 p = p * polevl( z, RP , 3)/p1evl( z, RQ, 8 );78 p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 ); 168 79 return( p ); 169 80 } … … 171 82 w = 5.0/x; 172 83 q = 25.0/(x*x); 173 p = polevl( q, PP , 6)/polevl( q, PQ, 6 );174 q = polevl( q, QP , 7)/p1evl( q, QQ, 7 );84 p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 ); 85 q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 ); 175 86 xn = x - PIO4; 176 87 … … 184 95 double xx, w, z, p, q, xn; 185 96 186 const double YZ1 = 0.43221455686510834878;187 const double YZ2 = 22.401876406482861405;188 const double YZ3 = 64.130620282338755553;97 //const double YZ1 = 0.43221455686510834878; 98 //const double YZ2 = 22.401876406482861405; 99 //const double YZ3 = 64.130620282338755553; 189 100 const double DR1 = 5.78318596294678452118; 190 101 const double PIO4F = 0.7853981633974483096; 191 192 double MO[8] = {193 -6.838999669318810E-002,194 1.864949361379502E-001,195 -2.145007480346739E-001,196 1.197549369473540E-001,197 -3.560281861530129E-003,198 -4.969382655296620E-002,199 -3.355424622293709E-006,200 7.978845717621440E-001201 };202 203 double PH[8] = {204 3.242077816988247E+001,205 -3.630592630518434E+001,206 1.756221482109099E+001,207 -4.974978466280903E+000,208 1.001973420681837E+000,209 -1.939906941791308E-001,210 6.490598792654666E-002,211 -1.249992184872738E-001212 };213 214 double YP[8] = {215 9.454583683980369E-008,216 -9.413212653797057E-006,217 5.344486707214273E-004,218 -1.584289289821316E-002,219 1.707584643733568E-001,220 0.0,221 0.0,222 0.0223 };224 225 double JP[8] = {226 -6.068350350393235E-008,227 6.388945720783375E-006,228 -3.969646342510940E-004,229 1.332913422519003E-002,230 -1.729150680240724E-001,231 0.0,232 0.0,233 0.0234 };235 102 236 103 if( x < 0 ) … … 244 111 return( 1.0 - 0.25*z ); 245 112 246 p = (z-DR1) * polevl( z, JP , 4);113 p = (z-DR1) * polevl( z, JPJ0, 4); 247 114 return( p ); 248 115 } … … 251 118 w = sqrt(q); 252 119 253 p = w * polevl( q, MO , 7);120 p = w * polevl( q, MOJ0, 7); 254 121 w = q*q; 255 xn = q * polevl( w, PH , 7) - PIO4F;122 xn = q * polevl( w, PHJ0, 7) - PIO4F; 256 123 p = p * cos(xn + xx); 257 124 return(p); -
sasmodels/models/lib/j1_cephes.c
rfad5dc1 r95ce773 40 40 */ 41 41 42 constant double RP[8] = {43 -8.99971225705559398224E8,44 4.52228297998194034323E11,45 -7.27494245221818276015E13,46 3.68295732863852883286E15,47 0.0,48 0.0,49 0.0,50 0.0 };51 52 constant double RQ[8] = {53 6.20836478118054335476E2,54 2.56987256757748830383E5,55 8.35146791431949253037E7,56 2.21511595479792499675E10,57 4.74914122079991414898E12,58 7.84369607876235854894E14,59 8.95222336184627338078E16,60 5.32278620332680085395E1861 };62 63 constant double PP[8] = {64 7.62125616208173112003E-4,65 7.31397056940917570436E-2,66 1.12719608129684925192E0,67 5.11207951146807644818E0,68 8.42404590141772420927E0,69 5.21451598682361504063E0,70 1.00000000000000000254E0,71 0.0} ;72 73 74 constant double PQ[8] = {75 5.71323128072548699714E-4,76 6.88455908754495404082E-2,77 1.10514232634061696926E0,78 5.07386386128601488557E0,79 8.39985554327604159757E0,80 5.20982848682361821619E0,81 9.99999999999999997461E-1,82 0.0 };83 84 constant double QP[8] = {85 5.10862594750176621635E-2,86 4.98213872951233449420E0,87 7.58238284132545283818E1,88 3.66779609360150777800E2,89 7.10856304998926107277E2,90 5.97489612400613639965E2,91 2.11688757100572135698E2,92 2.52070205858023719784E1 };93 94 constant double QQ[8] = {95 7.42373277035675149943E1,96 1.05644886038262816351E3,97 4.98641058337653607651E3,98 9.56231892404756170795E3,99 7.99704160447350683650E3,100 2.82619278517639096600E3,101 3.36093607810698293419E2,102 0.0 };103 104 constant double JP[8] = {105 -4.878788132172128E-009,106 6.009061827883699E-007,107 -4.541343896997497E-005,108 1.937383947804541E-003,109 -3.405537384615824E-002,110 0.0,111 0.0,112 0.0113 };114 115 constant double MO1[8] = {116 6.913942741265801E-002,117 -2.284801500053359E-001,118 3.138238455499697E-001,119 -2.102302420403875E-001,120 5.435364690523026E-003,121 1.493389585089498E-001,122 4.976029650847191E-006,123 7.978845453073848E-001124 };125 126 constant double PH1[8] = {127 -4.497014141919556E+001,128 5.073465654089319E+001,129 -2.485774108720340E+001,130 7.222973196770240E+000,131 -1.544842782180211E+000,132 3.503787691653334E-001,133 -1.637986776941202E-001,134 3.749989509080821E-001135 };136 137 42 double J1(double x) { 138 43 … … 154 59 { 155 60 z = x * x; 156 w = polevl( z, RP , 3 ) / p1evl( z, RQ, 8 );61 w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 157 62 w = w * x * (z - Z1) * (z - Z2); 158 63 return( w ); … … 162 67 z = w * w; 163 68 164 p = polevl( z, PP , 6)/polevl( z, PQ, 6 );165 q = polevl( z, QP , 7)/p1evl( z, QQ, 7 );69 p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 70 q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 166 71 167 72 xn = x - THPIO4; … … 189 94 { 190 95 z = xx * xx; 191 p = (z-Z1) * xx * polevl( z, JP , 4 );96 p = (z-Z1) * xx * polevl( z, JPJ1, 4 ); 192 97 return( p ); 193 98 } … … 196 101 w = sqrt(q); 197 102 198 p = w * polevl( q, MO1 , 7);103 p = w * polevl( q, MO1J1, 7); 199 104 w = q*q; 200 xn = q * polevl( w, PH1 , 7) - THPIO4F;105 xn = q * polevl( w, PH1J1, 7) - THPIO4F; 201 106 p = p * cos(xn + xx); 202 107 … … 205 110 } 206 111 207 208 /*double J1c(double x) {209 return 210 } */112 //Finally J1c function that equals 2*J1(x)/x 113 double J1c(double x) { 114 return (x != 0.0 ) ? 2.0*J1(x)/x : 1.0; 115 } -
sasmodels/models/lib/polevl.c
rfad5dc1 r95ce773 51 51 */ 52 52 53 constant double RPJ1[8] = { 54 -8.99971225705559398224E8, 55 4.52228297998194034323E11, 56 -7.27494245221818276015E13, 57 3.68295732863852883286E15, 58 0.0, 59 0.0, 60 0.0, 61 0.0 }; 62 63 constant double RQJ1[8] = { 64 6.20836478118054335476E2, 65 2.56987256757748830383E5, 66 8.35146791431949253037E7, 67 2.21511595479792499675E10, 68 4.74914122079991414898E12, 69 7.84369607876235854894E14, 70 8.95222336184627338078E16, 71 5.32278620332680085395E18 72 }; 73 74 constant double PPJ1[8] = { 75 7.62125616208173112003E-4, 76 7.31397056940917570436E-2, 77 1.12719608129684925192E0, 78 5.11207951146807644818E0, 79 8.42404590141772420927E0, 80 5.21451598682361504063E0, 81 1.00000000000000000254E0, 82 0.0} ; 83 84 85 constant double PQJ1[8] = { 86 5.71323128072548699714E-4, 87 6.88455908754495404082E-2, 88 1.10514232634061696926E0, 89 5.07386386128601488557E0, 90 8.39985554327604159757E0, 91 5.20982848682361821619E0, 92 9.99999999999999997461E-1, 93 0.0 }; 94 95 constant double QPJ1[8] = { 96 5.10862594750176621635E-2, 97 4.98213872951233449420E0, 98 7.58238284132545283818E1, 99 3.66779609360150777800E2, 100 7.10856304998926107277E2, 101 5.97489612400613639965E2, 102 2.11688757100572135698E2, 103 2.52070205858023719784E1 }; 104 105 constant double QQJ1[8] = { 106 7.42373277035675149943E1, 107 1.05644886038262816351E3, 108 4.98641058337653607651E3, 109 9.56231892404756170795E3, 110 7.99704160447350683650E3, 111 2.82619278517639096600E3, 112 3.36093607810698293419E2, 113 0.0 }; 114 115 constant double JPJ1[8] = { 116 -4.878788132172128E-009, 117 6.009061827883699E-007, 118 -4.541343896997497E-005, 119 1.937383947804541E-003, 120 -3.405537384615824E-002, 121 0.0, 122 0.0, 123 0.0 124 }; 125 126 constant double MO1J1[8] = { 127 6.913942741265801E-002, 128 -2.284801500053359E-001, 129 3.138238455499697E-001, 130 -2.102302420403875E-001, 131 5.435364690523026E-003, 132 1.493389585089498E-001, 133 4.976029650847191E-006, 134 7.978845453073848E-001 135 }; 136 137 constant double PH1J1[8] = { 138 -4.497014141919556E+001, 139 5.073465654089319E+001, 140 -2.485774108720340E+001, 141 7.222973196770240E+000, 142 -1.544842782180211E+000, 143 3.503787691653334E-001, 144 -1.637986776941202E-001, 145 3.749989509080821E-001 146 }; 147 148 constant double PPJ0[8] = { 149 7.96936729297347051624E-4, 150 8.28352392107440799803E-2, 151 1.23953371646414299388E0, 152 5.44725003058768775090E0, 153 8.74716500199817011941E0, 154 5.30324038235394892183E0, 155 9.99999999999999997821E-1, 156 0.0 157 }; 158 159 constant double PQJ0[8] = { 160 9.24408810558863637013E-4, 161 8.56288474354474431428E-2, 162 1.25352743901058953537E0, 163 5.47097740330417105182E0, 164 8.76190883237069594232E0, 165 5.30605288235394617618E0, 166 1.00000000000000000218E0, 167 0.0 168 }; 169 170 constant double QPJ0[8] = { 171 -1.13663838898469149931E-2, 172 -1.28252718670509318512E0, 173 -1.95539544257735972385E1, 174 -9.32060152123768231369E1, 175 -1.77681167980488050595E2, 176 -1.47077505154951170175E2, 177 -5.14105326766599330220E1, 178 -6.05014350600728481186E0, 179 }; 180 181 constant double QQJ0[8] = { 182 /* 1.00000000000000000000E0,*/ 183 6.43178256118178023184E1, 184 8.56430025976980587198E2, 185 3.88240183605401609683E3, 186 7.24046774195652478189E3, 187 5.93072701187316984827E3, 188 2.06209331660327847417E3, 189 2.42005740240291393179E2, 190 }; 191 192 constant double YPJ0[8] = { 193 1.55924367855235737965E4, 194 -1.46639295903971606143E7, 195 5.43526477051876500413E9, 196 -9.82136065717911466409E11, 197 8.75906394395366999549E13, 198 -3.46628303384729719441E15, 199 4.42733268572569800351E16, 200 -1.84950800436986690637E16, 201 }; 202 203 204 constant double YQJ0[7] = { 205 /* 1.00000000000000000000E0,*/ 206 1.04128353664259848412E3, 207 6.26107330137134956842E5, 208 2.68919633393814121987E8, 209 8.64002487103935000337E10, 210 2.02979612750105546709E13, 211 3.17157752842975028269E15, 212 2.50596256172653059228E17, 213 }; 214 215 constant double RPJ0[8] = { 216 -4.79443220978201773821E9, 217 1.95617491946556577543E12, 218 -2.49248344360967716204E14, 219 9.70862251047306323952E15, 220 0.0, 221 0.0, 222 0.0, 223 0.0 224 }; 225 226 constant double RQJ0[8] = { 227 /* 1.00000000000000000000E0,*/ 228 4.99563147152651017219E2, 229 1.73785401676374683123E5, 230 4.84409658339962045305E7, 231 1.11855537045356834862E10, 232 2.11277520115489217587E12, 233 3.10518229857422583814E14, 234 3.18121955943204943306E16, 235 1.71086294081043136091E18, 236 }; 237 238 constant double MOJ0[8] = { 239 -6.838999669318810E-002, 240 1.864949361379502E-001, 241 -2.145007480346739E-001, 242 1.197549369473540E-001, 243 -3.560281861530129E-003, 244 -4.969382655296620E-002, 245 -3.355424622293709E-006, 246 7.978845717621440E-001 247 }; 248 249 constant double PHJ0[8] = { 250 3.242077816988247E+001, 251 -3.630592630518434E+001, 252 1.756221482109099E+001, 253 -4.974978466280903E+000, 254 1.001973420681837E+000, 255 -1.939906941791308E-001, 256 6.490598792654666E-002, 257 -1.249992184872738E-001 258 }; 259 260 constant double JPJ0[8] = { 261 -6.068350350393235E-008, 262 6.388945720783375E-006, 263 -3.969646342510940E-004, 264 1.332913422519003E-002, 265 -1.729150680240724E-001, 266 0.0, 267 0.0, 268 0.0 269 }; 53 270 54 271 double polevl( double x, constant double *coef, int N ) { -
doc/genmodel.py
r044a9c1 r824465b 29 29 'q_max' : 1.0, 30 30 'nq' : 1000, 31 'nq2d' : 100 ,31 'nq2d' : 1000, 32 32 'vmin' : 1e-3, # floor for the 2D data results 33 33 'qx_max' : 0.5, 34 #'colormap' : 'gist_ncar', 35 'colormap' : 'jet', 34 36 } 35 37 … … 61 63 Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf)) 62 64 ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='lower', 63 extent=[-qx_max, qx_max, -qx_max, qx_max], cmap= pylab.cm.jet)65 extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=opts['colormap']) 64 66 ax.set_xlabel(r'$Q_x \/(\AA^{-1})$') 65 67 ax.set_ylabel(r'$Q_y \/(\AA^{-1})$') 66 67 # im = self.subplot.imshow(output, interpolation='nearest',68 # origin='lower',69 # vmin=zmin_temp, vmax=self.zmax_2D,70 # cmap=self.cmap,71 # extent=(self.xmin_2D, self.xmax_2D,72 # self.ymin_2D, self.ymax_2D))73 74 def ice_cm():75 from matplotlib._cm import _Blues_data76 from matplotlib import colors77 from matplotlib import rcParams78 def from_white(segments):79 scale = 1.0/segments[0][1]80 return [(k, v*scale, w*scale) for k, v, w in segments]81 ice_data = dict((k,from_white(v)) for k,v in _Blues_data.items())82 ice = colors.LinearSegmentedColormap("ice", ice_data, rcParams['image.lut'])83 return ice84 85 68 86 69 # Generate image -
sasmodels/compare.py
rd6850fa raf92b73 539 539 try: 540 540 base_value, base_time = time_calculation(base, pars, Nbase) 541 print("%s t=%. 1f ms, intensity=%.0f"541 print("%s t=%.2f ms, intensity=%.0f" 542 542 % (base.engine, base_time, sum(base_value))) 543 543 except ImportError: … … 550 550 try: 551 551 comp_value, comp_time = time_calculation(comp, pars, Ncomp) 552 print("%s t=%. 1f ms, intensity=%.0f"552 print("%s t=%.2f ms, intensity=%.0f" 553 553 % (comp.engine, comp_time, sum(comp_value))) 554 554 except ImportError: … … 584 584 if Ncomp > 0: plt.subplot(131) 585 585 plot_theory(data, base_value, view=view, use_data=False, limits=limits) 586 plt.title("%s t=%. 1f ms"%(base.engine, base_time))586 plt.title("%s t=%.2f ms"%(base.engine, base_time)) 587 587 #cbar_title = "log I" 588 588 if Ncomp > 0: 589 589 if Nbase > 0: plt.subplot(132) 590 590 plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 591 plt.title("%s t=%. 1f ms"%(comp.engine, comp_time))591 plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 592 592 #cbar_title = "log I" 593 593 if Ncomp > 0 and Nbase > 0: -
sasmodels/data.py
rc094758 r092cb3c 389 389 if view is 'log': 390 390 mdata[mdata <= 0] = masked 391 plt.errorbar(data.x /10, scale*mdata, yerr=data.dy, fmt='.')391 plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.') 392 392 all_positive = all_positive and (mdata > 0).all() 393 393 some_present = some_present or (mdata.count() > 0) … … 412 412 if view == 'q4' or not some_present or not all_positive 413 413 else view) 414 plt.xlabel("$q$/ nm$^{-1}$")414 plt.xlabel("$q$/A$^{-1}$") 415 415 plt.ylabel('$I(q)$') 416 416 … … 422 422 if num_plots > 1: 423 423 plt.subplot(1, num_plots, (use_data or use_theory) + 1) 424 plt.plot(data.x /10, mresid, '-')425 plt.xlabel("$q$/ nm$^{-1}$")424 plt.plot(data.x, mresid, '-') 425 plt.xlabel("$q$/A$^{-1}$") 426 426 plt.ylabel('residuals') 427 427 plt.xscale('linear' if not some_present else view)
Note: See TracChangeset
for help on using the changeset viewer.