LORENE
eos_mag.C
1 /* Methods of class Eos_mag
2 *
3 * (see file eos_mag.h for documentation).
4 *
5 */
6
7/*
8 * Copyright (c) 2011 Thomas Elghozi & Jerome Novak
9 * Copyright (c) 2013 Debarati Chatterjee
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30char eos_mag_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $" ;
31
32/*
33 * $Id: eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $
34 * $Log: eos_mag.C,v $
35 * Revision 1.13 2014/10/13 08:52:53 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.12 2014/09/22 16:13:10 j_novak
39 * Minor modif.
40 *
41 * Revision 1.11 2014/05/27 12:32:28 j_novak
42 * Added possibility to converge to a given magnetic moment.
43 *
44 * Revision 1.10 2014/05/13 15:37:12 j_novak
45 * Updated to new magnetic units.
46 *
47 * Revision 1.9 2014/04/28 12:48:13 j_novak
48 * Minor modifications.
49 *
50 * Revision 1.8 2014/03/11 14:27:26 j_novak
51 * Corrected a missing 4pi term.
52 *
53 * Revision 1.7 2014/03/03 16:23:08 j_novak
54 * Updated error message
55 *
56 * Revision 1.6 2013/12/12 16:07:30 j_novak
57 * interpol_herm_2d outputs df/dx, used to get the magnetization.
58 *
59 * Revision 1.5 2013/11/25 15:00:52 j_novak
60 * Correction of memory error.
61 *
62 * Revision 1.4 2013/11/14 16:12:55 j_novak
63 * Corrected a mistake in the units.
64 *
65 * Revision 1.2 2011/10/04 16:05:19 j_novak
66 * Update of Eos_mag class. Suppression of loge, re-definition of the derivatives
67 * and use of interpol_herm_2d.
68 *
69 * Revision 1.1 2011/06/16 10:49:18 j_novak
70 * New class Eos_mag for EOSs depending on density and magnetic field.
71 *
72 *
73 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $
74 *
75 */
76
77// headers C
78#include <cmath>
79
80// Headers Lorene
81#include "eos.h"
82#include "cmp.h"
83#include "param.h"
84#include "utilitaires.h"
85#include "unites.h"
86
87
88namespace Lorene {
89void interpol_herm_2d(const Tbl&, const Tbl&, const Tbl&, const Tbl&, const Tbl&, const Tbl&, double, double, double&, double&, double&) ;
90
91
92 //----------------------------//
93 // Constructors //
94 //----------------------------//
95
96// Standard constructor
97// --------------------
98Eos_mag::Eos_mag(const char* name_i, const char* table,
99 const char* path) : Eos(name_i), tablename(path) {
100
101 tablename += '/' ;
102 tablename += table ;
103
104 read_table() ;
105
106}
107
108// Standard constructor with full filename
109// ---------------------------------------
110Eos_mag::Eos_mag(const char* name_i, const char* file_name)
111 : Eos(name_i), tablename(file_name) {
112
113 read_table() ;
114
115}
116
117
118// Constructor from binary file
119// ----------------------------
120Eos_mag::Eos_mag(FILE* fich) : Eos(fich) {
121
122 char tmp_name[160] ;
123
124 fread(tmp_name, sizeof(char), 160, fich) ;
125 tablename += tmp_name ;
126
127 read_table() ;
128
129}
130
131
132
133// Constructor from a formatted file
134// ---------------------------------
135Eos_mag::Eos_mag(ifstream& fich, const char* table) : Eos(fich) {
136
137 fich >> tablename ;
138 tablename += '/' ;
139 tablename += table ;
140
141 read_table() ;
142
143}
144
145Eos_mag::Eos_mag(ifstream& fich) : Eos(fich) {
146
147 fich >> tablename ;
148
149 read_table() ;
150
151}
152
153
154 //--------------//
155 // Destructor //
156 //--------------//
157
159 delete d2lp ;
160 delete dlpsdB ;
161 delete dlpsdlh ;
162 delete Bfield ;
163 delete logh ;
164 delete logp ;
165}
166
167 //------------//
168 // Outputs //
169 //------------//
170
171void Eos_mag::sauve(FILE* fich) const {
172
173 Eos::sauve(fich) ;
174
175 fwrite(tablename.c_str(), sizeof(char), tablename.size(), fich) ;
176
177}
178 //------------------------//
179 // Comparison operators //
180 //------------------------//
181
182
183bool Eos_mag::operator==(const Eos& eos_i) const {
184
185 bool resu = true ;
186
187 if ( eos_i.identify() != identify() ) {
188 cerr << "The second EOS is not of type Eos_mag !" << endl ;
189 resu = false ;
190 }
191
192 return resu ;
193
194}
195
196bool Eos_mag::operator!=(const Eos& eos_i) const {
197
198 return !(operator==(eos_i)) ;
199
200}
201
202 //------------//
203 // Outputs //
204 //------------//
205
206
207ostream& Eos_mag::operator>>(ostream & ost) const {
208
209 ost <<
210 "EOS of class Eos_mag : tabulated EOS depending on two variables: enthalpy and magnetic field."
211 << '\n' ;
212 ost << "Table read from " << tablename << endl ;
213
214 return ost ;
215
216}
217
218
219
220 //------------------------//
221 // Reading of the table //
222 //------------------------//
223
225
226 using namespace Unites_mag ;
227
228 ifstream fich(tablename.data()) ;
229
230 if (!fich) {
231 cerr << "Eos_mag::read_table(): " << endl ;
232 cerr << "Problem in opening the EOS file!" << endl ;
233 cerr << "Aborting..." << endl ;
234 abort() ;
235 }
236
237 for (int i=0; i<5; i++) { // jump over the file
238 fich.ignore(1000, '\n') ; // header
239 } //
240
241 int nbp1, nbp2 ;
242 fich >> nbp1 >> nbp2 ; fich.ignore(1000, '\n') ; // number of data
243 if ((nbp1<=0) || (nbp2<=0)) {
244 cerr << "Eos_mag::read_table(): " << endl ;
245 cerr << "Wrong value for the number of lines!" << endl ;
246 cerr << "nbp1 = " << nbp1 << ", nbp2 = " << nbp2 << endl ;
247 cerr << "Aborting..." << endl ;
248 abort() ;
249 }
250
251 for (int i=0; i<3; i++) { // jump over the table
252 fich.ignore(1000, '\n') ;
253 }
254
255 logp = new Tbl(nbp2, nbp1) ;
256 logh = new Tbl(nbp2, nbp1) ;
257 Bfield = new Tbl(nbp2, nbp1) ;
258 dlpsdlh = new Tbl(nbp2, nbp1) ;
259 dlpsdB = new Tbl(nbp2, nbp1) ;
260 d2lp = new Tbl(nbp2, nbp1) ;
261
262 logp->set_etat_qcq() ;
263 logh->set_etat_qcq() ;
267 d2lp->set_etat_qcq() ;
268
269 double c2 = c_si * c_si ;
270 double B_unit = mag_unit / 1.e11 ;
271 double M_unit = B_unit*mu0/(4*M_PI) ;
272
273 int no1 ;
274 double nb_fm3, rho_cgs, p_cgs, mu_MeV, magB_PG, magM_PG, chi_PGpMeV ;
275
276 double ww = 0. ;
277
278 for (int j=0; j<nbp2; j++) {
279
280 for (int i=0; i<nbp1; i++) {
281 fich >> no1 >> nb_fm3 >> rho_cgs >> p_cgs >> mu_MeV
282 >> magB_PG >> magM_PG >> chi_PGpMeV ;
283 fich.ignore(1000,'\n') ;
284 if ( (nb_fm3<0) || (rho_cgs<0)) { // || (p_cgs < 0) ){
285 cerr << "Eos_mag::read_table(): " << endl ;
286 cerr << "Negative value in table!" << endl ;
287 cerr << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
288 ", p = " << p_cgs << endl ;
289 cerr << "Aborting..." << endl ;
290 abort() ;
291 }
292
293 double psc2 = 0.1*p_cgs/c2 ; // in kg m^-3
294 double rho_si = rho_cgs*1000. ; // in kg m^-3
295
296 double h_read = log(mu_MeV) ;
297 if ( (i==0) && (j==0) ) ww = h_read ;
298 double h_new = h_read - ww ;
299
300 logp->set(j, i) = psc2/rhonuc_si ;
301 logh->set(j, i) = h_new ;
302 Bfield->set(j, i) = magB_PG / B_unit ; // in Lorene units
303 dlpsdlh->set(j, i) = (rho_si + psc2)/rhonuc_si ;
304 dlpsdB->set(j, i) = magM_PG / M_unit ;
305 d2lp->set(j, i) = mu_MeV*chi_PGpMeV / (M_unit) ;
306
307 }
308 }
309
310 hmin = (*logh)(0, 0) ;
311 hmax = (*logh)(0, nbp1-1) ;
312
313 Bmax = (*Bfield)(nbp2-1, 0) ;
314
315 // cout << "hmin: " << hmin << ", hmax: " << hmax << endl ;
316 // cout << "Bmax: " << Bmax << endl ;
317
318 fich.close();
319
320}
321
322
323 //------------------------------//
324 // Computational routines //
325 //------------------------------//
326
327// Baryon density from enthalpy
328//------------------------------
329
330double Eos_mag::nbar_ent_p(double ent, const Param* par ) const {
331
332 using namespace Unites_mag ;
333
334 if ((ent > hmin - 1.e-12) && (ent < hmin))
335 ent = hmin ;
336
337 if ( ent >= hmin) {
338 if (ent > hmax) {
339 cerr << "Eos_tabul::nbar_ent_p : ent > hmax !" << endl ;
340 abort() ;
341 }
342 // recuperer magB0 (input)
343 double magB0 = 0. ;
344 if (par != 0x0)
345 if (par->get_n_double_mod() > 0) {
346 magB0 = par->get_double_mod() ;
347 }
348
349 if ( magB0 > Bmax) {
350 cerr << "Eos_tabul::nbar_ent_p : B > Bmax !" << endl ;
351 cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
352 abort() ;
353 }
354
355 double p_int, dpdb_int, dpdh_int ;
356
357 interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
358 p_int, dpdb_int, dpdh_int) ;
359
360 double nbar_int = dpdh_int * exp(-ent) ;
361
362 return nbar_int ;
363 }
364 else{
365 return 0 ;
366 }
367}
368
369
370// Energy density from enthalpy
371//------------------------------
372
373double Eos_mag::ener_ent_p(double ent, const Param* par ) const {
374
375 using namespace Unites_mag ;
376
377 if ((ent > hmin - 1.e-12) && (ent < hmin))
378 ent = hmin ;
379
380 if ( ent >= hmin ) {
381 if (ent > hmax) {
382 cerr << "Eos_tabul::ener_ent_p : ent > hmax !" << endl ;
383 abort() ;
384 }
385
386 // recuperer magB0 (input)
387 double magB0 = 0. ;
388 if (par != 0x0)
389 if (par->get_n_double_mod() > 0) {
390 magB0 = par->get_double_mod() ;
391 }
392
393 if ( magB0 > Bmax) {
394 cerr << "Eos_tabul::ener_ent_p : B > Bmax !" << endl ;
395 cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
396 abort() ;
397 }
398
399
400 double p_int, dpdb_int, dpdh_int ;
401
402 interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
403 p_int, dpdb_int, dpdh_int) ;
404
405 double nbar_int = dpdh_int * exp(-ent) ;
406
407 double f_int = - p_int + exp(ent) * nbar_int;
408
409 return f_int ;
410 }
411 else{
412 return 0 ;
413 }
414}
415
416// Pressure from enthalpy
417//------------------------
418
419double Eos_mag::press_ent_p(double ent, const Param* par ) const {
420
421 using namespace Unites_mag ;
422
423 if ((ent > hmin - 1.e-12) && (ent < hmin))
424 ent = hmin ;
425
426 if ( ent >= hmin ) {
427 if (ent > hmax) {
428 cout << "Eos_mag::press_ent_p : ent > hmax !" << endl ;
429 abort() ;
430 }
431
432 // recuperer magB0 (input)
433 double magB0 = 0. ;
434 if (par != 0x0)
435 if (par->get_n_double_mod() > 0) {
436 magB0 = par->get_double_mod() ;
437 }
438
439 if ( magB0 > Bmax) {
440 cerr << "Eos_tabul::press_ent_p : B > Bmax !" << endl ;
441 cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
442 abort() ;
443 }
444
445 double p_int, dpdb_int, dpdh_int ;
446
447 interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
448 p_int, dpdb_int, dpdh_int) ;
449
450 return p_int;
451 }
452 else{
453 return 0 ;
454 }
455}
456
457// Magnetization from enthalpy
458//----------------------------
459
460double Eos_mag::mag_ent_p(double ent, const Param* par) const {
461
462 using namespace Unites_mag ;
463
464 if ((ent > hmin - 1.e-12) && (ent < hmin))
465 ent = hmin ;
466
467 if ( ent >= hmin ) {
468 if (ent > hmax) {
469 cout << "Eos_mag::mag_ent_p : ent > hmax !" << endl ;
470 abort() ;
471 }
472
473 // recuperer magB0 (input)
474 double magB0 = 0. ;
475 if (par != 0x0)
476 if (par->get_n_double_mod() > 0) {
477 magB0 = par->get_double_mod() ;
478 }
479
480 if ( magB0 > Bmax) {
481 cerr << "Eos_tabul::mag_ent_p : B > Bmax !" << endl ;
482 cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
483 abort() ;
484 }
485
486 double p_int, dpdb_int, dpdh_int ;
487
488 interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
489 p_int, dpdb_int, dpdh_int) ;
490
491 double magnetization = dpdb_int ;
492
493 if (magB0 == 0.)
494 return 0 ;
495 else
496 return mu0*magnetization / magB0 ;
497 }
498
499 else
500 return 0. ;
501
502}
503
504
505// dln(n)/ln(H) from enthalpy
506//---------------------------
507
508double Eos_mag::der_nbar_ent_p(double ent, const Param* ) const {
509
510 c_est_pas_fait("Eos_mag::der_nbar_ent_p" ) ;
511
512 return ent ;
513
514}
515
516
517// dln(e)/ln(H) from enthalpy
518//---------------------------
519
520double Eos_mag::der_ener_ent_p(double ent, const Param* ) const {
521
522
523 c_est_pas_fait("Eos_mag::der_ener_enr_p" ) ;
524
525 return ent ;
526
527}
528
529
530// dln(p)/ln(H) from enthalpy
531//---------------------------
532
533double Eos_mag::der_press_ent_p(double ent, const Param* ) const {
534
535 c_est_pas_fait("Eos_mag::" ) ;
536
537 return ent ;
538
539}
540
541
542// dln(p)/dln(n) from enthalpy
543//---------------------------
544
545double Eos_mag::der_press_nbar_p(double ent, const Param*) const {
546
547 c_est_pas_fait("Eos_mag::der_press_nbar_p" ) ;
548
549 return ent ;
550
551}
552}
Tbl * dlpsdB
Table of .
Definition eos_mag.h:109
double mag_ent_p(double ent, const Param *par=0x0) const
Computes the magnetisation.
Definition eos_mag.C:460
Tbl * d2lp
Table of .
Definition eos_mag.h:112
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Definition eos_mag.C:196
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Tbl * logp
Table of .
Definition eos_mag.h:97
Tbl * logh
Table of .
Definition eos_mag.h:100
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:533
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition eos_mag.C:419
virtual void sauve(FILE *) const
Save in a file.
Definition eos_mag.C:171
double Bmax
Upper boundary of the magnetic field interval.
Definition eos_mag.h:94
virtual ~Eos_mag()
Destructor.
Definition eos_mag.C:158
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition eos_mag.C:330
Tbl * Bfield
Table of .
Definition eos_mag.h:103
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:508
virtual bool operator==(const Eos &) const
Comparison operator (egality)
Definition eos_mag.C:183
double hmin
Lower boundary of the log-enthalpy interval.
Definition eos_mag.h:88
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:520
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition eos_mag.C:373
string tablename
Name of the file containing the tabulated data.
Definition eos_mag.h:85
void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh .
Definition eos_mag.C:224
Tbl * dlpsdlh
Table of .
Definition eos_mag.h:106
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition eos_mag.C:545
Eos_mag(const char *name_i, const char *table, const char *path)
Standard constructor.
Definition eos_mag.C:98
double hmax
Upper boundary of the log-enthalpy interval.
Definition eos_mag.h:91
virtual ostream & operator>>(ostream &) const
Operator >>
Definition eos_mag.C:207
Equation of state base class.
Definition eos.h:190
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual void sauve(FILE *) const
Save in a file.
Definition eos.C:179
Parameter storage.
Definition param.h:125
int get_n_double_mod() const
Returns the number of stored modifiable double 's addresses.
Definition param.C:446
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition param.C:498
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Lorene prototypes.
Definition app_hor.h:64
Standard electro-magnetic units.