LORENE
des_domaine.C
1/*
2 * Basic routines for drawing the boundaries between domains.
3 */
4
5/*
6 * Copyright (c) 1999-2001 Eric Gourgoulhon
7 *
8 * This file is part of LORENE.
9 *
10 * LORENE is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * LORENE is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with LORENE; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
26
27char des_domaine_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_domaine.C,v 1.6 2014/10/13 08:53:22 j_novak Exp $" ;
28
29/*
30 * $Id: des_domaine.C,v 1.6 2014/10/13 08:53:22 j_novak Exp $
31 * $Log: des_domaine.C,v $
32 * Revision 1.6 2014/10/13 08:53:22 j_novak
33 * Lorene classes and functions now belong to the namespace Lorene.
34 *
35 * Revision 1.5 2014/10/06 15:16:05 j_novak
36 * Modified #include directives to use c++ syntax.
37 *
38 * Revision 1.4 2008/08/19 06:42:00 j_novak
39 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
40 * cast-type operations, and constant strings that must be defined as const char*
41 *
42 * Revision 1.3 2005/03/24 14:33:03 e_gourgoulhon
43 * Ponderation of precis by rhomax-rhomax (to draw domains with
44 * large r).
45 *
46 * Revision 1.2 2004/03/25 10:29:25 j_novak
47 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
48 *
49 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
50 * LORENE
51 *
52 * Revision 1.4 2001/02/28 09:45:02 eric
53 * Correction erreur affichage "des_domaine_y" dans des_domaine_z.
54 *
55 * Revision 1.3 2000/02/11 16:53:50 eric
56 * Utilisation des coordonnees cartesiennes abolues (X,Y,Z) et non plus
57 * des coordonnees relatives (x,y,z).
58 *
59 * Revision 1.2 1999/12/27 12:20:59 eric
60 * *** empty log message ***
61 *
62 * Revision 1.1 1999/12/27 12:19:03 eric
63 * Initial revision
64 *
65 *
66 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_domaine.C,v 1.6 2014/10/13 08:53:22 j_novak Exp $
67 *
68 */
69
70// C headers:
71#include <cmath>
72
73// PGPLOT headers:
74#include <cpgplot.h>
75
76// Lorene headers
77#include "map.h"
78#include "param.h"
79#include "utilitaires.h"
80#include "unites.h"
81
82// Local prototypes
83namespace Lorene {
84double fonc_des_domaine_x(double, const Param&) ;
85double fonc_des_domaine_y(double, const Param&) ;
86double fonc_des_domaine_z(double, const Param&) ;
87
88//******************************************************************************
89
90void des_domaine_x(const Map& mp, int l0, double x0, const char* device, int newgraph,
91 double y_min, double y_max, double z_min, double z_max,
92 const char* nomy, const char* nomz, const char* title, int nxpage, int nypage)
93{
94 using namespace Unites ;
95
96 double khi ;
97
98 Param parzerosec ;
99 parzerosec.add_int(l0, 0) ;
100 parzerosec.add_double_mod(x0, 0) ;
101 parzerosec.add_double_mod(khi, 1) ;
102 parzerosec.add_map(mp, 0) ;
103
104 double rhomin = 0 ;
105 double rhomax = 2 *
106 mp.val_r(mp.get_mg()->get_nzone() - 1, -1., 0., 0.) ;
107 double precis = 1.e-14 * (rhomax - rhomin) ;
108 int nitermax = 100 ;
109 int niter ;
110
111 const int np = 101 ;
112 float yg[np] ;
113 float zg[np] ;
114
115 double hkhi = 2 * M_PI / (np-1) ;
116
117 bool coupe_surface = true ;
118
119 for (int i=0; i< np; i++) {
120
121 khi = hkhi * i ;
122
123 // Search for the interval [rhomin0, rhomax0] which contains
124 // the first zero of des_surf:
125
126 double rhomin0 ;
127 double rhomax0 ;
128
129 if ( zero_premier(fonc_des_domaine_x, parzerosec, rhomin, rhomax, 100,
130 rhomin0, rhomax0) == false ) {
131 cout <<
132 "des_domaine_x : WARNING : no crossing with the domain boundary"
133 << endl ;
134 cout << " has been found for khi = " << khi << " !" << endl ;
135
136 coupe_surface = false ;
137 break ;
138
139 }
140
141
142 // Search for the zero in the interval [rhomin0, rhomax0] :
143
144 double rho = zerosec(fonc_des_domaine_x, parzerosec, rhomin0, rhomax0,
145 precis, nitermax, niter) ;
146
147 yg[i] = float(( rho * cos(khi) + mp.get_ori_y() ) / km) ;
148 zg[i] = float(( rho * sin(khi) + mp.get_ori_z() ) / km) ;
149
150 }
151
152 // Graphics display
153 // ----------------
154
155 if ( (newgraph == 1) || (newgraph == 3) ) {
156
157 if (device == 0x0) device = "?" ;
158
159 int ier = cpgbeg(0, device, nxpage, nypage) ;
160 if (ier != 1) {
161 cout << "des_domaine_x: problem in opening PGPLOT display !" << endl ;
162 }
163
164 // Taille des caracteres:
165 float size = float(1.3) ;
166 cpgsch(size) ;
167
168 // Epaisseur des traits:
169 int lepais = 1 ;
170 cpgslw(lepais) ;
171
172 cpgscf(2) ; // Fonte axes: caracteres romains
173
174 float ymin1 = float(y_min / km) ;
175 float ymax1 = float(y_max / km) ;
176 float zmin1 = float(z_min / km) ;
177 float zmax1 = float(z_max / km) ;
178
179 cpgenv(ymin1, ymax1, zmin1, zmax1, 1, 0 ) ;
180
181 if (nomy == 0x0) nomy = "y [km]" ;
182 if (nomz == 0x0) nomz = "z [km]" ;
183 if (title == 0x0) title = " " ;
184 cpglab(nomy,nomz,title) ;
185
186 }
187
188 if (coupe_surface) {
189 cpgsls(3) ; // lignes en trait mixte
190 cpgsci(3) ; // couleur verte
191 cpgline(np, yg, zg) ;
192 cpgsls(1) ; // retour aux lignes en trait plein
193 cpgsci(1) ; // couleur noire
194 }
195
196
197 // Closing graphic display
198 // -----------------------
199
200 if ( (newgraph == 2) || (newgraph == 3) ) {
201 cpgend() ;
202 }
203
204}
205
206
207//******************************************************************************
208
209void des_domaine_y(const Map& mp, int l0, double y0, const char* device, int newgraph,
210 double x_min, double x_max, double z_min, double z_max,
211 const char* nomx, const char* nomz, const char* title, int nxpage, int nypage)
212{
213 using namespace Unites ;
214
215 double khi ;
216
217 Param parzerosec ;
218 parzerosec.add_int(l0, 0) ;
219 parzerosec.add_double_mod(y0, 0) ;
220 parzerosec.add_double_mod(khi, 1) ;
221 parzerosec.add_map(mp, 0) ;
222
223 double rhomin = 0 ;
224 double rhomax = 2 *
225 mp.val_r(mp.get_mg()->get_nzone() - 1, -1., 0., 0.) ;
226 double precis = 1.e-14 * (rhomax - rhomin) ;
227 int nitermax = 100 ;
228 int niter ;
229
230 const int np = 101 ;
231 float xg[np] ;
232 float zg[np] ;
233
234 double hkhi = 2 * M_PI / (np-1) ;
235
236 bool coupe_surface = true ;
237
238 for (int i=0; i< np; i++) {
239
240 khi = hkhi * i ;
241
242 // Search for the interval [rhomin0, rhomax0] which contains
243 // the first zero of des_surf:
244
245 double rhomin0 ;
246 double rhomax0 ;
247
248 if ( zero_premier(fonc_des_domaine_y, parzerosec, rhomin, rhomax, 100,
249 rhomin0, rhomax0) == false ) {
250 cout <<
251 "des_domaine_y : WARNING : no crossing with the domain boundary"
252 << endl ;
253 cout << " has been found for khi = " << khi << " !" << endl ;
254
255 coupe_surface = false ;
256 break ;
257
258 }
259
260
261 // Search for the zero in the interval [rhomin0, rhomax0] :
262
263 double rho = zerosec(fonc_des_domaine_y, parzerosec, rhomin0, rhomax0,
264 precis, nitermax, niter) ;
265
266 xg[i] = float(( rho * cos(khi) + mp.get_ori_x() ) / km) ;
267 zg[i] = float(( rho * sin(khi) + mp.get_ori_z() ) / km) ;
268
269 }
270
271 // Graphics display
272 // ----------------
273
274 if ( (newgraph == 1) || (newgraph == 3) ) {
275
276 if (device == 0x0) device = "?" ;
277
278 int ier = cpgbeg(0, device, nxpage, nypage) ;
279 if (ier != 1) {
280 cout << "des_domaine_y: problem in opening PGPLOT display !" << endl ;
281 }
282
283 // Taille des caracteres:
284 float size = float(1.3) ;
285 cpgsch(size) ;
286
287 // Epaisseur des traits:
288 int lepais = 1 ;
289 cpgslw(lepais) ;
290
291 cpgscf(2) ; // Fonte axes: caracteres romains
292
293 float xmin1 = float(x_min / km) ;
294 float xmax1 = float(x_max / km) ;
295 float zmin1 = float(z_min / km) ;
296 float zmax1 = float(z_max / km) ;
297
298 cpgenv(xmin1, xmax1, zmin1, zmax1, 1, 0 ) ;
299
300 if (nomx == 0x0) nomx = "x [km]" ;
301 if (nomz == 0x0) nomz = "z [km]" ;
302 if (title == 0x0) title = " " ;
303 cpglab(nomx,nomz,title) ;
304
305 }
306
307 if (coupe_surface) {
308 cpgsls(3) ; // lignes en trait mixte
309 cpgsci(3) ; // couleur verte
310 cpgline(np, xg, zg) ;
311 cpgsls(1) ; // retour aux lignes en trait plein
312 cpgsci(1) ; // couleur noire
313 }
314
315
316 // Closing graphic display
317 // -----------------------
318
319 if ( (newgraph == 2) || (newgraph == 3) ) {
320 cpgend() ;
321 }
322
323}
324
325//******************************************************************************
326
327void des_domaine_z(const Map& mp, int l0, double z0, const char* device, int newgraph,
328 double x_min, double x_max, double y_min, double y_max,
329 const char* nomx, const char* nomy, const char* title, int nxpage, int nypage)
330{
331 using namespace Unites ;
332
333 double khi ;
334
335 Param parzerosec ;
336 parzerosec.add_int(l0, 0) ;
337 parzerosec.add_double_mod(z0, 0) ;
338 parzerosec.add_double_mod(khi, 1) ;
339 parzerosec.add_map(mp, 0) ;
340
341 double rhomin = 0 ;
342 double rhomax = 2 *
343 mp.val_r(mp.get_mg()->get_nzone() - 1, -1., 0., 0.) ;
344 double precis = 1.e-14 * (rhomax - rhomin) ;
345 int nitermax = 100 ;
346 int niter ;
347
348 const int np = 101 ;
349 float xg[np] ;
350 float yg[np] ;
351
352 double hkhi = 2 * M_PI / (np-1) ;
353
354 bool coupe_surface = true ;
355
356 for (int i=0; i< np; i++) {
357
358 khi = hkhi * i ;
359
360 // Search for the interval [rhomin0, rhomax0] which contains
361 // the first zero of des_surf:
362
363 double rhomin0 ;
364 double rhomax0 ;
365
366 if ( zero_premier(fonc_des_domaine_z, parzerosec, rhomin, rhomax, 100,
367 rhomin0, rhomax0) == false ) {
368 cout <<
369 "des_domaine_z : WARNING : no crossing with the domain boundary"
370 << endl ;
371 cout << " has been found for khi = " << khi << " !" << endl ;
372
373 coupe_surface = false ;
374 break ;
375
376 }
377
378
379 // Search for the zero in the interval [rhomin0, rhomax0] :
380
381 double rho = zerosec(fonc_des_domaine_z, parzerosec, rhomin0, rhomax0,
382 precis, nitermax, niter) ;
383
384 xg[i] = float(( rho * cos(khi) + mp.get_ori_x() ) / km) ;
385 yg[i] = float(( rho * sin(khi) + mp.get_ori_y() ) / km) ;
386
387 }
388
389 // Graphics display
390 // ----------------
391
392 if ( (newgraph == 1) || (newgraph == 3) ) {
393
394 if (device == 0x0) device = "?" ;
395
396 int ier = cpgbeg(0, device, nxpage, nypage) ;
397 if (ier != 1) {
398 cout << "des_domaine_z: problem in opening PGPLOT display !" << endl ;
399 }
400
401 // Taille des caracteres:
402 float size = float(1.3) ;
403 cpgsch(size) ;
404
405 // Epaisseur des traits:
406 int lepais = 1 ;
407 cpgslw(lepais) ;
408
409 cpgscf(2) ; // Fonte axes: caracteres romains
410
411 float xmin1 = float(x_min / km) ;
412 float xmax1 = float(x_max / km) ;
413 float ymin1 = float(y_min / km) ;
414 float ymax1 = float(y_max / km) ;
415
416 cpgenv(xmin1, xmax1, ymin1, ymax1, 1, 0 ) ;
417
418 if (nomx == 0x0) nomx = "x [km]" ;
419 if (nomy == 0x0) nomy = "y [km]" ;
420 if (title == 0x0) title = " " ;
421 cpglab(nomx,nomy,title) ;
422
423 }
424
425 if (coupe_surface) {
426 cpgsls(3) ; // lignes en trait mixte
427 cpgsci(3) ; // couleur verte
428 cpgline(np, xg, yg) ;
429 cpgsls(1) ; // retour aux lignes en trait plein
430 cpgsci(1) ; // couleur noire
431 }
432
433
434 // Closing graphic display
435 // -----------------------
436
437 if ( (newgraph == 2) || (newgraph == 3) ) {
438 cpgend() ;
439 }
440
441}
442
443
444
445
446//*****************************************************************************
447
448double fonc_des_domaine_x(double vrho, const Param& par) {
449
450 int l = par.get_int(0) ;
451 double x = par.get_double_mod(0) ;
452 double khi = par.get_double_mod(1) ;
453 const Map& mp = par.get_map(0) ;
454
455 // Absolute Cartesian coordinates:
456 double y = vrho * cos(khi) + mp.get_ori_y() ;
457 double z = vrho * sin(khi) + mp.get_ori_z() ;
458
459 // Spherical coordinates of the mapping:
460 double r, theta, phi ;
461 mp.convert_absolute(x, y, z, r, theta, phi) ;
462
463 return r - mp.val_r(l, 1., theta, phi) ;
464
465}
466
467//*****************************************************************************
468
469double fonc_des_domaine_y(double vrho, const Param& par) {
470
471 int l = par.get_int(0) ;
472 double y = par.get_double_mod(0) ;
473 double khi = par.get_double_mod(1) ;
474 const Map& mp = par.get_map(0) ;
475
476 // Absolute Cartesian coordinates:
477 double x = vrho * cos(khi) + mp.get_ori_x() ;
478 double z = vrho * sin(khi) + mp.get_ori_z() ;
479
480 // Spherical coordinates of the mapping:
481 double r, theta, phi ;
482 mp.convert_absolute(x, y, z, r, theta, phi) ;
483
484 return r - mp.val_r(l, 1., theta, phi) ;
485
486}
487
488//*****************************************************************************
489
490double fonc_des_domaine_z(double vrho, const Param& par) {
491
492 int l = par.get_int(0) ;
493 double z = par.get_double_mod(0) ;
494 double khi = par.get_double_mod(1) ;
495 const Map& mp = par.get_map(0) ;
496
497 // Absolute Cartesian coordinates:
498 double x = vrho * cos(khi) + mp.get_ori_x() ;
499 double y = vrho * sin(khi) + mp.get_ori_y() ;
500
501 // Spherical coordinates of the mapping:
502 double r, theta, phi ;
503 mp.convert_absolute(x, y, z, r, theta, phi) ;
504
505 return r - mp.val_r(l, 1., theta, phi) ;
506
507}
508
509}
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
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.
bool zero_premier(double(*f)(double, const Param &), const Param &par, double a, double b, int n, double &a0, double &b0)
Locates the sub-interval containing the first zero of a function in a given interval.
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.
Definition zerosec.C:89
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.