From c3e8dc23f1a0bce0310f3414aeeea0279a78a7f3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek.necada@aalto.fi>
Date: Mon, 6 Feb 2017 09:39:43 +0200
Subject: [PATCH] Function to precalculate and save hex lattice translation
 coefficients.

Former-commit-id: e43dfa30224a4a7d79ee92b0c5ff1856d4f53da9
---
 qpms/hexpoints.py | 79 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 79 insertions(+)

diff --git a/qpms/hexpoints.py b/qpms/hexpoints.py
index 7a04db5..7780513 100644
--- a/qpms/hexpoints.py
+++ b/qpms/hexpoints.py
@@ -167,3 +167,82 @@ def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True,
                 'mi' : mi if mirrorindices else None,
                 'nmi' : nmi if mirrorindices else None
                }
+
+
+from qpms_c import get_mn_y
+from .qpms_p import Ã, B̃, cart2sph
+
+def hexlattice_precalc_AB_save(file, lMax, k_hexside, maxlayer, circular=True, savepointinfo = False, J_scat=3):
+    params = {
+        'lMax' : lMax,
+        'k_hexside' : k_hexside,
+        'maxlayer' : maxlayer,
+        'circular' : circular,
+        'savepointinfo' : savepointinfo,
+        'J_scat' : J_scat
+    }
+    tpdict = generate_trianglepoints(maxlayer, v3d=True, circular=circular, sixthindices=True, mirrorindices=True)
+    tphcdict = generate_trianglepoints_hexcomplement(maxlayer, v3d=True, circular=circular, thirdindices=True, mirrorindices=True)
+    my, ny = get_mn_y(lMax)
+    nelem = len(my)
+    a_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex)
+    b_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex)
+    a_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex)
+    b_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex)
+    a_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex)
+    b_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex)
+    a_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex)
+    b_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex)
+    
+    k_0 = k_hexside*_s3 # not really a wave vector here because of the normalisation!
+    
+    points = tpdict['points'][tpdict['nmi']]
+    for j in range(points.shape[0]):
+        d_i2j = cart2sph(points[j])
+        for yi in range(nelem):
+            for yj in range(nelem):
+                a_self_nm[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+                b_self_nm[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+    points = tpdict['points'][tpdict['mi'][0]]
+    for j in range(points.shape[0]):
+        d_i2j = cart2sph(points[j])
+        for yi in range(nelem):
+            for yj in range(nelem):
+                a_self_m0[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+                b_self_m0[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+    points = tphcdict['points'][tphcdict['nmi']]
+    for j in range(points.shape[0]):
+        d_i2j = cart2sph(points[j])
+        for yi in range(nelem):
+            for yj in range(nelem):
+                a_d2u_nm[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+                b_d2u_nm[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+    points = tphcdict['points'][tphcdict['mi'][0]]
+    for j in range(points.shape[0]):
+        d_i2j = cart2sph(points[j])
+        for yi in range(nelem):
+            for yj in range(nelem):
+                a_d2u_m0[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+                b_d2u_m0[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat)
+    tosave = {
+        'a_self_nm' : a_self_nm,
+        'a_self_m0' : a_self_m0,
+        'b_self_nm' : b_self_nm,
+        'b_self_m0' : b_self_m0,
+        'a_d2u_nm' : a_d2u_nm,
+        'a_d2u_m0' : a_d2u_m0,
+        'b_d2u_nm' : b_d2u_nm,
+        'b_d2u_m0' : b_d2u_m0,
+        'precalc_params' : params
+    }
+    if savepointinfo:
+        tosave['tp_points'] = tpdict['points'],
+        tosave['tp_si'] = tpdict['si'],
+        tosave['tp_mi'] = tpdict['mi'],
+        tosave['tp_nmi'] = tpdict['nmi']
+        tosave['tphc_points'] = tphcdict['points'],
+        tosave['tphc_ti'] = tphcdict['ti'],
+        tosave['tphc_mi'] = tphcdict['mi'],
+        tosave['tphc_nmi'] = tphcdict['nmi']
+    np.savez(file, **tosave)
+