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 char et_rot_mag_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag.C,v 1.24 2014/10/13 08:52:58 j_novak Exp $" ;
31 
32 /*
33  * $Id: et_rot_mag.C,v 1.24 2014/10/13 08:52:58 j_novak Exp $
34  * $Log: et_rot_mag.C,v $
35  * Revision 1.24 2014/10/13 08:52:58 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.23 2014/05/13 15:36:54 j_novak
39  * *** empty log message ***
40  *
41  * Revision 1.22 2014/05/13 10:06:13 j_novak
42  * 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.
43  *
44  * Revision 1.21 2013/11/25 13:52:11 j_novak
45  * New class Et_magnetisation to include magnetization terms in the stress energy tensor.
46  *
47  * Revision 1.20 2013/11/14 16:12:55 j_novak
48  * Corrected a mistake in the units.
49  *
50  * Revision 1.19 2012/08/12 17:48:35 p_cerda
51  * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
52  *
53  * Revision 1.18 2011/10/06 14:55:36 j_novak
54  * equation_of_state() is now virtual to be able to call to the magnetized
55  * Eos_mag.
56  *
57  * Revision 1.17 2005/06/02 11:35:30 j_novak
58  * Added members for sving to a file and reading from it.
59  *
60  * Revision 1.16 2004/03/25 10:29:06 j_novak
61  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
62  *
63  * Revision 1.15 2002/09/30 14:21:21 j_novak
64  * *** empty log message ***
65  *
66  * Revision 1.14 2002/09/30 08:50:37 j_novak
67  * Output for central magnetic field value
68  *
69  * Revision 1.13 2002/06/05 15:15:59 j_novak
70  * The case of non-adapted mapping is treated.
71  * parmag.d and parrot.d have been merged.
72  *
73  * Revision 1.12 2002/06/03 13:00:45 e_marcq
74  *
75  * conduc parameter read in parmag.d
76  *
77  * Revision 1.10 2002/05/22 12:20:17 j_novak
78  * *** empty log message ***
79  *
80  * Revision 1.9 2002/05/20 15:44:55 e_marcq
81  *
82  * Dimension errors corrected, parmag.d input file created and read
83  *
84  * Revision 1.8 2002/05/20 08:27:59 j_novak
85  * *** empty log message ***
86  *
87  * Revision 1.7 2002/05/17 15:08:01 e_marcq
88  *
89  * Rotation progressive plug-in, units corrected, Q and a_j new member data
90  *
91  * Revision 1.6 2002/05/16 13:27:11 j_novak
92  * *** empty log message ***
93  *
94  * Revision 1.5 2002/05/16 11:54:11 j_novak
95  * Fixed output pbs
96  *
97  * Revision 1.4 2002/05/15 09:53:59 j_novak
98  * First operational version
99  *
100  * Revision 1.3 2002/05/14 13:38:36 e_marcq
101  *
102  *
103  * Unit update, new outputs
104  *
105  * Revision 1.1 2002/05/10 09:26:52 j_novak
106  * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
107  *
108  *
109  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag.C,v 1.24 2014/10/13 08:52:58 j_novak Exp $
110  *
111  */
112 
113 // Headers C
114 #include "cmath"
115 
116 // Headers Lorene
117 #include "et_rot_mag.h"
118 #include "utilitaires.h"
119 #include "unites.h"
120 #include "param.h"
121 #include "eos.h"
122 
123  //--------------//
124  // Constructors //
125  //--------------//
126 // Standard constructor
127 // --------------------
128 
129 
130 namespace Lorene {
131 Et_rot_mag::Et_rot_mag(Map& mp_i, int nzet_i, bool relat, const Eos& eos_i,
132  const int cond)
133  : Etoile_rot(mp_i, nzet_i, relat, eos_i),
134  A_t(mp_i),
135  A_phi(mp_i),
136  B_phi(mp_i),
137  j_t(mp_i),
138  j_phi(mp_i),
139  E_em(mp_i),
140  Jp_em(mp_i),
141  Srr_em(mp_i),
142  Spp_em(mp_i)
143 
144 {
145 
146  A_t = 0;
147  A_phi = 0;
148  B_phi = 0;
149  j_t = 0 ;
150  j_phi = 0 ;
151 
152  Q = 0 ;
153  a_j = 0 ;
154  conduc = cond ;
155 
156 set_der_0x0() ;
157 }
158 
159 
160 
161 Et_rot_mag::Et_rot_mag(Map& mp_i, const Eos& eos_i, FILE* fich, int withbphi)
162  : Etoile_rot(mp_i, eos_i, fich),
163  A_t(mp_i),
164  A_phi(mp_i),
165  B_phi(mp_i),
166  j_t(mp_i),
167  j_phi(mp_i),
168  E_em(mp_i),
169  Jp_em(mp_i),
170  Srr_em(mp_i),
171  Spp_em(mp_i)
172 {
173 
174  // Etoile parameters
175  // -----------------
176 
177  fread_be(&conduc, sizeof(int), 1, fich) ;
178  fread_be(&Q, sizeof(double), 1, fich) ;
179  fread_be(&a_j, sizeof(double), 1, fich) ;
180 
181  // Read of the saved fields:
182  // ------------------------
183 
184  Cmp A_t_file(mp, *mp.get_mg(), fich) ;
185  A_t = A_t_file ;
186 
187  Cmp A_phi_file(mp, *mp.get_mg(), fich) ;
188  A_phi = A_phi_file ;
189 
190  Cmp j_t_file(mp, *mp.get_mg(), fich) ;
191  j_t = j_t_file ;
192 
193  Cmp j_phi_file(mp, *mp.get_mg(), fich) ;
194  j_phi = j_phi_file ;
195 
196  Tenseur E_em_file(mp, fich) ;
197  E_em = E_em_file ;
198 
199  Tenseur Jp_em_file(mp, fich) ;
200  Jp_em = Jp_em_file ;
201 
202  Tenseur Srr_em_file(mp, fich) ;
203  Srr_em = Srr_em_file ;
204 
205  Tenseur Spp_em_file(mp, fich) ;
206  Spp_em = Spp_em_file ;
207 
208  if ( withbphi == 0 ) {
209  B_phi = 0;
210  }else{
211  Cmp B_phi_file(mp, *mp.get_mg(), fich) ;
212  B_phi = B_phi_file ;
213  }
214 
215  // Pointers of derived quantities initialized to zero
216  // --------------------------------------------------
217  set_der_0x0() ;
218 
219 }
220 
221 
222 // Copy constructor
223 // ----------------
224 
226  : Etoile_rot(et),
227  A_t(et.A_t),
228  A_phi(et.A_phi),
229  B_phi(et.B_phi),
230  j_t(et.j_t),
231  j_phi(et.j_phi),
232  E_em(et.E_em),
233  Jp_em(et.Jp_em),
234  Srr_em(et.Srr_em),
235  Spp_em(et.Spp_em)
236 
237 {
238  Q = et.Q ;
239  a_j = et.a_j ;
240  conduc = et.conduc ;
241  set_der_0x0() ;
242 }
243 
244 
245  //------------//
246  // Destructor //
247  //------------//
248 
250  del_deriv() ;
251 }
252 
253 
254  //----------------------------------//
255  // Management of derived quantities //
256  //----------------------------------//
257 
258 void Et_rot_mag::del_deriv() const {
259 
261 
262  set_der_0x0() ;
263 
264 }
265 
266 
269 
270 }
271 
272 
275 
276  del_deriv() ;
277 }
278 
279 
280 // Assignment to another Et_rot_mag
281 // --------------------------------
282 
284 
285  // Assignement of quantities common to all the derived classes of Etoile
287  A_t = et.A_t ;
288  A_phi = et.A_phi ;
289  B_phi = et.B_phi ;
290  j_t = et.j_t ;
291  j_phi = et.j_phi ;
292  E_em = et.E_em ;
293  Jp_em = et.Jp_em ;
294  Srr_em = et.Srr_em ;
295  Spp_em = et.Spp_em ;
296  Q = et.Q ;
297  a_j = et.a_j ;
298  conduc = et.conduc ;
299 
300  del_deriv() ;
301 
302 }
303 
304 
305  //--------------//
306  // Outputs //
307  //--------------//
308 
309 // Save in a file
310 // --------------
311 void Et_rot_mag::sauve(FILE* fich) const {
312 
313  Etoile_rot::sauve(fich) ;
314 
315  fwrite_be(&conduc, sizeof(int), 1, fich) ;
316  fwrite_be(&Q, sizeof(double), 1, fich) ;
317  fwrite_be(&a_j, sizeof(double), 1, fich) ;
318 
319  A_t.sauve(fich) ;
320  A_phi.sauve(fich) ;
321  j_t.sauve(fich) ;
322  j_phi.sauve(fich) ;
323  E_em.sauve(fich) ;
324  Jp_em.sauve(fich) ;
325  Srr_em.sauve(fich) ;
326  Spp_em.sauve(fich) ;
327  B_phi.sauve(fich) ;
328 
329 }
330 
331 
332 // Printing
333 // --------
334 
335 
336 ostream& Et_rot_mag::operator>>(ostream& ost) const {
337 
338  using namespace Unites_mag ;
339 
341  int theta_eq = mp.get_mg()->get_nt(nzet-1)-1 ;
342  ost << endl ;
343  ost << "Electromagnetic quantities" << endl ;
344  ost << "----------------------" << endl ;
345  ost << endl ;
346  if (is_conduct()) {
347  ost << "***************************" << endl ;
348  ost << "**** perfect conductor ****" << endl ;
349  ost << "***************************" << endl ;
350  ost << "Prescribed charge : " << Q*j_unit*pow(r_unit,3)/v_unit
351  << " [C]" << endl;
352  }
353  else {
354  ost << "***************************" << endl ;
355  ost << "**** isolator ****" << endl ;
356  ost << "***************************" << endl ;
357  }
358  ost << "Prescribed current amplitude : " << a_j*j_unit
359  << " [A/m2]" << endl ;
360  ost << "Magnetic Momentum : " << MagMom()/1.e32
361  << " [10^32 Am^2]" << endl ;
362  ost << "Radial magnetic field polar value : " <<
363  Magn()(0).va.val_point(l_surf()(0,0),xi_surf()(0,0),0.,0.)*mag_unit/1.e9
364  << " [GT]" << endl;
365 
366  ost << "Tangent magnetic field equatorial value : " <<
367  Magn()(1).va.val_point(l_surf()(0,theta_eq),xi_surf()(0,theta_eq),M_PI_2,0.)
368  *mag_unit/1.e9
369  << " [GT]" << endl;
370 
371  ost << "Central magnetic field values : " <<
372  Magn()(0)(0,0,0,0)*mag_unit/1.e9
373  << " [GT]" << endl;
374 
375  ost << "Radial electric field polar value : " <<
376  Elec()(0).va.val_point(l_surf()(0,0),xi_surf()(0,0),0.,0.)
377  *elec_unit/1.e12
378  << " [TV]" << endl;
379 
380  ost << "Radial electric field equatorial value : " <<
381  Elec()(0).va.val_point(l_surf()(0,theta_eq),xi_surf()(0,theta_eq),M_PI_2,0.)
382  *elec_unit/1.e12
383  << " [TV]" << endl;
384 
385  ost << "Magnetic/fluid pressure : "
386  << 1/(2*mu_si)*(pow(Magn()(0)(0,0,0,0),2) +
387  pow(Magn()(1)(0,0,0,0),2) +
388  pow(Magn()(2)(0,0,0,0),2))*mag_unit*mag_unit
389  / (press()(0,0,0,0)*rho_unit*pow(v_unit,2)) << endl ;
390  ost << "Computed charge : " << Q_comput() << " [C]" << endl ;
391  ost << "Interior charge : " << Q_int() << " [C]" << endl ;
392  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)
393  << endl ;
394  ost << "Gyromagnetic ratio : " << GyroMag() << endl ;
395 
396  return ost ;
397 }
398 
399 
400 
401 
402 
403 
404 
405 
406 
407 
408 
409 
410 
411 
412 
413 
414 }
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:64
Equation of state base class.
Definition: eos.h:190
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Base class for coordinate mappings.
Definition: map.h:670
Class for isolated rotating stars *** DEPRECATED : use class Star_rot instead ***.
Definition: etoile.h:1496
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:404
virtual ~Et_rot_mag()
Destructor.
Definition: et_rot_mag.C:249
Tenseur press
Fluid pressure.
Definition: etoile.h:461
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...
Definition: et_rot_global.C:95
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
Definition: et_rot_mag.C:273
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:335
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:466
void operator=(const Et_rot_mag &)
Assignment to another Et_rot_mag.
Definition: et_rot_mag.C:283
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:1325
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:336
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile_rot.C:441
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
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:131
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: et_rot_mag.C:258
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
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
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:432
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
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:311
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: et_rot_mag.C:267
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:364
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
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:561
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:389
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:301
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