[27953d1] | 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.1). 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_avg1.txt" % model.__class__.__name__,'w') |
---|
| 84 | output_f.write("<q_average> <2d_average> <1d_average>\n") |
---|
| 85 | |
---|
| 86 | for i_q in range(1, 30): |
---|
| 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_passed =validator(TriaxialEllipsoidModel, points=76) |
---|
| 118 | pp_passed = validator(ParallelepipedModel, points=76) |
---|
| 119 | ell_passed = validator(EllipticalCylinderModel, points=76) |
---|
| 120 | |
---|
| 121 | print "" |
---|
| 122 | print "Model Passed" |
---|
| 123 | print "TriaxialEllipsoid %s" % te_passed |
---|
| 124 | print "EllipticalCylinder %s" % ell_passed |
---|
| 125 | print "ParallelepipedModel %s" % pp_passed |
---|
| 126 | |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | |
---|
| 130 | |
---|