In [93]:

    

import numpy as np
import pandas as pd
from scipy import integrate


    

In [80]:

    

stellar_freqs = pd.DataFrame(np.loadtxt('Sun-freqs.dat', skiprows=1), 
                     columns=['n', 'l', 'nu_s', 'sigma'])
stellar_freqs.head()


    

    
Out[80]:






  

    

      

      
n
      
l
      
nu_s
      
sigma
    
  
  

    

      
0
      
6
      
0
      
972.613
      
0.002
    
    

      
1
      
7
      
1
      
1185.592
      
0.004
    
    

      
2
      
8
      
0
      
1263.162
      
0.012
    
    

      
3
      
8
      
1
      
1329.629
      
0.004
    
    

      
4
      
8
      
2
      
1394.680
      
0.011
    
  

In [81]:

    

model_freqs = pd.DataFrame(np.loadtxt('modelS/modelS.freq'),
                   columns=['l', 'n', 'nu_m', 'Q'])
model_freqs.head()


    

    
Out[81]:






  

    

      

      
l
      
n
      
nu_m
      
Q
    
  
  

    

      
0
      
0
      
2
      
404.4800
      
1.005840e-04
    
    

      
1
      
0
      
3
      
535.9387
      
2.491420e-05
    
    

      
2
      
0
      
4
      
680.5718
      
7.507410e-06
    
    

      
3
      
0
      
5
      
825.3639
      
2.596760e-06
    
    

      
4
      
0
      
6
      
972.7445
      
9.922740e-07
    
  

In [31]:

    

nu = pd.merge(nu_s, nu_m)
nu.head()


    

    
Out[31]:






  

    

      

      
n
      
l
      
nu_s
      
sigma
      
nu_m
      
Q
    
  
  

    

      
0
      
6
      
0
      
972.613
      
0.002
      
972.7445
      
9.922740e-07
    
    

      
1
      
7
      
1
      
1185.592
      
0.004
      
1185.6151
      
2.831820e-07
    
    

      
2
      
8
      
0
      
1263.162
      
0.012
      
1263.5229
      
1.797740e-07
    
    

      
3
      
8
      
1
      
1329.629
      
0.004
      
1329.6956
      
1.238890e-07
    
    

      
4
      
8
      
2
      
1394.680
      
0.011
      
1394.7049
      
8.553300e-08
    
  

In [84]:

    

nls = [(int(n), int(l)) 
       for (n,l) in np.array(nu[['n', 'l']])]


    

In [85]:

    

K = {(int(n),int(l)) :
     pd.DataFrame(np.loadtxt('modelS_ker/c^2-\\rho_l='+\
                      str(int(l))+'_n='+str(int(n))+'.dat'),
                  columns=['r', 'c2_m', 'rho_m'])
     for (n,l) in nls}


    

In [86]:

    

K[6,0].head()


    

    
Out[86]:






  

    

      

      
r
      
c2_m
      
rho_m
    
  
  

    

      
0
      
0.000000
      
0.000000e+00
      
-6.696783e-129
    
    

      
1
      
0.008292
      
1.700927e-13
      
-5.187623e-14
    
    

      
2
      
0.008359
      
1.728319e-13
      
-5.270772e-14
    
    

      
3
      
0.008426
      
1.756144e-13
      
-5.354725e-14
    
    

      
4
      
0.008494
      
1.784437e-13
      
-5.439765e-14
    
  

In [90]:

    

modelS = pd.DataFrame(np.loadtxt('c2_rho.dat'), 
                      columns=['r', 'c2', 'rho'])
c2_m = modelS['c2']
rho_m = modelS['rho']
modelS.head()


    

    
Out[90]:






  

    

      

      
r
      
c2
      
rho
    
  
  

    

      
0
      
1.436801e-60
      
2.541536e+15
      
154.236483
    
    

      
1
      
8.291977e-03
      
2.544074e+15
      
153.356344
    
    

      
2
      
8.358816e-03
      
2.544077e+15
      
153.344451
    
    

      
3
      
8.426194e-03
      
2.544148e+15
      
153.328338
    
    

      
4
      
8.494117e-03
      
2.544221e+15
      
153.311917
    
  

In [ ]:

    

def A(n, l, c2_s, rho_s):
    kernel = K[n,l]
    fx = kernel['c2_m'] * (c2_m - c2_s) / c2_m + \
         kernel['rho_m'] * (rho_m - rho_s) / rho_m
    return sp.integrate.simps(fx, x=kernel['r'])

def freq_diff(n, l):
    idx = nls.index((n,l))
    nu_m = nu['nu_m'][idx]
    nu_s = nu['nu_s'][idx]
    return (nu_m - nu_s) / nu_m

def chi_sq(c2_s, rho_s, F_surf):
    return sum(
        (freq_diff(n,l) - A(n,l,c2_s,rho_s) - \
         F_surf[nls.index((n,l))] / nu['Q'][nls.index((n,l))]
        )**2 / nu['sigma'][nls.index((n,l))]
               for (n,l) in nls)

def L_2(c2_s, rho_s):
    c2 = (c2_m - c2_s) / c2_m
    rho = (rho_m - rho_s) / rho_m
    d2c2_dr2 = np.diff(c2, n=2) / np.diff(modelS['r'], 2)
    d2rho_dr2 = np.diff(rho, n=2) / np.diff(modelS['r'], 2)
    return integrate.simps(d2c2_dr2**2 + d2rho_dr2**2, modelS['r'][2:])


    

In [119]:

    

A(nls[0][0], nls[0][1], c2_m, rho_m) # should be 0


    

    
Out[119]:





0.0

In [129]:

    

chi_sq(c2_m, rho_m, np.zeros(len(nls))) # should be close to 0


    

    
Out[129]:





0.0045089943862229465

In [148]:

    

L_2(c2_m, rho_m) # should be 0


    

    
Out[148]:





0.0