# -*- coding: utf-8 -*-
"""
Created on Fri Apr 21 09:59:40 2017

Routines to read in a standard solar atmosphere model computed with Phoenix V15 and
calculate associated quantities.  Needed for physical treatments normalized with
solar qantities

@author: ishort
"""

import math
import ToolBox
import Useful

def phxSunTeff():
    return 5777.0

def phxSunLogEg():
    return math.log(10.0) * 4.44  #//base e!

def getPhxSunTau64():
    
    phxSunTau64 = [
    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.07710505603676912e-05, 1.44974067037263168e-05, 1.95129342263596224e-05,
    2.62636352765333536e-05, 3.53498110503010944e-05, 4.75794431400941376e-05,
    6.40400427119728256e-05, 8.61953566475303296e-05, 1.16015530173997152e-04,
    1.56152300600049664e-04, 2.10174801133248704e-04, 2.82886943462596928e-04,
    3.80754602122237184e-04, 5.12480587696093120e-04, 6.89778537938765824e-04,
    9.28414544519474432e-04, 1.24960914129198688e-03, 1.68192432488086880e-03,
    2.26380340952144672e-03, 3.04698957090350784e-03, 4.10112707055130048e-03,
    5.51995432128156800e-03, 7.42963950759494912e-03, 1.00000000000000002e-02,
    1.34596032415536416e-02, 1.81160919420041312e-02, 2.43835409826882656e-02,
    3.28192787251147072e-02, 4.41734470314006464e-02, 5.94557070854439424e-02,
    8.00250227816105216e-02, 1.07710505603676912e-01, 1.44974067037263136e-01,
    1.95129342263596224e-01, 2.62636352765332992e-01, 3.53498110503010240e-01,
    4.75794431400941440e-01, 6.40400427119728384e-01, 8.61953566475303168e-01,
    1.16015530173997152e+00, 1.56152300600049664e+00, 2.10174801133248704e+00,
    2.82886943462596640e+00, 3.80754602122236800e+00, 5.12480587696092608e+00,
    6.89778537938765824e+00, 9.28414544519474432e+00, 1.24960914129198672e+01,
    1.68192432488086880e+01, 2.26380340952144640e+01, 3.04698957090350528e+01,
    4.10112707055129536e+01, 5.51995432128157312e+01, 7.42963950759495040e+01,
    1.00000000000000000e+02]
    
    return phxSunTau64

def getLogPhxSunTau64():

    logE = math.log10(math.e)

    phxSunTau64 = getPhxSunTau64()
    numPhxDep = len(phxSunTau64)
    #print("numPhxDep ", numPhxDep)
    logPhxSunTau64 = [0.0 for i in range(numPhxDep)]
    for i in range(1, numPhxDep):
        #print("i ", i, " phxSunTau64[i] ", phxSunTau64[i])
        logPhxSunTau64[i] = math.log(phxSunTau64[i])
        
    logPhxSunTau64[0] = logPhxSunTau64[1] - (logPhxSunTau64[numPhxDep - 1] - logPhxSunTau64[1]) / numPhxDep
    return logPhxSunTau64


