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