##Geometric properties of a rectangular patch
##Orthogonal geometry with line roll=1
###Hablemos de s¡smica - Jaime Checa J.



import array
import numpy
import matplotlib.pyplot as plt
import math
# Variables defining the acquisition geometry:
#   nl:Number of active lines in the pacth (even number)
#   ncl:number of active stations per line
#   dr, ds: Receiver and source interval
#   dl,dls: Receiver and source line interval
nl=28
ncl=140
dr=40
ds=80
dl=240
dls=400

#Inicialization of variables

salvos=int(dl/ds)
o = numpy.zeros(salvos*ncl*nl+1)
az= numpy.zeros(salvos*ncl*nl+1)

omil=dr*(ncl-1)/2
omxl=nl/2*dl-ds/2

#Draw receiver lines
m=0
plt.gca().set_aspect('equal', adjustable='box')
for j in range(1,nl+1):
 for k in range(1,ncl+1):
      xr=(k-1)*dr
      yr=(j-1)*dl
      plt.plot(xr,yr,marker='.',color='blue')
      m=m+1


#Function to calculate azimut in degrees
def azimuth(x1,y1,x2,y2):
  """Calcula el azimut en grados de una linea entre dos puntos"""
  dx=x2-x1
  dy=y2-y1
  azi=abs(math.degrees(math.atan(dx/dy)))
  if dy >0 and dx >0:
    azi=azi
  if dy<0 and dx>0:
    azi =180-azi
  if dy<0 and dx<0:
    azi+=180
  if dy>0 and dx<0:
    azi=360-azi
  return azi

# Draw salvo
# Calculates offset and azimut of each trace and save them 
# in o and az vectors.
# xs,ys: Source coordinates, xr,yr: Receiver coordinates

m=0
for i in range(1, salvos+1):
  xs=(ncl-1)/2*dr
  ys=(nl/2-1)*dl+ds/2+(i-1)*ds
  plt.plot(xs,ys,marker='.',color='red')
for i in range(1, salvos+1):
  xs=(ncl-1)/2*dr
  ys=(nl/2-1)*dl+ds/2+(i-1)*ds

  for j in range(1,nl+1):
    for k in range(1,ncl+1):
      xr=(k-1)*dr
      yr=(j-1)*dl
      o[m]=int(((xr-xs)**2+(yr-ys)**2))**0.5
      az[m]=azimuth(xs,ys,xr,yr)
      m=m+1

# Prints patch properties
print ('Geometrya of the Patch')
print ('-------------------')
print ('Number of lines: ',nl)
print ('Receiver line interval: ',dl)
print ('Receiver interval: ',dr)
print ('Source line interval: ',dls)
print ('Source interval: ',ds)
print ('Active stations per line: ',ncl)
print ('Number of channels in the patch: ',ncl*nl)
print ('-------------------')
print ('Nominal Fold:',ncl/2/dls*dr*nl/2,' in a bin of ',int(dr/2),' x ',int(ds/2))
print ('Nominal trace density:',int(ncl/2/dls*dr*nl/2/dr/ds*4*1000000) )
print ('Sources per salvo:',salvos)
print ('Source density:',int(1000/ds*1000/dls),' Sources/km2')
print ('Maximum inline offset: ',int(omil))
print ('Maximum xline offset:',int(omxl))
print ('Total maximum offset:',int((omil**2+omxl**2)**0.5))
print ('Maximum minimum offset:',int(((dls**2+dl**2))**0.5))
print ('Aspect ratio xl/il: %5.2f' % (omxl/omil))
print ('Lenght of receiver lines %5.1f'%(1000/dl), 'km/km2')
print ('Length of source lines: %5.1f'%(1000/dls), 'km/km2')
print ('Roll on xline: %5.1f' %(dl*(nl/4-0.5)))
print ('Roll on inline: %5.1f'% (dls*(ncl/2/dls*dr*0.5-0.5)))

# Offset and Azimut histogram
# 52 bins by default

# Plots cummulative offset histogram.
# Cand be also used to check fold vs offset by 
#multiplying density times area of subsurface bin.
x=plt.hist(o, bins=52, cumulative=True)
plt.xlabel('Offset (m)')
plt.ylabel('Number of traces')
plt.title("Cummulative offset histogram")
plt.grid(axis='y')
plt.grid(axis='x')
plt.show()
a=x[0]
b=x[1]
b=b[0:52]
densidad=ncl/2/dls*dr*nl/2/dr/ds*4*1000000
escala=densidad/numpy.max(a)
sal=(b,a*escala)
numpy.savetxt('cummulativehistogram.csv',numpy.transpose(sal),fmt='%d', delimiter=",")

plt.plot (b,a*escala)
plt.xlabel('Offset (m)')
plt.ylabel('Trace density. tr/km2')
plt.title("Cummulative trace density")
plt.grid(axis='y')
plt.grid(axis='x')
plt.show()

# Plots offset histogram showing the relative importance of offset values.
x=plt.hist(o, bins=52, cumulative=False)
a=x[0]
b=x[1]
b=b[0:52]
sal=(b,a)
numpy.savetxt('histograma.csv',numpy.transpose(sal),fmt='%d', delimiter=",")
plt.xlabel('Offset (m)')
plt.ylabel('Cantidad de trazas')
plt.title("Histograma de offsets")
plt.grid(axis='y')
plt.grid(axis='x')
plt.show()

# Grafica la dstribuci¢n de offsets vs Azimut
plt.plot (az,o,'*', color='orange')
plt.xlabel('Azimut in degrees')
plt.ylabel('Offset in m')
plt.title("Azimut vs offset distribution")
plt.grid(axis='y')
plt.grid(axis='x')
plt.xlabel('Angle (degrees)')
plt.ylabel('Number of traces')
plt.title("Offset histogram")
plt.grid(axis='y')
plt.grid(axis='x')
plt.show()