def plotmdp():
    spin00_pol_degree_spline = buildspline(0.5)
    spin00_mcube =  xBinnedModulationCube(fetch_mcubepath(0.5))

    spin998_pol_degree_spline = buildspline(0.998)
    spin998_mcube =  xBinnedModulationCube(fetch_mcubepath(0.998))

    spin00_mcube.fit()
    spin998_mcube.fit()
    
    spin00_fit_results = spin00_mcube.fit_results[0]
    spin998_fit_results = spin998_mcube.fit_results[0]

    
    plt.figure('MDP')
    
    spin00_mdp = spin00_mcube.mdp99[:-1]
    spin998_mdp = spin998_mcube.mdp99[:-1]
    emean = spin00_mcube.emean[:-1]
    emin =  spin00_mcube.emin[:-1]
    emax =  spin00_mcube.emax[:-1]
    width = (emax-emin)/2.
    plt.errorbar(emean,spin00_mdp,xerr=width, label='Spin 0.5',marker='o',linestyle='--')
    
    plt.errorbar(emean,spin998_mdp,xerr=width, label='Spin 0.998',marker='o',linestyle='--')
            
    plt.figtext(0.2, 0.85,'XIPE %s ks'%((SIM_DURATION*NUM_RUNS)/1000.),size=18)
    plt.xlim([1,10])
    
    plt.ylabel('MPD 99\%')
    plt.xlabel('Energy (keV)')
    plt.legend()
    plt.show()

def plot_doc():
    """
    """
    img = xFITSImage(os.path.join(XIMPOL_DATA, 'tycho_cmap.fits'))
    fig = img.plot(show=False)
    xFITSImage.add_label(fig, 'XIPE 300 ks')
    plt.show()

def view():
      
    _mcube = xBinnedModulationCube(MCUBE_FILE_PATH)
    _mcube.fit()
    _fit_results = _mcube.fit_results[0]
    plt.figure('Polarization degree')
    _mcube.plot_polarization_degree(show=False, color='blue')
    pol_degree_spline.plot(color='lightgray',label='Spin %s'%spindegree, show=False)
    plt.figtext(0.2, 0.85,'XIPE %s ks'%(SIM_DURATION/1000.),size=18)
    #plt.errorbar(_energy_mean, _pol_deg, yerr=_pol_deg_err, color='blue',marker='o')
    
    plt.legend()

    plt.figure('Polarization angle')
    _mcube.plot_polarization_angle(show=False, color='blue', degree=False)
    pol_angle_spline.plot(color='lightgray',label='Spin %s'%spindegree, show=False)
    plt.figtext(0.2, 0.85,'XIPE %s ks'%(SIM_DURATION/1000.),size=18)
    #plt.errorbar(_energy_mean,_pol_angle, yerr= _pol_angle_err,color='blue',marker='o')
    plt.xlim([1,10])
    plt.legend()
    plt.figure('MDP %s'%base_name)
    mdp = _mcube.mdp99[:-1]
    emean = _mcube.emean[:-1]
    emin =  _mcube.emin[:-1]
    emax =  _mcube.emax[:-1]
    width = (emax-emin)/2.
    plt.errorbar(emean,mdp,xerr=width, label='MDP99',marker='o',linestyle='--')
    plt.figtext(0.2, 0.85,'XIPE %s ks'%(SIM_DURATION/1000.),size=18)
    plt.xlim([1,10])
    plt.ylabel('MPD 99\%')
    plt.xlabel('Energy (keV)')
    #plt.legend()
    plt.show()