def phxSunTemp(teff, numDeps, tauRos): 

    logE = math.log10(math.e)

    #//Theoretical radiative/convective model from Phoenix V15:
    phxSunTemp64 = [
    3.75778887392339840e+03, 3.75778887392339840e+03, 3.78480175327941504e+03,
    3.81385432525541760e+03, 3.84360130602512768e+03, 3.87340585446516608e+03,
    3.90300184305606656e+03, 3.93231689265254528e+03, 3.96137919852984000e+03,
    3.99027119028325824e+03, 4.01910484194699648e+03, 4.04798292490651008e+03,
    4.07699548886169152e+03, 4.10623218035810816e+03, 4.13574364539801920e+03,
    4.16548101060783104e+03, 4.19541371831173824e+03, 4.22551121760088000e+03,
    4.25571229065970624e+03, 4.28594188575783232e+03, 4.31613168919769152e+03,
    4.34620698440244928e+03, 4.37603327507328960e+03, 4.40564394765877952e+03,
    4.43507740841559296e+03, 4.46439148496796224e+03, 4.49375530130093952e+03,
    4.52341166116436480e+03, 4.55357281866347264e+03, 4.58446079852491520e+03,
    4.61663974201107520e+03, 4.65052341797810624e+03, 4.68623381803595456e+03,
    4.72408924142126144e+03, 4.76494152329308416e+03, 4.80984310271200128e+03,
    4.85897778977827584e+03, 4.91315894280032960e+03, 4.97390461818851328e+03,
    5.04531167969494336e+03, 5.12680296183560704e+03, 5.22061204180252480e+03,
    5.32918534350649152e+03, 5.46202432323604352e+03, 5.61966782651567040e+03,
    5.80986721241013376e+03, 6.03911828822760320e+03, 6.23433005487621120e+03,
    6.53458311644527488e+03, 6.87429103746811904e+03, 7.29999981509928192e+03,
    7.66682942009826304e+03, 7.94223816217841024e+03, 8.16133659245977728e+03,
    8.35020013757955200e+03, 8.52047273964030720e+03, 8.67812135633704064e+03,
    8.82687568743616768e+03, 8.96926538519515648e+03, 9.10706359999037824e+03,
    9.24154121553023488e+03, 9.37363000902155008e+03, 9.50427569030960000e+03,
    9.63219702937432192e+03]

    #// interpolate onto gS3 tauRos grid and re-scale with Teff:
       
    phxSunTemp = [0.0 for i in range(numDeps)]
    scaleTemp = [ [0.0 for i in range(numDeps)] for j in range(2)]
    logPhxSunTau64 = getLogPhxSunTau64()
    #for i in range(64):
       #print("i ", i, " logPhxSunTau64 ", logPhxSunTau64[i], " phxSunTemp64 ", phxSunTemp64[i]) 
        
    #print("phxSunTemp: numDeps ", numDeps)
    for i in range(numDeps):
        phxSunTemp[i] = ToolBox.interpol(logPhxSunTau64, phxSunTemp64, tauRos[1][i])
        #print("i ", i, " tauRos[1][i] ", tauRos[1][i], " phxSunTemp[i] ", phxSunTemp[i])
        scaleTemp[0][i] = teff * phxSunTemp[i] / phxSunTeff()
        scaleTemp[1][i] = math.log(scaleTemp[0][i])
        #//System.out.println("tauRos[1][i] " + logE * tauRos[1][i] + " scaleTemp[1][i] " + logE * scaleTemp[1][i]);
        
    return scaleTemp
    
def phxSunPGas(grav, numDeps, tauRos): 

    logE = math.log10(math.e)
    logEg = math.log(grav) #//base e!

    #//Theoretical radiative/convective model from Phoenix V15:
    phxSunPGas64 = [
    1.00000000000000005e-04, 7.28828683006412544e+01, 8.61732126528505984e+01,
    1.01843641855932976e+02, 1.20317369304629504e+02, 1.42093296011949696e+02,
    1.67758727999644384e+02, 1.98004769223716256e+02, 2.33644726494082176e+02,
    2.75635953319757664e+02, 3.25104809120938880e+02, 3.83378880706399168e+02,
    4.52022443862726592e+02, 5.32877321364649344e+02, 6.28113128741022208e+02,
    7.40284569989930496e+02, 8.72399144145001216e+02, 1.02799724165148560e+03,
    1.21124571517496000e+03, 1.42704756928025427e+03, 1.68117132827309248e+03,
    1.98040330055171296e+03, 2.33272402439094176e+03, 2.74752260927171264e+03,
    3.23584954067544384e+03, 3.81071167175796544e+03, 4.48742481283848128e+03,
    5.28403135449994368e+03, 6.22178478543013120e+03, 7.32571484052561408e+03,
    8.62531818498740864e+03, 1.01553497268350960e+04, 1.19567104697253520e+04,
    1.40775115384306991e+04, 1.65743702896828832e+04, 1.95139034178162464e+04,
    2.29742653550211872e+04, 2.70468752817440448e+04, 3.18381441391788864e+04,
    3.74704748898233472e+04, 4.40799661582952512e+04, 5.18080650391892096e+04,
    6.07793647633492224e+04, 7.10351288049853440e+04, 8.24259773567987968e+04,
    9.44866985169806080e+04, 1.06329924298695632e+05, 1.17862219382348656e+05,
    1.28295128203359424e+05, 1.36933948396180352e+05, 1.43493910023715958e+05,
    1.48487688700034048e+05, 1.52795575243316608e+05, 1.56932489940248512e+05,
    1.61140965195830048e+05, 1.65564070780028256e+05, 1.70312554701480352e+05,
    1.75486986284790656e+05, 1.81187218697219744e+05, 1.87518050413513344e+05,
    1.94593473563783808e+05, 2.02540389901047584e+05, 2.11500759107428064e+05,
    2.21643078023966592e+05]
                
    numPhxDeps = len(phxSunPGas64) #//yeah, I know, 64, but that could change!
    logPhxSunPGas64 = [0.0 for i in range(numPhxDeps)]
    for i in range(len(phxSunPGas64)):
        logPhxSunPGas64[i] = math.log(phxSunPGas64[i]);
     
    logPhxSunTau64 = getLogPhxSunTau64()    
