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
28
29/*
30 * $Id: des_surface.C,v 1.6 2016/12/05 16:18:06 j_novak Exp $
31 * $Log: des_surface.C,v $
32 * Revision 1.6 2016/12/05 16:18:06 j_novak
33 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
34 *
35 * Revision 1.5 2014/10/13 08:53:23 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.4 2014/10/06 15:16:05 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.3 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.2 2004/03/25 10:29:25 j_novak
46 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
47 *
48 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
49 * LORENE
50 *
51 * Revision 1.3 2000/03/21 11:14:45 eric
52 * Le parametre precis est mis a 1e-8.
53 *
54 * Revision 1.2 2000/02/11 16:54:00 eric
55 * Utilisation des coordonnees cartesiennes abolues (X,Y,Z) et non plus
56 * des coordonnees relatives (x,y,z).
57 *
58 * Revision 1.1 1999/12/24 13:01:21 eric
59 * Initial revision
60 *
61 *
62 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_surface.C,v 1.6 2016/12/05 16:18:06 j_novak Exp $
63 *
64 */
65
66// C headers:
67#include <cmath>
68
69// PGPLOT headers:
70#include <cpgplot.h>
71
72// Lorene headers
73#include "cmp.h"
74#include "param.h"
75#include "utilitaires.h"
76#include "unites.h"
77
78// Local prototypes
79namespace Lorene {
80double fonc_des_surface_x(double, const Param&) ;
81double fonc_des_surface_y(double, const Param&) ;
82double fonc_des_surface_z(double, const Param&) ;
83
84//******************************************************************************
85
86void des_surface_x(const Cmp& defsurf, double x0, const char* device, int newgraph,
87 double y_min, double y_max, double z_min, double z_max,
88 const char* nomy, const char* nomz, const char* title, int nxpage, int nypage)
89{
90 using namespace Unites ;
91
92 assert(defsurf.get_etat() == ETATQCQ) ;
93
94
95 const Map* mp = defsurf.get_mp();
96
97 double khi ;
98
99 Param parzerosec ;
100 parzerosec.add_double_mod(x0, 0) ;
101 parzerosec.add_double_mod(khi, 1) ;
102 parzerosec.add_cmp(defsurf) ;
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-8 ;
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_surface_x, parzerosec, rhomin, rhomax, 100,
130 rhomin0, rhomax0) == false ) {
131 cout <<
132 "des_surface_x : WARNING : no interval containing a zero of defsurf"
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_surface_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 // Graphics display
152 // ----------------
153
154 if ( (newgraph == 1) || (newgraph == 3) ) {
155
156 if (device == 0x0) device = "?" ;
157
158 int ier = cpgbeg(0, device, nxpage, nypage) ;
159 if (ier != 1) {
160 cout << "des_surface_x: problem in opening PGPLOT display !" << endl ;
161 }
162
163 // Taille des caracteres:
164 float size = float(1.3) ;
165 cpgsch(size) ;
166
167 // Epaisseur des traits:
168 int lepais = 1 ;
169 cpgslw(lepais) ;
170
171 cpgscf(2) ; // Fonte axes: caracteres romains
172
173 float ymin1 = float(y_min / km) ;
174 float ymax1 = float(y_max / km) ;
175 float zmin1 = float(z_min / km) ;
176 float zmax1 = float(z_max / km) ;
177
178 cpgenv(ymin1, ymax1, zmin1, zmax1, 1, 0 ) ;
179
180 if (nomy == 0x0) nomy = "y [km]" ;
181 if (nomz == 0x0) nomz = "z [km]" ;
182 if (title == 0x0) title = " " ;
183 cpglab(nomy,nomz,title) ;
184
185 }
186
187 if (coupe_surface) {
188 cpgsls(1) ; // lignes en trait plein
189 cpgslw(6) ; // traits gras
190 cpgline(np, yg, zg) ;
191 cpgslw(1) ; // traits normaux
192 }
193
194
195 // Closing graphic display
196 // -----------------------
197
198 if ( (newgraph == 2) || (newgraph == 3) ) {
199 cpgend() ;
200 }
201
202}
203
204//******************************************************************************
205
206void des_surface_y(const Cmp& defsurf, double y0, const char* device, int newgraph,
207 double x_min, double x_max, double z_min, double z_max,
208 const char* nomx, const char* nomz, const char* title, int nxpage, int nypage)
209{
210 using namespace Unites ;
211
212 assert(defsurf.get_etat() == ETATQCQ) ;
213
214
215 const Map* mp = defsurf.get_mp();
216
217 double khi ;
218
219 Param parzerosec ;
220 parzerosec.add_double_mod(y0, 0) ;
221 parzerosec.add_double_mod(khi, 1) ;
222 parzerosec.add_cmp(defsurf) ;
223
224 double rhomin = 0 ;
225 double rhomax = 2 *
226 mp->val_r(mp->get_mg()->get_nzone() - 1, -1., 0., 0.) ;
227 double precis = 1.e-8 ;
228 int nitermax = 100 ;
229 int niter ;
230
231 const int np = 101 ;
232 float xg[np] ;
233 float zg[np] ;
234
235 double hkhi = 2 * M_PI / (np-1) ;
236
237 bool coupe_surface = true ;
238
239 for (int i=0; i< np; i++) {
240
241 khi = hkhi * i ;
242
243 // Search for the interval [rhomin0, rhomax0] which contains
244 // the first zero of des_surf:
245
246 double rhomin0 ;
247 double rhomax0 ;
248
249 if ( zero_premier(fonc_des_surface_y, parzerosec, rhomin, rhomax, 100,
250 rhomin0, rhomax0) == false ) {
251 cout <<
252 "des_surface_y : WARNING : no interval containing a zero of defsurf"
253 << endl ;
254 cout << " has been found for khi = " << khi << " !" << endl ;
255
256 coupe_surface = false ;
257 break ;
258
259 }
260
261
262 // Search for the zero in the interval [rhomin0, rhomax0] :
263
264 double rho = zerosec(fonc_des_surface_y, parzerosec, rhomin0, rhomax0,
265 precis, nitermax, niter) ;
266
267 xg[i] = float(( rho * cos(khi) + mp->get_ori_x() ) / km) ;
268 zg[i] = float(( rho * sin(khi) + mp->get_ori_z() ) / km) ;
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_surface_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(1) ; // lignes en trait plein
309 cpgslw(6) ; // traits gras
310 cpgline(np, xg, zg) ;
311 cpgslw(1) ; // traits normaux
312 }
313
314
315 // Closing graphic display
316 // -----------------------
317
318 if ( (newgraph == 2) || (newgraph == 3) ) {
319 cpgend() ;
320 }
321
322}
323
324//******************************************************************************
325
326void des_surface_z(const Cmp& defsurf, double z0, const char* device, int newgraph,
327 double x_min, double x_max, double y_min, double y_max,
328 const char* nomx, const char* nomy, const char* title, int nxpage, int nypage)
329{
330 using namespace Unites ;
331
332 assert(defsurf.get_etat() == ETATQCQ) ;
333
334
335 const Map* mp = defsurf.get_mp();
336
337 double khi ;
338
339 Param parzerosec ;
340 parzerosec.add_double_mod(z0, 0) ;
341 parzerosec.add_double_mod(khi, 1) ;
342 parzerosec.add_cmp(defsurf) ;
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-8 ;
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_surface_z, parzerosec, rhomin, rhomax, 100,
370 rhomin0, rhomax0) == false ) {
371 cout <<
372 "des_surface_z : WARNING : no interval containing a zero of defsurf"
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_surface_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 // Graphics display
392 // ----------------
393
394 if ( (newgraph == 1) || (newgraph == 3) ) {
395
396 if (device == 0x0) device = "?" ;
397
398 int ier = cpgbeg(0, device, nxpage, nypage) ;
399 if (ier != 1) {
400 cout << "des_surface_z: problem in opening PGPLOT display !" << endl ;
401 }
402
403 // Taille des caracteres:
404 float size = float(1.3) ;
405 cpgsch(size) ;
406
407 // Epaisseur des traits:
408 int lepais = 1 ;
409 cpgslw(lepais) ;
410
411 cpgscf(2) ; // Fonte axes: caracteres romains
412
413 float xmin1 = float(x_min / km) ;
414 float xmax1 = float(x_max / km) ;
415 float ymin1 = float(y_min / km) ;
416 float ymax1 = float(y_max / km) ;
417
418 cpgenv(xmin1, xmax1, ymin1, ymax1, 1, 0 ) ;
419
420 if (nomx == 0x0) nomx = "x [km]" ;
421 if (nomy == 0x0) nomy = "y [km]" ;
422 if (title == 0x0) title = " " ;
423 cpglab(nomx,nomy,title) ;
424
425 }
426
427 if (coupe_surface) {
428 cpgsls(1) ; // lignes en trait plein
429 cpgslw(6) ; // traits gras
430 cpgline(np, xg, yg) ;
431 cpgslw(1) ; // traits normaux
432 }
433
434
435 // Closing graphic display
436 // -----------------------
437
438 if ( (newgraph == 2) || (newgraph == 3) ) {
439 cpgend() ;
440 }
441
442}
443
444
445
446
447
448
449
450
451
452
453//*****************************************************************************
454
455double fonc_des_surface_x(double vrho, const Param& par) {
456
457 double x = par.get_double_mod(0) ;
458 double khi = par.get_double_mod(1) ;
459 const Cmp& defsurf = par.get_cmp() ;
460 const Map& mp = *(defsurf.get_mp()) ;
461
462 // Absolute Cartesian coordinates:
463 double y = vrho * cos(khi) + mp.get_ori_y() ;
464 double z = vrho * sin(khi) + mp.get_ori_z() ;
465
466 // Spherical coordinates of the mapping:
467 double r, theta, phi ;
468 mp.convert_absolute(x, y, z, r, theta, phi) ;
469
470 return defsurf.val_point(r, theta, phi) ;
471
472}
473
474//*****************************************************************************
475
476double fonc_des_surface_y(double vrho, const Param& par) {
477
478 double y = par.get_double_mod(0) ;
479 double khi = par.get_double_mod(1) ;
480 const Cmp& defsurf = par.get_cmp() ;
481 const Map& mp = *(defsurf.get_mp()) ;
482
483 // Absolute Cartesian coordinates:
484 double x = vrho * cos(khi) + mp.get_ori_x() ;
485 double z = vrho * sin(khi) + mp.get_ori_z() ;
486
487 // Spherical coordinates of the mapping:
488 double r, theta, phi ;
489 mp.convert_absolute(x, y, z, r, theta, phi) ;
490
491 return defsurf.val_point(r, theta, phi) ;
492
493}
494
495//*****************************************************************************
496
497double fonc_des_surface_z(double vrho, const Param& par) {
498
499 double z = par.get_double_mod(0) ;
500 double khi = par.get_double_mod(1) ;
501 const Cmp& defsurf = par.get_cmp() ;
502 const Map& mp = *(defsurf.get_mp()) ;
503
504 // Absolute Cartesian coordinates:
505 double x = vrho * cos(khi) + mp.get_ori_x() ;
506 double y = vrho * sin(khi) + mp.get_ori_y() ;
507
508 // Spherical coordinates of the mapping:
509 double r, theta, phi ;
510 mp.convert_absolute(x, y, z, r, theta, phi) ;
511
512 return defsurf.val_point(r, theta, phi) ;
513
514}
515}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_double_mod
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition param.C:456
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