import math
import sys
import glob
import os
import pexpect
import hull
import psyco
import fileManager

psyco.full()

def getAvgAlpha(depth):
    """Note that the depth is positive in meters
    This uses a polynomial fit to the change in the RF loss coefficient
    at 300 MHz versus depth for the data from bogorodsky's book"""

    alpha0_300 = 6.6e-4
    coeff = 6.24e-17
    dD = depth / 100.0
    losstot = 1.0

    # note for range the first argument is incluse
    # and the second is not
    # so range(10,1, -1)
    # gives 10,9,8,7,6,5,4,3,2
    depth1 = depth
    while depth1>=0:
        losstot = losstot * math.exp(-(alpha0_300 + coeff * math.pow(depth1, 4.0))*dD)
        depth1 = depth1 - dD
        
    avgloss = -1.0*math.log(losstot)/depth

    return avgloss


def getDndz(depth):

    if(depth>=1.0):
        K = 0.438
        alpha = 0.0132

        dndz = K * alpha * math.exp(-1.0*alpha*depth)
    else:
        dndz = 0.35

    return dndz

def getIndexOfRefraction(depth):

    if(depth>=1.0):
        n0 = 1.35
        K = 0.438
        alpha = 0.0132

        idx = n0 + K * ( 1.0 - math.exp(-1.0*alpha*depth))
    else:
        idx = 1.0 + 0.35 * depth

    return idx


def rayTrace(dir,initialDepth, initialTheta,
             thermalNoiseFloor, showerEnergy,
             highFreq, lowFreq,
             showFirnExit=False, numBins=19):

    print "Starting ray trace depth = %f\ttheta=%f\n" % ( initialDepth, initialTheta)

    totalPathLength = 0
    notAttenuated = True
    depth = initialDepth
    xpos = 0

    # shower energy passed to this routine is log(10)(eV)
    # convert it to plain eV
    showerEnergy = math.pow(10.0, showerEnergy)

    # convert to TeV
    showerEnergy = showerEnergy * 10**-12

    # thermal noise floor was in microV/m now in V/m
    thermalNoiseFloor = thermalNoiseFloor * 10**-6
            
    print "shower energy: %f" % showerEnergy

    binWidth = (highFreq-lowFreq)/numBins
    eField0 = []
    intensity = []

    for bin in range(numBins):
        binLowFreq = bin*binWidth + lowFreq
        binHighFreq = (bin+1)*binWidth + lowFreq

        deltaF = binHighFreq-binLowFreq
        barF = (binHighFreq+binLowFreq)/2.0
    
        # this is in volts per meter
        e0 = (1.1*10.0**-7)*(showerEnergy)*(barF/500.0)*(1.0/(1+0.4*math.pow(barF/500.0, 2.0)))*deltaF

        eField0.append(e0)
        intensity.append(1.0)

    print "SUM EFIELD: %f" % sum(eField0)
    
    totalEField = sum(eField0)
    if(totalEField<thermalNoiseFloor):
        print "initial ray energy too small, not used"
        print "Electric Field: %f" % totalEField
        print "Thermal Noise Floor: %f" % thermalNoiseFloor
        return

    theta = initialTheta

    outFd = open('%s/ray-%04.1f-%04.1f.dat' % ( dir, initialDepth, initialTheta ), 'a')

    n = getIndexOfRefraction(depth)
    dndz = getDndz(depth)
    
    outFd.write('%f\t%f\t%f\t%f\t%f\t%f\n' % ( depth, xpos,
                                               totalPathLength, theta,
                                               n, dndz))

    eField = []
    stepSize = 1.0
    while(depth>=0 and notAttenuated):
        # use path length steps of 0.5 meters
        x = stepSize * math.sin(math.radians(theta))
        z = stepSize * math.cos(math.radians(theta))

        depth = depth - z
        xpos = xpos + x
        totalPathLength = totalPathLength + stepSize

        eField = []
        for bin in range(numBins):
            eField.append(eField0[bin]*1.0/(totalPathLength)*intensity[bin])
            intensity[bin] = intensity[bin] * math.exp(-1.0*getAvgAlpha(depth)*stepSize)

        n = getIndexOfRefraction(depth)
        dndz = getDndz(depth)
        deltaTheta = (1.0/n) * math.sin(math.radians(theta)) * dndz * stepSize
        theta = theta + math.degrees(deltaTheta)
        
        if(sum(eField)<thermalNoiseFloor):
            print "Electric field below noise floor ( initial Theta = %f )" % initialTheta
            print "Total EField at attenuation: %f" % sum(eField)
            notAttenuated=False
            
        outFd.write('%f\t%f\t%f\t%f\t%f\t%f\n' % ( depth, xpos,
                                                   totalPathLength, theta,
                                                   n, dndz))
        if(depth<1.0):
            stepSize = 0.01
        else:
            if( (depth-2.0)<1.0):
                stepSize = 0.01
            else:
                stepSize = 2

        if(depth>=2800.0):
            # hit the bedrock..  should probably reflect but don't know how yet
            break
        
    if(depth<=0.0 and showFirnExit):
        print "Showing the exit from the firn!"
        # draw in a last section
        stepSize = 500.0
        x = stepSize * math.sin(math.radians(theta))
        z = stepSize * math.cos(math.radians(theta))
        
        depth = depth - z
        xpos = xpos + x
        totalPathLength = totalPathLength + stepSize

        outFd.write('%f\t%f\t%f\t%f\t%f\t%f\n' % ( depth, xpos,
                                                   totalPathLength, 0.0,
                                                   n, dndz))
        outFd.close()
    else:
        outFd.close()

