ex_eos_had_apr.cpp
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file ex_eos_had_apr.cpp
24  \brief File defining \ref ex_eos_had_apr class
25 */
26 /* Example: ex_eos_had_apr.cpp
27  -------------------------------------------------------------------
28 */
29 #include <iostream>
30 #include <cmath>
31 #include <boost/numeric/ublas/vector.hpp>
32 #include <o2scl/eos_had_apr.h>
33 #include <o2scl/table.h>
34 #include <o2scl/constants.h>
35 #include <o2scl/fermion_nonrel.h>
36 #include <o2scl/tov_solve.h>
37 #include <o2scl/deriv_cern.h>
38 #include <o2scl/mroot_cern.h>
39 #include <o2scl/test_mgr.h>
40 #include <o2scl/hdf_io.h>
41 
42 using namespace std;
43 using namespace o2scl;
44 using namespace o2scl_hdf;
45 using namespace o2scl_const;
46 using namespace o2scl_cgs;
47 
50 
51 /** \brief Compute the APR EOS with a Gibbs construction and the mass
52  versus radius curve [Example class]
53 
54  In succession, calculates nuclear matter, neutron matter, and then
55  neutron star matter with Maxwell and Gibbs constructions.
56 
57  We could use the more accurate masses in
58  <tt>o2scl/constants.h</tt> here, but APR appears to have been
59  designed to be used with neutron and protons masses equal
60  to 939 MeV.
61 */
63 
64 protected:
65 
66  /// \name Fermions
67  //@{
68  /// Compute zero-temperature thermodynamics
70  /// Neutron for low-density phase
72  /// Proton for low-density phase
74  /// Neutron for high-density phase
76  /// Proton for high-density phase
78  /// Electron for low-density phase
80  /// Muon for low-density phase
82  /// Electron for high-density phase
84  /// Muon for high-density phase
86  //@}
87 
88  /// \name 'Thermo' objects
89  //@{
90  /// Baryon thermodynamics for low-density phase
92  /// Leptonic thermodynamics for low-density phase
94  /// Baryon thermodynamics for high-density phase
96  /// Total thermodynamics
98  /// Leptonic thermodynamics for high-density phase
100  //@}
101 
102  /// \name Numerical methods
103  //@{
104  /// General solver
106  /// Solver for transition densities (lower tolerances)
108  /// Derivative object
110  //@}
111 
112  /// Baryon density
113  double nb;
114  /// Volume fraction of low-density phase
115  double chi;
116  /// Baryon chemical potential
117  double mub;
118  /// Charge chemical potential
119  double muq;
120  /// Proton fraction for Fig. 7
121  double f7x;
122  /// Choice of model from APR
123  int choice;
124  /// \name Phase specification
125  //@{
126  int phase;
127  static const int low_phase=1;
128  static const int mixed_phase=2;
129  static const int high_phase=3;
130  //@}
131 
132  /// Base APR EOS
134 
135  /// Table for output
137  /// HDF file for output
139 
140  /// Function for the Maxwell construction in Fig. 7
141  int maxwell_fig7(size_t nv, const ubvector &x, ubvector &y) {
142 
143  n.n=x[0];
144  p.n=x[1];
145  n2.n=x[2];
146  p2.n=x[3];
147 
148  if (x[0]<0.0 || x[1]<0.0 || x[2]<0.0 || x[3]<0.0) return 1;
149 
150  // Low-density phase
151  ap.pion=eos_had_apr::ldp;
152  ap.calc_e(n,p,hb);
153 
154  e.mu=n.mu-p.mu;
155  mu.mu=e.mu;
156  fzt.calc_mu_zerot(e);
157  fzt.calc_mu_zerot(mu);
158  l.ed=e.ed+mu.ed;
159  l.pr=e.pr+mu.pr;
160  l.en=e.en+mu.en;
161 
162  // High-density phase
163  ap.pion=eos_had_apr::hdp;
164  ap.calc_e(n2,p2,hb2);
165 
166  e2.mu=n2.mu-p2.mu;
167  mu2.mu=e2.mu;
168  fzt.calc_mu_zerot(e2);
169  fzt.calc_mu_zerot(mu2);
170  l2.ed=e2.ed+mu2.ed;
171  l2.pr=e2.pr+mu2.pr;
172  l2.en=e2.en+mu2.en;
173 
174  // Charge neutrality for low-density phase
175  y[0]=p.n-e.n-mu.n;
176  // Equal neutron chemical potentials
177  y[1]=n.mu-n2.mu;
178  // Equal pressures
179  y[2]=hb.pr-hb2.pr;
180  // Charge neutrality for high-density phase
181  y[3]=p2.n-e2.n-mu2.n;
182 
183  return 0;
184  }
185 
186  /// Maxwell construction of the nuclear matter mixed phase
187  int mixedmaxwell(size_t nv, const ubvector &x, ubvector &y) {
188  n.n=x[0];
189  p.n=x[1];
190  n2.n=x[2];
191  p2.n=x[3];
192  chi=x[4];
193 
194  ap.pion=eos_had_apr::ldp;
195  ap.calc_e(n,p,hb);
196 
197  e.mu=n.mu-p.mu;
198  mu.mu=e.mu;
199  fzt.calc_mu_zerot(e);
200  fzt.calc_mu_zerot(mu);
201  l.ed=e.ed+mu.ed;
202  l.pr=e.pr+mu.pr;
203  l.en=e.en+mu.en;
204 
205  ap.pion=eos_had_apr::hdp;
206  ap.calc_e(n2,p2,hb2);
207 
208  e2.mu=n2.mu-p2.mu;
209  mu2.mu=e2.mu;
210  fzt.calc_mu_zerot(e2);
211  fzt.calc_mu_zerot(mu2);
212  l2.ed=e2.ed+mu2.ed;
213  l2.pr=e2.pr+mu2.pr;
214  l2.en=e2.en+mu2.en;
215 
216  y[0]=nb-chi*(n.n+p.n)-(1.0-chi)*(n2.n+p2.n);
217  y[1]=n.mu-n2.mu;
218  y[2]=hb.pr-hb2.pr;
219  y[3]=p.n-e.n-mu.n;
220  y[4]=p2.n-e2.n-mu2.n;
221 
222  return 0;
223  }
224 
225  /// Function to construct Fig. 7
226  int fig7fun(size_t nv, const ubvector &x, ubvector &y) {
227 
228  nb=x[0];
229  p.n=f7x*nb;
230  n.n=nb-p.n;
231 
232  p2.n=f7x*nb;
233  n2.n=nb-p2.n;
234 
235  ap.pion=eos_had_apr::ldp;
236  ap.calc_e(n,p,hb);
237 
238  ap.pion=eos_had_apr::hdp;
239  ap.calc_e(n2,p2,hb2);
240 
241  y[0]=hb.ed-hb2.ed;
242 
243  return 0;
244  }
245 
246  /// Solve for neutron star matter (low-density phase)
247  int nstar_low(size_t nv, const ubvector &x, ubvector &y) {
248 
249  n.n=x[0];
250  p.n=nb-n.n;
251 
252  if (n.n<0.0 || p.n<0.0) {
253  return 1;
254  }
255 
256  ap.pion=eos_had_apr::ldp;
257  ap.calc_e(n,p,hb);
258 
259  e.mu=n.mu-p.mu;
260  mu.mu=e.mu;
261  fzt.calc_mu_zerot(e);
262  fzt.calc_mu_zerot(mu);
263  l.ed=e.ed+mu.ed;
264  l.pr=e.pr+mu.pr;
265  l.en=e.en+mu.en;
266 
267  y[0]=p.n-e.n-mu.n;
268 
269  tot.pr=hb.pr+l.pr;
270  tot.ed=hb.ed+l.ed;
271 
272  return 0;
273  }
274 
275  /// Solve for neutron star matter (high-density phase)
276  int nstar_high(size_t nv, const ubvector &x, ubvector &y) {
277 
278  n2.n=x[0];
279  p2.n=nb-n2.n;
280 
281  if (n2.n<0.0 || p2.n<0.0) return 1;
282 
283  ap.pion=eos_had_apr::hdp;
284  ap.calc_e(n2,p2,hb2);
285 
286  e.mu=n2.mu-p2.mu;
287  mu.mu=e.mu;
288  fzt.calc_mu_zerot(e);
289  fzt.calc_mu_zerot(mu);
290  l.ed=e.ed+mu.ed;
291  l.pr=e.pr+mu.pr;
292  l.en=e.en+mu.en;
293 
294  y[0]=p2.n-e.n-mu.n;
295 
296  tot.pr=hb2.pr+l.pr;
297  tot.ed=hb2.ed+l.ed;
298 
299  return 0;
300  }
301 
302  /// Solve for neutron star matter (mixed phase)
303  int nstar_mixed(size_t nv, const ubvector &x, ubvector &y) {
304 
305  n.n=x[0];
306  p.n=x[1];
307  e.mu=x[2];
308  n2.n=x[3];
309  p2.n=x[4];
310  mu.mu=e.mu;
311 
312  if (phase==low_phase) chi=1.0;
313  else if (phase==high_phase) chi=0.0;
314  else chi=x[5];
315 
316  if (n.n<0.0 || n2.n<0.0) return 1;
317  if (p.n<0.0 || p2.n<0.0) return 1;
318 
319  ap.pion=eos_had_apr::ldp;
320  ap.calc_e(n,p,hb);
321 
322  ap.pion=eos_had_apr::hdp;
323  ap.calc_e(n2,p2,hb2);
324 
325  fzt.calc_mu_zerot(e);
326  fzt.calc_mu_zerot(mu);
327  l.ed=e.ed+mu.ed;
328  l.pr=e.pr+mu.pr;
329  l.en=e.en+mu.en;
330 
331  y[0]=nb-chi*(n.n+p.n)-(1.0-chi)*(n2.n+p2.n);
332  y[1]=chi*p.n+(1.0-chi)*p2.n-e.n-mu.n;
333  y[2]=n.mu-p.mu-e.mu;
334  y[3]=p2.mu-p.mu;
335  y[4]=n.mu-n2.mu;
336 
337  if (phase==mixed_phase) y[5]=hb.pr-hb2.pr;
338 
339  if (phase==low_phase) {
340  tot.pr=hb.pr+l.pr;
341  tot.ed=hb.ed+l.ed;
342  } else if (phase==mixed_phase) {
343  tot.pr=hb.pr+l.pr;
344  tot.ed=hb.ed*chi+hb2.ed*(1.0-chi)+l.ed;
345  } else {
346  tot.pr=hb2.pr+l.pr;
347  tot.ed=hb2.ed+l.ed;
348  }
349 
350  return 0;
351  }
352 
353  /// Write a line of data to the table
354  void store_data() {
355 
356  std::vector<double> line;
357  line.push_back(tot.ed);
358  line.push_back(tot.pr);
359  line.push_back(nb);
360  line.push_back((chi*n.n+(1.0-chi)*n2.n));
361  line.push_back((chi*p.n+(1.0-chi)*p2.n));
362  line.push_back(n.n);
363  line.push_back(p.n);
364  line.push_back(n2.n);
365  line.push_back(p2.n);
366  line.push_back(chi);
367  line.push_back(e.n);
368  line.push_back(mu.n);
369  line.push_back(n.mu);
370  line.push_back(p.mu);
371  line.push_back(e.mu);
372  line.push_back(mu.mu);
373  line.push_back(n.ms);
374  line.push_back(p.ms);
375  line.push_back(n2.ms);
376  line.push_back(p2.ms);
377  line.push_back(n.kf);
378  line.push_back(p.kf);
379  line.push_back(n2.kf);
380  line.push_back(p2.kf);
381  line.push_back(e.kf);
382  line.push_back(mu.kf);
383  if (line.size()!=at.get_ncolumns()) {
384  O2SCL_ERR("Table misaligned.",exc_efailed);
385  }
386 
387  at.line_of_data(line.size(),line);
388  return;
389  }
390 
391  /// Solve for nuclear matter (mixed phase)
392  int nucmixed(size_t nv, const ubvector &x, ubvector &y) {
393 
394  n.n=x[0];
395  n2.n=x[1];
396 
397  if (phase==low_phase) chi=1.0;
398  else if (phase==high_phase) chi=0.0;
399  else chi=x[2];
400 
401  if (n.n<0.0 || n2.n<0.0) return 1;
402 
403  p.n=n.n;
404  p2.n=n2.n;
405 
406  ap.pion=eos_had_apr::ldp;
407  ap.calc_e(n,p,hb);
408 
409  ap.pion=eos_had_apr::hdp;
410  ap.calc_e(n2,p2,hb2);
411 
412  y[0]=n.mu-n2.mu;
413  y[1]=nb-chi*(n.n+p.n)-(1.0-chi)*(n2.n+p2.n);
414 
415  if (phase==mixed_phase) y[2]=hb.pr-hb2.pr;
416 
417  if (phase==low_phase) {
418  tot.pr=hb.pr;
419  tot.ed=hb.ed;
420  } else if (phase==mixed_phase) {
421  tot.pr=hb.pr;
422  tot.ed=hb.ed*chi+hb2.ed*(1.0-chi);
423  } else {
424  tot.pr=hb2.pr;
425  tot.ed=hb2.ed;
426  }
427 
428 
429  return 0;
430  }
431 
432  /// Solve for neutron matter (mixed phase)
433  int neutmixed(size_t nv, const ubvector &x, ubvector &y) {
434 
435  n.n=x[0];
436  n2.n=x[1];
437 
438  if (phase==low_phase) chi=1.0;
439  else if (phase==high_phase) chi=0.0;
440  else chi=x[2];
441 
442  if (n.n<0.0 || n2.n<0.0) return 1;
443 
444  p.n=0.0;
445  p2.n=0.0;
446 
447  ap.pion=eos_had_apr::ldp;
448  ap.calc_e(n,p,hb);
449 
450  ap.pion=eos_had_apr::hdp;
451  ap.calc_e(n2,p2,hb2);
452 
453  y[0]=n.mu-n2.mu;
454  y[1]=nb-chi*(n.n+p.n)-(1.0-chi)*(n2.n+p2.n);
455 
456  if (phase==mixed_phase) y[2]=hb.pr-hb2.pr;
457 
458  if (phase==low_phase) {
459  tot.pr=hb.pr;
460  tot.ed=hb.ed;
461  } else if (phase==mixed_phase) {
462  tot.pr=hb.pr;
463  tot.ed=hb.ed*chi+hb2.ed*(1.0-chi);
464  } else {
465  tot.pr=hb2.pr;
466  tot.ed=hb2.ed;
467  }
468 
469  return 0;
470  }
471 
472  /// Solve for phase transition to nuclei
473  int nucleimat(size_t nv, const ubvector &ex, ubvector &ey) {
474 
475  double nn1,np1,mun1,mup1;
476  double nn2,np2,mun2,mup2;
477  double u=0.0;
478  thermo th1, th2;
479 
480  n.n=ex[0];
481  p.n=ex[1];
482  ap.pion=eos_had_apr::ldp;
483  ap.calc_e(n,p,th1);
484  nn1=n.n;
485  np1=p.n;
486  mun1=n.mu;
487  mup1=p.mu;
488 
489  n.n=ex[2];
490  p.n=0.0;
491  ap.pion=eos_had_apr::ldp;
492  ap.calc_e(n,p,th2);
493  nn2=n.n;
494  np2=p.n;
495  mun2=n.mu;
496  mup2=p.mu;
497 
498  e.mu=mun1-mup1;
499  fzt.calc_mu_zerot(e);
500  l.ed=e.ed;
501  l.pr=e.pr;
502  l.en=e.en;
503  if (nv>3) u=ex[3];
504  else u=1.0;
505 
506  ey[0]=(mun1-mun2)/1.0;
507  ey[1]=(th1.pr-th2.pr)/1.0;
508  ey[2]=u*np1-e.n;
509  if (nv>3)ey[3]= barn-u*(np1+nn1)-(1.0-u)*nn2;
510 
511  tot.pr=th1.pr+e.pr;
512  hb.ed=u*th1.ed+(1.0-u)*th2.ed;
513  tot.ed=u*th1.ed+(1.0-u)*th2.ed+e.ed;
514 
515  n.n=ex[0];
516  p.n=ex[1];
517  ap.pion=eos_had_apr::ldp;
518  ap.calc_e(n,p,th1);
519 
520  return 0;
521  }
522 
523  /// Solve for phase transition to nuclei with a proton drip
524  int nucleimat_pdrip(size_t nv, const ubvector &ex, ubvector &ey) {
525 
526  double nn1,np1,mun1,mup1;
527  double nn2,np2,mun2,mup2;
528  double u=0.0;
529  thermo th1, th2;
530 
531  if (ex[0]<0.0 || ex[1]<0.0 || ex[2]<0.0 || ex[3]<0.0) {
532  return 1;
533  }
534 
535  n.n=ex[0];
536  p.n=ex[1];
537  ap.pion=eos_had_apr::ldp;
538  ap.calc_e(n,p,th1);
539  nn1=n.n;
540  np1=p.n;
541  mun1=n.mu;
542  mup1=p.mu;
543 
544  n.n=ex[2];
545  p.n=ex[3];
546  ap.pion=eos_had_apr::ldp;
547  ap.calc_e(n,p,th2);
548  nn2=n.n;
549  np2=p.n;
550  mun2=n.mu;
551  mup2=p.mu;
552 
553  e.mu=mun1-mup1;
554  fzt.calc_mu_zerot(e);
555  l.ed=e.ed;
556  l.pr=e.pr;
557  l.en=e.en;
558  if (nv>4) u=ex[4];
559  else u=1.0;
560 
561  ey[0]=(mun1-mun2)/1.0;
562  ey[1]=(th1.pr-th2.pr)/1.0;
563  ey[2]=e.n-u*np1-(1.0-u)*np2;
564  ey[3]=(mup1-mup2)/1.0;
565  if (nv>4) ey[4]=barn-u*(np1+nn1)-(1.0-u)*(nn2+np2);
566 
567  tot.pr=th1.pr+e.pr;
568  hb.ed=u*th1.ed+(1.0-u)*th2.ed;
569  tot.ed=u*th1.ed+(1.0-u)*th2.ed+e.ed;
570 
571  n.n=ex[0];
572  p.n=ex[1];
573  ap.pion=eos_had_apr::ldp;
574  ap.calc_e(n,p,th1);
575 
576  return 0;
577  }
578 
579 public:
580 
581  ex_eos_had_apr() {
582 
583  n.init(939.0/hc_mev_fm,2.0);
584  p.init(939.0/hc_mev_fm,2.0);
585  n2.init(939.0/hc_mev_fm,2.0);
586  p2.init(939.0/hc_mev_fm,2.0);
587 
588  // Ensure that this works without GNU units
591  ("kg","1/fm",o2scl_mks::mass_electron),2.0);
593  ("kg","1/fm",o2scl_mks::mass_muon),2.0);
595  ("kg","1/fm",o2scl_mks::mass_electron),2.0);
597  ("kg","1/fm",o2scl_mks::mass_muon),2.0);
598 
599  n.non_interacting=false;
600  p.non_interacting=false;
601  n2.non_interacting=false;
602  p2.non_interacting=false;
603 
604  at.inc_maxlines(2000);
605  }
606 
607  /// Main driver, computing the APR EOS and the associated M vs. R curve
608  void run() {
609 
610  test_mgr t;
611  t.set_output_level(1);
612 
613  // Output directory
614  std::string prefix="ex_eos_had_apr_";
615  cout << "Set output prefix to '" << prefix << "' ." << endl;
616  cout << endl;
617 
618  // Density grid
619  double nbstart, nb_end, dnb;
620 
621  // Density at which to start looking for a mixed phase
622  double nb_switch;
623 
624  // Temporary filename
625  string filenm;
626 
627  // Error code value
628  int ret;
629 
630  // Output file stream
631  ofstream fout;
632 
633  // Variables and function values for solvers
634  ubvector x(7), y(7), xx(3), yy(3);
635 
636  // Function objects
637  mm_funct11 f_nucmixed=std::bind
638  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
639  (&ex_eos_had_apr::nucmixed),this,
640  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
641  mm_funct11 f_neutmixed=std::bind
642  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
644  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
645  mm_funct11 f_nstar_mixed=std::bind
646  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
648  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
649  mm_funct11 f_nstar_low=std::bind
650  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
652  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
653  mm_funct11 f_nstar_high=std::bind
654  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
656  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
657  mm_funct11 f_fig7fun=std::bind
658  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
659  (&ex_eos_had_apr::fig7fun),this,
660  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
661  mm_funct11 f_maxwell_fig7=std::bind
662  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
664  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
665  mm_funct11 f_mixedmaxwell=std::bind
666  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
668  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
669 
670  // Init density grid
671  nbstart=0.005;
672  dnb=0.002;
673  nb_end=2.0;
674 
675  // Select APR EOS
676  choice=1;
677  ap.select(choice);
678 
679  // Init solver tolerances
680  solver.tol_abs=1.0e-10;
681  solver.tol_rel=1.0e-12;
682  solver_trans_density.tol_abs=1.0e-10;
683  solver_trans_density.tol_rel=1.0e-12;
684 
685  // Initialize table
686  at.line_of_names(((string)("ed pr nb "))+
687  "nn np nn1 np1 nn2 np2 chi ne nmu "+
688  "mun mup mue mumu "+
689  "msn msp msn2 msp2 kfn kfp kfn2 kfp2 kfe kfmu");
690  at.set_unit("ed","1/fm^4");
691  at.set_unit("pr","1/fm^4");
692  at.set_unit("nb","1/fm^3");
693  at.set_unit("nn","1/fm^3");
694  at.set_unit("np","1/fm^3");
695  at.set_unit("nn1","1/fm^3");
696  at.set_unit("np1","1/fm^3");
697  at.set_unit("nn2","1/fm^3");
698  at.set_unit("np2","1/fm^3");
699  at.set_unit("ne","1/fm^3");
700  at.set_unit("nmu","1/fm^3");
701  at.set_unit("mun","1/fm");
702  at.set_unit("mup","1/fm");
703  at.set_unit("mue","1/fm");
704  at.set_unit("mumu","1/fm");
705  at.set_unit("msn","1/fm");
706  at.set_unit("msp","1/fm");
707  at.set_unit("msn2","1/fm");
708  at.set_unit("msp2","1/fm");
709  at.set_unit("kfn","1/fm");
710  at.set_unit("kfp","1/fm");
711  at.set_unit("kfn2","1/fm");
712  at.set_unit("kfp2","1/fm");
713  at.set_unit("kfe","1/fm");
714  at.set_unit("kfmu","1/fm");
715 
716  //--------------------------------------------------------------------
717  // Saturation properties
718 
719  ap.set_thermo(hb);
720  ap.set_mroot(solver);
721  ap.set_n_and_p(n,p);
722  ap.saturation();
723 
724  cout << "--------- Saturation properties --------------------\n" << endl;
725  cout << "n_0: " << ap.n0 << " fm^{-3}\nE_B: " << ap.eoa*hc_mev_fm
726  << " MeV\nK: " << ap.comp*hc_mev_fm << " MeV\nS: "
727  << ap.esym*hc_mev_fm << " MeV\nM^{*}/M: " << ap.msom << endl;
728  t.test_rel(ap.n0,0.1598371,1.0e-5,"n0");
729  t.test_rel(ap.eoa*hc_mev_fm,-16.00066,1.0e-5,"eoa");
730  t.test_rel(ap.esym*hc_mev_fm,32.5887,1.0e-5,"esym");
731  cout << endl;
732 
733  //--------------------------------------------------------------------
734  // Init lepton fields to zero to start
735 
736  e.mu=0.0;
737  e.n=0.0;
738  e.kf=0.0;
739  mu.mu=0.0;
740  mu.n=0.0;
741  mu.kf=0.0;
742 
743  //--------------------------------------------------------------------
744  // Nuclear matter
745 
746  cout << "--------- Nuclear matter ---------------------------\n" << endl;
747 
748  // Begin with just the low density phase
749 
750  nb_switch=0.19;
751 
752  ap.pion=eos_had_apr::ldp;
753  chi=1.0;
754  at.clear_data();
755  for(nb=nbstart;nb<nb_switch;nb+=dnb) {
756 
757  n.n=nb/2.0;
758  p.n=n.n;
759 
760  ret=ap.calc_e(n,p,hb);
761  if (ret!=0) {
762  O2SCL_ERR("Failed to compute nuclear matter.",exc_efailed);
763  }
764  tot.pr=hb.pr;
765  tot.ed=hb.ed;
766 
767  store_data();
768  }
769 
770  phase=low_phase;
771  x[0]=nb_switch/2.0;
772  x[1]=x[0];
773 
774  // Now start searching for the mixed phase
775 
776  for(nb=nb_switch;nb<=nb_end;nb+=dnb) {
777 
778  if (phase!=mixed_phase) ret=solver.msolve(2,x,f_nucmixed);
779  else ret=solver.msolve(3,x,f_nucmixed);
780  if (ret!=0) {
781  cout << nb << endl;
782  O2SCL_ERR("Solving nuclear matter failed.",
783  exc_efailed);
784  }
785 
786  if (hb.pr<hb2.pr && phase==low_phase) {
787  cout << "Mixed phase begins near nb=" << nb << " fm^{-3}." << endl;
788  phase=mixed_phase;
789  // Pick a value of chi close to, but less than, one
790  x[2]=0.90;
791  nb-=dnb;
792  } else if (phase==mixed_phase && x[2]<0.0) {
793  cout << "Mixed phase ends near nb=" << nb << " fm^{-3}." << endl;
794  phase=high_phase;
795  nb-=dnb;
796  } else {
797  store_data();
798  }
799  }
800 
801  string fn1=((string)prefix)+"nuc.o2";
802  hf.open_or_create(fn1);
803  hdf_output(hf,at,"nuc");
804  hf.close();
805  cout << "Generated file '" << fn1 << "'." << endl;
806 
807  cout << endl;
808 
809  //--------------------------------------------------------------------
810  // Neutron matter
811 
812  cout << "--------- Neutron matter ---------------------------\n" << endl;
813 
814  ap.pion=eos_had_apr::ldp;
815  chi=1.0;
816  at.clear_data();
817 
818  nb_switch=0.16;
819 
820  // Begin with just the low density phase
821 
822  for(nb=nbstart;nb<nb_switch;nb+=dnb) {
823 
824  n.n=nb;
825  p.n=0.0;
826 
827  ret=ap.calc_e(n,p,hb);
828  if (ret!=0) {
829  O2SCL_ERR("Failed to compute neutron matter.",exc_efailed);
830  }
831 
832  tot.pr=hb.pr;
833  tot.ed=hb.ed;
834 
835  store_data();
836  }
837 
838  phase=low_phase;
839  x[0]=nb_switch;
840  x[1]=nb_switch;
841 
842  // Now start searching for the mixed phase
843 
844  for(nb=nb_switch;nb<=nb_end;nb+=dnb) {
845 
846  if (phase!=mixed_phase) {
847  ret=solver.msolve(2,x,f_neutmixed);
848  } else {
849  ret=solver.msolve(3,x,f_neutmixed);
850  }
851  if (ret!=0) {
852  cout << nb << endl;
853  O2SCL_ERR("Solving neutron matter failed.",
854  exc_efailed);
855  }
856 
857  if (hb.pr<hb2.pr && phase==low_phase) {
858  cout << "Mixed phase begins near nb=" << nb << " fm^{-3}." << endl;
859  phase=mixed_phase;
860  // Pick a value of chi close to, but less than, one
861  x[2]=0.90;
862  nb-=dnb;
863  } else if (phase==mixed_phase && x[2]<0.0) {
864  cout << "Mixed phase ends near nb=" << nb << " fm^{-3}." << endl;
865  phase=high_phase;
866  nb-=dnb;
867  } else {
868  store_data();
869  }
870  }
871 
872  string fn2=((string)prefix)+"neut.o2";
873  hf.open_or_create(fn2);
874  hdf_output(hf,at,"neut");
875  hf.close();
876  cout << "Generated file '" << fn2 << "'." << endl;
877 
878  cout << endl;
879 
880  //--------------------------------------------------------------------
881  // Neutron star matter - Maxwell construction
882 
883  cout << "--------- Neutron star matter ----------------------\n" << endl;
884 
885  // Removing this for now, as it caused negative densities
886  if (false) {
887 
888  cout << "Maxwell construction." << endl;
889 
890  // Verify Figure 7
891 
892  filenm=prefix;
893  filenm+="fig7.1.txt";
894  fout.open(filenm.c_str());
895  fout.setf(ios::scientific);
896  x[0]=0.32;
897  for(f7x=0.5;f7x>=-0.001;f7x-=0.025) {
898  ret=solver.msolve(1,x,f_fig7fun);
899  if (ret!=0) {
900  O2SCL_ERR("Failed to perform maxwell.",exc_efailed);
901  }
902  fout << f7x << " " << x[0] << endl;
903  }
904  fout.close();
905  cout << "Generated file '" << filenm << "'." << endl;
906 
907  // Compute the beginning and ending densities
908  // for the Maxwell construction
909 
910  x[0]=0.94*0.2;
911  x[1]=0.06*0.2;
912  x[2]=0.95*0.24;
913  x[3]=0.05*0.24;
914  solver.msolve(4,x,f_maxwell_fig7);
915 
916  double nb_low=x[0]+x[1], nb_high=x[2]+x[3];
917  cout << "Mixed phase begins at nb=" << nb_low << " fm^{-3}." << endl;
918  cout << "Mixed phase ends at nb=" << nb_high << " fm^{-3}." << endl;
919  x[0]=nbstart/1.1;
920 
921  filenm=prefix;
922  filenm+="fig7.2.txt";
923  fout.open(filenm.c_str());
924 
925  /// Compute matter at densities below the maxwell construction
926 
927  for(nb=0.02;nb<nb_low*1.00001;nb+=(nb_low-nbstart)/20.0) {
928  ret=solver.msolve(1,x,f_nstar_low);
929  if (ret!=0) {
930  cout << nb << endl;
931  f_nstar_low(1,x,y);
932  cout << x[0] << " " << y[0] << endl;
933  O2SCL_ERR("Solving Maxwell construction failed.",
934  exc_efailed);
935  }
936  fout << p.n/nb << " " << nb << endl;
937  }
938 
939  // Compute matter at densities inside the Maxwell construction
940 
941  x[0]=(1.0-0.07)*nb_low;
942  x[1]=0.07*nb_low;
943  x[2]=(1.0-0.06)*nb_high;
944  x[3]=0.06*nb_high;
945  x[4]=1.0;
946  dnb=(nb_high-nb_low)/20.0;
947  for(nb=nb_low+dnb;nb<=nb_high*1.00001;nb+=dnb) {
948  ret=solver.msolve(5,x,f_mixedmaxwell);
949  if (ret!=0) {
950  cout << nb << endl;
951  O2SCL_ERR("Solving Maxwell construction (part 2) failed.",
952  exc_efailed);
953  }
954  fout << (chi*p.n+(1.0-chi)*p2.n)/nb << " " << nb << endl;
955  }
956 
957  // Compute matter at densities above the Maxwell construction
958 
959  x[0]=0.23;
960  for(nb=nb_high;nb<nb_end;nb+=(nb_end-nb_high)/40.0) {
961  ret=solver.msolve(1,x,f_nstar_high);
962  if (ret!=0) {
963  cout << nb << endl;
964  O2SCL_ERR("Solving Maxwell construction (part 3) failed.",
965  exc_efailed);
966  }
967  fout << p2.n/nb << " " << nb << endl;
968  }
969  fout.close();
970  cout << "Generated file '" << filenm << "'." << endl;
971  }
972 
973  //--------------------------------------------------------------------
974  // Neutron star matter - Gibbs construction
975 
976  cout << "\nGibbs construction." << endl;
977 
978  dnb=0.002;
979  ap.pion=eos_had_apr::ldp;
980  chi=1.0;
981  at.clear_data();
982 
983  x[0]=nbstart/1.1;
984  for(nb=nbstart;nb<=0.1701;nb+=dnb) {
985  ret=solver.msolve(1,x,f_nstar_low);
986  if (ret!=0) {
987  cout << nb << endl;
988  O2SCL_ERR("Solving Gibbs construction failed.",
989  exc_efailed);
990  }
991  store_data();
992  }
993  nb-=dnb;
994 
995  phase=low_phase;
996  x[0]=n.n;
997  x[1]=p.n;
998  x[2]=n.mu-p.mu;
999  x[3]=n.n;
1000  x[4]=p.n;
1001 
1002  for(nb=0.17;nb<=nb_end;nb+=dnb) {
1003 
1004  if (phase!=mixed_phase) {
1005  ret=solver.msolve(5,x,f_nstar_mixed);
1006  nstar_mixed(5,x,y);
1007  } else {
1008  ret=solver.msolve(6,x,f_nstar_mixed);
1009  }
1010  if (ret!=0) {
1011  O2SCL_ERR("Gibbs construction (part 2) failed.",exc_esanity);
1012  }
1013 
1014  if (hb.pr<hb2.pr && phase==low_phase) {
1015  cout << "Mixed phase begins at nb=" << nb << " fm^{-3}." << endl;
1016  phase=mixed_phase;
1017  x[5]=0.90;
1018  nb-=dnb;
1019  } else if (phase==mixed_phase && x[5]<0.0) {
1020  cout << "Mixed phase ends at nb=" << nb << " fm^{-3}." << endl;
1021  phase=high_phase;
1022  nb-=dnb;
1023  } else {
1024  store_data();
1025  }
1026  }
1027 
1028  //--------------------------------------------------------------------
1029  // Estimate transition density
1030 
1031  cout << "\nEstimate of transition density." << endl;
1032 
1033  ubvector newx(12), newy(12);
1034  mm_funct11 nuclei_f=std::bind
1035  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
1036  (&ex_eos_had_apr::nucleimat),this,
1037  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
1038  mm_funct11 nucleip_f=std::bind
1039  (std::mem_fn<int(size_t,const ubvector &,ubvector &)>
1041  std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
1042 
1043  solver_trans_density.tol_abs/=1.0e4;
1044  solver_trans_density.tol_rel/=1.0e4;
1045 
1046  newx[0]=0.1;
1047  newx[1]=0.001;
1048  newx[2]=0.1;
1049  ret=solver_trans_density.msolve(3,newx,nuclei_f);
1050  if (ret!=0) {
1051  nucleimat(3,newx,newy);
1052  if (fabs(newy[0])>1.0e-8 || fabs(newy[1])>1.0e-8 ||
1053  fabs(newy[2])>1.0e-8) {
1054  cout << endl;
1055  cout << "Problem in transition density (1)." << endl;
1056  cout << err_hnd->get_str() << endl;
1057  cout << newx[0] << " " << newy[0] << endl;
1058  cout << newx[1] << " " << newy[1] << endl;
1059  cout << newx[2] << " " << newy[2] << endl;
1060  cout << endl;
1061  }
1062  }
1063  cout << "Without proton drip: density: " << newx[0]+newx[1]
1064  << " fm^{-3},\n pressure: "
1065  << at.interp("nb",newx[0]+newx[1],"pr") << " fm^{-4}." << endl;
1066 
1067  solver_trans_density.tol_abs=5.0e-8;
1068  solver_trans_density.tol_rel=5.0e-8;
1069 
1070  newx[3]=0.001;
1071  ret=solver_trans_density.msolve(4,newx,nucleip_f);
1072  if (ret!=0) {
1073  nucleimat_pdrip(4,newx,newy);
1074  if (fabs(newy[0])>1.0e-8 || fabs(newy[1])>1.0e-8 ||
1075  fabs(newy[2])>1.0e-8 || fabs(newy[3])>1.0e-8) {
1076  cout << endl;
1077  cout << "Problem in transition density (2)." << endl;
1078  cout << err_hnd->get_str() << endl;
1079  cout << newx[0] << " " << newy[0] << endl;
1080  cout << newx[1] << " " << newy[1] << endl;
1081  cout << newx[2] << " " << newy[2] << endl;
1082  cout << newx[3] << " " << newy[3] << endl;
1083  cout << endl;
1084  }
1085  }
1086  nucleimat_pdrip(4,newx,newy);
1087  cout << newx[0] << " " << newy[0] << endl;
1088  cout << newx[1] << " " << newy[1] << endl;
1089  cout << newx[2] << " " << newy[2] << endl;
1090  cout << newx[3] << " " << newy[3] << endl;
1091  cout << "With proton drip: density: " << newx[0]+newx[1]
1092  << " fm^{-3},\n pressure: "
1093  << at.interp("nb",newx[0]+newx[1],"pr") << " fm^{-4}." << endl;
1094  cout << endl;
1095 
1096  solver_trans_density.tol_abs=1.0e-16;
1097  solver_trans_density.tol_rel=1.0e-16;
1098 
1099  // Output to file
1100 
1101  string fn3=((string)prefix)+"nstar.o2";
1102  hf.open_or_create(fn3);
1103  hdf_output(hf,at,"nstar");
1104  hf.close();
1105  cout << "Generated file '" << fn3 << "'." << endl;
1106 
1107  cout << endl;
1108 
1109  //--------------------------------------------------------------------
1110  // Integrate TOV equations
1111 
1112  cout << "--------- TOV solver, M vs. R. ---------------------\n" << endl;
1113 
1114  tov_solve atov;
1115  eos_tov_interp teos;
1116  atov.verbose=0;
1117  teos.verbose=0;
1118  atov.set_units("1/fm^4","1/fm^4","1/fm^3");
1119  teos.default_low_dens_eos();
1120  teos.read_table(at,"ed","pr","nb");
1121  atov.set_eos(teos);
1122  atov.calc_gpot=true;
1123 
1124  // Get results
1125 
1126  atov.mvsr();
1127  std::shared_ptr<table_units<> > tov_tmp=atov.get_results();
1128  cout << "Maximum mass is " << tov_tmp->max("gm")
1129  << " solar masses." << endl;
1130  t.test_rel(tov_tmp->max("gm"),2.191338,5.0e-4,"max mass.");
1131 
1132  // Output to file
1133 
1134  string fn4=((string)prefix)+"mvsr.o2";
1135  hf.open_or_create(fn4);
1136  hdf_output(hf,*tov_tmp,"mvsr");
1137  hf.close();
1138  cout << "Generated file '" << fn4 << "'." << endl;
1139  cout << endl;
1140 
1141  cout << "--------- TOV solver, 1.4 Msun ---------------------\n" << endl;
1142 
1143  atov.fixed(1.4);
1144  std::shared_ptr<table_units<> > tov_tmp2=atov.get_results();
1145  cout << "Aprpoximate radius of a 1.4 solar mass neutron star is "
1146  << tov_tmp2->get("r",tov_tmp2->lookup("gm",1.4)) << " km." << endl;
1147  t.test_rel(tov_tmp2->get("r",tov_tmp2->lookup("gm",1.4)),11.46630,
1148  5.0e-3,"r14");
1149 
1150  string fn5=((string)prefix)+"m14.o2";
1151  hf.open_or_create(fn5);
1152  hdf_output(hf,*tov_tmp2,"m14");
1153  hf.close();
1154  cout << "Generated file '" << fn5 << "'." << endl;
1155  cout << endl;
1156 
1157  t.report();
1158 
1159  return;
1160  }
1161 
1162 };
1163 
1164 int main(void) {
1165 
1166  cout.setf(ios::scientific);
1167 
1168  ex_eos_had_apr ac;
1169  ac.run();
1170 
1171  return 0;
1172 }
1173 // End of example
void set_eos(eos_tov &ter)
Set the EOS to use.
Definition: tov_solve.h:534
virtual int mvsr()
Calculate the mass vs. radius curve.
thermo hb
Baryon thermodynamics for low-density phase.
double chi
Volume fraction of low-density phase.
int nstar_mixed(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron star matter (mixed phase)
int fig7fun(size_t nv, const ubvector &x, ubvector &y)
Function to construct Fig. 7.
double comp
Compression modulus in .
Definition: eos_had_base.h:340
bool calc_gpot
calculate the gravitational potential and the enclosed baryon mass (default false) ...
Definition: tov_solve.h:463
virtual int fixed(double target_mass, double pmax=1.0e20)
Calculate the profile of a star with fixed mass.
void default_low_dens_eos()
Default crust EOS from Negele73 and Baym71tg.
fermion e2
Electron for high-density phase.
void run()
Main driver, computing the APR EOS and the associated M vs. R curve.
int choice
Choice of model from APR.
err_hnd_type * err_hnd
fermion p2
Proton for high-density phase.
std::shared_ptr< table_units<> > get_results()
Return the results data table.
Definition: tov_solve.h:588
thermo l2
Leptonic thermodynamics for high-density phase.
const double barn
virtual void calc_mu_zerot(fermion &f)
virtual double convert(std::string from, std::string to, double val)
exc_esanity
fermion mu
Muon for low-density phase.
fermion p
Proton for low-density phase.
int verbose
control for output (default 1)
Definition: tov_solve.h:471
lib_settings_class o2scl_settings
int pion
Choice of phase (default best)
Definition: eos_had_apr.h:164
thermo tot
Total thermodynamics.
void store_data()
Write a line of data to the table.
double muq
Charge chemical potential.
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct11
int nstar_low(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron star matter (low-density phase)
void hdf_output(hdf_file &hf, o2scl::table3d &h, std::string name)
fermion mu2
Muon for high-density phase.
void set_units(double s_efactor=1.0, double s_pfactor=1.0, double s_nbfactor=1.0)
Set output units for the table.
eos_had_apr ap
Base APR EOS.
bool test_rel(data_t result, data_t expected, data_t rel_error, std::string description)
int nucmixed(size_t nv, const ubvector &x, ubvector &y)
Solve for nuclear matter (mixed phase)
double esym
Symmetry energy in .
Definition: eos_had_base.h:343
const double hc_mev_fm
void read_table(table_units<> &eosat, std::string s_cole, std::string s_colp, std::string s_colnb="")
Specify the EOS through a table.
fermion n2
Neutron for high-density phase.
table_units at
Table for output.
fermion e
Electron for low-density phase.
bool report() const
double n0
Saturation density in .
Definition: eos_had_base.h:346
exc_efailed
fermion n
Neutron for low-density phase.
double nb
Baryon density.
double f7x
Proton fraction for Fig. 7.
virtual int calc_e(fermion &n, fermion &p, thermo &th)
Equation of state as a function of density.
mroot_hybrids solver
General solver.
void clear_data()
int nucleimat_pdrip(size_t nv, const ubvector &ex, ubvector &ey)
Solve for phase transition to nuclei with a proton drip.
void set_unit(std::string scol, std::string unit)
int nucleimat(size_t nv, const ubvector &ex, ubvector &ey)
Solve for phase transition to nuclei.
double interp(std::string sx, double x0, std::string sy)
void open_or_create(std::string fname)
virtual const char * get_str()=0
int mixedmaxwell(size_t nv, const ubvector &x, ubvector &y)
Maxwell construction of the nuclear matter mixed phase.
void line_of_names(std::string newheads)
An EOS for the TOV solver using simple linear interpolation and an optional crust EOS...
Definition: eos_tov.h:675
void select(int model_index=1)
Select model.
Additional functions to read and write EOS data to HDF5 files.
convert_units & get_convert_units()
Compute the APR EOS with a Gibbs construction and the mass versus radius curve [Example class]...
int maxwell_fig7(size_t nv, const ubvector &x, ubvector &y)
Function for the Maxwell construction in Fig. 7.
#define O2SCL_ERR(d, n)
void line_of_data(size_t nv, const vec2_t &v)
thermo hb2
Baryon thermodynamics for high-density phase.
const double mass_muon
virtual void init(double mass, double dof)
Class to solve the Tolman-Oppenheimer-Volkov equations.
Definition: tov_solve.h:265
size_t get_ncolumns() const
thermo l
Leptonic thermodynamics for low-density phase.
int nstar_high(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron star matter (high-density phase)
double mub
Baryon chemical potential.
int neutmixed(size_t nv, const ubvector &x, ubvector &y)
Solve for neutron matter (mixed phase)
mroot_hybrids solver_trans_density
Solver for transition densities (lower tolerances)
hdf_file hf
HDF file for output.
virtual int msolve(size_t nn, vec_t &xx, func_t &ufunc)
double tol_abs
EOS from Akmal, Pandharipande, and Ravenhall.
Definition: eos_had_apr.h:124
bool non_interacting
void set_output_level(int l)
virtual void set_thermo(thermo &th)
Set class thermo object.
int verbose
Control for output (default 1)
Definition: eos_tov.h:63
double tol_rel
double msom
Effective mass (neutron)
Definition: eos_had_base.h:349
const double mass_electron
fermion_zerot fzt
Compute zero-temperature thermodynamics.
virtual void saturation()
Calculates some of the EOS properties at the saturation density.
double eoa
Binding energy (without the rest mass) in .
Definition: eos_had_base.h:332
virtual void set_mroot(mroot<> &mr)
Set class mroot object for use in calculating chemical potentials from densities. ...
virtual void set_n_and_p(fermion &n, fermion &p)
Set neutron and proton.
void inc_maxlines(size_t llines)
deriv_cern cd
Derivative object.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).