#-------------------------------------------------------------------------------
# Name:		    PlateModelCalculation.py
# Purpose:      Calculate depth of depth-to-basement based on Plate Model and
#               EB08 age od the oceanic crusts.
#
# Author:		Dr. Arghya Goswami
# Contact:		arghya.goswami@gmail.com
#
# Created:		14/03/2014
#
# Copyright:	(c) Dr. Arghya Goswami 2014
#
# Licence:		This program is free software; you can redistribute it and/or
#			    modify it under the terms of the GNU General Public License,
#			    version 2, as published by the Free Software Foundation.
#
#			    This program is distributed in the hope that it will be useful,
#			    but WITHOUT ANY WARRANTY; without even the implied warranty of
#			    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#			    GNU General Public License for more details.
#
#			    You should have received a copy of the GNU General Public License
#			    along with this program; if not, write to:
#
#			    Free Software Foundation, Inc.,
#			    51 Franklin St, Fifth Floor,
#			    Boston, MA, 02110-1301  USA
#-------------------------------------------------------------------------------
#!/usr/bin/env python
#-------------------------------------------------------------------------------
import arcpy, sets, numpy, sys, os, math, shutil, glob, time, fileinput, string, Tkinter, tkFileDialog
from Tkinter import *
##from arcpy.sa import *
if arcpy.CheckExtension("3D") == "Available":
	arcpy.CheckOutExtension("3D")
if arcpy.CheckExtension("GeoStats") == "Available":
	arcpy.CheckOutExtension("GeoStats")
if arcpy.CheckExtension("Spatial") == "Available":
	arcpy.CheckOutExtension("Spatial")
"""-----------------------------------------------------------------------------
# Start Time
-----------------------------------------------------------------------------"""
StartTime = time.asctime(time.localtime(time.time()))
print "Start Time: ", StartTime,"\n"
"""-----------------------------------------------------------------------------
# Check Script name
-----------------------------------------------------------------------------"""
scriptname = os.path.basename(__file__)
if "module" in scriptname:
	print "Please SAVE the script with a Meaningful NAME and RE-RUN!\n"
"""-----------------------------------------------------------------------------
# Definitions used in this Script
-----------------------------------------------------------------------------"""
# Browse a file function
def BrowseFile(msg1,indir,ext):
	#Creating Main Window
	root = Tk()
	if ext == "csv":
		desc1 = "Comma Seperated Value Files"
		desc2 = "*." + ext
	elif ext == "txt":
		desc1 = "Text Files"
		desc2 = "*." + ext
	elif ext == "shp":
		desc1 = "Shape Files"
		desc2 = "*." + ext
	elif ext == "cdf":
		desc1 = "Common Data Format Files"
		desc2 = "*." + ext
	# Creating Secondary Window for File Browsing
	filename = tkFileDialog.askopenfilename(filetypes = [(desc1,desc2),('All files','*')], title = msg1, initialdir = indir)
	# Close the Main Window
	root.withdraw()
	# Return Filename
	return filename

# Browse a directory function
def BrowseDir(msg2,indir):
	#Creating Main Window
	root = Tk()
	# Creating Secondary Window for Folder Browsing
	dirname = tkFileDialog.askdirectory(title = msg2, initialdir=indir)
	# Close the Main Window
	root.withdraw()
	# Return Dirname
	return dirname

# Check the existance of a directory function
def DirCheck(filename):
	# Check if the string is a file name or directory name
	if filename[len(filename)-4] == ".":
		dname = os.path.dirname(filename)
	else:
		dname = filename
	# Check if Dir exists, else create
	if os.path.exists(dname):
		pass
	else:
		os.mkdir(dname)
	# Check if the path name if Dir ends with a '/' else add it
	if dname[(len(dname)-1):]=="\\" or dname[(len(dname)-1):] == "/":
		dpath = dname
	else:
		dpath = dname + "/"
	return dpath

def AllfFields(infile):
	global FieldNnames
	fields = arcpy.ListFields(infile)
	FieldNnames = []
	for field in fields:
		FieldNnames.append(field.name)
	return FieldNnames

# Create the list of files in a folder
def getFileList(root,ext):
	# for folder path either use "Browse"/"browse" or provide the path to a folder
	# for all filetypes send extention as an empty string (""), else use like ".csv"
	global filelist
	if root == "Browse" or root == "browse":
		root = BrowseDir("Please select the Directory", "C:/")
		# select the type of file, for instance *.jpg or all files *.*
		for folder in glob.glob(root):
			print "Files in", folder, "is being populated to a list"
			if ext == "":
				ext = "/*.*"
			else:
				ext = "/*" + ext
		for file in glob.glob(folder + ext):
			# retrieves the stats for the current file as a tuple
			# (mode, ino, dev, nlink, uid, gid, size, atime, mtime, ctime)
			# the tuple element mtime at index 8 is the last-modified-date
			filelist.append(file)
	print "\nFile List Prepared, starting to list them......\n"
def PMDepth(MAge):
	RidgeDepth = -2639.796587
	DepthDiff = 0.00
	EqDepth = -5875.0
	Age = float(MAge)
##	Age = 57.528880999999998
	sumTerm = 0.0
	m=0.0
	while m<26:
	##    print math.pi
	##    print math.sqrt(math.pi)
	    KappaByYsq = ((-329.4618772*math.sqrt(math.pi))/(2*EqDepth))**2
	##    print KappaByYsq
	    expTerm = pow(math.e,(-KappaByYsq*pow((1+2*m),2)*pow(math.pi,2)*Age))
	##    print expTerm
	##    print "2m: ", 2*m
	##    print "1+2m: ", (1+2*m)
	##    print "(1+2m)^2: ", pow((1+2*m),2)
	##    print "1/((1+2m)^2): ", 1/pow((1+2*m),2)
	    Term = (1/(pow((1+2*m),2)))*expTerm
	    sumTerm = sumTerm + Term
##	    print "m: ",m, "Term: ", Term
	    m = m + 1
	SigmaTerm = (4/pow(math.pi,2))*sumTerm
##	print SigmaTerm
	DepthDiff = EqDepth*(0.5 - SigmaTerm)
##	print DepthDiff
	RajD2B = RidgeDepth + DepthDiff
##	print RajD2B
"""-----------------------------------------------------------------------------
# Main Script
-----------------------------------------------------------------------------"""
fIn = BrowseFile("Open Muller Age File","C:\\","csv")
fOut = fIn.replace("MAge","RajPM")

f1 = open(fIn,'r')
ftext = f1.readlines()
f1.close()

f2 = open(fOut,'w')
for lines in ftext:
	line = lines.replace("\n","")
	ages = lines.split(",")
	PM = []
	for items in ages:
		if items == "nan":
			pm = "nan"
		else:
			pm = PMDepth(items)
		PM.append(pm)
	newline = ""
	for things in PM:
		if newline == "":
			newline = str(things)
		else:
			newline = newline + "," + str(things)
	newline = newline + "\n"
	f2.write(newline)
f2.close()
"""-----------------------------------------------------------------------------
# End Time
-----------------------------------------------------------------------------"""
EndTime = time.asctime(time.localtime(time.time()))
print "\nEnd Time: ", EndTime,"\n"