import numpy as np

def calc_ZL(ZS):

    ZL=np.log10(50e0*ZS)
    if(ZS >= 0.03):
        ZL=np.log10(50e0*0.03)
    if(ZS <=0.001):
        ZL=np.log10(50e0*0.001)
    return ZL

def calc_ZAMS_coeffs_HURLEY(ZS):

    dA1=(0.39704170,-0.32913574 ,0.34776688 ,0.37470851 ,0.09011915)
    dA2=(8.52762600,-24.41225973,56.43597107,37.06152575,5.45624060)
    dA3=(0.00025546,-0.00123461 ,-0.00023246,0.00045519 ,0.00016176)
    dA4=(5.43288900,-8.62157806 ,13.44202049,14.51584135,3.39793084)
    dA5=(5.56357900,-10.32345224,19.44322980,18.97361347,4.16903097)
    dA6=(0.78866060,-2.90870942 ,6.54713531 ,4.05606657 ,0.53287322)
    dA7=(0.00586685,-0.01704237 ,0.03872348 ,0.02570041 ,0.00383376)

    dB1=(1.71535900,0.62246212 ,-0.92557761 ,-1.16996966 ,-0.30631491)
    dB2=(6.59778800,-0.42450044,-12.13339427,-10.73509484,-2.51487077)
    dB3=(10.0885500,-7.11727086,-31.67119479,-24.24848322,-5.33608972)
    dB4=(1.01249500,0.32699690 ,-0.00923418 ,-0.03876858 ,-0.00412750)
    dB5=(0.07490166,0.02410413 ,0.07233664  ,0.03040467  ,0.00197741)
    dB6=(0.01077422,0.00000000 ,0.00000000  ,0.00000000  ,0.00000000)
    dB7=(3.08223400,0.94472050 ,-2.15200882 ,-2.49219496 ,-0.63848738)
    dB8=(17.8477800,-7.45345690,-48.96066856,-40.05386135,-9.09331816)
    dB9=(0.00022582,-0.00186899,0.00388783  ,0.00142402  ,-0.00007671)

    ZL=np.log10(50e0*ZS)
    #print(ZS,ZL)
    if(ZS >= 0.03):
        ZL=np.log10(50e0*0.03)
    if(ZS <= 0.001):
        ZL=np.log10(50e0*0.001)
    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
    
    return A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7, B8, B9