def gnuplotGenerate(dir, pattern, postfix, startAngle, stopAngle,
                    highFreq, lowFreq, thermalFloor, showerEnergy,
                    vertFile, vol):

    files = glob.glob("%s/%s" % (dir, pattern))
    
    outFd = open('%s/plot.m' % (dir), 'w')
    initialized = False

    if(vertFile!=None and vol!=0):
        titleStr = "Version A, Depth %s meters Theta %3.2f to %3.2f degrees\\nShowerEnergy %2.2f (log10 eV) Freq %5.1f to %5.1f MHz\\nNoise Floor %3.1f uV/m Volume: %f" % (postfix, startAngle, stopAngle, showerEnergy, highFreq, lowFreq, thermalFloor, vol/(10.**9))
    else:
        titleStr = "Depth %s meters Theta %3.2f to %3.2f degrees\\nShowerEnergy %2.2f (log10 eV) Freq %5.1f to %5.1f MHz\\nNoise Floor %3.1f uV/m" % (postfix, startAngle, stopAngle, showerEnergy, highFreq, lowFreq, thermalFloor)

    if len(files)==0:
        print "<H3>Settings produced no useful raytrace!  Please edit settings and try again!</H3>"
        return -1
    
    for fName in files:
        if(not initialized):
            outFd.write("#!/usr/local/bin/gnuplot -persist\n")
            outFd.write("set output '/dev/null'\n")
            outFd.write('set nokey\n')
            outFd.write('set pointsize 0.2\n')
            outFd.write('set yrange [ 2800:-500 ]\n')
            outFd.write('set xrange [ -1000:5500 ]\n')

            #outFd.write("set title '%s'\n" % titleStr)
            outFd.write("set xlabel 'X Displacement ( Meters )'\n")
            outFd.write("set ylabel 'Depth ( Meters )'\n")
            outFd.write("plot 0 with linespoints\n")
            outFd.write("set style line 1 lt 1 linecolor rgb \"blue\"\n")
            outFd.write("replot '%s' using 2:1 with lines ls 1\n" % fName)
            initialized=True
        else:
            outFd.write("replot '%s' using 2:1 with lines ls 1\n" % fName)

    if(vertFile!=None):
        outFd.write("replot '%s/%s' using 1:2 with lines title 'Volume Area'\n" %  (dir, vertFile))
    
    outFd.write("set terminal png\n")
    outFd.write("set output '%s/plot.png'\n"% dir)
    outFd.write("set multiplot layout 1,1 title \"%s\"\n" % titleStr)
    outFd.write("replot\n")
    outFd.write("unset multiplot\n")

    outFd.write("set terminal postscript enhanced color\n")
    outFd.write("set output '%s/plot.ps'\n" % dir)
    outFd.write("set multiplot layout 1,1 title \"%s\"\n" % titleStr)
    outFd.write("replot\n")
    outFd.write("unset multiplot\n")
    
    outFd.close()

    os.chmod('%s/plot.m' % (dir), 0755)

    print "About to generate the postscript version of the plots"
    p = pexpect.spawn('%s/plot.m' % dir)
    p.expect(pexpect.EOF)
    print "Done..  About to generate the pdf file"
    
    p = pexpect.spawn('/usr/bin/ps2pdf %s/plot.ps %s/plot.pdf' % ( dir,dir ), logfile=sys.stdout)
    p.expect(pexpect.EOF)
    return 0

def generatePlot(depth, startAngle, stopAngle, angleStep,
                 thermalNoiseFloor, showerEnergy,
                 highFreq, lowFreq, showFirnExt):
    fm = fileManager.fileManager('/var/www/html/rayguiA/data', fileTimeOut=3600)
    dir = fm.getDir()
    print "generate plot.. dir = %s" % dir
    
    depths = [ depth ]

    for depth in depths:
        for theta in range(startAngle, stopAngle+angleStep, angleStep):
            rayTrace(dir, depth, theta, thermalNoiseFloor, showerEnergy,
                     highFreq, lowFreq, showFirnExit=showFirnExt)

    print "startAngle: %f, stopAngle: %f" % ( startAngle, stopAngle )
    if(startAngle==0.0 and stopAngle==180.0):
        a = hull.hull(dir, 'ray-%04.1f*.dat' % depth)
        a.load()
        verts = a.calculate()
        a.save('%s/verts.dat' % dir, verts)
        vol = a.volume(depth)
        print "VOLUME: %f" % vol
    
        returnVal = gnuplotGenerate(dir, 'ray-%04.1f*.dat' % depth, '%04.1f' % depth, startAngle, stopAngle, highFreq, lowFreq, thermalNoiseFloor, showerEnergy, 'verts.dat', vol)
    else:
        returnVal = gnuplotGenerate(dir, 'ray-%04.1f*.dat' % depth, '%04.1f' % depth, startAngle, stopAngle, highFreq, lowFreq, thermalNoiseFloor, showerEnergy, None, 0)

    if(returnVal==-1):
        return ""
    else:
        return dir