#!/usr/bin/env python 

'''
License
    This file is part of foam-extend.

    foam-extend is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation, either version 3 of the License, or (at your
    option) any later version.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Author
    Matthias Rauter matthias.rauter@uibk.ac.at
'''

# these 3 lines make python2 compatible with python3
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
from scipy.spatial import Voronoi
from scipy import interpolate
import re
import copy
import argparse
import sys

class vertex(object):
    def __init__(self, x, y, z):
        self.x=float(x); self.y=float(y); self.z=float(z);
    def __abs__(self):
        return (self.x**2.+self.y**2.+self.z**2.)**0.5
    def __add__(self, o):
        return vertex(self.x+o.x, self.y+o.y, self.z+o.z) 
    def __sub__(self, o):
        return vertex(self.x-o.x, self.y-o.y, self.z-o.z) 
    def __mul__(self, o):
        if type(self) == type(o):
            return self.x*o.x+self.y*o.y+self.z*o.z
        else:
            return vertex(self.x*o, self.y*o, self.z*o)
    def __truediv__(self, o):
        return vertex(self.x*o, self.y*o, self.z*o)
    def __pow__(self, o):
        return vertex(self.y*o.z-self.z*o.y, self.z*o.x-self.x*o.z, self.x*o.y-self.y*o.x)
    def __neg__(self, o):
        return vertex(-self.x, -self.y, -self.z)
    def __str__(self):
        return 'vertex %f %f %f'%(self.x, self.y, self.z)
            
class surf(object):
    def __init__(self, v0, v1, v2):
        self.v0 = v0; self.v1 = v1; self.v2 = v2;
        self.updateNormal()
    def updateNormal(self):
        self.n = (self.v1-self.v0)**(self.v2-self.v0)
        self.n = self.n/abs(self.n)
    def __str__(self):
        self.updateNormal()
        s = 'facet normal %f %f %f\n'%(self.n.x, self.n.y, self.n.z)
        s += 'outer loop\n'
        s += '%s\n'%(self.v0)
        s += '%s\n'%(self.v1)
        s += '%s\n'%(self.v2)
        s += 'endloop\n'
        s += 'endfacet\n'
        return s



def header(objectName, className='scalarField', locationName="constant/polyMesh"):
    return ("/*--------------------------------*- C++ -*----------------------------------*\\\n"
            "| =========                 |                                                 |\n"
            "| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n"
            "|  \\\\    /   O peration     | Version:  4.1                                   |\n"
            "|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |\n"
            "|    \\\\/     M anipulation  |                                                 |\n"
            "\*---------------------------------------------------------------------------*/\n"
            "FoamFile\n"
            "{}\n"
            "    version     2.0;\n"
            "    format      ascii;\n"
            "    class       {};\n"
            "    location    \"{}\";\n"
            "    object      {};\n"
            "{}\n"
            "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n"
            ).format("{", className, locationName, objectName, "}");

def footer():
    return ("\n"
            "// ************************************************************************* //\n"
            "\n")


def N(xis):
    return 1/4.*(1+xis[0])*(1+xis[1])

def getX(xi, x0):
    xi0 = np.array([1, -1])
    xi1 = np.array([1, 1])
    xi2 = np.array([-1, 1])
    xi3 = np.array([-1, -1])

    xg = x0[0]*N(xi*xi0)+x0[1]*N(xi*xi1)+x0[2]*N(xi*xi2)+x0[3]*N(xi*xi3)
    return xg
    
pattern = r"""
    \-*[0-9]+\.*[0-9]*
    """
rx = re.compile(pattern, re.VERBOSE)


def coords(s):
    try:
        x, y = map(float, s.split(','))
        return x, y
    except:
        raise argparse.ArgumentTypeError("Coordinates must be 'x,y'")