#// interpolate onto gS3 tauRos grid and re-scale with Teff:
        
    phxSunPGas = [0.0 for i in range(numDeps)]
    logPhxSunPGas = [0.0 for i in range(numDeps)]
    scalePGas = [ [0.0 for i in range(numDeps)] for j in range(2)]
    for i in range(numDeps):
        logPhxSunPGas[i] = ToolBox.interpol(logPhxSunTau64, logPhxSunPGas64, tauRos[1][i])
        scalePGas[1][i] = logEg + logPhxSunPGas[i] - phxSunLogEg()
        #print("i ", i, " scalePGas[1][i] ", scalePGas[1][i])
        scalePGas[0][i] = math.exp(scalePGas[1][i])
        #//System.out.println("scalePGas[1][i] " + logE * scalePGas[1][i])
        
    return scalePGas

def phxSunNe(grav, numDeps, tauRos, scaleTemp, kappaScale): 

    logE = math.log10(math.e)
    logEg = math.log(grav) #//base e!
    logEkappaScale = math.log(kappaScale);

#//Theoretical radiative/convective model from Phoenix V15:
    phxSunPe64 = [
    1.53086468021591745e-07, 5.66518458165471424e-03, 6.72808433760886656e-03,
    8.00271552708326656e-03, 9.51809762875982208e-03, 1.13117438884935648e-02,
    1.34299756939525680e-02, 1.59287848014678144e-02, 1.88751877391284448e-02,
    2.23491173128862976e-02, 2.64457686695698400e-02, 3.12779350532322240e-02,
    3.69791374171045888e-02, 4.37078139287801024e-02, 5.16503829681397248e-02,
    6.10221573903118336e-02, 7.20768505868849536e-02, 8.51123959415642752e-02,
    1.00475763241309840e-01, 1.18571138726675232e-01, 1.39870552376136714e-01,
    1.64923053015554560e-01, 1.94357063774820192e-01, 2.28928720249475840e-01,
    2.69525262128246720e-01, 3.17192228891198592e-01, 3.73192988074577856e-01,
    4.39058414038311360e-01, 5.16615873984964544e-01, 6.08066526878471680e-01,
    7.16264581324812416e-01, 8.44657163125294336e-01, 9.97267452897639808e-01,
    1.17915717019238848e+00, 1.39715732004723136e+00, 1.66026825646718432e+00,
    1.97886823850223904e+00, 2.36716912384854112e+00, 2.84540915928013805e+00,
    3.44853013665125120e+00, 4.21529199485384704e+00, 5.21488490421314560e+00,
    6.56660005867586432e+00, 8.55643059606379776e+00, 1.16931723772200080e+01,
    1.71629079266534368e+01, 2.75152019254691616e+01, 4.18720694941323264e+01,
    7.66283674228108288e+01, 1.45995186997127872e+02, 3.04766672331673792e+02,
    5.44151864837275328e+02, 8.17181982032739072e+02, 1.11216222784450608e+03,
    1.43633935534913856e+03, 1.79603721463325728e+03, 2.19692608617747040e+03,
    2.64548745663525184e+03, 3.14931730610757952e+03, 3.71721361233669376e+03,
    4.35932065708395904e+03, 5.08736399892079296e+03, 5.91634943413070720e+03,
    6.85104524590000384e+03]

    numPhxDeps = len(phxSunPe64)  #//yeah, I know, 64, but that could change!
    logPhxSunPe64 = [0.0 for i in range(numPhxDeps)]
    for i in range(len(phxSunPe64)):
        logPhxSunPe64[i] = math.log(phxSunPe64[i]);

    logPhxSunTau64 = getLogPhxSunTau64()