def makeMDPComparisonPlot(imaging_file_path):
    
    non_imaging_file_path = imaging_file_path.replace('_imaging.txt','_non_imaging.txt')
    print "Using %s for non imaging file path"%non_imaging_file_path
    
    _phase_ave, _phase_err, mdp_imaging = numpy.loadtxt(imaging_file_path, unpack=True)
    _phase_ave, _phase_err, mdp_nonimaging = numpy.loadtxt(non_imaging_file_path, unpack=True)
    scale_factor = 10.
    print "Improvement with imaging"
    for i, phase in enumerate(_phase_ave):
        imaging = 100*mdp_imaging[i]*(1/numpy.sqrt(scale_factor))
        non_imaging = 100*mdp_nonimaging[i]*(1/numpy.sqrt(scale_factor))
        print "%s\t Imaging (non):%s (%s) %s"%(phase,imaging,non_imaging,non_imaging/imaging)
        
   
    sim_label_imaging = 'XIPE %s ks\n Imaging 15"' % (SIM_DURATION*scale_factor/1000.)
    sim_label_nonimaging = 'Non imaging'
    lc_label = 'Light curve'
    plt.errorbar(_phase_ave, 100*mdp_imaging*(1/numpy.sqrt(scale_factor)),xerr=_phase_err, label=sim_label_imaging,fmt='o',markersize=6)
    
    plt.errorbar(_phase_ave, 100*mdp_nonimaging*(1/numpy.sqrt(scale_factor)),xerr=_phase_err, label=sim_label_nonimaging,fmt='v',markersize=6)
    
    pl_normalization_spline.plot(scale=10., show=False,
                                color='darkgray',label=lc_label)
    #on_phase = 0.25, 0.45
    #off_phase = 0.45,0.9
    plt.axvspan(0.25, 0.45, color='r', alpha=0.45, lw=0)
    plt.axvspan(0.45, 0.9, color='gray', alpha=0.25, lw=0)
    
    plt.ylabel('MDP 99\%')
    plt.legend()
    plt.savefig('crab_complex_mdp_imaging_vs_nonimaging_%i_shaded.png'%(SIM_DURATION*scale_factor/1000.))
    plt.show()

def plot_swift_lc(grb_list,show=True):
    """Plots Swift GRB light curves.
    """
    plt.figure(figsize=(10, 8), dpi=80)
    plt.title('Swift XRT light curves')
    num_grb = 0
    for grb_name in grb_list:
        flux_outfile = download_swift_grb_lc_file(grb_name, min_obs_time=21600)
        if flux_outfile is not None:
            integral_flux_spline = parse_light_curve(flux_outfile)
            if integral_flux_spline is not None:
                if grb_name == 'GRB 130427A':
                    integral_flux_spline.plot(num_points=1000,logx=True,\
                                              logy=True,show=False,\
                                              color="red",linewidth=1.0)
                    num_grb += 1
                else:
                    c = random.uniform(0.4,0.8)
                    integral_flux_spline.plot(num_points=1000,logx=True,\
                                              logy=True,show=False,\
                                              color='%f'%c,linewidth=1.0)
                    num_grb += 1
        else:
            continue
    logger.info('%i GRBs included in the plot.'%num_grb)
    if show:
        plt.show()

 def plot_polarization_angle(self, show=True, degree=False, **kwargs):
     """Plot the polarization angle as a function of energy.
     """
     if self.fit_results == []:
         self.fit()
     _x = self.emean
     _dx = [self.emean - self.emin, self.emax - self.emean]
     if degree:
         _y = [numpy.degrees(r.phase) for r in self.fit_results]
         _dy = [numpy.degrees(r.phase_error) for r in self.fit_results]
     else:
         _y = [(r.phase) for r in self.fit_results]
         _dy = [(r.phase_error) for r in self.fit_results]
     # If there's more than one energy binning we also fit the entire
     # energy interval, but we don't want the corresponding data point to
     # appear in the plot, se we brutally get rid of it.
     if len(self.fit_results) > 1:
         _x = _x[:-1]
         _dx = [_x - self.emin[:-1], self.emax[:-1] - _x]
         _y = _y[:-1]
         _dy = _dy[:-1]
     plt.errorbar(_x, _y, _dy, _dx, fmt='o', **kwargs)
     plt.xlabel('Energy [keV]')
     if degree:
         plt.ylabel('Polarization angle [$^\circ$]')
     else:
         plt.ylabel('Polarization angle [rad]')
     if show:
         plt.show()

def plot_doc():
    """
    """
    img = xFITSImage(os.path.join(XIMPOL_DATA, 'crab_complex_cmap.fits'))
    fig = img.plot(show=False)
    xFITSImage.add_label(fig, 'XIPE %d ks' %DURATION/1000.)
    plt.show()

