LORENE
des_vect_bin.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_vect_bin_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_vect_bin.C,v 1.6 2014/10/13 08:53:23 j_novak Exp $" ;
24
25/*
26 * $Id: des_vect_bin.C,v 1.6 2014/10/13 08:53:23 j_novak Exp $
27 * $Log: des_vect_bin.C,v $
28 * Revision 1.6 2014/10/13 08:53:23 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.5 2014/10/06 15:16:05 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:25 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.1 2000/10/09 12:27:20 eric
51 * Ajout du test sur l'identite des triades de vv1 et vv2.
52 *
53 * Revision 2.0 2000/03/02 10:34:24 eric
54 * *** empty log message ***
55 *
56 *
57 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_vect_bin.C,v 1.6 2014/10/13 08:53:23 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_vect_bin_x(const Tenseur& vv1, const Tenseur& vv2, double x0,
76 double scale, double sizefl, double y_min, double y_max,
77 double z_min, double z_max, const char* title,
78 const Cmp* defsurf1, const Cmp* defsurf2,
79 bool draw_bound, int ny, int nz) {
80
81 using namespace Unites ;
82
83 const Map& mp1 = *(vv1.get_mp()) ;
84 const Map& mp2 = *(vv2.get_mp()) ;
85
86 // Protections
87 // -----------
88
89 if (vv1.get_valence() != 1) {
90 cout <<
91 "des_vect_bin_x: the Tenseur vv1 must be of valence 1 (vector) !" << endl ;
92 abort() ;
93 }
94
95 if ( vv1.get_triad()->identify() != mp1.get_bvect_cart().identify() ) {
96 cout <<
97 "des_vect_bin_x: the vector vv1 must be given in Cartesian components !"
98 << endl ;
99 abort() ;
100 }
101
102 if (vv2.get_valence() != 1) {
103 cout <<
104 "des_vect_bin_x: the Tenseur vv2 must be of valence 1 (vector) !" << endl ;
105 abort() ;
106 }
107
108 if ( vv2.get_triad()->identify() != mp2.get_bvect_cart().identify() ) {
109 cout <<
110 "des_vect_bin_x: the vector vv2 must be given in Cartesian components !"
111 << endl ;
112 abort() ;
113 }
114
115 if ( *(vv1.get_triad()) != *(vv2.get_triad()) ) {
116 cout <<
117 "des_vect_bin_x: the components of the two vectors are not defined"
118 << endl << "on the same triad !" << endl ;
119 abort() ;
120 }
121
122
123 // Plot of the vector field
124 // ------------------------
125
126 float* vvy = new float[ny*nz] ;
127 float* vvz = new float[ny*nz] ;
128
129 double hy = (y_max - y_min) / double(ny-1) ;
130 double hza = (z_max - z_min) / double(nz-1) ;
131
132 for (int j=0; j<nz; j++) {
133
134 double z = z_min + hza * j ;
135
136 for (int i=0; i<ny; i++) {
137
138 double y = y_min + hy * i ;
139
140 double r, theta, phi ;
141
142 mp1.convert_absolute(x0, y, z, r, theta, phi) ;
143 double vv1y = vv1(1).val_point(r, theta, phi) ;
144 double vv1z = vv1(2).val_point(r, theta, phi) ;
145
146 mp2.convert_absolute(x0, y, z, r, theta, phi) ;
147 double vv2y = vv2(1).val_point(r, theta, phi) ;
148 double vv2z = vv2(2).val_point(r, theta, phi) ;
149
150 vvy[ny*j+i] = float(vv1y + vv2y) ;
151 vvz[ny*j+i] = float(vv1z + vv2z) ;
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 = ( (defsurf1 != 0x0) || (defsurf2 != 0x0) || draw_bound ) ?
170 1 : 3 ;
171
172 des_vect(vvy, vvz, ny, nz, ymin1, ymax1, zmin1, zmax1,
173 scale, sizefl, nomy, nomz, title, device, newgraph) ;
174
175 delete [] vvy ;
176 delete [] vvz ;
177
178 // Plot of the surfaces
179 // --------------------
180
181 if (defsurf1 != 0x0) {
182
183 assert(defsurf1->get_mp() == vv1.get_mp()) ;
184 newgraph = ( (defsurf2 != 0x0) || draw_bound ) ? 0 : 2 ;
185 des_surface_x(*defsurf1, x0, device, newgraph) ;
186 }
187
188 if (defsurf2 != 0x0) {
189
190 assert(defsurf2->get_mp() == vv2.get_mp()) ;
191 newgraph = draw_bound ? 0 : 2 ;
192 des_surface_x(*defsurf2, x0, device, newgraph) ;
193 }
194
195
196 // Plot of the domains outer boundaries
197 // ------------------------------------
198
199 if (draw_bound) {
200
201 int ndom1 = mp1.get_mg()->get_nzone() ;
202 int ndom2 = mp2.get_mg()->get_nzone() ;
203
204 for (int l=0; l<ndom1-1; l++) { // loop on the domains (except the
205 // last one)
206 newgraph = 0 ;
207 des_domaine_x(mp1, l, x0, device, newgraph) ;
208 }
209
210 for (int l=0; l<ndom2-1; l++) { // loop on the domains (except the
211 // last one)
212
213 newgraph = (l == ndom2-2) ? 2 : 0 ;
214
215 des_domaine_x(mp2, l, x0, device, newgraph) ;
216 }
217
218 }
219}
220
221
222//******************************************************************************
223
224void des_vect_bin_y(const Tenseur& vv1, const Tenseur& vv2, double y0,
225 double scale, double sizefl, double x_min, double x_max,
226 double z_min, double z_max, const char* title,
227 const Cmp* defsurf1, const Cmp* defsurf2,
228 bool draw_bound, int nx, int nz) {
229
230 using namespace Unites ;
231
232 const Map& mp1 = *(vv1.get_mp()) ;
233 const Map& mp2 = *(vv2.get_mp()) ;
234
235 // Protections
236 // -----------
237
238 if (vv1.get_valence() != 1) {
239 cout <<
240 "des_vect_bin_y: the Tenseur vv1 must be of valence 1 (vector) !" << endl ;
241 abort() ;
242 }
243
244 if ( vv1.get_triad()->identify() != mp1.get_bvect_cart().identify() ) {
245 cout <<
246 "des_vect_bin_y: the vector vv1 must be given in Cartesian components !"
247 << endl ;
248 abort() ;
249 }
250
251 if (vv2.get_valence() != 1) {
252 cout <<
253 "des_vect_bin_y: the Tenseur vv2 must be of valence 1 (vector) !" << endl ;
254 abort() ;
255 }
256
257 if ( vv2.get_triad()->identify() != mp2.get_bvect_cart().identify() ) {
258 cout <<
259 "des_vect_bin_y: the vector vv2 must be given in Cartesian components !"
260 << endl ;
261 abort() ;
262 }
263
264 if ( *(vv1.get_triad()) != *(vv2.get_triad()) ) {
265 cout <<
266 "des_vect_bin_y: the components of the two vectors are not defined"
267 << endl << "on the same triad !" << endl ;
268 abort() ;
269 }
270
271 // Plot of the vector field
272 // ------------------------
273
274 float* vvx = new float[nx*nz] ;
275 float* vvz = new float[nx*nz] ;
276
277 double hx = (x_max - x_min) / double(nx-1) ;
278 double hza = (z_max - z_min) / double(nz-1) ;
279
280 for (int j=0; j<nz; j++) {
281
282 double z = z_min + hza * j ;
283
284 for (int i=0; i<nx; i++) {
285
286 double x = x_min + hx * i ;
287
288 double r, theta, phi ;
289
290 mp1.convert_absolute(x, y0, z, r, theta, phi) ;
291 double vv1x = vv1(0).val_point(r, theta, phi) ;
292 double vv1z = vv1(2).val_point(r, theta, phi) ;
293
294 mp2.convert_absolute(x, y0, z, r, theta, phi) ;
295 double vv2x = vv2(0).val_point(r, theta, phi) ;
296 double vv2z = vv2(2).val_point(r, theta, phi) ;
297
298 vvx[nx*j+i] = float(vv1x + vv2x) ;
299 vvz[nx*j+i] = float(vv1z + vv2z) ;
300 }
301 }
302
303 float xmin1 = float(x_min / km) ;
304 float xmax1 = float(x_max / km) ;
305 float zmin1 = float(z_min / km) ;
306 float zmax1 = float(z_max / km) ;
307
308 const char* nomx = "x [km]" ;
309 const char* nomz = "z [km]" ;
310
311 if (title == 0x0) {
312 title = "" ;
313 }
314
315 const char* device = 0x0 ;
316 int newgraph = ( (defsurf1 != 0x0) || (defsurf2 != 0x0) || draw_bound ) ?
317 1 : 3 ;
318
319 des_vect(vvx, vvz, nx, nz, xmin1, xmax1, zmin1, zmax1,
320 scale, sizefl, nomx, nomz, title, device, newgraph) ;
321
322 delete [] vvx ;
323 delete [] vvz ;
324
325 // Plot of the surfaces
326 // --------------------
327
328 if (defsurf1 != 0x0) {
329
330 assert(defsurf1->get_mp() == vv1.get_mp()) ;
331 newgraph = ( (defsurf2 != 0x0) || draw_bound ) ? 0 : 2 ;
332 des_surface_y(*defsurf1, y0, device, newgraph) ;
333 }
334
335 if (defsurf2 != 0x0) {
336
337 assert(defsurf2->get_mp() == vv2.get_mp()) ;
338 newgraph = draw_bound ? 0 : 2 ;
339 des_surface_y(*defsurf2, y0, device, newgraph) ;
340 }
341
342
343 // Plot of the domains outer boundaries
344 // ------------------------------------
345
346 if (draw_bound) {
347
348 int ndom1 = mp1.get_mg()->get_nzone() ;
349 int ndom2 = mp2.get_mg()->get_nzone() ;
350
351 for (int l=0; l<ndom1-1; l++) { // loop on the domains (except the
352 // last one)
353 newgraph = 0 ;
354 des_domaine_y(mp1, l, y0, device, newgraph) ;
355 }
356
357 for (int l=0; l<ndom2-1; l++) { // loop on the domains (except the
358 // last one)
359
360 newgraph = (l == ndom2-2) ? 2 : 0 ;
361
362 des_domaine_y(mp2, l, y0, device, newgraph) ;
363 }
364
365 }
366}
367
368
369//******************************************************************************
370
371void des_vect_bin_z(const Tenseur& vv1, const Tenseur& vv2, double z0,
372 double scale, double sizefl, double x_min, double x_max,
373 double y_min, double y_max, const char* title,
374 const Cmp* defsurf1, const Cmp* defsurf2,
375 bool draw_bound, int nx, int ny) {
376
377 using namespace Unites ;
378
379 const Map& mp1 = *(vv1.get_mp()) ;
380 const Map& mp2 = *(vv2.get_mp()) ;
381
382 // Protections
383 // -----------
384
385 if (vv1.get_valence() != 1) {
386 cout <<
387 "des_vect_bin_z: the Tenseur vv1 must be of valence 1 (vector) !" << endl ;
388 abort() ;
389 }
390
391 if ( vv1.get_triad()->identify() != mp1.get_bvect_cart().identify() ) {
392 cout <<
393 "des_vect_bin_z: the vector vv1 must be given in Cartesian components !"
394 << endl ;
395 abort() ;
396 }
397
398 if (vv2.get_valence() != 1) {
399 cout <<
400 "des_vect_bin_z: the Tenseur vv2 must be of valence 1 (vector) !" << endl ;
401 abort() ;
402 }
403
404 if ( vv2.get_triad()->identify() != mp2.get_bvect_cart().identify() ) {
405 cout <<
406 "des_vect_bin_z: the vector vv2 must be given in Cartesian components !"
407 << endl ;
408 abort() ;
409 }
410
411 if ( *(vv1.get_triad()) != *(vv2.get_triad()) ) {
412 cout <<
413 "des_vect_bin_z: the components of the two vectors are not defined"
414 << endl << "on the same triad !" << endl ;
415 abort() ;
416 }
417
418 // Plot of the vector field
419 // ------------------------
420
421 float* vvx = new float[nx*ny] ;
422 float* vvy = new float[nx*ny] ;
423
424 double hx = (x_max - x_min) / double(nx-1) ;
425 double hy = (y_max - y_min) / double(ny-1) ;
426
427 for (int j=0; j<ny; j++) {
428
429 double y = y_min + hy * j ;
430
431 for (int i=0; i<nx; i++) {
432
433 double x = x_min + hx * i ;
434
435 double r, theta, phi ;
436
437 mp1.convert_absolute(x, y, z0, r, theta, phi) ;
438 double vv1x = vv1(0).val_point(r, theta, phi) ;
439 double vv1y = vv1(1).val_point(r, theta, phi) ;
440
441 mp2.convert_absolute(x, y, z0, r, theta, phi) ;
442 double vv2x = vv2(0).val_point(r, theta, phi) ;
443 double vv2y = vv2(1).val_point(r, theta, phi) ;
444
445 vvx[nx*j+i] = float(vv1x + vv2x) ;
446 vvy[nx*j+i] = float(vv1y + vv2y) ;
447 }
448 }
449
450 float xmin1 = float(x_min / km) ;
451 float xmax1 = float(x_max / km) ;
452 float ymin1 = float(y_min / km) ;
453 float ymax1 = float(y_max / km) ;
454
455 const char* nomx = "x [km]" ;
456 const char* nomy = "y [km]" ;
457
458 if (title == 0x0) {
459 title = "" ;
460 }
461
462 const char* device = 0x0 ;
463 int newgraph = ( (defsurf1 != 0x0) || (defsurf2 != 0x0) || draw_bound ) ?
464 1 : 3 ;
465
466 des_vect(vvx, vvy, nx, ny, xmin1, xmax1, ymin1, ymax1,
467 scale, sizefl, nomx, nomy, title, device, newgraph) ;
468
469 delete [] vvx ;
470 delete [] vvy ;
471
472 // Plot of the surfaces
473 // --------------------
474
475 if (defsurf1 != 0x0) {
476
477 assert(defsurf1->get_mp() == vv1.get_mp()) ;
478 newgraph = ( (defsurf2 != 0x0) || draw_bound ) ? 0 : 2 ;
479 des_surface_z(*defsurf1, z0, device, newgraph) ;
480 }
481
482 if (defsurf2 != 0x0) {
483
484 assert(defsurf2->get_mp() == vv2.get_mp()) ;
485 newgraph = draw_bound ? 0 : 2 ;
486 des_surface_z(*defsurf2, z0, device, newgraph) ;
487 }
488
489
490 // Plot of the domains outer boundaries
491 // ------------------------------------
492
493 if (draw_bound) {
494
495 int ndom1 = mp1.get_mg()->get_nzone() ;
496 int ndom2 = mp2.get_mg()->get_nzone() ;
497
498 for (int l=0; l<ndom1-1; l++) { // loop on the domains (except the
499 // last one)
500 newgraph = 0 ;
501 des_domaine_z(mp1, l, z0, device, newgraph) ;
502 }
503
504 for (int l=0; l<ndom2-1; l++) { // loop on the domains (except the
505 // last one)
506
507 newgraph = (l == ndom2-2) ? 2 : 0 ;
508
509 des_domaine_z(mp2, l, z0, device, newgraph) ;
510 }
511
512 }
513}
514}
void des_domaine_x(const Map &mp, int l0, double x0, const char *device=0x0, int newgraph=3, double y_min=-1, double y_max=1, double z_min=-1, double z_max=1, const char *nomy=0x0, const char *nomz=0x0, const char *title=0x0, int nxpage=1, int nypage=1)
Basic routine for drawing the outer boundary of a given domain in a plane X=constant.
void des_domaine_y(const Map &mp, int l0, double y0, const char *device=0x0, int newgraph=3, double x_min=-1, double x_max=1, double z_min=-1, double z_max=1, const char *nomx=0x0, const char *nomz=0x0, const char *title=0x0, int nxpage=1, int nypage=1)
Basic routine for drawing the outer boundary of a given domain in a plane Y=constant.
void des_domaine_z(const Map &mp, int l0, double z0, const char *device=0x0, int newgraph=3, double x_min=-1, double x_max=1, double y_min=-1, double y_max=1, const char *nomx=0x0, const char *nomz=0x0, const char *title=0x0, int nxpage=1, int nypage=1)
Basic routine for drawing the outer boundary of a given domain in a plane Z=constant.
void des_vect(float *vvx, float *vvy, int nx, int ny, float xmin, float xmax, float ymin, float ymax, double scale, double sizefl, const char *nomx, const char *nomy, const char *title, const char *device=0x0, int newgraph=3, int nxpage=1, int nypage=1)
Basic routine for plotting vector field.
void des_surface_x(const Scalar &defsurf, double x0, const char *device=0x0, int newgraph=3, double y_min=-1, double y_max=1, double z_min=-1, double z_max=1, const char *nomy=0x0, const char *nomz=0x0, const char *title=0x0, int nxpage=1, int nypage=1)
Basic routine for drawing a stellar surface in a plane X=constant.
void des_surface_y(const Scalar &defsurf, double y0, const char *device=0x0, int newgraph=3, double x_min=-1, double x_max=1, double z_min=-1, double z_max=1, const char *nomx=0x0, const char *nomz=0x0, const char *title=0x0, int nxpage=1, int nypage=1)
Basic routine for drawing a stellar surface in a plane Y=constant.
void des_surface_z(const Scalar &defsurf, double z0, const char *device=0x0, int newgraph=3, double x_min=-1, double x_max=1, double y_min=-1, double y_max=1, const char *nomx=0x0, const char *nomz=0x0, const char *title=0x0, int nxpage=1, int nypage=1)
Basic routine for drawing a stellar surface in a plane Z=constant.
void des_vect_bin_x(const Tenseur &vv1, const Tenseur &vv2, double x0, double scale, double sizefl, double y_min, double y_max, double z_min, double z_max, const char *title, const Cmp *defsurf1=0x0, const Cmp *defsurf2=0x0, bool draw_bound=true, int ny=20, int nz=20)
Plots the sum of two vectors in a plane X=constant.
void des_vect_bin_y(const Tenseur &vv1, const Tenseur &vv2, double x0, double scale, double sizefl, double x_min, double x_max, double z_min, double z_max, const char *title, const Cmp *defsurf1=0x0, const Cmp *defsurf2=0x0, bool draw_bound=true, int nx=20, int nz=20)
Plots the sum of two vectors in a plane Y=constant.
void des_vect_bin_z(const Tenseur &vv1, const Tenseur &vv2, double x0, double scale, double sizefl, double x_min, double x_max, double y_min, double y_max, const char *title, const Cmp *defsurf1=0x0, const Cmp *defsurf2=0x0, bool draw_bound=true, int nx=20, int ny=20)
Plots the sum of two vectors in a plane Z=constant.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.