LORENE
des_coupe_vect.C
1/*
2 * Copyright (c) 2000-2001 Eric Gourgoulhon
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: des_coupe_vect.C,v 1.7 2016/12/05 16:18:06 j_novak Exp $
27 * $Log: des_coupe_vect.C,v $
28 * Revision 1.7 2016/12/05 16:18:06 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.6 2014/10/13 08:53:22 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.5 2014/10/06 15:16:04 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.4 2008/08/19 06:42:00 j_novak
38 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
39 * cast-type operations, and constant strings that must be defined as const char*
40 *
41 * Revision 1.3 2004/03/25 10:29:24 j_novak
42 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
43 *
44 * Revision 1.2 2002/09/06 15:18:52 e_gourgoulhon
45 * Changement du nom de la variable "hz" en "hza"
46 * pour assurer la compatibilite avec le compilateur xlC_r
47 * sur IBM Regatta (le preprocesseur de ce compilateur remplace
48 * "hz" par "100").
49 *
50 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
51 * LORENE
52 *
53 * Revision 2.0 2000/03/01 16:11:40 eric
54 * *** empty log message ***
55 *
56 *
57 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_coupe_vect.C,v 1.7 2016/12/05 16:18:06 j_novak Exp $
58 *
59 */
60
61
62// Header C
63#include <cmath>
64
65// Header Lorene
66#include "tenseur.h"
67#include "graphique.h"
68#include "param.h"
69#include "utilitaires.h"
70#include "unites.h"
71
72namespace Lorene {
73//******************************************************************************
74
75void des_coupe_vect_x(const Tenseur& vv, double x0, double scale, double sizefl,
76 int nzdes, const char* title, const Cmp* defsurf, double zoom,
77 bool draw_bound, int ny, int nz) {
78
79 const Map& mp = *(vv.get_mp()) ;
80
81 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
82 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
83 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
84 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
85
86 ray = ( a1 > ray ) ? a1 : ray ;
87 ray = ( a2 > ray ) ? a2 : ray ;
88 ray = ( a3 > ray ) ? a3 : ray ;
89
90 ray *= zoom ;
91
92 double y_min = mp.get_ori_y() - ray ;
93 double y_max = mp.get_ori_y() + ray ;
94 double z_min = mp.get_ori_z() - ray ;
95 double z_max = mp.get_ori_z() + ray ;
96
97 des_coupe_vect_x(vv, x0, scale, sizefl, y_min, y_max, z_min, z_max, title,
98 defsurf, draw_bound, ny, nz) ;
99
100}
101
102
103
104//******************************************************************************
105
106void des_coupe_vect_x(const Tenseur& vv, double x0, double scale, double
107 sizefl, double y_min, double y_max, double z_min,
108 double z_max, const char* title, const Cmp* defsurf,
109 bool draw_bound, int ny, int nz) {
110
111 using namespace Unites ;
112
113 const Map& mp = *(vv.get_mp()) ;
114
115 if (vv.get_valence() != 1) {
116 cout <<
117 "des_coupe_vect_x: the Tenseur must be of valence 1 (vector) !" << endl ;
118 abort() ;
119 }
120
121 if ( vv.get_triad()->identify() != mp.get_bvect_cart().identify() ) {
122 cout <<
123 "des_coupe_vect_x: the vector must be given in Cartesian components !"
124 << endl ;
125 abort() ;
126 }
127
128
129 // Plot of the vector field
130 // ------------------------
131
132 float* vvy = new float[ny*nz] ;
133 float* vvz = new float[ny*nz] ;
134
135 double hy = (y_max - y_min) / double(ny-1) ;
136 double hza = (z_max - z_min) / double(nz-1) ;
137
138 for (int j=0; j<nz; j++) {
139
140 double z = z_min + hza * j ;
141
142 for (int i=0; i<ny; i++) {
143
144 double y = y_min + hy * i ;
145
146 // Computation of (r,theta,phi) :
147 double r, theta, phi ;
148 mp.convert_absolute(x0, y, z, r, theta, phi) ;
149
150 vvy[ny*j+i] = float(vv(1).val_point(r, theta, phi)) ;
151 vvz[ny*j+i] = float(vv(2).val_point(r, theta, phi)) ;
152
153 }
154 }
155
156 float ymin1 = float(y_min / km) ;
157 float ymax1 = float(y_max / km) ;
158 float zmin1 = float(z_min / km) ;
159 float zmax1 = float(z_max / km) ;
160
161 const char* nomy = "y [km]" ;
162 const char* nomz = "z [km]" ;
163
164 if (title == 0x0) {
165 title = "" ;
166 }
167
168 const char* device = 0x0 ;
169 int newgraph = ( (defsurf != 0x0) || draw_bound ) ? 1 : 3 ;
170
171 des_vect(vvy, vvz, ny, nz, ymin1, ymax1, zmin1, zmax1,
172 scale, sizefl, nomy, nomz, title, device, newgraph) ;
173
174
175 delete [] vvy ;
176 delete [] vvz ;
177
178 // Plot of the surface
179 // -------------------
180
181 if (defsurf != 0x0) {
182
183 assert(defsurf->get_mp() == vv.get_mp()) ;
184
185 newgraph = draw_bound ? 0 : 2 ;
186
187 des_surface_x(*defsurf, x0, device, newgraph) ;
188
189 } // End of the surface drawing
190
191
192 // Plot of the domains outer boundaries
193 // ------------------------------------
194
195 if (draw_bound) {
196
197 int ndom = mp.get_mg()->get_nzone() ; // total number of domains
198
199 for (int l=0; l<ndom-1; l++) { // loop on the domains (except the
200 // last one)
201
202 newgraph = (l == ndom-2) ? 2 : 0 ;
203
204 des_domaine_x(mp, l, x0, device, newgraph) ;
205 }
206 }
207
208
209}
210
211
212
213//******************************************************************************
214
215void des_coupe_vect_y(const Tenseur& vv, double y0, double scale, double sizefl,
216 int nzdes, const char* title, const Cmp* defsurf, double zoom,
217 bool draw_bound, int nx, int nz) {
218
219 const Map& mp = *(vv.get_mp()) ;
220
221 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
222 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
223 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
224 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
225
226 ray = ( a1 > ray ) ? a1 : ray ;
227 ray = ( a2 > ray ) ? a2 : ray ;
228 ray = ( a3 > ray ) ? a3 : ray ;
229
230 ray *= zoom ;
231
232 double x_min = mp.get_ori_x() - ray ;
233 double x_max = mp.get_ori_x() + ray ;
234 double z_min = mp.get_ori_z() - ray ;
235 double z_max = mp.get_ori_z() + ray ;
236
237
238 des_coupe_vect_y(vv, y0, scale, sizefl, x_min, x_max, z_min, z_max, title,
239 defsurf, draw_bound, nx, nz) ;
240
241}
242
243
244
245//******************************************************************************
246
247void des_coupe_vect_y(const Tenseur& vv, double y0, double scale, double
248 sizefl, double x_min, double x_max, double z_min,
249 double z_max, const char* title, const Cmp* defsurf,
250 bool draw_bound, int nx, int nz) {
251
252 using namespace Unites ;
253
254 const Map& mp = *(vv.get_mp()) ;
255
256 if (vv.get_valence() != 1) {
257 cout <<
258 "des_coupe_vect_y: the Tenseur must be of valence 1 (vector) !" << endl ;
259 abort() ;
260 }
261
262 if ( vv.get_triad()->identify() != mp.get_bvect_cart().identify() ) {
263 cout <<
264 "des_coupe_vect_y: the vector must be given in Cartesian components !"
265 << endl ;
266 abort() ;
267 }
268
269
270 // Plot of the vector field
271 // ------------------------
272
273 float* vvx = new float[nx*nz] ;
274 float* vvz = new float[nx*nz] ;
275
276 double hx = (x_max - x_min) / double(nx-1) ;
277 double hza = (z_max - z_min) / double(nz-1) ;
278
279 for (int j=0; j<nz; j++) {
280
281 double z = z_min + hza * j ;
282
283 for (int i=0; i<nx; i++) {
284
285 double x = x_min + hx * i ;
286
287 // Computation of (r,theta,phi) :
288 double r, theta, phi ;
289 mp.convert_absolute(x, y0, z, r, theta, phi) ;
290
291 vvx[nx*j+i] = float(vv(0).val_point(r, theta, phi)) ;
292 vvz[nx*j+i] = float(vv(2).val_point(r, theta, phi)) ;
293
294 }
295 }
296
297 float xmin1 = float(x_min / km) ;
298 float xmax1 = float(x_max / km) ;
299 float zmin1 = float(z_min / km) ;
300 float zmax1 = float(z_max / km) ;
301
302 const char* nomx = "x [km]" ;
303 const char* nomz = "z [km]" ;
304
305 if (title == 0x0) {
306 title = "" ;
307 }
308
309
310 const char* device = 0x0 ;
311 int newgraph = ( (defsurf != 0x0) || draw_bound ) ? 1 : 3 ;
312
313 des_vect(vvx, vvz, nx, nz, xmin1, xmax1, zmin1, zmax1,
314 scale, sizefl, nomx, nomz, title, device, newgraph) ;
315
316
317 delete [] vvx ;
318 delete [] vvz ;
319
320 // Plot of the surface
321 // -------------------
322
323 if (defsurf != 0x0) {
324
325 assert(defsurf->get_mp() == vv.get_mp()) ;
326
327 newgraph = draw_bound ? 0 : 2 ;
328
329 des_surface_y(*defsurf, y0, device, newgraph) ;
330
331 } // End of the surface drawing
332
333
334 // Plot of the domains outer boundaries
335 // ------------------------------------
336
337 if (draw_bound) {
338
339 int ndom = mp.get_mg()->get_nzone() ; // total number of domains
340
341 for (int l=0; l<ndom-1; l++) { // loop on the domains (except the
342 // last one)
343
344 newgraph = (l == ndom-2) ? 2 : 0 ;
345
346 des_domaine_y(mp, l, y0, device, newgraph) ;
347 }
348 }
349
350
351}
352
353
354//******************************************************************************
355
356void des_coupe_vect_z(const Tenseur& vv, double z0, double scale, double sizefl,
357 int nzdes, const char* title, const Cmp* defsurf, double zoom,
358 bool draw_bound, int nx, int ny) {
359
360 const Map& mp = *(vv.get_mp()) ;
361
362 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
363 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
364 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
365 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
366
367 ray = ( a1 > ray ) ? a1 : ray ;
368 ray = ( a2 > ray ) ? a2 : ray ;
369 ray = ( a3 > ray ) ? a3 : ray ;
370
371 ray *= zoom ;
372
373 double x_min = mp.get_ori_x() - ray ;
374 double x_max = mp.get_ori_x() + ray ;
375 double y_min = mp.get_ori_y() - ray ;
376 double y_max = mp.get_ori_y() + ray ;
377
378 des_coupe_vect_z(vv, z0, scale, sizefl, x_min, x_max, y_min, y_max, title,
379 defsurf, draw_bound, nx, ny) ;
380
381}
382
383
384
385//******************************************************************************
386
387void des_coupe_vect_z(const Tenseur& vv, double z0, double scale, double
388 sizefl, double x_min, double x_max, double y_min,
389 double y_max, const char* title, const Cmp* defsurf,
390 bool draw_bound, int nx, int ny) {
391
392 using namespace Unites ;
393
394 const Map& mp = *(vv.get_mp()) ;
395
396 if (vv.get_valence() != 1) {
397 cout <<
398 "des_coupe_vect_y: the Tenseur must be of valence 1 (vector) !" << endl ;
399 abort() ;
400 }
401
402 if ( vv.get_triad()->identify() != mp.get_bvect_cart().identify() ) {
403 cout <<
404 "des_coupe_vect_y: the vector must be given in Cartesian components !"
405 << endl ;
406 abort() ;
407 }
408
409
410 // Plot of the vector field
411 // ------------------------
412
413 float* vvx = new float[nx*ny] ;
414 float* vvy = new float[nx*ny] ;
415
416 double hy = (y_max - y_min) / double(ny-1) ;
417 double hx = (x_max - x_min) / double(nx-1) ;
418
419 for (int j=0; j<ny; j++) {
420
421 double y = y_min + hy * j ;
422
423 for (int i=0; i<nx; i++) {
424
425 double x = x_min + hx * i ;
426
427 // Computation of (r,theta,phi) :
428 double r, theta, phi ;
429 mp.convert_absolute(x, y, z0, r, theta, phi) ;
430
431 vvx[nx*j+i] = float(vv(0).val_point(r, theta, phi)) ;
432 vvy[nx*j+i] = float(vv(1).val_point(r, theta, phi)) ;
433
434 }
435 }
436
437 float ymin1 = float(y_min / km) ;
438 float ymax1 = float(y_max / km) ;
439 float xmin1 = float(x_min / km) ;
440 float xmax1 = float(x_max / km) ;
441
442 const char* nomy = "y [km]" ;
443 const char* nomx = "x [km]" ;
444
445 if (title == 0x0) {
446 title = "" ;
447 }
448
449 const char* device = 0x0 ;
450 int newgraph = ( (defsurf != 0x0) || draw_bound ) ? 1 : 3 ;
451
452 des_vect(vvx, vvy, nx, ny, xmin1, xmax1, ymin1, ymax1,
453 scale, sizefl, nomx, nomy, title, device, newgraph) ;
454
455
456 delete [] vvx ;
457 delete [] vvy ;
458
459 // Plot of the surface
460 // -------------------
461
462 if (defsurf != 0x0) {
463
464 assert(defsurf->get_mp() == vv.get_mp()) ;
465
466 newgraph = draw_bound ? 0 : 2 ;
467
468 des_surface_z(*defsurf, z0, device, newgraph) ;
469
470 } // End of the surface drawing
471
472 // Plot of the domains outer boundaries
473 // ------------------------------------
474
475 if (draw_bound) {
476
477 int ndom = mp.get_mg()->get_nzone() ; // total number of domains
478
479 for (int l=0; l<ndom-1; l++) { // loop on the domains (except the
480 // last one)
481
482 newgraph = (l == ndom-2) ? 2 : 0 ;
483
484 des_domaine_z(mp, l, z0, device, newgraph) ;
485 }
486 }
487
488}
489}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:304
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::z
Coord z
z coordinate centered on the grid
Definition map.h:740
Lorene::y
Coord y
y coordinate centered on the grid
Definition map.h:739
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:732
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:730