Vlasiator ebf0dd394 on dev (v5.4.0 + 1054 commits)
Loading...
Searching...
No Matches
vectorpotentialdipole_verify2.py
Go to the documentation of this file.
1#!/usr/bin/env python import matplotlib.pyplot as plt
2
3# /*
4# * This file is part of Vlasiator.
5# * Copyright 2010-2016 Finnish Meteorological Institute
6# * Copyright 2017-2019 University of Helsinki
7# *
8# * For details of usage, see the COPYING file and read the "Rules of the Road"
9# * at http://www.physics.helsinki.fi/vlasiator/
10# *
11# * This program is free software; you can redistribute it and/or modify
12# * it under the terms of the GNU General Public License as published by
13# * the Free Software Foundation; either version 2 of the License, or
14# * (at your option) any later version.
15# *
16# * This program is distributed in the hope that it will be useful,
17# * but WITHOUT ANY WARRANTY; without even the implied warranty of
18# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19# * GNU General Public License for more details.
20# *
21# * You should have received a copy of the GNU General Public License along
22# * with this program; if not, write to the Free Software Foundation, Inc.,
23# * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24# */
25import numpy as np
26import math
27import sys,os
28import pytools as pt
29import matplotlib.pyplot as plt
30import fieldmodels
31
32''' Testing routine for different dipole formulations
33
34 run using "python vectorpotentialdipole_verify.py [arg1] [arg2]" where arg1 is a number from 0 to 4
35 for different test profile starting positions, directions, and dipole tilts.
36
37 If arg2 is present, the code also calculates verification of derivative terms.
38
39 Generates plots of magnetic field components for different dipole models along cuts through the
40 simulation domain. For derivative analysis, calculates the derivatives along said cuts analytically
41 and numerically and outputs the ratio.
42
43'''
44
45if len(sys.argv)!=1:
46 testset = int(sys.argv[1])
47else:
48 testset = 0
49
50if len(sys.argv)!=2:
51 calcderivatives=True
52else:
53 calcderivatives=False
54
55plotmagnitude=False
56
57plt.switch_backend('Agg')
58
59outfilename = "./vecpotdip_verify2_"+str(testset)+".png"
60
61RE=6371000.
62#RE=1
63epsilon=1.e-15
64
65BGB=[0.,0.,0.]
66if testset==0:
67 tilt_angle_phi = 0.
68 tilt_angle_theta = 0.
69 line_theta = 0.
70 line_start = np.array([0,0,0])
71elif testset==1:
72 tilt_angle_phi = 0.
73 tilt_angle_theta = 0.
74 line_theta = 45.
75 line_start = np.array([0,0,0])
76elif testset==2:
77 tilt_angle_phi = 0.
78 tilt_angle_theta = 0.
79 line_theta = 0.
80 line_start = np.array([-3,-3,-3])
81elif testset==3:
82 tilt_angle_phi = 0.
83 tilt_angle_theta = 0.
84 line_theta = 45.
85 line_start = np.array([-3,-3,-3])
86elif testset==4:
87 tilt_angle_phi = 10
88 tilt_angle_theta = 0.
89 line_theta = 0.
90 line_start = np.array([0,0,0])
91elif testset==5:
92 tilt_angle_phi = 10
93 tilt_angle_theta = 45.
94 line_theta = 0.
95 line_start = np.array([0,0,0])
96elif testset==6:
97 tilt_angle_phi = 0
98 tilt_angle_theta = 5.
99 line_theta = 0.
100 line_start = np.array([0,0,0])
101 BGB=[0,0,-5e-9]
102else: # Same as 0
103 print("Default")
104 tilt_angle_phi = 0.
105 tilt_angle_theta = 0.
106 line_theta = 0.
107 line_start = np.array([0,0,0])
108
109print("Test set "+str(testset)+" line start "+str(line_start)+" tilt phi "+str(tilt_angle_phi)+" tilt theta "+str(tilt_angle_theta)+" line theta "+str(line_theta))
110
111line_phi = np.array([0,45,80,90,110,135])*math.pi/180.
112line_theta = np.zeros(len(line_phi)) + line_theta * math.pi/180.
113#line_start = np.array([-5,-5,-5])
114step = 0.1
115
116linewidth=2
117linthresh=1.e-10
118fontsize=20
119
120#fieldmodels.dipole.set_dipole(centerx, centery, centerz, tilt_phi, tilt_theta, mult=1.0, radius_f=None, radius_z=None):
121dip = fieldmodels.dipole(0,0,0,tilt_angle_phi,tilt_angle_theta)
122mdip = fieldmodels.dipole(80*RE,0,0,tilt_angle_phi,180.-tilt_angle_theta)
123imfpot = fieldmodels.IMFpotential(radius_z=10, radius_f=40, IMF=BGB)
124
125# Create figure
126fig = plt.figure()
127fig.set_size_inches(20,30)
128nsubplots=len(line_theta)
129for i in range(nsubplots):
130 fig.add_subplot(nsubplots,1,i+1)
131axes = fig.get_axes()
132
133radii = np.arange(0.1,100,step)*RE
134nr=len(radii)
135radiiRE = radii/RE
136
137fig.suptitle(r"Profiles starting from ("+str(line_start[0])+","+str(line_start[1])+","+str(line_start[2])+") [RE] with dipole tilt $\Phi="+str(int(tilt_angle_phi))+"$, $\Theta="+str(int(tilt_angle_theta))+"$", fontsize=fontsize)
138
139for i in range(nsubplots):
140 print("subplot ",i)
141 ax = axes[i]
142
143 ax.text(0.2,0.08,r"profile with $\theta="+str(int(line_theta[i]*180./math.pi))+"$, $\phi="+str(int(line_phi[i]*180./math.pi))+"$",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
144
145 xv = line_start[0]*RE + radii*np.sin(line_phi[i])*np.cos(line_theta[i])
146 yv = line_start[1]*RE + radii*np.sin(line_phi[i])*np.sin(line_theta[i])
147 zv = line_start[2]*RE + radii*np.cos(line_phi[i])
148
149 B1 = np.zeros([nr,4]) # X-scaled vector dipole
150 B2 = np.zeros([nr,4]) # regular dipole
151 B3 = np.zeros([nr,4]) # regular dipole + mirror dipole
152 B4 = np.zeros([nr,4]) # line dipole + mirror dipole
153
154 for j in range(nr):
155 for k in range(3):
156 B1[j,k] = dip.getX(xv[j],yv[j],zv[j],0,k,0) + imfpot.get(xv[j],yv[j],zv[j],0,k,0)
157 B2[j,k] = dip.get_old(xv[j],yv[j],zv[j],0,k,0)
158 B3[j,k] = B2[j,k] + mdip.get_old(xv[j],yv[j],zv[j],0,k,0)
159 B4[j,k] = dip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
160 B4[j,k] = B4[j,k] + mdip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
161 if plotmagnitude is True:
162 B1[j,3] = np.linalg.norm(B1[j,0:3])
163 B2[j,3] = np.linalg.norm(B2[j,0:3])
164 B3[j,3] = np.linalg.norm(B3[j,0:3])
165 B4[j,3] = np.linalg.norm(B4[j,0:3])
166
167 colors=['r','k','b','magenta']
168 coords = ['x','y','z','mag']
169
170 plotrange = range(3)
171 if plotmagnitude is True:
172 plotrange = range(4)
173 for k in plotrange:
174 ax.plot(radiiRE, B1[:,k], c=colors[k], linestyle='-', linewidth=linewidth, label='vectorpot B'+coords[k], zorder=-10)
175 ax.plot(radiiRE, B2[:,k], c=colors[k], linestyle='--', linewidth=linewidth, label='regular B'+coords[k])
176 ax.plot(radiiRE, B3[:,k], c=colors[k], linestyle=':', linewidth=linewidth, label='reg+mirror B'+coords[k])
177 if tilt_angle_phi<epsilon:
178 ax.plot(radiiRE, B4[:,k], c=colors[k], linestyle='-.', linewidth=linewidth, label='line+mirror B'+coords[k])
179
180 ax.set_xlabel(r"$r$ [$r_\mathrm{E}$]", fontsize=fontsize)
181 ax.set_xlim([0,70])
182 #ax.set_yscale('log', nonposy='clip')
183 ax.set_yscale('symlog', linthreshy=linthresh)
184 for item in ax.get_xticklabels():
185 item.set_fontsize(fontsize)
186 for item in ax.get_yticklabels():
187 item.set_fontsize(fontsize)
188
189 ylims = np.array(ax.get_ylim())
190 if ylims[0] < -1e-4:
191 ylims[0] = -1e-4
192 if ylims[1] > 1e-4:
193 ylims[1] = 1e-4
194 ax.set_ylim(ylims)
195
196handles, labels = axes[-1].get_legend_handles_labels()
197axes[-1].legend(handles, labels, fontsize=fontsize)
198
199fig.savefig(outfilename)
200plt.close()
201
202
203
204if calcderivatives:
205 # Derivatives
206 step2=0.00001 # distance in each direction for calculating numerical derivative
207 for kkk in range(3):
208 print("derivatives d"+coords[kkk])
209 dB1 = np.zeros([nr,3,3])
210 dB2 = np.zeros([nr,3,3])
211 dB3 = np.zeros([nr,3,3])
212 dB4 = np.zeros([nr,3,3])
213
214 # Create figure
215 fig = plt.figure()
216 fig.set_size_inches(20,30)
217 for i in range(nsubplots):
218 fig.add_subplot(nsubplots,1,i+1)
219 axes = fig.get_axes()
220
221 fig.suptitle(r"Numerical and analytical derivative ratios, profiles starting from ("+str(line_start[0])+","+str(line_start[1])+","+str(line_start[2])+") [RE] with dipole tilt $\Phi="+str(int(tilt_angle_phi))+"$, $\Theta="+str(int(tilt_angle_theta))+"$", fontsize=fontsize)
222
223 for i in range(nsubplots):
224 print("derivatives subplot ",i)
225 ax = axes[i]
226
227 xv = line_start[0]*RE + radii*np.sin(line_phi[i])*np.cos(line_theta[i])
228 yv = line_start[1]*RE + radii*np.sin(line_phi[i])*np.sin(line_theta[i])
229 zv = line_start[2]*RE + radii*np.cos(line_phi[i])
230
231 for j in range(nr):
232 for k in range(3):
233 B1[j,k] = dip.getX(xv[j],yv[j],zv[j],0,k,0) + imfpot.get(xv[j],yv[j],zv[j],0,k,0)
234 B2[j,k] = dip.get_old(xv[j],yv[j],zv[j],0,k,0)
235 B3[j,k] = B2[j,k] + mdip.get_old(xv[j],yv[j],zv[j],0,k,0)
236 B4[j,k] = dip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
237 B4[j,k] = B4[j,k] + mdip.get_ldp(xv[j],yv[j],zv[j],0,k,0)
238 #for kk in range(3):
239 kk=kkk
240 dB1[j,k,kk] = dip.getX(xv[j],yv[j],zv[j],1,k,kk) + imfpot.get(xv[j],yv[j],zv[j],1,k,kk)
241 dB2[j,k,kk] = dip.get_old(xv[j],yv[j],zv[j],1,k,kk)
242 dB3[j,k,kk] = dB2[j,k,kk] + mdip.get_old(xv[j],yv[j],zv[j],1,k,kk)
243 dB4[j,k,kk] = dip.get_ldp(xv[j],yv[j],zv[j],1,k,kk)
244 dB4[j,k,kk] = dB4[j,k,kk] + mdip.get_ldp(xv[j],yv[j],zv[j],1,k,kk)
245
246 # analytical derivative vs numerical derivative
247 for j in np.arange(1,nr-1):
248 for k in range(3):
249
250 # d/dx
251 if kkk==0:
252 #cdbx=(dip.getX(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.getX(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
253 cdbx=(dip.getX(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.getX(xv[j]-step2*RE,yv[j],zv[j],0,k,0) + imfpot.get(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - imfpot.get(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
254 if abs(cdbx) > epsilon*B1[j,k]:
255 dB1[j,k,0] = dB1[j,k,0]/cdbx
256 elif (abs(cdbx)<epsilon*B1[j,k]) and (abs(dB1[j,k,0])<epsilon*B1[j,k]):
257 dB1[j,k,0] = None
258 else:
259 dB1[j,k,0] = 0
260 if abs(B1[j,k])<epsilon:
261 dB1[j,k,0] = None
262
263 cdbx=(dip.get_old(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.get_old(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
264 if abs(cdbx) > epsilon*B2[j,k]:
265 dB2[j,k,0] = dB2[j,k,0]/cdbx
266 elif (abs(cdbx)<epsilon*B2[j,k]) and (abs(dB2[j,k,0])<epsilon*B2[j,k]):
267 dB2[j,k,0] = None
268 else:
269 dB2[j,k,0] = 0
270 if abs(B2[j,k])<epsilon:
271 dB2[j,k,0] = None
272
273 cdbx=(dip.get_old(xv[j]+step2*RE,yv[j],zv[j],0,k,0)+mdip.get_old(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.get_old(xv[j]-step2*RE,yv[j],zv[j],0,k,0) - mdip.get_old(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
274 if abs(cdbx) > epsilon*B3[j,k]:
275 dB3[j,k,0] = dB3[j,k,0]/cdbx
276 elif (abs(cdbx)<epsilon*B3[j,k]) and (abs(dB3[j,k,0])<epsilon*B3[j,k]):
277 dB3[j,k,0] = None
278 else:
279 dB3[j,k,0] = 0
280 if abs(B3[j,k])<epsilon:
281 dB3[j,k,0] = None
282
283 cdbx=(dip.get_ldp(xv[j]+step2*RE,yv[j],zv[j],0,k,0)+mdip.get_ldp(xv[j]+step2*RE,yv[j],zv[j],0,k,0) - dip.get_ldp(xv[j]-step2*RE,yv[j],zv[j],0,k,0) - mdip.get_ldp(xv[j]-step2*RE,yv[j],zv[j],0,k,0))/(2*step2*RE)
284 if abs(cdbx) > epsilon*B4[j,k]:
285 dB4[j,k,0] = dB4[j,k,0]/cdbx
286 elif (abs(cdbx)<epsilon*B4[j,k]) and (abs(dB4[j,k,0])<epsilon*B4[j,k]):
287 dB4[j,k,0] = None
288 else:
289 dB4[j,k,0] = 0
290 if abs(B4[j,k])<epsilon:
291 dB4[j,k,0] = None
292
293 # d/dy
294 if kkk==1:
295 cdby=(dip.getX(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.getX(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
296 if abs(cdby) > epsilon*B1[j,k]:
297 dB1[j,k,1] = dB1[j,k,1]/cdby
298 elif (abs(cdby)<epsilon*B1[j,k]) and (abs(dB1[j,k,1])<epsilon*B1[j,k]):
299 dB1[j,k,1] = None
300 else:
301 dB1[j,k,1] = 0
302 if abs(B1[j,k])<epsilon:
303 dB1[j,k,1] = None
304
305 cdby=(dip.get_old(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.get_old(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
306 if abs(cdby) > epsilon*B2[j,k]:
307 dB2[j,k,1] = dB2[j,k,1]/cdby
308 elif (abs(cdby)<epsilon*B2[j,k]) and (abs(dB2[j,k,1])<epsilon*B2[j,k]):
309 dB2[j,k,1] = None
310 else:
311 dB2[j,k,1] = 0
312 if abs(B2[j,k])<epsilon:
313 dB2[j,k,1] = None
314
315 cdby=(dip.get_old(xv[j],yv[j]+step2*RE,zv[j],0,k,0)+mdip.get_old(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.get_old(xv[j],yv[j]-step2*RE,zv[j],0,k,0) - mdip.get_old(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
316 if abs(cdby) > epsilon*B3[j,k]:
317 dB3[j,k,1] = dB3[j,k,1]/cdby
318 elif (abs(cdby)<epsilon*B3[j,k]) and (abs(dB3[j,k,1])<epsilon*B3[j,k]):
319 dB3[j,k,1] = None
320 else:
321 dB3[j,k,1] = 0
322 if abs(B3[j,k])<epsilon:
323 dB3[j,k,1] = None
324
325 cdby=(dip.get_ldp(xv[j],yv[j]+step2*RE,zv[j],0,k,0)+mdip.get_ldp(xv[j],yv[j]+step2*RE,zv[j],0,k,0) - dip.get_ldp(xv[j],yv[j]-step2*RE,zv[j],0,k,0) - mdip.get_ldp(xv[j],yv[j]-step2*RE,zv[j],0,k,0))/(2*step2*RE)
326 if abs(cdby) > epsilon*B4[j,k]:
327 dB4[j,k,1] = dB4[j,k,1]/cdby
328 elif (abs(cdby)<epsilon*B4[j,k]) and (abs(dB4[j,k,1])<epsilon*B4[j,k]):
329 dB4[j,k,1] = None
330 else:
331 dB4[j,k,1] = 0
332 if abs(B4[j,k])<epsilon:
333 dB4[j,k,1] = None
334
335 # d/dz
336 if kkk==2:
337 cdbz=(dip.getX(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.getX(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
338 if abs(cdbz) > epsilon*B1[j,k]:
339 dB1[j,k,2] = dB1[j,k,2]/cdbz
340 elif (abs(cdbz)<epsilon*B1[j,k]) and (abs(dB1[j,k,2])<epsilon*B1[j,k]):
341 dB1[j,k,2] = None
342 else:
343 dB1[j,k,2] = 0
344 if abs(B1[j,k])<epsilon:
345 dB1[j,k,2] = None
346
347 cdbz=(dip.get_old(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.get_old(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
348 if abs(cdbz) > epsilon*B2[j,k]:
349 dB2[j,k,2] = dB2[j,k,2]/cdbz
350 elif (abs(cdbz)<epsilon*B2[j,k]) and (abs(dB2[j,k,2])<epsilon*B2[j,k]):
351 dB2[j,k,2] = None
352 else:
353 dB2[j,k,2] = 0
354 if abs(B2[j,k])<epsilon:
355 dB2[j,k,2] = None
356
357 cdbz=(dip.get_old(xv[j],yv[j],zv[j]+step2*RE,0,k,0)+mdip.get_old(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.get_old(xv[j],yv[j],zv[j]-step2*RE,0,k,0) - mdip.get_old(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
358 if abs(cdbz) > epsilon*B3[j,k]:
359 dB3[j,k,2] = dB3[j,k,2]/cdbz
360 elif (abs(cdbz)<epsilon*B3[j,k]) and (abs(dB3[j,k,2])<epsilon*B3[j,k]):
361 dB3[j,k,2] = None
362 else:
363 dB3[j,k,2] = 0
364 if abs(B3[j,k])<epsilon:
365 dB3[j,k,2] = None
366
367 cdbz=(dip.get_ldp(xv[j],yv[j],zv[j]+step2*RE,0,k,0)+mdip.get_ldp(xv[j],yv[j],zv[j]+step2*RE,0,k,0) - dip.get_ldp(xv[j],yv[j],zv[j]-step2*RE,0,k,0) - mdip.get_ldp(xv[j],yv[j],zv[j]-step2*RE,0,k,0))/(2*step2*RE)
368 if abs(cdbz) > epsilon*B4[j,k]:
369 dB4[j,k,2] = dB4[j,k,2]/cdbz
370 elif (abs(cdbz)<epsilon*B4[j,k]) and (abs(dB4[j,k,2])<epsilon*B4[j,k]):
371 dB4[j,k,2] = None
372 else:
373 dB4[j,k,2] = 0
374 if abs(B4[j,k])<epsilon:
375 dB4[j,k,2] = None
376
377 print(np.ma.amin(np.ma.masked_invalid(dB1)),np.ma.amax(np.ma.masked_invalid(dB1)))
378 print(np.ma.amin(np.ma.masked_invalid(dB2)),np.ma.amax(np.ma.masked_invalid(dB2)))
379 print(np.ma.amin(np.ma.masked_invalid(dB3)),np.ma.amax(np.ma.masked_invalid(dB3)))
380 print(np.ma.amin(np.ma.masked_invalid(dB4)),np.ma.amax(np.ma.masked_invalid(dB4)))
381
382 # print(np.amin(B1),np.amax(B1))
383 # print(np.amin(B2),np.amax(B2))
384 # print(np.amin(B3),np.amax(B3))
385 # print(np.amin(B4),np.amax(B4))
386
387 colors=['r','k','b']
388 coords = ['x','y','z']
389
390 # for k in range(3):
391 # B1[:,k] = abs(B1[:,k])
392 # B2[:,k] = abs(B2[:,k])
393 # B3[:,k] = abs(B3[:,k])
394 # B4[:,k] = abs(B4[:,k])
395 # for kk in range(3):
396 # dB1[:,k,kk] = abs(dB1[:,k,kk])
397 # dB2[:,k,kk] = abs(dB2[:,k,kk])
398 # dB3[:,k,kk] = abs(dB3[:,k,kk])
399 # dB4[:,k,kk] = abs(dB4[:,k,kk])
400
401
402 for k in range(3):
403 ax.plot(radiiRE, dB1[:,k,kkk], c=colors[k], linestyle='-', linewidth=linewidth, label='vectorpot dB'+coords[k]+'/d'+coords[kkk])
404 ax.plot(radiiRE, dB2[:,k,kkk], c=colors[k], linestyle='--', linewidth=linewidth, label='regular dB'+coords[k]+'/d'+coords[kkk])
405 ax.plot(radiiRE, dB3[:,k,kkk], c=colors[k], linestyle=':', linewidth=linewidth, label='reg+mirror dB'+coords[k]+'/d'+coords[kkk])
406 if tilt_angle_phi<epsilon:
407 ax.plot(radiiRE, dB4[:,k,kkk], c=colors[k], linestyle='-.', linewidth=linewidth, label='line+mirror dB'+coords[k]+'/d'+coords[kkk])
408
409 ax.text(0.2,0.08,r"profile with $\theta="+str(int(line_theta[i]*180./math.pi))+"$, $\phi="+str(int(line_phi[i]*180./math.pi))+"$",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
410
411
412 ax.set_xlabel(r"$r$ [$r_\mathrm{E}$]", fontsize=fontsize)
413 ax.set_xlim([1,70])
414 #ax.set_yscale('log', nonposy='clip')
415 #ax.set_ylim([1.e-6,1.e-1])
416 #ax.set_ylim([1.e-26,1.e-21])
417 ax.set_ylim([-.1,1.1])
418
419 for item in ax.get_xticklabels():
420 item.set_fontsize(fontsize)
421 for item in ax.get_yticklabels():
422 item.set_fontsize(fontsize)
423
424 handles, labels = axes[-1].get_legend_handles_labels()
425 axes[-1].legend(handles, labels, fontsize=fontsize).set_zorder(10)
426
427 fig.savefig(outfilename[:-4]+"_d"+coords[kkk]+outfilename[-4:])
428 plt.close()
429
430
static ARCH_HOSTDEV VecSimple< T > abs(const VecSimple< T > &l)