import math

def tempModel(depth):
    """The equation below came from a fit to temperature versus depth data provided by kurt w."""
    return 1.83415e-09*math.pow(depth,3) + (-1.59061e-08*math.pow(depth,2)) + 0.00267687*depth + (-51.0696 )


def densityModel(temp):
    """The equation below came from fitting a curve through some density versus temperature data in Ice(h) from the CRC handbook"""
    return 0.916899+0.0227655*(1-math.exp(0.00761904*temp))


def arconne(density):
    """This formula came from bob morse.  Given a density it should return an index of refraction"""

    index = 1.0+0.845*density

    return index


if __name__=="__main__":
    fd = open('calculated-index.dat', 'w')
    fd.write('# depth ( in meters )  density  calculated index of refraction\n')

    for depth in range(0,2500):
        temp = tempModel(depth)
        density = densityModel(temp)
        index = arconne(density)

        fd.write('%f\t%f\t%f\n' % ( depth, density, index))
    fd.close()