def main():
    """
    """
    plt.figure(figsize=(10, 6), dpi=80)
    plt.title('Swift XRT light curves of GRBs up to 130427A')
    #get_all_swift_grb_names = ['GRB 130427A','GRB 041223']
    num_grb = 0
    for i,grb_name in enumerate(get_all_swift_grb_names()):
        flux_outfile = download_swift_grb_lc_file(grb_name)
        if type(flux_outfile) is str:
            integral_flux_spline = parse_light_curve(flux_outfile)
            if integral_flux_spline != 0:
                if grb_name == 'GRB 130427A':
                    integral_flux_spline.plot(num_points=1000,logx=True,\
                                              logy=True,show=False,\
                                              color="red",linewidth=1.0)
                    num_grb += 1
                    break
                    plt.title('Swift XRT light curves of GRBs up to now')
                else:
                    c = random.uniform(0.4,0.8)
                    integral_flux_spline.plot(num_points=1000,logx=True,\
                                              logy=True,show=False,\
                                              color='%f'%c,linewidth=1.0)
                    num_grb += 1
    print num_grb
    plt.show()

 def test_plots(self):
     """
     """
     model = xpeInterstellarAbsorptionModel()
     plt.figure()
     model.xsection.plot(logx=True, logy=True, show=False)
     save_current_figure('gabs_xsection.png', clear=False)
     plt.figure()
     model.xsection_ecube().plot(logx=True, show=False)
     save_current_figure('gabs_xsection_ecube.png', clear=False)
     plt.figure()
     ra, dec = 10.684, 41.269
     column_density = mapped_column_density_HI(ra, dec, 'LAB')
     trans = model.transmission_factor(column_density)
     trans.plot(logx=True, show=False, label='$n_H = $%.3e' % column_density)
     plt.legend(loc='upper left')
     save_current_figure('gabs_trans_lab.png', clear=False)
     plt.figure()
     for column_density in [1.e20, 1.e21, 1.e22, 1.e23]:
         trans = model.transmission_factor(column_density)
         trans.plot(logx=True, show=False, label='$n_H = $%.1e' %\
                    column_density)
     plt.legend(loc='upper left')
     save_current_figure('gabs_trans_samples.png', clear=False)
     if sys.flags.interactive:
         plt.show()

    def plot(self, show=True):
        """Plot the energy dispersion.
        """
        from ximpol.utils.matplotlib_ import pyplot as plt
        from ximpol.utils.matplotlib_ import context_two_by_two
        emin = self.matrix.xmin()
        emax = self.matrix.xmax()

        def _plot_vslice(energy, position):
            """Convenience function to plot a generic vertical slice of the
            energy dispersion.
            """
            ax = plt.subplot(2, 2, position)
            vslice = self.matrix.vslice(energy)
            vslice.plot(overlay=False, show=False)
            plt.text(0.1, 0.9, '$E = %.2f\\ \\rm{keV}$' % energy,
                     transform=ax.transAxes)
            ppf = vslice.build_ppf()
            eres = 0.5*(ppf(0.8413) - ppf(0.1586))/ppf(0.5)
            plt.text(0.1, 0.85, '$\sigma_E/E = %.3f$' % eres,
                     transform=ax.transAxes)

        with context_two_by_two():
            plt.figure(1)
            ax = plt.subplot(2, 2, 1)
            self.matrix.plot(show=False)
            ax = plt.subplot(2, 2, 2)
            self.ebounds.plot(overlay=False, show=False)
            _plot_vslice(emin + 0.333*(emax - emin), 3)
            _plot_vslice(emin + 0.666*(emax - emin), 4)
        if show:
            plt.show()

 def plot(self, show=True):
     """Overloaded plot method.
     """
     plt.errorbar(self.time, self.counts/self.timedel,
                  yerr=self.error/self.timedel, fmt='o')
     plt.xlabel('Time [s]')
     plt.ylabel('Rate [Hz]')
     if show:
         plt.show()

 def plot(self, show=True):
     """Overloaded plot method.
     """
     fig = plt.figure('Phasogram')
     plt.errorbar(self.phase, self.counts, yerr=self.error, fmt='o')
     plt.xlabel('Phase')
     plt.ylabel('Counts/bin')
     if show:
         plt.show()