#// interpolate onto gS3 tauRos grid and re-scale with Teff:
    phxSunPe = [0.0 for i in range(numDeps)]
    logPhxSunPe = [0.0 for i in range(numDeps)]
    logScalePe = [0.0 for i in range(numDeps)]
    scaleNe = [ [0.0 for i in range(numDeps)] for j in range(2) ]

    for i in range(numDeps):
        logPhxSunPe[i] = ToolBox.interpol(logPhxSunTau64, logPhxSunPe64, tauRos[1][i])
        logScalePe[i] = logEg + logPhxSunPe[i] - phxSunLogEg() - logEkappaScale
        scaleNe[1][i] = logScalePe[i] - scaleTemp[1][i] - Useful.logK()
        scaleNe[0][i] = math.exp(scaleNe[1][i])
        #//System.out.println("scaleNe[1][i] " + logE * scaleNe[1][i]);
        

    return scaleNe

#//Try to recover the opacity as lambda_0 = 1200 nm:
def phxSunKappa(numDeps, tauRos, kappaScale):

    logEkappaScale = math.log(kappaScale)

    #//Theoretical radiative/convective model from Phoenix V15:
    phxSunRho64 = [
    4.13782346832222649e-16, 3.02095569469690462e-10, 3.54633225055968270e-10,
    4.15928280610231993e-10, 4.87569895799879155e-10, 5.71381142733345291e-10,
    6.69468927495419999e-10, 7.84278468388299388e-10, 9.18654436245877140e-10,
    1.07590983297567878e-09, 1.25990158939278389e-09, 1.47513757382262481e-09,
    1.72688539188771193e-09, 2.02128936476074103e-09, 2.36554000030610158e-09,
    2.76809615861929229e-09, 3.23884396019102352e-09, 3.78934920783997866e-09,
    4.43317360103421215e-09, 5.18621173362546736e-09, 6.06707380164391496e-09,
    7.09757215466433105e-09, 8.30337600953291647e-09, 9.71426731449415417e-09,
    1.13650770268615465e-08, 1.32964932176367733e-08, 1.55557163673284530e-08,
    1.81974840999693492e-08, 2.12855768344032029e-08, 2.48940684847852482e-08,
    2.91068454381155637e-08, 3.40213170202104799e-08, 3.97519122004400661e-08,
    4.64290866159173997e-08, 5.41967343519845744e-08, 6.32144869975830899e-08,
    7.36729431582295057e-08, 8.57774421976652924e-08, 9.97399445761737017e-08,
    1.15721981027072251e-07, 1.33967659681056212e-07, 1.54620178670780798e-07,
    1.77690495649821781e-07, 2.02608223525831620e-07, 2.28481547026651195e-07,
    2.53309018291389784e-07, 2.74195019891415717e-07, 2.94373976046088894e-07,
    3.05614181338722779e-07, 3.09912387277346887e-07, 3.05484245799381785e-07,
    3.00519445088246902e-07, 2.98007120264342719e-07, 2.97336159154754909e-07,
    2.97854109132361140e-07, 2.99327766949861546e-07, 3.01691329467384893e-07,
    3.04944348605014908e-07, 3.09125225055924192e-07, 3.14302162196028050e-07,
    3.20569231575000568e-07, 3.28044919674719785e-07, 3.36858977566225440e-07,
    3.47271781807407172e-07]

    phxSunRadius64 = [
    9.98760000000000000e+10, 9.98660572490945152e+10, 9.98645871807186304e+10,
    9.98631098643980160e+10, 9.98616245003269760e+10, 9.98601306458076032e+10,
    9.98586280682994048e+10, 9.98571166428681216e+10, 9.98555962828737792e+10,
    9.98540668955362944e+10, 9.98525283799022080e+10, 9.98509805586940416e+10,
    9.98494232096872704e+10, 9.98478561022866944e+10, 9.98462789875034752e+10,
    9.98446916403608064e+10, 9.98430938763377024e+10, 9.98414855440511616e+10,
    9.98398665329129600e+10, 9.98382367818977152e+10, 9.98365962762478464e+10,
    9.98349450434777856e+10, 9.98332831693342848e+10, 9.98316107506358144e+10,
    9.98299278514395904e+10, 9.98282344996977408e+10, 9.98265306530218624e+10,
    9.98248161610832000e+10, 9.98230907740896512e+10, 9.98213541602841472e+10,
    9.98196058171550848e+10, 9.98178450629064448e+10, 9.98160711682936960e+10,
    9.98142833645708160e+10, 9.98124806655167488e+10, 9.98106617425436544e+10,
    9.98088251712680448e+10, 9.98069695137641728e+10, 9.98050931816160256e+10,
    9.98031941655960192e+10, 9.98012715884401664e+10, 9.97993275763000704e+10,
    9.97973698841808128e+10, 9.97954186071983616e+10, 9.97935139683624704e+10,
    9.97917197632915456e+10, 9.97901090800493440e+10, 9.97886636105590528e+10,
    9.97874361487011200e+10, 9.97864521731274880e+10, 9.97857037165240960e+10,
    9.97851119909964288e+10, 9.97845890321633152e+10, 9.97840832513957504e+10,
    9.97835682848038784e+10, 9.97830286360697344e+10, 9.97824528501662336e+10,
    9.97818311493811456e+10, 9.97811545315540224e+10, 9.97804143368400768e+10,
    9.97796020159347072e+10, 9.97787090120484864e+10, 9.97777268638511104e+10,
    9.97766460582020224e+10]

    numPhxDeps = len(phxSunRadius64)
    phxSunKappa64 = [0.0 for i in range(numPhxDeps)]
    logPhxSunKappa64 = [0.0 for i in range(numPhxDeps)]
    #//double[] logPhxSunRho64 = new double[numPhxDeps];
    #//double[] logPhxSunRadius64 = new double[numPhxDeps];
    phxSunTau64 = getPhxSunTau64()
    logPhxSunTau64 = getLogPhxSunTau64()

