#!/usr/bin/python


import numpy
import readData
import pyfits
import math
from scipy import interpolate


print "****************************************************"
print "*                                                  *"
print "*  Convert the TIGRE spectrum into an Fits file    *"
print "*                                                  *"
print "*                readable in IRAF                  *"
print "*                                                  *"
print "****************************************************"
print ""

print "TIGRE Spektrum: ",
file = raw_input()
print "Output file: ", 
output_file = raw_input()


spec = readData.readheros(file)

him = spec[0]
hs = spec[1]
fields = hs['TFIELDS']
object = him['OBJECT']
exptime = him['EXPTIME']
if fields == 12:
    wave = spec[2]
    spec_blz = spec[3]
    blaze = spec[4]
    sig = spec[5]
    spec_norm = spec[6]
    wave_merg = spec[7]
    spec_merg = spec[8]
    f_spec_merg = spec[9]
    spec_merg_norm = spec[10]
    f_spec_merg_norm = spec[11]
    global_res_curve= spec[12]
    day_var = spec[13]
if fields == 13:
    wave = spec[2]
    spec_blz = spec[3]
    blaze = spec[4]
    sig = spec[5]
    spec_norm = spec[6]
    single_spec = spec[7]
    wave_merg = spec[8]
    spec_merg = spec[9]
    f_spec_merg = spec[10]
    spec_merg_norm = spec[11]
    f_spec_merg_norm = spec[12]
    global_res_curve= spec[13]
    day_var = spec[14]