def plot(save_plots=False):
    """
    """
    sim_label = 'XIPE %s ks' % (SIM_DURATION/1000.)
    mod_label = 'Input model'
    _phase, _phase_err, _pol_deg, _pol_deg_err, _pol_angle,\
        _pol_angle_err = numpy.loadtxt(ANALYSIS_FILE_PATH, unpack=True)
    _colors = ['blue']*len(_pol_deg)
    plt.figure('Polarization degree')
    _good_fit = _pol_deg > 2*_pol_deg_err
    _bad_fit = numpy.logical_not(_good_fit)
    plt.errorbar(_phase[_good_fit], _pol_deg[_good_fit],
                 xerr=_phase_err[_good_fit], yerr=_pol_deg_err[_good_fit],
                 fmt='o', label=sim_label, color='blue')
    plt.errorbar(_phase[_bad_fit], _pol_deg[_bad_fit],
                 xerr=_phase_err[_bad_fit], yerr=_pol_deg_err[_bad_fit],
                 fmt='o', color='gray')
    pol_degree_spline.plot(show=False, label=mod_label, color='green')
    plt.axis([0., 1., 0., 0.1])
    plt.legend(bbox_to_anchor=(0.37, 0.95))
    plt.figtext(0.6, 0.8, '%.2f--%.2f keV' %\
                (E_BINNING[0], E_BINNING[-1]), size=16)
    if save_plots:
        plt.savefig('gk_per_polarization_degree.png')
    plt.figure('Polarization angle')
    plt.errorbar(_phase[_good_fit], _pol_angle[_good_fit],
                 xerr=_phase_err[_good_fit], yerr=_pol_angle_err[_good_fit],
                 fmt='o', label=sim_label, color='blue')
    plt.errorbar(_phase[_bad_fit], _pol_angle[_bad_fit],
                 xerr=_phase_err[_bad_fit], yerr=_pol_angle_err[_bad_fit],
                 fmt='o', color='gray')
    pol_angle_spline.plot(show=False, label=mod_label, color='green',
                          scale=numpy.radians(1.))
    plt.axis([0., 1., -0.1, 1.5])
    plt.xlabel('Rotational phase')
    plt.ylabel('Polarization angle [rad]')
    plt.legend(bbox_to_anchor=(0.37, 0.95))
    plt.figtext(0.6, 0.8, '%.2f--%.2f keV' %\
                (E_BINNING[0], E_BINNING[-1]), size=16)
    if save_plots:
        plt.savefig('gk_per_polarization_angle.png')
    _ebinning = zip(E_BINNING[:-1], E_BINNING[1:])
    if len(_ebinning) > 1:
        _ebinning.append((E_BINNING[0], E_BINNING[-1]))
    for i, (_emin, _emax) in enumerate(_ebinning):
        plt.figure('Phasogram %d' % i)
        phasogram = xBinnedPhasogram(_phasg_file_path(i))
        _scale = phasogram.counts.sum()/phasogram_spline.norm()/\
                 len(phasogram.counts)
        phasogram_spline.plot(show=False, label=mod_label, scale=_scale,
                              color='green')
        phasogram.plot(show=False, color='blue', label=sim_label )
        plt.legend(bbox_to_anchor=(0.37, 0.95))
        plt.figtext(0.65, 0.8, '%.2f--%.2f keV' % (_emin, _emax), size=16)
        if save_plots:
            plt.savefig('gk_per_phasogram_%d.png' % i)
    plt.show()

def display():
    """Display the source model.
    """
    from ximpol.utils.matplotlib_ import pyplot as plt
    from ximpol.srcmodel.img import xFITSImage

    print(ROI_MODEL)
    fig = plt.figure('Energy spectrum')
    spectral_model_spline.plot(logy=True, show=False, label='Total')
    plt.show()

