1 | #!/usr/bin/env python |
---|
2 | """ |
---|
3 | Class to validate a given 2D model w/ 3 axes by averaging it |
---|
4 | and comparing to 1D prediction. |
---|
5 | |
---|
6 | The equation used for averaging is: |
---|
7 | |
---|
8 | (integral dphi from 0 to 2pi)(integral dtheta from 0 to pi) |
---|
9 | p(theta, phi) I(q) sin(theta) dtheta |
---|
10 | |
---|
11 | = (1/N_phi) (1/N_theta) (pi/2) (sum over N_phi) (sum over N_theta) |
---|
12 | p(theta_i, phi_i) I(q) sin(theta_i) |
---|
13 | |
---|
14 | where p(theta, phi) is the probability distribution normalized to 4pi. |
---|
15 | In the current case, we put p(theta, phi) = 1. |
---|
16 | |
---|
17 | The normalization factor results from: |
---|
18 | 2pi/N_phi for the phi sum |
---|
19 | x pi/N_theta for the theta sum |
---|
20 | x 1/(4pi) because p is normalized to 4pi |
---|
21 | -------------- |
---|
22 | = (1/N_phi) (1/N_theta) (pi/2) |
---|
23 | |
---|
24 | Note: |
---|
25 | Averaging the 3-axes 2D scattering intensity give a slightly |
---|
26 | different output than the 1D function |
---|
27 | at hight Q (Q>~0.2). This is due(?) to the way the IGOR library |
---|
28 | averages, taking only 76 points in alpha, the angle between |
---|
29 | the axis of the ellipsoid and the q vector. |
---|
30 | |
---|
31 | Note: Core-shell and sphere models are symmetric around |
---|
32 | all axes and don't need to be tested in the following way. |
---|
33 | """ |
---|
34 | import sys, math |
---|
35 | from sans.models.EllipticalCylinderModel import EllipticalCylinderModel |
---|
36 | from sans.models.ParallelepipedModel import ParallelepipedModel |
---|
37 | from sans.models.TriaxialEllipsoidModel import TriaxialEllipsoidModel |
---|
38 | |
---|
39 | class Validate2D: |
---|
40 | """ |
---|
41 | Class to validate a given 2D model by averaging it |
---|
42 | and comparing to 1D prediction. |
---|
43 | """ |
---|
44 | |
---|
45 | def __init__(self): |
---|
46 | """ Initialization """ |
---|
47 | # Precision for the result comparison |
---|
48 | self.precision = 0.000001 |
---|
49 | # Flag for end result |
---|
50 | self.passed = True |
---|
51 | |
---|
52 | def __call__(self, model_class=EllipticalCylinderModel, points = 101): |
---|
53 | """ |
---|
54 | Perform test and produce output file |
---|
55 | @param model_class: python class of the model to test |
---|
56 | """ |
---|
57 | print "Averaging %s" % model_class.__name__ |
---|
58 | passed = True |
---|
59 | |
---|
60 | npts =points |
---|
61 | model = model_class() |
---|
62 | #model.setParam('scale', 1.0) |
---|
63 | #model.setParam('contrast', 1.0) |
---|
64 | |
---|
65 | theta_label = 'cyl_theta' |
---|
66 | if not model.params.has_key(theta_label): |
---|
67 | theta_label = 'parallel_theta' |
---|
68 | if not model.params.has_key(theta_label): |
---|
69 | theta_label = 'axis_theta' |
---|
70 | |
---|
71 | phi_label = 'cyl_phi' |
---|
72 | if not model.params.has_key(phi_label): |
---|
73 | phi_label = 'parallel_phi' |
---|
74 | if not model.params.has_key(phi_label): |
---|
75 | phi_label = 'axis_phi' |
---|
76 | |
---|
77 | psi_label = 'cyl_psi' |
---|
78 | if not model.params.has_key(psi_label): |
---|
79 | psi_label = 'parallel_psi' |
---|
80 | if not model.params.has_key(psi_label): |
---|
81 | psi_label = 'axis_psi' |
---|
82 | |
---|
83 | output_f = open("%s_avg.txt" % model.__class__.__name__,'w') |
---|
84 | output_f.write("<q_average> <2d_average> <1d_average>\n") |
---|
85 | |
---|
86 | for i_q in range(1, 40): |
---|
87 | q = 0.01*i_q |
---|
88 | sum = 0.0 |
---|
89 | for i_theta in range(npts): |
---|
90 | theta = math.pi/npts*(i_theta+1) |
---|
91 | |
---|
92 | model.setParam(theta_label, theta) |
---|
93 | |
---|
94 | for j in range(npts): |
---|
95 | model.setParam(phi_label, math.pi * 2.0 / npts * j) |
---|
96 | for k in range(npts): |
---|
97 | model.setParam(psi_label, math.pi * 2.0 / npts * k) |
---|
98 | if str(model.run([q, 0])).count("INF")>0: |
---|
99 | print "ERROR", q, theta, math.pi * 2.0 / npts * k |
---|
100 | sum += model.run([q, 0])*math.sin(theta) |
---|
101 | #sum += model.run([q, 0]) |
---|
102 | value = sum/npts/(npts)/npts*math.pi/2.0 |
---|
103 | ana = model.run(q) |
---|
104 | if q<0.2 and (value-ana)/ana>0.05: |
---|
105 | passed = False |
---|
106 | output_f.write("%10g %10g %10g\n" % (q, value, ana)) |
---|
107 | print "Q=%g: %10g %10g %10g %10g" % (q, value, ana, value-ana, value/ana) |
---|
108 | |
---|
109 | output_f.close() |
---|
110 | return passed |
---|
111 | |
---|
112 | |
---|
113 | |
---|
114 | if __name__ == '__main__': |
---|
115 | validator = Validate2D() |
---|
116 | |
---|
117 | #te was not passed. |
---|
118 | te_passed =validator(TriaxialEllipsoidModel, points=76) |
---|
119 | pp_passed = validator(ParallelepipedModel, points=76) |
---|
120 | ell_passed = validator(EllipticalCylinderModel, points=76) |
---|
121 | |
---|
122 | print "" |
---|
123 | print "Model Passed" |
---|
124 | print "TriaxialEllipsoid %s" % te_passed |
---|
125 | print "ParallelepipedModel %s" % pp_passed |
---|
126 | print "EllipticalCylinder %s" % ell_passed |
---|
127 | |
---|
128 | |
---|
129 | |
---|
130 | |
---|
131 | |
---|