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 
23 char 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 
72 namespace Lorene {
73 //******************************************************************************
74 
75 void 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 
224 void 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 
371 void 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.