def plot(save=False):
    """Plot the stuff in the analysis file.
    """
    sim_label = 'XIPE %s ks' % (SIM_DURATION/1000.)
    sim_label = 'OBS: %s ks' % (SIM_DURATION/1000.)
    mod_label = 'Input model'
    lc_label = 'Light curve'
    _phase, _phase_err, _pol_deg, _pol_deg_err, _pol_angle,\
        _pol_angle_err, _index, _index_err, _norm,\
        _norm_err = numpy.loadtxt(ANALYSIS_FILE_PATH, unpack=True)
    plt.figure('Polarization degree')
    pl_normalization_spline.plot(scale=0.12, show=False, color='lightgray',
                                 label=lc_label)

    _pol_deg_err_p=numpy.where((_pol_deg+_pol_deg_err)<1.0,_pol_deg_err,1.0-_pol_deg)
    _pol_deg_err_m=numpy.where((_pol_deg-_pol_deg_err)>0.0,_pol_deg_err,_pol_deg)
    #if (_pol_deg+_pol_deg_err)>1.0: _pol_deg_err=1.0-yerr_p
    
    plt.errorbar(_phase, _pol_deg, xerr=_phase_err, yerr=[_pol_deg_err_m,_pol_deg_err_p], fmt='o',
                 label=sim_label)
    pol_degree_spline.plot(show=False, label=mod_label)
    plt.axis([0., 1., 0., 0.5])
    plt.legend(bbox_to_anchor=(0.45, 0.95))
    if save:
        save_current_figure('crab_polarization_degree', OUTPUT_FOLDER, False)
    plt.figure('Polarization angle')
    pl_normalization_spline.plot(scale=0.4, offset=1.25, show=False,
                                 color='lightgray', label=lc_label)
    plt.errorbar(_phase, _pol_angle, xerr=_phase_err, yerr=_pol_angle_err,
                 fmt='o', label=sim_label)
    pol_angle_spline.plot(show=False, label=mod_label)
    plt.axis([0., 1., 1.25, 3.])
    plt.legend(bbox_to_anchor=(0.45, 0.95))
    if save:
        save_current_figure('crab_polarization_angle', OUTPUT_FOLDER, False)
    plt.figure('PL normalization')
    plt.errorbar(_phase, _norm, xerr=_phase_err, yerr=_norm_err, fmt='o',
                 label=sim_label)
    pl_normalization_spline.plot(show=False, label=mod_label)
    plt.axis([0., 1., None, None])
    plt.legend(bbox_to_anchor=(0.45, 0.95))
    if save:
        save_current_figure('crab_pl_norm', OUTPUT_FOLDER, False)
    plt.figure('PL index')
    pl_normalization_spline.plot(scale=0.18, offset=1.3, show=False,
                                 color='lightgray', label=lc_label)
    plt.errorbar(_phase, _index, xerr=_phase_err, yerr=_index_err, fmt='o',
                 label=sim_label)
    pl_index_spline.plot(show=False, label=mod_label)
    plt.axis([0., 1., 1.3, 2.1])
    plt.legend(bbox_to_anchor=(0.45, 0.95))
    if save:
        save_current_figure('crab_pl_index', OUTPUT_FOLDER, False)
    plt.show()

def plot(save=False):
    """Plot the stuff in the analysis file.
    """
    for j in range(len(E_BINNING)):
        if (j==len(E_BINNING)-1):
            analysis_file = ANALYSIS_FILE_PATH
            emin=E_BINNING[0]
            emax=E_BINNING[-1]
        else:
            analysis_file = '%s_%d' % (ANALYSIS_FILE_PATH,j)
            emin=E_BINNING[j]
            emax=E_BINNING[j+1]
            
        sim_label = 'XIPE %s ks' % (SIM_DURATION/1000.)
        mod_label = 'Input model'
        lc_label = 'Light curve'
        _phase, _phase_err, _pol_deg, _pol_deg_err, _pol_angle,\
            _pol_angle_err, _index, _index_err, _norm,\
            _norm_err = numpy.loadtxt(ANALYSIS_FILE_PATH, unpack=True)
        plt.figure('Polarization degree (%g-%g keV)' % (emin,emax))
        plt.title('Her X-1 (%g-%g keV)' % (emin,emax))
        pl_normalization_spline.plot(scale=20, show=False, color='lightgray',
                                     label=lc_label)
        plt.errorbar(_phase, _pol_deg, xerr=_phase_err, yerr=_pol_deg_err, fmt='o',
                     label=sim_label)
        pol_degree_spline.plot(show=False, label=mod_label)