def calc_coeffs_HURLEY(ZS):

    if(ZS >= 0.03):
        ZS=0.03
    if(ZS <=0.001):
        ZS=0.001

    Z=np.log10(50e0*ZS)
    Z2=Z*Z
    Z3=Z2*Z
    Z4=Z3*Z
    sigma=np.log10(ZS)
    rho=Z+1e0

    a=np.zeros(82)

    a[0]=0
    a[1]=1.593890e+3+Z*2.053038e+3+Z2*1.231226e+3+Z3*2.327785e+2
    a[2]=2.706708e+3+Z*1.483131e+3+Z2*5.772723e+2+Z3*7.411230e+1
    a[3]=1.466143e+2-Z*1.048442e+2-Z2*6.795374e+1-Z3*1.391127e+1
    a[4]=4.141960e-2+Z*4.564888e-2+Z2*2.958542e-2+Z3*5.571483e-3
    a[5]=3.426349e-1
    a[6]=1.949814e+1 +Z*1.758178e+0+Z2*-6.008212e+0+Z3*-4.470533e+0
    a[7]=4.903830e+0
    a[8]=5.212154e-2 +Z*3.166411e-2+Z2*-2.750074e-3+Z3*-2.271549e-3
    a[9]=1.312179e+0 +Z*-3.294936e-1+Z2*9.231860e-2+Z3*2.610989e-2
    a[10]=8.073972e-1
    a11=1.031538e+0+Z*-2.434480e-1+Z2*7.732821e+0+Z3*6.460705e+0+Z4*1.374484e+0
    a12=1.043715e+0+Z*-1.577474e+0+Z2*-5.168234e+0+Z3*-5.596506e+0+Z4*-1.299394e+0
    a[13]=7.859573e+2+Z*-8.542048e+0+Z2*-2.642511e+1+Z3*-9.585707e+0
    a[14]=3.858911e+3+Z*2.459681e+3+Z2*-7.630093e+1+Z3*-3.486057e+2+Z4*-4.861703e+1
    a[11]=a11*a[14]
    a[12]=a12*a[14]
    a[15]=2.888720e+2+Z*2.952979e+2+Z2*1.850341e+2+Z3*3.797254e+1
    a[16]=7.196580e+0+Z*5.613746e-1+Z2*3.805871e-1+Z3*8.398728e-2
    d1=np.max([0.097-0.1072*(sigma+3.0),np.max([0.097,np.min([0.1461,0.1461+0.1237*(sigma+2.0)])])])
    a[17]=10e0**d1
    a18=2.187715e-1+Z*-2.154437e+0+Z2*-3.768678e+0+Z3*-1.975518e+0+Z4*-3.021475e-1
    a19=1.466440e+0+Z*1.839725e+0+Z2*6.442199e+0+Z3*4.023635e+0+Z4*6.957529e-1
    a[20]=2.652091e+1+Z*8.178458e+1+Z2*1.156058e+2+Z3*7.633811e+1+Z4*1.950698e+1
    a[18]=a18*a[20]
    a[19]=a19*a[20]
    a[21]=1.472103e+0+Z*-2.947609e+0+Z2*-3.312828e+0+Z3*-9.945065e-1
    a[22]=3.071048e+0+Z*-5.679941e+0+Z2*-9.745523e+0+Z3*-3.594543e+0
    a[23]=2.617890e+0+Z*1.019135e+0+Z2*-3.292551e-2+Z3*-7.445123e-2
    a[24]=1.075567e-2+Z*1.773287e-2+Z2*9.610479e-3+Z3*1.732469e-3
    a[25]=1.476246e+0+Z*1.899331e+0+Z2*1.195010e+0+Z3*3.035051e-1
    a[26]=5.502535e+0+Z*-6.601663e-2+Z2*9.968707e-2+Z3*3.599801e-2
    a[27]=9.511033e+1+Z*6.819618e+1+Z2*-1.045625e+1+Z3*-1.474939e+1
    a[28]=3.113458e+1+Z*1.012033e+1+Z2*-4.650511e+0+Z3*-2.463185e+0
    a29=1.413057e+0+Z*4.578814e-1+Z2*-6.850581e-2+Z3*-5.588658e-2
    a[30]=3.910862e+1+Z*5.196646e+1+Z2*2.264970e+1+Z3*2.873680e+0
    a[31]=4.597479e+0+Z*-2.855179e-1+Z2*2.709724e-1
    a[32]=6.682518e+0+Z*2.827718e-1+Z2*-7.294429e-2
    a[29]=a29**a[32]
    a33=np.min([1.4,1.5135+0.3769*Z])
    a[33]=np.max([0.6355-0.4192*Z,np.max([1.25,a33])])
    a[34]=1.910302e-1+Z*1.158624e-1+Z2*3.348990e-2+Z3*2.599706e-3
    a[35]=3.931056e-1+Z*7.277637e-2+Z2*-1.366593e-1+Z3*-4.508946e-2
    a[36]=3.267776e-1+Z*1.204424e-1+Z2*9.988332e-2+Z3*2.455361e-2
    a[37]=5.990212e-1+Z*5.570264e-2+Z2*6.207626e-2+Z3*1.777283e-2
    a[38]=7.330122e-1+Z*5.192827e-1+Z2*2.316416e-1+Z3*8.346941e-3
    a[39]=1.172768e+0+Z*-1.209262e-1+Z2*-1.193023e-1+Z3*-2.859837e-2
    a[40]=3.982622e-1+Z*-2.296279e-1+Z2*-2.262539e-1+Z3*-5.219837e-2
    a[41]=3.571038e+0+Z*-2.223625e-2+Z2*-2.611794e-2+Z3*-6.359648e-3
    a42=1.9848e+0  +Z*1.1386e+0+Z2*3.5640e-1
    a[42]=np.min([1.25,np.max([1.1,a42])])
    a[43]=6.300e-2   +Z*4.810e-2+Z2*9.840e-3
    a44=1.200e+0   +Z*2.450e+0
    a[44]=np.min([1.3,np.max([0.45,a44])])
    a[45]=2.321400e-1+Z*1.828075e-3+Z2*-2.232007e-2+Z3*-3.378734e-3
    a[46]=1.163659e-2+Z*3.427682e-3+Z2*1.421393e-3+Z3*-3.710666e-3
    a[47]=1.048020e-2+Z*-1.231921e-2+Z2*-1.686860e-2+Z3*-4.234354e-3
    a[48]=1.555590e+0+Z*-3.223927e-1+Z2*-5.197429e-1+Z3*-1.066441e-1
    a49=9.7700e-2  +Z*-2.3100e-1+Z2*-7.5300e-2
    a[49]=np.max([a49,0.145])
    a50=2.4000e-1  +Z*1.8000e-1+Z2*5.9500e-1
    a[50]=np.min([a50,0.306+0.053*Z])
    a51=3.3000e-1  +Z*1.3200e-1+Z2*2.1800e-1
    a[51]=np.min([a51,0.3625+0.062*Z])
    a52=1.1064e+0  +Z*4.1500e-1+Z2*1.8000e-1
    a[52]=np.min([a52,0.9])
    if(ZS>0.01): a[52]=np.min([a52,1.0])
    a53=1.1900e+0  +Z*3.7700e-1+Z2*1.7600e-1
    a[53]=np.min([a53,1.0])
    if(ZS>0.01): a[53]=np.min([a53,1.1])
    a[54]=3.855707e-1+Z*-6.104166e-1+Z2*5.676742e+0+Z3*1.060894e+1+Z4*5.284014e+0
    a[55]=3.579064e-1+Z*-6.442936e-1+Z2*5.494644e+0+Z3*1.054952e+1+Z4*5.280991e+0
    a[56]=9.587587e-1+Z*8.777464e-1+Z2*2.017321e-1
    a57=1.5135e+0  +Z*3.7690e-1
    a57b=np.min([1.4,a57])
    a[57]=np.max([0.6355-0.4192*Z,np.max([1.25,a57b])])
    a[58]=4.907546e-1+Z*-1.683928e-1+Z2*-3.108742e-1+Z3*-7.202918e-2
    a[59]=4.537070e+0+Z*-4.465455e+0+Z2*-1.612690e+0+Z3*-1.623246e+0
    a[60]=1.796220e+0+Z*2.814020e-1+Z2*1.423325e+0+Z3*3.421036e-1
    a[61]=2.256216e+0+Z*3.773400e-1+Z2*1.537867e+0+Z3*4.396373e-1
    a62=8.4300e-2  +Z*-4.7500e-2+Z2*-3.5200e-2
    a[62]=np.max([0.065,a62])
    a[63]=7.3600e-2  +Z*7.4900e-2+Z2*4.4260e-2
    if(ZS<0.004): a[63]=np.min([0.055,a[63]])
    a64=1.3600e-1  +Z*3.5200e-2
    a[64]=np.max([0.091,np.min([0.121,a64])])
    a[65]=1.564231e-3+Z*1.653042e-3+Z2*-4.439786e-3+Z3*-4.951011e-3+Z4*-1.216530e-03
    a[66]=1.4770e+0  +Z*2.9600e-1
    a[66]=np.max([a[66],np.min([1.6,-0.308-1.046*Z])])
    a[66]=np.max([0.8,np.min([0.8-2.0*Z,a[66]])])
    a[67]=5.210157e+0+Z*-4.143695e+0+Z2*-2.120870e+0
    a[68]=1.1160e+0  +Z*1.6600e-1
    a[68]=np.max([0.9,np.min([a[68],1.0])])
    if(a[68]>a[66]): a[64]=a[58]*a[66]**a[60]/(a[59]+a[66]**a[61])
    a[68]=np.min([a[68],a[66]])
    a[69]=1.071489e+0+Z*-1.164852e-1+Z2*-8.623831e-2+Z3*-1.582349e-2
    a[70]=7.108492e-1+Z*7.935927e-1+Z2*3.926983e-1+Z3*3.622146e-2
    a[71]=3.478514e+0+Z*-2.585474e-2+Z2*-1.512955e-2+Z3*-2.833691e-3
    a[72]=9.132108e-1+Z*-1.653695e-1+Z2*3.636784e-2
    if(ZS>0.01): a[72]=np.max([a[72],0.95])
    a[73]=3.969331e-3+Z*4.539076e-3+Z2*1.720906e-3+Z3*1.897857e-4
    a[74]=1.600e+0   +Z*7.640e-1+Z2*3.322e-1
    a[74]=np.max([1.4,np.min([a[74],1.6])])
    a[75]=8.109e-1   +Z*-6.282e-1
    a[75]=np.max([1.0,np.min([a[75],1.27])])
    a[75]=np.max([a[75],0.6355-0.4192*Z])
    a[76]=1.192334e-2+Z*1.083057e-2+Z2*1.230969e+0+Z3*1.551656e+0
    a[76]=np.max([a[76],-0.1015564-0.2161264*Z-0.05182516*Z2])
    a[77]=-1.668868e-1+Z*5.818123e-1+Z2*-1.105027e+1+Z3*-1.668070e+1
    a[77]=np.max([-0.3868776-0.5457078*Z-0.1463472*Z2,np.min([0.0,a[77]])])
    a[78]=7.615495e-1+Z*1.068243e-1+Z2*-2.011333e-1+Z3*-9.371415e-2
    a[78]=np.max([0.0,np.min([a[78],7.454+9.046*Z])])
    a[79]=9.409838e+0+Z*1.522928e+0
    a[79]=np.min([a[79],np.max([2.0,-13.3-18.6*Z])])
    a[80]=-2.7110e-1 +Z*-5.7560e-1+Z2*-8.3800e-2
    a[80]=np.max([0.0585542,a[80]])
    a[81]=2.4930e+0  +Z*1.1475e+0
    a[81]=np.min([1.5,np.max([0.4,a[81]])])
    
    return a

