# PlaneTrendSurf: program that calculates an ASCII z-ccordinate ESRI raster for a known planar
# orientation.
# Input: ASCII text file "plane_in.txt".
# Line 1. trend and plunge degrees of the true dip vector (trend in azimuth format). (ex. 135.0, 1.43)
# Line 2. Elevation (Z value) of the measured attitude. (ex. 191.0)
# Line 3. "X,Y" coordinates of the position of the attitude measurement (ex. "282883.069 , 251535.329").
# Line 4. SW corner (lower left) "X,Y" coordinates of the calculated raster image.
# Line 5. "Columns,Rows" in calulated raster (ex. "251, 252").
# Line 6. Grid spacing between rows and columns of raster grid (ex. "21.6").
# ----------------------------------------------------------------------------------------------       
# Output: ESRI raster file with following structure:
# line 1: "ncols=" {number of columns in raster grid}
# line 2: "nrows=" {number of rows in raster grid}
# line 3: "xllcorner=" {x coordinate of the lower left (SW) corner of the grid}
# line 4: "yllcorner=" {y coordinate of the lower left (SW) corner of the grid}
# line 5: "cellsize=" {cell size spacing between columns and rows}
# line 6: "nodata value=" {value that represents no data at grid node}
#-----------------------------------------------------------------------------------------------

import sys, numpy, math
#import win32gui    #enable WIN32 MessageBox dialog
#import win32ui
#import win32con


# Function deg2rad(d): converts "d" from degrees to radians
def deg2rad(d):
    return d * math.pi / 180.0
#Function rad2deg(r): converts "r" from radians to degrees
def rad2deg(r):
    return r * 180.0 / math.pi

# Function dprod(x,y,z,a,b,c): returns the dot-product of the directional components (x,y,z) and
# (a,b,c).

def dprod(x,y,z,a,b,c):
   return x*a + y*b + z*c
# ------------------ Trend Surface Program ---------------------------------------------------

# initialize constants
tolerance = 1.0E-9
dl = " "
eoln = "\n"
# subdirectory path for output
path = "C:/QGIS_Data/OutcropPrediction/Prob2/"
# file name for output
ofn = "Prob2_bot_grd.txt"
rfn = "Prob2_bot_in.txt"
# open trend output file in write-mode
of = open(path+ofn,"w")
rf = open(path+rfn,"r")

# Enter azimuth, plunge of true dip:
att_str = rf.readline()
att_list = att_str.split(",")
dipazstr = att_list[0]
dipaz = float(dipazstr)
plungestr = att_list[1]
plunge = float(plungestr)
# Calculate the pole to the plane directional components, and the directional components to the
# X-Z plane
if dipaz < 180.0 :
    pole_az = dipaz + 180.0
else:
    pole_az = dipaz - 180.0
pole_pl = 90.0 - plunge
# directional components to pole to plane = (x1,y1,z1)
x1 = math.sin(deg2rad(pole_az)) * math.sin(deg2rad(90.0-pole_pl))
y1 = math.cos(deg2rad(pole_az)) * math.sin(deg2rad(90.0-pole_pl))
z1 = math.cos(deg2rad(90.0-pole_pl))
# directional components of pole to X-Z plane = (x2,y2,z2)
x2 = 0.0
y2 = 1.0
z2 = 0.0

# dot product yields theta angle between pole to plane and X-Z plane pole
theta1 = math.acos(dprod(x1,y1,z1,x2,y2,z2))
if abs(theta1 - deg2rad(90.0)) < tolerance :
    theta1 = deg2rad(90.0)
# cross product of pole to plane and X-Z plane pole yields the line of intersection
# between the 2 planes (x3,y3,z3), and this line is the slope of the plane in the X-Z plane
c1_x = (y1*z2-z1*y2)/math.sin(theta1)
c1_y = -(x1*z2-z1*x2)/math.sin(theta1)    
c1_z = (x1*y2-y1*x2)/math.sin(theta1)
    
if c1_z < 0.0 :
    c1_x = c1_x * -1.0
    c1_y = c1_y * -1.0
    c1_z = c1_z * -1.0
# check for round off errors    
if c1_x > 1.0 :
    c1_x = 1.0
if c1_x < -1.0 :
    c1_x = -1.0
c1 = -math.tan(math.acos(c1_x))

# directional components of pole to Y-Z plane = (x2,y2,z2)
x2 = -1.0
y2 = 0.0
z2 = 0.0
# dot product yields theta angle between pole to plane and Y-Z plane pole
theta2 = math.acos(dprod(x1,y1,z1,x2,y2,z2))
# cross product of pole to plane and Y-Z plane pole yields the line of intersection
# between the 2 planes (x3,y3,z3), and this line is the slope of the plane in the Y-Z plane
c2_x = (y1*z2-z1*y2)/math.sin(theta2)
c2_y = -(x1*z2-z1*x2)/math.sin(theta2)
c2_z = (x1*y2-y1*x2)/math.sin(theta2)
if c2_z < 0.0 :
    c2_x = c2_x * -1.0
    c2_y = c2_y * -1.0
    c2_z = c2_z * -1.0
# check for round off errors
if c2_y > 1.0 :
    c2_y = 1.0
if c2_y < -1.0 :
    c2_y = -1.0
 
c2 = -math.tan(math.acos(c2_y))

# the elevation (z) value of the measured attitude:
elev_str = rf.readline()
elev = float(elev_str)
# the (x,y) coordinates of the measured attitude:
xy_str = rf.readline()
crd_list = xy_str.split(",")
x_crd = float(crd_list[0])
y_crd = float(crd_list[1])
c0=elev-c1*x_crd-c2*y_crd
# (x,y) coordinates of grid SW corner:
sw_xy_str = rf.readline()
crd_list = sw_xy_str.split(",")
sw_x_crd = float(crd_list[0])
sw_y_crd = float(crd_list[1])
# number of columns,rows for grid:
grid_str = rf.readline()
grid_list = grid_str.split(",")
ncols_str = grid_list[0]
ncols = int(ncols_str)
nrows_str = grid_list[1]
nrows = int(nrows_str)
# cell size of raster grid:
cellsz_str = rf.readline()
cellsize = float(cellsz_str)

nodatavalue = 1e+30
# write header info
of.write("ncols" + dl + str(ncols) + eoln)
of.write("nrows" + dl + str(nrows) + eoln)
of.write("xllcorner" + dl + str(sw_x_crd) + eoln)
of.write("yllcorner" + dl + str(sw_y_crd) + eoln)
of.write("cellsize" + dl + str(cellsize) + eoln)
of.write("nodata_value" + dl + str(nodatavalue) + eoln)
linect = 6
i = 0
while i < nrows:
  y = sw_y_crd + (nrows - i) * cellsize
  j = 0
  while j < ncols:
      x = sw_x_crd + (j * cellsize)
      z = c0 + c1*x + c2*y  
      of.write(str(z)+eoln)
      linect = linect + 1
      j = j + 1
  i = i + 1
titlestr = "Results Message" + eoln
print titlestr
msgstr = "Planar surface raster written to: " + path + ofn + eoln
print msgstr
# win32ui.MessageBox(msgstr,titlestr,win32con.MB_OK)
titlestr = "z = c0 + c1*x + c2*y" + eoln
print titlestr
msgstr = "c0=" + str(c0) + eoln + "c1=" + str(c1) + eoln + "c2=" + str(c2) + eoln
print msgstr
# win32ui.MessageBox(msgstr,titlestr,win32con.MB_OK)

# close the files
of.close()
rf.close()