#!/usr/bin/python
###################################
# Author: Julien Peloton, j.peloton@sussex.ac.uk
# Pointing module: thanks Neil Goeckner-Wald!
###################################
import ConfigParser
import sys
import ephem
import copy
import os

import numpy as np
import healpy as hp
from . import util_CMB
from . import InputScan

from util_CMB import benchmark
try:
    from scipy import weave
except:
    ## weave has been removed from scipy version > 0.18
    import weave

try:
    from InsideMe import profiler
except:
    ## InsideMe will be released soon ;-)
    ## Stay tuned for a new great package
    class profiler():
        @staticmethod
        def benchmark(field=''):
            def outer_wrapper(func):
                def inner_wrapper(*args, **kwargs):
                    return func(*args, **kwargs)
                return inner_wrapper
            return outer_wrapper

# ## numerical constants
radToDeg = 180 / np.pi
sidDayToSec = 86164.0905

def generate_scans(config_file, env=None):
    """
    Main script to generate the scanning strategy.

    Parameters
    ----------
        * config_file: ini file, parser containing scanning parameters
            (see the setup_scanning.ini provided)
        * env: object, contain paths and debugging tools [Optional]
    """
    if env is None:
        env = lambda: -1
        env.verbose = False
        env.plot = False
        env.out_name = 'outputs'
        env.out_path = os.path.join('./', env.out_name)
        env.outpath_noise = os.path.join('./', env.out_name, 'noise')
        env.outpath_masks = os.path.join('./', env.out_name, 'masks')

    ## Load the parameters from the parser
    Config = ConfigParser.ConfigParser()
    Config.read(config_file)
    scanning = util_CMB.normalise_scanning_parser(
        Config._sections['ScanningParameters'])
    checks_param(scanning)

    ## Save ini file for later comparison
    path = os.path.join(env.outpath_masks, 'setup_scanning.ini')
    with open(path, 'w') as configfile:
        Config.write(configfile)

    ## Load experiment technical details
    experience = instrument(scanning.name_instrument)

    ## Initialise hdf5 file
    scan_info = InputScan.HealpixMapInfo(
        npix=hp.nside2npix(scanning.nside),
        obspix=np.zeros(hp.nside2npix(scanning.nside)),
        nside=scanning.nside,
        source=scanning.name_instrument)
    scan = InputScan.InputScan(scan_info)

    ## Loop over CESes and generate scans
    for n in range(scanning.length_of_cycle):
        ## Initialise the starting date of observation
        ## It will be updated then automatically for the run
        ## Do not run CES longer than 24h though!
        experience.inst.date = scanning.starting_date
        experience.inst.date += n * 24 * ephem.second * 3600
        position_in_the_cycle = n % scanning.length_of_cycle

        ## Grab the details of this CES only
        myargs = copy.copy(scanning)
        for myarg in myargs.__dict__.keys():
            if type(myargs.__dict__[myarg]) == np.ndarray:
                myargs.__dict__[myarg] = \
                    myargs.__dict__[myarg][position_in_the_cycle]

        ## Create the scan strategy
        create_scan_strategy(myargs, experience, scan)

    ## Convolve with the focal plane and extrapolate
    ## to obtain the full dataset (total number of days)
   
 scan.nhit, scan.cc, scan.ss, scan.cs = c
onvolve_focalplane(
        scan.nhit, scan.cc, scan.ss, scan.cs,
        nbolos=myargs.nbolos,
        fp_radius_amin=myargs.fp_radius_amin,
        boost=myargs.number_of_days/float(myargs.length_of_cycle))

    if env.plot is True:
        import pylab as pl
        hp.mollview(scan.nhit, sub=221)
        hp.mollview(scan.cc, sub=222)
        hp.mollview(scan.ss, sub=223)
        hp.mollview(scan.cs, sub=224)
        pl.show()

    scan.mapinfo.obspix = np.where(scan.nhit != 0)[0]
    scan.mapinfo.npix = len(scan.mapinfo.obspix)
    scan.nhit = scan.nhit[scan.mapinfo.obspix]
    scan.w = np.array(scan.nhit, dtype=int)
    scan.cc = scan.cc[scan.mapinfo.obspix]
    scan.ss = scan.ss[scan.mapinfo.obspix]
    scan.cs = scan.cs[scan.mapinfo.obspix]
    scan.mapinfo.source = scanning.name_instrument

    scan.save(os.path.join(env.outpath_masks,env.out_name+'.hdf5'))

