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
27char 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
76namespace Lorene {
77double fonc_des_surface_x(double, const Param&) ;
78double fonc_des_surface_y(double, const Param&) ;
79double fonc_des_surface_z(double, const Param&) ;
80
81//******************************************************************************
82
83void 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
203void 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
323void 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
452double 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
473double 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
494double 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}
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.