41 from subprocess
import Popen, PIPE
42 from sys
import stdout
43 from time
import sleep
46 from matplotlib
import cm
47 import matplotlib.pyplot
as plt
54 clears the current figure (no arg) or specified window
66 def system(cmd, output=False):
68 Execute the external command
72 out = system("ls", out=True)
73 out = system("ls -l", out=True)
78 and get its stdout exitcode and stderr.
80 args = shlex.split(cmd)
81 proc = Popen(args, stdout=PIPE, stderr=PIPE)
82 out, err = proc.communicate()
83 exitcode = proc.returncode
86 out = out.split(
'\n')[:-1]
88 for i
in range(len(out)):
93 return out, exitcode, err
96 def pli(data, color='gist_earth', cmin=9998, cmax=9998, win=1, origin=None,
99 plots the transpose of the data
101 color maps keywords can be found in
102 http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
107 exec(
'options += ",vmin=cmin"')
110 exec(
'options += ",vmax=cmax"')
112 if (color == b
'yorick'):
116 if (aspect !=
'auto'):
117 aspect =
"\'" + aspect +
"\'"
121 exec(
'plt.matshow(data, aspect=' + aspect +
', fignum=win, cmap=color' + options +
129 binned = np.zeros(w.shape[0] * w.shape[1]).reshape(w.shape[0], w.shape[1])
131 for i
in range(w.shape[0]):
132 for j
in range(w.shape[1]):
133 binned[i, j] = w[i, j].sum() / (footprint * footprint + 0.0)
139 tabVect = np.reshape(tab, tab.size)
140 return [np.min(tabVect), np.max(tabVect)]
155 color = "0.71" [0-1] gray scale
159 http://matplotlib.org/api/colors_api.html
162 fig = plt.figure(win)
163 ax = fig.add_subplot(1, 1, 1)
167 print((
"Warning %dD dimensions. Cannot plot data. Use pli instead. " %
172 ax.plot(data, color=color)
174 ax.plot(x, data, color=color)
179 ax.set_xscale(
'linear')
184 ax.set_yscale(
'linear')
190 data = np.array(data)
191 if (len(data.shape) > 1):
192 print(
"oups zcen with dims > 1 not coded yet...")
195 for i
in range(len(data) - 1):
196 tmp = (float(data[i]) + float(data[i + 1])) / 2.
197 tmp2 = np.append(tmp2, tmp)
207 if ((nssp == 7)
and (rint > 0.285
and rint < 0.29)):
209 print(
"cas particulier")
210 x =
zcen(np.linspace(-1, 1, num=nssp + 1))
212 for i
in range(nssp):
213 xx = np.hstack((xx, x))
214 x = np.reshape(xx, (nssp, nssp))
216 r = np.sqrt(x * x + y * y)
217 valid2dext = ((r < rext)) * 1
218 valid2dint = ((r >= rint)) * 1
219 valid2d = valid2dint * valid2dext
224 valid = np.reshape(valid2d, [nssp * nssp])
226 return valid.tolist()
230 def plsh(slopesvector, nssp=14, rmax=0.98, obs=0, win=1, invertxy=False):
232 tmp = getValidSubapArray( nssp, rmax, obs);
233 X,Y = meshgrid(np.linspace(-1, 1, nssp), np.linspace(-1, 1, nssp))
234 vx = np.zeros([nssp*nssp])
235 vy = np.zeros([nssp*nssp])
237 vx.flat[hart] = slopesvector.flat[0:len(slopesvector)/2]
238 vy.flat[hart] = slopesvector.flat[len(slopesvector)/2+1:]
239 vx = vx.reshape([nssp, nssp])
240 vy = vy.reshape([nssp, nssp])
244 Q = quiver(X,Y, vy, vx)
246 Q = quiver(X,Y, vx, vy)
247 #qk = quiverkey(Q, 0.5, 0.92, 2, r'$2 \frac{m}{s}$', labelpos='W', fontproperties={'weight': 'bold'})
250 axis([l-0.05*dx, r+0.05*dx, b-0.05*dy, t+0.05*dy]) # MUST DO OTHERWISE THE AUTOSCALE CAN MISS SOME ARROWS
251 #title('Minimal arguments, no kwargs')
255 def plpyr(slopesvector, validArray):
257 wao.config.p_wfss[0]._isvalid
259 nslopes = slopesvector.shape[0] / 2
260 x, y = np.where(validArray.T)
261 plt.quiver(x, y, slopesvector[0:nslopes], slopesvector[nslopes:])
264 def plsh(slopesvector, nssp, validint, sparta=False, invertxy=False, returnquiver=False):
266 <slopesvector> is the input vector of slopes
267 <nssp> is the number of subapertures in the diameter of the pupil
268 <validint> is the normalized diameter of central obscuration (between 0 and 1.00)
269 <sparta> when==1, slopes are ordered xyxyxyxy...
270 when==0, slopes are xxxxxxyyyyyyy
271 <xy> when==1, swap x and y. Does nothing special when xy==0.
273 The routine plots a field vector of subaperture gradients defined in
274 vector <slopesvector>.
275 The routine automatically adjusts/finds what are the valid subapertures
276 for plotting, depending on the number of elements in <slopesvector>. Only the
277 devalidated subapertures inside the central obscuration cannot be
278 known, that’s why <validint> has to be passed in the argument list.
281 nsub = slopesvector.shape[0] // 2
282 x = np.linspace(-1, 1, nssp)
283 x, y = np.meshgrid(x, x)
284 r = np.sqrt(x * x + y * y)
289 rorder = np.sort(r.reshape(nssp * nssp))
291 ncentral = nssp * nssp - np.sum(r >= validint, dtype=np.int32)
294 validext = rorder[ncentral + nsub]
296 valid = (r < validext) & (r >= validint)
297 ivalid = np.where(valid)
299 vx = np.zeros([nssp, nssp])
300 vy = np.zeros([nssp, nssp])
301 if (sparta
is False):
303 vy[ivalid] = slopesvector[0:nsub]
304 vx[ivalid] = slopesvector[nsub:]
307 vx[ivalid] = slopesvector[0::2]
308 vy[ivalid] = slopesvector[1::2]
309 if (invertxy
is True):
317 plt.quiver(x, y, vx, vy, pivot=
'mid')
322 ir = pyfits.get_data("/home/fvidal/data/Run2015/June2015_27_onsky/ir/ir_2015-06-28_06h27m40s_script44_gain.fits")
324 JAMAIS TESTEE !!!!!!!!!!!!!!
327 X = np.arange(-5, 5, 0.25)
328 Y = np.arange(-5, 5, 0.25)
329 X, Y = np.meshgrid(X, Y)
331 plt.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet)
335 def FFThz(signal, fe, freq=0):
336 """ PSD = FFThz( signal, fe ) OU f = FFThz( 1024, fe, freq=1 )
337 On the first form, returns the power spectral density of signal.
338 If signal has units 'u', the PSD has units 'u^2/Hz'.
339 The frequency axis can be get by using the keyword freq=1."""
342 d = np.linspace(0, fe, n + 1)[0:n / 2 + 1]
346 d = np.abs(np.fft.fft(signal))[0:n / 2 + 1]
347 d = d**2 / (fe * n / 2)
353 if np.isscalar(wfsNum):
357 PSD =
FFThz(zerall[ii][izerNum, :], fe)
360 if (len(wfsNum) > 1):
361 ff =
FFThz(zerall[wfsNum][izerNum, :], fe, freq=1)
363 ff =
FFThz(zerall[wfsNum[0]][izerNum, :], fe, freq=1)
369 for i
in range(1, int(seconds)):
370 stdout.write(
"\r%d" % i)
377 fig = plt.matshow(pup)
378 pdiam = istart[1] - istart[0]
381 if (isvalid[i // pdiam, j // pdiam]):
386 plt.Rectangle((i - 0.5, j - 0.5), pdiam, pdiam, fill=
False,