class instrument(object):
    """ Class to handle different instruments and their specificities """
    def __init__(self, instrument='so_deep'):
        """ Update here with new informations from TWGs. """

        self.instrument_ok = ['so_deep', 'so_shallow', 'polarbear', 'actpol']

        if instrument.lower() in self.instrument_ok:
            self.inst = self.chili_atacama()
        else:
            print 'Instrument not yet implemented\n'
            print 'Choose among ', self.instrument_ok
            sys.exit()

    def chili_atacama(self):
        """
        Put here instrument constants
        """
        ## Polarbear observatory location
        pb = ephem.Observer()
        pb.long = '-67:46.816'
        pb.lat = '-22:56.396'
        pb.elevation = 5200

        return pb

@profiler.benchmark(field='Convolve focal plane')
def convolve_focalplane(
        bore_nHits, bore_cos, bore_sin, bore_cs,
        nbolos, fp_radius_amin, boost):
    """
    Given a nHits and bore_cos and bore_sin map,
    perform the focal plane convolution.
    Original author: Neil Goeckner-Wald.
    Modifications by Julien Peloton.

    Parameters
    ----------
        * bore_nHits: 1D array, number of hits for the reference detector
        * bore_cos: 1D array, cumulative cos**2 for the reference detector
        * bore_sin: 1D array, cumulative sin**2 for the reference detector
        * bore_cs: 1D array, cumulative cos*sin for the reference detector
        * nbolos: int, total number of bolometers desired
        * fp_radius_amin: float, radius of the focal plane in arcmin
        * boost: float, boost factor to artificially
            increase the number of hits.
            It doesn't change the shape of the survey (just the amplitude)

    Outputs
    ----------
        * focalplane_nHits: 1D array, number of hits
            for the all the detectors
        * focalplane_cos: 1D array, cumulative cos**2
            for the all the detectors
        * focalplane_sin: 1D array, cumulative sin**2
            for the all the detectors
        * focalplane_cossin: 1D array, cumulative cos*sin
            for the all the detectors

    """
    #Now we want to make the focalplane maps
    focalplane_nHits = np.zeros(bore_nHits.shape)
    focalplane_sin = np.zeros(bore_sin.shape)
    focalplane_cos = np.zeros(bore_cos.shape)
    focalplane_cossin = np.zeros(bore_cs.shape)

    #Resolution of our healpix map
    nside = hp.npix2nside(focalplane_nHits.shape[0])
    resol_amin = hp.nside2resol(nside, arcmin=True)
    fp_rad_bins = int(fp_radius_amin * 2. / resol_amin)
    fp_diam_bins = (fp_rad_bins * 2) + 1

    #Build the focal plane model and a list of offsets
    (x_fp, y_fp) = np.array(
        np.unravel_index(
            range(0, fp_diam_bins**2),
            (fp_diam_bins, fp_diam_bins))).reshape(
                2, fp_diam_bins, fp_diam_bins) - (fp_rad_bins)
    fp_map = ((x_fp**2 + y_fp**2) < (fp_rad_bins)**2)

    bolo_per_pix = nbolos / float(np.sum(fp_map))

    dRA = np.ndarray.flatten(
        (x_fp[fp_map].astype(float) * fp_radius_amin) / (
            fp_rad_bins * 60. * (180. / (np.pi))))
    dDec = np.ndarray.flatten(
        (y_fp[fp_map].astype(float) * fp_radius_amin) / (
            fp_rad_bins * 60. * (180. / (np.pi))))

    pixels_global = np.array(np.where(bore_nHits != 0)[0], dtype=int)
    for n in pixels_global:
        n = int(n)

        #compute pointing offsets
        (theta_bore, phi_bore) = hp.pix2ang(nside, n)
        phi = phi_bore + dRA * np.sin(theta_bore)
        theta = theta_bore + dDec

        pixels = hp.ang2pix(nside, theta, phi)

        ## Necessary because the values in pixels aren't necessarily unique!
        ## This is a poor design choice and should probably be fixed
        c_code = r"""
        int i, pix;
        for (i=0;i<npix_loc;i++)
        {
            pix = pixel