diff --git a/common/lib/capa/capa/capa_problem.py b/common/lib/capa/capa/capa_problem.py index 451891d067..db42fb698a 100644 --- a/common/lib/capa/capa/capa_problem.py +++ b/common/lib/capa/capa/capa_problem.py @@ -33,6 +33,7 @@ from xml.sax.saxutils import unescape import chem import chem.chemcalc import chem.chemtools +import chem.miller import calc from correctmap import CorrectMap @@ -67,7 +68,8 @@ global_context = {'random': random, 'calc': calc, 'eia': eia, 'chemcalc': chem.chemcalc, - 'chemtools': chem.chemtools} + 'chemtools': chem.chemtools, + 'miller': chem.miller} # These should be removed from HTML output, including all subelements html_problem_semantics = ["codeparam", "responseparam", "answer", "script", "hintgroup"] diff --git a/common/lib/capa/capa/chem/miller.py b/common/lib/capa/capa/chem/miller.py new file mode 100644 index 0000000000..38b2ea88e2 --- /dev/null +++ b/common/lib/capa/capa/chem/miller.py @@ -0,0 +1,193 @@ +# 1) Calculate miller indices by points coordinates +# 2) Grader for miller indeces and lattice type + + +import numpy as np +import math +import fractions as fr +import decimal +import unittest + + +def lcm(a, b): + """ return lcm of a,b """ + return a * b / fr.gcd(a, b) + + +def section_to_fraction(distance): + """ Convert float distance, that plane cut on axis + to fraction. Return inverted fraction + + """ + if np.isnan(distance): # plane || to axis (or contains axis) + print distance, 0 + # return inverted fration to a == nan == 1/0 => 0 / 1 + return fr.Fraction(0, 1) + elif distance == 0: # plane goes through origin + return fr.Fraction(1, 1) # ERROR, need shift of coordinates + else: + # limit_denominator to closest nicest fraction + fract = fr.Fraction(distance).limit_denominator() # 5 / 2 : numerator / denominator + print 'Distance', distance, 'Inverted fraction', fract + # return inverted fraction + return fr.Fraction(fract.denominator, fract.numerator) + + +def sub_miller(sections): + ''' Calculate miller indices. + Plane does not intersect origin + ''' + fracts = [section_to_fraction(section) for section in sections] + + print sections, fracts + + common_denominator = reduce(lcm, [fract.denominator for fract in fracts]) + print 'common_denominator:', common_denominator + + # 2) lead to a common denominator + # 3) throw away denominator + miller = [fract.numerator * math.fabs(common_denominator) / fract.denominator for fract in fracts] + + # import ipdb; ipdb.set_trace() + # nice output: + output = '(' + ''.join(map(str, map(decimal.Decimal, miller))) + ')' + print 'Miller indices:', output + return output + + +def miller(points): + """Calculate miller indices of plane + """ + + print "\nCalculating miller indices:" + print 'Points:\n', points + # calculate normal to plane + N = np.cross(points[1] - points[0], points[2] - points[0]) + print "Normal:", N + + # origin + O = np.array([0, 0, 0]) + + # point of plane + P = points[0] + + # equation of a line for axes: O + (B-O) * t + # t - parameters, B = [Bx, By, Bz]: + B = map(np.array, [[1.0, 0, 0], [0, 1.0, 0], [0, 0, 1.0]]) + + # coordinates of intersections with axis: + sections = [np.dot(P - O, N) / np.dot(B_axis, N) if np.dot(B_axis, N) != 0 else np.nan for B_axis in B] + # import ipdb; ipdb.set_trace() + + if any(x == 0 for x in sections): # Plane goes through origin. + # Need to shift plane out of origin. + # For this change origin position + + # 1) find cube vertex, not crossed by plane + vertices = [ + # top + np.array([1.0, 1.0, 1.0]), + np.array([0.0, 0.0, 1.0]), + np.array([1.0, 0.0, 1.0]), + np.array([0.0, 1.0, 1.0]), + # bottom, except 0,0,0 + np.array([1.0, 0.0, 0.0]), + np.array([0.0, 1.0, 0.0]), + np.array([1.0, 1.0, 1.0]), + ] + + for vertex in vertices: + if np.dot(vertex - O, N) != 0: # vertex in plane + new_origin = vertex + break + + # get axis with center in new origin + X = np.array([1 - new_origin[0], new_origin[1], new_origin[2]]) + Y = np.array([new_origin[0], 1 - new_origin[1], new_origin[2]]) + Z = np.array([new_origin[0], new_origin[1], 1 - new_origin[2]]) + + # 2) calculate miller indexes by new origin + new_axis = [X - new_origin, Y - new_origin, Z - new_origin] + # coordinates of intersections with axis: + sections = [np.dot(P - new_origin, N) / np.dot(B_axis, N) if np.dot(B_axis, N) != 0 else np.nan for B_axis in new_axis] + # 3) fix signs of indexes + # 0 -> 1, 1 -> -1 + sections = (1 - 2 * new_origin) * sections + return sub_miller(sections) + + +class Test_Crystallography_Grader(unittest.TestCase): + ''' test crystallography grade function ''' + + def test_1(self): + x = np.array([0.5, 0, 0]) + y = np.array([0, 0.5, 0]) + z = np.array([0, 0, 0.5]) + self.assertEqual(miller(np.array([x, y, z])), '(222)') + + def test_2(self): + x = np.array([1, 0, 0]) + y = np.array([0, 1, 0]) + z = np.array([0, 0, 1]) + self.assertEqual(miller(np.array([x, y, z])), '(111)') + + def test_3(self): + x = np.array([1, 0.5, 1]) + y = np.array([1, 1, 0.5]) + z = np.array([0.5, 1, 1]) + self.assertEqual(miller(np.array([x, y, z])), '(222)') + + def test_4(self): + x = np.array([1. / 3, 1., 0]) + y = np.array([0, 2. / 3., 0]) + z = np.array([0, 1, 1. / 3]) + self.assertEqual(miller(np.array([x, y, z])), '(-33-3)') + + def test_5(self): + x = np.array([1. / 3, 1., 0]) + y = np.array([0, 1. / 3., 0]) + z = np.array([0, 1, 1. / 3]) + self.assertEqual(miller(np.array([x, y, z])), '(-63-6)') + + def test_6(self): + x = np.array([0, 1. / 4., 0]) + y = np.array([1. / 4, 0, 0]) + z = np.array([0, 0, 1. / 4]) + self.assertEqual(miller(np.array([x, y, z])), '(444)') + + def test_7(self): # goes throug origin + x = np.array([0, 1., 0]) + y = np.array([1., 0, 0]) + z = np.array([0.5, 1., 0]) + self.assertEqual(miller(np.array([x, y, z])), '(00-1)') + + def test_8(self): + x = np.array([0, 1., 0.5]) + y = np.array([1., 0, 0.5]) + z = np.array([0.5, 1., 0.5]) + self.assertEqual(miller(np.array([x, y, z])), '(002)') + + def test_9(self): + x = np.array([0, 1. / 4., 0]) + y = np.array([1. / 4, 0, 0]) + z = np.array([0, 0, 1. / 4]) + self.assertEqual(miller(np.array([x, y, z])), '(444)') + + def test_10(self): + x = np.array([0, 1. / 4., 0]) + y = np.array([1. / 4, 0, 0]) + z = np.array([0, 0, 1. / 4]) + self.assertEqual(miller(np.array([x, y, z])), '(444)') + + +def suite(): + + testcases = [Test_Crystallography_Grader] + suites = [] + for testcase in testcases: + suites.append(unittest.TestLoader().loadTestsFromTestCase(testcase)) + return unittest.TestSuite(suites) + +if __name__ == "__main__": + unittest.TextTestRunner(verbosity=2).run(suite()) +