From f7878cdb4aca4982e93bc5ca46ace32bbe52b2cc Mon Sep 17 00:00:00 2001
From: Marcel Paffrath <marcel@minos02.geophysik.ruhr-uni-bochum.de>
Date: Mon, 28 Sep 2015 10:20:26 +0200
Subject: [PATCH] implemented vgrids2VTK to surveyUtils

---
 pylot/core/active/surveyUtils.py | 112 ++++++++++++++++++++++++++++++-
 1 file changed, 111 insertions(+), 1 deletion(-)

diff --git a/pylot/core/active/surveyUtils.py b/pylot/core/active/surveyUtils.py
index 2f1c990b..a7086102 100644
--- a/pylot/core/active/surveyUtils.py
+++ b/pylot/core/active/surveyUtils.py
@@ -193,6 +193,116 @@ def findTracesInRanges(shot_dict, distancebin, pickbin):
                     shots_found[shot.getShotnumber()].append(traceID)
 
     return shots_found
-            
         
+def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'):
+    def getDistance(angle):
+        PI = np.pi
+        R = 6371.
+        distance = angle / 180 * (PI * R)
+        return distance
+
+    def readNumberOfPoints(filename):
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        nR = int(vglines[1].split()[0])
+        nTheta = int(vglines[1].split()[1])
+        nPhi = int(vglines[1].split()[2])
+
+        fin.close()
+        return nR, nTheta, nPhi
+
+    def readDelta(filename):
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        dR = float(vglines[2].split()[0])
+        dTheta = float(vglines[2].split()[1])
+        dPhi = float(vglines[2].split()[2])
+
+        fin.close()
+        return dR, dTheta, dPhi
+
+    def readStartpoints(filename):
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        sR = float(vglines[3].split()[0])
+        sTheta = float(vglines[3].split()[1])
+        sPhi = float(vglines[3].split()[2])
+
+        fin.close()
+        return sR, sTheta, sPhi
+
+    def readVelocity(filename):
+        vel = []; count = 0
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        for line in vglines:
+            count += 1
+            if count > 4:
+                vel.append(float(line.split()[0]))
+
+        print("Read %d points out of file: %s" %(count - 4, filename))
+        return vel
+
+    R = 6371 # earth radius
+    outfile = open(outputfile, 'w')
+
+    # Theta, Phi in radians, R in km
+    nR, nTheta, nPhi = readNumberOfPoints(inputfile)
+    dR, dTheta, dPhi = readDelta(inputfile)
+    sR, sTheta, sPhi = readStartpoints(inputfile)
+    vel = readVelocity(inputfile)
+
+    nX = nPhi; nY = nTheta; nZ = nR
+
+    sZ = sR - R
+    sX = getDistance(np.rad2deg(sPhi))
+    sY = getDistance(np.rad2deg(sTheta))
+
+    dX = getDistance(np.rad2deg(dPhi))
+    dY = getDistance(np.rad2deg(dTheta))
+    
+    xGrid = np.linspace(sX, sX + (dX * nX), nX)
+    yGrid = np.linspace(sZ, sZ + (nY * dY), nY)
+    zGrid = np.linspace(sZ, sZ + (nR * dR), nR)
+    nPoints = len(xGrid) * len(yGrid) * len(zGrid)
+
+    # write header
+    print("Writing header for VTK file...")
+    outfile.writelines('# vtk DataFile Version 3.1\n')
+    outfile.writelines('Velocity on FMTOMO vgrids.in points\n')
+    outfile.writelines('ASCII\n')
+    outfile.writelines('DATASET POLYDATA\n')
+    outfile.writelines('POINTS %15d float\n' %(nPoints))
+
+    # write coordinates
+    print("Writing coordinates to VTK file...")
+    for z in zGrid:
+        for y in yGrid:
+            for x in xGrid:
+                outfile.writelines('%10f %10f %10f \n' %(x, y, z))
+
+
+    outfile.writelines('VERTICES %15d %15d\n' %(nPoints, 2 * nPoints))
+
+    # write indices
+    print("Writing indices to VTK file...")
+    for index in range(nPoints):
+        outfile.writelines('%10d %10d\n' %(1, index))
+
+    outfile.writelines('POINT_DATA %15d\n' %(nPoints))
+    outfile.writelines('SCALARS velocity float %d\n' %(1))
+    outfile.writelines('LOOKUP_TABLE default\n')
+
+    # write velocity
+    print("Writing velocity values to VTK file...")
+    for velocity in vel:
+        outfile.writelines('%10f\n' %velocity)
+
+    outfile.close()
+    print("Wrote velocity grid for %d points to file: %s" %(nPoints, outputfile))
+    return