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