LORENE
des_prof_scalar.C
1 /*
2  * Draws the profile of a {\tt Scalar} along some radial axis determined by
3  * a fixed value of $(\theta, \phi)$.
4  */
5 
6 /*
7  * Copyright (c) 2004-2005 Eric Gourgoulhon & Philippe Grandclement
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE 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  * LORENE 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
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 
28 char des_prof_scalar_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_prof_scalar.C,v 1.12 2014/10/13 08:53:22 j_novak Exp $" ;
29 
30 
31 /*
32  * $Id: des_prof_scalar.C,v 1.12 2014/10/13 08:53:22 j_novak Exp $
33  * $Log: des_prof_scalar.C,v $
34  * Revision 1.12 2014/10/13 08:53:22 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.11 2014/10/06 15:16:05 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.10 2012/01/17 10:35:40 j_penner
41  * added point plot
42  *
43  * Revision 1.9 2008/08/19 06:42:00 j_novak
44  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
45  * cast-type operations, and constant strings that must be defined as const char*
46  *
47  * Revision 1.8 2005/03/25 19:57:04 e_gourgoulhon
48  * Added plot of domain boundaries (new argument draw_bound).
49  *
50  * Revision 1.7 2004/05/17 19:47:25 e_gourgoulhon
51  * -- Function des_profile_mult(const Scalar**,...): added argument
52  * device.
53  * -- Functions des_meridian: added arguments device and closeit.
54  *
55  * Revision 1.6 2004/04/06 07:47:29 j_novak
56  * Added a #include "string.h"
57  *
58  * Revision 1.5 2004/04/05 14:42:02 e_gourgoulhon
59  * Added functions des_meridian.
60  *
61  * Revision 1.4 2004/02/17 22:18:00 e_gourgoulhon
62  * Changed prototype of des_profile_mult (added radial_scale, theta and
63  * phi can now vary from one profile to the other, etc...)
64  *
65  * Revision 1.3 2004/02/16 13:23:33 e_gourgoulhon
66  * Function des_profile_mult: added delete [] uutab at the end.
67  *
68  * Revision 1.2 2004/02/15 21:57:45 e_gourgoulhon
69  * des_profile_mult: changed argument Scalar* to Scalar**.
70  *
71  * Revision 1.1 2004/02/12 16:21:28 e_gourgoulhon
72  * Functions des_profile for Scalar's transfered from file des_prof_cmp.C.
73  * Added new function des_profile_mult.
74  *
75  *
76  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_prof_scalar.C,v 1.12 2014/10/13 08:53:22 j_novak Exp $
77  *
78  */
79 
80 // Header C
81 #include <cmath>
82 #include <cstring>
83 
84 // Header Lorene
85 #include "scalar.h"
86 #include "graphique.h"
87 
88 #include <vector>
89 namespace Lorene {
90 //******************************************************************************
91 // VERSION SCALAR SANS UNITES
92 
93 void des_profile(const Scalar& uu, double r_min, double r_max,
94  double theta, double phi, const char* nomy, const char* title,
95  bool draw_bound) {
96 
97 
98  const int npt = 400 ; // Number of points along the axis
99 
100  float uutab[npt] ; // Value of uu at the npt points
101 
102  double hr = (r_max - r_min) / double(npt-1) ;
103 
104  for (int i=0; i<npt; i++) {
105 
106  double r = hr * i + r_min ;
107 
108  uutab[i] = float(uu.val_point(r, theta, phi)) ;
109 
110  }
111 
112  float xmin = float(r_min) ;
113  float xmax = float(r_max) ;
114 
115  const char* nomx = "r" ;
116 
117  if (title == 0x0) {
118  title = "" ;
119  }
120 
121  if (nomy == 0x0) {
122  nomy = "" ;
123  }
124 
125  // Preparations for the drawing of boundaries
126  // ------------------------------------------
127  const Map& mp = uu.get_mp() ;
128  int nz = mp.get_mg()->get_nzone() ;
129  int l_max = (mp.get_mg()->get_type_r(nz-1) == UNSURR) ? nz-2 : nz-1 ;
130 
131  float* xbound = new float[l_max+1] ;
132  int nbound = 0 ;
133 
134  if (draw_bound) {
135  const double xi_max = 1. ;
136  for (int l=0; l<=l_max; l++) {
137 
138  double rb = mp.val_r(l, xi_max, theta, phi) ;
139 
140  if ((rb >= r_min) && (rb <= r_max)) {
141  xbound[nbound] = float(rb) ;
142  nbound++ ;
143  }
144  }
145  }
146 
147  des_profile(uutab, npt, xmin, xmax, nomx, nomy, title, 0x0,
148  nbound, xbound) ;
149 
150  delete [] xbound ;
151 
152 }
153 
154 //******************************************************************************
155 
156 void des_profile(const Scalar& uu, double r_min, double r_max, double scale,
157  double theta, double phi, const char* nomx, const char* nomy,
158  const char* title, bool draw_bound) {
159 
160 
161  const int npt = 400 ; // Number of points along the axis
162 
163  float uutab[npt] ; // Value of uu at the npt points
164 
165  double hr = (r_max - r_min) / double(npt-1) ;
166 
167  for (int i=0; i<npt; i++) {
168 
169  double r = hr * i + r_min ;
170 
171  uutab[i] = float(uu.val_point(r, theta, phi)) ;
172 
173  }
174 
175  float xmin = float(r_min * scale) ;
176  float xmax = float(r_max * scale) ;
177 
178 
179  if (title == 0x0) {
180  title = "" ;
181  }
182 
183  if (nomx == 0x0) {
184  nomx = "" ;
185  }
186 
187  if (nomy == 0x0) {
188  nomy = "" ;
189  }
190 
191  // Preparations for the drawing of boundaries
192  // ------------------------------------------
193  const Map& mp = uu.get_mp() ;
194  int nz = mp.get_mg()->get_nzone() ;
195  int l_max = (mp.get_mg()->get_type_r(nz-1) == UNSURR) ? nz-2 : nz-1 ;
196 
197  float* xbound = new float[l_max+1] ;
198  int nbound = 0 ;
199 
200  if (draw_bound) {
201  const double xi_max = 1. ;
202  for (int l=0; l<=l_max; l++) {
203 
204  double rb = mp.val_r(l, xi_max, theta, phi) ;
205 
206  if ((rb >= r_min) && (rb <= r_max)) {
207  xbound[nbound] = float(rb) ;
208  nbound++ ;
209  }
210  }
211  }
212 
213  // Call to the low level routine
214  // -----------------------------
215  des_profile(uutab, npt, xmin, xmax, nomx, nomy, title, 0x0,
216  nbound, xbound) ;
217 
218  delete [] xbound ;
219 
220 }
221 
222 
223 //******************************************************************************
224 
225 void des_profile_mult(const Scalar** uu, int nprof, double r_min, double r_max,
226  const double* theta, const double* phi, double radial_scale,
227  bool closeit, const char* nomy, const char* title, int ngraph,
228  const char* nomx, const int* line_style, const char* device,
229  bool draw_bound) {
230 
231  // Special case of no graphical output:
232  if (device != 0x0) {
233  if ((device[0] == '/') && (device[1] == 'n')) return ;
234  }
235 
236  const int npt = 400 ; // Number of points along the axis
237  double rr[npt] ;
238 
239  float* uutab = new float[npt*nprof] ; // Value of uu at the npt points
240  // for each of the nprof profiles
241 
242  double hr = (r_max - r_min) / double(npt-1) ;
243 
244  for (int i=0; i<npt; i++) {
245  rr[i] = hr * i + r_min ;
246  }
247 
248 
249  for (int j=0; j<nprof; j++) {
250 
251  const Scalar& vv = *(uu[j]) ;
252 
253  for (int i=0; i<npt; i++) {
254  uutab[j*npt+i] = float(vv.val_point(rr[i], theta[j], phi[j])) ;
255  }
256  }
257 
258 
259  float xmin = float(radial_scale * r_min) ;
260  float xmax = float(radial_scale * r_max) ;
261 
262  if (nomx == 0x0) nomx = "r" ;
263 
264  if (nomy == 0x0) nomy = "" ;
265 
266  if (title == 0x0) title = "" ;
267 
268  // Preparations for the drawing of boundaries
269  // ------------------------------------------
270 
271  int nbound_max = 100 * nprof ;
272  float* xbound = new float[nbound_max] ;
273  int nbound = 0 ;
274 
275  if (draw_bound) {
276  const double xi_max = 1. ;
277  for (int j=0; j<nprof; j++) {
278 
279  const Map& mp = uu[j]->get_mp() ;
280  int nz = mp.get_mg()->get_nzone() ;
281  int l_max = (mp.get_mg()->get_type_r(nz-1) == UNSURR) ? nz-2 : nz-1 ;
282 
283  for (int l=0; l<=l_max; l++) {
284 
285  double rb = mp.val_r(l, xi_max, theta[j], phi[j]) ;
286 
287  if ((rb >= r_min) && (rb <= r_max)) {
288  xbound[nbound] = float(rb * radial_scale) ;
289  nbound++ ;
290  if (nbound > nbound_max-1) {
291  cout << "des_profile_mult : nbound too large !" << endl ;
292  abort() ;
293  }
294  }
295  }
296  }
297  }
298 
299  // Call to the low level routine
300  // -----------------------------
301 
302  des_profile_mult(uutab, nprof, npt, xmin, xmax, nomx, nomy, title,
303  line_style, ngraph, closeit, device, nbound, xbound) ;
304 
305 
306  delete [] uutab ;
307  delete [] xbound ;
308 
309 }
310 
311 //******************************************************************************
312 
313 void des_meridian(const Scalar& uu, double r_min, double r_max,
314  const char* nomy, int ngraph, const char* device,
315  bool closeit, bool draw_bound) {
316 
317  // Special case of no graphical output:
318  if (device != 0x0) {
319  if ((device[0] == '/') && (device[1] == 'n')) return ;
320  }
321 
322  const Scalar* des[] = {&uu, &uu, &uu, &uu, &uu} ;
323  double phi1[] = {0., 0., 0., 0.25*M_PI, 0.25*M_PI} ;
324  double theta1[] = {0., 0.25*M_PI, 0.5*M_PI, 0., 0.25*M_PI} ;
325 
326  des_profile_mult(des, 5, r_min, r_max, theta1, phi1, 1., closeit,
327  nomy,
328  "phi=0: th=0, pi/4, pi/2, phi=pi/4: th=0, pi/4",
329  ngraph, 0x0, 0x0, device, draw_bound) ;
330 
331 }
332 
333 //******************************************************************************
334 
335 
336 void des_meridian(const Sym_tensor& hh, double r_min, double r_max,
337  const char* name, int ngraph0, const char* device,
338  bool closeit) {
339 
340  // Special case of no graphical output:
341  if (device != 0x0) {
342  if ((device[0] == '/') && (device[1] == 'n')) return ;
343  }
344 
345  char nomy[80] ;
346 
347  int k = 0 ;
348  for (int i=1; i<=3; i++) {
349  for (int j=i; j<=3; j++) {
350 
351  char nom_i[3] ;
352  sprintf(nom_i, "%d", i) ;
353  char nom_j[3] ;
354  sprintf(nom_j, "%d", j) ;
355  strncpy(nomy, name, 40) ;
356  strcat(nomy, " comp. ") ;
357  strcat(nomy, nom_i) ;
358  strcat(nomy, nom_j) ;
359 
360  des_meridian(hh(i,j), r_min, r_max, nomy, ngraph0+k, device,
361  closeit) ;
362  k++ ;
363 
364  }
365  }
366 
367 }
368 
369 
370 //******************************************************************************
371 // VERSION SCALAR SANS UNITES
372 
373 void des_points(const Scalar& uu,
374  double theta, double phi, const char* nomy, const char* title,
375  bool draw_bound) {
376 
377  const Map& mp = uu.get_mp() ;
378  int nz = mp.get_mg()->get_nzone() ;
379  int nt = mp.get_mg()->get_nt(nz-1) ;
380  int np = mp.get_mg()->get_np(nz-1) ;
381 
382 // const int npt = *(uu.get_mp().get_mg())->get_nzone() ; // Number of points along the axis
383 
384 
385  int npt=0;
386 
387  for(int ii = 0; ii<nz; ii++)
388  npt += (uu.get_mp().get_mg())->get_nr(ii) ;
389 
390  float *uutab = new float[npt] ; // define a dynamic array
391  float *xtab = new float[npt] ; // define a dynamic array
392 
393  Mtbl r = *(mp.r.c);
394 
395  for(int ii = 0; ii<nz; ii++){
396  int nr = (uu.get_mp().get_mg())->get_nr(ii) ;
397  for(int ij=0; ij<nr; ij++){
398  uutab[ii*nr+ij] = float(uu.val_grid_point(ii,np-1,nt-1,ij)) ;
399  xtab[ii*nr+ij] = float(r(ii,np-1,nt-1,ij)) ;
400  }
401  }
402 
403  float xmin = float(totalmin(r)) ;
404  float xmax = float(totalmax(r)) ;
405 
406  const char* nomx = "r" ;
407 
408  if (title == 0x0) {
409  title = "" ;
410  }
411 
412  if (nomy == 0x0) {
413  nomy = "" ;
414  }
415 
416  // Preparations for the drawing of boundaries
417  // ------------------------------------------
418  int l_max = (mp.get_mg()->get_type_r(nz-1) == UNSURR) ? nz-2 : nz-1 ;
419 
420  float* xbound = new float[l_max+1] ;
421  int nbound = 0 ;
422 
423  if (draw_bound) {
424  const double xi_max = 1. ;
425  for (int l=0; l<=l_max; l++) {
426 
427  double rb = mp.val_r(l, xi_max, theta, phi) ;
428 
429  if ((rb >= xmin) && (rb <= xmax)) {
430  xbound[nbound] = float(rb) ;
431  nbound++ ;
432  }
433  }
434  }
435 
436  des_profile(uutab, npt, xtab, nomx, nomy, title, 0x0,
437  nbound, xbound) ;
438 
439  delete [] xbound ;
440 
441 }
442 
443 //******************************************************************************
444 
445 void des_points(const Scalar& uu, double scale,
446  double theta, double phi, const char* nomx, const char* nomy,
447  const char* title, bool draw_bound) {
448 
449  const Map& mp = uu.get_mp() ;
450  int nz = mp.get_mg()->get_nzone() ;
451  int nt = mp.get_mg()->get_nt(nz-1) ;
452  int np = mp.get_mg()->get_np(nz-1) ;
453 
454 // const int npt = *(uu.get_mp().get_mg())->get_nzone() ; // Number of points along the axis
455 
456 
457  int npt=0;
458 
459  for(int ii = 0; ii<nz; ii++)
460  npt += (uu.get_mp().get_mg())->get_nr(ii) ;
461 
462  float *uutab = new float[npt] ; // define a dynamic array
463  float *xtab = new float[npt] ; // define a dynamic array
464 
465  Mtbl r = *(mp.r.c);
466 
467  for(int ii = 0; ii<nz; ii++){
468  int nr = (uu.get_mp().get_mg())->get_nr(ii) ;
469  for(int ij=0; ij<nr; ij++){
470  uutab[ii*nr+ij] = float(uu.val_grid_point(ii,np-1,nt-1,ij)) ;
471  xtab[ii*nr+ij] = float(r(ii,np-1,nt-1,ij)) ;
472  }
473  }
474 
475  float xmin = float(totalmin(r) * scale) ;
476  float xmax = float(totalmax(r) * scale) ;
477 
478 
479  if (title == 0x0) {
480  title = "" ;
481  }
482 
483  if (nomx == 0x0) {
484  nomx = "" ;
485  }
486 
487  if (nomy == 0x0) {
488  nomy = "" ;
489  }
490 
491  // Preparations for the drawing of boundaries
492  // ------------------------------------------
493  int l_max = (mp.get_mg()->get_type_r(nz-1) == UNSURR) ? nz-2 : nz-1 ;
494 
495  float* xbound = new float[l_max+1] ;
496  int nbound = 0 ;
497 
498  if (draw_bound) {
499  const double xi_max = 1. ;
500  for (int l=0; l<=l_max; l++) {
501 
502  double rb = mp.val_r(l, xi_max, theta, phi) ;
503 
504  if ((rb >= xmin/scale) && (rb <= xmax/scale)) {
505  xbound[nbound] = float(rb) ;
506  nbound++ ;
507  }
508  }
509  }
510 
511  // Call to the low level routine
512  // -----------------------------
513  des_profile(uutab, npt, xtab, nomx, nomy, title, 0x0,
514  nbound, xbound) ;
515 
516  delete [] xbound ;
517 
518 }
519 }
des_profile
void des_profile(const float *uutab, int nx, float xmin, float xmax, const char *nomx, const char *nomy, const char *title, const char *device=0x0, int nbound=0, float *xbound=0x0)
Basic routine for drawing a single profile with uniform x sampling.
Definition: des_profile.C:91
Lorene::totalmin
double totalmin(const Mtbl &)
Minimum value of the Mtbl elements in all domain.
Definition: mtbl_math.C:522
Lorene
Lorene prototypes.
Definition: app_hor.h:64
des_meridian
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and .
des_profile_mult
void des_profile_mult(const float *uutab, int nprof, int nx, float xmin, float xmax, const char *nomx, const char *nomy, const char *title, const int *line_style, int ngraph, bool closeit, const char *device=0x0, int nbound=0, float *xbound=0x0)
Basic routine for drawing multiple profiles with uniform x sampling.
Definition: des_profile.C:183
Lorene::totalmax
double totalmax(const Mtbl &)
Maximum value of the Mtbl elements in all domains.
Definition: mtbl_math.C:494
des_points
void des_points(const float *uutab, int nx, float xmin, float xmax, const char *nomx=0x0, const char *nomy=0x0, const char *title=0x0, const char *device=0x0, int nbound=0, float *xbound=0x0)
Basic routine for plotting points using grid locations.