# -*- coding: utf-8 -*-
"""
Created on Fri Apr 21 17:28:32 2017

Initializes and re-scales a Phoenix LTE spherical reference model of 
 * Teff=10000K, log(g)=4.0, [Fe/H]=0.0, xi=1.0 km/s, l=1.0H_p, M=1M_Sun, R=1.1613E+11cm

@author: ishort
"""

import math
import ToolBox
import Useful

def phxRefTeff():
    return 10000.0

def phxRefLogEg(): 
    return math.log(10.0) * 4.0  #//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 = [
    6.07574016685149309E+03, 6.07574016685149309E+03, 6.13264671606194861E+03,
    6.20030362747541585E+03, 6.27534705504544127E+03, 6.35396254937768026E+03,
    6.43299900128272293E+03, 6.51018808525609893E+03, 6.58411555606889124E+03,
    6.65406717610081068E+03, 6.71983498258185136E+03, 6.78154367852633823E+03,
    6.83954193198123903E+03, 6.89437231818902364E+03, 6.94676889243451842E+03,
    6.99759489202792247E+03, 7.04769490055547158E+03, 7.09773520027041195E+03,
    7.14812062339764907E+03, 7.19901426577775601E+03, 7.25041827414427917E+03,
    7.30225171801659872E+03, 7.35440093819652611E+03, 7.40675066225539558E+03,
    7.45920456139609178E+03, 7.51166464185182758E+03, 7.56404228766520191E+03,
    7.61627005664532771E+03, 7.66833575187113820E+03, 7.72034173334201841E+03,
    7.77258785750414881E+03, 7.82555139374063583E+03, 7.87986936059489017E+03,
    7.93639246968124371E+03, 7.99620846303960116E+03, 8.06052820253916161E+03,
    8.13047124123426238E+03, 8.20741189262034641E+03, 8.29307358429898159E+03,
    8.38980788216330802E+03, 8.49906053657168923E+03, 8.62314483632361771E+03,
    8.76456384216990409E+03, 8.92693370905029224E+03, 9.11177170396923248E+03,
    9.32167977041711492E+03, 9.56236981551314602E+03, 9.82432656703466455E+03,
    1.01311427939962559E+04, 1.04299661074183350E+04, 1.08355089220389909E+04,
    1.12094886773674716E+04, 1.16360710406256258E+04, 1.20991237739366334E+04,
    1.25891111265208237E+04, 1.31070008299570563E+04, 1.36522498965801387E+04,
    1.42233473670298790E+04, 1.48188302103200131E+04, 1.54423659243804523E+04,
    1.60892587452310745E+04, 1.67828517694842230E+04, 1.74930217234773954E+04,
    1.82922661949382236E+04]

    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):

    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, 8.32127743125684882E-02, 1.29584527404206007E-01,
    1.94435381478779895E-01, 2.81524759872055830E-01, 3.94850766488002047E-01,
    5.39098197994885120E-01, 7.20109114447812781E-01, 9.45331395103965466E-01,
    1.22424260721948497E+00, 1.56877812718506826E+00, 1.99379948180689048E+00,
    2.51761637911653136E+00, 3.16251087800302599E+00, 3.95513966878889667E+00,
    4.92671520309637767E+00, 6.11303768406991388E+00, 7.55464145673977505E+00,
    9.29736005428628154E+00, 1.13934670418806956E+01, 1.39033471883101818E+01,
    1.68975909460311797E+01, 2.04594801623940121E+01, 2.46878880919212804E+01,
    2.97005964646718432E+01, 3.56383534114781142E+01, 4.26698208468708984E+01,
    5.09974403334007107E+01, 6.08640463074419387E+01, 7.25594340179816726E+01,
    8.64248329112294158E+01, 1.02854593091977520E+02, 1.22294652156180661E+02,
    1.45234045163109670E+02, 1.72184927273123520E+02, 2.03652334583264832E+02,
    2.40105656346438934E+02, 2.81936164286554344E+02, 3.29393094590693863E+02,
    3.82482413201705356E+02, 4.40963324580460835E+02, 5.04333229685725428E+02,
    5.71827998329611432E+02, 6.42424030136117835E+02, 7.15115448265608620E+02,
    7.89188190751975185E+02, 8.64179477829598227E+02, 9.41037808653716070E+02,
    1.02093026109089942E+03, 1.10808816566702853E+03, 1.20591338801728261E+03,
    1.32157321934523725E+03, 1.46400967396971282E+03, 1.64395527893530380E+03,
    1.87431044562489683E+03, 2.16986659968736876E+03, 2.54753164223200429E+03,
    3.02667796755900645E+03, 3.62964225373483487E+03, 4.38288420138537458E+03,
    5.31730879832813844E+03, 6.47251190142057658E+03, 7.89413608165941059E+03,
    9.64747840003540659E+03]

    logPhxRefTau64 =getLogPhxRefTau64()

    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 ]
        

    #// interpolate onto gS3 tauRos grid and re-scale with Teff:
    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.53 #//top of model
    gexpBottom = 0.85 #//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):
        logPhxRefPGas[i] = ToolBox.interpol(logPhxRefTau64, logPhxRefPGas64, tauRos[1][i])
        thisGexp = gexpTop + gexpRange * (tauRos[1][i] - tauRos[1][0]) / tauLogRange
        #//scaling with g
        scalePGas[1][i] = thisGexp*logEg + logPhxRefPGas[i] - thisGexp*phxRefLogEg()
        #//scaling with zscl:
        #scalePGas[1][i] = -0.5*logZScale + scalePGas[1][i]
        ##//scaling with A_He:
        #scalePGas[1][i] = 0.666667 * math.log(1.0 + 4.0*AHe) + scalePGas[1][i] - logHeDenom 
        #scalePGas[0][i] = math.exp(scalePGas[1][i])
        #//System.out.println("scalePGas[1][i] " + logE * scalePGas[1][i])   
    scalePGas[1] = [ -0.5*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 = [
    4.77258390479251340E-05, 1.54333794509103339E-02, 2.24384775218179552E-02,
    3.24056217848841463E-02, 4.62639509784656192E-02, 6.49897301016105072E-02,
    8.96001972148401798E-02, 1.21161157265374353E-01, 1.60825358340301261E-01,
    2.09891146620685975E-01, 2.69867426146356171E-01, 3.42538888354808724E-01,
    4.30045384007358256E-01, 5.35006986797593842E-01, 6.60704782988379868E-01,
    8.11262305821688567E-01, 9.91741961224463009E-01, 1.20813527252446407E+00,
    1.46731521914247520E+00, 1.77705126262480850E+00, 2.14614122290851617E+00,
    2.58462667298359561E+00, 3.10405210627260297E+00, 3.71777653138435804E+00,
    4.44135288803457673E+00, 5.29279499891786465E+00, 6.29303772366266312E+00,
    7.46652989782078702E+00, 8.84221515332682451E+00, 1.04552216626003140E+01,
    1.23496848557054300E+01, 1.45813048229500843E+01, 1.72206436663779385E+01,
    2.03589457441922157E+01, 2.41156208954111868E+01, 2.86442876094033458E+01,
    3.41355927487861948E+01, 4.08398462152914732E+01, 4.90908766488638761E+01,
    5.93486059459067832E+01, 7.21405304518226842E+01, 8.81824952094146681E+01,
    1.08367129768339950E+02, 1.33856171619767082E+02, 1.65693080738235807E+02,
    2.04943252558813072E+02, 2.52705001145053956E+02, 3.07224623951268654E+02,
    3.70334137141753217E+02, 4.33722318385145741E+02, 5.08910395587106336E+02,
    5.82220694357564639E+02, 6.65278728107771599E+02, 7.62124991657425880E+02,
    8.79654481582760809E+02, 1.02622262715821921E+03, 1.21099204341081804E+03,
    1.44432886438589208E+03, 1.73838904022049860E+03, 2.10808802008476914E+03,
    2.57102379769462232E+03, 3.14976025581092108E+03, 3.86645770963505538E+03,
    4.75493678618616923E+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.53 #//top of model
    gexpBottom = 0.82 #//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):
        #logPhxRefPe[i] = ToolBox.interpol(logPhxRefTau64, logPhxRefPe64, tauRos[1][i])
        thisGexp = gexpTop + gexpRange * (tauRos[1][i] - tauRos[1][0]) / tauLogRange
        #//scaling with g
        scalePe[1][i] = thisGexp*logEg + logPhxRefPe[i] - thisGexp*phxRefLogEg()
        #//scale with Teff:
        if (teff < 10000.0):
            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
               
            scalePe[1][i] = thisOmega*teff + scalePe[1][i] - thisOmega*phxRefTeff()
            
        #//scaling with zscl:
        #scalePe[1][i] = 0.5*logZScale + scalePe[1][i]
        #//scaling with A_He:
        #scalePe[1][i] = 0.333333 * math.log(1.0 + 4.0*AHe) + scalePe[1][i] - logHeDenom
        #scalePe[1][i] = logEg + logPhxRefPe[i] - phxRefLogEg()
        #scalePe[0][i] = math.exp(scalePe[1][i])
    scalePe[1] = [ 0.5*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[1] = [ logEg + x - phxRefLogEg() for x in logPhxRefPe ]
    scalePe[0] = [ math.exp(x) for x in scalePe[1] ]
        #//System.out.println("scaleNe[1][i] " + logE * scaleNe[1][i]);
        
    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