LORENE
des_surface.C
1 /*
2  * Basic routine for drawing the surface of a star
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_surface_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_surface.C,v 1.5 2014/10/13 08:53:23 j_novak Exp $" ;
28 
29 /*
30  * $Id: des_surface.C,v 1.5 2014/10/13 08:53:23 j_novak Exp $
31  * $Log: des_surface.C,v $
32  * Revision 1.5 2014/10/13 08:53:23 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.4 2014/10/06 15:16:05 j_novak
36  * Modified #include directives to use c++ syntax.
37  *
38  * Revision 1.3 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.2 2004/03/25 10:29:25 j_novak
43  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
44  *
45  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
46  * LORENE
47  *
48  * Revision 1.3 2000/03/21 11:14:45 eric
49  * Le parametre precis est mis a 1e-8.
50  *
51  * Revision 1.2 2000/02/11 16:54:00 eric
52  * Utilisation des coordonnees cartesiennes abolues (X,Y,Z) et non plus
53  * des coordonnees relatives (x,y,z).
54  *
55  * Revision 1.1 1999/12/24 13:01:21 eric
56  * Initial revision
57  *
58  *
59  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_surface.C,v 1.5 2014/10/13 08:53:23 j_novak Exp $
60  *
61  */
62 
63 // C headers:
64 #include <cmath>
65 
66 // PGPLOT headers:
67 #include <cpgplot.h>
68 
69 // Lorene headers
70 #include "cmp.h"
71 #include "param.h"
72 #include "utilitaires.h"
73 #include "unites.h"
74 
75 // Local prototypes
76 namespace Lorene {
77 double fonc_des_surface_x(double, const Param&) ;
78 double fonc_des_surface_y(double, const Param&) ;
79 double fonc_des_surface_z(double, const Param&) ;
80 
81 //******************************************************************************
82 
83 void des_surface_x(const Cmp& defsurf, double x0, const char* device, int newgraph,
84  double y_min, double y_max, double z_min, double z_max,
85  const char* nomy, const char* nomz, const char* title, int nxpage, int nypage)
86 {
87  using namespace Unites ;
88 
89  assert(defsurf.get_etat() == ETATQCQ) ;
90 
91 
92  const Map* mp = defsurf.get_mp();
93 
94  double khi ;
95 
96  Param parzerosec ;
97  parzerosec.add_double_mod(x0, 0) ;
98  parzerosec.add_double_mod(khi, 1) ;
99  parzerosec.add_cmp(defsurf) ;
100 
101  double rhomin = 0 ;
102  double rhomax = 2 *
103  mp->val_r(mp->get_mg()->get_nzone() - 1, -1., 0., 0.) ;
104  double precis = 1.e-8 ;
105  int nitermax = 100 ;
106  int niter ;
107 
108  const int np = 101 ;
109  float yg[np] ;
110  float zg[np] ;
111 
112  double hkhi = 2 * M_PI / (np-1) ;
113 
114  bool coupe_surface = true ;
115 
116  for (int i=0; i< np; i++) {
117 
118  khi = hkhi * i ;
119 
120  // Search for the interval [rhomin0, rhomax0] which contains
121  // the first zero of des_surf:
122 
123  double rhomin0 ;
124  double rhomax0 ;
125 
126  if ( zero_premier(fonc_des_surface_x, parzerosec, rhomin, rhomax, 100,
127  rhomin0, rhomax0) == false ) {
128  cout <<
129  "des_surface_x : WARNING : no interval containing a zero of defsurf"
130  << endl ;
131  cout << " has been found for khi = " << khi << " !" << endl ;
132 
133  coupe_surface = false ;
134  break ;
135 
136  }
137 
138 
139  // Search for the zero in the interval [rhomin0, rhomax0] :
140 
141  double rho = zerosec(fonc_des_surface_x, parzerosec, rhomin0, rhomax0,
142  precis, nitermax, niter) ;
143 
144  yg[i] = float(( rho * cos(khi) + mp->get_ori_y() ) / km) ;
145  zg[i] = float(( rho * sin(khi) + mp->get_ori_z() ) / km) ;
146  }
147 
148  // Graphics display
149  // ----------------
150 
151  if ( (newgraph == 1) || (newgraph == 3) ) {
152 
153  if (device == 0x0) device = "?" ;
154 
155  int ier = cpgbeg(0, device, nxpage, nypage) ;
156  if (ier != 1) {
157  cout << "des_surface_x: problem in opening PGPLOT display !" << endl ;
158  }
159 
160  // Taille des caracteres:
161  float size = float(1.3) ;
162  cpgsch(size) ;
163 
164  // Epaisseur des traits:
165  int lepais = 1 ;
166  cpgslw(lepais) ;
167 
168  cpgscf(2) ; // Fonte axes: caracteres romains
169 
170  float ymin1 = float(y_min / km) ;
171  float ymax1 = float(y_max / km) ;
172  float zmin1 = float(z_min / km) ;
173  float zmax1 = float(z_max / km) ;
174 
175  cpgenv(ymin1, ymax1, zmin1, zmax1, 1, 0 ) ;
176 
177  if (nomy == 0x0) nomy = "y [km]" ;
178  if (nomz == 0x0) nomz = "z [km]" ;
179  if (title == 0x0) title = " " ;
180  cpglab(nomy,nomz,title) ;
181 
182  }
183 
184  if (coupe_surface) {
185  cpgsls(1) ; // lignes en trait plein
186  cpgslw(6) ; // traits gras
187  cpgline(np, yg, zg) ;
188  cpgslw(1) ; // traits normaux
189  }
190 
191 
192  // Closing graphic display
193  // -----------------------
194 
195  if ( (newgraph == 2) || (newgraph == 3) ) {
196  cpgend() ;
197  }
198 
199 }
200 
201 //******************************************************************************
202 
203 void des_surface_y(const Cmp& defsurf, double y0, const char* device, int newgraph,
204  double x_min, double x_max, double z_min, double z_max,
205  const char* nomx, const char* nomz, const char* title, int nxpage, int nypage)
206 {
207  using namespace Unites ;
208 
209  assert(defsurf.get_etat() == ETATQCQ) ;
210 
211 
212  const Map* mp = defsurf.get_mp();
213 
214  double khi ;
215 
216  Param parzerosec ;
217  parzerosec.add_double_mod(y0, 0) ;
218  parzerosec.add_double_mod(khi, 1) ;
219  parzerosec.add_cmp(defsurf) ;
220 
221  double rhomin = 0 ;
222  double rhomax = 2 *
223  mp->val_r(mp->get_mg()->get_nzone() - 1, -1., 0., 0.) ;
224  double precis = 1.e-8 ;
225  int nitermax = 100 ;
226  int niter ;
227 
228  const int np = 101 ;
229  float xg[np] ;
230  float zg[np] ;
231 
232  double hkhi = 2 * M_PI / (np-1) ;
233 
234  bool coupe_surface = true ;
235 
236  for (int i=0; i< np; i++) {
237 
238  khi = hkhi * i ;
239 
240  // Search for the interval [rhomin0, rhomax0] which contains
241  // the first zero of des_surf:
242 
243  double rhomin0 ;
244  double rhomax0 ;
245 
246  if ( zero_premier(fonc_des_surface_y, parzerosec, rhomin, rhomax, 100,
247  rhomin0, rhomax0) == false ) {
248  cout <<
249  "des_surface_y : WARNING : no interval containing a zero of defsurf"
250  << endl ;
251  cout << " has been found for khi = " << khi << " !" << endl ;
252 
253  coupe_surface = false ;
254  break ;
255 
256  }
257 
258 
259  // Search for the zero in the interval [rhomin0, rhomax0] :
260 
261  double rho = zerosec(fonc_des_surface_y, parzerosec, rhomin0, rhomax0,
262  precis, nitermax, niter) ;
263 
264  xg[i] = float(( rho * cos(khi) + mp->get_ori_x() ) / km) ;
265  zg[i] = float(( rho * sin(khi) + mp->get_ori_z() ) / km) ;
266  }
267 
268  // Graphics display
269  // ----------------
270 
271  if ( (newgraph == 1) || (newgraph == 3) ) {
272 
273  if (device == 0x0) device = "?" ;
274 
275  int ier = cpgbeg(0, device, nxpage, nypage) ;
276  if (ier != 1) {
277  cout << "des_surface_y: problem in opening PGPLOT display !" << endl ;
278  }
279 
280  // Taille des caracteres:
281  float size = float(1.3) ;
282  cpgsch(size) ;
283 
284  // Epaisseur des traits:
285  int lepais = 1 ;
286  cpgslw(lepais) ;
287 
288  cpgscf(2) ; // Fonte axes: caracteres romains
289 
290  float xmin1 = float(x_min / km) ;
291  float xmax1 = float(x_max / km) ;
292  float zmin1 = float(z_min / km) ;
293  float zmax1 = float(z_max / km) ;
294 
295  cpgenv(xmin1, xmax1, zmin1, zmax1, 1, 0 ) ;
296 
297  if (nomx == 0x0) nomx = "x [km]" ;
298  if (nomz == 0x0) nomz = "z [km]" ;
299  if (title == 0x0) title = " " ;
300  cpglab(nomx,nomz,title) ;
301 
302  }
303 
304  if (coupe_surface) {
305  cpgsls(1) ; // lignes en trait plein
306  cpgslw(6) ; // traits gras
307  cpgline(np, xg, zg) ;
308  cpgslw(1) ; // traits normaux
309  }
310 
311 
312  // Closing graphic display
313  // -----------------------
314 
315  if ( (newgraph == 2) || (newgraph == 3) ) {
316  cpgend() ;
317  }
318 
319 }
320 
321 //******************************************************************************
322 
323 void des_surface_z(const Cmp& defsurf, double z0, const char* device, int newgraph,
324  double x_min, double x_max, double y_min, double y_max,
325  const char* nomx, const char* nomy, const char* title, int nxpage, int nypage)
326 {
327  using namespace Unites ;
328 
329  assert(defsurf.get_etat() == ETATQCQ) ;
330 
331 
332  const Map* mp = defsurf.get_mp();
333 
334  double khi ;
335 
336  Param parzerosec ;
337  parzerosec.add_double_mod(z0, 0) ;
338  parzerosec.add_double_mod(khi, 1) ;
339  parzerosec.add_cmp(defsurf) ;
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-8 ;
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_surface_z, parzerosec, rhomin, rhomax, 100,
367  rhomin0, rhomax0) == false ) {
368  cout <<
369  "des_surface_z : WARNING : no interval containing a zero of defsurf"
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_surface_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  // Graphics display
389  // ----------------
390 
391  if ( (newgraph == 1) || (newgraph == 3) ) {
392 
393  if (device == 0x0) device = "?" ;
394 
395  int ier = cpgbeg(0, device, nxpage, nypage) ;
396  if (ier != 1) {
397  cout << "des_surface_z: problem in opening PGPLOT display !" << endl ;
398  }
399 
400  // Taille des caracteres:
401  float size = float(1.3) ;
402  cpgsch(size) ;
403 
404  // Epaisseur des traits:
405  int lepais = 1 ;
406  cpgslw(lepais) ;
407 
408  cpgscf(2) ; // Fonte axes: caracteres romains
409 
410  float xmin1 = float(x_min / km) ;
411  float xmax1 = float(x_max / km) ;
412  float ymin1 = float(y_min / km) ;
413  float ymax1 = float(y_max / km) ;
414 
415  cpgenv(xmin1, xmax1, ymin1, ymax1, 1, 0 ) ;
416 
417  if (nomx == 0x0) nomx = "x [km]" ;
418  if (nomy == 0x0) nomy = "y [km]" ;
419  if (title == 0x0) title = " " ;
420  cpglab(nomx,nomy,title) ;
421 
422  }
423 
424  if (coupe_surface) {
425  cpgsls(1) ; // lignes en trait plein
426  cpgslw(6) ; // traits gras
427  cpgline(np, xg, yg) ;
428  cpgslw(1) ; // traits normaux
429  }
430 
431 
432  // Closing graphic display
433  // -----------------------
434 
435  if ( (newgraph == 2) || (newgraph == 3) ) {
436  cpgend() ;
437  }
438 
439 }
440 
441 
442 
443 
444 
445 
446 
447 
448 
449 
450 //*****************************************************************************
451 
452 double fonc_des_surface_x(double vrho, const Param& par) {
453 
454  double x = par.get_double_mod(0) ;
455  double khi = par.get_double_mod(1) ;
456  const Cmp& defsurf = par.get_cmp() ;
457  const Map& mp = *(defsurf.get_mp()) ;
458 
459  // Absolute Cartesian coordinates:
460  double y = vrho * cos(khi) + mp.get_ori_y() ;
461  double z = vrho * sin(khi) + mp.get_ori_z() ;
462 
463  // Spherical coordinates of the mapping:
464  double r, theta, phi ;
465  mp.convert_absolute(x, y, z, r, theta, phi) ;
466 
467  return defsurf.val_point(r, theta, phi) ;
468 
469 }
470 
471 //*****************************************************************************
472 
473 double fonc_des_surface_y(double vrho, const Param& par) {
474 
475  double y = par.get_double_mod(0) ;
476  double khi = par.get_double_mod(1) ;
477  const Cmp& defsurf = par.get_cmp() ;
478  const Map& mp = *(defsurf.get_mp()) ;
479 
480  // Absolute Cartesian coordinates:
481  double x = vrho * cos(khi) + mp.get_ori_x() ;
482  double z = vrho * sin(khi) + mp.get_ori_z() ;
483 
484  // Spherical coordinates of the mapping:
485  double r, theta, phi ;
486  mp.convert_absolute(x, y, z, r, theta, phi) ;
487 
488  return defsurf.val_point(r, theta, phi) ;
489 
490 }
491 
492 //*****************************************************************************
493 
494 double fonc_des_surface_z(double vrho, const Param& par) {
495 
496  double z = par.get_double_mod(0) ;
497  double khi = par.get_double_mod(1) ;
498  const Cmp& defsurf = par.get_cmp() ;
499  const Map& mp = *(defsurf.get_mp()) ;
500 
501  // Absolute Cartesian coordinates:
502  double x = vrho * cos(khi) + mp.get_ori_x() ;
503  double y = vrho * sin(khi) + mp.get_ori_y() ;
504 
505  // Spherical coordinates of the mapping:
506  double r, theta, phi ;
507  mp.convert_absolute(x, y, z, r, theta, phi) ;
508 
509  return defsurf.val_point(r, theta, phi) ;
510 
511 }
512 }
des_surface_y
void des_surface_y(const Scalar &defsurf, 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 a stellar surface in a plane Y=constant.
Lorene
Lorene prototypes.
Definition: app_hor.h:64
des_surface_z
void des_surface_z(const Scalar &defsurf, 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 a stellar surface in a plane Z=constant.
Lorene::Cmp::val_point
double val_point(double r, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point , by means of the spectral...
Definition: cmp.C:732
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:89
Unites
Standard units of space, time and mass.
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:61
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
des_surface_x
void des_surface_x(const Scalar &defsurf, 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 a stellar surface in a plane X=constant.
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69