forked from Probe-Particle/ppafm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplotZ.py
executable file
·158 lines (119 loc) · 6.32 KB
/
plotZ.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#!/usr/bin/python
# This is a sead of simple plotting script which should get AFM frequency delta 'df.xsf' and generate 2D plots for different 'z'
import os
import sys
import __main__ as main
import numpy as np
import matplotlib.pyplot as plt
#import GridUtils as GU
import pyProbeParticle as PPU
import pyProbeParticle.GridUtils as GU
from scipy.interpolate import interp1d
from optparse import OptionParser
def find_minimum(array,precision=0.0001):
i=1
while i< len(array):
if (array[i-1] - array[i]) > precision and ( array[i+1] - array[i]) > precision:
return i
i+=1
HELP_MSG="""Use this program in the following way:
"""+os.path.basename(main.__file__) +""" -p "X1xY1" [-p "X2xY2" ...] """
parser = OptionParser()
parser.add_option( "-k", action="store", type="float", help="tip stiffenss [N/m]" )
parser.add_option( "--krange", action="store", type="float", help="tip stiffenss range (min,max,n) [N/m]", nargs=3)
parser.add_option( "-q", action="store", type="float", help="tip charge [e]" )
parser.add_option( "--qrange", action="store", type="float", help="tip charge range (min,max,n) [e]", nargs=3)
parser.add_option( "-a", action="store", type="float", help="oscilation amplitude [A]" )
parser.add_option( "--arange", action="store", type="float", help="oscilation amplitude range (min,max,n) [A]", nargs=3)
parser.add_option("-p", "--points", default=[], type=str, help="Point where to perform Z-scan", action="append")
parser.add_option( "--npy" , action="store_true" , help="load and save fields in npy instead of xsf" , default=False)
#parser.add_option( "-y", action="store", type="float", help="format of input file")
#parser.add_option( "--yrange", action="store", type="float", help="y positions of the tip range (min,max,n) [A]", nargs=3)
(options, args) = parser.parse_args()
opt_dict = vars(options)
print options
if options.npy:
data_format ="npy"
else:
data_format ="xsf"
if options.points==[]:
sys.exit(HELP_MSG)
print " >> OVEWRITING SETTINGS by params.ini "
PPU.loadParams( 'params.ini' )
dz = PPU.params['scanStep'][2]
Amp = [ PPU.params['Amplitude'] ]
scan_max=PPU.params['scanMax'][2]
scan_min=PPU.params['scanMin'][2]
scan_step=PPU.params['scanStep'][2]
print " >> OVEWRITING SETTINGS by command line arguments "
if opt_dict['krange'] is not None:
Ks = np.linspace( opt_dict['krange'][0], opt_dict['krange'][1], opt_dict['krange'][2] )
elif opt_dict['k'] is not None:
Ks = [ opt_dict['k'] ]
else:
Ks = [ PPU.params['stiffness'][0] ]
# Qs
if opt_dict['qrange'] is not None:
Qs = np.linspace( opt_dict['qrange'][0], opt_dict['qrange'][1], opt_dict['qrange'][2] )
elif opt_dict['q'] is not None:
Qs = [ opt_dict['q'] ]
else:
Qs = [ PPU.params['charge'] ]
# Amps
if opt_dict['arange'] is not None:
Amps = np.linspace( opt_dict['arange'][0], opt_dict['arange'][1], opt_dict['arange'][2] )
elif opt_dict['a'] is not None:
Amps = [ opt_dict['a'] ]
else:
Amps = [ PPU.params['Amplitude'] ]
for iq,Q in enumerate( Qs ):
for ik,K in enumerate( Ks ):
dirname = "Q%1.2fK%1.2f" %(Q,K)
print "Working in {} directory".format(dirname)
fzs,lvec,nDim,head=GU.load_scal_field(dirname+'/OutFz', data_format=data_format)
dfs = PPU.Fz2df( fzs, dz = dz, k0 = PPU.params['kCantilever'], f0=PPU.params['f0Cantilever'], n=Amp/dz )
for p in options.points:
x=float(p.split('x')[0])
y=float(p.split('x')[1])
x_pos=int(x/scan_step)
y_pos=int(y/scan_step)
Zplot=np.zeros(fzs.shape[0])
Fplot=np.zeros(fzs.shape[0])
DFplot=np.zeros(fzs.shape[0])
for k in range(0,fzs.shape[0]):
Fplot[k]=fzs[-k-1][y_pos][x_pos]
Zplot[k]=scan_max-scan_step*k
# shifting the df plot
for k in range(0,dfs.shape[0]-1):
DFplot[k+(int)(Amp/scan_step/2)]=dfs[-k-1][y_pos][x_pos]
xnew = np.linspace(Zplot[0], Zplot[-1], num=41, endpoint=True)
F_interp=interp1d(Zplot, Fplot,kind='cubic')
fig,ax1 = plt.subplots()
ax1.plot(Zplot, Fplot, 'ko', xnew, F_interp(xnew), 'k--')
ax1.set_xlabel('Z coordinate of the tip ($\AA$)')
ax1.set_ylabel('Force (eV/$\AA$)', color='black')
for tl in ax1.get_yticklabels():
tl.set_color('black')
F_interp=interp1d(Zplot, DFplot,kind='cubic')
ax2=ax1.twinx()
# min_index= np.argmin(DFplot)
min_index= find_minimum(DFplot)
# print "MIN", min_index
# print DFplot
ax2.plot(Zplot, DFplot,'bo', xnew, F_interp(xnew), 'b--')
axes = plt.gca()
ax2.set_ylabel('Frequency shift (Hz)', color='b')
for tl in ax2.get_yticklabels():
tl.set_color('b')
ax2.text(Zplot[min_index]+0.02, DFplot[min_index]-1.0,
'x:{:4.2f} ($\AA$); y:{:4.2f} (Hz)'.format(Zplot[min_index],
DFplot[min_index]), style='italic',
bbox={'facecolor':'blue', 'alpha':0.5, 'pad':0})
plt.axhline(y=0, color='black', ls='-.')
perplane=fig.add_axes([0.65, 0.6, 0.25, 0.25])
# perplane.imshow(dfs[min_index+int(0.5/scan_step),:, :], origin='image', cmap='gray')
# perplane.imshow(dfs[len(Zplot)-min_index-(int)(Amp/scan_step/2)-5, :, :], origin='image', cmap='gray')
perplane.imshow(dfs[len(Zplot)-min_index-(int)(Amp/scan_step/2)-int(1.0/scan_step), :, :], origin='image', cmap='gray')
perplane.scatter(x=x_pos, y=y_pos, s=50, c='red', alpha=0.8)
perplane.axis('off')
plt.show()