#CELDA 1
import array
import numpy
import matplotlib.pyplot as plt
import math
# Entrada de las 6 variables que definen la geometria del patch rectangular.
#   nl:Nœmero de l’neas activas en el pacth (se asume nœmero par)
#   ncl:nœmero de estaciones activas en cada l’nea
#   dr, ds: distancia entre receptores y entre fuentes rspectivamente
#   dl,dls: distancia entre lineas de receptores y lineas de fuentes resp.
nl=10
ncl=120
dr=40
ds=80
dl=320
dls=400

#Inicializaci—n de 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-1)*dl+ds/2

#Dibuja las lineas de receptores
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

#Funci—n para calcular azimut en grados
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

# Dibuja el salvo
# Calcula el offset y azimut de cada traza que genera el Patch y las almacena
# en los vectores o y az respectivamente
# xs,ys: Coordenadas de la fuente, xr,yr: Coordenadas del receptor

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

# Reporta las caracter’sticas del Patch y sus propiedades principales
print ('Geometria del Patch')
print ('-------------------')
print ('Numero de lineas: ',nl)
print ('Distancia entre lineas de receptores: ',dl)
print ('Distancia entre receptores: ',dr)
print ('Distancia entre lineas de disparo: ',dls)
print ('Distancia entre disparos: ',ds)
print ('Numero de estaciones activas por linea: ',ncl)
print ('Numero de canales en el patch: ',ncl*nl)
print ('-------------------')
print ('Fold nominal:',ncl/2/dls*dr*nl/2,' en bin de ',int(dr/2),' x ',int(ds/2))
print ('Densidad de trazas nominal:',int(ncl/2/dls*dr*nl/2/dr/ds*4*1000000) )
print ('Disparos por salvo:',salvos)
print ('Densidad de disparos:',int(1000/ds*1000/dls),' Disp/km2')
print ('Offset maximo inline: ',int(omil))
print ('Offset maximo xline:',int(omxl))
print ('Offset maximo total:',int((omil**2+omxl**2)**0.5))
print ('Maximo offset minimo:',int(((dls**2+dl**2))**0.5))
print ('Relacion de aspecto xl/il: %5.2f' % (omxl/omil))
print ('Kilometros de linea de receptores %5.1f'%(1000/dl), 'km/km2')
print ('Kilometros de linea de fuentes: %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)))


#CELDA 2
# Histograma de offsets y azimuts
# El histograma esta dividido en 52 celdas

# Grafica el histograma acumulado de offsets.
# Con un factor de escala representa la densidad acumulada de trazas 
# con respecto al offset o el cubrimiento en el subsuelo. 
x=plt.hist(o, bins=52, cumulative=True)
plt.xlabel('Offset (m)')
plt.ylabel('Cantidad de trazas')
plt.title("Histograma acumulado de offsets")
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('histograma acumulado.csv',numpy.transpose(sal),fmt='%d', delimiter=",")

plt.plot (b,a)
plt.xlabel('Offset (m)')
plt.ylabel('Densidad de trazas. tr/km2')
plt.title("Densidad de trazas acumulada vs offset")
plt.grid(axis='y')
plt.grid(axis='x')
plt.show()

# Grafica el histograma de offsets que muestra la importancia relativa de los
# diferentes valores de offset que aporta este patch.
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 en grados')
plt.ylabel('Offset en m')
plt.title("Distribucion de azimut vs offset")
plt.grid(axis='y')
plt.grid(axis='x')
plt.xlabel('Offset (m)')
plt.ylabel('Cantidad de trazas')
plt.title("Histograma de offsets")
plt.grid(axis='y')
plt.grid(axis='x')
plt.show()