@@ -860,3 +860,134 @@ def test_square_lattice_gives_two_perpendicular_peaks(self):
860860 f"or { expected2 :.1f} ° within 5°. "
861861 "One lattice direction may be misidentified."
862862 )
863+
864+
865+ # ── FFT lattice symmetry ──────────────────────────────────────────────────────
866+
867+ class TestFftLatticeSymmetry :
868+ """FFT power spectrum symmetry for square and hexagonal atomic corrugations."""
869+
870+ @staticmethod
871+ def _power_peaks (arr : np .ndarray , n_peaks : int ) -> list [tuple [int , int ]]:
872+ """Hanning-windowed FFT power; return (row_offset, col_offset) tuples
873+ of the n_peaks strongest peaks, excluding a DC disc of radius 2 px."""
874+ Ny , Nx = arr .shape
875+ win2d = np .outer (np .hanning (Ny ), np .hanning (Nx ))
876+ power = np .abs (np .fft .fftshift (np .fft .fft2 (arr * win2d ))) ** 2
877+
878+ cy , cx = Ny // 2 , Nx // 2
879+ Y , X = np .mgrid [:Ny , :Nx ]
880+ power [(Y - cy ) ** 2 + (X - cx ) ** 2 < 4.0 ] = 0.0 # DC disc r < 2
881+
882+ peaks = []
883+ suppress = max (3 , min (Ny , Nx ) // 20 )
884+ p = power .copy ()
885+ for _ in range (n_peaks ):
886+ idx = int (np .argmax (p ))
887+ py , px = divmod (idx , Nx )
888+ if p [py , px ] <= 0 :
889+ break
890+ peaks .append ((py - cy , px - cx ))
891+ yl = max (0 , py - suppress ); yh = min (Ny , py + suppress + 1 )
892+ xl = max (0 , px - suppress ); xh = min (Nx , px + suppress + 1 )
893+ p [yl :yh , xl :xh ] = 0.0
894+ return peaks
895+
896+ def test_square_lattice_gives_fourfold_fft_symmetry (self ):
897+ """z = cos(2πx/a) + cos(2πy/a) must give 4 FFT peaks at 90° intervals.
898+
899+ Physical context: square surface reconstructions (Si(001), oxide
900+ interfaces) display two orthogonal corrugations; their STM FFT has
901+ fourfold rotational symmetry.
902+
903+ Period 16 px on 128×128 → k_bin = 8 (exact integer FFT bin, zero
904+ discretisation error). Expected peak positions in (row_offset,
905+ col_offset) coordinates from DC: (0,±8) and (±8,0).
906+
907+ Each angle must fall within 5° of a multiple of 90°. A missing
908+ cosine term (only one direction) would give 2 peaks (180° spacing),
909+ not 4 (90° spacing), failing this test.
910+ """
911+ Ny , Nx = 128 , 128
912+ period_px = 16
913+ k_expected = float (Nx // period_px ) # 8.0
914+
915+ Y , X = np .mgrid [:Ny , :Nx ]
916+ arr = (
917+ np .cos (2.0 * np .pi * X / period_px )
918+ + np .cos (2.0 * np .pi * Y / period_px )
919+ ).astype (np .float64 )
920+
921+ peaks = self ._power_peaks (arr , n_peaks = 4 )
922+ assert len (peaks ) == 4 , f"Expected 4 peaks, found { len (peaks )} ."
923+
924+ for i , (dy , dx ) in enumerate (peaks ):
925+ radius = float (np .hypot (dy , dx ))
926+ assert abs (radius - k_expected ) < 0.5 , (
927+ f"Peak { i } at offset ({ dy } ,{ dx } ): radius { radius :.2f} px, "
928+ f"expected { k_expected :.1f} ± 0.5 px."
929+ )
930+ angle = float (np .degrees (np .arctan2 (dy , dx )))
931+ dev = float (abs (angle % 90.0 ))
932+ dev = min (dev , 90.0 - dev )
933+ assert dev < 5.0 , (
934+ f"Peak { i } at ({ dy } ,{ dx } ), angle { angle :.1f} °: { dev :.1f} ° from "
935+ "nearest 90° multiple; expected fourfold (square) symmetry. "
936+ "A missing cosine direction gives 2-fold (180°) symmetry."
937+ )
938+
939+ def test_hexagonal_lattice_gives_sixfold_fft_symmetry (self ):
940+ """Three equal-amplitude plane waves at 120° intervals must give 6
941+ FFT peaks at 60° intervals, all on the same reciprocal circle.
942+
943+ Physical context: hexagonal close-packed surfaces (Cu(111), Au(111),
944+ graphene, MoS₂) display sixfold symmetry in the STM FFT.
945+
946+ Wave vectors (magnitude 1/a, at 0°/120°/240°):
947+ k₁ = (1/a, 0)
948+ k₂ = (−1/(2a), √3/(2a))
949+ k₃ = (−1/(2a), −√3/(2a))
950+ Each cosine contributes two FFT peaks (±k), giving 6 peaks total.
951+ All six lie on the circle |k| = 1/a → FFT radius = Nx/a = 8 bins.
952+
953+ The ky bins for k₂/k₃ fall at ±4√3 ≈ ±6.928 (non-integer); the
954+ nearest integer bin is 7, so detected radius = √(7²+4²) ≈ 8.06 px.
955+ Radius tolerance ±1 px and angle tolerance ±10° accommodate this.
956+
957+ Adjacent-peak separation ≥ 8 px > suppress_r = 6 px, so all six
958+ peaks are found independently.
959+ """
960+ Ny , Nx = 128 , 128
961+ period_px = 16.0
962+ k_expected = float (Nx ) / period_px # 8.0
963+
964+ Y , X = np .mgrid [:Ny , :Nx ]
965+ sq3_half = np .sqrt (3.0 ) / 2.0
966+ arr = (
967+ np .cos (2.0 * np .pi * X / period_px )
968+ + np .cos (2.0 * np .pi * (- 0.5 * X + sq3_half * Y ) / period_px )
969+ + np .cos (2.0 * np .pi * (- 0.5 * X - sq3_half * Y ) / period_px )
970+ ).astype (np .float64 )
971+
972+ peaks = self ._power_peaks (arr , n_peaks = 6 )
973+ assert len (peaks ) == 6 , f"Expected 6 peaks, found { len (peaks )} ."
974+
975+ for i , (dy , dx ) in enumerate (peaks ):
976+ radius = float (np .hypot (dy , dx ))
977+ assert abs (radius - k_expected ) < 1.0 , (
978+ f"Peak { i } at ({ dy } ,{ dx } ): radius { radius :.3f} px, "
979+ f"expected { k_expected :.1f} ± 1 px. "
980+ "All 6 peaks should lie on the same reciprocal circle."
981+ )
982+
983+ angles = sorted (
984+ float (np .degrees (np .arctan2 (dy , dx ))) % 360.0
985+ for dy , dx in peaks
986+ )
987+ gaps = [angles [i + 1 ] - angles [i ] for i in range (5 )]
988+ gaps .append (360.0 - angles [5 ] + angles [0 ])
989+ for i , gap in enumerate (gaps ):
990+ assert abs (gap - 60.0 ) < 10.0 , (
991+ f"Angular gap { i } : { gap :.1f} °, expected 60° ± 10°. "
992+ f"Sorted peak angles: { [f'{ a :.0f} °' for a in angles ]} ."
993+ )
0 commit comments