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