nucleus_rmf.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 nucleus_rmf.h
24  \brief File defining \ref o2scl::nucleus_rmf
25 */
26 #ifndef RMF_NUCLEUS_H
27 #define RMF_NUCLEUS_H
28 
29 #include <iostream>
30 #include <string>
31 #include <vector>
32 #include <o2scl/interp.h>
33 #include <o2scl/constants.h>
34 #include <o2scl/part.h>
35 #include <o2scl/eos_had_rmf.h>
36 #include <o2scl/table_units.h>
37 #include <o2scl/ode_rkck_gsl.h>
38 #include <o2scl/ode_funct.h>
39 
40 #ifndef DOXYGEN_NO_O2NS
41 namespace o2scl {
42 #endif
43 
44  /** \brief Spherical closed-shell nuclei with a relativistic
45  mean-field model in the Hartree approximation
46 
47  This code is very experimental.
48 
49  This class is based on a code developed by C.J. Horowitz and
50  B.D. Serot, and used in \ref Horowitz81 which was then adapted
51  by P.J. Ellis and used in \ref Heide94 and \ref Prakash94. Ellis
52  and A.W. Steiner adapted it for the parameterization in in \ref
53  eos_had_rmf for \ref Steiner05b, and then converted to C++ by
54  Steiner afterwards.
55 
56  The standard usage is something like:
57  \code
58  nucleus_rmf rn;
59  o2scl_hdf::rmf_load(rn.rmf,"NL4");
60  rn.run_nucleus(82,208,0,0);
61  cout << rn.rnrp << endl;
62  \endcode
63  which computes the structure of \f$ ^{208}\mathrm{Pb} \f$ and
64  outputs the neutron skin thickness using the model \c 'NL4'.
65 
66  Potential exceptions are
67  - Failed to converge
68  - Failed to solve meson field equations
69  - Energy not finite (usually a problem in the equation of
70  state)
71  - Energy not finite in final calculation
72  - Function \ref iterate() called before \ref init_run()
73  - Not a closed-shell nucleus
74 
75  The initial level pattern is
76  \verbatim
77  1 S 1/2
78  // 2 nucleons
79  1 P 3/2
80  1 P 1/2
81  // 8 nucleus
82  1 D 5/2
83  1 D 3/2
84  2 S 1/2
85  // 20 nucleons
86  1 F 7/2
87  // 28 nucleons
88  1 F 5/2
89  2 P 3/2
90  2 P 1/2
91  // 40 nucleons
92  1 G 9/2
93  // 50 nucleus
94  1 G 7/2
95  2 D 5/2
96  1 H 11/2
97  2 D 3/2
98  3 S 1/2
99  // 82 nucleons
100  1 H 9/2
101  2 F 7/2
102  1 I 13/2
103  2 F 5/2
104  3 P 3/2
105  3 P 1/2
106  // 126 nucleons
107  2 G 9/2
108  1 I 11/2
109  1 J 15/2
110  3 D 5/2
111  4 S 1/2
112  2 G 7/2
113  3 D 3/2
114  // 184 nucleons
115  \endverbatim
116 
117  Below, \f$ \alpha \f$ is a generic index for the isospin, the
118  radial quantum number \f$ n \f$ and the angular quantum numbers
119  \f$ \kappa \f$ and \f$ m \f$. The meson fields are \f$
120  \sigma(r), \omega(r) \f$ and \f$ \rho(r) \f$. The baryon density
121  is \f$ n(r) \f$, the neutron and proton densities are \f$ n_n(r)
122  \f$ and \f$ n_p(r) \f$, and the baryon scalar density is \f$
123  n_s(r) \f$.
124  The nucleon field equations are
125  \f{eqnarray*}
126  F^{\prime}_{\alpha}(r)- \frac{\kappa}{r} F_{\alpha}(r)
127  + \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r)
128  - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2} \right)
129  e A(r) - M + g_{\sigma} \sigma(r) \right] G_{\alpha}(r) &=& 0 \\
130  G^{\prime}_{\alpha}(r)+ \frac{\kappa}{r} G_{\alpha}(r)
131  - \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r)
132  - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2}
133  \right) e A(r) + M - g_{\sigma} \sigma(r) \right] F_{\alpha}(r) &=& 0
134  \f}
135  where the isospin number, \f$ t_{\alpha} \f$ is \f$ 1/2 \f$ for
136  protons and \f$ -1/2 \f$ for neutrons.
137  The meson field equations are
138  \f{eqnarray*}
139  \sigma^{\prime \prime}(r) + \frac{2}{r} \sigma^{\prime}(r)
140  - m_{\sigma}^2 \sigma &=& - g_{\sigma} n_s(r) + b M g_{\sigma}^3
141  \sigma^2 + c g_{\sigma}^4 \sigma^3 - g_{\rho}^2 \rho^2
142  \frac{\partial f}{\partial \sigma} \\
143  \omega^{\prime \prime}(r) + \frac{2}{r} \omega^{\prime}(r)
144  - m_{\omega}^2 \omega &=& - g_{\omega} n(r) +
145  \frac{\zeta}{6} g_{\omega}^4 \omega^3 + g_{\rho}^2 \rho^2
146  \frac{\partial f}{\partial \omega} \\
147  \rho^{\prime \prime}(r) + \frac{2}{r} \rho^{\prime}(r)
148  - m_{\rho}^2 \rho &=& - \frac{g_{\rho}}{2}
149  \left[n_n(r)-n_p(r)\right] + 2 g_{\rho}^2 \rho f + \frac{\xi}{6}
150  g_{\rho}^4 \rho^3
151  \f}
152  and the Coulomb field equation is
153  \f[
154  A^{\prime \prime}(r) + \frac{2}{r} A^{\prime}(r) = - e n_p(r)
155  \f]
156  The meson field equations plus a pair of Dirac-like nucleon
157  field equations for each index \f$ \alpha \f$ must be solved to
158  compute the structure of a given nucleus.
159 
160  The densities (scalar, baryon, isospin, and charge) are
161  \f{eqnarray*}
162  n_s(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2-F(r)^2
163  \right] \right\} \\
164  n(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2+F(r)^2
165  \right] \right\} \\
166  n_i(r) &=& \sum_{\alpha} \left\{ t_{\alpha} \int d^3 r
167  \left[ G(r)^2-F(r)^2 \right] \right\} \\
168  n_c(r) &=& \sum_{\alpha} \left\{ \left[t_{\alpha}+\frac{1}{2}\right]
169  \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\}
170  \f}
171 
172  The total energy is
173  \f[
174  E = \sum_{\alpha} \varepsilon_{\alpha} (2 j_{\alpha}+1)
175  - \frac{1}{2} \int d^{3} r
176  \left[- g_{\sigma} \sigma(r) \rho_s(r) + g_{\omega} \omega(r)
177  \rho(r) + \frac{1}{2} g_{\rho} \rho(r) + e A(r) n_p(r) \right]
178  \f]
179 
180  The charge density is the proton density folded with the
181  charge density of the proton, i.e.
182  \f[
183  \rho_{\mathrm{ch}}(r) = \int d^{3} r^{\prime}
184  \rho_{\mathrm{prot}}(r-r^{\prime}) \rho_p(r)
185  \f]
186  where the proton charge density is assumed to be of the form
187  \f[
188  \rho_{\mathrm{prot}}(r) = \frac{\mu^3}{8 \pi} \exp \left(
189  - \mu |r|\right)
190  \f]
191  and the parameter \f$ \mu = (0.71)^{1/2}~\mathrm{GeV} \f$ (see
192  Eq. 20b in \ref Horowitz81). The default value of \ref a_proton
193  is the value of \f$ \mu \f$ converted into \f$ \mathrm{fm}^{-1}
194  \f$.
195 
196  Generally, the first array index associated with a function
197  is not the value at \f$ r=0 \f$, but at \f$ r=\Delta r \f$
198  (stored in \ref step_size) and so the \f$ r=0 \f$ part of the
199  algorithm is handled separately.
200 
201  \todo Better documentation
202  \todo Convert energies() to use EOS and possibly
203  replace sigma_rhs() and related functions by the associated
204  field equation method of eos_had_rmf.
205 
206  \todo I believe currently the center of mass correction
207  for the binding energy per nucleon is not done and has
208  to be added in after the fact
209 
210  \future Sort energy levels at the end by energy
211  \future Improve the numerical methods
212  \future Make the neutron and proton orbitals more configurable
213  \future Generalize to \f$ m_n \neq m_p \f$ .
214  \future Allow more freedom in the integrations
215  \future Consider converting everything to inverse fermis.
216  \future Convert to zero-indexed arrays (mostly done)
217  \future Warn when the level ordering is wrong, and unoccupied levels
218  are lower energy than occupied levels
219  \future Connect with \ref o2scl::nucmass ?
220  */
221  class nucleus_rmf {
222 
225 
226 #ifndef DOXYGEN_INTERNAL
227 
228  protected:
229 
230  /// A convenient struct for the solution of the Dirac equations
231  typedef struct {
232  /// Eigenvalue
233  double eigen;
234  /// Quantum number \f$ \kappa \f$ .
235  double kappa;
236  /// The meson fields
238  /// Desc
240  } odparms;
241 
242  /// The total number of shells stored internally
243  static const int n_internal_levels=29;
244 
245 #endif
246 
247  public:
248 
249  nucleus_rmf();
250 
251  ~nucleus_rmf();
252 
253  /** \name Basic operation
254  */
255  //@{
256  /// A shell of nucleons for \ref nucleus_rmf
257  typedef struct {
258  /// Degeneracy \f$ 2 j+1 \f$ .
259  int twojp1;
260  /// \f$ \kappa \f$
261  int kappa;
262  /// Energy eigenvalue
263  double energy;
264  /// Isospin ( \f$ +1/2 \f$ or \f$ -1/2 \f$ .
265  double isospin;
266  /// Angular momentum-spin state \f$ ^{2s+1} \ell_{j} \f$
267  std::string state;
268  /// Matching radius (in fm)
269  double match_point;
270  /// Desc.
271  double eigen;
272  /// Desc.
273  double eigenc;
274  /// Number of nodes in the wave function
275  int nodes;
276  } shell;
277 
278  /** \brief Computes the structure of a nucleus with the specified
279  number of levels
280 
281  Note that \ref rmf must be set before run_nucleus() is called.
282 
283  This calls init_run(), and then iterate() until \c iconverged is
284  1, and then post_converge().
285  */
286  int run_nucleus(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
287 
288  /// Set output level
289  void set_verbose(int v) { verbose=v; };
290  //@}
291 
292  /** \name Lower-level interface
293  */
294  //@{
295  /** \brief Initialize a run
296 
297  Note that \ref rmf must be set before run_nucleus() is called.
298  */
299  void init_run(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
300 
301  /// Perform an iteration
302  int iterate(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N,
303  int &iconverged, int &dirac_converged, int &meson_converged);
304 
305  /// After convergence, make CM corrections, etc.
306  int post_converge(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
307 
308  //@}
309 
310  /// \name Results
311  //@{
312  /** \brief Get the radial profiles
313 
314  The profiles are calculated each iteration by iterate().
315  */
316  std::shared_ptr<table_units<> > get_profiles() { return profiles; };
317 
318  /** \brief The final charge densities
319  */
320  std::shared_ptr<table_units<> > get_chden() { return chden_table; };
321 
322  /// The number of levels
323  int nlevels;
324 
325  /** \brief The levels (protons first, then neutrons)
326 
327  An array of size \ref nlevels
328  */
329  std::vector<shell> levels;
330 
331  /// The number of unoccupied levels (equal to \c unocc_Z + \c unocc_N)
333 
334  /** \brief The unoccupied levels (protons first, then neutrons)
335 
336  An array of size \ref nuolevels
337  */
338  std::vector<shell> unocc_levels;
339 
340  /** \brief Surface tension (in \f$ \mathrm{fm}^{-3} \f$ )
341 
342  Computed in post_converge() or automatically in run_nucleus()
343  */
344  double stens;
345 
346  /** \brief Skin thickness (in fm)
347 
348  Computed every iteration in iterate()
349  or automatically in run_nucleus()
350  */
351  double rnrp;
352 
353  /** \brief Neutron RMS radius (in fm)
354 
355  Computed every iteration in iterate() or automatically in
356  run_nucleus()
357  */
358  double rnrms;
359 
360  /** \brief Proton RMS radius (in fm)
361 
362  Computed every iteration in iterate() or automatically in
363  run_nucleus()
364  */
365  double rprms;
366 
367  /** \brief Total energy (in MeV)
368 
369  Computed every iteration in iterate() or automatically in
370  run_nucleus()
371  */
372  double etot;
373 
374  /** \brief Charge radius (in fm)
375 
376  Computed in post_converge() or automatically in run_nucleus()
377  */
378  double r_charge;
379 
380  /** \brief Charge radius corrected by the center of mass (in fm)
381 
382  Computed in post_converge() or automatically in run_nucleus()
383  */
384  double r_charge_cm;
385  //@}
386 
387  /** \name Equation of state
388  */
389  //@{
390  /** \brief The default equation of state (default NL3)
391 
392  This is set in the constructor to be the default
393  model, NL3, using the function \ref load_nl3().
394  */
396 
397  /** \brief Set the base EOS to be used
398 
399  The equation of state must be set before run_nucleus() or
400  init_fun() are called, including the value of eos_had_rmf::mnuc.
401  */
403  rmf=&r;
404  return 0;
405  }
406 
407  /** \brief \ref thermo object for the EOS
408 
409  This is just used as temporary storage.
410  */
412 
413  /** \brief The neutron
414 
415  The mass of the neutron is ignored and set by init_run()
416  to be eos_had_rmf::mnuc from \ref rmf.
417  */
418  fermion n;
419 
420  /** \brief The proton
421 
422  The mass of the proton is ignored and set by init_run()
423  to be eos_had_rmf::mnuc from \ref rmf.
424  */
425  fermion p;
426  //@}
427 
428  /** \name Numeric configuration
429  */
430  //@{
431  /** \brief If true, call the error handler if the routine does not
432  converge or reach the desired tolerance (default true)
433  */
435 
436  /// Set the stepper for the Dirac differential equation
438  ode_funct> &step) {
439  ostep=&step;
440  return;
441  }
442 
443  /// Maximum number of total iterations (default 70)
444  int itmax;
445 
446  /** \brief Maximum number of iterations for solving the meson
447  field equations (default 10000)
448  */
450 
451  /** \brief Maximum number of iterations for solving the Dirac
452  equations (default 100)
453  */
455 
456  /// Tolerance for Dirac equations (default \f$ 5 \times 10^{-3} \f$ ).
457  double dirac_tol;
458 
459  /** \brief Second tolerance for Dirac equations
460  (default \f$ 5 \times 10^{-4} \f$ ).
461  */
462  double dirac_tol2;
463 
464  /// Tolerance for meson field equations (default \f$ 10^{-6} \f$ ).
465  double meson_tol;
466  //@}
467 
468  /** \brief Initial guess structure
469 
470  The initial guess for the meson field profiles is
471  a set of Fermi-Dirac functions, i.e.
472  \f[
473  \sigma(r)=\mathrm{sigma0}/
474  [1+\exp((r-\mathrm{fermi\_radius})/\mathrm{fermi\_width})]
475  \f]
476  */
477  typedef struct {
478  /// Scalar field at r=0
479  double sigma0;
480  /// Vector field at r=0
481  double omega0;
482  /// Isubvector field at r=0
483  double rho0;
484  /// Coulomb field at r=0
485  double A0;
486  /// The radius for which the fields are half their central value
487  double fermi_radius;
488  /// The "width" of the Fermi-Dirac function
489  double fermi_width;
490  } initial_guess;
491 
492  /** \brief Parameters for initial guess
493 
494  Default is {310,240,-6,25.9,6.85,0.6}
495  */
497 
498  /** \brief If true, use the generic ODE solver instead of the
499  internal 4th order Runge-Kutta
500  */
502 
503 #ifndef DOXYGEN_INTERNAL
504 
505  protected:
506 
507  /** \brief The parameter for the charge density of the proton
508  (default is about 4.27073)
509  */
510  double a_proton;
511 
512  /// Load the default model NL3 into the given \ref eos_had_rmf object
513  int load_nl3(eos_had_rmf &r);
514 
515  /// The base EOS
517 
518  /** \brief The radial profiles
519  */
520  std::shared_ptr<table_units<> > profiles;
521 
522  /** \brief The final charge densities
523  */
524  std::shared_ptr<table_units<> > chden_table;
525 
526  /** \brief A pointer to the current vector of levels
527  (either \ref levels or \ref unocc_levels)
528  */
529  std::vector<shell> *levp;
530 
531  /** \brief Control output (default 1)
532  */
533  int verbose;
534 
535  /// The starting neutron levels
537 
538  /// The starting proton levels
540 
541  /// The grid size
542  static const int grid_size=300;
543 
544  /// The grid step size (default 0.04)
545  double step_size;
546 
547  /// The nucleon mass (automatically set in init_fun())
548  double mnuc;
549 
550  /** \name The meson and photon fields and field equations (protected)
551  */
552  //@{
553  /// Values of the fields from the last iteration
555 
556  /// The values of the fields
558 
559  /// The Green's functions inside
561 
562  /// The Green's functions outside
564 
565  /// Scalar density RHS
566  double sigma_rhs(double sig, double ome, double rho);
567 
568  /// Vector density RHS
569  double omega_rhs(double sig, double ome, double rho);
570 
571  /// Isubvector density RHS
572  double rho_rhs(double sig, double ome, double rho);
573 
574  /** \brief Calculate meson and photon
575  Green's functions \ref gin and \ref gout
576  */
577  void meson_init();
578 
579  /** \brief Calculate meson and photon fields
580 
581  The meson and photon field equations are of the
582  Sturm-Liouville form, e.g.
583  \f[
584  \left[r^2 \sigma^{\prime}(r) \right]^{\prime} -
585  r^2 m_{\sigma}^2 \sigma(r) = f(r)
586  \f]
587  where \f$ \sigma(0) = \sigma_0 \f$ and \f$ \sigma(+\infty) = 0
588  \f$. A solution of the homogeneous equation with \f$ f(r) =0
589  \f$ is \f$ \sigma(r) = \sigma_0 \sinh( m_{\sigma} r)/
590  (m_{\sigma} r) \f$. The associated Green's function is
591  \f[
592  D(r,r^{\prime},m_{\sigma}) = \frac{-1}{m_{\sigma} r r^{\prime}}
593  \sinh (m_{\sigma} r_{<}) \exp (-m_{\sigma} r_{>})
594  \f]
595  where \f$ r_{<} \f$ is the smaller of \f$ r \f$ and
596  \f$ r^{\prime} \f$ and \f$ r_{>} \f$ is the larger.
597  One can then write the solution of the meson field
598  equation given the density
599  \f[
600  \sigma(r) = \int_0^{\infty} r^{\prime 2}~d r^{\prime}
601  \left[ - g_{\sigma} n_{s}(r) \right]
602  D\left(r,r^{\prime},m_{\sigma}\right)
603  \f]
604 
605  Since radii are positive, \f$ \sinh (m r) \approx
606  e^{m r}/2 \f$ and
607  \f[
608  D(r,r^{\prime},m_{\sigma}) \approx \left[
609  \frac{-1}{2 m_{\sigma} r_{<}}
610  \exp (m_{\sigma} r_{<})
611  \right] \left[
612  \frac{1}{r_{>}}
613  \exp (-m_{\sigma} r_{>})
614  \right]
615  \f]
616  The two terms in the Green's function are separated into
617  \f[
618  \mathrm{gin}(r) = \frac{e^{m_{\sigma} r}}{2 m_{\sigma} r}
619  \f]
620  and
621  \f[
622  \mathrm{gout}(r) = \frac{e^{-m_{\sigma} r}}{r} \, .
623  \f]
624  These functions are computed in \ref meson_init() . Then the
625  field is given by
626  \f[
627  \sigma(r)= \left(\frac{e^{-m_{\sigma} r}}{r}\right)
628  \int_0^{r} r^{\prime 2} g_{\sigma} n_{s}
629  \left(\frac{e^{m_{\sigma} r^{\prime}}}{2 m_{\sigma}
630  r^{\prime}} \right)~d r^{\prime} +
631  \left(\frac{e^{m_{\sigma} r}}{2 m_{\sigma} r} \right)
632  \int_r^{\infty} r^{\prime 2} g_{\sigma} n_{s}
633  \left(\frac{e^{-m_{\sigma} r^{\prime}}}
634  {r^{\prime}}\right)~d r^{\prime}
635  \f]
636  where the first integral is stored in <tt>xi2</tt> and
637  the second is in <tt>xi1</tt> in the function \ref meson_iter() .
638  The part of <tt>xi2</tt> at \f$ r^{\prime}=0 \f$ is
639  stored in <tt>xi20</tt>.
640  */
641  void meson_iter(double ic);
642 
643  /** \brief Solve for the meson profiles
644  */
645  int meson_solve();
646 
647  /** \brief The grid index corresponding to the nuclear surface
648  (computed by init_run())
649  */
651  //@}
652 
653  /** \name Density information (protected)
654  */
655  //@{
656  /// The densities times radius squared
658 
659  /// The proton scalar density times radius squared
661 
662  /// The scalar field RHS
664 
665  /// The vector field RHS
667 
668  /// The isubvector field RHS
670 
671  /// Charge density
673 
674  /// Charge density
676 
677  /// Baryon density
679  //@}
680 
681  /// Energy integrand
683 
684  /// Initialize the meson and photon fields, the densities, etc.
685  void init_meson_density();
686 
687  /// Calculate the energy profile
688  int energy_radii(double xpro, double xnu, double e);
689 
690  /// Compute the center of mass correction
691  void center_mass_corr(double atot);
692 
693  /** \name Calculating the form factor, etc. (protected)
694  */
695  //@{
696  /// Fold in proton form factor
697  void pfold(double x, double &xrhof);
698 
699  /// Function representing proton form factor
700  double xpform(double x, double xp, double a);
701 
702  /// Perform integrations for form factor
703  void gauss(double xmin, double xmax, double x, double &xi);
704 
705  /// Desc.
706  double xrhop(double x1);
707 
708  /// Interpolation object
710  //@}
711 
712  /** \name Solving the Dirac equations (protected)
713  */
714  //@{
715  /** \brief Solve the Dirac equations
716 
717  Solves the Dirac equation in from 12 fm to the match point and
718  then out from .04 fm and adjusts eigenvalue with
719  \f[
720  \Delta \varepsilon = -g(r=\mathrm{match\_point})
721  \times (f^{+}-f^{-})
722  \f]
723  */
724  int dirac(int ilevel);
725 
726  /// Take a step in the Dirac equations
727  void dirac_step(double &x, double h, double eigen,
728  double kappa, ubvector &varr);
729 
730  /// The form of the Dirac equations for the ODE stepper
731  int odefun(double x, size_t nv, const ubvector &y,
732  ubvector &dydx, odparms &op);
733 
734  /// Compute the fields for the Dirac equations
735  void field(double x, double &s, double &v, ubvector &varr);
736 
737  /// The default stepper
740 
741  /// The ODE stepper
744 
745  /** \brief Integrate the Dirac equations using a simple
746  inline 4th order Runge-Kutta
747  */
748  double dirac_rk4(double x, double g1, double f1, double &funt,
749  double eigen, double kappa, ubvector &varr);
750  //@}
751 
752  /// True if init() has been called
754 
755  /// ODE functions
757 
758  /// ODE derivatives
760 
761  /// ODE errors
763 
764  /// \name Gauss-Legendre integration points and weights
765  //@{
766  double x12[6], w12[6];
767  double x100[50], w100[50];
768  //@}
769 
770 #endif
771 
772  };
773 
774 #ifndef DOXYGEN_NO_O2NS
775 }
776 #endif
777 
778 #endif
boost::numeric::ublas::matrix
o2scl::nucleus_rmf::set_step
void set_step(ode_step< ubvector, ubvector, ubvector, ode_funct > &step)
Set the stepper for the Dirac differential equation.
Definition: nucleus_rmf.h:437
o2scl::nucleus_rmf::stens
double stens
Surface tension (in )
Definition: nucleus_rmf.h:344
o2scl::nucleus_rmf::pfold
void pfold(double x, double &xrhof)
Fold in proton form factor.
o2scl::nucleus_rmf::get_chden
std::shared_ptr< table_units<> > get_chden()
The final charge densities.
Definition: nucleus_rmf.h:320
o2scl::nucleus_rmf::shell::eigenc
double eigenc
Desc.
Definition: nucleus_rmf.h:273
o2scl::nucleus_rmf::initial_guess::A0
double A0
Coulomb field at r=0.
Definition: nucleus_rmf.h:485
boost::numeric::ublas::vector< double >
o2scl::nucleus_rmf::rnrms
double rnrms
Neutron RMS radius (in fm)
Definition: nucleus_rmf.h:358
o2scl::nucleus_rmf::energy
ubvector energy
Energy integrand.
Definition: nucleus_rmf.h:682
o2scl::nucleus_rmf::hb
thermo hb
thermo object for the EOS
Definition: nucleus_rmf.h:411
o2scl::nucleus_rmf::field
void field(double x, double &s, double &v, ubvector &varr)
Compute the fields for the Dirac equations.
o2scl::nucleus_rmf::def_step
ode_rkck_gsl< ubvector, ubvector, ubvector, ode_funct > def_step
The default stepper.
Definition: nucleus_rmf.h:739
o2scl::nucleus_rmf::meson_iter
void meson_iter(double ic)
Calculate meson and photon fields.
o2scl::ode_step
o2scl::nucleus_rmf::field0
ubmatrix field0
Values of the fields from the last iteration.
Definition: nucleus_rmf.h:554
o2scl::nucleus_rmf::shell::match_point
double match_point
Matching radius (in fm)
Definition: nucleus_rmf.h:269
o2scl::eos_had_rmf
Relativistic mean field theory EOS.
Definition: eos_had_rmf.h:298
o2scl::nucleus_rmf::err_nonconv
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: nucleus_rmf.h:434
o2scl::nucleus_rmf::chden1
ubvector chden1
Charge density.
Definition: nucleus_rmf.h:672
o2scl::nucleus_rmf::dirac
int dirac(int ilevel)
Solve the Dirac equations.
o2scl::ode_rkck_gsl
o2scl::nucleus_rmf::initial_guess::sigma0
double sigma0
Scalar field at r=0.
Definition: nucleus_rmf.h:479
o2scl::nucleus_rmf::run_nucleus
int run_nucleus(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
Computes the structure of a nucleus with the specified number of levels.
o2scl::nucleus_rmf::rho_rhs
double rho_rhs(double sig, double ome, double rho)
Isubvector density RHS.
o2scl::nucleus_rmf::odparms
A convenient struct for the solution of the Dirac equations.
Definition: nucleus_rmf.h:231
o2scl::nucleus_rmf::shell::state
std::string state
Angular momentum-spin state .
Definition: nucleus_rmf.h:267
o2scl::nucleus_rmf::chdenc
ubvector chdenc
Charge density.
Definition: nucleus_rmf.h:675
o2scl::nucleus_rmf::shell::nodes
int nodes
Number of nodes in the wave function.
Definition: nucleus_rmf.h:275
o2scl::nucleus_rmf::rnrp
double rnrp
Skin thickness (in fm)
Definition: nucleus_rmf.h:351
o2scl::nucleus_rmf::shell::energy
double energy
Energy eigenvalue.
Definition: nucleus_rmf.h:263
o2scl::nucleus_rmf::meson_solve
int meson_solve()
Solve for the meson profiles.
thermo_tl< double >
o2scl::nucleus_rmf::chden_table
std::shared_ptr< table_units<> > chden_table
The final charge densities.
Definition: nucleus_rmf.h:524
o2scl::nucleus_rmf::verbose
int verbose
Control output (default 1)
Definition: nucleus_rmf.h:533
o2scl::nucleus_rmf::initial_guess::fermi_radius
double fermi_radius
The radius for which the fields are half their central value.
Definition: nucleus_rmf.h:487
o2scl::nucleus_rmf::a_proton
double a_proton
The parameter for the charge density of the proton (default is about 4.27073)
Definition: nucleus_rmf.h:510
o2scl::nucleus_rmf::odparms::varr
ubvector * varr
Desc.
Definition: nucleus_rmf.h:239
o2scl::nucleus_rmf::odparms::fields
ubmatrix * fields
The meson fields.
Definition: nucleus_rmf.h:237
o2scl::nucleus_rmf::center_mass_corr
void center_mass_corr(double atot)
Compute the center of mass correction.
o2scl::nucleus_rmf::energy_radii
int energy_radii(double xpro, double xnu, double e)
Calculate the energy profile.
o2scl::nucleus_rmf::unocc_levels
std::vector< shell > unocc_levels
The unoccupied levels (protons first, then neutrons)
Definition: nucleus_rmf.h:338
o2scl::nucleus_rmf::levels
std::vector< shell > levels
The levels (protons first, then neutrons)
Definition: nucleus_rmf.h:329
o2scl::interp_vec
o2scl::nucleus_rmf::dirac_tol2
double dirac_tol2
Second tolerance for Dirac equations (default ).
Definition: nucleus_rmf.h:462
o2scl::nucleus_rmf::xrho
ubmatrix xrho
The densities times radius squared.
Definition: nucleus_rmf.h:657
o2scl::nucleus_rmf
Spherical closed-shell nuclei with a relativistic mean-field model in the Hartree approximation.
Definition: nucleus_rmf.h:221
o2scl::nucleus_rmf::fields
ubmatrix fields
The values of the fields.
Definition: nucleus_rmf.h:557
o2scl::nucleus_rmf::grid_size
static const int grid_size
The grid size.
Definition: nucleus_rmf.h:542
o2scl::nucleus_rmf::gauss
void gauss(double xmin, double xmax, double x, double &xi)
Perform integrations for form factor.
o2scl::nucleus_rmf::shell::eigen
double eigen
Desc.
Definition: nucleus_rmf.h:271
o2scl::nucleus_rmf::mnuc
double mnuc
The nucleon mass (automatically set in init_fun())
Definition: nucleus_rmf.h:548
o2scl::nucleus_rmf::rmf
eos_had_rmf * rmf
The base EOS.
Definition: nucleus_rmf.h:516
o2scl::nucleus_rmf::ode_y
ubvector ode_y
ODE functions.
Definition: nucleus_rmf.h:756
o2scl::nucleus_rmf::n_internal_levels
static const int n_internal_levels
The total number of shells stored internally.
Definition: nucleus_rmf.h:243
o2scl::nucleus_rmf::generic_ode
bool generic_ode
If true, use the generic ODE solver instead of the internal 4th order Runge-Kutta.
Definition: nucleus_rmf.h:501
o2scl::nucleus_rmf::dirac_step
void dirac_step(double &x, double h, double eigen, double kappa, ubvector &varr)
Take a step in the Dirac equations.
o2scl::nucleus_rmf::odparms::kappa
double kappa
Quantum number .
Definition: nucleus_rmf.h:235
o2scl::nucleus_rmf::init_meson_density
void init_meson_density()
Initialize the meson and photon fields, the densities, etc.
o2scl::ode_funct
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct
o2scl::nucleus_rmf::sigma_rhs
double sigma_rhs(double sig, double ome, double rho)
Scalar density RHS.
o2scl::nucleus_rmf::set_verbose
void set_verbose(int v)
Set output level.
Definition: nucleus_rmf.h:289
o2scl::nucleus_rmf::ig
initial_guess ig
Parameters for initial guess.
Definition: nucleus_rmf.h:496
o2scl::nucleus_rmf::init_run
void init_run(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
Initialize a run.
o2scl::nucleus_rmf::set_eos
int set_eos(eos_had_rmf &r)
Set the base EOS to be used.
Definition: nucleus_rmf.h:402
o2scl::nucleus_rmf::ode_yerr
ubvector ode_yerr
ODE errors.
Definition: nucleus_rmf.h:762
o2scl::nucleus_rmf::shell::twojp1
int twojp1
Degeneracy .
Definition: nucleus_rmf.h:259
o2scl::nucleus_rmf::profiles
std::shared_ptr< table_units<> > profiles
The radial profiles.
Definition: nucleus_rmf.h:520
o2scl::nucleus_rmf::xrhosp
ubvector xrhosp
The proton scalar density times radius squared.
Definition: nucleus_rmf.h:660
o2scl::nucleus_rmf::shell
A shell of nucleons for nucleus_rmf.
Definition: nucleus_rmf.h:257
o2scl::nucleus_rmf::step_size
double step_size
The grid step size (default 0.04)
Definition: nucleus_rmf.h:545
o2scl::nucleus_rmf::meson_init
void meson_init()
Calculate meson and photon Green's functions gin and gout.
o2scl::nucleus_rmf::gout
ubmatrix gout
The Green's functions outside.
Definition: nucleus_rmf.h:563
o2scl::nucleus_rmf::proton_shells
shell proton_shells[n_internal_levels]
The starting proton levels.
Definition: nucleus_rmf.h:539
o2scl::nucleus_rmf::gi
interp_vec< ubvector > * gi
Interpolation object.
Definition: nucleus_rmf.h:709
o2scl::nucleus_rmf::rprms
double rprms
Proton RMS radius (in fm)
Definition: nucleus_rmf.h:365
o2scl::nucleus_rmf::dirac_tol
double dirac_tol
Tolerance for Dirac equations (default ).
Definition: nucleus_rmf.h:457
o2scl::nucleus_rmf::omega_rhs
double omega_rhs(double sig, double ome, double rho)
Vector density RHS.
o2scl::nucleus_rmf::init_called
bool init_called
True if init() has been called.
Definition: nucleus_rmf.h:753
o2scl::nucleus_rmf::levp
std::vector< shell > * levp
A pointer to the current vector of levels (either levels or unocc_levels)
Definition: nucleus_rmf.h:529
o2scl::nucleus_rmf::def_rmf
eos_had_rmf def_rmf
The default equation of state (default NL3)
Definition: nucleus_rmf.h:395
o2scl::nucleus_rmf::r_charge
double r_charge
Charge radius (in fm)
Definition: nucleus_rmf.h:378
o2scl::nucleus_rmf::neutron_shells
shell neutron_shells[n_internal_levels]
The starting neutron levels.
Definition: nucleus_rmf.h:536
o2scl::nucleus_rmf::p
fermion p
The proton.
Definition: nucleus_rmf.h:425
o2scl::nucleus_rmf::ode_dydx
ubvector ode_dydx
ODE derivatives.
Definition: nucleus_rmf.h:759
o2scl::nucleus_rmf::shell::kappa
int kappa
Definition: nucleus_rmf.h:261
o2scl::nucleus_rmf::xrhop
double xrhop(double x1)
Desc.
o2scl::nucleus_rmf::initial_guess::rho0
double rho0
Isubvector field at r=0.
Definition: nucleus_rmf.h:483
o2scl::nucleus_rmf::n
fermion n
The neutron.
Definition: nucleus_rmf.h:418
o2scl::nucleus_rmf::shell::isospin
double isospin
Isospin ( or .
Definition: nucleus_rmf.h:265
o2scl::nucleus_rmf::ostep
ode_step< ubvector, ubvector, ubvector, ode_funct > * ostep
The ODE stepper.
Definition: nucleus_rmf.h:743
o2scl::nucleus_rmf::odefun
int odefun(double x, size_t nv, const ubvector &y, ubvector &dydx, odparms &op)
The form of the Dirac equations for the ODE stepper.
o2scl::nucleus_rmf::initial_guess::fermi_width
double fermi_width
The "width" of the Fermi-Dirac function.
Definition: nucleus_rmf.h:489
o2scl::nucleus_rmf::meson_tol
double meson_tol
Tolerance for meson field equations (default ).
Definition: nucleus_rmf.h:465
o2scl::nucleus_rmf::meson_itmax
int meson_itmax
Maximum number of iterations for solving the meson field equations (default 10000)
Definition: nucleus_rmf.h:449
o2scl::nucleus_rmf::r_charge_cm
double r_charge_cm
Charge radius corrected by the center of mass (in fm)
Definition: nucleus_rmf.h:384
o2scl::nucleus_rmf::dirac_itmax
int dirac_itmax
Maximum number of iterations for solving the Dirac equations (default 100)
Definition: nucleus_rmf.h:454
o2scl::nucleus_rmf::arho
ubvector arho
Baryon density.
Definition: nucleus_rmf.h:678
o2scl::nucleus_rmf::xpform
double xpform(double x, double xp, double a)
Function representing proton form factor.
o2scl::nucleus_rmf::itmax
int itmax
Maximum number of total iterations (default 70)
Definition: nucleus_rmf.h:444
o2scl::nucleus_rmf::dirac_rk4
double dirac_rk4(double x, double g1, double f1, double &funt, double eigen, double kappa, ubvector &varr)
Integrate the Dirac equations using a simple inline 4th order Runge-Kutta.
o2scl::nucleus_rmf::surf_index
int surf_index
The grid index corresponding to the nuclear surface (computed by init_run())
Definition: nucleus_rmf.h:650
o2scl::nucleus_rmf::xrhor
ubvector xrhor
The isubvector field RHS.
Definition: nucleus_rmf.h:669
o2scl::nucleus_rmf::get_profiles
std::shared_ptr< table_units<> > get_profiles()
Get the radial profiles.
Definition: nucleus_rmf.h:316
o2scl::nucleus_rmf::nlevels
int nlevels
The number of levels.
Definition: nucleus_rmf.h:320
o2scl::nucleus_rmf::initial_guess::omega0
double omega0
Vector field at r=0.
Definition: nucleus_rmf.h:481
o2scl::nucleus_rmf::xrhov
ubvector xrhov
The vector field RHS.
Definition: nucleus_rmf.h:666
o2scl::nucleus_rmf::xrhos
ubvector xrhos
The scalar field RHS.
Definition: nucleus_rmf.h:663
o2scl::nucleus_rmf::etot
double etot
Total energy (in MeV)
Definition: nucleus_rmf.h:372
o2scl::nucleus_rmf::load_nl3
int load_nl3(eos_had_rmf &r)
Load the default model NL3 into the given eos_had_rmf object.
o2scl::nucleus_rmf::post_converge
int post_converge(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
After convergence, make CM corrections, etc.
o2scl::nucleus_rmf::odparms::eigen
double eigen
Eigenvalue.
Definition: nucleus_rmf.h:233
o2scl::nucleus_rmf::nuolevels
int nuolevels
The number of unoccupied levels (equal to unocc_Z + unocc_N)
Definition: nucleus_rmf.h:332
o2scl::nucleus_rmf::gin
ubmatrix gin
The Green's functions inside.
Definition: nucleus_rmf.h:560
o2scl::nucleus_rmf::initial_guess
Initial guess structure.
Definition: nucleus_rmf.h:477
o2scl::nucleus_rmf::iterate
int iterate(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N, int &iconverged, int &dirac_converged, int &meson_converged)
Perform an iteration.

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