#        pol_degree_spline_noqed.plot(show=False, label='No QED')
        plt.axis([0., 1., 0., 1.1])
        plt.legend(bbox_to_anchor=(0.4, 0.5))
        if save:
            save_current_figure('herx1_polarization_degree_%d' % j, OUTPUT_FOLDER, False)
        plt.figure('Polarization angle (%g-%g keV)' % (emin,emax))
        plt.title('Her X-1 (%g-%g keV)' % (emin,emax))
        pl_normalization_spline.plot(scale=60, offset=0, show=False,
                                     color='lightgray', label=lc_label)
        plt.errorbar(_phase, _pol_angle, xerr=_phase_err, yerr=_pol_angle_err,
            fmt='o', label=sim_label)
        pol_angle_spline.plot(show=False, label=mod_label)
        plt.axis([0., 1., 0, 3.15])
        plt.legend(bbox_to_anchor=(0.4, 0.3))
        if save:
            save_current_figure('herx1_polarization_angle_%d' %j, OUTPUT_FOLDER, False)
        plt.figure('Flux (%g-%g keV)' % (emin,emax))
        plt.title('Her X-1 (%g-%g keV)' % (emin,emax))
        plt.errorbar(_phase,
                     0.21*_norm/(2-_index)*(10.**(2.0-_index)-2.**(2.-_index)), xerr=_phase_err,
                     yerr=0.21*_norm_err/(2-_index)*(10.**(2-_index)-2.**(2.-_index)), fmt='o',
                     label=sim_label)
        pl_normalization_spline.plot(scale=1,show=False, label=mod_label)
        plt.axis([0., 1., None, None])
        plt.legend(bbox_to_anchor=(0.75, 0.95))
        if save:
            save_current_figure('herx1_flux_%d' % j, OUTPUT_FOLDER, False)
        plt.show()

 def plot_ymap(self, overlay=True, show=False):
     """Plot the y polarization map.
     """
     if self.y_img is None:
         self.y_img = xFITSImage(self.ymap_file_path, build_cdf=False)
     fig = self.y_img.plot(show=False, zlabel="Polarization degree (y)")
     if overlay:
         self.overlay_arrows(fig)
     if show:
         plt.show()
     return fig

def plot_doc():
    """
    """
    img = xFITSImage(os.path.join(XIMPOL_DATA, 'casa_cmap.fits'))
    fig = img.plot(show=False)
    #Option to draw the psf circle on the count map
    RAD_PSF = 11/60.
    fig.show_circles(350.769, 58.873, RAD_PSF/60., lw=2, color='white')
    fig.add_label(0.73,0.90, 'PSF', relative=True, size='x-large',
                                    color='white', horizontalalignment='left') 
    xFITSImage.add_label(fig, 'XIPE 250 ks')
    plt.show()

def main(interactive=False):
    """Test the script plotting the light curve og GRB 130427A
    """
    #If you want all the GRBs use the following line:    
    grb_list = get_all_swift_grb_names()
    #grb_list = ['GRB 130427A','GRB 050124']
    plot_swift_lc(grb_list,show=False)
    overlay_tag()
    save_current_figure('Swift_XRT_light_curves',
                        clear=False)
    if interactive:
        plt.show()

def display():
    """Display the source model.
    """
    print(ROI_MODEL)
    from ximpol.utils.matplotlib_ import pyplot as plt
    fig = plt.figure('Energy spectrum')
    spectral_model.plot(show=False, logy=True)
    fig = plt.figure('Polarization degree')
    pol_degree.plot(show=False)
    fig = plt.figure('Polarization angle')
    pol_angle.plot(show=False)
    plt.show()