c Paramètres pour etat_irwin c dans le fichier de donnees coder: NOM_ETAT='etat_irwin1', autres extensions c 1a, 2, 3, 4, 5 l'absence d'extension correspond à 'etat_irwin1' c Les détails pour personnalisation se trouvent dans c SUN_STAR_DAT/ETAT_irwin/free_eos-2.2.1/src c Extraits du README de la source téléchargeable: c The EOS1, EOS1a, EOS2, and EOS3 option suites respectively take about 80, c 30, 8, and 2 times as long to compute as EOS4 for the locus of points in a c solar model. They represent various compromises between speed and quality c (see paper). For solar vibrational frequencies or generating EOS tables to c be subsequently interpolated I recommend using EOS1. For most stellar work c involving calculations of luminosities and effective temperatures, EOS4 c provides excellent results when called directly from the stellar interiors' c code, but it might be worthwhile to try EOS1 or EOS1a for a few test cases c (especially along the LMS where differences are the maximum) to make sure c the changes are within the observational errors. Also, note that c interpolating EOS tables is a non-trivial task since the tables have to be c quite large to reduce interpolation errors to negligible proportions. C*******input flags: C ifoption = 1 (pteh style) C ifoption = 2 (fjs style) C ifoption = 3 (mdh style) C ifoption = 4 (variant of mdh style to fit SCVH) C ifoption = 5 (no explicit pteh, fjs, or mdh pressure ionization, C just_use_Planck-Larkin) C ifmodified = 0 (original formulation suitable for stellar interiors C calculation in particular style) C ifmodified = 1 (modified version suitable for stellar interiors calculation. C ifmodified = 2 (modify style for best fit of opal tables + extension) C ifmodified = many other values for a number of test cases décrits dans c free_eps.f dans SUN_STAR_DAT/ETAT_irwin/free_eos-2.2.1/src C ifion is a flag that controls the way that ionization is done. In C general, the lower ifion, the slower the code, and the more C ionization details that are calculated. C ifion = -2 sets lw=4 and ifmtrace = 0, which implies all 295 ionization C stages of the 20 elements are treated in detail. C ifion = -1 sets lw=3 and ifmtrace = 0, which implies minor metals C are treated as fully ionized while H, He, and the major C metals are treated in detail. C ifion = 0 (recommended) sets lw = 3, and ifmtrace = 0 below T = 1.d6 C and ifmtrace = 1 (both major and minor metals treated as C fully ionized) above T = 1.d6. C ifion = 1 sets lw = 3, and ifmtrace = 1, which implies all major C and minor metals are always treated as fully ionized C ifion = 2 sets lw = 2, and ifmtrace = 1, which implies all elements C are treated as fully ionized. C n.b. major metals controlled by array iftracemetal in free_eos_detailed. C currently list includes C, N, O, Ne, Mg, Si, S, Fe C kif_in = 0, ln f = match_variable and tl are independent variables. C kif_in = 1, ln p = match_variable and tl are independent variables. C kif_in = -1 signal kif = 1, and ifrad = 2 (see calculated flags below). C kif_in = 2, ln rho = match_variable and tl are independent variables. C n.b. match_variable and tl are *always* the independent variables, C and it is the calling routines responsibility to place the correct C value consistent with kif_in in match_variable for each call. C the fl, pl, and rl values in the argument list are used *strictly for C output* C n.b. if kif_in is not equal to 0, then much more computer time is C required to compute the EOS because an fl iteration must C be used to to match the match_variable. However, the initial guess C for fl is improved by a Taylor series approach in this case to C reduce these fl iterations to a minimum. C*******input real quantities which are not flags: C eps(nelements_in) is an array of neps = nelements = 20 values of relative C abundance by weight divided by the appropriate atomic weight. C We_use_the atomic weight scale where (un-ionized) C(12) C has a weight of 12.00000000.... All weights are for the C un-ionized element. The eps value for an element C should be the sum of the individual isotopic eps values for C that element. The eps array refers to the C elements in the following order: C H,He,C,N,O,Ne,Na,Mg,Al,Si,P,S,Cl,A,Ca,Ti,Cr,Mn,Fe,Ni C match_variable is matched by iterative adjustment of fl when kif > 0. C for kif = 0, match_variable is interpreted as fl C tl = ln t C*******calculated flags used in free_eos_detailed call: C kif is ordinarily equal to kif_in, but the kif_in = -1 case is C used to signal kif = 1, and ifrad = 2. C ifh2 = 0 no h2 C ifh2 = 1 vdb h2 C ifh2 = 2 st h2 C ifh2 = 3 irwin h2 C ifh2plus = 0 no h2plus C ifh2plus = 1 st h2plus C ifh2plus = 2 irwin h2plus C morder = 3, 5, or 8 (or 13, 15, or 18) means 3rd, 5th, or 8th order C fermi_dirac integral approximation following EFF fit (or modified C version of EFF fit which reduces to Cody-Thacher approximation for C low relativistic correction). C morder = -3, -5, or -8 (or -13, -15, or -18) means_use_above C approximations in non-relativistic limit. C morder = 1 uses Cody-Thacher approximation directly. C morder = 21 calculate Fermi-Dirac integrals with slow, but precise C (~1.d-9 relative errors) numerical integration. C morder = -21 is same as 21 in non-relativistic limit. C morder = 23 means original 3rd order eff result C morder = -23 is same as 23 in non-relativistic limit. C if abs(ifexchange_in) > 100 then_use_linear transform approximation C for exchange treatement. Otherwise,_use_numerical transform as C described in Paper IV. C ifexchange = mod(ifexchange_in,100): C One other wrinkle on deciding ifexchange is that large degeneracy C approximations (mod(ifexchange,10) = 2) only allowed above psi_lim. C ifexchange: C To understand ifexchange, must summarize free-energy model of fex used C in research note. In general from kapusta relation, C fex is proportional to I - J + 2^1.5 pi^2/3 beta^2 K, where C K and J and known integrals which are functions of psi and beta, C and I = K^2. (see Kovetz et al 1972, ApJ 174, 109, hereafter KLVH) C 0 < ifexchange < 10 --> KLVH treatment (i.e., drop K term). C 10 <= ifexchange < 20 --> Kapusta treatment (i.e., retain K term). C 20 <= ifexchange < 30 --> special test case,_use_I term alone. C 30 <= ifexchange < 40 --> special test case,_use_J term alone. C 40 <= ifexchange < 50 --> special test case,_use_K term alone. C N.B. if ifexchange is negative, then_use_non-relativistic limit C of corresponding positive ifexchange option. C ifexchange details: C ifexchange = 1 --> G(psi) + weak relativistic correction from KLVH C ifexchange = 2 --> I-J degenerate expression from KLVH (with corrected C sign error on a2 and high numerical precision a2 and a3, see C paper IV) C ifexchange = 11 or 12 K term added (both series from CG). C ifexchange = 21 or 22 I term alone (both series from KLVH) C ifexchange = 31 or 32 -J term alone (both series from KLVH). C ifexchange = 41 or 42 K term alone. C mod(ifexchange,10) = 4 is lowest order fit of J, K C mod(ifexchange,10) = 5 is next higher order fit of J, K C mod(ifexchange,10) = 6 is highest order fit of J, K C ifcoulomb > 9, do diffraction correction C ifcoulomb > 0 means do metal Coulomb contribution to sum0 and sum2 exactly C -10 < ifcoulomb < 0 means do metal Coulomb contribution to sum0 and C sum2 using the "metal Coulomb" approximation. C n.b. this approximation is internally replaced by PTEH C approximation (next line) if either eps(1) or eps(2) are zero. C ifcoulomb < -9 means_use_PTEH approximation for sum0 and sum2. C n.b. if this approximation is combined with the PTEH C approximation for pressure ionization, then the C routine is much faster because no ionization fraction C iterations are required. C mod(|ifcoulomb|,10) = 0 ignore Coulomb interaction. C mod(|ifcoulomb|,10) = 1_use_Debye-Huckel Coulomb approximation. C mod(|ifcoulomb|,10) = 2_use_Debye-Huckel Coulomb approximation with tau(x) correction C mod(|ifcoulomb|,10) = 3_use_PTEH Coulomb approximation with their theta_e C mod(|ifcoulomb|,10) = 4_use_PTEH Coulomb approximation with fermi-dirac theta_e C mod(|ifcoulomb|,10) = 9 same as 4 with DeWitt definition of lambda C (using sum0a = sum0 + ne*theta_e). C mod(|ifcoulomb|,10) = 5_use_DH smoothly connected to modified OCP C with DeWitt definition of lambda. C mod(|ifcoulomb|,10) = 6 same as 5 with alternative smooth connection C mod(|ifcoulomb|,10) = 7 DH (Gamma < 1) or OCP using new DeWitt lambda. C mod(|ifcoulomb|,10) = 8 same as 7 with theta_e = 0. C C ifpi contains meaning of two flags: C ifpi > 0 means_use_Planck-Larkin occupation probability, otherwise not. C remaining meaning in absolute value of ifpi C |ifpi| = 0,_use_no pressure ionization C |ifpi| = 1,_use_pteh pressure ionization C |ifpi| = 2,_use_fjs pressure ionization C |ifpi| = 3,_use_MDH pressure ionization C |ifpi| = 4,_use_Saumon-like variation of MDH pressure ionization C |ifpi| > 4 same as zero, i.e.,_use_no pressure ionization. C ifrad = 0, no radiation pressure, C input match_variable is consistent (ln P excluding radiation C pressure) for kif = 1. C ifrad = 1, radiation pressure included, C input match_variable is consistent (ln P including radiation C pressure) for kif = 1. C ifrad = 2, radiation pressure included, C input match_variable is ln(ptotal-prad) for kif = 1, but C all output quantities are calculated with radiation included. C this feature is used to reduce significance loss in regions which C are dominated by radiation pressure. C for kif = 2, the input match_variable is ln rho as per normal, but C the output pressure(1) is ln(ptotal-prad). C n.b. this latter case is only used for some tables which have C rho and T as the independent variable and all quantities including C pressure derivatives calculated with radiation pressure *except for* C the pressure itself. Also note for this latter case C (ifrad = 2, kif =2) that pressure(1) = ln (ptotal-prad) will be C inconsistent with pressure(2) and pressure(3) which will be C partials of ln ptotal wrt ln rho and ln t. C lw < 3 every element treated as fully ionized. C lw = 3 trace metals treated as fully ionized. C lw > 3 very slow option with all elements treated as partially ionized. C iftc =1 only used for thermodynamic consistency tests on entropy(2) C (in which case entropy(2) returned via rtp). C*******output quantities: C fl = EFF degeneracy parameter. For kif=0 this is merely determined C from fl = match_variable. For kif>0 this is determined C from a Taylor series approach as an initial guess, then C refined through iteration inside free_eos_detailed. C t = temperature C rho = density C rl = ln rho C p = pressure (total if ifrad = 1, prad subtracted if ifrad = 0) C pl = ln p C cf = cp (- partial ln rho(T,P)/partial T)^{1/2} C cp = specific heat at constant pressure C qp and qf are legacy parameters which were used to calculate an C approximate gravitational energy generation rate when C abundances were changing. These parameters have now been C disabled (set to 1.d300) since the exact gravitational energy C generation rate for the case of changing abundances can be C calculated using the precepts in Strittmatter, P.A., C Faulkner, J.. Robertson, J.W., and Faulkner, D.J. 1970, C "A Question of Entropy", Ap.J. 161, 369-373. C http://adsbit.harvard.edu/cgi-bin/ (join URL to next line) C nph-iarticle_query?bibcode=1970ApJ...161..369S C From eq. 10a of that paper the preferred rate form is eps_grav = C -(eps_nuclear - d L/dm - eps_neutrino) = - du/dt - P d (1/rho)/dt, C where p, rho, and energy(1) (=u) are variables already returned by C free_eos. That paper also gives an alternative form (eq. 10b) of C the rate equation which is equivalent to Cox and Giuli eq. C 17.75''', but that form is deprecated since it will be invalid for C discontinuous abundances and is much more complicated than eq. C 10a. C sf = partial entropy(fl, tl)/partial fl C st = partial entropy(fl, tl)/partial tl C grada = the adiabatic temperature gradient C rtp = - partial ln rho(T,P) wrt ln T. n.b. note negative sign so C quantity should be positive for normal EOS. C rmue = rho/mu_e, where mu_e is the mean molecular C weight per electron. C fh2 = n(H+)/n(all H) C fhe2 = n(He+)/n(all He) C fhe3 = n(He++)/n(all He) C xmu1 = 1/mu, where mu is the mean molecular C weight per particle. C xmu3 = 1/mu_e, where mu is the mean molecular C weight per electron. C eta = degeneracy parameter (Cox and Guili eta), related C to EFF ln f = fl by C wf = sqrt(1.d0 + f) = d eta/d fl C eta = fl+2.d0*(wf-log(1.d0+wf)) C gamma1-gamma3 are the gammas as defined in Cox and Guili C h2rat = 2*n(H2)/n(all H) C h2plusrat = 2*n(H2+)/n(all H) C lambda = Coulomb interaction parameter (two definitions depending on C ifcoulomb). C gamma_e = Coulomb diffraction parameter. C iteration_count (if positive) = C the total number of ionization fraction loops completed C to get the solution C iteration_count (if negative) = C negative of info status variable returned by free_eos_detailed C (only if that info is non-zero, i.e., it signals an abormal end). C degeneracy, pressure, density, energy, enthalpy, and entropy C are all 3-vectors, with the second component being the derivative C of the first component wrt match_variable (except for the case where C kif_in = -1), and the third component being the derivative of C the first component wrt tl. C definitions: C degeneracy(1) is EFF degeneracy parameter ln f defined above. C pressure(1) = ln pressure (except for kif=2, ifrad=2). C density(1) = ln density. C energy(1) = internal energy per unit mass. C enthalpy(1) = enthalpy per unit mass = energy(1) + p/rho. C entropy(1) = entropy per unit mass. C*******internal quantities: C fm = partial fl(tl, match_variable)/partial match_variable (kif>0) C ft = partial fl(tl, match_variable)/partial tl (kif>0) C*******