LORENE
et_rot_mag.C
1 /*
2  * Methods for magnetized axisymmetric rotating neutron stars.
3  *
4  * See the file et_rot_mag.h for documentation
5  *
6  */
7 
8 /*
9  * Copyright (c) 2002 Emmanuel Marcq
10  * Copyright (c) 2002 Jerome Novak
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 /*
33  * $Id: et_rot_mag.C,v 1.25 2016/12/05 16:17:54 j_novak Exp $
34  * $Log: et_rot_mag.C,v $
35  * Revision 1.25 2016/12/05 16:17:54 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.24 2014/10/13 08:52:58 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.23 2014/05/13 15:36:54 j_novak
42  * *** empty log message ***
43  *
44  * Revision 1.22 2014/05/13 10:06:13 j_novak
45  * Change of magnetic units, to make the Lorene unit system coherent. Magnetic field is now expressed in Lorene units. Improvement on the comments on units.
46  *
47  * Revision 1.21 2013/11/25 13:52:11 j_novak
48  * New class Et_magnetisation to include magnetization terms in the stress energy tensor.
49  *
50  * Revision 1.20 2013/11/14 16:12:55 j_novak
51  * Corrected a mistake in the units.
52  *
53  * Revision 1.19 2012/08/12 17:48:35 p_cerda
54  * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
55  *
56  * Revision 1.18 2011/10/06 14:55:36 j_novak
57  * equation_of_state() is now virtual to be able to call to the magnetized
58  * Eos_mag.
59  *
60  * Revision 1.17 2005/06/02 11:35:30 j_novak
61  * Added members for sving to a file and reading from it.
62  *
63  * Revision 1.16 2004/03/25 10:29:06 j_novak
64  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
65  *
66  * Revision 1.15 2002/09/30 14:21:21 j_novak
67  * *** empty log message ***
68  *
69  * Revision 1.14 2002/09/30 08:50:37 j_novak
70  * Output for central magnetic field value
71  *
72  * Revision 1.13 2002/06/05 15:15:59 j_novak
73  * The case of non-adapted mapping is treated.
74  * parmag.d and parrot.d have been merged.
75  *
76  * Revision 1.12 2002/06/03 13:00:45 e_marcq
77  *
78  * conduc parameter read in parmag.d
79  *
80  * Revision 1.10 2002/05/22 12:20:17 j_novak
81  * *** empty log message ***
82  *
83  * Revision 1.9 2002/05/20 15:44:55 e_marcq
84  *
85  * Dimension errors corrected, parmag.d input file created and read
86  *
87  * Revision 1.8 2002/05/20 08:27:59 j_novak
88  * *** empty log message ***
89  *
90  * Revision 1.7 2002/05/17 15:08:01 e_marcq
91  *
92  * Rotation progressive plug-in, units corrected, Q and a_j new member data
93  *
94  * Revision 1.6 2002/05/16 13:27:11 j_novak
95  * *** empty log message ***
96  *
97  * Revision 1.5 2002/05/16 11:54:11 j_novak
98  * Fixed output pbs
99  *
100  * Revision 1.4 2002/05/15 09:53:59 j_novak
101  * First operational version
102  *
103  * Revision 1.3 2002/05/14 13:38:36 e_marcq
104  *
105  *
106  * Unit update, new outputs
107  *
108  * Revision 1.1 2002/05/10 09:26:52 j_novak
109  * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
110  *
111  *
112  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag.C,v 1.25 2016/12/05 16:17:54 j_novak Exp $
113  *
114  */
115 
116 // Headers C
117 #include "cmath"
118 
119 // Headers Lorene
120 #include "et_rot_mag.h"
121 #include "utilitaires.h"
122 #include "unites.h"
123 #include "param.h"
124 #include "eos.h"
125 
126  //--------------//
127  // Constructors //
128  //--------------//
129 // Standard constructor
130 // --------------------
131 
132 
133 namespace Lorene {
134 Et_rot_mag::Et_rot_mag(Map& mp_i, int nzet_i, bool relat, const Eos& eos_i,
135  const int cond)
136  : Etoile_rot(mp_i, nzet_i, relat, eos_i),
137  A_t(mp_i),
138  A_phi(mp_i),
139  B_phi(mp_i),
140  j_t(mp_i),
141  j_phi(mp_i),
142  E_em(mp_i),
143  Jp_em(mp_i),
144  Srr_em(mp_i),
145  Spp_em(mp_i)
146 
147 {
148 
149  A_t = 0;
150  A_phi = 0;
151  B_phi = 0;
152  j_t = 0 ;
153  j_phi = 0 ;
154 
155  Q = 0 ;
156  a_j = 0 ;
157  conduc = cond ;
158 
159 set_der_0x0() ;
160 }
161 
162 
163 
164 Et_rot_mag::Et_rot_mag(Map& mp_i, const Eos& eos_i, FILE* fich, int withbphi)
165  : Etoile_rot(mp_i, eos_i, fich),
166  A_t(mp_i),
167  A_phi(mp_i),
168  B_phi(mp_i),
169  j_t(mp_i),
170  j_phi(mp_i),
171  E_em(mp_i),
172  Jp_em(mp_i),
173  Srr_em(mp_i),
174  Spp_em(mp_i)
175 {
176 
177  // Etoile parameters
178  // -----------------
179 
180  fread_be(&conduc, sizeof(int), 1, fich) ;
181  fread_be(&Q, sizeof(double), 1, fich) ;
182  fread_be(&a_j, sizeof(double), 1, fich) ;
183 
184  // Read of the saved fields:
185  // ------------------------
186 
187  Cmp A_t_file(mp, *mp.get_mg(), fich) ;
188  A_t = A_t_file ;
189 
190  Cmp A_phi_file(mp, *mp.get_mg(), fich) ;
191  A_phi = A_phi_file ;
192 
193  Cmp j_t_file(mp, *mp.get_mg(), fich) ;
194  j_t = j_t_file ;
195 
196  Cmp j_phi_file(mp, *mp.get_mg(), fich) ;
197  j_phi = j_phi_file ;
198 
199  Tenseur E_em_file(mp, fich) ;
200  E_em = E_em_file ;
201 
202  Tenseur Jp_em_file(mp, fich) ;
203  Jp_em = Jp_em_file ;
204 
205  Tenseur Srr_em_file(mp, fich) ;
206  Srr_em = Srr_em_file ;
207 
208  Tenseur Spp_em_file(mp, fich) ;
209  Spp_em = Spp_em_file ;
210 
211  if ( withbphi == 0 ) {
212  B_phi = 0;
213  }else{
214  Cmp B_phi_file(mp, *mp.get_mg(), fich) ;
215  B_phi = B_phi_file ;
216  }
217 
218  // Pointers of derived quantities initialized to zero
219  // --------------------------------------------------
220  set_der_0x0() ;
221 
222 }
223 
224 
225 // Copy constructor
226 // ----------------
227 
229  : Etoile_rot(et),
230  A_t(et.A_t),
231  A_phi(et.A_phi),
232  B_phi(et.B_phi),
233  j_t(et.j_t),
234  j_phi(et.j_phi),
235  E_em(et.E_em),
236  Jp_em(et.Jp_em),
237  Srr_em(et.Srr_em),
238  Spp_em(et.Spp_em)
239 
240 {
241  Q = et.Q ;
242  a_j = et.a_j ;
243  conduc = et.conduc ;
244  set_der_0x0() ;
245 }
246 
247 
248  //------------//
249  // Destructor //
250  //------------//
251 
253  del_deriv() ;
254 }
255 
256 
257  //----------------------------------//
258  // Management of derived quantities //
259  //----------------------------------//
260 
261 void Et_rot_mag::del_deriv() const {
262 
264 
265  set_der_0x0() ;
266 
267 }
268 
269 
272 
273 }
274 
275 
278 
279  del_deriv() ;
280 }
281 
282 
283 // Assignment to another Et_rot_mag
284 // --------------------------------
285 
287 
288  // Assignement of quantities common to all the derived classes of Etoile
290  A_t = et.A_t ;
291  A_phi = et.A_phi ;
292  B_phi = et.B_phi ;
293  j_t = et.j_t ;
294  j_phi = et.j_phi ;
295  E_em = et.E_em ;
296  Jp_em = et.Jp_em ;
297  Srr_em = et.Srr_em ;
298  Spp_em = et.Spp_em ;
299  Q = et.Q ;
300  a_j = et.a_j ;
301  conduc = et.conduc ;
302 
303  del_deriv() ;
304 
305 }
306 
307 
308  //--------------//
309  // Outputs //
310  //--------------//
311 
312 // Save in a file
313 // --------------
314 void Et_rot_mag::sauve(FILE* fich) const {
315 
316  Etoile_rot::sauve(fich) ;
317 
318  fwrite_be(&conduc, sizeof(int), 1, fich) ;
319  fwrite_be(&Q, sizeof(double), 1, fich) ;
320  fwrite_be(&a_j, sizeof(double), 1, fich) ;
321 
322  A_t.sauve(fich) ;
323  A_phi.sauve(fich) ;
324  j_t.sauve(fich) ;
325  j_phi.sauve(fich) ;
326  E_em.sauve(fich) ;
327  Jp_em.sauve(fich) ;
328  Srr_em.sauve(fich) ;
329  Spp_em.sauve(fich) ;
330  B_phi.sauve(fich) ;
331 
332 }
333 
334 
335 // Printing
336 // --------
337 
338 
339 ostream& Et_rot_mag::operator>>(ostream& ost) const {
340 
341  using namespace Unites_mag ;
342 
344  int theta_eq = mp.get_mg()->get_nt(nzet-1)-1 ;
345  ost << endl ;
346  ost << "Electromagnetic quantities" << endl ;
347  ost << "----------------------" << endl ;
348  ost << endl ;
349  if (is_conduct()) {
350  ost << "***************************" << endl ;
351  ost << "**** perfect conductor ****" << endl ;
352  ost << "***************************" << endl ;
353  ost << "Prescribed charge : " << Q*j_unit*pow(r_unit,3)/v_unit
354  << " [C]" << endl;
355  }
356  else {
357  ost << "***************************" << endl ;
358  ost << "**** isolator ****" << endl ;
359  ost << "***************************" << endl ;
360  }
361  ost << "Prescribed current amplitude : " << a_j*j_unit
362  << " [A/m2]" << endl ;
363  ost << "Magnetic Momentum : " << MagMom()/1.e32
364  << " [10^32 Am^2]" << endl ;
365  ost << "Radial magnetic field polar value : " <<
366  Magn()(0).va.val_point(l_surf()(0,0),xi_surf()(0,0),0.,0.)*mag_unit/1.e9
367  << " [GT]" << endl;
368 
369  ost << "Tangent magnetic field equatorial value : " <<
370  Magn()(1).va.val_point(l_surf()(0,theta_eq),xi_surf()(0,theta_eq),M_PI_2,0.)
371  *mag_unit/1.e9
372  << " [GT]" << endl;
373 
374  ost << "Central magnetic field values : " <<
375  Magn()(0)(0,0,0,0)*mag_unit/1.e9
376  << " [GT]" << endl;
377 
378  ost << "Radial electric field polar value : " <<
379  Elec()(0).va.val_point(l_surf()(0,0),xi_surf()(0,0),0.,0.)
380  *elec_unit/1.e12
381  << " [TV]" << endl;
382 
383  ost << "Radial electric field equatorial value : " <<
384  Elec()(0).va.val_point(l_surf()(0,theta_eq),xi_surf()(0,theta_eq),M_PI_2,0.)
385  *elec_unit/1.e12
386  << " [TV]" << endl;
387 
388  ost << "Magnetic/fluid pressure : "
389  << 1/(2*mu_si)*(pow(Magn()(0)(0,0,0,0),2) +
390  pow(Magn()(1)(0,0,0,0),2) +
391  pow(Magn()(2)(0,0,0,0),2))*mag_unit*mag_unit
392  / (press()(0,0,0,0)*rho_unit*pow(v_unit,2)) << endl ;
393  ost << "Computed charge : " << Q_comput() << " [C]" << endl ;
394  ost << "Interior charge : " << Q_int() << " [C]" << endl ;
395  ost << "Q^2/M^2 :" << mu_si*c_si*c_si*Q_comput()*Q_comput()/(4*M_PI*g_si*mass_g()*mass_g()*m_unit*m_unit)
396  << endl ;
397  ost << "Gyromagnetic ratio : " << GyroMag() << endl ;
398 
399  return ost ;
400 }
401 
402 
403 
404 
405 
406 
407 
408 
409 
410 
411 
412 
413 
414 
415 
416 
417 }
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinate...
Tenseur E_em
electromagnetic energy density in the Eulerian frame
Definition: et_rot_mag.h:161
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
double GyroMag() const
Gyromagnetic ratio .
double Q
In the case of a perfect conductor, the requated baryonic charge.
Definition: et_rot_mag.h:179
Lorene prototypes.
Definition: app_hor.h:67
Equation of state base class.
Definition: eos.h:209
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Base class for coordinate mappings.
Definition: map.h:688
Class for isolated rotating stars *** DEPRECATED : use class Star_rot instead ***.
Definition: etoile.h:1499
virtual double mass_g() const
Gravitational mass.
Tenseur Magn() const
Computes the magnetic field spherical components in Lorene&#39;s units.
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
Definition: etoile_rot.C:415
virtual ~Et_rot_mag()
Destructor.
Definition: et_rot_mag.C:252
Tenseur press
Fluid pressure.
Definition: etoile.h:464
virtual const Itbl & l_surf() const
Description of the stellar surface: returns a 2-D Itbl containing the values of the domain index l on...
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
Definition: et_rot_mag.C:276
Cmp A_phi
-component of the electromagnetic potential 1-form divided by .
Definition: et_rot_mag.h:155
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: etoile_rot.C:344
Tenseur Srr_em
rr component of the electromagnetic stress 3-tensor, as measured in the Eulerian frame. (not used and set to 0, should be supressed)
Definition: et_rot_mag.h:170
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: etoile_rot.C:477
void operator=(const Et_rot_mag &)
Assignment to another Et_rot_mag.
Definition: et_rot_mag.C:286
double a_j
Amplitude of the curent/charge function.
Definition: et_rot_mag.h:180
void sauve(FILE *) const
Save in a file.
Definition: tenseur.C:1341
bool is_conduct() const
Tells if the star is made of conducting or isolating material.
Definition: et_rot_mag.h:241
Class for magnetized (isolator or perfect conductor), rigidly rotating stars.
Definition: et_rot_mag.h:143
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: et_rot_mag.C:339
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile_rot.C:452
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
Et_rot_mag(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i, const int cond)
Standard constructor.
Definition: et_rot_mag.C:134
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: et_rot_mag.C:261
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:73
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Cmp j_phi
-component of the current 4-vector
Definition: et_rot_mag.h:159
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:435
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:72
Cmp j_t
t-component of the current 4-vector
Definition: et_rot_mag.h:158
double Q_int() const
Computed charge from the integration of charge density over the star (i.e.
virtual void sauve(FILE *) const
Save in a file.
Definition: et_rot_mag.C:314
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: et_rot_mag.C:270
Tenseur Elec() const
Computes the electric field spherical components in Lorene&#39;s units.
Cmp A_t
t-component of the elecctromagnetic potential 1-form, divided by .
Definition: et_rot_mag.h:150
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: etoile_rot.C:374
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
int conduc
Flag: conduc=0->isolator, 1->perfect conductor.
Definition: et_rot_mag.h:181
void sauve(FILE *) const
Save in a file.
Definition: cmp.C:564
double Q_comput() const
Computed charge deduced from the asymptotic behaviour of At [SI units].
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
Definition: etoile_rot.C:400
double MagMom() const
Magnetic Momentum in SI units.
Cmp B_phi
-component of the magnetic field
Definition: et_rot_mag.h:157
Standard electro-magnetic units.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Tenseur Jp_em
component of the electromagnetic momentum density 3-vector, as measured in the Eulerian frame...
Definition: et_rot_mag.h:167
Tenseur Spp_em
component of the electromagnetic stress 3-tensor, as measured in the Eulerian frame.
Definition: et_rot_mag.h:173