def calc_coeffs_BYRNE(ZS,afe):

    if(ZS >= 0.04):
        ZS=0.03
    if(ZS <=0.001):
        ZS=0.001

    ZL = calc_ZL(ZS)
    
    z = ZL
    z2 = z*z
    z3 = z2*z
    z4 = z3*z

    AFE = afe
    AFE2 = AFE*AFE

    b = np.zeros(69)
    b[0] =  8.97539e00
    #R0_low mass stars
    b[1]  = 5.33807e02
    b[2]  =-3.24394e00
    b[3]  = 2.60967e01 + 1.54326e00*z + 1.05658e-1*AFE*z2
    b[4]  = 1.00967e01
    b[5]  =-2.18839e01 - 7.27841e00*z + 3.26805e-1*AFE
    b[6]  = 1.00965e01
    b[7]  = 5.80086e02 + 4.96437e00*z
    b[8]  =-4.29736e00
    b[9]  = 3.74971e00 - 6.36795e00*z + 1.79430e00*z2
    b[10] = 9.54465e00
    #R0_high mass stars
    b[11] =-5.6305e-1 + 4.1528e-1*z + 3.0894e-1*z2 
    b[12] = 1.5503e00 + 3.1469e-1*z + 3.2019e-2*AFE2*z2 - 1.0484e-3*AFE2*z3
    b[13] = 4.7372e-1 - 2.7101e-2*z + 1.1697e-2*AFE2*z2
    b[14] = 4.2709e-4 + 4.0553e-5*z + 1.2158e-3*AFE*z   + 4.9584e-4*AFE*z3
    b[15] = 1.9321e00 + 4.8018e-1*z + 4.7041e-2*AFE*z
    #L0
    b[16] = 1.7047e00 - 8.4516e00*z + 1.3890e01*z2 - 9.8473e00*z3 + 1.3510e01*z4 + 5.76784e-1*AFE + 1.40516e01*AFE2
    b[17] = 3.6911e01 - 7.6272e01*z + 9.9889e01*z2 - 2.4105e01*z3 + 2.4279e01*z4 + 2.83841e01*AFE + 4.53143e01*AFE2
    b[18] = 6.9398e-1 - 5.0715e00*z + 5.2277e00*z2 - 8.3235e00*z3 + 7.9322e-1*z4 - 1.25317e-1*AFE + 1.39085e01*AFE2
    b[19] = 2.8543e01 - 2.8499e01*z + 2.7650e01*z2 + 2.7932e01*z3 + 2.8138e01*z4 + 1.92527e00*AFE + 1.23523e01*AFE2
    b[20] = 2.5531e01 - 2.1217e01*z + 2.5429e01*z2 + 8.6413e00*z3 + 1.5680e01*z4 + 1.71860e01*AFE + 9.70280e01*AFE2
    b[21] = 3.3885e00 - 8.9833e00*z + 1.1551e01*z2 - 4.6203e00*z3 + 3.1806e-1*z4 + 9.96020e-1*AFE + 2.67763e00*AFE2
    b[22] = 2.6312e-2 - 5.4357e-2*z + 6.9621e-2*z2 - 1.6750e-2*z3 + 1.8662e-2*z4 + 1.85043e-2*AFE + 2.70854e-2*AFE2
    #AGE_TMS
    b[23] = 6.45547e2*((z+2.)**(9.02878e-1)+(1.47738e-3)*(z+2.)**(6.51418e0+AFE*-1.25468e1+AFE2*1.55575e-2))
    b[24] = 1.53973e3*((z+2.)**(6.25289e-2)+(5.24213e-2)*(z+2.)**(4.22274e0+AFE*-1.04958e0+AFE2*-3.02614e0))
    b[25] = 6.17196e-2 - 2.04701e-1*z - 4.87759e-2*z2 - 4.12218e-2*AFE + 4.38856e-2*AFE2
    b[26] = 6.87345e01 - 7.13544e01*z - 7.35245e01*z2 - 4.71213e00*z3
    b[27] = 3.01514e-1 + 2.02105e-1*z + 2.62031e-1*z2 - 4.92494e-2*AFE - 2.35150e-2*AFE2
    b[28] = 7.45886e-1 + 1.91214e-1*z + 3.04747e-1*z2 + 2.02850e-1*z3
    #LTMS_low mass stars, plus high mass stars at max R
    b29   = 1.40136e00 - 2.33766e-1*z + 2.32278e00*z2 + 1.06889e00*z3
    b30   = 6.71190e-1 - 1.90987e00*z - 1.39804e00*z2 - 2.60405e00*z3 - 8.73777e-1*z4
    b[31] = 6.63154e02 - 6.72684e01*z + 5.35932e-1*z2 - 1.00923e01*z3
    b[32] = 8.78872e00 - 1.27360e00*z + 6.95135e-1*z2 + 4.56210e-2*z3 + 7.96724e-3*z*AFE
    b[33] = 1.28464e00 + 5.15597e-2*z + 1.23244e-1*z2
    b[34] = 3.87440e03 - 4.01073e-4*z - 7.63010e01*z2 - 3.48610e02*z3 - 4.86180e01*z4
    b[29] = b29*b[34]
    b[30] = b30*b[34]
    b[35] = 3.68307e02 + 4.54792e02*z + 1.11320e02*z2 + 4.13110e-1*z3
    b[36] = 1.57492e-1 + 8.86317e-2*z + 1.12835e-1*z2 + 4.03034e-3*z3 + 1.63748e-2*z4
    b[37] = 1.35375e02 + 1.99643e00*z - 1.02517e00*z2 + 1.10427e00*z3 - 1.24576e00*z4
    b[38] = 6.48550e00 - 1.37040e00*z + 6.56024e-1*z2
    #LTMS_massive stars, at TMS
    b[39] = 10**1.61220e01
    b[40] =-5.56835e00 - 2.83833e00*z - 2.59356e-1*z2 - 1.17485e-1*AFE
    b[41] = 10**8.36172e00
    b[42] =-1.30895e00 - 8.99849e-1*z - 2.77480e-1*z2 - 4.34803e-2*AFE
    b[43] = 10**4.8956600
    b[44] = 1.95814e-1 - 1.37772e-1*z + 1.06041e-2*z2 + 2.41298e-2*AFE
    #RTMS_low mass stars
    b45   = 1.11300e00 - 1.30650e-1*z - 5.35450e-1*z2
    b46   = 4.93212e-1 - 3.74170e-1*z + 2.84740e-1*z2
    b[47] = 3.51950e00 + 4.26910e-1*z - 1.76940e00*z2
    b[48] = 3.33410e01 + 6.99640e01*z + 1.22530e02*z2
    b[45] = b45*b[48]
    b[46] = b46*b[48] 
    b[49] =-5.57740e00 - 8.67930e00*z - 1.03770e01*z2
    #RTMS_intermediate stars
    b[50] = 3.54080e00 + 7.07320e-1*z
    b[51] = 5.39990e00 
    b[52] = 1.96210e-1 + 2.26140e-5*z
    b[53] = 6.19300e00 
    b[54] = 5.75789e-6 + 2.70629e-5*z + 4.73741e-5*z2 + 3.62919e-5*z3 + 1.02411e-5*z4 + 1.64681e-9*AFE
    b[55] = 9.26580e00 + 2.67420e-3*z - 2.08340e-3*z2 + 1.59660e-3*z3 - 1.93750e-9*z4
    b[56] = 1.47630e00 - 1.34180e-3*z 
    b[57] = 1.30210e00 - 7.52720e-4*z 
    b[58] = 5.00750e00 - 2.70760e-3*z - 5.00400e-8*z2 - 1.66280e-8*z3   
    #RTMS_massive stars, maximum radius
    b[59] = 1.13330e01 + 9.38750e00*z + 7.66066e00*z2 + 1.02193e00*z3 - 1.35006e-1*AFE2 
    b[60] =-5.31271e00 - 6.62489e00*z - 4.74497e00*z2 - 2.97616e00*z3 - 2.05366e-1*AFE
    b[61] = 5.81595e00 - 8.69598e00*z - 1.07450e00*z2 - 2.06903e-1*z3
    b[62] = 5.38100e-1 + 1.84685e-1*z + 4.47447e-1*z2 + 5.76257e-1*z3
    #RTMS_massive stars, end of MS
    b[63] = 1.54080e01 + 2.21700e00*z + 8.20890e00*z2
    b[64] =-7.44410e00 - 3.90170e00*z - 3.96250e00*z2
    b[65] = 2.56010e01 + 3.91800e00*z + 1.15320e00*z2
    b[66] =-2.81160e-1 - 8.91678e-2*z + 7.88068e-2*z2
    #AGE_TMS, deltaT for massive stars, TMS minus Rmax
    b[67] =-1.05390e-1 + 2.01633e-1*z + 5.53316e-1*z2 + 1.90185e-1*z3 + 2.97669e-2*AFE
    b[68] = 3.89793e-3 + 8.72170e-3*z + 9.33672e-3*z2 + 7.09167e-3*z3 + 3.89436e-4*AFE
    
    return b

