From 0f2e8e700e846569bca685e136aaab5937c9785b Mon Sep 17 00:00:00 2001
From: Beale John Henry <beale_j@saresc-cons-04.psi.ch>
Date: Wed, 12 Jul 2023 15:56:12 +0200
Subject: [PATCH] raw code from Thomas Barends for visulaising preferential
 orientations from CrystFEL .stream files

---
 cool_tools/plot_orientations_orthographic.py | 191 +++++++++++++++++++
 1 file changed, 191 insertions(+)
 create mode 100644 cool_tools/plot_orientations_orthographic.py

diff --git a/cool_tools/plot_orientations_orthographic.py b/cool_tools/plot_orientations_orthographic.py
new file mode 100644
index 0000000..b8231e5
--- /dev/null
+++ b/cool_tools/plot_orientations_orthographic.py
@@ -0,0 +1,191 @@
+# IMPORTant stuff....
+import string
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D  
+
+
+def ProjectionPlot(X,Y,Z,C,title):
+    f=plt.figure()
+    ax1=f.add_subplot(221)
+    ax1.scatter(X,Z,s=0.5,c=C)
+    ax1.set_title(title+',XZ plane')
+    ax1.set_xlabel('X ')
+    ax1.set_ylabel('Z (along beam)')
+    ax1.set_aspect('equal', 'box')
+    
+    ax2=f.add_subplot(222,projection='3d')
+    ax2.scatter(X,Z,-Y,s=0.5,c=C)
+    ax2.set_xlabel('X ')
+    ax2.set_ylabel('Z (along beam)')
+    ax2.set_zlabel('Y (along jet)')
+    ax2.set_aspect('equal', 'box')
+    
+    ax3=f.add_subplot(223)
+    ax3.scatter(X,-Y,c=C,s=0.5)
+    ax3.set_title(title+',XY plane')
+    ax3.set_xlabel('X ')
+    ax3.set_ylabel('Y (along jet)')
+    ax3.set_aspect('equal', 'box')
+    
+    ax4=f.add_subplot(224)
+    ax4.scatter(Z,-Y,c=C,s=0.5)
+    ax4.set_title(title+',ZY plane')
+    ax4.set_xlabel('Z (along beam)')
+    ax4.set_ylabel('Y (along jet)')
+    ax4.set_aspect('equal', 'box')
+
+    f.canvas.set_window_title(title)
+    f.set_size_inches(9,7)
+
+    plt.tight_layout()
+    plt.savefig(title+'.png')
+    plt.show()
+    
+    
+
+# declare lists to read reciprocal vectors
+astar=[]
+bstar=[]
+cstar=[]
+
+# set path and filenames
+#pathname='/lfs/l1/agorel/2018-03-20-shock-wave-8ns-cxim1816/2-indexing-pump-probe-zaef-nomulti/'
+#filename='indexed_65T_pumpprobe_SNR4_THR100.stream.rogue.stream'
+
+#pathname='/lfs/l1/agorel/2018-03-20-shock-wave-8ns-cxim1816/2-indexing-probe-zaef-nomulti/'
+#filename='indexed_65T_probe_SNR4_THR100.stream'
+
+pathname='./'
+filename='FAP-dark.stream'
+
+#pathname='./'
+#filename='smalla.stream'
+
+#reindex=True
+reindex=False
+
+# open stream file, read all lines and close again
+infile=open(pathname+filename,'r')
+inlines=infile.readlines()
+infile.close()
+
+
+
+# loop over lines and extract reciprocal cell axes 
+for line in inlines:
+    if 'astar' in line:
+        elements=line.split()
+        vector=[float(elements[2]),float(elements[3]),float(elements[4])]
+        astar.append(vector)
+    
+    if 'bstar' in line:
+        elements=line.split()
+        vector=[float(elements[2]),float(elements[3]),float(elements[4])]
+        bstar.append(vector)
+
+    if 'cstar' in line:
+        elements=line.split()
+        vector=[float(elements[2]),float(elements[3]),float(elements[4])]
+        cstar.append(vector)
+
+# turn all into numpy arrays
+astararray=np.asarray(astar)
+bstararray=np.asarray(bstar)
+cstararray=np.asarray(cstar)
+
+# get number of cells from shape of array
+s=np.shape(astararray)
+numcells=s[0]
+
+# declare lists for real-space cell vector directions (i.e. normalized vectors)
+adir=[]
+bdir=[]
+cdir=[]
+
+# declare lists for vector lengths
+alen=[]
+blen=[]
+clen=[]
+
+# loop over all cells
+for n in range(numcells):
+    astar=astararray[n,:]
+    bstar=bstararray[n,:]
+    cstar=cstararray[n,:]
+    Vstar=np.dot( astar , np.cross( bstar , cstar )  )    # reciprocal cell volume
+    a=(np.cross(bstar,cstar) / Vstar )                    # calculate real-space a
+    b=(np.cross(cstar,astar) / Vstar )                    # calculate real-space b
+    c=(np.cross(astar,bstar) / Vstar )                    # calculate real-space c
+
+    al=np.linalg.norm(a)                                  # calculate lengths 
+    bl=np.linalg.norm(b)
+    cl=np.linalg.norm(c)
+
+    if al<bl and reindex==True:
+        newa=-1.0*b
+        newb=a
+        newc=c
+         
+        newal=bl
+        newbl=al
+        newcl=cl
+        
+        a=newa
+        b=newb
+        c=newc
+        al=newal
+        bl=newbl
+        cl=newcl 
+
+    adir.append(a/al)                                     # put normalized vectors in lists
+    bdir.append(b/bl)
+    cdir.append(c/cl)
+    
+    alen.append(al)                                       # put lengths in lists
+    blen.append(bl)
+    clen.append(cl)
+    
+
+adir=np.asarray(adir)                                     # turn lists into numpy arrays
+bdir=np.asarray(bdir)
+cdir=np.asarray(cdir)
+
+alen=np.asarray(alen)
+blen=np.asarray(blen)
+clen=np.asarray(clen)
+
+f,((ax1,ax2,ax3))=plt.subplots(1,3)
+
+ax1.hist(10.*alen,int(5*np.sqrt(len(alen))),normed=1,facecolor='green',edgecolor='green')
+atitle=r'$\mathit{a}$'
+ax1.set_title(atitle)
+
+ax2.hist(10.*blen,int(5*np.sqrt(len(blen))),normed=1,facecolor='green',edgecolor='green')
+btitle=r'$\mathit{b}$'
+ax2.set_title(btitle)
+
+ax3.hist(10.*clen,int(5*np.sqrt(len(clen))),normed=1,facecolor='green',edgecolor='green')
+ctitle=r'$\mathit{c}$'
+ax3.set_title(ctitle)
+
+
+
+f.canvas.set_window_title('Unit Cell Parameter Histograms')
+f.set_size_inches(10,5)
+plt.tight_layout()
+plt.show()
+
+step=50                                                    # plot every nth point only
+
+# plot a, coloured by length of a (normalized to maximum length observed)
+print 'I will plot every',step,'th point,',len(adir[::step,0]),'points in total...'
+
+ProjectionPlot(adir[::step,0],adir[::step,1],adir[::step,2],(alen[::step]/np.max(alen)),'a_axis_direction')
+
+ProjectionPlot(bdir[::step,0],bdir[::step,1],bdir[::step,2],(blen[::step]/np.max(blen)),'b_axis_direction')
+
+ProjectionPlot(cdir[::step,0],cdir[::step,1],cdir[::step,2],(clen[::step]/np.max(clen)),'c_axis_direction')
+
+
+# FIN....
\ No newline at end of file