#!/usr/bin/env python
# ----------------------------------------------------------------------
#
# Brad T. Aagaard, U.S. Geological Survey
# Charles A. Williams, GNS Science
# Matthew G. Knepley, University of Chicago
#
# This code was developed as part of the Computational Infrastructure
# for Geodynamics (http://geodynamics.org).
#
# Copyright (c) 2010-2015 University of California, Davis
#
# See COPYING for license information.
#
# ----------------------------------------------------------------------
#
# PREREQUISITES: numpy

# ======================================================================
import numpy
import sys

for i in range(len(sys.argv)):                        # read parameter from shell(contain the prefix)
  if(sys.argv[i].find("-p")==0):
    P1=sys.argv[i+1]
    P2=sys.argv[i+2]
    P3=sys.argv[i+3]
    P4=sys.argv[i+4]
    P5=sys.argv[i+5]

print P1,P2,P3,P4,P5
# Define parameters.
interfaceFile = P1
numContours = int(P3)
pointsPerContour = int(P2)
journalFile = P4
acisFile = P5

# Journal file formatting, etc.

# Read coordinates and reshape them.
intCoords = numpy.loadtxt(interfaceFile, dtype=numpy.float64).reshape(numContours, pointsPerContour, 3)

j = open(journalFile, 'w')
j.write("# CUBIT journal file generated by interface_netsurf.py.\n" + \
            "#\n" + \
            "# Create an ACIS NURBS surface from intersecting lines.\n" + \
            "#\n" \
            "# ----------------------------------------------------------\n" + \
            "reset\n")

# Loop over contours (u-lines).
for iContour in range(numContours):
    points = intCoords[iContour,:,:]
    for iPoint in range(pointsPerContour):
        point = points[iPoint,:]
        j.write("create vertex x %10.5e y %10.5e z %10.5e\n" % \
                    (point[0], point[1], point[2]))
        if 0 == iPoint:
            j.write("${idBeg=Id('vertex')}\n")
    j.write("${idEnd=Id('vertex')}\n")
    j.write("create curve spline vertex {idBeg} to {idEnd} delete\n")
    if 0 == iContour:
        j.write("${idCBeg=Id('curve')}\n")
j.write("${idCEnd=Id('curve')}\n\n")

# Loop over profiles (v-lines).
for iProfile in range(pointsPerContour):
    points = intCoords[:,iProfile,:]
    for iPoint in range(numContours):
        point = points[iPoint,:]
        j.write("create vertex x %10.5e y %10.5e z %10.5e\n" % \
                    (point[0], point[1], point[2]))
        if 0 == iPoint:
            j.write("${idBeg=Id('vertex')}\n")
    j.write("${idEnd=Id('vertex')}\n")
    j.write("create curve spline vertex {idBeg} to {idEnd} delete\n")
    if 0 == iProfile:
        j.write("${idPBeg=Id('curve')}\n")
j.write("${idPEnd=Id('curve')}\n\n")

# Create net surface.
j.write("create surface net u curve {idCBeg} to {idCEnd} v curve {idPBeg} to {idPEnd}\n")

## Rotate the surface
#j.write("\n")
#length = numpy.max(intCoords[:,:,0]) - numpy.min(intCoords[:,:,0])
#width  = numpy.max(intCoords[:,:,1]) - numpy.min(intCoords[:,:,1])
#offset_x = (numpy.max(intCoords[:,:,0]) + numpy.min(intCoords[:,:,0]))/2.0
#offset_y = (numpy.max(intCoords[:,:,1]) + numpy.min(intCoords[:,:,1]))/2.0
#print length,width, offset_x, offset_y
#
#
#create surface rectangle width length  height width
#surface 2 move x offset_x y offset_y z 0
#surface 2 move x 0 y 0 z -3000 copy
#surface 2 move x 0 y 0 z -30000
#

# Delete spline curves and export Acis file.
j.write("delete curve all\n")
j.write("export Acis '%s' overwrite\n" % acisFile)
j.close()

# End of file