def calc_R0_BYRNE(ZS,afe,M):

    # ZS = Linear, Iron-scaled Metallicity Mass Fraction
    # ZL = log-scaled metallicity log(ZS/0.02)

    if isinstance(ZS,float):  ZS  = np.zeros(1)+ZS
    if isinstance(afe,float): afe = np.zeros(1)+afe
    if isinstance(M,float):   M   = np.zeros(1)+M
    if isinstance(M,int):     M   = np.zeros(1)+M

    if len(ZS) == 1:  ZS =np.zeros(len(M))+ZS
    if len(afe) == 1: afe=np.zeros(len(M))+afe

    R0 = np.zeros(len(M))
    for i in range(0,len(M)):
    
        b = calc_coeffs_BYRNE(ZS[i],afe[i])
        
        R0a= (b[0]+(b[1]*M[i]**b[2])+(b[3]*M[i]**b[4])+(b[5]*M[i]**b[6]))/((b[7]*M[i]**b[8])+(b[9]*M[i]**b[10]))
        R0b= (b[11])+(b[12]*M[i]**b[13])+(b[14]*M[i]**b[15])
        # R0lm=(b[0]+(b[1]*M**b[2]+(b[3]*M**b[4]+(b[5]*M**b[6])/((b[7]*M**b[8]+(b[9]*M**b[10])
    
        if(M[i]<=3.4): R0[i]=R0a
        if(M[i]>=3.5): R0[i]=R0b
        if((M[i]>3.4) & (M[i]<3.5)):
            R0a= (b[0]+(b[1]*3.4**b[2])+(b[3]*3.4**b[4])+(b[5]*3.4**b[6]))/((b[7]*3.4**b[8])+(b[9]*3.4**b[10]))
            R0b= (b[11])+(b[12]*3.5**b[13])+(b[14]*3.5**b[15])
            R0[i]=((3.5-M[i])*R0a+(M[i]-3.4)*R0b)*10.0
    
    return R0

