LORENE
eos_bf_poly.C
1/*
2 * Methods of the class Eos_bf_poly.
3 *
4 * (see file eos_bifluid.h for documentation).
5 */
6
7/*
8 * Copyright (c) 2001-2002 Jerome Novak
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28
29char eos_bf_poly_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_poly.C,v 1.22 2014/10/13 08:52:52 j_novak Exp $" ;
30
31/*
32 * $Id: eos_bf_poly.C,v 1.22 2014/10/13 08:52:52 j_novak Exp $
33 * $Log: eos_bf_poly.C,v $
34 * Revision 1.22 2014/10/13 08:52:52 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.21 2014/10/06 15:13:06 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.20 2014/04/25 10:43:51 j_novak
41 * The member 'name' is of type string now. Correction of a few const-related issues.
42 *
43 * Revision 1.19 2008/08/19 06:42:00 j_novak
44 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
45 * cast-type operations, and constant strings that must be defined as const char*
46 *
47 * Revision 1.18 2004/05/13 15:27:42 r_prix
48 * fixed a little eos-bug also in the relativistic case (same as done already in Newt)
49 *
50 * Revision 1.17 2003/12/17 23:12:32 r_prix
51 * replaced use of C++ <string> by standard ANSI char* to be backwards compatible
52 * with broken compilers like MIPSpro Compiler 7.2 on SGI Origin200. ;-)
53 *
54 * Revision 1.16 2003/12/10 08:58:20 r_prix
55 * - added new Eos_bifluid paramter for eos-file: bool slow_rot_style
56 * to indicate if we want this particular kind of EOS-inversion (only works for
57 * the Newtonian 'analytic' polytropes) --> replaces former dirty hack with gamma1<0
58 *
59 * Revision 1.15 2003/12/05 15:09:47 r_prix
60 * adapted Eos_bifluid class and subclasses to use read_variable() for
61 * (formatted) file-reading.
62 *
63 * Revision 1.14 2003/12/04 14:24:41 r_prix
64 * a really dirty hack: gam1 < 0 signals to use 'slow-rot-style' EOS
65 * inversion (i.e. typeos=5). This only works for the 'analytic EOS'.
66 *
67 * Revision 1.13 2003/11/18 18:28:38 r_prix
68 * moved particle-masses m_1, m_2 of the two fluids into class eos_bifluid (from eos_bf_poly)
69 *
70 * Revision 1.12 2003/03/17 10:28:04 j_novak
71 * Removed too strong asserts on beta and kappa
72 *
73 * Revision 1.11 2003/02/07 10:47:43 j_novak
74 * The possibility of having gamma5 xor gamma6 =0 has been introduced for
75 * tests. The corresponding parameter files have been added.
76 *
77 * Revision 1.10 2002/10/16 14:36:34 j_novak
78 * Reorganization of #include instructions of standard C++, in order to
79 * use experimental version 3 of gcc.
80 *
81 * Revision 1.9 2002/05/31 16:13:36 j_novak
82 * better inversion for eos_bifluid
83 *
84 * Revision 1.8 2002/05/10 09:26:52 j_novak
85 * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
86 *
87 * Revision 1.7 2002/05/02 15:16:22 j_novak
88 * Added functions for more general bi-fluid EOS
89 *
90 * Revision 1.6 2002/04/05 09:09:36 j_novak
91 * The inversion of the EOS for 2-fluids polytrope has been modified.
92 * Some errors in the determination of the surface were corrected.
93 *
94 * Revision 1.5 2002/01/16 15:03:28 j_novak
95 * *** empty log message ***
96 *
97 * Revision 1.4 2002/01/11 14:09:34 j_novak
98 * Added newtonian version for 2-fluid stars
99 *
100 * Revision 1.3 2001/12/04 21:27:53 e_gourgoulhon
101 *
102 * All writing/reading to a binary file are now performed according to
103 * the big endian convention, whatever the system is big endian or
104 * small endian, thanks to the functions fwrite_be and fread_be
105 *
106 * Revision 1.2 2001/11/29 15:05:26 j_novak
107 * The entrainment term in 2-fluid eos is modified
108 *
109 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
110 * LORENE
111 *
112 * Revision 1.4 2001/08/31 15:48:50 novak
113 * The flag tronc has been added to nbar_ent_p
114 *
115 * Revision 1.3 2001/08/27 09:52:21 novak
116 * New version of formulas
117 *
118 * Revision 1.2 2001/06/22 15:36:46 novak
119 * Modification de trans2Eos
120 *
121 * Revision 1.1 2001/06/21 15:24:46 novak
122 * Initial revision
123 *
124 *
125 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_poly.C,v 1.22 2014/10/13 08:52:52 j_novak Exp $
126 *
127 */
128
129// Headers C
130#include <cstdlib>
131#include <cmath>
132
133// Headers Lorene
134#include "eos_bifluid.h"
135#include "param.h"
136#include "eos.h"
137#include "cmp.h"
138#include "utilitaires.h"
139
140namespace Lorene {
141
142//****************************************************************************
143// local prototypes
144double puis(double, double) ;
145double enthal1(const double x, const Param& parent) ;
146double enthal23(const double x, const Param& parent) ;
147double enthal(const double x, const Param& parent) ;
148//****************************************************************************
149
150 //--------------//
151 // Constructors //
152 //--------------//
153
154// Standard constructor with gam1 = gam2 = 2,
155// gam3 = gam4 = gam5 = gam6 = 1, m_1 = 1 and m_2 =1
156// -------------------------------------------------
157Eos_bf_poly::Eos_bf_poly(double kappa1, double kappa2, double kappa3, double bet) :
158 Eos_bifluid("bi-fluid polytropic EOS", 1, 1),
159 gam1(2), gam2(2),gam3(1),gam4(1),gam5(1),
160 gam6(1),kap1(kappa1), kap2(kappa2), kap3(kappa3),beta(bet),
161 relax(0.5), precis(1.e-9), ecart(1.e-8)
162{
163
165 set_auxiliary() ;
166
167}
168
169// Standard constructor with everything specified
170// -----------------------------------------------
171Eos_bf_poly::Eos_bf_poly(double gamma1, double gamma2, double gamma3,
172 double gamma4, double gamma5, double gamma6,
173 double kappa1, double kappa2, double kappa3,
174 double bet, double mass1, double mass2,
175 double rel, double prec, double ec) :
176 Eos_bifluid("bi-fluid polytropic EOS", mass1, mass2),
177 gam1(gamma1),gam2(gamma2),gam3(gamma3),gam4(gamma4),gam5(gamma5),
178 gam6(gamma6),kap1(kappa1),kap2(kappa2),kap3(kappa3),beta(bet),
179 relax(rel), precis(prec), ecart(ec)
180{
181
183 set_auxiliary() ;
184
185}
186
187// Copy constructor
188// ----------------
190 Eos_bifluid(eosi),
191 gam1(eosi.gam1), gam2(eosi.gam2), gam3(eosi.gam3),
192 gam4(eosi.gam4), gam5(eosi.gam5), gam6(eosi.gam6),
193 kap1(eosi.kap1), kap2(eosi.kap2), kap3(eosi.kap3),
194 beta(eosi.beta),
195 typeos(eosi.typeos), relax(eosi.relax), precis(eosi.precis),
196 ecart(eosi.ecart) {
197
198 set_auxiliary() ;
199
200}
201
202
203// Constructor from binary file
204// ----------------------------
206 Eos_bifluid(fich) {
207
208 fread_be(&gam1, sizeof(double), 1, fich) ;
209 fread_be(&gam2, sizeof(double), 1, fich) ;
210 fread_be(&gam3, sizeof(double), 1, fich) ;
211 fread_be(&gam4, sizeof(double), 1, fich) ;
212 fread_be(&gam5, sizeof(double), 1, fich) ;
213 fread_be(&gam6, sizeof(double), 1, fich) ;
214 fread_be(&kap1, sizeof(double), 1, fich) ;
215 fread_be(&kap2, sizeof(double), 1, fich) ;
216 fread_be(&kap3, sizeof(double), 1, fich) ;
217 fread_be(&beta, sizeof(double), 1, fich) ;
218 fread_be(&relax, sizeof(double), 1, fich) ;
219 fread_be(&precis, sizeof(double), 1, fich) ;
220 fread_be(&ecart, sizeof(double), 1, fich) ;
221
223 set_auxiliary() ;
224
225
226}
227
228
229// Constructor from a formatted file
230// ---------------------------------
231Eos_bf_poly::Eos_bf_poly(const char *fname ) :
232 Eos_bifluid(fname), relax(0.5), precis(1.e-8), ecart(1.e-7)
233{
234 int res = 0;
235
236 res += read_variable (fname, const_cast<char*>("gamma1"), gam1);
237 res += read_variable (fname, const_cast<char*>("gamma2"), gam2);
238 res += read_variable (fname, const_cast<char*>("gamma3"), gam3);
239 res += read_variable (fname, const_cast<char*>("gamma4"), gam4);
240 res += read_variable (fname, const_cast<char*>("gamma5"), gam5);
241 res += read_variable (fname, const_cast<char*>("gamma6"), gam6);
242 res += read_variable (fname, const_cast<char*>("kappa1"), kap1);
243 res += read_variable (fname, const_cast<char*>("kappa2"), kap2);
244 res += read_variable (fname, const_cast<char*>("kappa3"), kap3);
245 res += read_variable (fname, const_cast<char*>("beta"), beta);
246
247
248 if (res != 0)
249 {
250 cerr << "ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
251 exit (-1);
252 }
253
255
256 if (get_typeos() == 4)
257 {
258 res += read_variable (fname, const_cast<char*>("relax"), relax);
259 res += read_variable (fname, const_cast<char*>("precis"), precis);
260 res += read_variable (fname, const_cast<char*>("ecart"), ecart);
261 }
262 else if (get_typeos() == 0) // analytic EOS: check if we want slow-rot-style inversion
263 {
264 bool slowrot = false;
265 read_variable (fname, const_cast<char*>("slow_rot_style"), slowrot); // dont require this variable!
266 if (slowrot)
267 typeos = 5; // type=5 is reserved for (type0 + slow-rot-style)
268 }
269
270
271 if (res != 0)
272 {
273 cerr << "ERROR: could not read all required eos-paramters for Eos_bf_poly()" << endl;
274 exit (-1);
275 }
276
277
278 set_auxiliary() ;
279
280}
281 //--------------//
282 // Destructor //
283 //--------------//
284
286
287 //Does nothing ...
288
289}
290 //--------------//
291 // Assignment //
292 //--------------//
293
295
296 // Assignment of quantities common to all the derived classes of Eos_bifluid
298
299 gam1 = eosi.gam1 ;
300 gam2 = eosi.gam2 ;
301 gam3 = eosi.gam3 ;
302 kap1 = eosi.kap1 ;
303 kap2 = eosi.kap2 ;
304 kap3 = eosi.kap3 ;
305 beta = eosi.beta ;
306 typeos = eosi.typeos ;
307 relax = eosi.relax ;
308 precis = eosi.precis ;
309 ecart = eosi.ecart ;
310
311 set_auxiliary() ;
312
313}
314
315
316 //-----------------------//
317 // Miscellaneous //
318 //-----------------------//
319
321
322 gam1m1 = gam1 - double(1) ;
323 gam2m1 = gam2 - double(1) ;
324 gam34m1 = gam3 + gam4 - double(1) ;
325 gam56m1 = gam5 + gam6 - double(1) ;
326
327 if (fabs(kap3*kap3-kap2*kap1) < 1.e-15) {
328 cout << "WARNING!: Eos_bf_poly: the parameters are degenerate!" << endl ;
329 abort() ;
330 }
331
332}
333
335
336 if ((gam1 == double(2)) && (gam2 == double(2)) && (gam3 == double(1))
337 && (gam4 == double(1)) && (gam5 == double(1))
338 && (gam6 == double(0))) {
339 typeos = -1 ;
340 cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
341 cout << "WARNING!! The entrainment factor does not depend on" <<endl ;
342 cout << " density 2! This may be incorrect and should only be used"<<endl ;
343 cout << " for testing purposes..." << endl ;
344 cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
345 }
346 else if ((gam1 == double(2)) && (gam2 == double(2)) && (gam3 == double(1))
347 && (gam4 == double(1)) && (gam5 == double(0))
348 && (gam6 == double(1))) {
349 typeos = -2 ;
350 cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
351 cout << "WARNING!! The entrainment factor does not depend on" << endl ;
352 cout << " density 1! This may be incorrect and should only be used"<<endl ;
353 cout << " for testing purposes..." << endl ;
354 cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
355 }
356 else if ((gam1 == double(2)) && (gam2 == double(2)) && (gam3 == double(1))
357 && (gam4 == double(1)) && (gam5 == double(1))
358 && (gam6 == double(1))) {
359 typeos = 0 ;
360 }
361 else if ((gam3 == double(1)) && (gam4 == double(1)) && (gam5 == double(1))
362 && (gam6 == double(1))) {
363 typeos = 1 ;
364 }
365 else if ((gam3 == double(1)) && (gam5 == double(1))) {
366 typeos = 2 ;
367 }
368 else if ((gam4 == double(1)) && (gam6 == double(1))) {
369 typeos = 3 ;
370 return ;
371 }
372 else {
373 typeos = 4 ;
374 }
375
376 cout << "DEBUG: EOS-type was determined as typeos = " << typeos << endl;
377
378 return ;
379}
380
381 //------------------------//
382 // Comparison operators //
383 //------------------------//
384
385
386bool Eos_bf_poly::operator==(const Eos_bifluid& eos_i) const {
387
388 bool resu = true ;
389
390 if ( eos_i.identify() != identify() ) {
391 cout << "The second EOS is not of type Eos_bf_poly !" << endl ;
392 resu = false ;
393 }
394 else{
395
396 const Eos_bf_poly& eos = dynamic_cast<const Eos_bf_poly&>( eos_i ) ;
397
398 if ((eos.gam1 != gam1)||(eos.gam2 != gam2)||(eos.gam3 != gam3)
399 ||(eos.gam4 != gam4)||(eos.gam5 != gam5)||(eos.gam6 != gam6)) {
400 cout
401 << "The two Eos_bf_poly have different gammas : " << gam1 << " <-> "
402 << eos.gam1 << ", " << gam2 << " <-> "
403 << eos.gam2 << ", " << gam3 << " <-> "
404 << eos.gam3 << ", " << gam4 << " <-> "
405 << eos.gam4 << ", " << gam5 << " <-> "
406 << eos.gam5 << ", " << gam6 << " <-> "
407 << eos.gam6 << endl ;
408 resu = false ;
409 }
410
411 if ((eos.kap1 != kap1)||(eos.kap2 != kap2)|| (eos.kap3 != kap3)){
412 cout
413 << "The two Eos_bf_poly have different kappas : " << kap1 << " <-> "
414 << eos.kap1 << ", " << kap2 << " <-> "
415 << eos.kap2 << ", " << kap3 << " <-> "
416 << eos.kap3 << endl ;
417 resu = false ;
418 }
419
420 if (eos.beta != beta) {
421 cout
422 << "The two Eos_bf_poly have different betas : " << beta << " <-> "
423 << eos.beta << endl ;
424 resu = false ;
425 }
426
427 if ((eos.m_1 != m_1)||(eos.m_2 != m_2)) {
428 cout
429 << "The two Eos_bf_poly have different masses : " << m_1 << " <-> "
430 << eos.m_1 << ", " << m_2 << " <-> "
431 << eos.m_2 << endl ;
432 resu = false ;
433 }
434
435 }
436
437 return resu ;
438
439}
440
441bool Eos_bf_poly::operator!=(const Eos_bifluid& eos_i) const {
442
443 return !(operator==(eos_i)) ;
444
445}
446
447
448 //------------//
449 // Outputs //
450 //------------//
451
452void Eos_bf_poly::sauve(FILE* fich) const {
453
454 Eos_bifluid::sauve(fich) ;
455
456 fwrite_be(&gam1, sizeof(double), 1, fich) ;
457 fwrite_be(&gam2, sizeof(double), 1, fich) ;
458 fwrite_be(&gam3, sizeof(double), 1, fich) ;
459 fwrite_be(&gam4, sizeof(double), 1, fich) ;
460 fwrite_be(&gam5, sizeof(double), 1, fich) ;
461 fwrite_be(&gam6, sizeof(double), 1, fich) ;
462 fwrite_be(&kap1, sizeof(double), 1, fich) ;
463 fwrite_be(&kap2, sizeof(double), 1, fich) ;
464 fwrite_be(&kap3, sizeof(double), 1, fich) ;
465 fwrite_be(&beta, sizeof(double), 1, fich) ;
466 fwrite_be(&relax, sizeof(double), 1, fich) ;
467 fwrite_be(&precis, sizeof(double), 1, fich) ;
468 fwrite_be(&ecart, sizeof(double), 1, fich) ;
469
470}
471
472ostream& Eos_bf_poly::operator>>(ostream & ost) const {
473
474 ost << "EOS of class Eos_bf_poly (relativistic polytrope) : " << endl ;
475 ost << " Adiabatic index gamma1 : " << gam1 << endl ;
476 ost << " Adiabatic index gamma2 : " << gam2 << endl ;
477 ost << " Adiabatic index gamma3 : " << gam3 << endl ;
478 ost << " Adiabatic index gamma4 : " << gam4 << endl ;
479 ost << " Adiabatic index gamma5 : " << gam5 << endl ;
480 ost << " Adiabatic index gamma6 : " << gam6 << endl ;
481 ost << " Pressure coefficient kappa1 : " << kap1 <<
482 " rho_nuc c^2 / n_nuc^gamma" << endl ;
483 ost << " Pressure coefficient kappa2 : " << kap2 <<
484 " rho_nuc c^2 / n_nuc^gamma" << endl ;
485 ost << " Pressure coefficient kappa3 : " << kap3 <<
486 " rho_nuc c^2 / n_nuc^gamma" << endl ;
487 ost << " Coefficient beta : " << beta <<
488 " rho_nuc c^2 / n_nuc^gamma" << endl ;
489 ost << " Type for inversion : " << typeos << endl ;
490 ost << " Parameters for inversion (used if typeos = 4) : " << endl ;
491 ost << " relaxation : " << relax << endl ;
492 ost << " precision for zerosec_b : " << precis << endl ;
493 ost << " final discrepancy in the result : " << ecart << endl ;
494
495 return ost ;
496
497}
498
499
500 //------------------------------//
501 // Computational routines //
502 //------------------------------//
503
504// Baryon densities from enthalpies
505//---------------------------------
506
507bool Eos_bf_poly::nbar_ent_p(const double ent1, const double ent2,
508 const double delta2, double& nbar1,
509 double& nbar2) const {
510
511 // RP: follow Newtonian case in this: the following is wrong, I think
512 // bool one_fluid = ((ent1<=0.)||(ent2<=0.)) ;
513
514 bool one_fluid = false;
515
516 if (!one_fluid) {
517 switch (typeos) {
518 case 5: // same as typeos=0 but with slow-rot-style inversion
519 case 0: {
520 double kpd = kap3+beta*delta2 ;
521 double determ = kap1*kap2 - kpd*kpd ;
522
523 nbar1 = (kap2*(exp(ent1) - m_1) - kpd*(exp(ent2) - m_2)) / determ ;
524 nbar2 = (kap1*(exp(ent2) - m_2) - kpd*(exp(ent1) - m_1)) / determ ;
525 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
526 break ;
527 }
528 case 1: {
529 double b1 = exp(ent1) - m_1 ;
530 double b2 = exp(ent2) - m_2 ;
531 double denom = kap3 + beta*delta2 ;
532 if (fabs(denom) < 1.e-15) {
533 nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
534 nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
535 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
536 }
537 else {
538 double coef0 = 0.5*gam2*kap2/pow(denom, gam2m1) ;
539 double coef1 = 0.5*kap1*gam1 ;
540 Param parent ;
541 parent.add_double(coef0, 0) ;
542 parent.add_double(b1, 1) ;
543 parent.add_double(coef1, 2) ;
544 parent.add_double(gam1m1,3) ;
545 parent.add_double(gam2m1,4) ;
546 parent.add_double(denom, 5) ;
547 parent.add_double(b2, 6) ;
548
549 double xmin, xmax ;
550 double f0 = enthal1(0.,parent) ;
551 if (fabs(f0)<1.e-15) {
552 nbar1 = 0. ;}
553 else {
554 double cas = (gam1m1*gam2m1 - 1.)*f0;
555 assert (fabs(cas) > 1.e-15) ;
556 xmin = 0. ;
557 xmax = cas/fabs(cas) ;
558 do {
559 xmax *= 1.1 ;
560 if (fabs(xmax) > 1.e10) {
561 cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
562 cout << f0 << ", " << cas << endl ; // to be removed!!
563 cout << "typeos = 1" << endl ;
564 abort() ;
565 }
566 } while (enthal1(xmax,parent)*f0 > 0.) ;
567 double l_precis = 1.e-12 ;
568 int nitermax = 400 ;
569 int niter = 0 ;
570 nbar1 = zerosec_b(enthal1, parent, xmin, xmax, l_precis,
571 nitermax, niter) ;
572 }
573 nbar2 = (b1 - coef1*puis(nbar1, gam1m1))/denom ;
574 double resu1 = coef1*puis(nbar1,gam1m1) + denom*nbar2 ;
575 double resu2 = 0.5*gam2*kap2*puis(nbar2,gam2m1) + denom*nbar1 ;
576 double erreur = fabs((log(fabs(1+resu1))-ent1)/ent1) +
577 fabs((log(fabs(1+resu2))-ent2)/ent2) ;
578 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
579 }
580 break ;
581 }
582 case 2: {
583 double b1 = exp(ent1) - m_1 ;
584 double b2 = exp(ent2) - m_2 ;
585 double denom = kap3 + beta*delta2 ;
586 if (fabs(denom) < 1.e-15) {
587 nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
588 nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
589 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
590 }
591 else {
592 double coef0 = beta*delta2 ;
593 double coef1 = 0.5*kap1*gam1 ;
594 double coef2 = 0.5*kap2*gam2 ;
595 double coef3 = 1./gam1m1 ;
596 Param parent ;
597 parent.add_double(b1, 0) ;
598 parent.add_double(kap3, 1) ;
599 parent.add_double(gam4, 2) ;
600 parent.add_double(coef0, 3) ;
601 parent.add_double(gam6,4) ;
602 parent.add_double(coef1, 5) ;
603 parent.add_double(coef3, 6) ;
604 parent.add_double(coef2, 7) ;
605 parent.add_double(gam2m1, 8) ;
606 parent.add_double(b2, 9) ;
607
608 double xmin, xmax ;
609 double f0 = enthal23(0.,parent) ;
610 if (fabs(f0)<1.e-15) {
611 nbar2 = 0. ;}
612 else {
613 double pmax = (fabs(kap3) < 1.e-15 ? 0. : gam4*(gam4-1) ) ;
614 double ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0.
615 : gam6*(gam4-1) ) ;
616 pmax = (pmax>ptemp ? pmax : ptemp) ;
617 ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0. : gam4*(gam6-1) ) ;
618 pmax = (pmax>ptemp ? pmax : ptemp) ;
619 ptemp = (fabs(coef0) < 1.e-15 ? 0. : gam6*(gam6-1) ) ;
620 pmax = (pmax>ptemp ? pmax : ptemp) ;
621 double cas = (pmax - gam1m1*gam2m1)*f0;
622 // cout << "Pmax, cas: " << pmax << ", " << cas << endl ;
623 assert (fabs(cas) > 1.e-15) ;
624 xmin = 0. ;
625 xmax = cas/fabs(cas) ;
626 do {
627 xmax *= 1.1 ;
628 if (fabs(xmax) > 1.e10) {
629 cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
630 cout << "typeos = 2" << endl ;
631 abort() ;
632 }
633 } while (enthal23(xmax,parent)*f0 > 0.) ;
634 double l_precis = 1.e-12 ;
635 int nitermax = 400 ;
636 int niter = 0 ;
637 nbar2 = zerosec_b(enthal23, parent, xmin, xmax, l_precis,
638 nitermax, niter) ;
639 }
640 nbar1 = (b1 - kap3*puis(nbar2,gam4) - coef0*puis(nbar2,gam6))
641 / coef1 ;
642 nbar1 = puis(nbar1,coef3) ;
643 double resu1 = coef1*puis(nbar1,gam1m1) + kap3*puis(nbar2,gam4)
644 + coef0*puis(nbar2, gam6) ;
645 double resu2 = coef2*puis(nbar2,gam2m1)
646 + gam4*kap3*puis(nbar2, gam4-1)*nbar1
647 + gam6*coef0*puis(nbar2, gam6-1)*nbar1 ;
648 double erreur = fabs((log(fabs(1+resu1))-ent1)/ent1) +
649 fabs((log(fabs(1+resu2))-ent2)/ent2) ;
650 //cout << "Erreur d'inversion: " << erreur << endl ;
651 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
652 }
653 break ;
654 }
655 case 3: {
656 double b1 = exp(ent1) - m_1 ;
657 double b2 = exp(ent2) - m_2 ;
658 double denom = kap3 + beta*delta2 ;
659 if (fabs(denom) < 1.e-15) {
660 nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
661 nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
662 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
663 }
664 else {
665 double coef0 = beta*delta2 ;
666 double coef1 = 0.5*kap1*gam1 ;
667 double coef2 = 0.5*kap2*gam2 ;
668 double coef3 = 1./gam2m1 ;
669 Param parent ;
670 parent.add_double(b2, 0) ;
671 parent.add_double(kap3, 1) ;
672 parent.add_double(gam3, 2) ;
673 parent.add_double(coef0, 3) ;
674 parent.add_double(gam5, 4) ;
675 parent.add_double(coef2, 5) ;
676 parent.add_double(coef3, 6) ;
677 parent.add_double(coef1, 7) ;
678 parent.add_double(gam1m1, 8) ;
679 parent.add_double(b1, 9) ;
680
681 double xmin, xmax ;
682 double f0 = enthal23(0.,parent) ;
683 if (fabs(f0)<1.e-15) {
684 nbar1 = 0. ;}
685 else {
686 double pmax = (fabs(kap3) < 1.e-15 ? 0. : gam3*(gam3-1) ) ;
687 double ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0.
688 : gam5*(gam3-1) ) ;
689 pmax = (pmax>ptemp ? pmax : ptemp) ;
690 ptemp = (fabs(kap3*coef0) < 1.e-15 ? 0. : gam3*(gam5-1) ) ;
691 pmax = (pmax>ptemp ? pmax : ptemp) ;
692 ptemp = (fabs(coef0) < 1.e-15 ? 0. : gam5*(gam5-1) ) ;
693 pmax = (pmax>ptemp ? pmax : ptemp) ;
694 double cas = (pmax - gam1m1*gam2m1)*f0;
695 // cout << "Pmax, cas: " << pmax << ", " << cas << endl ;
696 assert (fabs(cas) > 1.e-15) ;
697 xmin = 0. ;
698 xmax = cas/fabs(cas) ;
699 do {
700 xmax *= 1.1 ;
701 if (fabs(xmax) > 1.e10) {
702 cout << "Problem in Eos_bf_poly::nbar_ent_p!" << endl ;
703 cout << "typeos = 3" << endl ;
704 abort() ;
705 }
706 } while (enthal23(xmax,parent)*f0 > 0.) ;
707 double l_precis = 1.e-12 ;
708 int nitermax = 400 ;
709 int niter = 0 ;
710 nbar1 = zerosec_b(enthal23, parent, xmin, xmax, l_precis,
711 nitermax, niter) ;
712 }
713 nbar2 = (b2 - kap3*puis(nbar1,gam3) - coef0*puis(nbar1,gam5))
714 / coef2 ;
715 nbar2 = puis(nbar2,coef3) ;
716 double resu1 = coef1*puis(nbar1,gam1m1)
717 + gam3*kap3*puis(nbar1,gam3-1)*nbar2
718 + coef0*gam5*puis(nbar1, gam5-1)*nbar2 ;
719 double resu2 = coef2*puis(nbar2,gam2m1)
720 + kap3*puis(nbar1, gam3) + coef0*puis(nbar1, gam5);
721 double erreur = fabs((log(fabs(1+resu1))-ent1)/ent1) +
722 fabs((log(fabs(1+resu2))-ent2)/ent2) ;
723 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)||(erreur > 1.e-4)) ;
724 }
725 break ;
726 }
727 case 4:{
728 double b1 = exp(ent1) - m_1 ;
729 double b2 = exp(ent2) - m_2 ;
730 double denom = kap3 + beta*delta2 ;
731 if (fabs(denom) < 1.e-15) {
732 nbar1 = puis(2*b1/(gam1*kap1), 1./gam1m1) ;
733 nbar2 = puis(2*b2/(gam2*kap2), 1./gam2m1) ;
734 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
735 }
736 else {
737 int nitermax = 200 ;
738 int niter = 0 ;
739 int nmermax = 800 ;
740
741 double a11 = 0.5*gam1*kap1 ;
742 double a13 = gam3*kap3 ;
743 double a14 = beta*gam5*delta2 ;
744
745 double a22 = 0.5*gam2*kap2 ;
746 double a23 = gam4*kap3 ;
747 double a24 = beta*gam6*delta2 ;
748
749 double n1l, n2l, n1s, n2s ;
750
751 double delta = a11*a22 - (a13+a14)*(a23+a24) ;
752 n1l = (a22*b1 - (a13+a14)*b2)/delta ;
753 n2l = (a11*b2 - (a23+a24)*b1)/delta ;
754 n1s = puis(b1/a11, 1./(gam1m1)) ;
755 n2s = puis(b2/a22, 1./(gam2m1)) ;
756
757 double n1m = (n1l<0. ? n1s : sqrt(n1l*n1s)) ;
758 double n2m = (n2l<0. ? n2s : sqrt(n2l*n2s)) ;
759
760 Param parf1 ;
761 Param parf2 ;
762 double c1, c2, c3, c4, c5, c6, c7 ;
763 c1 = gam1m1 ;
764 c2 = gam3 - 1. ;
765 c3 = gam5 - 1. ;
766 c4 = a11 ;
767 c5 = a13*puis(n2m,gam4) ;
768 c6 = a14*puis(n2m,gam6) ;
769 c7 = b1 ;
770 parf1.add_double_mod(c1,0) ;
771 parf1.add_double_mod(c2,1) ;
772 parf1.add_double_mod(c3,2) ;
773 parf1.add_double_mod(c4,3) ;
774 parf1.add_double_mod(c5,4) ;
775 parf1.add_double_mod(c6,5) ;
776 parf1.add_double_mod(c7,6) ;
777
778 double d1, d2, d3, d4, d5, d6, d7 ;
779 d1 = gam2m1 ;
780 d2 = gam4 - 1. ;
781 d3 = gam6 - 1. ;
782 d4 = a22 ;
783 d5 = a23*puis(n1m,gam3) ;
784 d6 = a24*puis(n1m,gam5) ;
785 d7 = b2 ;
786 parf2.add_double_mod(d1,0) ;
787 parf2.add_double_mod(d2,1) ;
788 parf2.add_double_mod(d3,2) ;
789 parf2.add_double_mod(d4,3) ;
790 parf2.add_double_mod(d5,4) ;
791 parf2.add_double_mod(d6,5) ;
792 parf2.add_double_mod(d7,6) ;
793
794 double xmin = -3*(n1s>n2s ? n1s : n2s) ;
795 double xmax = 3*(n1s>n2s ? n1s : n2s) ;
796
797 double n1 = n1m ;
798 double n2 = n2m ;
799 bool sortie = true ;
800 int mer = 0 ;
801
802 //cout << "Initial guess: " << n1 << ", " << n2 << endl ;
803 n1s *= 0.1 ;
804 n2s *= 0.1 ;
805 do {
806 //cout << "n1, n2: " << n1 << ", " << n2 << endl ;
807 n1 = zerosec_b(&enthal, parf1, xmin, xmax, precis, nitermax, niter) ;
808 n2 = zerosec_b(&enthal, parf2, xmin, xmax, precis, nitermax, niter) ;
809
810 sortie = (fabs((n1m-n1)/(n1m+n1)) + fabs((n2m-n2)/(n2m+n2)) > ecart) ;
811 n1m = relax*n1 + (1.-relax)*n1m ;
812 n2m = relax*n2 + (1.-relax)*n2m ;
813 if (n2m>-n2s) {
814 parf1.get_double_mod(4) = a13*puis(n2m,gam4) ;
815 parf1.get_double_mod(5) = a14*puis(n2m,gam6) ;
816 }
817 else {
818 parf1.get_double_mod(4) = a13*puis(-n2s,gam4) ;
819 parf1.get_double_mod(5) = a14*puis(-n2s,gam6) ;
820 }
821
822 if (n1m>-n1s) {
823 parf2.get_double_mod(4) = a23*puis(n1m,gam3) ;
824 parf2.get_double_mod(5) = a24*puis(n1m,gam5) ;
825 }
826 else {
827 parf2.get_double_mod(4) = a23*puis(-n1s,gam3) ;
828 parf2.get_double_mod(5) = a24*puis(-n1s,gam5) ;
829 }
830
831 mer++ ;
832 } while (sortie&&(mer<nmermax)) ;
833 nbar1 = n1m ;
834 nbar2 = n2m ;
835
836// double resu1 = a11*puis(n1,gam1m1) + a13*puis(n1,gam3-1.)*puis(n2,gam4)
837// +a14*puis(n1,gam5-1.)*puis(n2,gam6) ;
838// double resu2 = a22*puis(n2,gam2m1) + a23*puis(n1,gam3)*puis(n2,gam4-1.)
839// +a24*puis(n1,gam5)*puis(n2,gam6-1.) ;
840 //cout << "Nbre d'iterations: " << mer << endl ;
841 //cout << "Resus: " << n1m << ", " << n2m << endl ;
842 //cout << "Verification: " << log(fabs(1+resu1)) << ", "
843 // << log(fabs(1+resu2)) << endl ;
844 // cout << "Erreur: " << fabs(enthal(n1, parf1)/b1) +
845 // fabs(enthal(n2, parf2)/b2) << endl ;
846 //cout << "Erreur: " << fabs((log(fabs(1+resu1))-ent1)/ent1) +
847 //fabs((log(fabs(1+resu2))-ent2)/ent2) << endl ;
848 }
849 break ;
850 }
851 case -1: {
852 double determ = kap1*kap2 - kap3*kap3 ;
853
854 nbar1 = (kap2*(exp(ent1) - m_1 - beta*delta2)
855 - kap3*(exp(ent2) - m_2)) / determ ;
856 nbar2 = (kap1*(exp(ent2) - m_2) - kap3*(exp(ent1) - m_1 - beta*delta2))
857 / determ ;
858 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
859 break ;
860 }
861 case -2: {
862 double determ = kap1*kap2 - kap3*kap3 ;
863
864 nbar1 = (kap2*(exp(ent1) - m_1 )
865 - kap3*(exp(ent2) - m_2 - beta*delta2)) / determ ;
866 nbar2 = (kap1*(exp(ent2) - m_2 - beta*delta2)
867 - kap3*(exp(ent1) - m_1)) / determ ;
868 one_fluid = ((nbar1 < 0.)||(nbar2 < 0.)) ;
869 break ;
870 }
871 }
872 }
873 return one_fluid ;
874}
875// One fluid sub-EOSs
876//-------------------
877
878double Eos_bf_poly::nbar_ent_p1(const double ent1) const {
879 return puis(2*(exp(ent1) - m_1)/(gam1*kap1), 1./gam1m1) ;
880}
881
882double Eos_bf_poly::nbar_ent_p2(const double ent2) const {
883 return puis(2*(exp(ent2) - m_2)/(gam2*kap2), 1./gam2m1) ;
884}
885
886// Energy density from baryonic densities
887//---------------------------------------
888
889double Eos_bf_poly::ener_nbar_p(const double nbar1, const double nbar2,
890 const double delta2) const {
891
892 if (( nbar1 > double(0) ) || ( nbar2 > double(0))) {
893
894 double n1 = (nbar1>double(0) ? nbar1 : double(0)) ;
895 double n2 = (nbar2>double(0) ? nbar2 : double(0)) ;
896 double x2 = ((nbar1>double(0))&&(nbar2>double(0))) ? delta2 : 0 ;
897
898 double resu = 0.5*kap1*pow(n1, gam1) + 0.5*kap2*pow(n2,gam2)
899 + kap3*pow(n1,gam3)*pow(n2,gam4) + m_1*n1 + m_2*n2
900 + x2*beta*pow(n1,gam5)*pow(n2,gam6) ;
901 return resu ;
902 }
903 else return 0 ;
904}
905
906// Pressure from baryonic densities
907//---------------------------------
908
909double Eos_bf_poly::press_nbar_p(const double nbar1, const double nbar2,
910 const double delta2) const {
911
912 if (( nbar1 > double(0) ) || ( nbar2 > double(0))) {
913
914 double n1 = (nbar1>double(0) ? nbar1 : double(0)) ;
915 double n2 = (nbar2>double(0) ? nbar2 : double(0)) ;
916 double x2 = ((nbar1>double(0))&&(nbar2>double(0))) ? delta2 : 0 ;
917
918 double resu = 0.5*gam1m1*kap1*pow(n1,gam1) + 0.5*gam2m1*kap2*pow(n2,gam2)
919 + gam34m1*kap3*pow(n1,gam3)*pow(n2,gam4) +
920 x2*gam56m1*beta*pow(n1,gam5)*pow(n2,gam6) ;
921
922 return resu ;
923 }
924 else return 0 ;
925}
926
927// Derivatives of energy
928//----------------------
929
930double Eos_bf_poly::get_K11(const double n1, const double n2, const
931 double delta2) const
932{
933 double xx ;
934 if (n1 <= 0.) {
935 xx = 0. ;
936 }
937 else {
938 xx = 0.5*gam1*kap1 * pow(n1,gam1 - 2) + m_1/n1 +
939 gam3*kap3 * pow(n1,gam3 - 2) * pow(n2,gam4) +
940 (delta2*(gam5 + 2) - 2)*beta * pow(n1,gam5 - 2)*pow(n2, gam6) ;
941 }
942 return xx ;
943}
944
945double Eos_bf_poly::get_K22(const double n1, const double n2, const
946 double delta2) const
947{
948 double xx ;
949 if (n2 <= 0.) {
950 xx = 0. ;
951 }
952 else {
953 xx = 0.5*gam2*kap2 * pow(n2,gam2 - 2) + m_2/n2 +
954 gam4*kap3 * pow(n2,gam4 - 2) * pow(n1,gam3) +
955 (delta2*(gam6 + 2) - 2)*beta * pow(n1, gam5) * pow(n2,gam6 - 2) ;
956 }
957 return xx ;
958}
959
960double Eos_bf_poly::get_K12(const double n1, const double n2, const
961 double delta2) const
962{
963 double xx ;
964 if ((n1 <= 0.) || (n2 <= 0.)) { xx = 0.; }
965 else {
966 double gamma_delta3 = pow(1-delta2,-1.5) ;
967 xx = 2*beta*pow(n1,gam5-1)*pow(n2,gam6-1) / gamma_delta3 ;
968 }
969 return xx ;
970}
971
972// Conversion functions
973// ---------------------
974
976
977 Eos_poly* eos_simple = new Eos_poly(gam1, kap1, m_1) ;
978 return eos_simple ;
979}
980/***************************************************************************
981 *
982 * Non-class functions
983 *
984 ***************************************************************************/
985
986// New "pow"
987//-----------
988
989 double puis(double x, double p) {
990 assert(p>=0.) ;
991 if (p==0.) return (x>=0 ? 1 : -1) ;
992 //if (p==0.) return 1 ;
993 if (x<0.) return (-pow(-x,p)) ;
994 else return pow(x,p) ;
995}
996
997// Auxilliary functions for nbar_ent_p
998//------------------------------------
999
1000double enthal1(const double x, const Param& parent) {
1001 assert(parent.get_n_double() == 7) ;
1002
1003 return parent.get_double(0)*puis(parent.get_double(1) - parent.get_double(2)
1004 *puis(x,parent.get_double(3)), parent.get_double(4))
1005 + parent.get_double(5)*x - parent.get_double(6) ;
1006
1007}
1008
1009double enthal23(const double x, const Param& parent) {
1010 assert(parent.get_n_double() == 10) ;
1011
1012 double nx = (parent.get_double(0) - parent.get_double(1)*
1013 puis(x,parent.get_double(2)) - parent.get_double(3)*
1014 puis(x,parent.get_double(4)) )/parent.get_double(5) ;
1015 nx = puis(nx,parent.get_double(6)) ;
1016 return parent.get_double(7)*puis(x,parent.get_double(8))
1017 + parent.get_double(1)*parent.get_double(2)*nx*
1018 puis(x,parent.get_double(2) - 1)
1019 + parent.get_double(3)*parent.get_double(4)*nx*
1020 puis(x,parent.get_double(4) - 1)
1021 - parent.get_double(9) ;
1022
1023}
1024
1025double enthal(const double x, const Param& parent) {
1026 assert(parent.get_n_double_mod() == 7) ;
1027
1028 double alp1 = parent.get_double_mod(0) ;
1029 double alp2 = parent.get_double_mod(1) ;
1030 double alp3 = parent.get_double_mod(2) ;
1031 double cc1 = parent.get_double_mod(3) ;
1032 double cc2 = parent.get_double_mod(4) ;
1033 double cc3 = parent.get_double_mod(5) ;
1034 double cc4 = parent.get_double_mod(6) ;
1035
1036 return (cc1*puis(x,alp1) + cc2*puis(x,alp2) + cc3*puis(x,alp3) - cc4) ;
1037
1038}
1039
1040}
1041
Analytic equation of state for two fluids (relativistic case).
double gam5
Adiabatic indexes , see Eq.~eeosbfpolye}.
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
double kap1
Pressure coefficient , see Eq.
void determine_type()
Determines the type of the analytical EOS (see typeos )
double gam1
Adiabatic indexes , see Eq.~eeosbfpolye}.
double kap2
Pressure coefficient , see Eq.
void operator=(const Eos_bf_poly &)
Assignment to another Eos_bf_poly.
double gam4
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam3
Adiabatic indexes , see Eq.~eeosbfpolye}.
virtual double get_K11(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to (baryonic density 1) .
virtual double get_K12(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy with respect to .
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
virtual ostream & operator>>(ostream &) const
Operator >>
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
virtual double get_K22(const double n1, const double n2, const double delta2) const
Computes the derivative of the energy/(baryonic density 2) .
double gam6
Adiabatic indexes , see Eq.~eeosbfpolye}.
double gam2
Adiabatic indexes , see Eq.~eeosbfpolye}.
double precis
contains the precision required in zerosec_b
virtual int identify() const
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
Definition eos_bf_file.C:95
virtual double press_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the pressure from the baryonic densities and the relative velocity.
virtual void sauve(FILE *) const
Save in a file.
void set_auxiliary()
Computes the auxiliary quantities gam1m1 , gam2m1 and gam3m1.
virtual double ener_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the total energy density from the baryonic densities and the relative velocity.
virtual bool nbar_ent_p(const double ent1, const double ent2, const double delta2, double &nbar1, double &nbar2) const
Computes both baryon densities from the log-enthalpies.
int typeos
The bi-fluid analytical EOS type:
double beta
Coefficient , see Eq.
double kap3
Pressure coefficient , see Eq.
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present.
double ecart
contains the precision required in the relaxation nbar_ent_p
Eos_bf_poly(double kappa1, double kappa2, double kappa3, double beta)
Standard constructor.
double relax
Parameters needed for some inversions of the EOS.
virtual ~Eos_bf_poly()
Destructor.
2-fluids equation of state base class.
virtual void sauve(FILE *) const
Save in a file.
double m_1
Individual particle mass [unit: ].
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
double m_2
Individual particle mass [unit: ].
Polytropic equation of state (relativistic case).
Definition eos.h:757
Equation of state base class.
Definition eos.h:190
Parameter storage.
Definition param.h:125
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition param.C:453
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition param.C:498
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
double zerosec_b(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter)
Finding the zero a function on a bounded domain.
int read_variable(const char *fname, const char *var_name, char *fmt, void *varp)
Reads a variable from file.
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
Lorene prototypes.
Definition app_hor.h:64