1 | |
---|
2 | import VolumeCanvas |
---|
3 | from sans.models.SphereModel import SphereModel |
---|
4 | from sans.models.CoreShellModel import CoreShellModel |
---|
5 | |
---|
6 | import math, time |
---|
7 | |
---|
8 | def form_factor(q, r): |
---|
9 | qr = q*r |
---|
10 | f = 3*( math.sin(qr) - qr*math.cos(qr) ) / (qr*qr*qr) |
---|
11 | return f*f |
---|
12 | |
---|
13 | def test_1(): |
---|
14 | |
---|
15 | radius = 15 |
---|
16 | |
---|
17 | density = .1 |
---|
18 | vol = 4/3*math.pi*radius*radius*radius |
---|
19 | npts = vol*density |
---|
20 | |
---|
21 | canvas = VolumeCanvas.VolumeCanvas() |
---|
22 | canvas.setParam('lores_density', density) |
---|
23 | handle = canvas.add('sphere') |
---|
24 | canvas.setParam('%s.radius' % handle, radius) |
---|
25 | canvas.setParam('%s.contrast' % handle, 1.0) |
---|
26 | |
---|
27 | |
---|
28 | if False: |
---|
29 | # Time test |
---|
30 | t_0 = time.time() |
---|
31 | value_1 = 1.0e8*canvas.getIq(0.1) |
---|
32 | print "density = 0.1: output=%g time=%g" % (value_1, time.time()-t_0) |
---|
33 | |
---|
34 | t_0 = time.time() |
---|
35 | canvas.setParam('lores_density', 1) |
---|
36 | value_1 = 1.0e8*canvas.getIq(0.1) |
---|
37 | print "density = 1000: output=%g time=%g" % (value_1, time.time()-t_0) |
---|
38 | |
---|
39 | t_0 = time.time() |
---|
40 | canvas.setParam('lores_density', 0.01) |
---|
41 | value_1 = 1.0e8*canvas.getIq(0.1) |
---|
42 | print "density = 0.00001: output=%g time=%g" % (value_1, time.time()-t_0) |
---|
43 | print |
---|
44 | |
---|
45 | |
---|
46 | sphere = SphereModel() |
---|
47 | sphere.setParam('radius', radius) |
---|
48 | sphere.setParam('scale', 1.0) |
---|
49 | sphere.setParam('contrast', 1.0) |
---|
50 | |
---|
51 | |
---|
52 | # Simple sphere sum(Pr) = (rho*V)^2 |
---|
53 | # each p(r) point has a volume of 1/density |
---|
54 | |
---|
55 | for i in range(35): |
---|
56 | q = 0.001 + 0.01*i |
---|
57 | |
---|
58 | |
---|
59 | |
---|
60 | #sim_1 = 1.0e8*canvas.getIq(q)*4/3*math.pi/(density*density*density) |
---|
61 | sim_1 = canvas.getIq(q) |
---|
62 | ana_1 = sphere.run(q) |
---|
63 | #ana_1 = form_factor(q, radius) |
---|
64 | |
---|
65 | print "q=%g sim=%g ana=%g ratio=%g" % (q, sim_1, ana_1, sim_1/ana_1) |
---|
66 | |
---|
67 | def test_2(): |
---|
68 | radius = 15.0 |
---|
69 | thickness = 5.0 |
---|
70 | |
---|
71 | core_vol = 4.0/3.0*math.pi*radius*radius*radius |
---|
72 | outer_radius = radius+thickness |
---|
73 | shell_vol = 4.0/3.0*math.pi*outer_radius*outer_radius*outer_radius - core_vol |
---|
74 | shell_sld = -1.0*core_vol/shell_vol |
---|
75 | print "Shell SLD", shell_sld |
---|
76 | |
---|
77 | |
---|
78 | density = .1 |
---|
79 | vol = 4/3*math.pi*radius*radius*radius |
---|
80 | npts = vol*density |
---|
81 | |
---|
82 | canvas = VolumeCanvas.VolumeCanvas() |
---|
83 | canvas.setParam('lores_density', density) |
---|
84 | handle = canvas.add('sphere') |
---|
85 | canvas.setParam('%s.radius' % handle, outer_radius) |
---|
86 | canvas.setParam('%s.contrast' % handle, shell_sld) |
---|
87 | |
---|
88 | handle2 = canvas.add('sphere') |
---|
89 | canvas.setParam('%s.radius' % handle2, radius) |
---|
90 | canvas.setParam('%s.contrast' % handle2, 1.0) |
---|
91 | |
---|
92 | |
---|
93 | |
---|
94 | # Core-shell |
---|
95 | sphere = CoreShellModel() |
---|
96 | # Core radius |
---|
97 | sphere.setParam('radius', radius) |
---|
98 | # Shell thickness |
---|
99 | sphere.setParam('thickness', thickness) |
---|
100 | sphere.setParam('core_sld', 1.0) |
---|
101 | |
---|
102 | |
---|
103 | sphere.setParam('shell_sld', shell_sld) |
---|
104 | sphere.setParam('solvent_sld',0.0) |
---|
105 | sphere.setParam('background',0.0) |
---|
106 | sphere.setParam('scale',1.0) |
---|
107 | |
---|
108 | out = open("lores_test.txt",'w') |
---|
109 | out.write("<q> <sim> <ana>\n") |
---|
110 | |
---|
111 | for i in range(65): |
---|
112 | q = 0.001 + 0.01*i |
---|
113 | |
---|
114 | # For each volume integral that we change to a sum, |
---|
115 | # we must multiply by 1/density = V/N |
---|
116 | # Since we want P(r)/V, we will need to multiply |
---|
117 | # the sum by 1/(N*density), where N is the number of |
---|
118 | # points without overlap. Since we already divide |
---|
119 | # by N when calculating I(q), we only need to divide |
---|
120 | # by the density here. We divide by N in the |
---|
121 | # calculation because it is difficult to estimate it here. |
---|
122 | |
---|
123 | |
---|
124 | # Put the factor 2 in the simulation two... |
---|
125 | sim_1 = canvas.getIq(q) |
---|
126 | ana_1 = sphere.run(q) |
---|
127 | |
---|
128 | print "q=%g sim=%g ana=%g ratio=%g" % (q, sim_1, ana_1, sim_1/ana_1) |
---|
129 | out.write( "%g %g %g\n" % (q, sim_1, ana_1)) |
---|
130 | |
---|
131 | out.close() |
---|
132 | |
---|
133 | def test_4(): |
---|
134 | radius = 15 |
---|
135 | |
---|
136 | density = .1 |
---|
137 | vol = 4/3*math.pi*radius*radius*radius |
---|
138 | npts = vol*density |
---|
139 | |
---|
140 | |
---|
141 | canvas = VolumeCanvas.VolumeCanvas() |
---|
142 | canvas.setParam('lores_density', density) |
---|
143 | #handle = canvas.add('sphere') |
---|
144 | #canvas.setParam('%s.radius' % handle, radius) |
---|
145 | #canvas.setParam('%s.contrast' % handle, 1.0) |
---|
146 | |
---|
147 | pdb = canvas.add('test.pdb') |
---|
148 | |
---|
149 | |
---|
150 | |
---|
151 | sphere = SphereModel() |
---|
152 | sphere.setParam('radius', radius) |
---|
153 | sphere.setParam('scale', 1.0) |
---|
154 | sphere.setParam('contrast', 1.0) |
---|
155 | |
---|
156 | |
---|
157 | # Simple sphere sum(Pr) = (rho*V)^2 |
---|
158 | # each p(r) point has a volume of 1/density |
---|
159 | |
---|
160 | for i in range(35): |
---|
161 | q = 0.001 + 0.01*i |
---|
162 | |
---|
163 | |
---|
164 | |
---|
165 | #sim_1 = 1.0e8*canvas.getIq(q)*4/3*math.pi/(density*density*density) |
---|
166 | sim_1 = canvas.getIq(q) |
---|
167 | ana_1 = sphere.run(q) |
---|
168 | #ana_1 = form_factor(q, radius) |
---|
169 | |
---|
170 | print "q=%g sim=%g ana=%g ratio=%g" % (q, sim_1, ana_1, sim_1/ana_1) |
---|
171 | |
---|
172 | def test_5(): |
---|
173 | from sans.models.SphereModel import SphereModel |
---|
174 | model = VolumeCanvas.VolumeCanvas() |
---|
175 | |
---|
176 | handle = model.add('sphere') |
---|
177 | |
---|
178 | radius = 10 |
---|
179 | density = .1 |
---|
180 | |
---|
181 | ana = SphereModel() |
---|
182 | ana.setParam('scale', 1.0) |
---|
183 | ana.setParam('contrast', 1.0) |
---|
184 | ana.setParam('background', 0.0) |
---|
185 | ana.setParam('radius', radius) |
---|
186 | |
---|
187 | model.setParam('lores_density', density) |
---|
188 | model.setParam('%s.radius' % handle, radius) |
---|
189 | model.setParam('scale' , 1.0) |
---|
190 | model.setParam('%s.contrast' % handle, 1.0) |
---|
191 | model.setParam('background' , 0.0) |
---|
192 | |
---|
193 | ana = ana.runXY([0.1, 0.1]) |
---|
194 | sim = model.getIq2D(0.1, 0.1) |
---|
195 | print ana, sim, sim/ana, ana/sim |
---|
196 | |
---|
197 | def test_6(): |
---|
198 | from sans.models.CylinderModel import CylinderModel |
---|
199 | radius = 5 |
---|
200 | length = 40 |
---|
201 | density = 20 |
---|
202 | |
---|
203 | ana = CylinderModel() |
---|
204 | ana.setParam('scale', 1.0) |
---|
205 | ana.setParam('contrast', 1.0) |
---|
206 | ana.setParam('background', 0.0) |
---|
207 | ana.setParam('radius', radius) |
---|
208 | ana.setParam('length', length) |
---|
209 | |
---|
210 | # Along Y |
---|
211 | ana.setParam('cyl_theta', 1.57) |
---|
212 | ana.setParam('cyl_phi', 1.57) |
---|
213 | |
---|
214 | # Along Z |
---|
215 | #ana.setParam('cyl_theta', 0) |
---|
216 | #ana.setParam('cyl_phi', 0) |
---|
217 | |
---|
218 | model = VolumeCanvas.VolumeCanvas() |
---|
219 | handle = model.add('cylinder') |
---|
220 | model.setParam('lores_density', density) |
---|
221 | model.setParam('%s.radius' % handle, radius) |
---|
222 | model.setParam('%s.length' % handle, length) |
---|
223 | model.setParam('scale' , 1.0) |
---|
224 | model.setParam('%s.contrast' % handle, 1.0) |
---|
225 | model.setParam('background' , 0.0) |
---|
226 | |
---|
227 | # Along Y |
---|
228 | model.setParam('%s.orientation' % handle, [0,0,0]) |
---|
229 | |
---|
230 | # Along Z |
---|
231 | #model.setParam('%s.orientation' % handle, [1.57,0,0]) |
---|
232 | |
---|
233 | |
---|
234 | print model.npts |
---|
235 | for i in range(40): |
---|
236 | qmax = 0.5 |
---|
237 | anaX = ana.runXY([qmax*i/40.0, 0.0]) |
---|
238 | simX = model.getIq2D(qmax*i/40.0, 0.0) |
---|
239 | |
---|
240 | anaY = ana.runXY([0, qmax*i/40.0]) |
---|
241 | simY = model.getIq2D(0, qmax*i/40.0) |
---|
242 | print anaX, simX, simX/anaX, '|', anaY, simY, simY/anaY |
---|
243 | |
---|
244 | def test_7(): |
---|
245 | from sans.models.CoreShellModel import CoreShellModel |
---|
246 | print "Testing core-shell" |
---|
247 | radius = 15 |
---|
248 | thickness = 5 |
---|
249 | density = 5 |
---|
250 | |
---|
251 | core_vol = 4.0/3.0*math.pi*radius*radius*radius |
---|
252 | outer_radius = radius+thickness |
---|
253 | shell_vol = 4.0/3.0*math.pi*outer_radius*outer_radius*outer_radius - core_vol |
---|
254 | shell_sld = -1.0*core_vol/shell_vol |
---|
255 | |
---|
256 | # Core-shell |
---|
257 | sphere = CoreShellModel() |
---|
258 | # Core radius |
---|
259 | sphere.setParam('radius', radius) |
---|
260 | # Shell thickness |
---|
261 | sphere.setParam('thickness', thickness) |
---|
262 | sphere.setParam('core_sld', 1.0) |
---|
263 | sphere.setParam('shell_sld', shell_sld) |
---|
264 | sphere.setParam('solvent_sld', 0.0) |
---|
265 | sphere.setParam('background', 0.0) |
---|
266 | sphere.setParam('scale', 1.0) |
---|
267 | ana = sphere |
---|
268 | |
---|
269 | canvas = VolumeCanvas.VolumeCanvas() |
---|
270 | canvas.setParam('lores_density', density) |
---|
271 | |
---|
272 | handle = canvas.add('sphere') |
---|
273 | canvas.setParam('%s.radius' % handle, outer_radius) |
---|
274 | canvas.setParam('%s.contrast' % handle, shell_sld) |
---|
275 | |
---|
276 | handle2 = canvas.add('sphere') |
---|
277 | canvas.setParam('%s.radius' % handle2, radius) |
---|
278 | canvas.setParam('%s.contrast' % handle2, 1.0) |
---|
279 | |
---|
280 | canvas.setParam('scale' , 1.0) |
---|
281 | canvas.setParam('background' , 0.0) |
---|
282 | |
---|
283 | |
---|
284 | """ Testing default core-shell orientation """ |
---|
285 | qlist = [.0001, 0.002, .01, .1, 1.0, 5.] |
---|
286 | for q in qlist: |
---|
287 | ana_val = ana.runXY([q, 0.2]) |
---|
288 | sim_val, err = canvas.getIq2DError(q, 0.2) |
---|
289 | print ana_val, sim_val, sim_val/ana_val, err, (sim_val-ana_val)/err |
---|
290 | |
---|
291 | |
---|
292 | |
---|
293 | if __name__ == "__main__": |
---|
294 | test_6() |
---|