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
27char 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
63namespace Lorene {
64double fonc_des_surf_scal_x(double, const Param&) ;
65double fonc_des_surf_scal_y(double, const Param&) ;
66double fonc_des_surf_scal_z(double, const Param&) ;
67
68//******************************************************************************
69
70void 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
190void 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
310void 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
439double 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
460double 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
481double 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}
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
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
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.
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.
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.
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.
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
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.