def calc_L0_BYRNE(ZS,afe,M):

    # ZS = Linear, Iron-scaled Metallicity Mass Fraction
    # ZL = log-scaled metallicity log(ZS/0.02)

    if isinstance(ZS,float):  ZS  = np.zeros(1)+ZS
    if isinstance(afe,float): afe = np.zeros(1)+afe
    if isinstance(M,float):   M   = np.zeros(1)+M
    if isinstance(M,int):     M   = np.zeros(1)+M

    if len(ZS) == 1:  ZS =np.zeros(len(M))+ZS
    if len(afe) == 1: afe=np.zeros(len(M))+afe

    L0 = np.zeros(len(M))
    for i in range(0,len(M)):
    
        b = calc_coeffs_BYRNE(ZS[i],afe[i])

        L0[i]=(b[16]*M[i]**5.5+b[17]*M[i]**11.0)/(b[18]+M[i]**3.0+b[19]*M[i]**5.0+b[20]*M[i]**7.0+b[21]*M[i]**8.0+b[22]*M[i]**9.5)
        
    return L0

def calc_tMS_BYRNE(ZS,afe,M):

    # ZS = Linear, Iron-scaled Metallicity Mass Fraction
    # ZL = log-scaled metallicity log(ZS/0.02)

    if isinstance(ZS,float):  ZS  = np.zeros(1)+ZS
    if isinstance(afe,float): afe = np.zeros(1)+afe
    if isinstance(M,float):   M   = np.zeros(1)+M
    if isinstance(M,int):     M   = np.zeros(1)+M

    if len(ZS) == 1:  ZS =np.zeros(len(M))+ZS
    if len(afe) == 1: afe=np.zeros(len(M))+afe

    tMS = np.zeros(len(M))
    tBGB= np.zeros(len(M))
    tWR = np.zeros(len(M))
    
    for i in range(0,len(M)):

        ZL = calc_ZL(ZS[i])
        a = calc_coeffs_HURLEY(ZS[i])
        b = calc_coeffs_BYRNE(ZS[i],afe[i])
        
        tBGB[i]=(b[23]+b[24]*M[i]**(4.0+b[25])+b[26]*M[i]**(5.5+b[27])+b[28]*M[i]**7.0)/(a[4]*M[i]**2+a[5]*M[i]**7.0)
        tX=np.maximum(0.95,np.minimum(0.95-0.03*(ZL+0.30103),0.99))
        tY=np.maximum(0.5,1.0-0.01*np.maximum(a[6]/(M[i]**a[7]),a[8]+a[9]/(M[i]**a[10])))
        tHOOK=tY*tBGB[i]
        tMS[i]=np.maximum(tHOOK,tX*tBGB[i])

        deltaT = np.maximum(0,b[67]+M[i]*b[68])
        tWR[i] = tMS[i]-deltaT
        
    return tMS, tWR