def main(argv):
    
    infile = ''
    outfile = ''
    xres = 100
    yres = 100

    parser = argparse.ArgumentParser(description='Taking a Raster-File and creating a STL-File')

    parser.add_argument('-i', type=str, help='input filename')
    parser.add_argument('-o', type=str, help='output filename')
    parser.add_argument('-xres', type=int, help='x resolution')
    parser.add_argument('-yres', type=int, help='y resolution')
    parser.add_argument('-p1', type=coords, help='lower left edgepoint')
    parser.add_argument('-p2', type=coords, help='lower right edgepoint')
    parser.add_argument('-p3', type=coords, help='upper right edgepoint')
    parser.add_argument('-p4', type=coords, help='upper left edgepoint')
    parser.add_argument('-plot', help='show a plot of dem and boundaries', action="store_true")
    parser.add_argument('-walls', help='add walls to the stl (e.g. for cfmesh)', action="store_true")
    parser.add_argument('-mesh', help='create mesh instead of stl', action="store_true")
    parser.add_argument('-fillup', help='fill up nan values with neatrest neighbor', action="store_true")
    parser.add_argument('-quad', help='create mesh with quads', action="store_true")
    parser.add_argument('-trionly', help='create mesh with triangles only', action="store_true")
    parser.add_argument('-exactcopy', help='create a triangulation with the resolution of the dem', action="store_true")
    parser.add_argument('-butterflyangle', type=float, help='maximum allowed butterfly angle', default=5)
    parser.add_argument('-offsetx', type=float, help='offset mesh by x m', default=0)
    parser.add_argument('-offsety', type=float, help='offset mesh by y m', default=0)
    args = parser.parse_args()
    
    infile = args.i
    outfile = args.o
    if args.xres is not None:
        xres = args.xres
    if args.yres is not None:
        yres = args.yres
    p1in = args.p1
    p2in = args.p2
    p3in = args.p3
    p4in = args.p4
    walls = args.walls
    plot = args.plot
    mesh = args.mesh
    fillup = args.fillup
    quad = args.quad
    trionly = args.trionly
    exactcopy = args.exactcopy
    butterflyangle = args.butterflyangle
    if exactcopy:
        xres = 0
        yres = 0
        walls = False
        mesh = False


    if p1in is not None:
        p1in = np.array(p1in)
    if p2in is not None:
        p2in = np.array(p2in)
    if p3in is not None:
        p3in = np.array(p3in)
    if p4in is not None:
        p4in = np.array(p4in)


    print("Reading from {}".format(infile))
    print("Writing to {}".format(outfile))
    if walls:     print("adding walls to stl")
    if plot:      print("plotting overview")
    if mesh:      print("creating OpenFOAM mesh directly")
    if exactcopy: print("creating an exact copy of the ascii")
    if fillup:    print("filling up nodata values")
    

    print("moving by ({}, {})".format(args.offsetx, args.offsety))

    f = open(infile, 'r')
    ncols = int(rx.findall(f.readline())[0])
    nrows = int(rx.findall(f.readline())[0])
    xllcenter = float(rx.findall(f.readline())[0])
    yllcenter = float(rx.findall(f.readline())[0])
    cs = float(rx.findall(f.readline())[0])

    #print(xllcenter, yllcenter)
    xllcenter = xllcenter + args.offsetx
    yllcenter = yllcenter + args.offsety

    x = [float(xllcenter+i*cs) for i in range(ncols)]
    y = [float(yllcenter+(nrows-i)*cs) for i in range(nrows)]

    i = 0
    z = np.zeros((nrows, ncols))
    for line in f:
        st = rx.findall(line)
        if len(st) == ncols:
            z[i,:] = [float(s) for s in st]
            i = i+1
    #z = (np.loadtxt(infile, delimiter=' ', converters=None, skiprows=6))

    if exactcopy:
        xres = ncols
        yres = nrows

    print("xres = {}".format(xres))
    print("yres = {}".format(yres))

    if fillup:
        fn = True
        k = 1
        while fn and k < 10:
            print("filling up ({}) ...".format(k)); k = k+1
            fn = False
            fixed = 0
            unfixed = 0
            if k%2 == 0:
                irun = [z.shape[0]-1-i for i in range(z.shape[0])]
                jrun = [z.shape[1]-1-j for j in range(z.shape[1])]
            else:
                irun = range(z.shape[0])
                jrun = range(z.shape[1])
            jj = 0
            for i in irun:
                for j in jrun:
                    jj = jj+1
                    if z[i,j] < 0:
                        val = 0.; n = 0
                        if i > 0: 
                            if z[i-1,j] > 0: val = val+z[i-1,j]; n = n+1
                        if i < z.shape[0]-1: 
                            if z[i+1,j] > 0: val = val+z[i+1,j]; n = n+1
                        if j > 0: 
                            if z[i,j-1] > 0: val = val+z[i,j-1]; n = n+1
                        if j < z.shape[1]-1: 
                            if z[i,j+1] > 0: val = val+z[i,j+1]; n = n+1
                        if n>0: z[i,j] = val/n; fixed = fixed+1
                        else: unfixed = unfixed+1
                        fn = True
            print("...checked: {}, fixed: {}, unfixed: {}".format(jj,fixed,unfixed))       
       
    if exactcopy:
        f = lambda x,y: z[int((yllcenter-y)/cs+nrows), int((x-xllcenter)/cs)]    
    else: 
        ff = interpolate.interp2d(x, y, z, kind='cubic')
        dd = 5
        f = ff#lambda x,y: 0.25*ff(x-dd, y-dd)+0.25*ff(x-dd, y+dd)+0.25*ff(x+dd, y+dd)+0.25*ff(x+dd, y-dd)

    minx = min(x)
    maxx = max(x)
    miny = min(y)
    maxy = max(y)

    print("Limits of map:")
    print("    x = [{},{}]".format(minx, maxx))
    print("    y = [{},{}]".format(miny, maxy))


    p1 = [minx+10, miny+10]
    p2 = [minx+10, maxy+10]
    p3 = [maxx-10, maxy+10] 
    p4 = [maxx-10, miny+10]


    if p1in is not None:
        p1 = p1in+np.array([args.offsetx, args.offsety])
    if p2in is not None:
        p2 = p2in+np.array([args.offsetx, args.offsety])
    if p3in is not None:
        p3 = p3in+np.array([args.offsetx, args.offsety])
    if p4in is not None:
        p4 = p4in+np.array([args.offsetx, args.offsety])


    pe = np.array([p1, p2, p3, p4])

    print("Meshing between:")
    print("    p1 = ({},{})".format(pe[0,0],pe[0,1]))
    print("    p2 = ({},{})".format(pe[1,0],pe[1,1]))
    print("    p3 = ({},{})".format(pe[2,0],pe[2,1]))
    print("    p4 = ({},{})".format(pe[3,0],pe[3,1]))
    s = 0.5*np.abs(np.dot(pe[:,0],np.roll(pe[:,1],1))-np.dot(pe[:,1],np.roll(pe[:,0],1)))
    print("Surface = {}".format(s))
    pm = [0.5*pe[0,0]+0.5*pe[2,0], 0.5*pe[0,1]+0.5*pe[2,1]]
    print("Point inside = ({},{},{})".format(pm[0], pm[1], 10+f(pm[0], pm[1])))
    print("Suggested points for blockMesh:")
    norm = lambda v: v/np.linalg.norm(v)
    ps = [norm(pe[i+1]-pe[i])+norm(pe[i-1]-pe[i]) for i in range(3)]
    ps.append(norm(pe[0]-pe[3])+norm(pe[2]-pe[3]))
    ps = np.array([ps[i]*10.+pe[i] for i in range(4)])
    print("    p1' = ({},{})".format(ps[0,0], ps[0,1]))
    print("    p2' = ({},{})".format(ps[1,0], ps[1,1]))
    print("    p3' = ({},{})".format(ps[2,0], ps[2,1]))
    print("    p4' = ({},{})".format(ps[3,0], ps[3,1]))

    if plot:
        fig, ax = plt.subplots()
        z[z==-9999] = float('nan')
        plt.imshow(z, extent=[minx, maxx, miny, maxy], origin='upper', aspect='equal')

        ax.fill(pe[:,0], pe[:,1], fill=False, alpha=0.5)
        ax.scatter(pe[:,0], pe[:,1])
        ax.fill(ps[:,0], ps[:,1], fill=False, alpha=0.5)
        ax.scatter(ps[:,0], ps[:,1])
        ax.scatter([pm[0]], [pm[1]])
        fig.tight_layout()
        plt.show()


    ztop = [f(p1[0], p1[1])+200, 
            f(p2[0], p2[1])+200, 
            f(p3[0], p3[1])+200, 
            f(p4[0], p4[1])+200]

    f2 = interpolate.interp2d([p1[0], p2[0], p3[0], p4[0]],[p1[1], p2[1], p3[1], p4[1]], ztop)

    pe = np.array([p1, p2, p3, p4])
    p = []


    if exactcopy:
        yglob, xglob = np.meshgrid(y, x)
    else:
        xloc = np.linspace(-1, 1, xres)
        yloc = np.linspace(-1, 1, yres)

        yyloc, xxloc = np.meshgrid(yloc, xloc)

        xglob = []
        yglob = []

        for i in range(len(xxloc)):
            xglob.append([])
            yglob.append([])
            for j in range(len(xxloc[i])):
                p = getX(np.array([xxloc[i][j], yyloc[i][j]]), pe)
                xglob[-1].append(p[0])
                yglob[-1].append(p[1])

    if mesh:

        #1 -> quad
        #2 -> triangle: ul - lr
        #3 -> triangle: ll - ur

        p = np.array([[[xglob[i][j],yglob[i][j],f(xglob[i][j], yglob[i][j])[0]] for j in range(yres)] for i in range(xres)])
        z = np.array([[p[i,j,2] for j in range(yres)] for i in range(xres)])

        facetype = np.zeros(dtype=int, shape=(xres-1,yres-1))

        if quad:
            for i in range(xres-1):
                for j in range(yres-1):
                    facetype[i,j] = 1
        elif trionly:
            for i in range(xres-1):
                for j in range(yres-1):
                    if (i+j)%2 == 0:
                        facetype[i,j] = 2
                    else:
                        facetype[i,j] = 3
        else:
            #upper-left
            ul = [(0,0), (1,0), (0,1)]
            #lower-right
            lr = [(1,0), (1,1), (0,1)]

            #uperr-right
            ur = [(1,0), (1,1), (0,0)]
            #lower-left
            ll = [(0,0), (1,1), (0,1)]

            tri = np.array([ul,lr,ur,ll])

            ul = 0
            lr = 1
            ur = 2
            ll = 3

        n = np.zeros(shape=(4,xres-1,yres-1,3))

        for i in range(xres-1):
            for j in range(yres-1):
                for t in range(4):

                    p0 = p[i+tri[t][0][0],j+tri[t][0][1]]
                    p1 = p[i+tri[t][1][0],j+tri[t][1][1]]
                    p2 = p[i+tri[t][2][0],j+tri[t][2][1]]
                    n[t,i,j] = np.cross(p1-p0, p2-p0)
                    n[t,i,j] /= np.linalg.norm(n[t,i,j])

        kappa = np.zeros(shape=(2,xres-1,yres-1))
        zmid = np.zeros(shape=(xres+1,yres+11))
        ds = np.sqrt(2*(1.5*cs)**2)


        for i in range(xres-1):
            for j in range(yres-1):
                zmid[i,j] = (z[i,j]+z[i+1,j]+z[i,j+1]+z[i+1,j+1])/4.



        limit = np.cos(butterflyangle/180.*np.pi)

        for i in range(xres-1):
            for j in range(yres-1):

                if np.dot(n[ul,i,j], n[lr,i,j]) > limit and np.dot(n[ur,i,j], n[ll,i,j]) > limit:
                    facetype[i,j] = 1

                if i>0 and i<xres-2 and j > 0 and j<yres-2:
                    #ul-lr
                    dzds = (z[i+2,j+2]-z[i-1,j-1])/ds
                    dzds2 = (z[i+2,j+2]-2*zmid[i,j]+z[i-1,j-1])/2./ds
                    kappa[0,i,j] = dzds2/(1+dzds**2)**1.5

                    dzds = (z[i+2,j-1]-z[i-1,j+2])/ds
                    dzds2 = (z[i+2,j-1]-2*zmid[i,j]+z[i-1,j+2])/2./ds
                    kappa[1,i,j] = dzds2/(1+dzds**2)**1.5


                if facetype[i,j] == 1:
                    pass
                elif kappa[0,i,j] == 0 and kappa[1,i,j] == 0:
                    facetype[i,j] = 1
                else:
                    if abs(kappa[0,i,j]) > abs(kappa[1,i,j]):
                        facetype[i,j] = 3
                    else:
                        facetype[i,j] = 2


        l = 0
        ls = []
        npoints = xres*yres*2
    
        nfaces = (xres-1)*(yres-1)*2+(xres-1)*yres+(yres-1)*xres

        for j in range(yres-1):
            for i in range(xres-1):
                if facetype[i,j] != 1:
                    nfaces += 2
        
        print("mesh with {} points, {} faces".format(npoints, nfaces))
        
        fi = open(outfile+"/points", 'w')
        fi.write(header("points", "vectorField"))
        fi.write("{}\n".format(npoints))
        fi.write("(\n")
        for i in range(xres):
            ls.append([])
            for j in range(yres):
                ls[-1].append([])    
                for k in range(2):
                    ls[-1][-1].append(l)
                    l += 1
                    fi.write("({} {} {})\n".format(xglob[i][j], yglob[i][j], k*20.+f(xglob[i][j], yglob[i][j])[0]))
        fi.write(")\n")
        fi.write(footer())
        fi.close()


        owner = [-1 for i in range(nfaces)]
        neighbor = [-1 for i in range(nfaces)]

        fi = open(outfile+"/faces", 'w')
        fi.write(header("faces", "faceList"))
        fi.write("{}\n".format(nfaces))
        fi.write("(\n")

        l = 0
        innerXStart = l
        #Inner Xnormal
        for j in range(yres-1):
            for i in range(1, xres-1):
                fi.write("4({0} {3} {2} {1})//inner_x\n".format(ls[i][j][0], ls[i][j+1][0], ls[i][j+1][1], ls[i][j][1]))
                owner[l] = (i-1)+(xres-1)*j
                neighbor[l] = owner[l]+1
                l = l+1
        innerYStart = l
        #Inner Ynormal
        for j in range(1, yres-1):
            for i in range(xres-1):
                fi.write("4({0} {3} {2} {1})//inner_y\n".format(ls[i][j][0], ls[i][j][1], ls[i+1][j][1], ls[i+1][j][0]))
                owner[l] = i+(xres-1)*(j-1)
                neighbor[l] = owner[l]+(xres-1)
                l = l+1
        bottomStart = l
        #Bottom
        for j in range(yres-1):
            for i in range(xres-1):
                if facetype[i,j] == 1:
                    fi.write("4({0} {1} {2} {3})//bottom(q)\n".format(ls[i][j][0], ls[i+1][j][0], ls[i+1][j+1][0], ls[i][j+1][0]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                elif facetype[i,j] == 2:
                    fi.write("3({1} {0} {2})//bottom(t)\n".format(ls[i][j][0], ls[i+1][j+1][0], ls[i+1][j][0]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                    fi.write("3({1} {0} {2})//bottom(t)\n".format(ls[i][j][0], ls[i][j+1][0], ls[i+1][j+1][0]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                else:
                    fi.write("3({1} {0} {2})//bottom(t)\n".format(ls[i][j][0], ls[i][j+1][0], ls[i+1][j][0]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                    fi.write("3({1} {0} {2})//bottom(t)\n".format(ls[i+1][j][0], ls[i][j+1][0], ls[i+1][j+1][0]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
        topStart = l
        #Top
        for i in range(len(xglob)-1):
            for j in range(len(xglob[i])-1):
                if facetype[i,j] == 1:
                    fi.write("4({0} {1} {2} {3})//top(q)\n".format(ls[i][j][1], ls[i][j+1][1], ls[i+1][j+1][1], ls[i+1][j][1]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                elif facetype[i,j] == 2:
                    fi.write("3({1} {0} {2})//top(t)\n".format(ls[i][j][1], ls[i+1][j][1], ls[i+1][j+1][1]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                    fi.write("3({1} {0} {2})//top(t)\n".format(ls[i][j][1], ls[i+1][j+1][1], ls[i][j+1][1]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                else:
                    fi.write("3({1} {0} {2})//top(t)\n".format(ls[i][j][1], ls[i+1][j][1], ls[i][j+1][1]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
                    fi.write("3({1} {0} {2})//top(t)\n".format(ls[i+1][j][1], ls[i+1][j+1][1], ls[i][j+1][1]))
                    owner[l] = i+(xres-1)*j
                    l = l+1
        xminStart = l
        #Xmin
        for j in range(len(yglob[0])-1):
            fi.write("4({} {} {} {})//xmin\n".format(ls[0][j][0], ls[0][j+1][0], ls[0][j+1][1], ls[0][j][1]))
            owner[l] = (xres-1)*j
            l = l+1
        xmaxStart = l
        #Xmax
        for j in range(len(yglob[-1])-1):
            fi.write("4({} {} {} {})//xmax\n".format(ls[-1][j][0], ls[-1][j][1], ls[-1][j+1][1], ls[-1][j+1][0]))
            owner[l] = (xres-1)*j+xres-2
            l = l+1
        yminStart = l
        #Ymin
        for i in range(len(xglob)-1):
            fi.write("4({} {} {} {})//ymin\n".format(ls[i][0][0], ls[i][0][1], ls[i+1][0][1], ls[i+1][0][0]))
            owner[l] = i
            l = l+1
        ymaxStart = l
        #Ymax
        for i in range(len(xglob)-1):
            fi.write("4({} {} {} {})//ymax\n".format(ls[i][-1][0], ls[i+1][-1][0], ls[i+1][-1][1], ls[i][-1][1]))
            owner[l] = (xres-1)*(yres-2)+i
            l = l+1
        end = l

        fi.write(")\n")
        fi.write(footer())
        fi.close()


        fi = open(outfile+"/owner", 'w')
        fi.write(header("owner", "labelList"))
        fi.write("{}\n".format(nfaces))
        fi.write("(\n")

        for kk in owner:
            fi.write("{}\n".format(kk))
                    
        fi.write(")\n")
        fi.write(footer())
        fi.close()



        fi = open(outfile+"/neighbour", 'w')
        fi.write(header("neighbour", "labelList"))
        fi.write("{}\n".format(nfaces))
        fi.write("(\n")

        for kk in neighbor:
            fi.write("{}\n".format(kk))
                    
        fi.write(")\n")
        fi.write(footer())
        fi.close()



        fi = open(outfile+"/boundary", 'w')
        fi.write(header("boundary", "polyBoundaryMesh"))

        fi.write("6\n")
        fi.write("(\n")


        fi.write("    terrain\n")
        fi.write("    {\n")
        fi.write("        type            patch;\n")
        fi.write("        nFaces          {};\n".format(topStart-bottomStart))
        fi.write("        startFace       {};\n".format(bottomStart))
        fi.write("    }\n")

        fi.write("    maxZ\n")
        fi.write("    {\n")
        fi.write("        type            patch;\n")
        fi.write("        nFaces          {};\n".format(xminStart-topStart))
        fi.write("        startFace       {};\n".format(topStart))
        fi.write("    }\n")

        fi.write("    minX\n")
        fi.write("    {\n")
        fi.write("        type            patch;\n")
        fi.write("        nFaces          {};\n".format(xmaxStart-xminStart))
        fi.write("        startFace       {};\n".format(xminStart))
        fi.write("    }\n")

        fi.write("    maxX\n")
        fi.write("    {\n")
        fi.write("        type            patch;\n")
        fi.write("        nFaces          {};\n".format(yminStart-xmaxStart))
        fi.write("        startFace       {};\n".format(xmaxStart))
        fi.write("    }\n")

        fi.write("    minY\n")
        fi.write("    {\n")
        fi.write("        type            patch;\n")
        fi.write("        nFaces          {};\n".format(ymaxStart-yminStart))
        fi.write("        startFace       {};\n".format(yminStart))
        fi.write("    }\n")

        fi.write("    maxY\n")
        fi.write("    {\n")
        fi.write("        type            patch;\n")
        fi.write("        nFaces          {};\n".format(end-ymaxStart))
        fi.write("        startFace       {};\n".format(ymaxStart))
        fi.write("    }\n")

        fi.write(")\n")
        fi.write(footer())
        fi.close()


    else:
        stl = open(outfile, 'w')
        stl.write('solid terrain\n')
        print("writing terrain...");
        for i in range(len(xglob)-1):
            if i%100 == 0:
                print("...line {}".format(i))
            for j in range(len(xglob[i])-1):
                if exactcopy:
                    writeFlag = (f(xglob[i][j],yglob[i][j]) > -9990 and f(xglob[i+1][j],yglob[i+1][j]) > -9990 and f(xglob[i+1][j+1],yglob[i+1][j+1]) > -9990 and  f(xglob[i][j+1], yglob[i][j+1]) > -9990)
                else:
                    writeFlag = True
                if writeFlag:
                    vertices = [vertex(xglob[i][j], yglob[i][j], f(xglob[i][j], yglob[i][j])),
                                vertex(xglob[i][j+1], yglob[i][j+1], f(xglob[i][j+1], yglob[i][j+1])),
                                vertex(xglob[i+1][j+1], yglob[i+1][j+1], f(xglob[i+1][j+1], yglob[i+1][j+1])),
                                vertex(xglob[i+1][j], yglob[i+1][j], f(xglob[i+1][j], yglob[i+1][j]))]
                    for k in range(2):
                        if (i+j)%2 == 0:
                            stl.write(str(surf(vertices[0], vertices[k+1], vertices[k+2])))
                        else:
                            stl.write(str(surf(vertices[k+0], vertices[k+1], vertices[3])))
        stl.write('endsolid terrain\n')


        if walls:
            stl.write('solid maxZ\n')
            for i in range(len(xglob)-1):
                for j in range(len(xglob[i])-1):
                    vertices = [vertex(xglob[i][j], yglob[i][j], f2(xglob[i][j], yglob[i][j])),
                                vertex(xglob[i][j+1], yglob[i][j+1], f2(xglob[i][j+1], yglob[i][j+1])),
                                vertex(xglob[i+1][j+1], yglob[i+1][j+1], f2(xglob[i+1][j+1], yglob[i+1][j+1])),
                                vertex(xglob[i+1][j], yglob[i+1][j], f2(xglob[i+1][j], yglob[i+1][j]))]
                    for k in range(2):
                        stl.write(str(surf(vertices[0], vertices[k+1], vertices[k+2])))
            stl.write('endsolid maxZ\n')


            stl.write('solid minX\n')
            i = 0
            for j in range(len(xglob[i])-1):
                vertices = [vertex(xglob[i][j], yglob[i][j], f(xglob[i][j], yglob[i][j])),
                            vertex(xglob[i][j+1], yglob[i][j+1], f(xglob[i][j+1], yglob[i][j+1])),
                            vertex(xglob[i][j+1], yglob[i][j+1], f2(xglob[i][j+1], yglob[i][j+1])),
                            vertex(xglob[i][j], yglob[i][j], f2(xglob[i][j], yglob[i][j]))]
                for k in range(2):
                    stl.write(str(surf(vertices[0], vertices[k+1], vertices[k+2])))
            stl.write('endsolid minX\n')

            stl.write('solid maxX\n')
            i = -1
            for j in range(len(xglob[i])-1):
                vertices = [vertex(xglob[i][j], yglob[i][j], f(xglob[i][j], yglob[i][j])),
                            vertex(xglob[i][j+1], yglob[i][j+1], f(xglob[i][j+1], yglob[i][j+1])),
                            vertex(xglob[i][j+1], yglob[i][j+1], f2(xglob[i][j+1], yglob[i][j+1])),
                            vertex(xglob[i][j], yglob[i][j], f2(xglob[i][j], yglob[i][j]))]
                for k in range(2):
                    stl.write(str(surf(vertices[0], vertices[k+1], vertices[k+2])))
            stl.write('endsolid maxX\n')


            stl.write('solid minY\n')
            j = 0
            for i in range(len(xglob)-1):
                vertices = [vertex(xglob[i][j], yglob[i][j], f(xglob[i][j], yglob[i][j])),
                            vertex(xglob[i][j], yglob[i][j], f2(xglob[i][j], yglob[i][j])),
                            vertex(xglob[i+1][j], yglob[i+1][j], f2(xglob[i+1][j], yglob[i+1][j])),
                            vertex(xglob[i+1][j], yglob[i+1][j], f(xglob[i+1][j], yglob[i+1][j]))]
                for k in range(2):
                    stl.write(str(surf(vertices[0], vertices[k+1], vertices[k+2])))
            stl.write('endsolid minY\n')

            stl.write('solid maxY\n')
            j = -1
            for i in range(len(xglob)-1):
                vertices = [vertex(xglob[i][j], yglob[i][j], f(xglob[i][j], yglob[i][j])),
                            vertex(xglob[i][j], yglob[i][j], f2(xglob[i][j], yglob[i][j])),
                            vertex(xglob[i+1][j], yglob[i+1][j], f2(xglob[i+1][j], yglob[i+1][j])),
                            vertex(xglob[i+1][j], yglob[i+1][j], f(xglob[i+1][j], yglob[i+1][j]))]
                for k in range(2):
                    stl.write(str(surf(vertices[0], vertices[k+1], vertices[k+2])))
            stl.write('endsolid maxY\n')


        stl.close()


if __name__ == "__main__":
   main(sys.argv[1:])
