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