count = 0
while (count < 3):
    
    print ""
    print "If you want a single order (y/n): ",
    question_order = raw_input()
    
    
    if question_order == 'y':
        print "Select order: ",
        order_num = raw_input()
        
        wave_order = wave[::,order_num] 
        spec_order = spec_blz[::,order_num]
        sig_order = sig[::,order_num]
        blaze_order = blaze[::,order_num]
        nwave = len(wave_order)
        
        w0 = numpy.ceil(wave_order[0]*100.)/100.
        w_max = numpy.floor(wave_order[nwave-1]*100.)/100.
        
        dwave = numpy.around(((w_max - w0) / nwave)*100.)/100.
        
        npix_new_grid = numpy.floor( (w_max - w0) / dwave)
        
        new_grid = w0 + numpy.arange(npix_new_grid)*dwave
        data_dummy = interpolate.InterpolatedUnivariateSpline(wave_order,spec_order)
        err_data_dummy = interpolate.InterpolatedUnivariateSpline(wave_order,sig_order)
        blz_dummy = interpolate.InterpolatedUnivariateSpline(wave_order,blaze_order)
        
        data = data_dummy(new_grid)
        err_data = err_data_dummy(new_grid)
        blaze = blz_dummy(new_grid)
        
        nwave = len(new_grid)
        tab = numpy.zeros((3,nwave), float)
        
        tab[0,::] = data
        tab[1,::] = err_data
        tab[2,::] = blaze
        
        hdu = pyfits.PrimaryHDU(tab)
        hdu.header.update('CRVAL1', w0 ,'Coordinate at reference pixel')
        hdu.header.update('CDELT1', dwave, 'Corrdinate increment per pix')
        hdu.header.update('CTYPE1', 'WAVELENGTH', 'Units of corrdinate')
        hdu.header.update('CRVAL2', w0 ,'Coordinate at reference pixel')
        hdu.header.update('CDELT2', dwave, 'Corrdinate increment per pix')
        hdu.header.update('CTYPE2', 'WAVELENGTH', 'Units of corrdinate')
        hdu.header.update('CRVAL3', w0 ,'Coordinate at reference pixel')
        hdu.header.update('CDELT3', dwave, 'Corrdinate increment per pix')
        hdu.header.update('CTYPE3', 'WAVELENGTH', 'Units of corrdinate')
        
        hdu.header.update('OBJECT', object, 'TARGET NAME')
        hdu.header.update('EXPTIME', exptime, 'EXPOSURE TIME (seconds)')
        
        hdu.writeto(output_file)
        
        count=4
            
    elif question_order == 'n':
        
        wave_order = wave_merg
        spec_order = spec_merg
        sig_order = f_spec_merg
        spec_order_norm = spec_merg_norm
        sig_order_norm = f_spec_merg_norm

        glolal_corr = global_res_curve
        small_scale_corr = day_var
        cn = 0
        while (cn < 3):
            print "With global correction (y/n): ",
            global_corr_question = raw_input()

            if global_corr_question == 'y':
                spec_order = spec_order*glolal_corr
                sig_order = sig_order*glolal_corr
                cn=4
            elif global_corr_question == 'n':
                spec_order = spec_order
                sig_order = sig_order
                cn=4
            else:
                print ""
                print "You have to use y or n!"
            
                cn = cn+1
                if cn == 4:
                    print "No correction is performed"

        cn2 = 0
        while (cn2 < 3):
            print "With small scale correction (y/n): ",
            small_scale_corr_question = raw_input()

            if small_scale_corr_question == 'y':
                spec_order = spec_order*small_scale_corr
                sig_order = sig_order*small_scale_corr
                cn2=4
            elif small_scale_corr_question == 'n':
                spec_order = spec_order
                sig_order = sig_order
                cn2=4
            else:
                print ""
                print "You have to use y or n!"
            
                cn2 = cn2+1
                if cn2 == 4:
                    print "No correction is performed"

        nwave = len(wave_order)
        
        w0 = numpy.ceil(wave_order[0]*100.)/100.
        w_max = numpy.floor(wave_order[nwave-1]*100.)/100.
        
        dwave = numpy.around(((w_max - w0) / nwave)*100.)/100.
        
        npix_new_grid = numpy.floor( (w_max - w0) / dwave)
            
        new_grid = w0 + numpy.arange(npix_new_grid)*dwave
        
        data_dummy = interpolate.InterpolatedUnivariateSpline(wave_order,spec_order)
        err_data_dummy = interpolate.InterpolatedUnivariateSpline(wave_order,sig_order)
        data_norm_dummy = interpolate.InterpolatedUnivariateSpline(wave_order,spec_order_norm)
        err_data_norm_dummy = interpolate.InterpolatedUnivariateSpline(wave_order,sig_order_norm)
            
        data = data_dummy(new_grid)
        err_data = err_data_dummy(new_grid)
        data_norm = data_norm_dummy(new_grid)
        err_data_norm = err_data_norm_dummy(new_grid)
        
        nwave = len(new_grid)
        tab = numpy.zeros((4,nwave), float)
            
        tab[0,::] = data
        tab[1,::] = err_data
        tab[2,::] = data_norm
        tab[3,::] = err_data_norm
            
        hdu = pyfits.PrimaryHDU(tab)
        hdu.header.update('CRVAL1', w0 ,'Coordinate at reference pixel')
        hdu.header.update('CDELT1', dwave, 'Corrdinate increment per pix')
        hdu.header.update('CTYPE1', 'WAVELENGTH', 'Units of corrdinate')
        hdu.header.update('CRVAL2', w0 ,'Coordinate at reference pixel')
        hdu.header.update('CDELT2', dwave, 'Corrdinate increment per pix')
        hdu.header.update('CTYPE2', 'WAVELENGTH', 'Units of corrdinate')
        hdu.header.update('CRVAL3', w0 ,'Coordinate at reference pixel')
        hdu.header.update('CDELT3', dwave, 'Corrdinate increment per pix')
        hdu.header.update('CTYPE3', 'WAVELENGTH', 'Units of corrdinate')
        hdu.header.update('CRVAL4', w0 ,'Coordinate at reference pixel')
        hdu.header.update('CDELT4', dwave, 'Corrdinate increment per pix')
        hdu.header.update('CTYPE4', 'WAVELENGTH', 'Units of corrdinate')
        
        hdu.header.update('OBJECT', object, 'TARGET NAME')
        hdu.header.update('EXPTIME', exptime, 'EXPOSURE TIME (seconds)')
            
        hdu.writeto(output_file)
            
        count=4
        
    else:
        print ""
        print "You have to use y or n!"
        print ""

        count = count+1