def calc_RTMS_BYRNE(ZS,afe,M):

    # ZS = Linear, Iron-scaled Metallicity Mass Fraction
    # ZL = log-scaled metallicity log(ZS/0.02)

    if isinstance(ZS,float):  ZS  = np.zeros(1)+ZS
    if isinstance(afe,float): afe = np.zeros(1)+afe
    if isinstance(M,float):   M   = np.zeros(1)+M
    if isinstance(M,int):     M   = np.zeros(1)+M

    if len(ZS) == 1:  ZS =np.zeros(len(M))+ZS
    if len(afe) == 1: afe=np.zeros(len(M))+afe

    RTMS= np.zeros(len(M))
    RWR = np.zeros(len(M))
    lll = str(len(M))
    
    for i in range(0,len(M)):

        # print(str(i)+'/'+lll, end='\r')
        ZL = calc_ZL(ZS[i])
        a = calc_coeffs_HURLEY(ZS[i])
        b = calc_coeffs_BYRNE(ZS[i],afe[i])
        MWRR = 2164.9*(ZL+2)**(-7.0) + 64.947*(ZL+2)**(-0.8)

        RT1=(b[45]+b[46]*M[i]**b[47])/(b[48]+M[i]**b[49])
        RT2=(b[50]*M[i]**(b[51])+b[52]*M[i]**b[53]+b[54]*M[i]**b[55])/(b[56]+b[57]*M[i]**(b[58]))
        RT3=(10.**b[59])*M[i]**b[60]+b[61]*M[i]**b[62]
        RW4=(10.**b[63])*M[i]**b[64]+b[65]*M[i]**b[66]

        R0 =calc_R0_BYRNE(ZS,afe,M)

        if (M[i] < a[17]):
            RTMS[i] = RT1
            if RTMS[i] <= R0[i]:
                RTMS[i] = R0[i]

        if ((M[i]>=a[17]+0.1)&(M[i]<0.90*MWRR)):
            RTMS[i] = RT2
        if (M[i] >= 1.10*MWRR):
            RTMS[i] = RT3

        if ((M[i]>=a[17])&(M[i]<a[17]+0.1)):
            RT1=(b[45]+b[46]*a[17]**b[47])/(b[48]+a[17]**b[49])
            RT2=(b[50]*(a[17]+0.1)**(b[51])+b[52]*(a[17]+0.1)**b[53]+b[54]*(a[17]+0.1)**b[55])/(b[56]+b[57]*(a[17]+0.1)**(b[58]))
            RTMS[i]=(((a[17]+0.1)-M[i])*RT1+(M[i]-a[17])*RT2)*10.0

        if ((M[i]>=MWRR*0.90)&(M[i]<=MWRR*1.10)):
            RT2i=(b[50]*(MWRR*0.90)**(b[51])+b[52]*(MWRR*0.90)**b[53]+b[54]*(MWRR*0.90)**b[55])/(b[56]+b[57]*(MWRR*0.90)**(b[58]))
            RT3i=(10.**b[59])*(MWRR*1.10)**b[60]+b[61]*(MWRR*1.10)**b[62]
            ip_alpha = np.log10(RT3i / RT2i) / np.log10(1.10/0.9)
            x1 = MWRR*0.9
            # print("M =", M[i])
            # print("RT2(0.9MWRR) =", RT2i)
            # print("RT3(1.1MWRR) =", RT3i)
            # print("ratio RT3/RT2 =", RT3i/RT2i)
            RTMS[i] = RT2i * (M[i] / x1)**ip_alpha
            # RTMS[i]=((1.10*MWRR-M[i])*RT2+(M[i]-0.90*MWRR)*RT3)/(0.20*MWRR)
        
        if M[i] >= MWRR:
            RWR[i]  = RW4
            if (RWR[i] > RTMS[i]):
                RWR[i] = RTMS[i]
        
        if M[i] <= a[62]:
            RTMS[i] = np.maximum(1.5*R0[i], RTMS[i])

    return RTMS,RWR

