# -*- coding: utf-8 -*-
"""
Created on Sat Apr 29 14:26:42 2017
@author: Ian
"""
import math
import Useful
import ToolBox
"""/**
* 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 delta(linePoints, lam0In, numDeps, tauRos, massIn, xiTIn, teff):
"""//delta function line profile for initiali check of line strength"""
lam0 = lam0In #// * 1.0E-7; //nm to cm
logLam0 = math.log(lam0)
logE = math.log10(math.e) #// for debug output
#//System.out.println("LineProf: doppler, logDopp: " + doppler + " " + logE*logDopp);
#//Put input parameters into linear cgs units:
#//System.out.println("LINEGRID: Tau1: " + tau1);
#//logA = 2.0 * logLam0 + logGamma - ln4pi - logC - logDopp;
#//a = Math.exp(logA);
#//System.out.println("LINEGRID: logA: " + logE * logA);
#//Set up a half-profile Delta_lambda grid in Doppler width units
#//from line centre to wing
numPoints = 1
#//System.out.println("LineProf: numPoints: " + numPoints);
#// Return a 2D numPoints X numDeps array of normalized line profile points (phi)
lineProf = [ [ 0.0 for i in range(numDeps) ] for j in range(1) ]
c = Useful.c()
logC = Useful.logC()
logK = Useful.logK()
amu = Useful.amu()
ln10 = math.log(10.0)
ln2 = math.log(2.0)
lnSqRtPi = 0.5 * math.log(math.pi)
logTeff = math.log(teff)
xiT = xiTIn * 1.0E5 #//km/s to cm/s
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
#// Line profile points in Doppler widths - needed for Voigt function, H(a,v):
#double ii;
#// lineProf[0][0] = 0.0; v[0] = 0.0; //Line centre - cannot do logaritmically!
#double delta, core, logDelta;
#//int il0 = 36;
#//System.out.println("il0 " + il0 + " temp[il] " + temp[0][il0] + " press[il] " + logE*press[1][il0]);
for id in range(numDeps):
#//if (il <= numCore) {
#// - Gaussian ONLY - at line centre Lorentzian will diverge!
delta = 1.0
#//System.out.println("LINEGRID- CORE: core: " + core);
#//System.out.println("LINEGRID: il, v[il]: " + il + " " + v[il] + " lineProf[0][il]: " + lineProf[0][il]);
#//System.out.println("LINEGRID: il, Voigt, H(): " + il + " " + voigt);
#//Convert from H(a,v) in dimensionless Voigt units to physical phi((Delta lambda) profile:
logDelta = math.log(delta) + 2.0 * logLam0 - lnSqRtPi - logDopp - logC
lineProf[0][id] = math.exp(logDelta)
#//if (id == 36) {
#// System.out.println("il " + il + " linePoints " + 1.0e7 * linePoints[0][il] + " id " + id + " lineProf[il][id] " + lineProf[il][id]);
#//}
#//System.out.println("LineProf: il, id, lineProf[il][id]: " + il + " " + id + " " + lineProf[il][id]);
#// if (id == 20) {
#// for (int il = 0; il < numPoints; il++) {
#// System.out.format("Voigt: %20.16f %20.16f%n", linePoints[1][il], logE * Math.log(lineProf[il][id]));
#// }
#// }
#} //id loop
return lineProf
#} //end method delta()
def gauss(linePoints, lam0In, numDeps, teff, tauRos, temp, tempSun):
c = Useful.c()
logC = Useful.logC()
#//double k = Useful.k;
logK = Useful.logK()
#//double e = Useful.e;
#//double mE = Useful.mE;
lam0 = lam0In #// * 1.0E-7; //nm to cm
logLam0 = math.log(lam0)
ln10 = math.log(10.0)
ln2 = math.log(2.0)
ln4pi = math.log(4.0 * math.pi)
lnSqRtPi = 0.5 * math.log(math.pi)
sqPi = math.sqrt(math.pi)
#//double ln100 = 2.0*Math.log(10.0);
logE = math.log10(math.e) #// for debug output
doppler = linePoints[0][1] / linePoints[1][1]
logDopp = math.log(doppler)
tiny = 1.0e-19 #//??
#//System.out.println("LineProf: doppler, logDopp: " + doppler + " " + logE*logDopp);
#//Put input parameters into linear cgs units:
#//System.out.println("LINEGRID: Tau1: " + tau1);
#//logA = 2.0 * logLam0 + logGamma - ln4pi - logC - logDopp;
#//a = Math.exp(logA);
#//System.out.println("LINEGRID: logA: " + logE * logA);
#//Set up a half-profile Delta_lambda grid in Doppler width units
#//from line centre to wing
numPoints = len(linePoints[0])
#//System.out.println("LineProf: numPoints: " + numPoints);
#// Return a 2D numPoints X numDeps array of normalized line profile points (phi)
lineProf = [ [ 0.0 for i in range(numDeps) ] for j in range(numPoints) ]
#// Line profiel points in Doppler widths - needed for Voigt function, H(a,v):
v = [ 0.0 for i in range(numPoints) ]
#double logV, ii;
#// lineProf[0][0] = 0.0; v[0] = 0.0; //Line centre - cannot do logaritmically!
#double gauss, core, logGauss;
gauss = tiny #//default initialization
#//int il0 = 36;
#//System.out.println("il0 " + il0 + " temp[il] " + temp[0][il0] + " press[il] " + logE*press[1][il0]);
for id in range(numDeps):
for il in range(numPoints):
v[il] = linePoints[1][il]
#//System.out.println("LineProf: il, v[il]: " + il + " " + v[il]);
#//if (il <= numCore) {
if (v[il] <= 3.5 and v[il] >= -3.5):
#// - Gaussian ONLY - at line centre Lorentzian will diverge!
core = math.exp(-1.0 * (v[il] * v[il]))
gauss = core
#//System.out.println("LINEGRID- CORE: core: " + core);
#}
#//System.out.println("LINEGRID: il, v[il]: " + il + " " + v[il] + " lineProf[0][il]: " + lineProf[0][il]);
#//System.out.println("LINEGRID: il, Voigt, H(): " + il + " " + voigt);
#//Convert from H(a,v) in dimensionless Voigt units to physical phi((Delta lambda) profile:
logGauss = math.log(gauss) + 2.0 * logLam0 - lnSqRtPi - logDopp - logC
lineProf[il][id] = math.exp(logGauss)
#//if (id == 36) {
#// System.out.println("il " + il + " linePoints " + 1.0e7 * linePoints[0][il] + " id " + id + " lineProf[il][id] " + lineProf[il][id]);
#//}
#//System.out.println("LineProf: il, id, lineProf[il][id]: " + il + " " + id + " " + lineProf[il][id]);
#} // il lambda loop
#// if (id == 20) {
#// for (int il = 0; il < numPoints; il++) {
#// System.out.format("Voigt: %20.16f %20.16f%n", linePoints[1][il], logE * Math.log(lineProf[il][id]));
#// }
#// }
#} //id loop
""" /* Debug
// Check that line profile is area-normalized (it is NOT, area = 1.4845992503443734E-19!, but IS constant with depth - !?:
double delta;
for (int id = 0; id < numDeps; id++) {
double sum = 0.0;
for (int il = 1; il < numPoints2; il++) {
delta = lineProf2[0][il][id] - lineProf2[0][il - 1][id];
sum = sum + (lineProf2[1][il][id] * delta);
}
System.out.println("LineGrid: id, Profile area = " + id + " " + sum );
}
*/ """
return lineProf
#} //end method gauss()
def voigt(linePoints, lam0In, logAij, logGammaCol,
numDeps, teff, tauRos, temp, pGas,
tempSun, pGasSun, hjertComp, dbgHandle):
c = Useful.c()
logC = Useful.logC()
#//double k = Useful.k;
logK = Useful.logK()
#//double e = Useful.e;
#//double mE = Useful.mE;
lam0 = lam0In #// * 1.0E-7; //nm to cm
logLam0 = math.log(lam0)
ln10 = math.log(10.0)
ln2 = math.log(2.0);
ln4pi = math.log(4.0 * math.pi)
#lnSqRtPi = 0.5 * math.log(math.pi)
sqRtPi = math.sqrt(math.pi)
sqPi = math.sqrt(math.pi)
#//double ln100 = 2.0*Math.log(10.0);
logE = math.log10(math.e) #// for debug output
doppler = linePoints[0][1] / linePoints[1][1]
logDopp = math.log(doppler)
#//System.out.println("LineProf: doppler, logDopp: " + doppler + " " + logE*logDopp);
#//Put input parameters into linear cgs units:
#//double gammaCol = Math.pow(10.0, logGammaCol);
#// Lorentzian broadening:
#// Assumes Van der Waals dominates radiative damping
#// log_10 Gamma_6 for van der Waals damping around Tau_Cont = 1 in Sun
#// - p. 57 of Radiative Transfer in Stellar Atmospheres (Rutten)
logGammaSun = 9.0 * ln10 #// Convert to base e
#//double logFudge = Math.log(2.5); // Van der Waals enhancement factor
tau1 = ToolBox.tauPoint(numDeps, tauRos, 1.0)
#outline = ("tau1 "+ str(tau1) + "\n")
#dbgHandle.write(outline)
#//System.out.println("LINEGRID: Tau1: " + tau1);
#//logA = 2.0 * logLam0 + logGamma - ln4pi - logC - logDopp;
#//a = Math.exp(logA);
#//System.out.println("LINEGRID: logA: " + logE * logA);
#//Set up a half-profile Delta_lambda grid in Doppler width units
#//from line centre to wing
numPoints = len(linePoints[0])
#//System.out.println("LineProf: numPoints: " + numPoints);
#// Return a 2D numPoints X numDeps array of normalized line profile points (phi)
lineProf = [ [ 0.0 for i in range(numDeps) ] for j in range(numPoints) ]
#// Line profiel points in Doppler widths - needed for Voigt function, H(a,v):
v = [0.0 for i in range(numPoints)]
#double logV, ii;
#// lineProf[0][0] = 0.0; v[0] = 0.0; //Line centre - cannot do logaritmically!
#double gamma, logGamma, a, logA, voigt, core, wing, logWing, logVoigt;
Aij = math.pow(10.0, logAij)
il0 = 36
#// For Hjerting function approximation:
#double vSquare, vFourth, vAbs, a2, a3, a4, Hjert0, Hjert1, Hjert2, Hjert3, Hjert4, hjertFn;
#//System.out.println("il0 " + il0 + " temp[il] " + temp[0][il0] + " press[il] " + logE*press[1][il0]);
for id in range(numDeps):
#//Formula from p. 56 of Radiative Transfer in Stellar Atmospheres (Rutten),
#// logarithmically with respect to solar value:
logGamma = pGas[1][id] - pGasSun[1][tau1] + 0.7 * (tempSun[1][tau1] - temp[1][id]) + logGammaSun
#if (id%5 == 1):
# outline = ("id "+ str(id)+ " logGamma "+ str(logGamma) + "\n")
# dbgHandle.write(outline)
#//logGamma = logGamma + logFudge + logGammaCol;
logGamma = logGamma + logGammaCol
#//Add radiation (natural) broadning:
gamma = math.exp(logGamma) + Aij
logGamma = math.log(gamma)
#//
#//if (id == 12){
#//System.out.println("LineGrid: logGamma: " + id + " " + logE * logGamma + " press[1][id] " + press[1][id] + " pressSun[1][tau1] "
#// + pressSun[1][tau1] + " temp[1][id] " + temp[1][id] + " tempSun[1][tau1] " + tempSun[1][tau1]);
#// }
#//Voigt "a" parameter with line centre wavelength:
logA = 2.0 * logLam0 + logGamma - ln4pi - logC - logDopp
a = math.exp(logA)
a2 = math.exp(2.0*logA)
a3 = math.exp(3.0*logA)
a4 = math.exp(4.0*logA)
#// if (id == 12) {
#//System.out.println("LineGrid: lam0: " + lam0 + " logGam " + logE * logGamma + " logA " + logE * logA);
#// }
#//if (id == 30) {
#// //System.out.println("il v[il] iy y logNumerator logDenominator logInteg ");
#// System.out.println("voigt: v logVoigt: ");
#//}
for il in range(numPoints):
v[il] = linePoints[1][il]
vAbs = abs(v[il])
vSquare = vAbs * vAbs
vFourth = vSquare * vSquare
#//System.out.println("LineProf: il, v[il]: " + il + " " + v[il]);
#//Approximate Hjerting fn from tabulated expansion coefficients:
#// Interpolate in Hjerting table to exact "v" value for each expanstion coefficient:
#// Row 0 of Hjerting component table used for tabulated abscissae, Voigt "v" parameter
if (vAbs <= 12.0):
#//we are within abscissa domain of table
Hjert0 = ToolBox.interpol(hjertComp[0], hjertComp[1], vAbs)
Hjert1 = ToolBox.interpol(hjertComp[0], hjertComp[2], vAbs)
Hjert2 = ToolBox.interpol(hjertComp[0], hjertComp[3], vAbs)
Hjert3 = ToolBox.interpol(hjertComp[0], hjertComp[4], vAbs)
Hjert4 = ToolBox.interpol(hjertComp[0], hjertComp[5], vAbs)
else:
#// We use the analytic expansion
Hjert0 = 0.0
Hjert1 = (0.56419 / vSquare) + (0.846 / vFourth)
Hjert2 = 0.0
Hjert3 = -0.56 / vFourth
Hjert4 = 0.0
#//Approximate Hjerting fn with power expansion in Voigt "a" parameter
#// "Observation & Analysis of Stellar Photospeheres" (D. Gray), 3rd Ed., p. 258:
hjertFn = Hjert0 + a*Hjert1 + a2*Hjert2 + a3*Hjert3 + a4*Hjert4
#if ((id%5 == 1) and (il%2 == 0)):
# outline = ("il "+ str(il)+ " hjertFn "+ str(hjertFn) + "\n")
# dbgHandle.write(outline)
"""/* Gaussian + Lorentzian approximation:
//if (il <= numCore) {
if (v[il] <= 2.0 && v[il] >= -2.0) {
// - Gaussian ONLY - at line centre Lorentzian will diverge!
core = Math.exp(-1.0 * (v[il] * v[il]));
voigt = core;
//System.out.println("LINEGRID- CORE: core: " + core);
} else {
logV = Math.log(Math.abs(v[il]));
//Gaussian core:
core = Math.exp(-1.0 * (v[il] * v[il]));
// if (id == 12) {
// System.out.println("LINEGRID- WING: core: " + core);
// }
//Lorentzian wing:
logWing = logA - lnSqRtPi - (2.0 * logV);
wing = Math.exp(logWing);
voigt = core + wing;
// if (id == 12) {
// System.out.println("LINEGRID- WING: wing: " + wing + " logV " + logV);
// }
} // end else
*/"""
#//System.out.println("LINEGRID: il, v[il]: " + il + " " + v[il] + " lineProf[0][il]: " + lineProf[0][il]);
#//System.out.println("LINEGRID: il, Voigt, H(): " + il + " " + voigt);
#//Convert from H(a,v) in dimensionless Voigt units to physical phi((Delta lambda) profile:
#//logVoigt = Math.log(voigt) + 2.0 * logLam0 - lnSqRtPi - logDopp - logC;
#//System.out.println("voigt: Before log... id " + id + " il " + il + " hjertFn " + hjertFn);
#logVoigt = math.log(hjertFn) + 2.0 * logLam0 - lnSqRtPi - logDopp - logC
voigt = hjertFn * math.pow(lam0, 2) / sqRtPi / doppler / c
#logVoigt = math.log(voigt)
#lineProf[il][id] = math.exp(logVoigt)
lineProf[il][id] = voigt
if (lineProf[il][id] <= 0.0):
lineProf[il][id] = 1.0e-49
#// if (id == 12) {
#// System.out.println("il " + il + " linePoints " + 1.0e7 * linePoints[0][il] + " id " + id + " lineProf[il][id] " + lineProf[il][id]);
#// }
#//System.out.println("LineProf: il, id, lineProf[il][id]: " + il + " " + id + " " + lineProf[il][id]);
#} // il lambda loop
#// if (id == 20) {
#// for (int il = 0; il < numPoints; il++) {
#// System.out.format("Voigt: %20.16f %20.16f%n", linePoints[1][il], logE * Math.log(lineProf[il][id]));
#// }
#// }
#} //id loop
return lineProf
#} //end method voigt()
def stark(linePoints, lam0In, logAij, logGammaCol,
numDeps, teff, tauRos, temp, pGas, Ne,
tempSun, pGasSun, hjertComp, lineName):
c = Useful.c()
logC = Useful.logC()
#//double k = Useful.k;
logK = Useful.logK()
#//double e = Useful.e;
#//double mE = Useful.mE;
lam0 = lam0In #// * 1.0E-7; //nm to cm
logLam0 = math.log(lam0)
logLam0A = math.log(lam0) + 8.0*math.log(10.0) #//cm to A
ln10 = math.log(10.0)
ln2 = math.log(2.0)
ln4pi = math.log(4.0 * math.pi)
lnSqRtPi = 0.5 * math.log(math.pi)
sqRtPi = math.sqrt(math.pi)
sqPi = math.sqrt(math.pi)
#//double ln100 = 2.0*Math.log(10.0);
logE = math.log10(math.e) #// for debug output
doppler = linePoints[0][1] / linePoints[1][1]
logDopp = math.log(doppler)
#//System.out.println("LineProf: doppler, logDopp: " + doppler + " " + logE*logDopp);
#//Put input parameters into linear cgs units:
#//double gammaCol = Math.pow(10.0, logGammaCol);
#// Lorentzian broadening:
#// Assumes Van der Waals dominates radiative damping
#// log_10 Gamma_6 for van der Waals damping around Tau_Cont = 1 in Sun
#// - p. 57 of Radiative Transfer in Stellar Atmospheres (Rutten)
logGammaSun = 9.0 * ln10 #// Convert to base e
#//double logFudge = Math.log(2.5); // Van der Waals enhancement factor
tau1 = ToolBox.tauPoint(numDeps, tauRos, 1.0)
#//System.out.println("LINEGRID: Tau1: " + tau1);
#//logA = 2.0 * logLam0 + logGamma - ln4pi - logC - logDopp;
#//a = Math.exp(logA);
#//System.out.println("LINEGRID: logA: " + logE * logA);
#//Set up a half-profile Delta_lambda grid in Doppler width units
#//from line centre to wing
numPoints = len(linePoints[0])
#//System.out.println("LineProf: numPoints: " + numPoints);
#// Return a 2D numPoints X numDeps array of normalized line profile points (phi)
lineProf = [ [ 0.0 for i in range(numDeps) ] for j in range(numPoints) ]
#// Line profiel points in Doppler widths - needed for Voigt function, H(a,v):
v = [0.0 for i in range(numPoints)]
#double logV, ii;
#// lineProf[0][0] = 0.0; v[0] = 0.0; //Line centre - cannot do logaritmically!
#double gamma, logGamma, a, logA, voigt, core, wing, logWing, logVoigt;
Aij = math.pow(10.0, logAij)
il0 = 36
#// For Hjerting function approximation:
#double vSquare, vFourth, vAbs, a2, a3, a4, Hjert0, Hjert1, Hjert2, Hjert3, Hjert4, hjertFn;
#//Parameters for linear Stark broadening:
#//Assymptotic ("far wing") "K" parameters
#//Stehle & Hutcheon, 1999, A&A Supp Ser, 140, 93 and CDS data table
#//Assume K has something to do with "S" and proceed as in Observation and Analysis of
#// Stellar Photosphere, 3rd Ed. (D. Gray), Eq. 11.50,
#//
logTuneStark = math.log(3.16e7) #//convert DeltaI K parameters to deltaS STark profile parameters
logKStark = [0.0 for i in range(11)]
logKStark[0] = math.log(2.56e-03) + logTuneStark #//Halpha
logKStark[1] = math.log(7.06e-03) + logTuneStark #//Hbeta
logKStark[2] = math.log(1.19e-02) + logTuneStark #//Hgamma
logKStark[3] = math.log(1.94e-02) + logTuneStark #//Hdelta
logKStark[4] = math.log(2.95e-02) + logTuneStark #//Hepsilon
logKStark[5] = math.log(4.62e-02) + logTuneStark #//H8 JB
logKStark[6] = math.log(6.38e-02) + logTuneStark #//H9 JB
logKStark[7] = math.log(8.52e-02) + logTuneStark #//H10 JB
logKStark[8] = math.log(1.12e-01) + logTuneStark #//H11 JB
logKStark[9] = math.log(1.43e-01) + logTuneStark #//H12 JB
logKStark[10] = math.log(1.80e-01) + logTuneStark #//H13 JB
#//logKStark[11] = Math.log(2.11) + logTuneStark; //H30 JB
thisLogK = [0.0 for i in range(4)] #//default initialization
#//double thisLogK = logKStark[10]; //default initialization
#//which Balmer line are we? crude but effective:
if (lam0In > 650.0e-7):
thisLogK = logKStark[0] #//Halpha
#//System.out.println("Halpha")
#}
if ( (lam0In > 480.0e-7) and (lam0In < 650.0e-7) ):
#//System.out.println("Hbeta");
thisLogK = logKStark[1] #//Hbeta
#}
if ( (lam0In > 420.0e-7) and (lam0In < 470.0e-7) ):
#//System.out.println("Hgamma");
thisLogK = logKStark[2] #//Hgamma
#}
if ( (lam0In > 400.0e-7) and (lam0In < 450.0e-7) ):
#//System.out.println("Hdelta");
thisLogK = logKStark[3] #//Hdelta
if ( (lam0In < 400.0e-7) ):
#//System.out.println("Hepsilon");
thisLogK = logKStark[4] #//Hepsilon
#}
#// if ((lam0In < 390.0e-7)){
#//
#////This won't work here - "species" is always just "HI":
#// int numberInName = (int) lineName.substring("HI".length());
#// //console.log(numberInName);
#// thisLogK = logKStark[numberInName-3];
#// }
#//
#double F0, logF0, lamOverF0, logLamOverF0; //electrostatic field strength (e.s.u.)
#double deltaAlpha, logDeltaAlpha, logStark, logStarkTerm; //reduced wavelength de-tuning parameter (Angstroms/e.s.u.)
logF0Fac = math.log(1.249e-9)
#// log wavelength de-tunings in A:
#double logThisPoint, thisPoint;
#//System.out.println("il0 " + il0 + " temp[il] " + temp[0][il0] + " press[il] " + logE*press[1][il0]);
for id in range(numDeps):
#//linear Stark broadening stuff:
logF0 = logF0Fac + (0.666667)*Ne[1][id]
logLamOverF0 = logLam0A - logF0
lamOverF0 = math.exp(logLamOverF0)
#//System.out.println("id " + id + " logF0 " + logE*logF0 + " logLamOverF0 " + logE*logLamOverF0 + " lamOverF0 " + lamOverF0);
#//Formula from p. 56 of Radiative Transfer in Stellar Atmospheres (Rutten),
#// logarithmically with respect to solar value:
logGamma = pGas[1][id] - pGasSun[1][tau1] + 0.7 * (tempSun[1][tau1] - temp[1][id]) + logGammaSun
#//logGamma = logGamma + logFudge + logGammaCol
logGamma = logGamma + logGammaCol
#//Add radiation (natural) broadning:
gamma = math.exp(logGamma) + Aij
logGamma = math.log(gamma)
#//
#//if (id == 12){
#//System.out.println("LineGrid: logGamma: " + id + " " + logE * logGamma + " press[1][id] " + press[1][id] + " pressSun[1][tau1] "
#// + pressSun[1][tau1] + " temp[1][id] " + temp[1][id] + " tempSun[1][tau1] " + tempSun[1][tau1]);
#// }
#//Voigt "a" parameter with line centre wavelength:
logA = 2.0 * logLam0 + logGamma - ln4pi - logC - logDopp
a = math.exp(logA)
a2 = math.exp(2.0*logA)
a3 = math.exp(3.0*logA)
a4 = math.exp(4.0*logA)
#// if (id == 12) {
#//System.out.println("LineGrid: lam0: " + lam0 + " logGam " + logE * logGamma + " logA " + logE * logA);
#// }
#//if (id == 30) {
#// //System.out.println("il v[il] iy y logNumerator logDenominator logInteg ");
#// System.out.println("voigt: v logVoigt: ");
#//}
for il in range(numPoints):
v[il] = linePoints[1][il]
vAbs = abs(v[il])
vSquare = vAbs * vAbs
vFourth = vSquare * vSquare
#//System.out.println("LineProf: il, v[il]: " + il + " " + v[il]);
#//Approximate Hjerting fn from tabulated expansion coefficients:
#// Interpolate in Hjerting table to exact "v" value for each expanstion coefficient:
#// Row 0 of Hjerting component table used for tabulated abscissae, Voigt "v" parameter
if (vAbs <= 12.0):
#//we are within abscissa domain of table
Hjert0 = ToolBox.interpol(hjertComp[0], hjertComp[1], vAbs)
Hjert1 = ToolBox.interpol(hjertComp[0], hjertComp[2], vAbs)
Hjert2 = ToolBox.interpol(hjertComp[0], hjertComp[3], vAbs)
Hjert3 = ToolBox.interpol(hjertComp[0], hjertComp[4], vAbs)
Hjert4 = ToolBox.interpol(hjertComp[0], hjertComp[5], vAbs)
else:
#// We use the analytic expansion
Hjert0 = 0.0
Hjert1 = (0.56419 / vSquare) + (0.846 / vFourth)
Hjert2 = 0.0
Hjert3 = -0.56 / vFourth
Hjert4 = 0.0
#}
#//Approximate Hjerting fn with power expansion in Voigt "a" parameter
#// "Observation & Analysis of Stellar Photospeheres" (D. Gray), 3rd Ed., p. 258:
hjertFn = Hjert0 + a*Hjert1 + a2*Hjert2 + a3*Hjert3 + a4*Hjert4;
logStark = -49.0 #//re-initialize
if (vAbs > 2.0):
#//System.out.println("Adding in Stark wing");
thisPoint = 1.0e8 * abs(linePoints[0][il]) #//cm to A
logThisPoint = math.log(thisPoint)
logDeltaAlpha = logThisPoint - logF0
deltaAlpha = math.exp(logDeltaAlpha)
logStarkTerm = ( math.log(lamOverF0 + deltaAlpha) - logLamOverF0 )
logStark = thisLogK + 0.5*logStarkTerm - 2.5*logDeltaAlpha
#//System.out.println("il " + il + " logDeltaAlpha " + logE*logDeltaAlpha + " logStarkTerm " + logE*logStarkTerm + " logStark " + logE*logStark);
#//console.log("il " + il + " logDeltaAlpha " + logE*logDeltaAlpha + " logStarkTerm " + logE*logStarkTerm + " logStark " + logE*logStark);
#//System.out.println("id " + id + " il " + il + " v[il] " + v[il]
#// + " hjertFn " + hjertFn + " Math.exp(logStark) " + Math.exp(logStark));
#//not here! hjertFn = hjertFn + Math.exp(logStark);
#//System.out.println("LINEGRID: il, v[il]: " + il + " " + v[il] + " lineProf[0][il]: " + lineProf[0][il]);
#//System.out.println("LINEGRID: il, Voigt, H(): " + il + " " + voigt);
#//Convert from H(a,v) in dimensionless Voigt units to physical phi((Delta lambda) profile:
#//logVoigt = Math.log(voigt) + 2.0 * logLam0 - lnSqRtPi - logDopp - logC;
#//System.out.println("stark: Before log... id " + id + " il " + il + " hjertFn " + hjertFn);
#logVoigt = math.log(hjertFn) - lnSqRtPi - logDopp
voigt = hjertFn / sqRtPi / doppler
#//logVoigt = math.log(voigt)
logStark = logStark - logF0
if (vAbs > 2.0):
#//if (id == 24) {
#// System.out.println("il " + il + " v[il] " + v[il] + " logVoigt " + logE*logVoigt + " logStark " + logE*logStark);
#//}
#//voigt = math.exp(logVoigt) + math.exp(logStark)
voigt = voigt + math.exp(logStark)
#//logVoigt = math.log(voigt)
#logVoigt = logVoigt + 2.0 * logLam0 - logC
voigt = voigt * math.pow(lam0, 2) / c
#//lineProf[il][id] = math.exp(logVoigt)
lineProf[il][id] = voigt
if (lineProf[il][id] <= 0.0):
lineProf[il][id] = 1.0e-49
#//if (id == 24) {
#// System.out.println("lam0In " + lam0In);
#// System.out.println("il " + il + " linePoints " + 1.0e7 * linePoints[0][il] + " id " + id + " lineProf[il][id] " + lineProf[il][id]);
#//}
#//System.out.println("LineProf: il, id, lineProf[il][id]: " + il + " " + id + " " + lineProf[il][id]);
#} // il lambda loop
#// if (id == 20) {
#// for (int il = 0; il < numPoints; il++) {
#// System.out.format("Voigt: %20.16f %20.16f%n", linePoints[1][il], logE * Math.log(lineProf[il][id]));
#// }
#// }
#} //id loop
return lineProf
#} //end method stark()
def lineSource(numDeps, tau, temp, lambda2):
"""#// Make line source function:
#// Equivalenth two-level atom (ETLA) approx
#//CAUTION: input lambda in nm"""
lineSource = [0.0 for i in range(numDeps)]
#//thermal photon creation/destruction probability
epsilon = 0.01 #//should decrease with depth??
#//This is an artifact of jayBinner's original purpose:
grayLevel = 1.0
#//int iLam0 = numLams / 2; //+/- 1 deltaLambda
#//double lam0 = linePoints[0][iLam0]; //line centre lambda in cm - not needed:
#//double lamStart = lambda - 0.1; // nm
#//double lamStop = lambda + 0.1; // nm
#//double lamRange = (lamStop - lamStart); // * 1.0e-7; // line width in cm
#//System.out.println("lamStart " + lamStart + " lamStop " + lamStop + " lamRange " + lamRange);
jayLambda = [0.0 for i in range(numDeps)]
BLambda = [ [ 0.0 for i in range(numDeps) ] for j in range(2) ]
#double linSrc;
#// Dress up Blambda to look like what jayBinner expects:
for i in range(numDeps):
#//Planck.planck return log(B_lambda):
BLambda[0][i] = math.exp(Planck.planck(temp[0][i], lambda2))
BLambda[1][i] = 1.0 #//supposed to be dB/dT, but not needed.
#//CAUTION: planckBin Row 0 is linear lambda-integrated B_lambda; Row 1 is same for dB_lambda/dT
#//planckBin = MulGrayTCorr.planckBinner(numDeps, temp, lamStart, lamStop);
jayLambda = MulGrayTCorr.jayBinner(numDeps, tau, temp, BLambda, grayLevel)
#//To begin with, coherent scattering - we're not computing line profile-weighted average Js and Bs
for i in range(numDeps):
#//planckBin[0][i] = planckBin[0][i] / lamRange; //line average
#//jayBin[i] = jayBin[i];
linSrc = (1.0 - epsilon) * jayLambda[i] + epsilon * BLambda[0][i]
lineSource[i] = math.log(linSrc)
return lineSource
