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 char 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
83 namespace Lorene {
84 double fonc_des_domaine_x(double, const Param&) ;
85 double fonc_des_domaine_y(double, const Param&) ;
86 double fonc_des_domaine_z(double, const Param&) ;
87 
88 //******************************************************************************
89 
90 void 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 
209 void 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 
327 void 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 
448 double 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 
469 double 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 
490 double 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 }
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.
Lorene prototypes.
Definition: app_hor.h:64
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.
Standard units of space, time and mass.
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
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
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.
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:61
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69