def calc_LTMS_BYRNE(ZS,afe,M):

    # ZS = Linear, Iron-scaled Metallicity Mass Fraction
    # ZL = log-scaled metallicity log(ZS/0.02)

    if isinstance(ZS,float):  ZS  = np.zeros(1)+ZS
    if isinstance(afe,float): afe = np.zeros(1)+afe
    if isinstance(M,float):   M   = np.zeros(1)+M
    if isinstance(M,int):     M   = np.zeros(1)+M

    if len(ZS) == 1:  ZS =np.zeros(len(M))+ZS
    if len(afe) == 1: afe=np.zeros(len(M))+afe


    LTMS= np.zeros(len(M))
    LWR = np.zeros(len(M))
    
    for i in range(0,len(M)):

        ZL = calc_ZL(ZS[i])
        a = calc_coeffs_HURLEY(ZS[i])
        b = calc_coeffs_BYRNE(ZS[i],afe[i])
        MWRL = 694.252*(ZL+2)**(-3.0)

        LT1 = ((b[29]*M[i]**3.)+(b[30]*M[i]**4.)+(b[31]*M[i]**(b[32]+b[33])))/((b[34])+(b[35]*M[i]**(5.))+(b[36]*M[i]**(b[32]))+(b[37]*M[i]**b[38]))
        LT2 = ((b[39]*M[i]**b[40])+(b[41]*M[i]**b[42])+(b[43]*M[i]**b[44]))
        
        LTMS[i] = LT1
        if M[i] >= MWRL:
            LWR[i]  = LT2
        
    return LTMS,LWR