from astropy import units as u
from astropy.coordinates import SkyCoord
from scipy.ndimage.interpolation import map_coordinates
import os.path
from collections import Iterable
import pyfits as fits
def get_ebv_from_map(coordinates, interpolate=True, order=1):

    mapdir='/vol/Dust_Map/maps/'
    fname = os.path.join(mapdir, 'SFD_dust_4096_{0}.fits')

    # Parse input
    ra, dec = coordinates
    coordinates = SkyCoord(ra=ra, dec=dec, unit=(u.degree, u.degree))

    # Convert to galactic coordinates.
    coordinates = coordinates.galactic
    l = coordinates.l.radian
    b = coordinates.b.radian

    # Check if l, b are scalar
    return_scalar = False
    if not isinstance(l, Iterable):
        return_scalar = True
        l, b = np.array([l]), np.array([b])

    # Initialize return array
    ebv = np.empty_like(l)

    # Treat north (b>0) separately from south (b<0).
    for n, idx, ext in [(1, b >= 0, 'ngp'), (-1, b < 0, 'sgp')]:

        if not np.any(idx): continue
        hdulist = fits.open(fname.format(ext))
        mapd = hdulist[0].data

        # Project from galactic longitude/latitude to lambert pixels.
        # (See SFD98).
        npix = mapd.shape[0]
        x = (npix / 2 * np.cos(l[idx]) * np.sqrt(1. - n*np.sin(b[idx])) +
             npix / 2 - 0.5)
        y = (-npix / 2 * n * np.sin(l[idx]) * np.sqrt(1. - n*np.sin(b[idx])) +
             npix / 2 - 0.5)

        # Get map values at these pixel coordinates.
        if interpolate:
            ebv[idx] = map_coordinates(mapd, [y, x], order=order)
        else:
            x=np.round(x).astype(np.int)
            y=np.round(y).astype(np.int)
            ebv[idx] = mapd[y, x]

        hdulist.close()

    if return_scalar:
        return ebv[0]
    return ebv

ABtoVEGA=[+0.77,-0.13,-0.02,+0.19,+0.49,-0.19,-0.18,-0.06,-0.06,+0.04,+0.10,+0.22,+0.27,+0.36,+0.45,+0.56,+0.50]

ZP=np.loadtxt("zero-points.offsets.asc")
Gal=np.loadtxt("photometry_galaxies_photometry.cat", dtype=None)
ra=Gal[:,1]
dec=Gal[:,2]

Av=[3.836 ,3.070,2.542,2.058,1.313, 3.466 , 3.082 , 2.882 , 2.652 , 2.367 , 2.226 , 2.066 , 1.872 , 1.652 , 1.423 , 1.298 , 1.156]
count=0
for i in range(Gal.shape[0]):
    j=0
    m=39
    while (j<len(ABtoVEGA)):
         if ((isnan(float(Gal[i,j+5]))) or (isnan(float(Gal[i,j+22])))):
            Gal[i,j+5]='-99.0'
            Gal[i,j+22]='0.0'
            count+=1
         else:
            value=get_ebv_from_map((ra[i], dec[i]), interpolate=True, order=1)
            Gal[i,j+5]=str(float(Gal[i,j+5])-ABtoVEGA[j]-ZP[j]-value*Av[j])
         j+=1
    while (m<45):
          if (math.isnan(float( Gal[i,m]))):
              Gal[i,m] = '0.0'
          m+=1

for n, idx, ext in [(1, b >= 0, 'ngp'), (-1, b < 0, 'sgp')]:

    if not np.any(idx): continue
    hdulist = fits.open(fname.format(ext))  # this one
    mapd = hdulist[0].data

def main_loop(a):
    for n, idx, ext in [(1, b >= 0, 'ngp'), (-1, b < 0, 'sgp')]:
        # Rest of the code

for n, idx, ext in [(1, b >= 0, 'ngp'), (-1, b < 0, 'sgp')]:

    # Other necessary code
    if ext not in mapdcache:
        with fits.open(fname.format(ext)) as hdulist:
            mapdcache[ext] = hdulist[0].data
    else:
        mapd = mapdcache[ext]