#//Fix to get right depth scale and right line strengths:
#// Yeah - everywhere ya go - opacity fudge
    fudge = 0.25
    logFudge = math.log(fudge)
    #double deltaRho, deltaRadius, deltaTau, logDeltaRho, logDeltaRadius, logDeltaTau;
    logE = math.log10(math.e)
    for i in range(1, numPhxDeps):

        #//Renormalize radii before taking difference
        #//Caution: Radius *decreases* with increasing i (inward) and we'll be taking the log:
        deltaRadius = (1.0e-11 * phxSunRadius64[i - 1]) - (1.0e-11 * phxSunRadius64[i])
        deltaRadius = abs(deltaRadius)
        #//restore to cm:
        deltaRadius = 1.0e11 * deltaRadius
        #//Renormalize before taking rho difference
        deltaRho = (1.0e9 * phxSunRho64[i]) - (1.0e9 * phxSunRho64[i - 1])
        deltaRho = abs(deltaRho)
        #//Restore g/cm^3:
        deltaRho = 1.0e-9 * deltaRho
        #//Renormalize before taking rho difference
        deltaTau = (1.0e2 * phxSunTau64[i]) - (1.0e2 * phxSunTau64[i - 1])
        deltaTau = abs(deltaTau)
        deltaTau = 1.0e-2 * deltaTau

        logDeltaRadius = math.log(deltaRadius)
        logDeltaRho = math.log(deltaRho)
        logDeltaTau = math.log(deltaTau)

        logPhxSunKappa64[i] = logDeltaTau - logDeltaRho - logDeltaRadius - logEkappaScale + logFudge
        phxSunKappa64[i] = math.exp(logPhxSunKappa64[i])
        #//System.out.println("logPhxSunKappa64[i] " + logE*logPhxSunKappa64[i]);

    logPhxSunKappa64[0] = logPhxSunKappa64[1]
    phxSunKappa64[0] = phxSunKappa64[1]

    #// interpolate onto gS3 tauRos grid and re-scale with Teff:
    phxSunKappa = [ [0.0 for i in range(numDeps)] for j in range(2)]
    for i in range(numDeps):
        phxSunKappa[1][i] = ToolBox.interpol(logPhxSunTau64, logPhxSunKappa64, tauRos[1][i])
        phxSunKappa[0][i] = math.exp(phxSunKappa[1][i])
        #//System.out.println("phxSunKappa[1][i], i= " + i + " " + logE * phxSunKappa[1][i]);

    return phxSunKappa