https://github.com/sevenian3/ChromaStarPy
Tip revision: 103d3d0df6d9574c49818f149e1cae8100455d10 authored by Ian Short on 06 July 2023, 18:09:20 UTC
Create ReadMe
Create ReadMe
Tip revision: 103d3d0
ScaleT4250g20.py
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 21 16:39:41 2017
/**
* Initializes and re-scales a Phoenix LTE spherical reference model of
* Teff=4250K, log(g)=2.0, [Fe/H]=0.0, xi=1.0 km/s, l=1.0H_p, M=1M_Sun, R=6.4761e+10cm
*
* @author Ian
*/
@author: ishort
"""
import math
import ToolBox
import Useful
def phxRefTeff():
return 4250.0
def phxRefLogEg():
return math.log(10.0) * 2.0 #//base e!
#//He abundance from Grevesse Asplund et al 2010
def phxRefLogAHe():
return math.log(10.0) * (10.93 - 12.0) #//base e "A_12" logarithmic abundance scale!
def getphxRefTau64():
#//Corresponding Tau_12000 grid (ie. lambda_0 = 1200 nm):
phxRefTau64 = [
0.00000000000000000e+00, 9.99999999999999955e-07, 1.34596032415536424e-06,
1.81160919420041334e-06, 2.43835409826882661e-06, 3.28192787251147086e-06,
4.41734470314007309e-06, 5.94557070854439435e-06, 8.00250227816105150e-06,
1.07710505603676914e-05, 1.44974067037263169e-05, 1.95129342263596216e-05,
2.62636352765333530e-05, 3.53498110503010939e-05, 4.75794431400941383e-05,
6.40400427119728238e-05, 8.61953566475303262e-05, 1.16015530173997159e-04,
1.56152300600049659e-04, 2.10174801133248699e-04, 2.82886943462596935e-04,
3.80754602122237182e-04, 5.12480587696093125e-04, 6.89778537938765847e-04,
9.28414544519474451e-04, 1.24960914129198684e-03, 1.68192432488086874e-03,
2.26380340952144670e-03, 3.04698957090350801e-03, 4.10112707055130046e-03,
5.51995432128156785e-03, 7.42963950759494875e-03, 1.00000000000000002e-02,
1.34596032415536422e-02, 1.81160919420041318e-02, 2.43835409826882663e-02,
3.28192787251147047e-02, 4.41734470314006436e-02, 5.94557070854439401e-02,
8.00250227816105275e-02, 1.07710505603676912e-01, 1.44974067037263149e-01,
1.95129342263596212e-01, 2.62636352765332981e-01, 3.53498110503010221e-01,
4.75794431400941464e-01, 6.40400427119728333e-01, 8.61953566475303190e-01,
1.16015530173997150e+00, 1.56152300600049654e+00, 2.10174801133248712e+00,
2.82886943462596641e+00, 3.80754602122236818e+00, 5.12480587696092638e+00,
6.89778537938765801e+00, 9.28414544519474383e+00, 1.24960914129198670e+01,
1.68192432488086894e+01, 2.26380340952144650e+01, 3.04698957090350540e+01,
4.10112707055129562e+01, 5.51995432128157333e+01, 7.42963950759495049e+01,
1.00000000000000000e+02]
return phxRefTau64
def getLogPhxRefTau64():
logE = math.log10(math.e)
phxRefTau64 = getphxRefTau64()
numPhxDep = len(phxRefTau64)
logPhxRefTau64 = [0.0 for i in range(numPhxDep)]
#for i in range(1, numPhxDep):
# logPhxRefTau64[i] = math.log(phxRefTau64[i])
logPhxRefTau64[1: numPhxDep] = [ math.log(phxRefTau64[i]) for i in range(1, numPhxDep) ]
#print("*********")
#print(logPhxRefTau64)
#print("*********")
#stop
logPhxRefTau64[0] = logPhxRefTau64[1] - (logPhxRefTau64[numPhxDep - 1] - logPhxRefTau64[1]) / numPhxDep
return logPhxRefTau64
def phxRefTemp(teff, numDeps, tauRos):
logE = math.log10(math.e)
#//Theoretical radiative/convective model from Phoenix V15:
phxRefTemp64 = [
2.55177189467776134e+03, 2.55177189467776134e+03, 2.57792125655286191e+03,
2.60758002435901381e+03, 2.64012742898025908e+03, 2.67432875475946230e+03,
2.70893891273552163e+03, 2.74308870598832664e+03, 2.77632383320681129e+03,
2.80848811483166946e+03, 2.83959882227651724e+03, 2.86975658734314266e+03,
2.89908813206188461e+03, 2.92771359779851400e+03, 2.95573102630948051e+03,
2.98321198037901013e+03, 3.01020749755630777e+03, 3.03675606855582828e+03,
3.06288810300240721e+03, 3.08862982043339889e+03, 3.11400422160098924e+03,
3.13903037085205142e+03, 3.16372090154785838e+03, 3.18810905092703479e+03,
3.21223709737028685e+03, 3.23616013487906321e+03, 3.25997011697764810e+03,
3.28379847029234497e+03, 3.30779792957306154e+03, 3.33215641253420881e+03,
3.35718199939999613e+03, 3.38317056293436599e+03, 3.41035215459581877e+03,
3.43904544965103469e+03, 3.46977882782005645e+03, 3.50316432370656776e+03,
3.53952093749718961e+03, 3.57938756282554868e+03, 3.62358876591721355e+03,
3.67390373594475159e+03, 3.73057002905895024e+03, 3.79526239264127798e+03,
3.87014190881368677e+03, 3.96034605960104500e+03, 4.06785077337546363e+03,
4.19530126018311603e+03, 4.35807101922509446e+03, 4.53128152603078979e+03,
4.75647877892966608e+03, 4.99815140831592271e+03, 5.29640554819846147e+03,
5.60111278999269234e+03, 6.01099379170666271e+03, 6.36043009619938857e+03,
6.87521615935859882e+03, 7.26044412747461593e+03, 7.50448146936407466e+03,
7.68547260686038317e+03, 7.83660366095023619e+03, 7.97031461253075395e+03,
8.09288613904806061e+03, 8.20796577074222841e+03, 8.31800144917403668e+03,
8.42235419676150195e+03]
logPhxRefTau64 = getLogPhxRefTau64();
#// interpolate onto gS3 tauRos grid and re-scale with Teff:
phxRefTemp = [0.0 for i in range(numDeps)]
scaleTemp = [ [0.0 for i in range(numDeps)] for j in range(2)]
#for i in range(numDeps):
# phxRefTemp[i] = ToolBox.interpol(logPhxRefTau64, phxRefTemp64, tauRos[1][i])
# scaleTemp[0][i] = teff * phxRefTemp[i] / phxRefTeff()
# scaleTemp[1][i] = math.log(scaleTemp[0][i]);
phxRefTemp = [ ToolBox.interpol(logPhxRefTau64, phxRefTemp64, x) for x in tauRos[1] ]
scaleTemp[0] = [ teff * x / phxRefTeff() for x in phxRefTemp ]
scaleTemp[1] = [ math.log(x) for x in scaleTemp[0] ]
#//System.out.println("tauRos[1][i] " + logE * tauRos[1][i] + " scaleTemp[1][i] " + logE * scaleTemp[1][i]);
#print("*********")
#print(phxRefTemp)
#print("*********")
#print(scaleTemp[0])
#stop
return scaleTemp
def phxRefPGas(grav, zScale, logAHe, numDeps, tauRos):
#//System.out.println("ScaleT4250g20.phxRefPGas called");
logE = math.log10(math.e)
logEg = math.log(grav) #//base e!
AHe = math.exp(logAHe)
refAHe = math.exp(phxRefLogAHe())
logZScale = math.log(zScale)
#//Theoretical radiative/convective model from Phoenix V15:
phxRefPGas64 = [
1.00000000000000005e-04, 4.30797022881529035e+00, 5.23633209862413107e+00,
6.35110679837766412e+00, 7.68415255320091806e+00, 9.26944725328996455e+00,
1.11437944000671951e+01, 1.33486175933519231e+01, 1.59320732110696728e+01,
1.89510405786159843e+01, 2.24729535630686001e+01, 2.65776121702747155e+01,
3.13591521574494898e+01, 3.69283427891380711e+01, 4.34153596323810476e+01,
5.09731596421035889e+01, 5.97815115471393099e+01, 7.00515965533140843e+01,
8.20320988196306047e+01, 9.60157114466986172e+01, 1.12346894715709283e+02,
1.31431137235228107e+02, 1.53746059612369663e+02, 1.79853925611152988e+02,
2.10416068512255976e+02, 2.46209759085281831e+02, 2.88146638387011080e+02,
3.37292455489998588e+02, 3.94889455574428723e+02, 4.62381218232488322e+02,
5.41431239027500510e+02, 6.33944346680425951e+02, 7.42101436711750239e+02,
8.68388258638872003e+02, 1.01559408923251908e+03, 1.18679360231291594e+03,
1.38539626183849145e+03, 1.61516883770226173e+03, 1.88020120421177171e+03,
2.18463232976707104e+03, 2.53293658371589027e+03, 2.93000255565556563e+03,
3.38117288723951515e+03, 3.89194415900473678e+03, 4.46986439959823201e+03,
5.12776215162745939e+03, 5.88482548647403291e+03, 6.77685699601984379e+03,
7.85508152689431518e+03, 9.16793034379196797e+03, 1.06499829584269464e+04,
1.20573723930509423e+04, 1.31104385587741126e+04, 1.38632723386838443e+04,
1.43359336214159030e+04, 1.46332769212074945e+04, 1.48729928730357842e+04,
1.50999078995603995e+04, 1.53309455833121519e+04, 1.55753421735992888e+04,
1.58400487938745937e+04, 1.61313349003330495e+04, 1.64553234176194455e+04,
1.68190844909027946e+04]
numPhxDeps = len(phxRefPGas64) #//yeah, I know, 64, but that could change!
logPhxRefPGas64 = [0.0 for i in range(numPhxDeps)]
#for i in range(numPhxDeps):
# logPhxRefPGas64[i] = math.log(phxRefPGas64[i])
logPhxRefPGas64 = [ math.log(x) for x in phxRefPGas64 ]
logPhxRefTau64 = getLogPhxRefTau64();
#// interpolate onto gS3 tauRos grid and re-scale with grav, metallicity and He abundance
#// From Gray 3rd Ed. Ch.9, esp p. 189, 196:
phxRefPGas = [0.0 for i in range(numDeps)]
logPhxRefPGas = [0.0 for i in range(numDeps)]
scalePGas = [ [0.0 for i in range(numDeps)] for j in range(2) ]
#//exponents in scaling with g:
gexpTop = 0.54 #//top of model
gexpBottom = 0.64 #//bottom of model
gexpRange = (gexpBottom - gexpTop)
tauLogRange = tauRos[1][numDeps-1] - tauRos[1][0]
#double thisGexp;
#// factor for scaling with A_He:
logHeDenom = 0.666667 * math.log(1.0 + 4.0*refAHe)
logPhxRefPGas = [ ToolBox.interpol(logPhxRefTau64, logPhxRefPGas64, x) for x in tauRos[1] ]
for i in range(numDeps):
# #//if (i%10 == 0){
# #//System.out.println("i " + i);
# #//}
# logPhxRefPGas[i] = ToolBox.interpol(logPhxRefTau64, logPhxRefPGas64, tauRos[1][i])
# #//if (i%10 == 0){
# #//System.out.println("After tau interpolation: pGas " + logE*logPhxRefPGas[i]);
# #//}
thisGexp = gexpTop + gexpRange * (tauRos[1][i] - tauRos[1][0]) / tauLogRange
# #//scaling with g
# #//if (i%10 == 0){
# #//System.out.println("thisGexp " + thisGexp);
# #//}
scalePGas[1][i] = thisGexp*logEg + logPhxRefPGas[i] - thisGexp*phxRefLogEg()
# #//if (i%10 == 0){
# #//System.out.println("After scaling with g: pGas " + logE*scalePGas[1][i]);
# #//}
# #//scaling with zscl:
# #//if (i%10 == 0){
# #//System.out.println("logZScale " + logZScale);
# #//}
# scalePGas[1][i] = -0.333333*logZScale + scalePGas[1][i]
# #//if (i%10 == 0){
# #//System.out.println("After scaling with z: pGas " + logE*scalePGas[1][i]);
# #//}
# #//scaling with A_He:
# #//if (i%10 == 0){
# #//System.out.println("Math.log(1.0 + 4.0*AHe) - logHeDenom " + (0.666667*Math.log(1.0 + 4.0*AHe) - logHeDenom));
# #//}
# scalePGas[1][i] = 0.666667 * math.log(1.0 + 4.0*AHe) + scalePGas[1][i] - logHeDenom
# #//if (i%10 == 0){
# #//System.out.println("After scaling with AHe: pGas " + logE*scalePGas[1][i]);
# #//}
# scalePGas[0][i] = math.exp(scalePGas[1][i])
# #//if (i%10 == 0){
# #//System.out.println("logPhxRefPGas " + logE*logPhxRefPGas[i] + " scalePGas[1][i] " + logE * scalePGas[1][i]);
# #//}
#logPhxRefPGas = [ ToolBox.interpol(logPhxRefTau64, logPhxRefPGas64, x) for x in tauRos[1] ]
#No! scalePGas[1] = [ logEg*(gexpTop + gexpRange * (tauRos[1][i] - tauRos[1][0]) / tauLogRange)\
#No! + logPhxRefPGas[i] - thisGexp*phxRefLogEg() for i in range(numDeps) ]
scalePGas[1] = [ -0.333333*logZScale + x for x in scalePGas[1] ]
scalePGas[1] = [ 0.666667 * math.log(1.0 + 4.0*AHe) + x - logHeDenom for x in scalePGas[1] ]
scalePGas[0] = [ math.exp(x) for x in scalePGas[1] ]
#SOMETHING WRONG!!
#print("*********")
#print("logPhxRefPGas ", logPhxRefPGas)
#print("*********")
#print("scalePGas[1] ", scalePGas[1])
#print("*********")
#stop
#Carefull here - P at upper boundary can be an underestimate, but it must not be greater than value at next depth in!
if (scalePGas[0][0] >= scalePGas[0][1]):
scalePGas[0][0] = 0.5 * scalePGas[0][1];
scalePGas[1][0] = math.log(scalePGas[0][0]);
return scalePGas
def phxRefPe(teff, grav, numDeps, tauRos, zScale, logAHe):
logE = math.log10(math.e)
logEg = math.log(grav) #//base e!
AHe = math.exp(logAHe)
refAHe = math.exp(phxRefLogAHe())
logZScale = math.log(zScale)
#//Theoretical radiative/convective model from Phoenix V15:
phxRefPe64 = [
8.82107332460937786e-09, 2.78937385278448721e-05, 3.47576301038577686e-05,
4.35832572659463317e-05, 5.49376619037973015e-05, 6.94409656870532439e-05,
8.77707903792445917e-05, 1.10693651821575575e-04, 1.39104108185757799e-04,
1.74066374927776200e-04, 2.16859921028963076e-04, 2.69028731246404777e-04,
3.32432061837431244e-04, 4.09296056659801003e-04, 5.02267654306449442e-04,
6.14473664793957226e-04, 7.49597452349998888e-04, 9.11974871131562793e-04,
1.10671154087655646e-03, 1.33981929191927902e-03, 1.61837137445693482e-03,
1.95067571498008536e-03, 2.34645801822358068e-03, 2.81730972458459810e-03,
3.37701527100230971e-03, 4.04206883240371042e-03, 4.83264944868273382e-03,
5.77375902899801303e-03, 6.89635454260574543e-03, 8.23939377538898503e-03,
9.85478499694628605e-03, 1.18084450314643753e-02, 1.41821928408006545e-02,
1.70844502273588689e-02, 2.06675266544882504e-02, 2.51401461738746321e-02,
3.07601652353682656e-02, 3.78883555684116358e-02, 4.70434551509077148e-02,
5.90737609051253110e-02, 7.49278881219457016e-02, 9.61506501246208595e-02,
1.25001054449255633e-01, 1.65469764331727026e-01, 2.21877141054763416e-01,
2.99471993844184214e-01, 4.09463450124863737e-01, 5.44807309494293679e-01,
7.32622090968494066e-01, 1.00147011151172505e+00, 1.61234892801148755e+00,
3.05928587566573817e+00, 7.65951702544590596e+00, 1.63896549667325182e+01,
4.47471354140328827e+01, 8.76262151753557248e+01, 1.30080006341931238e+02,
1.72013432855285373e+02, 2.15462867458483203e+02, 2.61494582627386649e+02,
3.10959654478989705e+02, 3.64658419590852475e+02, 4.23480004219150715e+02,
4.87036243289867400e+02]
logPhxRefTau64 = getLogPhxRefTau64()
numPhxDeps = len(phxRefPe64) #//yeah, I know, 64, but that could change!
logPhxRefPe64 = [0.0 for i in range(numPhxDeps)]
#for i in range(numPhxDeps):
# logPhxRefPe64[i] = math.log(phxRefPe64[i])
logPhxRefPe64 = [ math.log(x) for x in phxRefPe64 ]
#// interpolate onto gS3 tauRos grid and re-scale with Teff:
phxRefPe = [0.0 for i in range(numDeps)]
logPhxRefPe = [0.0 for i in range(numDeps)]
scalePe = [ [0.0 for i in range(numDeps) ] for j in range(2) ]
#//exponents in scaling with Teff ONLY VALID FOR Teff < 10000K:
omegaTaum1 = 0.0012 #//log_10(tau) < 0.1
omegaTaup1 = 0.0015 #//log_10(tau) > 1.0
omegaRange = (omegaTaup1-omegaTaum1)
lonOfM1 = math.log(0.1)
#//exponents in scaling with g:
gexpTop = 0.48 #//top of model
gexpBottom = 0.33 #//bottom of model
gexpRange = (gexpBottom - gexpTop)
tauLogRange = tauRos[1][numDeps-1] - tauRos[1][0]
#double thisGexp
thisOmega = omegaTaum1 #//default initialization
#// factor for scaling with A_He:
logHeDenom = 0.333333 * math.log(1.0 + 4.0*refAHe)
logPhxRefPe = [ ToolBox.interpol(logPhxRefTau64, logPhxRefPe64, x) for x in tauRos[1] ]
for i in range(numDeps):
#//if (i%10 == 0){
#//System.out.println("i " + i);
#//}
#logPhxRefPe[i] = ToolBox.interpol(logPhxRefTau64, logPhxRefPe64, tauRos[1][i])
thisGexp = gexpTop + gexpRange * (tauRos[1][i] - tauRos[1][0]) / tauLogRange
if (tauRos[0][i] < 0.1):
thisOmega = omegaTaum1
if (tauRos[0][i] > 10.0):
thisOmega = omegaTaup1
if ( (tauRos[0][i] >= 0.1) and (tauRos[0][i] <= 10.0) ):
thisOmega = omegaTaum1 + omegaRange * (tauRos[1][i] - lonOfM1) / tauLogRange
#//if (i%10 == 0){
#//System.out.println("thisGexp " + thisGexp + " thisOmega " + thisOmega);
#//}
#//scaling with g
scalePe[1][i] = thisGexp*logEg + logPhxRefPe[i] - thisGexp*phxRefLogEg()
#//if (i%10 == 0){
#//System.out.println("After g scaling: pe " + logE*scalePe[1][i]);
#//}
#//scale with Teff:
scalePe[1][i] = thisOmega*teff + scalePe[1][i] - thisOmega*phxRefTeff()
#//if (i%10 == 0){
#//System.out.println("After Teff scaling: pe " + logE*scalePe[1][i]);
#//}
#//scaling with zscl:
#scalePe[1][i] = 0.333333*logZScale + scalePe[1][i]
#//if (i%10 == 0){
#//System.out.println("After z scaling: pe " + logE*scalePe[1][i]);
#//}
#//scaling with A_He:
#scalePe[1][i] = 0.333333 * math.log(1.0 + 4.0*AHe) + scalePe[1][i] - logHeDenom
#//if (i%10 == 0){
#//System.out.println(" logPhxRefPe " + logE*logPhxRefPe[i] + " After A_He scaling: pe " + logE*scalePe[1][i]);
#//}
#scalePe[0][i] = math.exp(scalePe[1][i]);
scalePe[1] = [ 0.333333*logZScale + x for x in scalePe[1] ]
scalePe[1] = [ 0.333333 * math.log(1.0 + 4.0*AHe) + x - logHeDenom for x in scalePe[1] ]
scalePe[0] = [ math.exp(x) for x in scalePe[1] ]
#print("*********")
#print("logPhxRefPe ", logPhxRefPe)
#print("*********")
#print("scalePe[1] ", scalePe[1])
#print("*********")
#print("scalePe[0] ", scalePe[0])
#print("*********")
#stop
return scalePe
def phxRefNe(numDeps, scaleTemp, scalePe):
logE = math.log10(math.e)
scaleNe = [ [0.0 for i in range(numDeps) ] for j in range(2) ]
#for i in range(numDeps):
# scaleNe[1][i] = scalePe[1][i] - scaleTemp[1][i] - Useful.logK()
# scaleNe[0][i] = math.exp(scaleNe[1][i])
scaleNe[1] = [ scalePe[1][i] - scaleTemp[1][i] - Useful.logK() for i in range(numDeps) ]
scaleNe[0] = [ math.exp(x) for x in scaleNe[1] ]
#print("*********")
#print("scaleNe[1] ", scaleNe[1])
#print("*********")
#print("scaleNe[0] ", scaleNe[0])
#print("*********")
#stop
return scaleNe
