module stellar_coeffs implicit none integer, parameter :: dp = kind(1.0d0) contains !--------------------------------------------------------------------------- function calc_ZL(ZS) result(ZL) real(dp), intent(in) :: ZS real(dp) :: ZL ZL = log10(50.0d0 * ZS) if (ZS >= 0.03d0) ZL = log10(50.0d0 * 0.03d0) if (ZS <= 0.001d0) ZL = log10(50.0d0 * 0.001d0) end function calc_ZL !--------------------------------------------------------------------------- subroutine calc_ZAMS_coeffs_HURLEY(ZS, A1, A2, A3, A4, A5, A6, A7, & B1, B2, B3, B4, B5, B6, B7, B8, B9) real(dp), intent(in) :: ZS real(dp), intent(out) :: A1, A2, A3, A4, A5, A6, A7 real(dp), intent(out) :: B1, B2, B3, B4, B5, B6, B7, B8, B9 real(dp) :: dA1(0:4), dA2(0:4), dA3(0:4), dA4(0:4) real(dp) :: dA5(0:4), dA6(0:4), dA7(0:4) real(dp) :: dB1(0:4), dB2(0:4), dB3(0:4), dB4(0:4) real(dp) :: dB5(0:4), dB6(0:4), dB7(0:4), dB8(0:4), dB9(0:4) real(dp) :: ZL dA1 = (/ 0.39704170d0, -0.32913574d0, 0.34776688d0, 0.37470851d0, 0.09011915d0 /) dA2 = (/ 8.52762600d0,-24.41225973d0, 56.43597107d0, 37.06152575d0, 5.45624060d0 /) dA3 = (/ 0.00025546d0, -0.00123461d0, -0.00023246d0, 0.00045519d0, 0.00016176d0 /) dA4 = (/ 5.43288900d0, -8.62157806d0, 13.44202049d0, 14.51584135d0, 3.39793084d0 /) dA5 = (/ 5.56357900d0,-10.32345224d0, 19.44322980d0, 18.97361347d0, 4.16903097d0 /) dA6 = (/ 0.78866060d0, -2.90870942d0, 6.54713531d0, 4.05606657d0, 0.53287322d0 /) dA7 = (/ 0.00586685d0, -0.01704237d0, 0.03872348d0, 0.02570041d0, 0.00383376d0 /) dB1 = (/ 1.71535900d0, 0.62246212d0, -0.92557761d0, -1.16996966d0, -0.30631491d0 /) dB2 = (/ 6.59778800d0, -0.42450044d0,-12.13339427d0,-10.73509484d0, -2.51487077d0 /) dB3 = (/ 10.08855000d0, -7.11727086d0,-31.67119479d0,-24.24848322d0, -5.33608972d0 /) dB4 = (/ 1.01249500d0, 0.32699690d0, -0.00923418d0, -0.03876858d0, -0.00412750d0 /) dB5 = (/ 0.07490166d0, 0.02410413d0, 0.07233664d0, 0.03040467d0, 0.00197741d0 /) dB6 = (/ 0.01077422d0, 0.00000000d0, 0.00000000d0, 0.00000000d0, 0.00000000d0 /) dB7 = (/ 3.08223400d0, 0.94472050d0, -2.15200882d0, -2.49219496d0, -0.63848738d0 /) dB8 = (/ 17.84778000d0, -7.45345690d0,-48.96066856d0,-40.05386135d0, -9.09331816d0 /) dB9 = (/ 0.00022582d0, -0.00186899d0, 0.00388783d0, 0.00142402d0, -0.00007671d0 /) ZL = log10(50.0d0 * ZS) if (ZS >= 0.03d0) ZL = log10(50.0d0 * 0.03d0) if (ZS <= 0.001d0) ZL = log10(50.0d0 * 0.001d0) A1=dA1(0)+dA1(1)*ZL+dA1(2)*ZL**2+dA1(3)*ZL**3+dA1(4)*ZL**4 A2=dA2(0)+dA2(1)*ZL+dA2(2)*ZL**2+dA2(3)*ZL**3+dA2(4)*ZL**4 A3=dA3(0)+dA3(1)*ZL+dA3(2)*ZL**2+dA3(3)*ZL**3+dA3(4)*ZL**4 A4=dA4(0)+dA4(1)*ZL+dA4(2)*ZL**2+dA4(3)*ZL**3+dA4(4)*ZL**4 A5=dA5(0)+dA5(1)*ZL+dA5(2)*ZL**2+dA5(3)*ZL**3+dA5(4)*ZL**4 A6=dA6(0)+dA6(1)*ZL+dA6(2)*ZL**2+dA6(3)*ZL**3+dA6(4)*ZL**4 A7=dA7(0)+dA7(1)*ZL+dA7(2)*ZL**2+dA7(3)*ZL**3+dA7(4)*ZL**4 B1=dB1(0)+dB1(1)*ZL+dB1(2)*ZL**2+dB1(3)*ZL**3+dB1(4)*ZL**4 B2=dB2(0)+dB2(1)*ZL+dB2(2)*ZL**2+dB2(3)*ZL**3+dB2(4)*ZL**4 B3=dB3(0)+dB3(1)*ZL+dB3(2)*ZL**2+dB3(3)*ZL**3+dB3(4)*ZL**4 B4=dB4(0)+dB4(1)*ZL+dB4(2)*ZL**2+dB4(3)*ZL**3+dB4(4)*ZL**4 B5=dB5(0)+dB5(1)*ZL+dB5(2)*ZL**2+dB5(3)*ZL**3+dB5(4)*ZL**4 B6=dB6(0)+dB6(1)*ZL+dB6(2)*ZL**2+dB6(3)*ZL**3+dB6(4)*ZL**4 B7=dB7(0)+dB7(1)*ZL+dB7(2)*ZL**2+dB7(3)*ZL**3+dB7(4)*ZL**4 B8=dB8(0)+dB8(1)*ZL+dB8(2)*ZL**2+dB8(3)*ZL**3+dB8(4)*ZL**4 B9=dB9(0)+dB9(1)*ZL+dB9(2)*ZL**2+dB9(3)*ZL**3+dB9(4)*ZL**4 end subroutine calc_ZAMS_coeffs_HURLEY !--------------------------------------------------------------------------- subroutine calc_coeffs_HURLEY(ZS_in, a) real(dp), intent(in) :: ZS_in real(dp), intent(out) :: a(0:81) real(dp) :: ZS, Z, Z2, Z3, Z4, sigma, rho real(dp) :: a11, a12, a18, a19, a29 real(dp) :: a33, a42, a44, a49, a50, a51, a52, a53, a57, a57b, a62, a64, d1 ZS = ZS_in if (ZS >= 0.03d0) ZS = 0.03d0 if (ZS <= 0.001d0) ZS = 0.001d0 Z = log10(50.0d0 * ZS) Z2 = Z * Z Z3 = Z2 * Z Z4 = Z3 * Z sigma = log10(ZS) rho = Z + 1.0d0 a(0) = 0.0d0 a(1) = 1.593890d3+Z*2.053038d3+Z2*1.231226d3+Z3*2.327785d2 a(2) = 2.706708d3+Z*1.483131d3+Z2*5.772723d2+Z3*7.411230d1 a(3) = 1.466143d2-Z*1.048442d2-Z2*6.795374d1-Z3*1.391127d1 a(4) = 4.141960d-2+Z*4.564888d-2+Z2*2.958542d-2+Z3*5.571483d-3 a(5) = 3.426349d-1 a(6) = 1.949814d1+Z*1.758178d0+Z2*(-6.008212d0)+Z3*(-4.470533d0) a(7) = 4.903830d0 a(8) = 5.212154d-2+Z*3.166411d-2+Z2*(-2.750074d-3)+Z3*(-2.271549d-3) a(9) = 1.312179d0+Z*(-3.294936d-1)+Z2*9.231860d-2+Z3*2.610989d-2 a(10) = 8.073972d-1 a11 = 1.031538d0+Z*(-2.434480d-1)+Z2*7.732821d0+Z3*6.460705d0+Z4*1.374484d0 a12 = 1.043715d0+Z*(-1.577474d0)+Z2*(-5.168234d0)+Z3*(-5.596506d0)+Z4*(-1.299394d0) a(13) = 7.859573d2+Z*(-8.542048d0)+Z2*(-2.642511d1)+Z3*(-9.585707d0) a(14) = 3.858911d3+Z*2.459681d3+Z2*(-7.630093d1)+Z3*(-3.486057d2)+Z4*(-4.861703d1) a(11) = a11 * a(14) a(12) = a12 * a(14) a(15) = 2.888720d2+Z*2.952979d2+Z2*1.850341d2+Z3*3.797254d1 a(16) = 7.196580d0+Z*5.613746d-1+Z2*3.805871d-1+Z3*8.398728d-2 d1 = max(0.097d0-0.1072d0*(sigma+3.0d0), max(0.097d0, min(0.1461d0, 0.1461d0+0.1237d0*(sigma+2.0d0)))) a(17) = 10.0d0**d1 a18 = 2.187715d-1+Z*(-2.154437d0)+Z2*(-3.768678d0)+Z3*(-1.975518d0)+Z4*(-3.021475d-1) a19 = 1.466440d0+Z*1.839725d0+Z2*6.442199d0+Z3*4.023635d0+Z4*6.957529d-1 a(20) = 2.652091d1+Z*8.178458d1+Z2*1.156058d2+Z3*7.633811d1+Z4*1.950698d1 a(18) = a18 * a(20) a(19) = a19 * a(20) a(21) = 1.472103d0+Z*(-2.947609d0)+Z2*(-3.312828d0)+Z3*(-9.945065d-1) a(22) = 3.071048d0+Z*(-5.679941d0)+Z2*(-9.745523d0)+Z3*(-3.594543d0) a(23) = 2.617890d0+Z*1.019135d0+Z2*(-3.292551d-2)+Z3*(-7.445123d-2) a(24) = 1.075567d-2+Z*1.773287d-2+Z2*9.610479d-3+Z3*1.732469d-3 a(25) = 1.476246d0+Z*1.899331d0+Z2*1.195010d0+Z3*3.035051d-1 a(26) = 5.502535d0+Z*(-6.601663d-2)+Z2*9.968707d-2+Z3*3.599801d-2 a(27) = 9.511033d1+Z*6.819618d1+Z2*(-1.045625d1)+Z3*(-1.474939d1) a(28) = 3.113458d1+Z*1.012033d1+Z2*(-4.650511d0)+Z3*(-2.463185d0) a29 = 1.413057d0+Z*4.578814d-1+Z2*(-6.850581d-2)+Z3*(-5.588658d-2) a(30) = 3.910862d1+Z*5.196646d1+Z2*2.264970d1+Z3*2.873680d0 a(31) = 4.597479d0+Z*(-2.855179d-1)+Z2*2.709724d-1 a(32) = 6.682518d0+Z*2.827718d-1+Z2*(-7.294429d-2) a(29) = a29**a(32) a33 = min(1.4d0, 1.5135d0+0.3769d0*Z) a(33) = max(0.6355d0-0.4192d0*Z, max(1.25d0, a33)) a(34) = 1.910302d-1+Z*1.158624d-1+Z2*3.348990d-2+Z3*2.599706d-3 a(35) = 3.931056d-1+Z*7.277637d-2+Z2*(-1.366593d-1)+Z3*(-4.508946d-2) a(36) = 3.267776d-1+Z*1.204424d-1+Z2*9.988332d-2+Z3*2.455361d-2 a(37) = 5.990212d-1+Z*5.570264d-2+Z2*6.207626d-2+Z3*1.777283d-2 a(38) = 7.330122d-1+Z*5.192827d-1+Z2*2.316416d-1+Z3*8.346941d-3 a(39) = 1.172768d0+Z*(-1.209262d-1)+Z2*(-1.193023d-1)+Z3*(-2.859837d-2) a(40) = 3.982622d-1+Z*(-2.296279d-1)+Z2*(-2.262539d-1)+Z3*(-5.219837d-2) a(41) = 3.571038d0+Z*(-2.223625d-2)+Z2*(-2.611794d-2)+Z3*(-6.359648d-3) a42 = 1.9848d0+Z*1.1386d0+Z2*3.5640d-1 a(42) = min(1.25d0, max(1.1d0, a42)) a(43) = 6.300d-2+Z*4.810d-2+Z2*9.840d-3 a44 = 1.200d0+Z*2.450d0 a(44) = min(1.3d0, max(0.45d0, a44)) a(45) = 2.321400d-1+Z*1.828075d-3+Z2*(-2.232007d-2)+Z3*(-3.378734d-3) a(46) = 1.163659d-2+Z*3.427682d-3+Z2*1.421393d-3+Z3*(-3.710666d-3) a(47) = 1.048020d-2+Z*(-1.231921d-2)+Z2*(-1.686860d-2)+Z3*(-4.234354d-3) a(48) = 1.555590d0+Z*(-3.223927d-1)+Z2*(-5.197429d-1)+Z3*(-1.066441d-1) a49 = 9.7700d-2+Z*(-2.3100d-1)+Z2*(-7.5300d-2) a(49) = max(a49, 0.145d0) a50 = 2.4000d-1+Z*1.8000d-1+Z2*5.9500d-1 a(50) = min(a50, 0.306d0+0.053d0*Z) a51 = 3.3000d-1+Z*1.3200d-1+Z2*2.1800d-1 a(51) = min(a51, 0.3625d0+0.062d0*Z) a52 = 1.1064d0+Z*4.1500d-1+Z2*1.8000d-1 a(52) = min(a52, 0.9d0) if (ZS > 0.01d0) a(52) = min(a52, 1.0d0) a53 = 1.1900d0+Z*3.7700d-1+Z2*1.7600d-1 a(53) = min(a53, 1.0d0) if (ZS > 0.01d0) a(53) = min(a53, 1.1d0) a(54) = 3.855707d-1+Z*(-6.104166d-1)+Z2*5.676742d0+Z3*1.060894d1+Z4*5.284014d0 a(55) = 3.579064d-1+Z*(-6.442936d-1)+Z2*5.494644d0+Z3*1.054952d1+Z4*5.280991d0 a(56) = 9.587587d-1+Z*8.777464d-1+Z2*2.017321d-1 a57 = 1.5135d0+Z*3.7690d-1 a57b = min(1.4d0, a57) a(57) = max(0.6355d0-0.4192d0*Z, max(1.25d0, a57b)) a(58) = 4.907546d-1+Z*(-1.683928d-1)+Z2*(-3.108742d-1)+Z3*(-7.202918d-2) a(59) = 4.537070d0+Z*(-4.465455d0)+Z2*(-1.612690d0)+Z3*(-1.623246d0) a(60) = 1.796220d0+Z*2.814020d-1+Z2*1.423325d0+Z3*3.421036d-1 a(61) = 2.256216d0+Z*3.773400d-1+Z2*1.537867d0+Z3*4.396373d-1 a62 = 8.4300d-2+Z*(-4.7500d-2)+Z2*(-3.5200d-2) a(62) = max(0.065d0, a62) a(63) = 7.3600d-2+Z*7.4900d-2+Z2*4.4260d-2 if (ZS < 0.004d0) a(63) = min(0.055d0, a(63)) a64 = 1.3600d-1+Z*3.5200d-2 a(64) = max(0.091d0, min(0.121d0, a64)) a(65) = 1.564231d-3+Z*1.653042d-3+Z2*(-4.439786d-3)+Z3*(-4.951011d-3)+Z4*(-1.216530d-3) a(66) = 1.4770d0+Z*2.9600d-1 a(66) = max(a(66), min(1.6d0, -0.308d0-1.046d0*Z)) a(66) = max(0.8d0, min(0.8d0-2.0d0*Z, a(66))) a(67) = 5.210157d0+Z*(-4.143695d0)+Z2*(-2.120870d0) a(68) = 1.1160d0+Z*1.6600d-1 a(68) = max(0.9d0, min(a(68), 1.0d0)) if (a(68) > a(66)) a(64) = a(58)*a(66)**a(60)/(a(59)+a(66)**a(61)) a(68) = min(a(68), a(66)) a(69) = 1.071489d0+Z*(-1.164852d-1)+Z2*(-8.623831d-2)+Z3*(-1.582349d-2) a(70) = 7.108492d-1+Z*7.935927d-1+Z2*3.926983d-1+Z3*3.622146d-2 a(71) = 3.478514d0+Z*(-2.585474d-2)+Z2*(-1.512955d-2)+Z3*(-2.833691d-3) a(72) = 9.132108d-1+Z*(-1.653695d-1)+Z2*3.636784d-2 if (ZS > 0.01d0) a(72) = max(a(72), 0.95d0) a(73) = 3.969331d-3+Z*4.539076d-3+Z2*1.720906d-3+Z3*1.897857d-4 a(74) = 1.600d0+Z*7.640d-1+Z2*3.322d-1 a(74) = max(1.4d0, min(a(74), 1.6d0)) a(75) = 8.109d-1+Z*(-6.282d-1) a(75) = max(1.0d0, min(a(75), 1.27d0)) a(75) = max(a(75), 0.6355d0-0.4192d0*Z) a(76) = 1.192334d-2+Z*1.083057d-2+Z2*1.230969d0+Z3*1.551656d0 a(76) = max(a(76), -0.1015564d0-0.2161264d0*Z-0.05182516d0*Z2) a(77) = -1.668868d-1+Z*5.818123d-1+Z2*(-1.105027d1)+Z3*(-1.668070d1) a(77) = max(-0.3868776d0-0.5457078d0*Z-0.1463472d0*Z2, min(0.0d0, a(77))) a(78) = 7.615495d-1+Z*1.068243d-1+Z2*(-2.011333d-1)+Z3*(-9.371415d-2) a(78) = max(0.0d0, min(a(78), 7.454d0+9.046d0*Z)) a(79) = 9.409838d0+Z*1.522928d0 a(79) = min(a(79), max(2.0d0, -13.3d0-18.6d0*Z)) a(80) = -2.7110d-1+Z*(-5.7560d-1)+Z2*(-8.3800d-2) a(80) = max(0.0585542d0, a(80)) a(81) = 2.4930d0+Z*1.1475d0 a(81) = min(1.5d0, max(0.4d0, a(81))) end subroutine calc_coeffs_HURLEY !--------------------------------------------------------------------------- subroutine calc_coeffs_BYRNE(ZS_in, afe, b) real(dp), intent(in) :: ZS_in, afe real(dp), intent(out) :: b(0:68) real(dp) :: ZS, ZL, z, z2, z3, z4, AFE, AFE2 real(dp) :: b29, b30, b45, b46 ZS = ZS_in if (ZS >= 0.04d0) ZS = 0.03d0 if (ZS <= 0.001d0) ZS = 0.001d0 ZL = calc_ZL(ZS) z = ZL z2 = z * z z3 = z2 * z z4 = z3 * z AFE = afe AFE2 = AFE * AFE b(0) = 8.97539d0 ! R0_low mass stars b(1) = 5.33807d2 b(2) = -3.24394d0 b(3) = 2.60967d1 + 1.54326d0*z + 1.05658d-1*AFE*z2 b(4) = 1.00967d1 b(5) = -2.18839d1 - 7.27841d0*z + 3.26805d-1*AFE b(6) = 1.00965d1 b(7) = 5.80086d2 + 4.96437d0*z b(8) = -4.29736d0 b(9) = 3.74971d0 - 6.36795d0*z + 1.79430d0*z2 b(10) = 9.54465d0 ! R0_high mass stars b(11) = -5.6305d-1 + 4.1528d-1*z + 3.0894d-1*z2 b(12) = 1.5503d0 + 3.1469d-1*z + 3.2019d-2*AFE2*z2 - 1.0484d-3*AFE2*z3 b(13) = 4.7372d-1 - 2.7101d-2*z + 1.1697d-2*AFE2*z2 b(14) = 4.2709d-4 + 4.0553d-5*z + 1.2158d-3*AFE*z + 4.9584d-4*AFE*z3 b(15) = 1.9321d0 + 4.8018d-1*z + 4.7041d-2*AFE*z ! L0 b(16) = 1.7047d0 - 8.4516d0*z + 1.3890d1*z2 - 9.8473d0*z3 + 1.3510d1*z4 + 5.76784d-1*AFE + 1.40516d1*AFE2 b(17) = 3.6911d1 - 7.6272d1*z + 9.9889d1*z2 - 2.4105d1*z3 + 2.4279d1*z4 + 2.83841d1*AFE + 4.53143d1*AFE2 b(18) = 6.9398d-1 - 5.0715d0*z + 5.2277d0*z2 - 8.3235d0*z3 + 7.9322d-1*z4 - 1.25317d-1*AFE + 1.39085d1*AFE2 b(19) = 2.8543d1 - 2.8499d1*z + 2.7650d1*z2 + 2.7932d1*z3 + 2.8138d1*z4 + 1.92527d0*AFE + 1.23523d1*AFE2 b(20) = 2.5531d1 - 2.1217d1*z + 2.5429d1*z2 + 8.6413d0*z3 + 1.5680d1*z4 + 1.71860d1*AFE + 9.70280d1*AFE2 b(21) = 3.3885d0 - 8.9833d0*z + 1.1551d1*z2 - 4.6203d0*z3 + 3.1806d-1*z4 + 9.96020d-1*AFE + 2.67763d0*AFE2 b(22) = 2.6312d-2 - 5.4357d-2*z + 6.9621d-2*z2 - 1.6750d-2*z3 + 1.8662d-2*z4 + 1.85043d-2*AFE + 2.70854d-2*AFE2 ! AGE_TMS b(23) = 6.45547d2*((z+2.0d0)**(9.02878d-1) + (1.47738d-3)*(z+2.0d0)**(6.51418d0+AFE*(-1.25468d1)+AFE2*1.55575d-2)) b(24) = 1.53973d3*((z+2.0d0)**(6.25289d-2) + (5.24213d-2)*(z+2.0d0)**(4.22274d0+AFE*(-1.04958d0)+AFE2*(-3.02614d0))) b(25) = 6.17196d-2 - 2.04701d-1*z - 4.87759d-2*z2 - 4.12218d-2*AFE + 4.38856d-2*AFE2 b(26) = 6.87345d1 - 7.13544d1*z - 7.35245d1*z2 - 4.71213d0*z3 b(27) = 3.01514d-1 + 2.02105d-1*z + 2.62031d-1*z2 - 4.92494d-2*AFE - 2.35150d-2*AFE2 b(28) = 7.45886d-1 + 1.91214d-1*z + 3.04747d-1*z2 + 2.02850d-1*z3 ! LTMS_low mass stars, plus high mass stars at max R b29 = 1.40136d0 - 2.33766d-1*z + 2.32278d0*z2 + 1.06889d0*z3 b30 = 6.71190d-1 - 1.90987d0*z - 1.39804d0*z2 - 2.60405d0*z3 - 8.73777d-1*z4 b(31) = 6.63154d2 - 6.72684d1*z + 5.35932d-1*z2 - 1.00923d1*z3 b(32) = 8.78872d0 - 1.27360d0*z + 6.95135d-1*z2 + 4.56210d-2*z3 + 7.96724d-3*z*AFE b(33) = 1.28464d0 + 5.15597d-2*z + 1.23244d-1*z2 b(34) = 3.87440d3 - 4.01073d-4*z - 7.63010d1*z2 - 3.48610d2*z3 - 4.86180d1*z4 b(29) = b29 * b(34) b(30) = b30 * b(34) b(35) = 3.68307d2 + 4.54792d2*z + 1.11320d2*z2 + 4.13110d-1*z3 b(36) = 1.57492d-1 + 8.86317d-2*z + 1.12835d-1*z2 + 4.03034d-3*z3 + 1.63748d-2*z4 b(37) = 1.35375d2 + 1.99643d0*z - 1.02517d0*z2 + 1.10427d0*z3 - 1.24576d0*z4 b(38) = 6.48550d0 - 1.37040d0*z + 6.56024d-1*z2 ! LTMS_massive stars, at TMS b(39) = 10.0d0**1.61220d1 b(40) = -5.56835d0 - 2.83833d0*z - 2.59356d-1*z2 - 1.17485d-1*AFE b(41) = 10.0d0**8.36172d0 b(42) = -1.30895d0 - 8.99849d-1*z - 2.77480d-1*z2 - 4.34803d-2*AFE b(43) = 10.0d0**4.89566d0 b(44) = 1.95814d-1 - 1.37772d-1*z + 1.06041d-2*z2 + 2.41298d-2*AFE ! RTMS_low mass stars b45 = 1.11300d0 - 1.30650d-1*z - 5.35450d-1*z2 b46 = 4.93212d-1 - 3.74170d-1*z + 2.84740d-1*z2 b(47) = 3.51950d0 + 4.26910d-1*z - 1.76940d0*z2 b(48) = 3.33410d1 + 6.99640d1*z + 1.22530d2*z2 b(45) = b45 * b(48) b(46) = b46 * b(48) b(49) = -5.57740d0 - 8.67930d0*z - 1.03770d1*z2 ! RTMS_intermediate stars b(50) = 3.54080d0 + 7.07320d-1*z b(51) = 5.39990d0 b(52) = 1.96210d-1 + 2.26140d-5*z b(53) = 6.19300d0 b(54) = 5.75789d-6 + 2.70629d-5*z + 4.73741d-5*z2 + 3.62919d-5*z3 + 1.02411d-5*z4 + 1.64681d-9*AFE b(55) = 9.26580d0 + 2.67420d-3*z - 2.08340d-3*z2 + 1.59660d-3*z3 - 1.93750d-9*z4 b(56) = 1.47630d0 - 1.34180d-3*z b(57) = 1.30210d0 - 7.52720d-4*z b(58) = 5.00750d0 - 2.70760d-3*z - 5.00400d-8*z2 - 1.66280d-8*z3 ! RTMS_massive stars, maximum radius b(59) = 1.13330d1 + 9.38750d0*z + 7.66066d0*z2 + 1.02193d0*z3 - 1.35006d-1*AFE2 b(60) = -5.31271d0 - 6.62489d0*z - 4.74497d0*z2 - 2.97616d0*z3 - 2.05366d-1*AFE b(61) = 5.81595d0 - 8.69598d0*z - 1.07450d0*z2 - 2.06903d-1*z3 b(62) = 5.38100d-1 + 1.84685d-1*z + 4.47447d-1*z2 + 5.76257d-1*z3 ! RTMS_massive stars, end of MS b(63) = 1.54080d1 + 2.21700d0*z + 8.20890d0*z2 b(64) = -7.44410d0 - 3.90170d0*z - 3.96250d0*z2 b(65) = 2.56010d1 + 3.91800d0*z + 1.15320d0*z2 b(66) = -2.81160d-1 - 8.91678d-2*z + 7.88068d-2*z2 ! AGE_TMS, deltaT for massive stars, TMS minus Rmax b(67) = -1.05390d-1 + 2.01633d-1*z + 5.53316d-1*z2 + 1.90185d-1*z3 + 2.97669d-2*AFE b(68) = 3.89793d-3 + 8.72170d-3*z + 9.33672d-3*z2 + 7.09167d-3*z3 + 3.89436d-4*AFE end subroutine calc_coeffs_BYRNE !--------------------------------------------------------------------------- subroutine calc_R0_BYRNE(ZS_arr, afe_arr, M_arr, N, R0) integer, intent(in) :: N real(dp), intent(in) :: ZS_arr(N), afe_arr(N), M_arr(N) real(dp), intent(out) :: R0(N) integer :: i real(dp) :: b(0:68) real(dp) :: R0a, R0b do i = 1, N call calc_coeffs_BYRNE(ZS_arr(i), afe_arr(i), b) R0a = (b(0)+(b(1)*M_arr(i)**b(2))+(b(3)*M_arr(i)**b(4))+(b(5)*M_arr(i)**b(6))) & /((b(7)*M_arr(i)**b(8))+(b(9)*M_arr(i)**b(10))) R0b = (b(11))+(b(12)*M_arr(i)**b(13))+(b(14)*M_arr(i)**b(15)) if (M_arr(i) <= 3.4d0) R0(i) = R0a if (M_arr(i) >= 3.5d0) R0(i) = R0b if ((M_arr(i) > 3.4d0) .and. (M_arr(i) < 3.5d0)) then R0a = (b(0)+(b(1)*3.4d0**b(2))+(b(3)*3.4d0**b(4))+(b(5)*3.4d0**b(6))) & /((b(7)*3.4d0**b(8))+(b(9)*3.4d0**b(10))) R0b = (b(11))+(b(12)*3.5d0**b(13))+(b(14)*3.5d0**b(15)) R0(i) = ((3.5d0-M_arr(i))*R0a+(M_arr(i)-3.4d0)*R0b)*10.0d0 end if end do end subroutine calc_R0_BYRNE !--------------------------------------------------------------------------- subroutine calc_L0_BYRNE(ZS_arr, afe_arr, M_arr, N, L0) integer, intent(in) :: N real(dp), intent(in) :: ZS_arr(N), afe_arr(N), M_arr(N) real(dp), intent(out) :: L0(N) integer :: i real(dp) :: b(0:68) do i = 1, N call calc_coeffs_BYRNE(ZS_arr(i), afe_arr(i), b) L0(i) = (b(16)*M_arr(i)**5.5d0+b(17)*M_arr(i)**11.0d0) & /(b(18)+M_arr(i)**3.0d0+b(19)*M_arr(i)**5.0d0 & +b(20)*M_arr(i)**7.0d0+b(21)*M_arr(i)**8.0d0+b(22)*M_arr(i)**9.5d0) end do end subroutine calc_L0_BYRNE !--------------------------------------------------------------------------- subroutine calc_tMS_BYRNE(ZS_arr, afe_arr, M_arr, N, tMS, tWR) integer, intent(in) :: N real(dp), intent(in) :: ZS_arr(N), afe_arr(N), M_arr(N) real(dp), intent(out) :: tMS(N), tWR(N) integer :: i real(dp) :: a(0:81), b(0:68) real(dp) :: ZL, tBGB, tX, tY, tHOOK, deltaT do i = 1, N ZL = calc_ZL(ZS_arr(i)) call calc_coeffs_HURLEY(ZS_arr(i), a) call calc_coeffs_BYRNE(ZS_arr(i), afe_arr(i), b) tBGB = (b(23)+b(24)*M_arr(i)**(4.0d0+b(25))+b(26)*M_arr(i)**(5.5d0+b(27))+b(28)*M_arr(i)**7.0d0) & /(a(4)*M_arr(i)**2.0d0+a(5)*M_arr(i)**7.0d0) tX = max(0.95d0, min(0.95d0-0.03d0*(ZL+0.30103d0), 0.99d0)) tY = max(0.5d0, 1.0d0-0.01d0*max(a(6)/(M_arr(i)**a(7)), a(8)+a(9)/(M_arr(i)**a(10)))) tHOOK = tY * tBGB tMS(i) = max(tHOOK, tX*tBGB) deltaT = max(0.0d0, b(67)+M_arr(i)*b(68)) tWR(i) = tMS(i) - deltaT end do end subroutine calc_tMS_BYRNE !--------------------------------------------------------------------------- subroutine calc_RTMS_BYRNE(ZS_arr, afe_arr, M_arr, N, RTMS, RWR) integer, intent(in) :: N real(dp), intent(in) :: ZS_arr(N), afe_arr(N), M_arr(N) real(dp), intent(out) :: RTMS(N), RWR(N) integer :: i real(dp) :: a(0:81), b(0:68) real(dp) :: ZL, MWRR real(dp) :: RT1, RT2, RT3, RW4 real(dp) :: RT2i, RT3i, ip_alpha, x1 real(dp) :: R0(N) call calc_R0_BYRNE(ZS_arr, afe_arr, M_arr, N, R0) RTMS = 0.0d0 RWR = 0.0d0 do i = 1, N ZL = calc_ZL(ZS_arr(i)) call calc_coeffs_HURLEY(ZS_arr(i), a) call calc_coeffs_BYRNE(ZS_arr(i), afe_arr(i), b) MWRR = 2164.9d0*(ZL+2.0d0)**(-7.0d0) + 64.947d0*(ZL+2.0d0)**(-0.8d0) RT1 = (b(45)+b(46)*M_arr(i)**b(47))/(b(48)+M_arr(i)**b(49)) RT2 = (b(50)*M_arr(i)**(b(51))+b(52)*M_arr(i)**b(53)+b(54)*M_arr(i)**b(55)) & /(b(56)+b(57)*M_arr(i)**(b(58))) RT3 = (10.0d0**b(59))*M_arr(i)**b(60)+b(61)*M_arr(i)**b(62) RW4 = (10.0d0**b(63))*M_arr(i)**b(64)+b(65)*M_arr(i)**b(66) if (M_arr(i) < a(17)) then RTMS(i) = RT1 if (RTMS(i) <= R0(i)) RTMS(i) = R0(i) end if if ((M_arr(i) >= a(17)+0.1d0) .and. (M_arr(i) < 0.90d0*MWRR)) RTMS(i) = RT2 if (M_arr(i) >= 1.10d0*MWRR) RTMS(i) = RT3 if ((M_arr(i) >= a(17)) .and. (M_arr(i) < a(17)+0.1d0)) then RT1 = (b(45)+b(46)*a(17)**b(47))/(b(48)+a(17)**b(49)) RT2 = (b(50)*(a(17)+0.1d0)**(b(51))+b(52)*(a(17)+0.1d0)**b(53)+b(54)*(a(17)+0.1d0)**b(55)) & /(b(56)+b(57)*(a(17)+0.1d0)**(b(58))) RTMS(i) = (((a(17)+0.1d0)-M_arr(i))*RT1+(M_arr(i)-a(17))*RT2)*10.0d0 end if if ((M_arr(i) >= MWRR*0.90d0) .and. (M_arr(i) <= MWRR*1.10d0)) then RT2i = (b(50)*(MWRR*0.90d0)**(b(51))+b(52)*(MWRR*0.90d0)**b(53)+b(54)*(MWRR*0.90d0)**b(55)) & /(b(56)+b(57)*(MWRR*0.90d0)**(b(58))) RT3i = (10.0d0**b(59))*(MWRR*1.10d0)**b(60)+b(61)*(MWRR*1.10d0)**b(62) ip_alpha = log10(RT3i / RT2i) / log10(1.10d0/0.9d0) x1 = MWRR * 0.9d0 RTMS(i) = RT2i * (M_arr(i) / x1)**ip_alpha end if if (M_arr(i) >= MWRR) then RWR(i) = RW4 if (RWR(i) > RTMS(i)) RWR(i) = RTMS(i) end if if (M_arr(i) <= a(62)) RTMS(i) = max(1.5d0*R0(i), RTMS(i)) end do end subroutine calc_RTMS_BYRNE !--------------------------------------------------------------------------- subroutine calc_LTMS_BYRNE(ZS_arr, afe_arr, M_arr, N, LTMS, LWR) integer, intent(in) :: N real(dp), intent(in) :: ZS_arr(N), afe_arr(N), M_arr(N) real(dp), intent(out) :: LTMS(N), LWR(N) integer :: i real(dp) :: a(0:81), b(0:68) real(dp) :: ZL, MWRL real(dp) :: LT1, LT2 LTMS = 0.0d0 LWR = 0.0d0 do i = 1, N ZL = calc_ZL(ZS_arr(i)) call calc_coeffs_HURLEY(ZS_arr(i), a) call calc_coeffs_BYRNE(ZS_arr(i), afe_arr(i), b) MWRL = 694.252d0*(ZL+2.0d0)**(-3.0d0) LT1 = ((b(29)*M_arr(i)**3.0d0)+(b(30)*M_arr(i)**4.0d0)+(b(31)*M_arr(i)**(b(32)+b(33)))) & /((b(34))+(b(35)*M_arr(i)**5.0d0)+(b(36)*M_arr(i)**(b(32)))+(b(37)*M_arr(i)**b(38))) LT2 = ((b(39)*M_arr(i)**b(40))+(b(41)*M_arr(i)**b(42))+(b(43)*M_arr(i)**b(44))) LTMS(i) = LT1 if (M_arr(i) >= MWRL) LWR(i) = LT2 end do end subroutine calc_LTMS_BYRNE end module stellar_coeffs