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