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
LineGrid.py
# -*- coding: utf-8 -*-
"""
Created on Sat Apr 29 13:31:10 2017
@author: Ian
"""
import math
import Useful
"""/**
* Line profile, phi_lambda(lambda): Assume Voigt function profile - need H(a,v)
* Assumes CRD, LTE, ??? Input parameters: lam0 - line center wavelength in nm
* mass - mass of absorbing particle (amu) logGammaCol - log_10(gamma) - base 10
* logarithmic collisional (pressure) damping co-efficient (s^-1) epsilon -
* convective microturbulence- non-thermal broadening parameter (km/s) Also
* needs atmospheric structure information: numDeps WON'T WORK - need observer's
* frame fixed lambda at all depths: temp structure for depth-dependent thermal
* line broadening Teff as typical temp instead of above pressure structure,
* pGas, if scaling gamma
*/"""
def lineGridDelta(lam0In, massIn, xiTIn, numDeps, teff):
c = Useful.c()
logC = Useful.logC()
#//double k = Useful.k;
logK = Useful.logK()
#//double e = Useful.e;
#//double mE = Useful.mE;
amu = Useful.amu()
ln10 = math.log(10.0)
ln2 = math.log(2.0)
logE = math.log10(math.e) #// for debug output
#//Put input parameters into linear cgs units:
#//double gammaCol = Math.pow(10.0, logGammaCol);
logTeff = math.log(teff)
xiT = xiTIn * 1.0E5 #//km/s to cm/s
lam0 = lam0In #// * 1.0E-7; //nm to cm
logLam0 = math.log(lam0)
logMass = math.log(massIn * amu) #//amu to g
#// Compute depth-independent Doppler width, Delta_lambda_D:
#double doppler, logDopp;
#double logHelp, help; //scratch
logHelp = ln2 + logK + logTeff - logMass #// M-B dist, square of v_mode
help = math.exp(logHelp) + xiT * xiT #// quadratic sum of thermal v and turbulent v
logHelp = 0.5 * math.log(help)
logDopp = logHelp + logLam0 - logC
doppler = math.exp(logDopp) #// cm
#//System.out.println("LineGrid: doppler, logDopp: " + doppler + " " + logE*logDopp);
#//Set up a half-profile Delta_lambda grid in Doppler width units
#//from line centre to wing
#//int numCore = 5;
#//int numWing = 5;
#//int numWing = 0; //debug
numPoints = 1
#// a 2D 2 X numPoints array of Delta Lambdas
#// Row 0 : Delta lambdas in cm - will need to be in nm for Planck and Rad Trans?
#// Row 1 : Delta lambdas in Doppler widths
linePoints = [ [ 0.0 for i in range(numPoints) ] for j in range(2) ]
#// Line profiel points in Doppler widths - needed for Voigt function, H(a,v):
v = [ 0.0 for i in range(numPoints) ]
#double logV, ii, jj;
il = 0
ii = float(il)
#// In core, space v points linearly:
#// Voigt "v" parameter
#// v > 0 --> This is the *red* wing:
v[il] = ii
linePoints[0][il] = doppler * v[il]
linePoints[1][il] = v[il]
#//System.out.println("LineGrid: il, lam, v: " + il + " " +
#// linePoints[0][il] + " " + linePoints[1][il]);
return linePoints
#} //end method lineGridDelta
#//
#//
#//
def lineGridGauss(lam0In, massIn, xiTIn, numDeps, teff, numCore):
c = Useful.c()
logC = Useful.logC()
#//double k = Useful.k;
logK = Useful.logK()
#//double e = Useful.e;
#//double mE = Useful.mE;
amu = Useful.amu()
dln10 = math.log(10.0)
ln2 = math.log(2.0)
logE = math.log10(math.e) #// for debug output
#//Put input parameters into linear cgs units:
#//double gammaCol = Math.pow(10.0, logGammaCol);
logTeff = math.log(teff)
xiT = xiTIn * 1.0E5 #//km/s to cm/s
lam0 = lam0In #// * 1.0E-7; //nm to cm
logLam0 = math.log(lam0)
logMass = math.log(massIn * amu) #//amu to g
#// Compute depth-independent Doppler width, Delta_lambda_D:
#double doppler, logDopp;
#double logHelp, help; //scratch
logHelp = ln2 + logK + logTeff - logMass #// M-B dist, square of v_mode
help = math.exp(logHelp) + xiT * xiT #// quadratic sum of thermal v and turbulent v
logHelp = 0.5 * math.log(help)
logDopp = logHelp + logLam0 - logC
doppler = math.exp(logDopp) #// cm
#//System.out.println("LineGrid: doppler, logDopp: " + doppler + " " + logE*logDopp);
#//Set up a half-profile Delta_lambda grid in Doppler width units
#//from line centre to wing
#//int numCore = 5;
#//int numWing = 5;
#//int numWing = 0; //debug
numPoints = numCore
#// a 2D 2 X numPoints array of Delta Lambdas
#// Row 0 : Delta lambdas in cm - will need to be in nm for Planck and Rad Trans?
#// Row 1 : Delta lambdas in Doppler widths
linePoints = [ [ 0.0 for i in range(numPoints) ] for j in range(2) ]
#// Line profiel points in Doppler widths - needed for Voigt function, H(a,v):
v = [0.0 for i in range(numPoints)]
maxCoreV = 3.5 #//core half-width ~ in Doppler widths
#//double maxWingDeltaLogV = 1.5 * ln10; //maximum base e logarithmic shift from line centre in Doppler widths
minWingDeltaLogV = math.log(maxCoreV + 1.5)
maxWingDeltaLogV = 9.0 + minWingDeltaLogV
#double logV, ii, jj;
for il in range(numPoints):
ii = float(il)
#// In core, space v points linearly:
#// Voigt "v" parameter
#// v > 0 --> This is the *red* wing:
v[il] = ii * maxCoreV / (numCore - 1)
linePoints[0][il] = doppler * v[il]
linePoints[1][il] = v[il]
#//System.out.println("LineGrid: il, lam, v: " + il + " " +
#// linePoints[0][il] + " " + linePoints[1][il]);
#} // il lambda loop
#// Add the negative DeltaLambda half of the line:
numPoints2 = (2 * numPoints) - 1
#//System.out.println("LineGrid: numpoints2: " + numPoints2);
#// Return a 2D 2 X (2xnumPoints-1) array of Delta Lambdas
#// Row 0 : Delta lambdas in cm - will need to be in nm for Planck and Rad Trans?
#// Row 1 : Delta lambdas in Doppler widths
linePoints2 = [ [ 0.0 for i in range(numPoints2) ] for j in range(2) ]
#//wavelengths are depth-independent - just put them in the 0th depth slot:
for il2 in range(numPoints2):
if (il2 < numPoints - 1):
il = (numPoints - 1) - il2
linePoints2[0][il2] = -1.0 * linePoints[0][il]
linePoints2[1][il2] = -1.0 * linePoints[1][il]
else:
#//Positive DelataLambda half:
il = il2 - (numPoints - 1)
linePoints2[0][il2] = linePoints[0][il]
linePoints2[1][il2] = linePoints[1][il]
#//System.out.println("LineGrid: il2, lam, v: " + il2 + " " +
#// linePoints2[0][il2] + " " + linePoints2[1][il2]);
#} //il2 loop
return linePoints2
#} //end method lineGridGauss
#//
#//
#//
def lineGridVoigt(lam0In, massIn, xiTIn, numDeps, teff, numCore, numWing, species):
c = Useful.c()
logC = Useful.logC()
#//double k = Useful.k;
logK = Useful.logK()
#//double e = Useful.e;
#//double mE = Useful.mE;
amu = Useful.amu()
ln10 = math.log(10.0)
ln2 = math.log(2.0)
logE = math.log10(math.e) #// for debug output
#//Put input parameters into linear cgs units:
#//double gammaCol = Math.pow(10.0, logGammaCol);
logTeff = math.log(teff)
xiT = xiTIn * 1.0E5 #//km/s to cm/s
lam0 = lam0In #// * 1.0E-7; //nm to cm
logLam0 = math.log(lam0)
logMass = math.log(massIn * amu) #//amu to g
#// Compute depth-independent Doppler width, Delta_lambda_D:
#double doppler, logDopp;
#double logHelp, help; //scratch
logHelp = ln2 + logK + logTeff - logMass #// M-B dist, square of v_mode
help = math.exp(logHelp) + xiT * xiT #// quadratic sum of thermal v and turbulent v
logHelp = 0.5 * math.log(help)
logDopp = logHelp + logLam0 - logC
doppler = math.exp(logDopp) #// cm
#//System.out.println("LineGrid: doppler, logDopp: " + doppler + " " + logE*logDopp);
#//Set up a half-profile Delta_lambda grid in Doppler width units
#//from line centre to wing
#//int numCore = 5;
#//int numWing = 5;
#//int numWing = 0; //debug
numPoints = numCore + numWing
#// a 2D 2 X numPoints array of Delta Lambdas
#// Row 0 : Delta lambdas in cm - will need to be in nm for Planck and Rad Trans?
#// Row 1 : Delta lambdas in Doppler widths
linePoints = [ [ 0.0 for i in range(numPoints) ] for j in range(2) ]
#// Line profiel points in Doppler widths - needed for Voigt function, H(a,v):
v = [0.0 for i in range(numPoints) ]
maxCoreV = 3.5 #//core half-width ~ in Doppler widths
#//double maxWingDeltaLogV = 1.5 * ln10; //maximum base e logarithmic shift from line centre in Doppler widths
minWingDeltaLogV = math.log(maxCoreV + 1.5)
maxWingDeltaLogV = 9.0 + minWingDeltaLogV
if(species=="HI" and teff>=7000):
maxCoreV = 3.5
minWingDeltaLogV = math.log(maxCoreV + 1.5)
maxWingDeltaLogV = 12.0 + minWingDeltaLogV
#//console.log("2)"+maxWingDeltaLogV);
#double logV, ii, jj;
for il in range(numPoints):
ii = float(il)
if (il < numCore):
#// In core, space v points linearly:
#// Voigt "v" parameter
#// v > 0 --> This is the *red* wing:
v[il] = ii * maxCoreV / (numCore - 1)
linePoints[0][il] = doppler * v[il]
linePoints[1][il] = v[il]
else:
#//Space v points logarithmically in wing
jj = ii - numCore
logV = (jj * (maxWingDeltaLogV - minWingDeltaLogV) / (numPoints - 1)) + minWingDeltaLogV
v[il] = math.exp(logV)
linePoints[0][il] = doppler * v[il]
linePoints[1][il] = v[il]
#} // end else
#//System.out.println("LineGrid: il, lam, v: " + il + " " +
#// linePoints[0][il] + " " + linePoints[1][il]);
#} // il lambda loop
#// Add the negative DeltaLambda half of the line:
numPoints2 = (2 * numPoints) - 1
#//System.out.println("LineGrid: numpoints2: " + numPoints2);
#// Return a 2D 2 X (2xnumPoints-1) array of Delta Lambdas
#// Row 0 : Delta lambdas in cm - will need to be in nm for Planck and Rad Trans?
#// Row 1 : Delta lambdas in Doppler widths
linePoints2 = [ [ 0.0 for i in range(numPoints2) ] for j in range(2) ]
#//wavelengths are depth-independent - just put them in the 0th depth slot:
for il2 in range(numPoints2):
if (il2 < numPoints - 1):
il = (numPoints - 1) - il2
linePoints2[0][il2] = -1.0 * linePoints[0][il]
linePoints2[1][il2] = -1.0 * linePoints[1][il]
else:
#//Positive DelataLambda half:
il = il2 - (numPoints - 1)
linePoints2[0][il2] = linePoints[0][il]
linePoints2[1][il2] = linePoints[1][il]
#//System.out.println("LineGrid: il2, lam, v: " + il2 + " " +
#// linePoints2[0][il2] + " " + linePoints2[1][il2]);
#} //il2 loop
return linePoints2