Isochrone interpolation yields incorrect results

[balbinot]

I have been trying to use the interpolation scheme fro the Dartmouth models, but this yield some strange results.

Tested in the following environment

1.0.0 >= scipy >= 0.12.0 numpy= 1.13.3 pandas= 0.20.3 python 2.7.3

from matplotlib import pyplot as plt
import numpy as np

from isochrones.dartmouth import Dartmouth_Isochrone, DartmouthModelGrid

dar = Dartmouth_Isochrone()
d = DartmouthModelGrid(['g', 'r']).df

g = dar.mag['g']
r = dar.mag['r']

## Original grid
plt.figure(figsize=(4,6))
j = (d['age']==9)&(d['feh']==-2.5)
plt.plot(d['g'][j]-d['r'][j], d['g'][j], 'k-', lw=2, label="1.00 Gyr")
j = (d['age']>9)&(d['age']< 9.097)&(d['feh']==-2.5)
plt.plot(d['g'][j]-d['r'][j], d['g'][j], 'b-', lw=2, label="1.25 Gyr")

## Interpolated model
mm = np.arange(0.1,3, 0.00001)
tg = g(mm, 9.05, -2.5)
tr = r(mm, 9.05, -2.5)

plt.plot(tg-tr,tg,'r-', lw=2, label='1.12Gyr (interp)')
plt.ylim(plt.ylim()[::-1])
plt.legend(loc='best')
plt.xlabel('g-r')
plt.ylabel('g')
plt.savefig('test.png')

Sorry for the chopped axis label.

[timothydmorton]

Yikes! That's pretty nasty. This appears to be a good illustration of what I've suspected, that the interpolation in mass-age-feh might yield questionable results off the main sequence, when things are changing rapidly. When I get some more time I'll take a look deeper into this. Thanks for pointing this out-- I think something like this would be a good test case for interpolation schemes.

I will note that a similar exercise (at feh=-2.0) using the MIST model grid (which uses a custom interpolation scheme rather than scipy) appears more sane, perhaps because the model grids are denser?

from matplotlib import pyplot as plt
import numpy as np

from isochrones.dartmouth import Dartmouth_Isochrone, Dartmouth_FastIsochrone, DartmouthModelGrid
from isochrones.mist import MIST_Isochrone, MISTModelGrid

# ic = Dartmouth_Isochrone()
# grid = DartmouthModelGrid
ic = MIST_Isochrone()
grid = MISTModelGrid

d = grid(['g', 'r']).df

# Give MISTModelGrid 'age' column
if 'age' not in d.columns:
    d = d.rename(columns={'log10_isochrone_age_yr':'age'})

g = ic.mag['g']
r = ic.mag['r']

feh = -2.0

## Original grid
plt.figure(figsize=(4,6))
j = (d['age']==9)&(d['feh']==feh)
plt.plot(d['g'][j]-d['r'][j], d['g'][j], 'k-', lw=2, label="1.00 Gyr")
j = (d['age']>9)&(d['age']< 9.097)&(d['feh']==feh)
plt.plot(d['g'][j]-d['r'][j], d['g'][j], 'b-', lw=2, label="1.25 Gyr")

## Interpolated model
mm = np.concatenate([np.arange(0.15,3, 0.0001)])
tg = g(mm, 9.05, feh)
tr = r(mm, 9.05, feh)

plt.plot(tg-tr,tg,'r.', ms=1, lw=2, label='1.12Gyr (interp)')
plt.ylim(plt.ylim()[::-1])
plt.legend(loc='best')
plt.xlabel('g-r')
plt.ylabel('g');

[balbinot]

So, for Dartmouth I found that they provide this EEP (equivalent evolutionary point) column. If you use that to interpolate, instead of mass, a sane result comes out. So basically in the Delaunay interpolation you need to change points[:,0] = eep

That produces something like this for Dartmouth.

I guess than the user could just produce an interpolation for the mass and use that to map to eep.

This is the code that generates this (sorry very messy; allinone3.dat is an ascii table containing all the darthmouth models with fixed afe and Y)

#!/usr/bin/env python
#-*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as p
from scipy.interpolate import LinearNDInterpolator as interpnd

a,Z,Zeff,feh,afe,Y,m,logg,V,I,eep = np.loadtxt('allinone3.dat',
                                               usecols=(0,1,2,3,4,5,7,9,16,21,6),
                                               unpack=True)

## Trim models for speed
j = (a>=5)*(feh<=-0.5)*(m>0.6)
N = len(a[j])
pts = np.zeros((N,3))
pts[:,0] = eep[j]
pts[:,1] = a[j]
pts[:,2] = feh[j]

jj = (a==9)*(feh==-2.5)
p.plot(V[jj]- I[jj], I[jj], '-', lw=2, label='age = 9 Gyr')
jj = (a==10)*(feh==-2.5)
p.plot(V[jj]- I[jj], I[jj], '-', lw=2, label='age = 10 Gyr')
jj = (a==9.5)*(feh==-2.5)
p.plot(V[jj]- I[jj], I[jj], '-', lw=2, label='age = 9.5 Gyr')

iV = interpnd(pts, V[j])
iI = interpnd(iV.tri, I[j])

np.save('asd.tri', iV.tri)

FEH = -2.5
AGE = 10.0
ms = np.arange(1,300, 0.1) ## actually eep
for age in [9.25, 9.75]:
    p.plot(iV(ms, age,FEH)-iI(ms, age,FEH), iI(ms, age,FEH), 'k:', label='age = {:.2f} Gyr'.format(age))

p.xlabel('V-I')
p.ylabel('V')
p.legend(loc='best')
p.ylim(p.ylim()[::-1])
p.show()
exit()

[timothydmorton]

Thanks for investigating this. I hope to have time soon to experiment with how to integrate EEP interpolation into this code.