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
ScaleT5000.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=5000K, log(g)=4.5, [Fe/H]=0.0, xi=1.0 km/s, l=1.0H_p, M=1M_Sun, R=6.4761D+10cm
*
* @author Ian
*/
@author: ishort
"""
import math
import ToolBox
import Useful
def phxRefTeff():
return 5000.0
def phxRefLogEg():
return math.log(10.0) * 4.5 #//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) ]
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 = [
3.15213572679982190e+03, 3.15213572679982190e+03, 3.17988621810632685e+03,
3.21012887128011243e+03, 3.24126626267038500e+03, 3.27276078893546673e+03,
3.30435725697820226e+03, 3.33589185632140106e+03, 3.36724151725549154e+03,
3.39831714195318273e+03, 3.42906935013664861e+03, 3.45949368388945595e+03,
3.48962758169505923e+03, 3.51953742647688796e+03, 3.54929791042697934e+03,
3.57896962155466872e+03, 3.60858205550851335e+03, 3.63812646699481775e+03,
3.66755983657917068e+03, 3.69681905522719444e+03, 3.72583932497757132e+03,
3.75457006928661031e+03, 3.78298372918123914e+03, 3.81109104721021231e+03,
3.83893072914395862e+03, 3.86656355962043835e+03, 3.89408059675027425e+03,
3.92160316230741546e+03, 3.94927225929978204e+03, 3.97726284805320847e+03,
4.00584847611869327e+03, 4.03531360317989993e+03, 4.06591896438200047e+03,
4.09802860937899732e+03, 4.13221207874272022e+03, 4.16915227717330799e+03,
4.20937593060261861e+03, 4.25369220113429128e+03, 4.30330739566306784e+03,
4.36035870964639616e+03, 4.42601579216115442e+03, 4.50281614584142153e+03,
4.59386420090837146e+03, 4.70448179136501403e+03, 4.83727710376560208e+03,
4.99516189027659129e+03, 5.19102132587796405e+03, 5.40505223548941285e+03,
5.67247302987449984e+03, 5.95695843497286933e+03, 6.27957483223234703e+03,
6.71365960956718118e+03, 7.06828382342861460e+03, 7.34157936910693206e+03,
7.56939938735570740e+03, 7.77138428264261165e+03, 7.95656000812699585e+03,
8.13006721530056711e+03, 8.29523535580475982e+03, 8.45429779465689171e+03,
8.60879260449185131e+03, 8.75981713693203528e+03, 8.90838141718757288e+03,
9.05361290415211806e+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]);
return scaleTemp
def phxRefPGas(grav, zScale, logAHe, numDeps, tauRos):
#//System.out.println("ScaleT5000.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, 1.03770217591881035e+02, 1.24242770084417913e+02,
1.47686628640383276e+02, 1.74578854906314291e+02, 2.05506972274478784e+02,
2.41168221287605292e+02, 2.82385081738383917e+02, 3.30127686150304896e+02,
3.85540773715381306e+02, 4.49974446823229414e+02, 5.25018679681323647e+02,
6.12542265074691159e+02, 7.14737800095933608e+02, 8.34175243666085407e+02,
9.73867213356324669e+02, 1.13734973870022168e+03, 1.32878148706864113e+03,
1.55306409432270971e+03, 1.81598529465124716e+03, 2.12438618583220841e+03,
2.48635477283421324e+03, 2.91145034581766595e+03, 3.41095942605562823e+03,
3.99819276314161607e+03, 4.68883438023894087e+03, 5.50134310662684311e+03,
6.45741052408807354e+03, 7.58249196327514983e+03, 8.90641248566333525e+03,
1.04639741154490002e+04, 1.22956502717452295e+04, 1.44484787849992390e+04,
1.69769301182948657e+04, 1.99435621814443475e+04, 2.34195796692420117e+04,
2.74860930366683497e+04, 3.22351125605895031e+04, 3.77699103578024442e+04,
4.42033085085744533e+04, 5.16616495136288213e+04, 6.02879692077906366e+04,
7.02475218656768702e+04, 8.17365047611011832e+04, 9.50146489805318997e+04,
1.10441316485543124e+05, 1.28451318144638804e+05, 1.49415613553191157e+05,
1.72877372164747008e+05, 1.96852852539717947e+05, 2.18808320050485723e+05,
2.35794833242603316e+05, 2.48716041541587241e+05, 2.59902150512206339e+05,
2.70560370352023339e+05, 2.81251297069544089e+05, 2.92310802132537181e+05,
3.03988239352240635e+05, 3.16495216131040419e+05, 3.30029076402488339e+05,
3.44786943994771456e+05, 3.60975297786138486e+05, 3.78815092131546407e+05,
3.98560549755298765e+05]
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]);
#//}
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] ]
#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 = [
1.17858427569630401e-08, 1.73073837795169436e-03, 2.13762360059438538e-03,
2.64586145846806451e-03, 3.26749020460433354e-03, 4.02219945676032288e-03,
4.93454747856481805e-03, 6.03357965110110344e-03, 7.35319802933484621e-03,
8.93306098318919460e-03, 1.08200092390451780e-02, 1.30700158082515377e-02,
1.57505131367194594e-02, 1.89428593874781982e-02, 2.27446519479000651e-02,
2.72716961646799864e-02, 3.26596927620770305e-02, 3.90659173672675136e-02,
4.66713907010225318e-02, 5.56843086932707065e-02, 6.63452384304821230e-02,
7.89341909634427297e-02, 9.37792909747245523e-02, 1.11270186635302790e-01,
1.31870014183696899e-01, 1.56130489360824298e-01, 1.84715397349025645e-01,
2.18428766543559472e-01, 2.58245610307223983e-01, 3.05363622257444900e-01,
3.61311333509324040e-01, 4.27990544717643029e-01, 5.07743853690445168e-01,
6.03604039632526179e-01, 7.19674246257567152e-01, 8.61422066803848585e-01,
1.03568172049434559e+00, 1.25187412720684454e+00, 1.52336996895144261e+00,
1.87078029858400652e+00, 2.31893413667797388e+00, 2.90597658045488094e+00,
3.68566481623166187e+00, 4.74110273402785865e+00, 6.16546324347510222e+00,
8.08486709272609971e+00, 1.07959796585076546e+01, 1.46390000057528482e+01,
2.17273927465764913e+01, 3.56194058574816239e+01, 6.57361652682183575e+01,
1.48468954779851543e+02, 2.80489497081349555e+02, 4.46587250419467807e+02,
6.46784311972032128e+02, 8.86744838282462069e+02, 1.17244960918767083e+03,
1.51089748714632174e+03, 1.91050957850908458e+03, 2.38115682377229541e+03,
2.93426662234414562e+03, 3.58305801646245618e+03, 4.34379670059742239e+03,
5.22642525609140284e+03]
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] ]
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] ]
return scaleNe
