import math, sys

from math import sin,cos,floor,pi,sqrt
from Wizardry.AIPSData import AIPSUVData

c_vacuum = 299792458.0

def create_scint_sn_table(sourcename, aipsclass, diskno, seqno, experiment, uv_aipsclass, uvdatatable_number, scinttime, exclude):
    pulsardata = AIPSUVData(sourcename, aipsclass, diskno, seqno)
    uvdata = AIPSUVData(experiment.upper(), uv_aipsclass, 1, 1)
    samplesntable = uvdata.table('SN', uvdatatable_number)
    keywords = samplesntable.keywords
    num_snterms = 8
    #num_snterms = 0
    #for term in keywords:
    #    num_snterms = num_snterms + 1
    samplerow = samplesntable[0]
    
    #print pulsardata.header
    snpoints = 0
    lasttime = -99999
    sntable = pulsardata.attach_table('SN', 1, no_term=num_snterms)
    sntable.keywords['NO_ANT'] = len(pulsardata.antennas)
    sntable.keywords['NO_IF'] = 4
    sntable.keywords['NO_NODES'] = samplesntable.keywords['NO_NODES']
    sntable.keywords['NO_POL'] = 2
    num_freqs = 4
    num_stokes = 2
    
    print samplerow
    gainlength = len(samplerow.real1)
    
    #set all the stuff in sample row to zero, except the time, antenna, IF and gain (real and weight)
    samplerow.i_far_rot = 0
    samplerow.time_interval = scinttime/1440.0
    samplerow.mbdelay2 = 0
    samplerow.source_id = 1
    samplerow.mbdelay1 = 0
    samplerow.node_no = 0
    samplerow.source_id = 1
    samplerow.freq_id = 1
    for i in range(gainlength):
        samplerow.imag1[i] = 0.0
        samplerow.imag2[i] = 0.0
        samplerow.delay_1[i] = 0.0
        samplerow.delay_2[i] = 0.0
        samplerow.rate_1[i] = 0.0
        samplerow.rate_2[i] = 0.0

    count = [0]*gainlength
    total = [0.0]*gainlength
    outtimes = []
    times = [None]*gainlength
    amps = [None]*gainlength
    outgains = [None]*gainlength
    for i in range(gainlength):
        times[i] = []
        amps[i] = []
        outgains[i] = []
    
    #Loop through the visibilities, averaging and working out amplitudes
    for visibility in pulsardata:
        skip = False
        for t in exclude:
            if visibility.baseline[0] == t or visibility.baseline[1] == t:
                skip = True
        if not skip:
            for i in range(gainlength):
                for j in range(2):
                    amp = math.sqrt(visibility.visibility[i][0][j][0]*visibility.visibility[i][0][j][0] + visibility.visibility[i][0][j][1]*visibility.visibility[i][0][j][1])
                    if amp > 0 and amp == amp:
                        total[i] = total[i] + amp
                        count[i] = count[i] + 1
                        times[i].append(visibility.time)
                        amps[i].append(amp)

    #Do a Hanning smooth, FWHM of scintillation time
    #time_interval = scinttime/1440.0

    #No Hanning smooth
    time_interval = scinttime/(1440.0*2.0) #This will ensure only one in range
    starttime = times[0][0]
    endtime = times[0][count[0]-1]
    avoidindices = []
    for i in range(1,gainlength):
        if times[i][0] < starttime:
            starttime = times[i][0]
        if times[i][count[i]-1] > endtime:
            endtime = times[i][count[i]-1]
    for i in range(gainlength):
        attime = starttime
        startat = 0
        averageamp = total[i]/float(count[i])
        while attime < endtime:
            num = 0
            denom = 0
            for j in range(startat, count[i]):
                if not times[i][j] < attime + time_interval:
                    break
                if times[i][j] > attime - time_interval:
                    scale = 0.5*(1 - math.cos(2*math.pi*(0.5*(times[i][j]-attime)/time_interval + 0.5)))
                    num = num + scale*amps[i][j]
                    denom = denom + scale
                else:
                    startat = i
            if i==0:
                outtimes.append(attime)
            if denom > 0.0:
                #Edited to correct scaling
                outgains[i].append(math.sqrt(averageamp/(num/denom)))
            else:
                outgains[i].append(0.0)
            attime = attime + time_interval

    for i in range(len(outtimes)):
        present = False
        for j in range(gainlength):
            if outgains[j][i] > 0.0:
                present = True
            else:
                outgains[j][i] = 1.0
        if not present:
            avoidindices.append(i)

    #Create the table
    avoidat = 0
    for i in range(len(outtimes)):
        if avoidat >= len(avoidindices) or (not i == avoidindices[avoidat]):
            for j in range(gainlength):
                samplerow.real1[j] = outgains[j][i]
                samplerow.real2[j] = outgains[j][i]
                samplerow.weight_1[j] = 1/(outgains[j][i]*outgains[j][i])
                samplerow.weight_2[j] = 1/(outgains[j][i]*outgains[j][i])
                samplerow.time = outtimes[i]
            for j in range(len(pulsardata.antennas)):
                samplerow.antenna_no = j + 1
                sntable.append(samplerow)
        else:
            avoidat = avoidat+1

    #Work out all the points
    #for visibility in pulsardata:
    #    skip = False
    #    for t in exclude:
    #        if visibility.baseline[0] == t or visibility.baseline[1] == t:
    #            skip = True
    #    if not skip:
    #        for i in range(gainlength):
    #            for j in range(2):
    #                if visibility.visibility[i][0][j][0] != 0.0 and visibility.visibility[i][0][j][1] != 0.0:
    #                    amp = total/(float(count)*math.sqrt(visibility.visibility[i][0][j][0]*visibility.visibility[i][0][j][0] + visibility.visibility[i][0][j][1]*visibility.visibility[i][0][j][1]))
    #                    if amp > 0 and amp == amp:
    #                        for k in range(gainlength):
    #                            samplerow.real1[k] = amp
    #                            samplerow.real2[k] = amp
    #                            samplerow.weight_1[k] = 1/(amp*amp)
    #                            samplerow.weight_2[k] = 1/(amp*amp)
    #                        samplerow.time = visibility.time
    #                        for k in range(len(pulsardata.antennas)):
    #                            samplerow.antenna_no = k + 1
    #                            sntable.append(samplerow)

    sntable.close()