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