https://github.com/sevenian3/ChromaStarPy
Raw File
Tip revision: 103d3d0df6d9574c49818f149e1cae8100455d10 authored by Ian Short on 06 July 2023, 18:09:20 UTC
Create ReadMe
Tip revision: 103d3d0
SpecSyn2.py
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 28 17:03:30 2017

@author: ishort
"""

#/**
# *
# * Create master kappa_lambda(lambda) and tau_lambda(lambda) for
# * FormalSoln.formalSoln()
# *
# * @author Ian
# */

import math
import ToolBox
#plotting:
import matplotlib
import pylab
import numpy

def masterLambda(numLams, numMaster, numNow, masterLams, numPoints, listLineLambdas):
        
    """//Merge continuum and line wavelength scales - for one line
//This expects *pure* line opacity - no continuum opacity pre-added!"""
    
    #//int numCnt = lambdaScale.length;
    #//skip the last wavelength point in the line lambda grid - it holds the line centre wavelength
    #//int numLine = lineLambdas.length - 1;

    numTot = numNow + numPoints #//current dynamic total

    #//System.out.println("numCnt " + numCnt + " numLine " + numLine + " numTot " + numTot);
    """/*
         for (int i = 0; i < numCnt; i++) {
         System.out.println("i " + i + " lambdaScale[i] " + lambdaScale[i]);
         }
         for (int i = 0; i < numLine; i++) {
         System.out.println("i " + i + " lineLambdas[i] " + lineLambdas[i]);
         }
         */ """
    #//Row 0 is merged lambda scale
    #//Row 1 is log of *total* (line plus continuum kappa
    masterLamsOut = [0.0 for i in range(numTot)]

    #// Merge wavelengths into a sorted master list
    #//initialize with first continuum lambda:
    lastLam = masterLams[0]
    masterLamsOut[0] = masterLams[0]
    nextCntPtr = 1
    nextLinePtr = 0
    for iL in range(1, numTot):
        if (nextCntPtr < numNow):
            #//System.out.println("nextCntPtr " + nextCntPtr + " lambdaScale[nextCntPtr] " + lambdaScale[nextCntPtr]);
            #//System.out.println("nextLinePtr " + nextLinePtr + " lineLambdas[nextLinePtr] " + lineLambdas[nextLinePtr]);
            if ((masterLams[nextCntPtr] <= listLineLambdas[nextLinePtr])
            or (nextLinePtr >= numPoints - 1)):
                #//Next point is a continuum point:
                masterLamsOut[iL] = masterLams[nextCntPtr]
                nextCntPtr+=1

            elif ((listLineLambdas[nextLinePtr] < masterLams[nextCntPtr])
                        and (nextLinePtr < numPoints - 1)):
                #//Next point is a line point:
                masterLamsOut[iL] = listLineLambdas[nextLinePtr]
                nextLinePtr+=1            
            
            #//System.out.println("iL " + iL + " masterLamsOut[iL] " + masterLamsOut[iL]);
    #} //iL loop
    #//Make sure final wavelength point in masterLams is secured:
    masterLamsOut[numTot-1] = masterLams[numNow-1]

    return masterLamsOut
#}

def masterKappa(numDeps, numLams, numMaster, numNow, masterLams, masterLamsOut, logMasterKaps, \
                numPoints, listLineLambdas, listLogKappaL):
#//                                          
    logE = math.log10(math.e) #// for debug output

    #//int numLams = masterLams.length;
    numTot = numNow + numPoints
        
    logMasterKapsOut = [ [ 0.0 for i in range(numDeps) ] for j in range(numTot) ]
    #//double[][] kappa2 = new double[2][numTot];
    #//double[][] lineKap2 = new double[2][numTot];
    #double kappa2, lineKap2, totKap;
    #lineKap2 = 1.0e-99 #//initialization
    #logLineKap2 = -49.0 #//initialization
    logKappa2 = [0.0 for i in range(numTot)]
    logLineKap2 = [-49.0 for i in range(numTot)]
    #//int numCnt = lambdaScale.length;
    #//int numLine = lineLambdas.length - 1;
    #kappa1D = [0.0 for i in range(numNow)]
    logKappa1D = [0.0 for i in range(numNow)]
    #thisMasterLams = [0.0 for i in range(numNow)]
    #lineKap1D = [0.0 for i in range(numPoints)]
    logLineKap1D = [0.0 for i in range(numPoints)]
    #//System.out.println("iL   masterLams    logMasterKappa");
    #print("numNow ", numNow, " numPoints ", numPoints)
    #print("iD ", iD, " len(masterLams) ", len(masterLams), " len(logKappa1D) ", len(logKappa1D))           
    
    #for k in range(numNow):
    #    thisMasterLams[k] = masterLams[k]
    thisMasterLams = [ masterLams[k] for k in range(numNow) ]
        
    for iD in range(numDeps):

        #//Extract 1D *linear* opacity vectors for interpol()
        #for k in range(numNow):
        #    #kappa1D[k] = math.exp(logMasterKaps[k][iD]) #//now wavelength dependent 
        #    logKappa1D[k] = logMasterKaps[k][iD] #//now wavelength dependent
        logKappa1D = [ logMasterKaps[k][iD] for k in range(numNow) ]

        #for k in range(numPoints):
        #    #lineKap1D[k] = math.exp(listLogKappaL[k][iD])
        #    logLineKap1D[k] = listLogKappaL[k][iD]
            #//     if (iD%10 == 1){
            #//        System.out.println("iD " + iD + " k " + k + " listLineLambdas " + listLineLambdas[k] + " lineKap1D " + lineKap1D[k]);
            #//     }
        logLineKap1D = [ listLogKappaL[k][iD] for k in range(numPoints) ]    

        #//Interpolate continuum and line opacity onto master lambda scale, and add them lambda-wise:
        #for iL in range(numTot):
        #    logLineKap2[iL] = -49.0 #//re-initialization
        logLineKap2 = [ -49.0 for iL in range(numTot) ]
 
        logKappa2 = numpy.interp(masterLamsOut, thisMasterLams, logKappa1D)
        logLineKap2 = numpy.interp(masterLamsOut, listLineLambdas, logLineKap1D)
        for iL in range(numTot):
            totKap = math.exp(logKappa2[iL]) + math.exp(logLineKap2[iL])
            logMasterKapsOut[iL][iD] = math.log(totKap)
        #logMasterKapsOut[:][iD] =\
        #[ math.log(math.exp(logKappa2[iL]) + math.exp(logLineKap2[iL])) for iL in range(numTot) ]
        
    #}  iD loop 
    #pylab.plot(masterLamsOut, [logMasterKaps[i][12] for i in range(numTot)]) 
    #pylab.plot(masterLamsOut, [logMasterKaps[i][12] for i in range(numTot)], '.')      

    return logMasterKapsOut;
#}
    
back to top