From 7aedc35c164f1eb38d2249ac48149008e70191f4 Mon Sep 17 00:00:00 2001
From: Marcel Paffrath <marcel@minos01.geophysik.ruhr-uni-bochum.de>
Date: Mon, 9 Nov 2015 16:12:12 +0100
Subject: [PATCH] first implementation of a checkerboard modification for a
 given vgrids.in file

---
 pylot/core/active/seismicArrayPreparation.py | 130 +++++++++++++++++++
 1 file changed, 130 insertions(+)

diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py
index a26b659c..2a22dbcb 100644
--- a/pylot/core/active/seismicArrayPreparation.py
+++ b/pylot/core/active/seismicArrayPreparation.py
@@ -511,6 +511,136 @@ class SeisArray(object):
         print('Wrote %d points to file %s for %d layers'%(count, outfilename, nlayers))
         outfile.close()
 
+    def addCheckerboard(self, spacing = 20, pertubation = 1, inputfile = 'vgrids.in', outputfile = 'vgrids_cb.in'):
+        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])
+
+            print('readNumberOf Points: Awaiting %d grid points in %s'
+                  %(nR*nTheta*nPhi, filename))
+            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):
+            '''            
+            Reads in velocity from vgrids file and returns a list containing all values in the same order
+            '''
+            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
+
+        def correctSpacing(spacing, delta, disttype = None):
+            if spacing > delta:
+                spacing_corr = round(spacing / delta) * delta
+            elif spacing < delta:
+                spacing_corr = delta
+            print('The spacing of the checkerboard of %s (%s) was corrected to '
+                  'a value of %s to fit the grid spacing of %s.' %(spacing, disttype, spacing_corr, delta))
+            return spacing_corr
+            
+
+        R = 6371. # earth radius
+        decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3)
+        outfile = open(outputfile, 'w')
+
+        # Theta, Phi in radians, R in km
+        nR, nTheta, nPhi = readNumberOfPoints(inputfile)
+        dR, dThetaRad, dPhiRad = readDelta(inputfile)
+        sR, sThetaRad, sPhiRad = readStartpoints(inputfile)
+        vel = readVelocity(inputfile)
+
+        dTheta, dPhi = np.rad2deg((dThetaRad, dPhiRad))
+        sTheta, sPhi = np.rad2deg((dThetaRad, dPhiRad))
+        
+        eR = sR + (nR - 1) * dR
+        ePhi = sPhi + (nPhi - 1) * dPhi
+        eTheta = sTheta + (nTheta - 1) * dTheta
+
+        nPoints = nR * nTheta * nPhi
+
+        thetaGrid = np.linspace(sTheta, eTheta, num = nTheta)
+        phiGrid = np.linspace(sPhi, ePhi, num = nPhi)
+        rGrid = np.linspace(sR, eR, num = nR)
+
+        # write header for velocity grid file (in RADIANS)
+        outfile.writelines('%10s %10s \n' %(1, 1))
+        outfile.writelines('%10s %10s %10s\n' %(nR, nTheta, nPhi))
+        outfile.writelines('%10s %10s %10s\n' %(dR, dThetaRad, dPhiRad))
+        outfile.writelines('%10s %10s %10s\n' %(sR, sThetaRad, sPhiRad))
+
+        spacR = correctSpacing(spacing, dR, '[meter], R')
+        spacTheta = correctSpacing(self._getAngle(spacing), dTheta, '[degree], Theta')
+        spacPhi = correctSpacing(self._getAngle(spacing), dPhi, '[degree], Phi')
+
+        count = 0
+        evenOdd = 1
+        even = 0; odd = 0
+
+        # In the following loop it is checked whether the positive distance from the border of the model
+        # for a point on the grid divided by the spacing is even or odd and then pertubated.
+        # The position is also shifted by half of the delta so that the position is directly on the point and
+        # not on the border between two points.
+        for radius in rGrid:
+            if np.floor((radius - sR - dR/2) / spacR) % 2:
+                evenOddR = 1
+            else:
+                evenOddR = -1
+            for theta in thetaGrid:
+                if np.floor((theta - sTheta - dTheta/2) / spacTheta) % 2:
+                    evenOddT = 1
+                else:
+                    evenOddT = -1
+                for phi in phiGrid:
+                    if np.floor((phi - sPhi - dPhi/2) / spacPhi) % 2:
+                        evenOddP = 1
+                    else:
+                        evenOddP = -1
+                    velocity = vel[count]
+                    evenOdd = evenOddR * evenOddT * evenOddP
+                    velocity += evenOdd * pertubation
+
+                    outfile.writelines('%10s %10s\n'%(velocity, decm))
+                    count += 1
+
+                    progress = float(count) / float(nPoints) * 100
+                    self._update_progress(progress)
+
+        outfile.close()
+
     def exportAll(self, filename = 'interpolated_receivers.out'):
         recfile_out = open(filename, 'w')
         count = 0