Raw File
# -*- 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

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
    #//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("SpecSyn: numNow ", numNow, " numPoints ", numPoints)
    #print("SpecSyn: len(thisMasterLams) ", len(thisMasterLams), " len(logKappa1D) ", len(logKappa1D))
    #print("SpecSyn: len(listLineLambdas) ", len(listLineLambdas), " len(logLineKap1D) ", len(logLineKap1D))
    for k in range(numNow):
        thisMasterLams[k] = masterLams[k]
        
    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 

        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]);
            #//     }
            

        #//Interpolate continuum and line opacity onto master lambda scale, and add them lambda-wise:
        for iL in range(numTot):
            #kappa2 = ToolBox.interpol(masterLams, kappa1D, masterLamsOut[iL])
            #logKappa2 = ToolBox.interpol(masterLams, logKappa1D, masterLamsOut[iL])
            logKappa2 = ToolBox.interpol(thisMasterLams, logKappa1D, masterLamsOut[iL])
            #lineKap2 = 1.0e-49 #//re-initialization
            logLineKap2 = -49.0 #//re-initialization
            if ( (masterLamsOut[iL] >= listLineLambdas[0]) and (masterLamsOut[iL] <= listLineLambdas[numPoints-1]) ):
                
                #lineKap2 = ToolBox.interpol(listLineLambdas, lineKap1D, masterLamsOut[iL])
                logLineKap2 = ToolBox.interpol(listLineLambdas, logLineKap1D, masterLamsOut[iL])
                #if (lineKap2 <= 0.0):
                #    lineKap2 = 1.0e-49
                #//lineKap2 = 1.0e-99;  //test
                
            #//test lineKap2 = 1.0e-99;  //test
            #// if (iD%10 == 1){
            #//   System.out.println("iD " + iD + " iL " + iL + " masterLamsOut " + masterLamsOut[iL] + " kappa2 " + kappa2 + " lineKap2 " + lineKap2);
            #//}
            #totKap = kappa2 + lineKap2
            totKap = math.exp(logKappa2) + math.exp(logLineKap2)
            logMasterKapsOut[iL][iD] = math.log(totKap)
            #//if (iD == 36) {
            #//    System.out.format("%02d   %12.8e   %12.8f%n", iL, masterLams[iL], logE * logMasterKappa[iL][iD]);
            #//}
        #} iL loop    
    #}  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